Journal of Applied Mathematics

Volume 2013 (2013), Article ID 454706, 11 pages

http://dx.doi.org/10.1155/2013/454706

Research Article

Image Reconstruction Based on Homotopy Perturbation Inversion Method for Electrical Impedance Tomography

Jing Wang and Bo Han

Department of Mathematics, Harbin Institute of Technology, Harbin 150001, China

Received 9 May 2013; Revised 26 November 2013; Accepted 27 November 2013

Academic Editor: Md. Sazzad Chowdhury

Copyright © 2013 Jing Wang and Bo Han. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The image reconstruction for electrical impedance tomography (EIT) mathematically is a typed nonlinear ill-posed inverse problem. In this paper, a novel iteration regularization scheme based on the homotopy perturbation technique, namely, homotopy perturbation inversion method, is applied to investigate the EIT image reconstruction problem. To verify the feasibility and effectiveness, simulations of image reconstruction have been performed in terms of considering different locations, sizes, and numbers of the inclusions, as well as robustness to data noise. Numerical results indicate that this method can overcome the numerical instability and is robust to data noise in the EIT image reconstruction. Moreover, compared with the classical Landweber iteration method, our approach improves the convergence rate. The results are promising.

1. Introduction

Electrical impedance tomography (EIT) is an imaging modality which seeks to reconstruct the electrical conductivity distribution inside the interested objection from the surface electrical measurements. The process of generating an image from the measurements is called image reconstruction. Being a noninvasive, portable, and inexpensive methodology, the EIT technique has been extensively applied in medical imaging, geophysical exploration, nondestructive testing of materials, and so forth [1–3].

However, the image reconstruction for EIT mathematically is by nature a nonlinear and ill-posed inverse problem, which makes the image reconstruction a challenging task. Nonlinearity increases the requirement for using iterative image reconstruction techniques rather than noniterative ones, while ill-posedness makes any numerical reconstruction method require some form of regularization technique to obtain stable solution. Image quality greatly depends on the performance of the forward solver and the measurement errors but also depends on the regularization technique and the converging nature of iterative algorithm used in inverse problem. At present, various reconstruction algorithms have been developed to improve the image quality and increase the speed of inversion algorithms [4–10]. In the past few years, iterative image reconstruction with an updated forward solution has been investigated. In particular, the Landweber iteration method is one of the most straightforward and popular iterative algorithms for image reconstruction problem [11–14]. The advantages of this algorithm include the simple implementation and low computational cost, due to the fact that only the gradient information of the data fitting is used in the process of implementing the image reconstruction. Nevertheless, the ill-conditioning of the sensitivity matrix makes the convergence of Landweber iteration algorithm relatively slow and need a lot of iterations to get a good-quality image. Landweber iteration has a semiconvergence characteristic; that is, reconstructed image error often starts to increase gradually after getting to local minimum. Some preconditioning methods have been proposed to speed up the convergence rate [15–17].

Homotopy perturbation method is a novel and effective iteration method, which has been successfully applied to solve various of nonlinear equations [18–20] in various branches of science and engineering. Recently, it has also been extended for dealing with nonlinear ill-posed inverse problems [21, 22]. The essential idea of this method is to introduce a homotopy parameter and combine the traditional perturbation method with the homotopy technique. One of the most notable features of homotopy perturbation method is that usually only a few perturbation terms can be sufficient to obtain a reasonably accurate solution. In this paper, a novel iteration regularization scheme based on the homotopy perturbation technique, namely, homotopy perturbation inversion method, is applied to investigate the EIT inverse problem for the first time. Numerical simulations are carried out to illustrate the feasibility and effectiveness of this method. Simulation results indicate that this method can overcome the numerical instability and is robust to data noise in the EIT image reconstruction. Moreover, in contrast to the classical Landweber iteration method, our approach improves the convergence rate. As a result an efficient method for EIT image reconstruction is introduced.

The outline of this paper is as follows. In Section 2, we briefly describe the mathematical model for the EIT forward and inverse problem. In Section 3, the homotopy perturbation inversion scheme is presented, which is an iteration regularization scheme. A stable reconstruction is guaranteed provided the step length parameter is chosen appropriately and a stopping rule is formulated. To illustrate the well performance of this method, numerical simulations are provided in Section 4. Finally, a short conclusion is drawn in Section 5.

