A Fidelityembedded Regularization Method for Robust Electrical Impedance Tomography
Abstract
Electrical impedance tomography (EIT) provides functional images of an electrical conductivity distribution inside the human body. Since the 1980s, many potential clinical applications have arisen using inexpensive portable EIT devices. EIT acquires multiple transimpedance measurements across the body from an array of surface electrodes around a chosen imaging slice. The conductivity image reconstruction from the measured data is a fundamentally illposed inverse problem notoriously vulnerable to measurement noise and artifacts. Most available methods invert the illconditioned sensitivity or Jacobian matrix using a regularized leastsquares datafitting technique. Their performances rely on the regularization parameter, which controls the tradeoff between fidelity and robustness. For clinical applications of EIT, it would be desirable to develop a method achieving consistent performance over various uncertain data, regardless of the choice of the regularization parameter. Based on the analysis of the structure of the Jacobian matrix, we propose a fidelityembedded regularization (FER) method and a motion artifact removal filter. Incorporating the Jacobian matrix in the regularization process, the new FER method with the motion artifact removal filter offers stable reconstructions of highfidelity images from noisy data by taking a very large regularization parameter value. The proposed method showed practical merits in experimental studies of chest EIT imaging.
1 Introduction
EIT is a noninvasive realtime functional imaging modality for the continuous monitoring of physiological functions such as lung ventilation and perfusion. The image contrast represents the time change of the electrical conductivity distribution inside the human body. Using an array of surface electrodes around a chosen imaging slice, the imaging device probes the internal conductivity distribution by injecting electrical currents at tens or hundreds of kHz. The injected currents (at safe levels) produce distributions of electric potentials that are noninvasively measured from the attached surface electrodes. A portable EIT system can provide functional images with an excellent temporal resolution of tens of frames per second. EIT was introduced in the late 1970s[1, 3, 2, 4], likely motivated by the success of Xray CT. Numerous image reconstruction methods and experimental validations have demonstrated its feasibility[5, 6, 7, 8] and clinical trials have begun especially in lung ventilation imaging and pulmonary function testing[9]. However, EIT images often suffer from measurement noise and artifacts especially in clinical environments and there still exist needs for new image reconstruction algorithms to achieve both high image quality and robustness.
The volume conduction or lead field theory indicates that a local perturbation of the internal conductivity distribution alters the measured current–voltage or transimpedance data, which provide core information for conductivity image reconstructions. The sensitivity of the data to conductivity changes varies significantly depending on the distance between the electrodes and the conductivity perturbation: the measured data are highly sensitive to conductivity changes near the currentinjection electrodes, whereas the sensitivity drops rapidly as the distance increases[10]. In addition, the boundary geometry and the electrode configuration also significantly affect the measured boundary data.
The sensitivity map between the measured data and the internal conductivity perturbation is the basis of the conductivity image reconstruction. A discretization of the imaging domain results in a sensitivity or Jacobian matrix, which is inverted to produce a conductivity image. The major difficulty arises from this inversion process, because the matrix is severely illconditioned. Conductivity image reconstruction in EIT is, therefore, known to be a fundamentally illposed inverse problem. Uncertainties in the body shape, body movements, and electrode positions are unavoidable in practice, and result in significant amounts of forward and inverse modeling errors. Measurement noise and these modeling errors, therefore, may deteriorate the quality of reconstructed images.
Numerous image reconstruction algorithms have been developed to tackle this illposed inverse problem with data uncertainties[11, 12, 13]. In the early 1980s, an EIT version of the Xray CT backprojection algorithm[3, 15, 16] was developed based on a careful understanding of the CT idea. The onestep GaussNewton method[17] is one of the widely used methods and often called the linearized sensitivity method. Several direct methods were also developed: the layer stripping method[18] recovered the conductivity distribution layer by layer, the Dbar method[19, 20] was motivated from the constructive uniqueness proof for the inverse conductivity problem[21], and the factorization method[22, 23] was originated as a shape reconstruction method in the inverse scattering problem[24]. Recently, the discrete cosine transform was adopted to reduce the number of unknowns of the inverse conductivity problem[25]. There is an open source software package, called EIDORS, for forward and inverse modelings of EIT[26]. There are also novel theoretical results showing a unique identification of the conductivity distribution under the ideal model of EIT[27, 28, 29, 21, 30, 31].
Most common EIT image reconstruction methods are based on some form of leastsquares inversion, minimizing the difference between the measured data and computed data provided by a forward model. Various regularization methods are adopted for the stable inversion of the illconditioned Jacobian matrix. These regularized leastsquares datafitting approaches adjust the degree of regularization by using a parameter, controlling the tradeoff between data fidelity and stability of reconstruction. Their performances, therefore, depend on the choice of the parameter.
In this paper, we propose a new regularization method, which is designed to achieve satisfactory performances in terms of both fidelity and stability regardless of the choice of the regularization parameter. Investigating the correlations among the column vectors of the Jacobian matrix, we developed a new regularization method in which the structure of data fidelity is incorporated. We also developed a motion artifact removal filter, that can be applied to the data before image reconstructions, by using a submatrix of the Jacobian matrix. After explaining the developed methods, we will describe experimental results showing that the proposed fidelityembedded regularization (FER) method combined with the motion artifact removal filter provides stable image reconstructions with satisfactory image quality even for very large regularization parameter values, thereby making the method irrelevant to the choice of the parameter value.
2 Method
2.1 Representation of Measured Data
To explain the EIT image reconstruction problem clearly and effectively, we restrict our description to the case of a 16channel EIT system for realtime timedifference imaging applications. The sixteen electrodes () are attached around a chosen imaging slice, denoted by . We adopt the neighboring data collection scheme, where the device injects current between a neighboring electrode pair and simultaneously measures the induced voltages between all neighboring pairs of electrodes for . Here, we denote . Let be the electrical conductivity distribution of , and let denote the boundary surface of . The electrical potential distribution corresponding to the th injection current, denoted by , is governed by the following equations:
(1) 
where is the outward unit normal vector to , is the electrode contact impedance of the th electrode , is the potential on subject to the th injection current, and is the amplitude of the injection current between and . Assuming that , the measured voltage between the electrode pair subject to the th injection current at time is expressed as:
(2) 
The voltage data in a clinical environment are seriously affected by the following unwanted factors[32, 33, 13]:

