Content area
A class of sparse optimization techniques that require solely matrix-vector products, rather than an explicit access to the forward matrix and its transpose, has been paid much attention in the recent decade for dealing with large-scale inverse problems. This study tailors application of the so-called Gradient Projection for Sparse Reconstruction (GPSR) to large-scale time-difference three-dimensional electrical impedance tomography (3D EIT). 3D EIT typically suffers from the need for a large number of voxels to cover the whole domain, so its application to real-time imaging, for example monitoring of lung function, remains scarce since the large number of degrees of freedom of the problem extremely increases storage space and reconstruction time. This study shows the great potential of the GPSR for large-size time-difference 3D EIT. Further studies are needed to improve its accuracy for imaging small-size anomalies.
http://crossmark.crossref.org/dialog/?doi=10.1007/s11517-015-1441-1&domain=pdf
Web End = Med Biol Eng Comput (2016) 54:12431255 DOI 10.1007/s11517-015-1441-1
ORIGINAL ARTICLE
http://crossmark.crossref.org/dialog/?doi=10.1007/s11517-015-1441-1&domain=pdf
Web End = A fast timedifference inverse solver for 3D EIT with application to lung imaging
Ashkan Javaherian1 Manuchehr Soleimani2 Knut Moeller1
Received: 18 December 2014 / Accepted: 20 November 2015 / Published online: 6 January 2016 International Federation for Medical and Biological Engineering 2015
1 Introduction
Electrical impedance tomography (EIT) is a diffusive imaging modality for reconstructing the conductivity eld inside an object from surface electrical measurements [15, 32]. This technique has many applications in medicine, e.g., real-time monitoring of lung [2, 3], detection of breast tumors [30], or imaging of brain activity [4]. Typically, a number of electrodes are attached to the skin of the subject, a small alternating current is injected through some of these electrodes successively, and the induced electrical potentials are measured on the remaining electrodes [15, 32].
The image reconstruction is done by iteratively updating the conductivity eld until 2 norm of discrepancy between simulated and real measured data is minimized. From a theoretical point of view, it involves alternatively a nonlinear forward problem of calculating the surface voltages from the conductivity eld via nite element method (FEM) and a severely ill-posed inverse problem for updating the conductivity eld from the surface voltages [40, 45, 56]. To cope with the nonlinear relationship between the conductivity eld and the datasets, a Jacobian (sensitivity) matrix is computed to linearize the problem around a homogenous conductivity [40, 56]. To deal with the high ill-posedness of the problem, the inverse problem is often regularized via assuming a priori assumptions about the conductivity eld [1, 14, 50].
EIT is an efcient tool for real-time monitoring of lung since there is a large contrast between the conductivity of air and that of the encompassing tissues [3]. There are, however, some sources of error during the measurement. Breathing action or posture changes may move the electrodes during the measurements, thus deleteriously affecting the recovered conductivity [5, 47].
Abstract A class of sparse optimization techniques that require solely matrixvector products, rather than an explicit access to the forward matrix and its transpose, has been paid much attention in the recent decade for dealing with large-scale inverse problems. This study tailors application of the so-called Gradient Projection for Sparse Reconstruction (GPSR) to large-scale time-difference three-dimensional electrical impedance tomography (3D EIT). 3D EIT typically suffers from the need for a large number of voxels to cover the whole domain, so its application to real-time imaging, for example monitoring of lung function, remains scarce since the large number of degrees of freedom of the problem extremely increases storage space and reconstruction time. This study shows the great potential of the GPSR for large-size time-difference 3D EIT. Further studies are needed to improve its accuracy for imaging small-size anomalies.
Keywords Three-dimensional electrical impedance tomography Sparse recovery Gradient Projection for Sparse Reconstruction Lung imaging
* Manuchehr Soleimani [email protected]
1 Institute of Technical Medicine, Faculty of Medical and Life
Sciences, Furtwangen University of Applied Sciences, Villingen-Schwenningen, Germany
2 Engineering Tomography Laboratory (ETL), Department of Electrical Engineering, University of Bath, Bath, UK
1 3
1244 Med Biol Eng Comput (2016) 54:12431255
To cope with such errors, time-difference reconstruction is often given precedence over absolute reconstruction. Employing the time-difference reconstruction, the objective is to infer conductivity changes from difference between two boundary datasets that are measured at different times [5].
Classical quadratic inverse solvers in EIT often consider some correlations between adjacent nite elements, thereby reducing the ill-posedness of the problem [1]. In this way, the problem is stabilized at the cost of imposing some smoothness on the reconstructed image, so detecting sharp discontinuities over the conductivity eld will be impossible [11, 36]. There are, however, many organs that have well-dened boundaries, and thus represent sharp variations over the conductivity prole, e.g., interfaces between collapsed and ventilated regions of lung [10, 11].
To overcome this effect, total variation (TV) regularization has been applied to EIT, thanks to its ability to better preserve sharp interfaces, and compared to the classical quadratic regularization [11, 17]. To the best of our knowledge, the TV reconstruction schemes that have been applied to EIT so far are based on Newtons method, e.g., primaldual interior point method or Lagged diffusivity method [11, 17]. These codes are available on EIDORS website [19]. Newtons methods intuitively require inverse Hessian, so their application to 3D EIT leads to very large computations.
EIT is inherently three-dimensional since electrical current cannot be conned to ow solely at the electrodes plane. As a result, 2D EIT is subject to artifacts produced by contrasts that are off the electrodes plane [27]. 3D EIT has thus received much attention with at least two planes of electrodes [9, 20, 24, 25, 27, 29, 33, 35, 42, 51, 55]. Among quadratic regularized solvers, Krylov subspace methods such as conjugate gradient (CG) best suit 3D EIT, as classical Newtons methods involve the costly inversion of Hessian [33, 35, 55]. Indeed, the very large size of forward operator in 3D EIT increases the ill-posedness of the problem, as well as computational time. As a result, the main advantage of EIT over other imaging modalities for real-time imaging will be lost [9, 20].
Finding sparse solutions to large-size linear systems of equations has attracted much interest recently [8, 21, 39]. The presence of an 1 norm as the regularization function encourages small components of the unknown parameters to become exactly zero, thus promoting sparse solutions [21].
To reduce the computational cost of calculating the conductivity over a large number of nite elements in 3D EIT, this study tailors the application of a class of sparsity
inverse solvers that does not require the Jacobian (J) to be stored explicitly, but only needs matrixvector products including J and JT [8, 21, 28, 39]. This study shows the great potential of the so-called Gradient Projection for Sparse Reconstruction (GPSR) method for time-difference 3D EIT. The GPSR was rst proposed in the context of signal processing [21] and has then been extended to other applications (See [12, 16, 31, 34, 44, 48, 54, 57]). Applying time-difference imaging, the GPSR splits the update at each iterate into its negative and positive parts and then enforces a nonnegativity constraint on each part so that the background conductivity is gradually set to zero with the progress of the solver [21]. Since for the time-difference reconstruction, the background is typically set to zero, and a solution around zero is expected, it benets notably from the splitting behavior of the GPSR, unlike other competing sparsity solvers, which fail to possess this advantage.
Furthermore, the superiority of the GPSR over other popular sparsity solvers, e.g., l1-ls [39], two-step iterative shrinkage thresholding (TWIST) [8], or xed-point continuation (FPC) [28], for large-scale inverse problems has already been demonstrated, e.g. [21].
The sparsity regularization for EIT has attracted much interest recently. To best of our knowledge, the most well-known algorithm for sparsity reconstruction in EIT was proposed in [38] and was then applied to real experiments [22, 23]. Although similar to the GPSR, this algorithm follows a gradient-based method which requires solely matrixvector products, it does not benet from the splitting scheme in the GPSR. In addition, a direct application of the gradient of the residual leads to numerical instability for this algorithm, even in two-dimensional cases, so a Sobolev smoothing step is applied to the gradient via solving an augmented Dirichlet boundary value problem at each iterate, which increases the computational cost [22, 23, 38].
Typically, the gradient-based sparsity reconstruction considerably reduces the computational cost in comparison with Newton-based 1 reconstructions such as [41]. In spite of the very fast nature of these solvers, which is highly demanded for 3D imaging of lung function, their application to 3D lung imaging has not been reported so far. To the best of our knowledge, this manuscript reports the rst application of the gradient-based sparsity reconstruction for this case.
Preconditioned conjugate gradient (PCG), the most popular algorithm for large-size 3D EIT, available on the EIDORS website [19], is considered as the rst benchmark [33, 35, 55]. The second benchmark is the sparsity algorithm specied for EIT in [22, 23, 38].
1 3
Med Biol Eng Comput (2016) 54:12431255
1245
2 Method
2.1 Forward and inverse models
The forward problem in EIT is to calculate the surface electrical potentials on the electrodes (U) as a function of the injected currents (I) and the conductivity distribution (~). To implement the forward problem, elliptic partial differential equations (PDEs) are dened over the mesh according to Ohms law. Neumann and Dirichlet boundary conditions are determined as functions of the boundary datasets (I, U). The resulting nonlinear systems of equations are written as U = F(~)I, where F(~) : I u denotes Neumann-to-Dirichlet (NtD) map, and : u U represents Dirichlet-to-observation map [15, 32]. The resulting problem is nonlinear with respect to the conductivity, so it is linearized around the background conductivity ~bg via computing the Jacobian (J) as follows [40].
Applying the time-difference imaging, the objective is to calculate conductivity changes ~~ from difference between two data frames that are measured at times t1 and t2(~V) [5, 47].
In a matrix notation, in light of = J~~, the inverse problem is to infer ~~ from the real difference measured data ~V in the form
The unconstrained Tikhonov form of (2) can be written as
where
2.2 Preconditioned conjugate gradient (PCG) inverse solver
The choice of r = 2 turns problem (3) to quadratic regularized form [1, 14, 40, 50, 56]. In this work, to evade the costly computation of inverse Hessian, the quadratic optimization problem is solved by the PCG, rather than Newtons methods. PCG method is a Krylov subspace techniques [46]. Generally, setting the derivative of least squares problem min~~ ~V J~~ 22 to zero yields
Applying the PCG to Eq. (5) yields the following optimization algorithm. Given J, ~V, initial guess ~~0, preconditioner M, and stopping tolerance , the algorithm is outlined as follows.
(1)
Algorithm 1. Preconditioned Conjugate Gradient [47]
0
~
~
0
r
0 ~
~ ~
J
V
0
d
0 r
M
1
0
0 )
( d
r
0
T
0
While )
/
(
k Do
0
q
k Jd
k
1
k
(
d )
k
k
T
q
k
k
U = (F(~) F(~bg))I J(~ ~bg).
~
~
k
k d
~
~
1
k
r
k
k q
r
1
k
k
s
k r
M
1
k
k
k s
r )
(
1
T
k
min
~~ J~~ ~V 22 s.t. a priori assumption on ~~
(2)
k
k
k
d
k
k d
s
1
k
k
min
~~ J~~ ~V 22 + [notdef]Rr(~~),
(3)
k
End Do
where dk and k denote the search direction and step length at iteration k, respectively. In this work, the PCG is implemented with the aid of the EIDORS software [19] and is regarded as the rst benchmark for evaluating the performance of the proposed inverse solver. The second benchmark is the most well-known sparsity algorithm in EIT [22, 23, 38], which is outlined in Appendix. For further theoretical details, the reader is referred to [22, 23, 38].
2.3 Sparse recovery
Many different approaches have been proposed to seek a sparse solution to large-size linear system y = Ax + n, where y is the observation data and n is the noise. Roughly they can be divided into two categories, i.e., constrained and unconstrained optimization problems. Assuming T and
k
1
Rr() = 1
r rr.
(4)
JTJ~~ = JT~V.
(5)
1 3
1246 Med Biol Eng Comput (2016) 54:12431255
to be nonnegative parameters, the constrained form leads to the following two formulations, i.e.,
the so-called quadratic program (QP), i.e., least absolute shrinkage and selection operator (LASSO) [49], and
namely quadratically constrained linear program (QCLP),
or Basis pursuit with = 0 [13].
There has also been much interest in solving the unconstrained form of the problem, i.e.,
where [notdef] is the regularization parameter. It was proved that a solution of (6) for T 0 is a minimizer of (8) for some
[notdef] > 0. Similarly, a solution of (7) is either X = 0, or a minimizer of (8) for some [notdef] > 0 (31).
2.4 Gradient Projection for Sparse Reconstruction (GPSR) inverse solver
The choice of r = 1 conducts problem (3) to
which is equivalent to Eq. (8).
The GPSR approach is applied in this work to infer a sparse solution of the conductivity changes from the difference data V. The base of this approach is to initially split the unknown vector ~~ into its positive and negative parts and then enforcing a nonnegativity constraint on each part [21], i.e.,
Accordingly, considering a mesh made up of n nite elements, for i = 1, 2, . . . , n,
where (x)+ = max{0, x}. Considering the penalty term in (9) in the form ~~ 1 = 1Tnu + 1Tnw with 1n = [1, 1, . . . , 1]T
yields [21]
Supposing z =
min
z cTz +
1
2zTBz (z) s.t. z 0.
(13)
Now gradient projection (GP) method is employed, which involves two stages at each iteration [7, 21]. First, given zk , k > 0 is chosen as a step length for searching along negative direction (zk) from zk in the feasible set.
(7) The nonnegativity constraint imposed to ~k iteratively nulls the background with the progress of the algorithm. k [0, 1] is then chosen to set
The step length k is chosen in two different ways, which is explained in the sequel.
2.4.1 Basic GPSR
Employing the basic variant of the GPSR, the gradient is dened as [21]
Applying the gradient in this way prevents the elements of z that were nulled by the constraint in the previous iterates from any further updates, reducing the number of degrees of freedom of the problem iteratively. An initial guess for is chosen so that is minimized along gk, i.e.,
An exact minimizer for the above problem is written as [21]
min
x y Ax 22 s.t. x 1 T,
(6)
min
x x 1 s.t. y Ax 22 ,
~k = (zk k(zk))+
(14)
zk+1 = zk + k(~k zk).
(15)
min
x y Ax 22 + [notdef] x 1,
(8)
gki =
[notdef]
(zk)i, if zki > 0 or ((zk))i < 0 0, otherwise ,
(16)
min
~~ ~V J~~ 2 + [notdef] ~~ 1,
(9)
k0 = arg min
(zk gk).
(17)
~~ = u w, u 0, w 0.
(10)
k0 =
(gk)Tgk
(gk)TBgk ,
(18)
ui = (~~i)+
wi = (~~i)+,
(11)
0 is then encountered by upper and lower bounds min and max. For the backtracking line search, considering scalar parameters (0, 1) and (0, 1/2), k is chosen to be the rst value in the sequence 0, 0, 20, . . . that satises
~((zk k~(zk))+) ~(zk) ~(zk)T(zk
(zk k~(zk))+)
(19)
min
u,w ~V J(u w) 22 + [notdef]1Tn u + [notdef]1Tn w s.t. u 0, w 0.
[notdef]
(12)
u w
, c = [notdef]12n +[notdef] JT~V JT~V
and
2.4.2 BarzilaiBorwein GPSR
Employing BarzilaiBorwein (BB) scheme,k is chosen such that kI approaches the inverse Hessian 2(x) over the latest step. Letting sk = zk zk1 and k = (zk) (zk1), the exact step length is computed as [6, 21]
B =
[notdef]
JTJ JTJ
JTJ JTJ
, (12) is rewritten as
1 3
Med Biol Eng Comput (2016) 54:12431255
1247
2.5 Analysis of computational cost
At rst glance, the problem (13) may appear more costly than the classical form of the problem in (12) since the dimension of the problem becomes twice, i.e., z R2n . However, considering (z) = cTz + 12zTBz, matrix
c = [notdef]12n + [notdef]
k0 = arg min
[notdef][notdef][notdef]
[notdef][notdef][notdef]
sk k
2 (sk)Tk
(sk)Tsk .
(20)
0 is then encountered by upper and lower bounds min and max, similar to the basic GPSR. In comparison with the BB GPSR in [21], a more sophisticated version of the
BB rule was employed in our study. This BB rule was similarly applied to the standard sparsity algorithm presented in Appendix, according to [22, 23, 38]. The step length computed by Eq. (20) is reduced in an inner iteration by a constant factor until the following criterion in Eq. (21) is satised. It turns out that enforcing monotonicity would deteriorate the behavior of the BB rule on the convergence. As a result, a globally convergent BarzilaiBorwein is employed, where zk+1 is accepted as a new iterate if [52]
where i denotes the Q previous iterations, and (0, 1) is a constant that is often chosen near zero. The stopping criterion is dened based on perturbation results for linear complementarity problems (LCP) as follows [21, 53]. There exists a scalar parameter such that
where s represents the solution to problem (13), dist() denotes the distance operator, and the minimizing operator is taken component-wise. In light of (22), the algorithm is terminated if
The GPSR algorithm is outlined as follows.
Algorithm 2. Gradient Projection for Sparse Reconstruction [31] Set 0
0
z , and choose max
min
JT~V
JT~V
is independent from the unknown
parameter ~~, so it can be computed once at the start of the algorithm. In addition, zTBz yields a scalar value, which can be easily calculated as
Therefore, the total cost for the calculation of (z) at each iteration involves a matrixvector product J(u w) and a vectorvector product cTz. To minimize in problem (13), the gradient of (z) is computed at each iterate in the form
Considering
the cost for computing (z) is two matrixvector products, as c is calculated independently from the iterates. In this study, the application of the GPSR to 3D EIT was tailored. To implement the GPSR, a MATLAB code in the context of signal processing, which is available online [26], was modied.
3 Numerical results
A 3D shape of an adult human thorax was simulated by using the EIDORS and NETGEN software as follows [19, 43]. The contour of a human thorax and lungs were plotted according to a CT image available on the EIDORS and were then mapped onto a 3D nite element mesh generated by the NETGEN software [43]. The created 3D mesh was made up of 161,021 tetrahedral elements with a height of 1. Thirty-two circular electrodes were installed around the chest in two rings aligned by axial planes 0.33 and 0.66. The electrodes were simulated based on complete electrode model with a contact impedance of 100 and a radius of 0.05 [51]. An electrical current with amplitude of 5 mA was successively injected through the electrodes, and the induced potentials were measured according to planar alignment protocol, the most well-known protocol for 3D EIT [27].
zTBz = (u w)TJTJ(u w) = J(u w) 22.
(24)
max
(k+Q+1ik (Zi) k
~zk+1
zk+1 zk
2 (21)
(z) = c + Bz.
(25)
Bz =
[notdef]
JTJ(u w)
JTJ(u w)
,
(26)
dist(z, s) min(z, (z) ,
(22)
min(z, (z) tol.
(23)
0 ,While stopping criterion is not satisfied Do
Choose based on basic (2.4.1) or BB (2.4.2) scheme
Set
k z
~ (
Calculate Set
k
z
k
k ~
z
1
k
k
End Do
k
1
1 3
1248 Med Biol Eng Comput (2016) 54:12431255
Fig. 1 The simulated human chest from: a a 3D view, b a top view, and c a lateral view
The noise contributed to the measurement data was an additive white Gaussian noise (AWGN). Considering the difference imaging, AWGN is simulated as
where NL is noise level, std is standard deviation of difference between two frames of data, and randn is a vector denoting pseudorandom values drawn over a standard Gaussian distribution. The simulated data were contaminated with a 20 db AWGN, i.e., NL = 1020/20 = 0.1 (see [19]).
To avoid the so-called inverse crime, the inverse solver was applied to a coarser mesh made up of 20,955 elements. The background conductivity was set to 1, while the lungs conductivity was set to 0.3 S m1. Figure 1a shows the simulated chest phantom from a 3D view. Figure 1b, c, respectively, exhibits the mesh from a top and a lateral view. The time-difference reconstruction was employed with the background conductivity as the reference data.
The PCG solver was implemented and terminated at threshold = 1e 2. This parameter was heuristically selected to produce the image best tting the simulated model over 100 iterates. The sparsity algorithm specied for EIT in [22, 23, 38], referred to as standard sparsity algorithm here, the basic GPSR and BB GPSR was implemented and was continued until the threshold tol = 1e 2. Heuristically, among a wide range of regularization parameters, [notdef] = 1e 4 produced the optimal image for both the basic and BB GPSR. The optimal image for the standard sparsity algorithm was produced by [notdef] = 5e 7.
Figure 2 displays the 3D images reconstructed by the solvers at equidistant axial cross sections. From the left side, the rst column corresponds to the 3D image reconstructed by the PCG solver, the second column represents the standard sparsity algorithm, and the third and fourth columns pertain to the basic GPSR, and BB GPSR, respectively. The slices were exhibited according to the
(27)
Noise = NL std(V) randn,
Fig. 2 From the left side, the 3D images reconstructed by the PCG, standard sparsity, basic GPSR, and BB GPSR, displayed at transverse planes written to the left of the gure
color bar shown to the right of each column, which was adjusted such that its minimal value represents discrepancy of the conductivity between the lungs and background in the simulated phantom, i.e., 0.7 S m1. The
vertical position of each transverse plane was written to the left of the gure.
3.1 Observations
As shown in the left column of Fig. 2, the PCG produced a great blurriness and failed to precisely detect the sharp conductivity jumps over the reconstructed image. These
1 3
Med Biol Eng Comput (2016) 54:12431255
1249
Fig. 3 The images reconstructed by a the standard sparsity (second column in Fig. 2), b basic GPSR (third column in Fig. 2), and c BB GPSR (fourth column in Fig. 2), from a 3D view. To better visualize the inclusions, maximal values of color brs were increased, compared to Fig. 2
Table 1 The RE of the images reconstructed from the simulated chest and the CPU time
PCG Standard sparse Basic GPSR BB GPSR
RE 0.61 0.35 0.24 0.20
Time 25.81 10.73 1.93 5.79
solution ~phantom
RE =
[notdef][notdef]~
[notdef][notdef]~
[notdef][notdef]
[notdef][notdef]
2
,
(28)
phantom
2
smoothing effects have given rise to a low spatial resolution aligning all dimensions. A fair comparison between Fig. 1b and the second column in Fig. 2 reveals that the standard sparsity algorithm was not tolerant enough to accurately determine the lung boundaries. Furthermore, comparing the second column to Fig. 1c shows that the spatial resolution aligning the longitudinal axis has been lost, and a great artifact has been produced in the top and bottom slices, where the lungs do not exist.
As displayed in the third and fourth columns of Fig. 2, the GPSR solver considerably improved the reconstruction in determining the sharp jumps over the conductivity prole. A comparison between Fig. 1 and the two last columns reveals that the GPSR determined the lung boundaries more accurately than the PCG and standard sparsity algorithms over all the transverse planes. According to the color bars, the GPSR determined the conductivity changes amplitude with a much higher accuracy as well. Furthermore, compared to the basic GPSR, conducting the GPSR through the BB scheme improved the solution regarding contrast, as well as the amount of artifact. Figure 3ac, respectively, exhibit the images reconstructed by the standard sparsity, basic GPSR and BB GPSR from a 3D view. In other words, Fig. 3ac represents 3D views of the images shown in the second, third and fourth columns in Fig. 2. The image reconstructed by the PCG was neglected, as it was covered by very large amount of artifacts around the lungs.
The performance of the solvers was evaluated in terms of Relative Error (RE), i.e.,
where ~phantom denotes the conductivity distribution over the simulated phantom, and ~solution denotes the computed absolute conductivity, i.e., the calculated conductivity changes plus the background. Accordingly, the rst row in Table 1 shows the RE for the reconstructed images. The RE of the images reconstructed by the basic and BB GPSR was, respectively, 31 and 43 % smaller than that of the standard sparsity algorithm. The PCG produced much greater RE than all the employed sparsity algorithms.
3.2 Computational cost
The second row in Table 1 shows the CPU time elapsed on implementing the solvers. All the computed CPU times involve 0.81 s for computation of the Jacobian. The processor that was employed in this work is an Intel Core i3-3220 Processor (3.30 GHz) with a RAM of 4 GB and a 64-bit operating system (Windows 7, Microsoft, Seattle, WA).
According to Table 1, the GPSR considerably reduced the computational time and compared to the competing algorithms. Since the BB scheme typically conducts the solution without forcing the objective function to decrease monotonically through all the iterations, the CPU time elapsed on the BB GPSR was three times more than the basic GPSR. In this way, the BB scheme reconstructed the image with an RE 17 % smaller than the basic GPSR.
4 Experimental results
To validate the proposed image reconstruction, the algorithms were also tested on a real data measured from a
1 3
1250 Med Biol Eng Comput (2016) 54:12431255
Fig. 4 The chest phantom created for reconstructing image from the experimental data pertaining to thirty-four frames of a breathing cycle
human chest. The data, which is available on the EIDORS website, pertain to thirty-four frames of a breathing cycle of a human subject [18]. To recover the real shape of lungs in a 3D representation, a very ne mesh was needed. Accordingly, the inverse problem was applied to a simulated chest made up of 143,119 voxels. Figure 4 displays the created mesh from a 3D view. The electrodes were placed aligning axial plane 0.5 and were represented by the green circles.
Employing the solvers, the time-difference imaging was applied such that the rst frame was used as the reference data. The PCG algorithm was rst employed, and the images were calculated at stopping threshold = 7e 3 , which produced the best image heuristically. Figure 5 shows the reconstructed images concerning all the thirty-four frames. Note that these frames are originally 3D images that are shown solely at axial plane 0.5, aligned by the electrodes plane, due to space constraints.
Since the simulated results showed the superiority of the GPSR over the standard sparsity algorithm proposed in [22, 23, 38] regarding both accuracy and time, the standard sparsity solver was neglected in this section. The GPSR solvers were implemented until tol = 1e 2, and 3D images were computed. The regularization parameter was heuristically chosen to be [notdef] = 1e 2. Figure 6a, b, respectively, shows the images reconstructed by the basic and BB GPSR at the electrodes plane.
4.1 Observations
The results show that PCG failed to accurately recover shape of the lungs over a 3D chest phantom with such large number of voxels. The images include great amounts of artifact as well. On the other hand, the GPSR considerably improved the reconstruction. Both the basic and BB GPSR better determine the lung shape during the breathing cycle, compared to the PCG. In addition, a fair comparison between Fig. 6a and b shows that the BB scheme provided
a better contrast over the frames than the basic scheme, as well as a smaller amount of artifact.
For all the solvers, the maximal mean of conductivity changes occurred for the 22nd frame of the breathing cycle. Figure 7 exhibits the images pertaining to the 14th, 22nd, and 30th frames from a 3D scene. Figure 7ac shows the images of the 14th frame, in the middle of inhalation stage, reconstructed by the PCG, basic and BB GPSR, respectively. Similarly, Fig. 7df shows 3D views of the images of the 22nd frame, maximum end-respiratory. Figure 7gi, respectively, represent the 3D images reconstructed by the PCG, basic and BB GPSR for the 30th frame of the breathing cycle, in the middle of exhalation stage. In other words, Fig. 7a, d, g represents the three aforementioned frames of Figs. 5, 7b, e, h show the 3D view of these frames in Figs. 6a, 7c, f, i represent the corresponding frames of Fig. 6b. To better visualize the 3D images, maximal values of color bars were increased, compared to Figs. 5 and 6 so that very small artifacts were neglected. As shown in this gure, compared to the PCG, the GPSR solvers better recovered shape of the lungs and reduced the amount of artifact.
4.2 Computational cost
In real-time imaging, e.g., monitoring of pulmonary function, the reconstruction time is vital, which limits the applicability of the 3D reconstruction as a result of the need for a large number of voxels. Table 2 shows the CPU time consumed by the employed solvers for reconstructing 3D images of all the thirty-four frames. According to this table, the GPSR noticeably reduced the reconstruction
Fig. 5 The 3D images reconstructed by the PCG solver at thirty-four frames of the breathing cycle. The slices were taken at the electrodes plane
1 3
Med Biol Eng Comput (2016) 54:12431255
1251
Fig. 6 The 3D images reconstructed by a the basic GPSR and b BB GPSR at thirty-four frames of the breathing cycle. The slices were taken at the electrodes plane
time, compared to the PCG. This table also afrms that the execution of the BB GPSR is more than the basic GPSR. Indeed, the BB scheme does not force the objective function to decrease monotonically through all iterates and adjusts the step length in a more sophisticated way, so the accuracy of the solution was improved at the cost of the reconstruction time. The results show that the BB and basic GPSR were, respectively, 4.83 and 9.09 times faster than the PCG in reconstructing all the frames.
5 Discussion
The results show that both the basic and BB GPSR are more tolerant than the standard sparsity and PCG algorithms for 3D EIT of human chest. To the best of our knowledge, the most valid sparsity solver for EIT is the algorithm proposed in [38], which was then tested on real data in [22, 23]. Namely standard sparsity here, we regard this algorithm as the benchmark for testing the performance of the GPSR inverse solvers. According to Appendix, the convergence of this algorithm is contingent on computation of a smoothed Sobolev gradient of the residual norm, which imposes some smoothness on the solution, as well as extra time for solving the corresponding Dirichlet boundary value problem at each iterate. As shown in the simulated results in Sect. 3, the basic and BB GPSR outperformed the standard sparsity solver regarding the accuracy and time.
In addition, the PCG is regarded as the most well-known solver for 3D EIT in both literature and EIDORS website
since it does not require the inverse Hessian, in contrast to classical Newtons methods [19, 33, 35, 55]. The PCG, however, deleteriously smoothed the solution and failed to accurately determine sharp conductivity jumps over lung boundary. As a result, the resulting images contain very little information for diagnostic purposes. In addition, although its computational cost is much lower than Newtons method, it is not yet tolerant enough to deal with a large number of nite elements, i.e., more than 20,000 voxels. Indeed, the 3D reconstruction for real cases suffers severely from the large number of degrees of freedom of the problem, which increases the ill-posedness, as well as the reconstruction time.
The GPSR was already proposed in the context of signal processing for sparse recovery of signals [21]. A modied version of the GPSR was tailored in this study for the time-difference 3D EIT. Both the numerical and experimental results reveal that the basic and BB GPSR noticeably improved 3D EIT for real-time lung imaging, compared to the PCG and the standard sparsity solvers. Applying the GPSR to the time-difference 3D EIT, the updated conductivity changes with positive and negative values are regarded as separate sparse vectors and subsequently enforcing a nonnegativity constraint to the gradient projection nulls the background conductivity through the iterates. The results show the high potential of the GPSR for time-difference 3D EIT. The computations only require matrix vector products, so the computational cost arising from explicitly storing J and JTJ will be removed. Both the basic and BB GPSR provided a more accurate conductivity
1 3
1252 Med Biol Eng Comput (2016) 54:12431255
Fig. 7 The images of the 14th frame of the breathing cycle from a 3D view, reconstructed by: a PCG, b basic GPSR, and c BB GPSR, the images of the 22nd frame of the breathing cycle from a 3D view,
reconstructed by: d PCG, e basic GPSR, and f BB GPSR, and the 3D images of the 30th frame, reconstructed by: g PCG, h basic GPSR, and i BB GPSR
Table 2 The CPU time elapsed on reconstructing images from the real human lung data over all the thirty-four frames (s)
PCG Basic GPSR BB GPSR
144.79 15.92 29.94
prole during the breathing cycle regarding determination of interfaces, as well as the conductivity amplitude. In addition, the results show that the GPSR better addressed the large number of voxels and appreciably reduced the CPU time against the competing solvers.
However, according to our previous numerical results, the main drawback of the GPSR for application to 3D EIT is that its accuracy is severely deteriorated in determining
very small inclusions [37]. This problem was addressed by adopting a compressive sensing scheme for sampling the nite elements covering the inclusions. By application of a preprocessing PCG step, the proposed scheme improved accuracy at the cost of speed. For further information, the reader is referred to [37].
6 Conclusion
3D EIT typically suffers from the need for a large number of nite elements to cover the whole domain, thus requiring very large computations [33, 35, 55]. Although the main advantage of 2D EIT for real-time imaging of lung is its high speed [3, 4], the applicability of 3D EIT to this
1 3
Med Biol Eng Comput (2016) 54:12431255
1253
case remains scarce, since the large number of degrees of freedom of the problem increases the computational cost notably. Indeed, recovery of inter-medium interfaces inside human organs will be erroneous over a coarse mesh. For example, for lung monitoring, sharp jumps over the conductivity prole cannot be detected suitably over a coarse mesh, thus providing misleading physiological information. This study showed that the GPSR algorithm best suits 3D EIT of chest. The GPSR suitably deals with a large
number of voxels and determined the conductivity eld more accurately than the PCG and standard sparsity solvers, at the same time considerably reducing the computational time. Further studies are still needed to cope with its poor performance in imaging small-size anomalies.
Acknowledgments This work was partially supported by the Federal Ministry of Education and Research (BMBF) in Germany under Grant No. 03FH038I3 (MOSES).
Appendix
Algorithm 3. The sparsity algorithm proposed for EIT in [43-45] (standard sparsity solver)
Set
While stopping criterion is not satisfied Do
Compute
Compute the Sobolev smoothed gradient via Dirchlet boundary value problem
s.t.
Deter r iteration based on the BB scheme Update the conductivity changes
End DO
where and represent the medium and its boundary, respectively, and Ds denotes Sobolev smoothed gradient of the residual (see [38]). S denotes the soft shrinkage operator with k, which sets small elements of each update to zero, thus promoting the sparse solution, denotes the Laplacian operator, and is a scalar parameter controlling the degree of smoothing of the gradient.
Similar to the GPSR algorithm, the stopping criterion was adopted based on LCP, according to Eqs. (22) and (23) with replacing zk by ~k.
The rationale behind employing the smoothed Sobolev gradient is that heuristically, the direct application of to this algorithm exhibits unappealing oscillations in the reconstruction, which often gives rise to numerical
instability. This is another superiority of the GPSR over this algorithm, as the GPSR suitably converges through the direct application of the gradient.
Applying the BB scheme, an initial guess for the step length is made according to Eq. (20). This initial guess is then iteratively reduced until the following equation is satised
where, i denotes the Q previous iterates, and (0, 1) is a constant that is often chosen near zero.
~ ~~0 + S
~~k kDs
~~k max (kQ+1ik)
~
~~i
k 2
S
~~k kDs
~~k ~k
2 (29)
1 3
1254 Med Biol Eng Comput (2016) 54:12431255
References
1. Adler A, Guardo R (1996) Electrical impedance tomography: regularized imaging and contrast detection. IEEE Trans Biomed Eng 15:170179
2. Adler A, Arnold JH, Bayford R, Borsic A, Brown B, Dixon P, Faes TJC, Frerichs I, Gagnon H, Grber Y, Grychtol B, Hahn G, Lion-heart WRB, Malik A, Patterson RP, Stocks J, Tizzard A, Weiler N, Wolf GK (2009) GREIT: a unied approach to 2D linear EIT reconstruction of lung images. Physiol Meas 30(6):3555
3. Adler A, Amato MB, Arnold JH, Bayford R, Bodenstein M, Bhm SH, Brown BH, Frerichs I, Stenqvist O, Weiler N, Wolf GK (2012) Whither Lung EIT: where are we, where do we want to go and what do we need to get there? Physiol Meas 33(5):679694
4. Bagshaw AP, Liston AD, Bayford RH, Tizzard A, Gibson AP, Tidswell T, Sparkes MK, Dehghani H, Binnie CD, Holder DS (2003) Electrical impedance tomography of human brain function using reconstruction algorithms based on the nite element method. Neuro Image 20:752764
5. Barber DC, Brown BH (1988) Errors in reconstruction of resistivity images using a linear reconstruction technique. Clin Phys Physiol Meas 9:101104
6. Barzilai J, Borwein J (1988) Two point step size gradient methods. IMA J Numer Anal 8:141148
7. Bertsekas DP (1999) Nonlinear programming, 2nd edn. Athena, Boston
8. Bioucas-Dias J, Figueiredo M (2007) A new TwIST: two-step iterative shrinkage/thresholding algorithms for image restoration. IEEE Trans Image Process 16(12):29923004
9. Blue RS, Isaacson D, Newell JC (2000) Real-time three-dimensional electrical impedance imaging. Physiol Meas 21:11210. Borsic A, Lionheart WRB, McLeod CN (2002) Generation of anisotropic-smoothness regularization lters for EIT. IEEE Trans Med Imag 21(6):579587
11. Borsic A, Graham BM, Adler A, Lionheart WRB (2010) In vivo impedance imaging with total variation regularization. IEEE Trans Med Imag 29(1):4454
12. Brunelli D, Caione C (2015) Sparse recovery optimization in wireless sensor networks with a sub-Nyquist sampling rate. Sensors 15:1665416673
13. Chen S, Donoho D, Saunders M (1998) Atomic decomposition by basis pursuit. SIAM J Sci Comput 20:3361
14. Cheney M, Isaacson D, Newell JC, Simske S, Goble J (1990) NOSER: an algorithm for solving the inverse conductivity problem. Int J Imaging Syst Technol 2:6675
15. Cheney M, Isaacson D, Newell JC (1999) Electrical impedance tomography. SIAM Rev 41:85101
16. Cuadros AP, Arce GR, Arguello H (2014) Coded aperture design in compressive X-ray tomography. In: IEEE global conference on signal and information processing, 35 Dec 2014
17. Dobson DC, Santosa F (1994) An image enhancement technique for electrical impedance tomography. Inverse Probl 10:317334
18. Electrical Impedance Tomography and Diffuse Optical Tomography Reconstruction Software (EIDORS). http://eidors3d.sourceforge.net/tutorial/lung_EIT/tutorial310-lung-images.shtml
Web End =http://eidors3d.source- http://eidors3d.sourceforge.net/tutorial/lung_EIT/tutorial310-lung-images.shtml
Web End =forge.net/tutorial/lung_EIT/tutorial310-lung-images.shtml
19. Electrical Impedance Tomography and Diffuse Optical Tomography Reconstruction Software (EIDORS), Released version: EIDORS 3.7.1 (29 May 2013), http://eidors3d.sourceforge.net/
Web End =http://eidors3d.sourceforge.net/
20. Fan WR, Wang HX (2010) 3D modeling of the human thorax for ventilation distribution measured through electrical impedance tomography. Meas Sci Technol 21:115801
21. Figueiredo M, Nowak R, Wright S (2007) Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems. IEEE J Sel Top Signal Process 1:586598
22. Gehre M, Kluth T, Lipponen A, Jin B, Seppnen A, Kaipio JP, Maass P (2012) Sparsity reconstruction in electrical impedance tomography: an experimental evaluation. J Comput Appl Math 236:21262136
23. Gehre M, Kluth T, Sebu C, Maass P (2014) Sparse 3D reconstructions in electrical impedance tomography using real data. Inverse Probl Sci Eng 22(1):3144
24. Gobel JC, Cheney M, Isaacson D (1992) Electrical impedance tomography in three dimensions. Appl Comput Electromagn Soc J 7:128147
25. Goharian M, Soleimani M, Moran G (2009) A trust region subproblem for 3D electrical impedance tomography inverse problem using experimental data. Prog Electromagn Res 94:1932
26. Gradient Projection for Sparse Reconstruction (GPSR v6.0). http://www.lx.it.pt/%7emtf/GPSR/
Web End =www.lx.it.pt/~mtf/GPSR/
27. Graham BM, Adler A (2007) Electrode placement congurations for 3D EIT. Physiol Meas 28:2944
28. Hale T, Yin W, Zhang Y (2007) A xed-point continuation method for l1-regularized minimization with applications to compressed sensing. Department of Computational and Applied Mathematics,
Rice University, Houston, TX, Technical Report TR07-0729. Halter RJ, Hartov A, Paulsen KD (2007) Experimental justication for using 3D conductivity reconstructions in electrical impedance tomography. Physiol Meas 28:115127
30. Halter RJ, Hartov A, Paulsen KD (2008) A broadband high-frequency electrical impedance tomography system for breast imaging. IEEE Trans Biomed Eng 55(2):650659
31. He X, Liang J, Wang X, Yu J, Qu X, Wang X, Hou Y, Chen D, Liu F, Tian J (2010) Sparse reconstruction for quantitative bioluminescence tomography based on the incomplete variables truncated conjugate gradient method. Opt Express 18(24):2482524841
32. Holder DS (2004) Electrical impedance tomography: methods, history and applications. Institute of Physics Publishing, pp 364
33. Horesh L, Bolhofer M, Schweiger M, Arridge SR, Holder DS (2007) Novel large-scale 3D electrical impedance tomography modeling of the human head. IFMBE Proc 14:38583861
34. Howland GA, Lum DJ, Howell JC (2014) Compressive wavefront sensing with weak values. Opt Express 22(16):1887018880
35. Javaherian A, Soleimani M (2013) Compressed sampling for boundary measurements in three-dimensional electrical impedance tomography. Physiol Meas 34:11331150
36. Javaherian A, Movafeghi A, Faghihi R (2013) Reducing negative effects of quadratic norm regularization on image reconstruction in electrical impedance tomography. Appl Math Model 37(8):56375652
37. Javaherian A, Soleimani M, Moeller K (2015) Sampling of nite elements for sparse recovery in large scale 3D electrical impedance tomography. Physiol Meas 36(1):4366
38. Jin B, Khan T, Maass P (2012) A reconstruction algorithm for electrical impedance tomography based on sparsity regularization. Int J Numer Methods Eng 89:337353
39. Kim SJ, Koh K, Lustig M, Boyd S, Gorinevsky D (2007) An interior-point method for large-scale l1-regularized least squares.
IEEE J Sel Top Signal Process 1(4):60661740. Lionheart WRB (2004) EIT reconstruction algorithms: pitfalls, challenges and recent developments. Physiol Meas 25:125142
41. Mamatjan Y, Borsic A, Gursoy D, Adler Andy (2013) An experimental clinical evaluation of EIT imaging with l1 data and image norms. Physiol Meas 34:10271039
42. Metherall P, Barber DC, Smallwood RH, Brown BH (1996) Three dimensional electrical impedance tomography. Nature 380:509512
43. Netgen Mesh Generator (v5.0). http://sourceforge.net/projects/netgen-mesher/files/netgen-mesher/5.0/
Web End =http://sourceforge.net/projects/ http://sourceforge.net/projects/netgen-mesher/files/netgen-mesher/5.0/
Web End =netgen-mesher/les/netgen-mesher/5.0/
44. Park JC, Song B, Kim JS, Park SH, Kim HK, Liu Z, Suh TS, Song WY (2012) Fast compressed sensing-based CBCT
1 3
Med Biol Eng Comput (2016) 54:12431255
1255
reconstruction using BarzilaiBorwein formulation for application to on-line IGRT. Med Phys 39(3):1207121745. Rezajoo S, Hossein-Zadeh G (2010) Reconstruction convergence and speed enhancement in electrical impedance tomography for domains with known internal boundaries. Physiol Meas 31:14991516
46. Shewchuk JR (1994) An introduction to the conjugate gradient method without the agonizing pain. Carnegie Mellon University, Pittsburgh
47. Soleimani M, Gomez-Laberge C, Adler A (2006) Imaging of conductivity changes and electrode movement in EIT. Physiol Meas 27:103113
48. Thompson D, Harmany Z, Marciat R (2011) Sparse video recovery using linearly constrained gradient projection. In: IEEE international conference on acoustic, speech and signal processing (ICASSP), Prague, Czech, 2227 May 2011
49. Tibshirani R (1996) Regression shrinkage and selection via the lasso. J R Stat Soc B 58:267288
50. Vauhkonen M, Vadasz D, Karjalainen PA, Somersalo E, Kaipio JP (1998) Tikhonov regularization and prior information in electrical impedance tomography. IEEE Trans Med Imag 17(2):285293
51. Vauhkonen PJ, Vauhkonen M, Savolainen T, Kaipio JP (1999) Three dimensional electrical impedance tomography based on the complete electrode model. IEEE Trans Biomed Eng 46:11501160
52. Wright SJ, Nowak RD, Figueiredo MAT (2009) Sparse reconstruction by separable approximation. IEEE Trans Signal Process 57:24792493
53. Xiao B, Harker PT (1989) Perturbation results for the linear complementarity problem. Appl Math Lett 2(4):401405
54. Xu X, Li E, Yu H, Gong W, Han S (2014) Morphology separation in ghost imaging via sparsity constraint. Opt Express 22(12):1437514381
55. Yang CL, Wei HY, Adler A, Soleimani M (2013) Reducing computational costs in large scale 3D EIT by using a sparse Jacobian matrix with block-wise CGLS reconstruction. Physiol Meas 34:645658
56. Yorkey TJ, Webster JG, Tompkins WJ (1987) Comparing reconstruction algorithms for electrical impedance tomography. IEEE Trans Biomed Eng 34:843852
57. Yu Y, Hong M, Liu F, Wang H, Crozier S (2010) Comparison and analysis of nonlinear algorithms for compressed sensing in MRI. In: 32nd annual international conference of the IEEE EMBS, Buenos Aires, Argentina, September 2010
Ashkan Javaherian received a B.Sc. in Biomedical Engineering, majoring in Bioelectric, and an M.Sc. in Medical Radiation Engineering from Shiraz University, Iran, in 2007 and 2010, respectively. From 2010 to 2013, he worked as a collaborative researcher with School of Mechanical Engineering, Shiraz University, Iran. Since 2013, he has worked as a junior scientist in Institute of Technical Medicine (ITEM), Faculty of Medical and Life Sciences, Furtwangen University of Applied Sciences, Germany. His research interest is mathematics of nonionizing imaging modalities. He has also worked in industry on software and electronics of medical imaging systems.
Manuchehr Soleimani received B.Sc. degree in electrical engineering and the M.Sc. degree in biomedical engineering, and the Ph.D. degree in inverse problems and electromagnetic tomography from the University of Manchester, Manchester, U.K., in 2005. From 2005 to 2007, he was a Research Associate with the School of Materials, University of Manchester.In 2007, he joined the Department of Electronic and Electrical Engineering, University of Bath, Bath, U.K., where he was a Research Associate and he then became a Lecturer in 2008, Senior Lecturer in 2013 and Reader in 2015. He has authored/co-authored over 100 journal papers and holds several patents. In 2011, he found the Engineering Tomography Laboratory (ETL) at the University of Bath, working in tomographic imaging.
Knut Moeller received a M.S. and a Ph.D. degree in computer science, and a M.D. in medicine from the University of Bonn, Bonn, Germany, in 1986, 1991, and 1996, respectively. From 1991 to 1997, he was an Assistant Professor in the Department of Computer Science, Bonn University, Bonn, Germany, where he was involved in the elds of machine learning, robotics, and image processing.In 1998, he became a Professor of Medical Informatics at Furtwangen University, Villingen-Schwenningen, Germany, where he currently is Director of the Institute of Technical Medicine (ITeM) and the Head of the Biomedical Engineering Division. His research interests include decision support systems, modeling, and signal analysis with the application to lung protective mechanical ventilation. Prof. Mller is a member of the German Society of Biomedical Engineering (DGBMT) and of the German Association for Electrical, Electronic & Information Technologies (VDE).
1 3
International Federation for Medical and Biological Engineering 2016