2. Mathematical Formulation of the EIT Problem

In this section we give a brief account on the mathematical model for EIT, which is composed of forward problem and inverse problem.

We first consider a typical acquisition experiment for EIT. Assume that electrodes with perfect conductor have been fixed around the surface of the object. Current is injected through some subset of these electrodes to excite the object under investigation, and the induced voltages are measured at all electrodes. This procedure is repeated many times with different electrodes until a sufficient amount of data has been obtained. Mathematical models for the measurement are described in [23, 24]. The complete electrode model, nowadays the standard model for medical applications, is presently the most accurate model for EIT since it takes into account the presence of the discrete electrodes and the effect of the contact impedance and can match the measurement precision of the experiment well [24]. The model is composed of the following equations:

where is the electrical conductivity, is the electrical potential, denotes the th electrode, is the effective contact impedance between the th electrode and tissue, and is the unit outward normal to and . Furthermore, and are the electrical potential and electrical current on the th electrode, respectively. In addition, in order to ensure the solvability of the elliptic boundary value problem (1) and the uniqueness of the solution, the following conditions are required to be satisfied:

We denote vectors consisting of the electrical potentials and currents on the electrodes by and ; that is, and , respectively.

2.1. The Forward Problem and Finite Element Approximation

For known conductivity distribution , injecting current , and contact impedance , the forward problem is to find the internal electrical potentials in and the voltages on the electrodes that satisfy the complete electrode model (1). In practice the system (1) cannot be solved analytically for any arbitrary distribution, and therefore numerical approximations must be used to find an approximate solution. Finite element method [25, 26] is particularly well suited to the discretization of this problem.

As we all know, the procedure in the finite element method starts with the so-called variational (Galerkin) formulation of the problem (1), also called the weak form. We first need define some notation. By we denote the -based Sobolev space with smoothness index 1, and . For convenience, let us denote the space of length by . The weak formulation of the forward problem can be stated as follows: given a current vector , a conductivity distribution , and positive contact impedance , to find the solution such that

where the strictly elliptic bilinear form is defined by

Indeed, the existence and uniqueness of a solution have been shown using the Lax-Milgram lemma in [24].

Now we consider the piecewise linear function as the finite element basis function, and the solution domain will be divided into small triangles elements. We approximate the potential distribution within as

and potentials on the electrodes as

where is the number of nodes in finite element mesh, , , and so forth and and are the coefficients to be determined. For notational reasons we identify and with their coordinates and write , . By substituting these two approximations into variational equation (3) and by choosing and , we obtain a discrete linear system

where the data vector with and , the determined coefficient vector , and the stiffness matrix

with

After solving the determined coefficients , we can use (6) to find the approximate potentials on the electrodes.

2.2. The Inverse Problem

In practice, the conductivity distribution is not known. What we know is merely all pairs of injected current data and resulting voltage data on the electrodes. The inverse problem under the complete electrode model (1) is to estimate the conductivity distribution from these data.

By Ohm’s law, the output is indeed a measurement (discretization) of the linear current-to-voltage operator that is, , which maps the applied electrode current data to the computed electrode voltage data or denoted by to emphasize the dependence on the conductivity parameter as follows:

and , where denotes the Dirichlet trace operator, which projects the solution of the forward problem to the electrode voltages .

Now we denote the solution operator by , such that

where relates the conductivity distribution with the corresponding virtual measurements , often referred to as the parameter-to-observation map, also called the forward operator. Notably, depends linearly on but nonlinearly on . For a fixed current vector , we can view as a function of only. It is important to distinguish between the PDE operator and the forward operator . Though both are associated with the model problem (1), their meaning is not the same. Indeed, describes the forward problem, while actually represents its solution.

However, in practical applications, actual observations for EIT are typically potential differences between certain electrodes, denoted by , which may not be known exactly and unavoidably noisy due to the measurement or model error. Thus the observation model response of EIT can be formulated as the following nonlinear operator equation:

where is the noise in the measurements.

Thus the inverse problem for EIT becomes that of finding the approximation conductivity distribution from the observed partial noisy knowledge , also called image reconstruction.

3. Homotopy Perturbation Inversion Method

One popular approach to deal with the inverse problem (12) is the nonlinear least-square method

The corresponding Euler equation of (13) is given by

By the homotopy perturbation technique, a fixed point homotopy function is constructed for (14), which satisfies

where is an embedding homotopy parameter and is an initial guess which incorporate a prior estimate of the true conductivity distribution . Obviously, from these definitions we will have

As the parameter changed continuously from zero to one, the trivial problem is continuously deformed to the original problem that is, changed from to the approximation solution . In topology, this is called deformation, and and are called homotopic.

Next, according to the homotopy perturbation technique, we can use the parameter as a small parameter and assume that the solution of (15) can be expressed as a power series in ,

As , the approximate solution of (14) is obtained with

We rewrite (15) with the Taylor expansion of at and ignore the higher order term, which yields

Substituting (17) into (19), we will have

Equating the coefficients for the different power of , we get the following formulations:

Then starting with an initial prior estimate and in turn solving (21), we can obtain

Thus by (18) we can get the approximate solution of (14); that is,

The above formulation will allow us to construct two types of iteration methods for solving the nonlinear inverse problem (12). The first iteration scheme is to use the first-order approximation, as follows, for a given prior ,

which is the well-known classical Landweber iteration method (CLIM) [27], viewed as one of the variations of the steepest gradient descent method, and the second iteration scheme is to use the second-order approximation, as follows, for a given prior ,

where is an appropriately chosen step length. Since the iteration-varying step length usually shows better performance than the constant value, the step length for iterations (24) and (25) was chosen with the same criterion; that is, , where is the Jacobian matrix of calculated at the th iteration and denotes the maximum eigenvalue of . We here call iteration scheme (25) the homotopy perturbation inversion method (HPIM). It has faster convergence than the classical Landweber iteration (24), which has been analyzed and proved in [28].

Considering the measurement errors of the observation data we employ the discrepancy principle as a stopping rule for HPIM, which determines the stopping index by

for some sufficiently large . In fact, this stopping rule renders this method as a regularization method.

4. Numerical Simulations

In this section, some numerical results for simulated data were presented to demonstrate the feasibility and effectiveness of the homotopy perturbation inversion method (HPIM). The quality of the images reconstructed by HPIM is compared with that of the classical Landweber iteration method (CLIM). All the simulations are 2D, as this allows an easier visualization of the results and faster simulation. A 16-electrode system was selected for simulations. We used constant injection current between adjacent electrodes and adjacent voltage measurement between all other electrodes. To avoid committing what is referred to as an inverse crime, two different finite element meshes were used for the forward inverse solvers, as shown in Figure 1, and the small green strips on the perimeter of the meshes show the position of the ideal electrodes. A mesh with 1968 elements and 1049 nodes was used for the forward simulations. For the inverse computations another mesh with 492 elements and 279 nodes was used to reduce the computational burden. Considering the simplicity’s sake of mesh refinement, we divided one element used in the inverse calculation into four elements to generate the forward mesh. In the reconstructions, contact impedances under the electrodes were assumed to be known and are set to be 0.05 for all electrodes.

(a)

(b)

Figure 1: The forward and inverse models. (a) the forward model with 1968 triangulation elements; (b) the inverse model with 492 elements.

Figure 2 shows three original phantom models used for our simulations. The phantoms from left to right have various inclusions at different locations, occupying 36 cells, 16 cells, and 28 cells on the total size (492 cells) of the unknown discrete conductivity vector, respectively. The higher (background) and lower (inclusions) conductivity distribution values are set to be 1.0 and 6.0, respectively.

(a)

(b)

(c)

Figure 2: Three simulation phantoms.

In the following, images for these three phantom models were reconstructed by the homotopy perturbation inversion method (HPIM) and classical Landweber iteration method (CLIM), respectively. For all simulations, the synthetic ideal data were generated using the finite element forward model based on the EIDORS 2D EIT software developed by Vauhkonen et al. [29]. There are a total of 256 data points being generated for the entire survey. The initial conductivity distribution values for both methods are set to the background values.