unknown and varying contact impedances making the data for the th injection current unreliable and

inaccuracies in the moving boundary geometry and electrode positions.
Here, we set and . Among those sixteen voltage data for each injection current, thirteen are measured between electrode pairs where no current is injected, that is, the normal components of the current density are zero. For these data, it is reasonable to assume to get in (1). Discarding the remaining three voltage data, which are sensitive to changes in the contact impedances, we obtain the following voltage data at each time:
(3) 
The total number of measured voltage data is at each time. The data vector reflects the conductivity distribution , body geometry , electrode positions , and data collection protocol.
Neglecting the contact impedances underneath the voltagesensing electrodes where no current is injected, the relation between and is expressed approximately as:
(4) 
where is a position in and is the area element. Note that since depends on time, also (implicitly) depends on time as well. For timedifference conductivity imaging, we take derivative with respect to time variable to obtain
(5) 
where and denote the timederivatives of and , respectively. The data depends nonlinearly on because of the nonlinear dependency of on . This can be linearized by replacing with a computed potential:
(6) 
where is the computed potential induced by the th injection current with a reference conductivity . We can compute by solving the following boundary value problem:
(7) 
For simplicity, the reference conductivity is set as . Simplifying (6) as
(8) 
denotes the timechange of the measured voltage data vector at time and denotes
(9) 
2.2 Sensitivity Matrix
Computerized image reconstructions require a crosssectional imaging plane (or electrode plane) to be discretized into finite elements (), where is the th element or pixel (Fig. 1). Recent advances in 3D scanner technology allow the boundary shape and electrode positions to be captured. If denotes the discretized version of the change of the conductivity distribution at time , the standard linearized reconstruction algorithm is based on the following linear approximation:
(10) 
where is the sensitivity matrix (or Jacobian matrix) given by
(11) 
The sensitivity matrix is precomputed assuming a homogeneous conductivity distribution in the imaging plane (Fig. 2). The numbers of rows and columns of are 208 and , respectively. For simplicity, we used a 2D forward model of the crosssection to compute . The th column of comprises the changes in the voltage data subject to a unit conductivity change in the th pixel .
Inevitable discrepancies exist between the forward model output and the measured data due to modeling errors and measurement noise: the real background conductivity is not homogeneous, the boundary shape and electrode positions change with body movements and there always exist electronic noise and interferences.
2.3 Motion Artifact Removal
Motion artifacts are inevitable in practice to produce errors in measured voltage data and deteriorate the quality of reconstructed images[34, 35]. To investigate how the motion artifacts influence the measured voltage data, we take timederivative to both sides of (4) assuming the domain varies with time. It follows from the Reynolds transport theorem that
(12) 
where is the outwardnormal directional velocity of . Note that the last term of (2.3) is the voltage change due to the boundary movement. Using the chain rule, the first term of the righthand side of (2.3) is expressed as
(13) 
It follows from the integration by parts and (1) that
(14)  
(15) 
Here, we neglected the contact impedances underneath the voltagesensing electrodes and approximated . From (13)(15), (2.3) can be expressed as:
(16) 
Note that (2.3) becomes (5) when the boundary does not vary with time (). We linearize (2.3) by replacing with the computed potential in (7):
(17) 
After discretization, (2.3) can be expressed as:
(18) 
where is given by
with .
Compared to (10), (18) has the additional term that is the (linearized) error caused by the boundary movement. This term multiplied by the strong sensitivities on the boundary becomes a serious troublemaker, and cannot be neglected in (18) because the vectors for have large magnitudes.
To filter out the uncertain data related with motion artifacts from , we introduce the boundary sensitive Jacobian matrix , which is a submatrix of consisting of all columns corresponding to the triangular elements located on the boundary. The boundary movement errors in the measured data are, then, assumed to be in the column space of . The boundary errors are extracted by
where is a regularization parameter and is the identity matrix. Then, the motion artifact is filtered out from data by subtraction . The proposed motion artifact removal is performed before image reconstruction using any image reconstruction method. In this paper, the filtered data were used in places of for all image reconstructions.
2.4 Main Result : Fidelityembedded Regularization (FER)
Severe instability arises in practice from the illconditioned structure of when some form of its inversion is tried. To deal with this fundamental difficulty, the regularized leastsquares datafitting approach is commonly adopted to compute
(19) 
with a suitably chosen regularization parameter and regularization operator . Such image reconstructions rely on the choice of (often empirically determined) and using a priori information, suffering from over or underregularization.
We propose the fidelityembedded regularization (FER) method:
(20) 
where the regularization operator is the diagonal matrix such that
(21) 
To explain the FER method, we closely examine the correlations among column vectors of the sensitivity matrix , described in Fig. 2. The correlation between and can be expressed as
(22) 
where denotes the standard inner product. Here, is a solution of
(23) 
where is the characteristic function having on and otherwise. The identity (22) follows from
(24) 
for [23]. This shows that the column vector is like an EEG (electroencephalography) data induced by dipole sources with directions at locations . Given that two dipole sources at distant locations produce mutually independent data, the correlation between and decreases with the distance between and . Fig. 2 shows a few images of the correlation as a function of for four different positions . The correlation decreases rapidly as the distance increases. In the green regions where the correlation is almost zero, is nearly orthogonal to .
Fig. 2 shows that if and are far from each other, the corresponding columns of the sensitivity matrix are nearly orthogonal. This somewhat orthogonal structure of the sensitivity matrix motivates an algebraic formula that directly computes the local ensemble average of conductivity changes at each point using the inner product between changes in the data and a scaled sensitivity vector at that point:
(25) 
where is the weighted average conductivity at the th element and the weight is expressed in terms of the correlations between columns of . It turns out that this simple formula shows a remarkable performance in terms of robustness, but requires a slight compromise in spatial resolution.
Substituting into (25), the relation between and can be expressed as the following convolution form:
(26) 
where . The nonzero scaling factor is designed for normalization. The kernel satisfies the following:

for each , due to the nonzero scaling factor.

decreases as the distance between and increases (except near boundary where strong sensitivity arises).
Hence, roughly behaves like a blurred version of the Dirac delta function. This is the reason why the formula (25) directly computes the local ensemble average of conductivity changes at each point.
The algebraic formula (25) can be seen as a regularized leastsquares datafitting method (19) when the regularization operator is . Then, the formula (25) can be expressed using in (21) as
(27) 
This can be formulated similarly as (19) for an extremely large value of () when :
(28) 
Here, is a devised scaling term to prevent the reconstructed image from becoming zero when goes to infinity.
The FER method in (20) was proposed based on (27) and (28). When the regularization parameter is small (), the FER method is equivalent to the regularized leastsquares datafitting method (19). When is large (), it converges to the algebraic formula (25) and directly recovers the weighted average conductivity . The regularization operation fully exploits the somewhat orthogonal structure of the sensitivity matrix, thereby embedding data fidelity in the regularization process. Adopting this carefully designed , the FER method provides stable conductivity image reconstructions with high fidelity even for very large regularization parameter values.
3 Results
We applied the FER method to experimental data to show its performance. We acquired the boundary geometry and electrode positions as accurate as possible to reduce forward modeling uncertainties[33]. A handheld 3D scanner was used to capture the boundary shape of the thorax and electrode positions (Fig. 1). Then, we set the electrode plane as the horizontal crosssection of the 3Dscanned thorax containing the attached electrodes (Fig. 1). The finite element method was employed to compute the sensitivity matrix by discretizing the imaging slice. Here, we used a mesh with 12,001 nodes and 23,320 triangular elements for subject A and a different mesh with 13,146 nodes and 25,610 triangular elements for subject B.
Figs. 3 and 4 compare the performance of the proposed FER method in (20) with the standard regularized leastsquares method ((19) when is the identity matrix). The regularization parameter of the standard method was heuristically chosen for its best performance, and the parameter of the FER method was set to be one of three different values . The injection current was 1 mA at 100 kHz, and the frame rate was 9 frames per second. The reference frame at was obtained from the maximum expiration state. The measured data, , represent the voltage differences between each time and . The blue regions, which denote where conductivity decreased by inhaled air, increased during inspiration and decreased during expiration. The FER method with was clearly more robust than the standard method that produced more artifacts originated from the inversion process.
Standard 

FER () 

FER () 

FER () 

Standard 

FER () 

FER () 

FER () 

Note that the degree of orthogonality of the columns of the sensitivity matrix depends on the current injection pattern. This makes the performance of the FER method depend on the current injection pattern since it incorporates the structure of the sensitivity matrix in the regularization process. For example, if we inject currents between diagonal pairs of electrodes, the corresponding sensitivity matrix becomes less orthogonal compared with that using the neighboring current injection protocol, thus producing more blurred images. In contrast, a more narrower injection angle, for example, using a 32channel EIT system with the adjacent injection pattern, can enhance the orthogonality of the corresponding sensitivity matrix. However, the narrower injection angle results in poor distinguishability[14] and may deteriorate the image quality. To maximize the performance of the FER method, balancing between orthogonality of the sensitivity matrix and distinguishability should be considered when designing a data collection protocol.
The direct algebraic formula (25) or (20) for can be expressed as a transpose of a scaled sensitivity matrix. This type of direct approach was suggested in the late 1980s by Kotre[36], but was soon abandoned owing to poor performance and lack of theoretical grounding[37]. Since then, regularized inversion of the sensitivity matrix has been the main approach for EIT image reconstruction. Kotre’s method using the normalized transpose of was regarded as an extreme version of the backprojection algorithm in EIT[33], in the sense of the Radon transform in CT; in the case when is the Radon transform of CT, its adjoint is known as the backprojector. With this interpretation, it seems that Kotre’s method is very sensitive to forward modeling uncertainties. In the FER method, is a scaled version of the adjoint of the Jacobian, which can be viewed as a blurred version of the Dirac delta function. The FER method with becomes a direct method for robust conductivity image reconstructions without inverting the sensitivity matrix .
4 Conclusion
Since the 1980s, many EIT image reconstruction methods have been developed to overcome difficulties in achieving robust and consistent images from patients in clinical environments. Recent clinical trials of applying EIT to mechanically ventilated patients have shown its feasibility as a new realtime bedside imaging modality. They also request, however, more robust image reconstructions from patients’ data contaminated by noise and artifacts. The proposed FER method achieves both robustness and fidelity by incorporating the structure of the sensitivity matrix in the regularization process. Unlike most other algorithms, the FER method also offers direct image reconstructions without matrix inversion. This has a practical advantage especially in clinical environments since the direct method does not require any adjustment of regularization parameters. In addition to timedifference conductivity imaging, the FER method enables robust spectroscopic admittivity imaging of both conductivity and permittivity, and can employ frequencydifference approaches as well.
Although we showed that in (26) provides a satisfactory approximation to , it is difficult to estimate its accuracy. This is related to the kernel and the structure of the sensitivity matrix in which the current injection pattern is incorporated. It is a challenging issue to mathematically characterize how precisely the kernel approximates the Dirac delta function with respect to the current inject pattern.
Acknowledgment
K.L. and J.K.S. were supported by the National Research Foundation of Korea (NRF) grant 2015R1A5A1009350. E.J.W. was supported by a grant of the Korean Health Technology R&D Project (HI14C0743).
References
 [1] R. P. Henderson and J. G. Webster, “An impedance camera for spatially specific measurements of the thorax”, IEEE Trans. Biomed. Eng., vol. BME25, no. 3, pp. 250254, 1978.
 [2] A. P. Calderón “On an inverse boundary value problem”, In seminar on numerical analysis and its applications to continuum Physics, Soc. Brasileira de Matemàtica, pp. 6573, 1980.
 [3] D. C. Barber and B. H. Brown “Applied potential tomography”, J. Phys. E. Sci. Instrum., vol. 17, no. 9, pp. 723733, 1984.
 [4] T. Yorkey, J. Webster, and W. Tompkins, “Comparing reconstruction algorithms for electrical impedance tomography”, IEEE Trans. Biomed. Engr., vol. 34, no. 11, pp. 843852, 1987.
 [5] P. Metherall, D. C. Barber, R. H. Smallwood, and B. H. Brown, “Threedimensional electrical impedance tomography”, Nature, vol. 380, pp. 509512, 1996.
 [6] M. Cheney, D. Isaacson, and J. C. Newell, “Electrical impedance tomography”, SIAM Review, vol. 41, no. 1, pp. 85101, 1999.
 [7] C. Putensen, H. Wrigge, J. Zinserling, “Electrical impedance tomography guided ventilation therapy”, Critical Care, vol. 13, no. 3, pp. 344350, 2007.
 [8] T. Meier et al., “Assessment of regional lung recruitment and derecruitment during a PEEP trial based on electrical impedance tomography”, Intensive Care Med., vol. 34, no. 3, pp. 543550, 2008.
 [9] I. Frerichs et al., “Chest electrical impedance tomography examination, data analysis, terminology, clinical use and recommendations: concesus statement of the TRanslational EIT developmeNt stuDy group”, Thorax at press, 2016.
 [10] D. C. Barber and B. H. Brown, “Errors in reconstruction of resistivity images using a linear reconstruction technique”, Clin. Phys. Physiol. Meas., vol. 9, no. Supple A, pp. 101104, 1988.
 [11] D. S. Holder, Electrical impedance tomography: Methods, History and Applications, Bristol and Philadelphia:IOP Publishing, 2005
 [12] J. L. Müller and S. Siltanen, Linear and nonlinear inverse problems with practical applications, USA:SIAM, 2012.
 [13] J. K. Seo and E. J. Woo, Nonlinear Inverse Problems in Imaging, Wiley, 2013.
 [14] D. Isaacson, “Distinguishability of conductivities by electric current computed tomography”, IEEE Trans. Med. Imaging, vol. MI5, nol. 2, pp. 9195, 1986.
 [15] B. H. Brown, D. C. Barber, and A. D. Seagar, “Applied potential tomography: possible clinical applications”, Clin. Phys. Physiol. Meas., vol. 6, no. 2, pp. 109121, 1985
 [16] F. Santosa and M. Vogelius, “A backprojection algorithm for electrical impedance imaging”, SIAM J. Appl. Math., vol. 50, no. 1, pp. 216243, 1990
 [17] M. Cheney, D. Isaacson, J. Newell, S. Simske, and J. Goble, “NOSER: An algorithm for solving the inverse conductivity problem”, Internat. J. Imaging Systems and Technology, vol. 2, no. 2, pp. 6675, 1990.
 [18] E. Somersalo, M. Cheney, D. Isaacson, and E. Isaacson, “Layer stripping: a direct numerical method for impedance imaging”, Inverse Problems, vol. 7, no. 6, pp. 899926, 1991.
 [19] S. Siltanen, J. Mueller, and D. Isaacson, “An implementation of the reconstruction algorithm of A Nachman for the 2D inverse conductivity problem”, Inverse Problems, vol. 16, no. 3, pp. 681699, 2000.
 [20] J. L. Mueller, S. Siltanen, D. Isaacson, “A direct reconstruction algorithm for electrical impedance tomography”, IEEE Trans. Med. Imag., vol. 21, no. 6, pp. 555559, 2002.
 [21] A. Nachman, “Global uniqueness for a twodimensional inverse boundary problem”, Ann. Math., vol. 143, no. 1, pp. 7196, 1996.
 [22] M. Hanke, M. Brühl, “Recent progress in electrical impedance tomography”, Inverse Problems, vol. 19, no. 6, pp. 126, 2003.
 [23] M. K. Choi, B. Bastian, and J. K. Seo, “Regularizing a linearized EIT reconstruction method using a sensitivitybased factorization method”, Inverse Probl. Sci. En., vol. 22, no. 7, pp. 10291044, 2014.
 [24] A. Kirsch, “Characterization of the shape of a scattering obstacle using the spectral data of the far field operator”, Inverse Problems, vol. 14, no. 6, pp. 14891512, 1998.
 [25] B. Schullcke et al., “Structuralfunctional lung imaging using a combined CTEIT and a Discrete Cosine Transformation reconstruction method”, Scientific Reports, vol. 6, pp. 112, 2016.
 [26] A. Adler and W. R. B. Lionheart, “Uses and abuses of EIDORS: An extensible software base for EIT”, Physiol. Meas., vol. 27, no. 5, pp. S25S42, 2006.
 [27] R. Kohn and M. Vogelius, “Determining conductivity by boundary measurements”, Comm. Pure Appl. Math., vol. 37, pp. 113123, 1984.
 [28] J. Sylvester and G. Uhlmann, “A global uniqueness theorem for an inverse boundary value problem”, Ann. Math., vol. 125, pp. 153169, 1987.
 [29] A. Nachman, “Reconstructions from boundary measurements”, Ann. Math., vol. 128, pp. 531576, 1988.
 [30] K. Astala, L. Päivärinta, “Calderon’s inverse conductivity problem in the plane”, Ann. Math., vol. 163, pp. 265299, 2006.
 [31] C. Kenig, J. Sjostrand, and G. Uhlmann, “The Calderon problem with partial data”, Ann. Math., vol. 165, pp. 567–591, 2007.
 [32] V. Kolehmainen, M. Vauhkonen, P. A. Karjalainen, J. P. Kaipio, “Assessment of errors in static electrical impedance tomography with adjacent and trigonometric current patterns”, Physiol. Meas., vol. 18, no. 4, pp. 289303, 1997.
 [33] W. R. B. Lionheart, “EIT reconstruction algorithms: pitfalls, challenges and recent developments”, Physiol. Meas., vol. 25, no. 1, pp. 125142, 2004.
 [34] A. Adler, R. Guardo, Y. Berthiaume, “Impedance Imaging of Lung Ventilation: Do We Need to Account for Chest Expansion?”, IEEE trans. Biomed. Eng., vol. 43, no. 4, pp. 414420, 2005.
 [35] J. Zhang, R. P. Patterson, “EIT images of ventilation: what contributes to the resistivity changes?”, Physiol. Meas., vol. 26, no. 2, pp. S81S92, 2005.
 [36] C. J. Kotre, “A sensitivity coefficient method for the reconstruction of electrical impedance tomograms”, Clin. Phys. Physiol. Meas., vol. 10, no. 3, pp. 275281, 1989.
 [37] D. C. Barber, “A sensitivity method for electrical impedance tomography”, Clin. Phys. Physiol. Meas., vol. 10, no. 4, pp. 368369, 1989.