In addition to the reconstructed images, in order to compare the reconstruction performance more quantitatively, the following criteria are used:

where represents the reconstructed conductivity distribution, is the true one, and DE and RE represent the discrepancy error and relative image error, respectively.

4.1. Noise-Free Case

The numerical performances of HPIM and CLIM were first tested by the noise-free data. Figures 3(a), 3(b), 4(a), 4(b), 5(a), and 5(b) show the reconstructed results obtained from HPIM and CLIM for three different test phantoms, respectively. Columns from left to right in each figure represent the reconstructions at iterations 10, 20, and 50, respectively. It can be seen that HPIM is superior to CLIM and can reconstruct the conductivity distribution even after 10-step iteration. After 20 iterations, the reconstruction performance with HPIM shows great improvement over CLIM. To compare the performance of the above two methods more quantitatively, the relative errors of reconstructed images in Figures 3(a), 3(b), 4(a), 4(b), 5(a), and 5(b), using (27), are calculated and listed in Table 1. It can be observed that HPIM has less relative error than CLIM at the same iterations. To visibly illustrate the convergence behavior of both methods, we also draw the curves from discrepancy errors and relative errors versus iteration time for the reconstructions of these three phantoms, as shown in Figures 3(c)–5(c). It is clear that the discrepancy errors and relative errors of HPIM are consistently lower than those of CLIM. As can be expected, the convergence rate of the proposed HPIM is accelerated significantly compared to that of CLIM.

Table 1: Comparison of relative errors for the reconstructions displayed in Figures 3(a), 3(b), 4(a), 4(b), 5(a), and 5(b).

Iterations
Figure 3(a)
Figure 3(b)
Figure 4(a)
Figure 4(b)
Figure 5(a)
Figure 5(b)

10 steps
0.2150
0.2239
0.1213
0.1414
0.1718
0.1895

20 steps
0.2072
0.2208
0.1133
0.1355
0.1593
0.1839

50 steps
0.1954
0.2147
0.1066
0.1216
0.1500
0.1713

(a)

(b)

(c)

Figure 3: Reconstruction results with HPIM (see (a)) and CLIM (see (b)) for the first phantom at different iteration number, as well as convergence behavior (see (c)). (c) Left: discrepancy error; right: relative error. In both plots, the horizontal axes denote the time of the iteration process.

(a)

(b)

(c)

Figure 4: Reconstruction results with HPIM (see (a)) and CLIM (see (b)) for the second phantom at different iteration number, as well as convergence behavior (see (c)). (c) Left: discrepancy error; right: relative error. In both plots, the horizontal axes denote the time of the iteration process.

(a)

(b)

(c)

Figure 5: Reconstruction results with HPIM (see (a)) and CLIM (see (b)) for the third phantom at different iteration number, as well as convergence behavior (see (c)). (c) Left: discrepancy error; right: relative error. In both plots, the horizontal axes denote the time of the iteration process.

4.2. Noise-Contaminated Case

Studies indicate that EIT image reconstruction process is a typical ill-posed problem. A small perturbation in the input data may yield a large fluctuation of the solution, which may make the inversion results meaningless. In the data acquisition processes, the measurement data unavoidably contain noise. The robustness of reconstruction technique with respect to noise is crucial. In order to simulate real measurement environment, the noise-contaminated input data are used to evaluate the numerical performance and robustness to noise of HPIM. Usually, the synthetic Gaussian noise was generated to simulate systematic and random errors existing in the data acquisition processes. The resulting noisy data was indicated by adding noise to the synthetic ideal data, defined as

where is the noise level and is a noise vector obeying a Gaussian distribution with zero mean and standard deviation of 1.

In the following, we illustrate the performance of HPIM for the noise-contaminated case by comparing against the common CLIM. According to the stopping rule (26) with , for each phantom, the corresponding reconstructed images by HPIM under four different noise levels of 5%, 3%, 1%, and 0.5% are listed in Figures 6(a)–6(d), respectively. As can be expected, HPIM shows the favorable robustness. Since the images reconstructed by CLIM have the similar qualities of those by HPIM, we do not list them in Figure 6. We summarized more detailed computational results in Table 2, which is against the second phantom with two inclusions, and the corresponding reconstructed images are partially displayed in the middle column of Figure 6. Here and denote the stooping steps and the CPU time (in seconds), respectively. We find from Table 2 that images reconstructed by HPIM have the similar error qualities to those by CLIM, but in less steps and computational time. Additionally, it also can be seen from Figure 6 and Table 2 that with the increase of the noise levels the reconstruction quality decreases, which indicates that the accuracy of the measurement data has an effect on the imaging quality.

Table 2: Comparison of reconstructions for the second phantom with at four different noise levels.

HPIM
CLIM

DE
RE
(s)

DE
RE
(s)

5%
20
0.9974
0.1199
8.6813
43
0.9997
0.1186
17.8413

3%
27
0.5999
0.1923
12.1581
59
0.5993
0.1912
25.1848

1%
91
0.1999
0.1051
38.3450
190
0.2000
0.1034
77.1894

0.5%
234
0.0993
0.0997
101.8921
479
0.1000
0.0980
221.6608

(a)

(b)

(c)

(d)

Figure 6: Reconstruction results by HPIM for three phantoms with different noise levels. (a) ; (b) ; (c) ; (d) .

5. Conclusion

This paper presents an application of homotopy perturbation inversion method to the image reconstruction for EIT. Simulation results indicate the feasibility and effectiveness of this method. Compared with the classical Landweber iteration method, our approach reduces the computational time and produces more satisfactory effect.

Acknowledgment

This work was supported by the National Natural Science Foundation of China under Grant no. 91230119.

References

M. Cheney, D. Isaacson, and J. C. Newell, “Electrical impedance tomography,” SIAM Review, vol. 41, no. 1, pp. 85–101, 1999.

L. Borcea, “Electrical impedance tomography,” Inverse Problems, vol. 18, no. 6, pp. R99–R136, 2002.

W. Lionheart, N. Polydorides, and A. Borsic, “The reconstruction problem in electrical impedance tomography,” in Methods, History and Applications, chapter 1, IOP, Philadelphia, Pa, USA, 2005.

T. Yorkey, J. Webster, and W. Tompkins, “Comparing reconstruction algorithms for electrical impedance tomography,” IEEE Transactions on Biomedical Engineering, vol. 34, no. 11, pp. 843–852, 1987.

W. L. Weng and F. J. Dickin, “Improved modified Newton-Raphson algorithm for electrical impedance tomography,” Electronics Letters, vol. 32, no. 3, pp. 206–207, 1996.

M. Wang, “Inverse solutions for electrical impedance tomography based on conjugate gradients methods,” Measurement Science and Technology, vol. 13, no. 1, pp. 101–117, 2002.

M. Vauhkonen, D. Vadasz, P. Karjalainen, E. Somersalo, and J. Kaipio, “Tikhonov regularization and prior information in electrical impedance tomography,” IEEE Transactions on Medical Imaging, vol. 17, no. 2, pp. 285–293, 2002.

E. Chung, T. Chan, and X. Tai, “Electrical impedance tomography using level set representation and total variational regularization,” Journal of Computational Physics, vol. 205, no. 1, pp. 357–372, 2005.

A. Borsic, B. Graham, A. Adler, and W. Lionheart, “In vivo impedance imaging with total variation regularization,” IEEE Transactions on Medical Imaging, vol. 29, no. 1, pp. 44–54, 2010.

B. T. Jin, T. Khan, and P. Maass, “A reconstruction algorithm for electrical impedance tomography based on sparsity regularization,” International Journal for Numerical Methods in Engineering, vol. 89, no. 3, pp. 337–353, 2011.

W. Q. Yang, D. M. Spink, T. A. York, and H. Mccann, “An image-reconstruction algorithm based on Landweber's iteration method for electrical-capacitance tomography,” Measurement Science and Technology, vol. 10, no. 11, pp. 1065–1069, 1999.

Y. Li and W. Q. Yang, “Image reconstruction by nonlinear Landweber iteration for complicated distributions,” Measurement Science and Technology, vol. 19, no. 9, Article ID 094014, 2008.

Z. H. Guo, F. Shao, and D. C. Lv, “Sensitivity matrix construction for electrical capacitance tomography based on the difference model,” Flow Measurement and Instrumentation, vol. 20, no. 3, pp. 95–102, 2009.

Z. Liu, G. Y. Yang, N. He, and X. Y. Tan, “Landweber iterative algorithm based on regularization in electromagnetic tomography for multiphase flow measurement,” Flow Measurement and Instrumentation, vol. 27, pp. 53–56, 2012.

H. X. Wang, C. Wang, and W. L. Yin, “A pre-iteration method for the inverse problem in electrical impedance tomography,” IEEE Transactions on Instrumentation and Measurement, vol. 53, no. 4, pp. 1093–1096, 2004.

G. Lu, L. H. Peng, B. F. Zhang, and Y. B. Liao, “Preconditioned Landweber iteration algorithm for electrical capacitance tomography,” Flow Measurement and Instrumentation, vol. 16, no. 2-3, pp. 163–167, 2005.

J. D. Jang, S. H. Lee, K. Y. Kim, and B. Y. Choi, “Modified iterative Landweber method in electrical capacitance tomography,” Measurement Science and Technology, vol. 17, no. 7, pp. 1909–1917, 2006.

J.-H. He, “A coupling method of a homotopy technique and a perturbation technique for non-linear problems,” International Journal of Non-Linear Mechanics, vol. 35, no. 1, pp. 37–43, 2000.

J.-H. He, “Homotopy perturbation method for bifurcation of nonlinear problems,” International Journal of Nonlinear Sciences and Numerical Simulation, vol. 6, no. 2, pp. 207–209, 2005.

H. Jafari and S. Momani, “Solving fractional diffusion and wave equations by modified homotopy perturbation method,” Physics Letters A, vol. 370, no. 5-6, pp. 388–396, 2007.

F. Shakeri and M. Dehghan, “Inverse problem of diffusion equation by He's homotopy perturbation method,” Physica Scripta, vol. 75, no. 4, pp. 551–556, 2007.

L. Cao, B. Han, and W. Wang, “Homotopy perturbation method for nonlinear ill-posed operator equations,” International Journal of Nonlinear Sciences and Numerical Simulation, vol. 10, no. 10, pp. 1319–1322, 2009.

K.-S. Cheng, D. Isaacson, J. Newell, and D. Gisser, “Electrode models for electric current computed tomography,” IEEE Transactions on Biomedical Engineering, vol. 36, no. 9, pp. 918–924, 1989.

E. Somersalo, M. Cheney, and D. Isaacson, “Existence and uniqueness for electrode models for electric current computed tomography,” SIAM Journal on Applied Mathematics, vol. 52, no. 4, pp. 1023–1040, 1992.

E. J. Woo, P. Hua, J. G. Webster, and W. J. Tompkins, “Finite-element method in electrical impedance tomography,” Medical and Biological Engineering and Computing, vol. 32, no. 5, pp. 530–536, 1994.

N. Polydorides and W. R. B. Lionheart, “A matlab toolkit for three-dimensional electrical impedance tomography: a contribution to the electrical impedance and diffuse optical reconstruction software project,” Measurement Science and Technology, vol. 13, no. 12, pp. 1871–1883, 2002.

M. Hanke, A. Neubauer, and O. Scherzer, “A convergence analysis of the Landweber iteration for nonlinear ill-posed problems,” Numerische Mathematik, vol. 72, no. 1, pp. 21–37, 1995.

L. Cao and B. Han, “Convergence analysis of the homotopy perturbation method for solving nonlinear ill-posed operator equations,” Computers & Mathematics with Applications, vol. 61, no. 8, pp. 2058–2061, 2011.

M. Vauhkonen, W. R. B. Lionheart, L. M. Heikkinen, P. J. Vauhkonen, and J. P. Kaipio, “A matlab package for the EIDORS project to reconstruct two-dimensional EIT images,” Physiological Measurement, vol. 22, no. 1, pp. 107–111, 2001.