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

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.


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][2][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][5][6][7][8][9][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][12][13][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 goodquality 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][16][17].
Homotopy perturbation method is a novel and effective iteration method, which has been successfully applied to solve various of nonlinear equations [18][19][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.

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,  = ( 1 ,  2 , . . .,   ) ∈ R  and  = ( 1 ,  2 , . . .,   ) ∈ R  , respectively.

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  1 (Ω) we denote the  2 (Ω)-based Sobolev space with smoothness index 1, and H1 (Ω) = { ∈  1 (Ω) : ∫ Ω   = 0}.For convenience, let us denote the space of length  by R  ⬦ = { ∈ R  : ∑  =1   = 0}.The weak formulation of the forward problem can be stated as follows: given a current vector  ∈ R  ⬦ , a conductivity distribution , and positive contact impedance   , to find the solution (, ) ∈  1 where the strictly elliptic bilinear form  is defined by Indeed, the existence and uniqueness of a solution (, ) ∈  1 (Ω) ⊕ R  ⬦ 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.

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,   :  ∈ R  ⬦ →  ∈ R  ⬦ , 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    = Υ   ∈ R  ⬦ , where Υ  : H1 (Ω) → R  ⬦ 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  obs , 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  obs , also called image reconstruction.

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 (V, ) :  × [0, 1] →  is constructed for (14), which satisfies where  ∈ [0, 1] is an embedding homotopy parameter and  0 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 (V, 0) = 0 is continuously deformed to the original problem (V, 1) = 0 that is, V changed from  0 to the approximation solution .In topology, this is called deformation, and  −  0 and   ()  (() −  obs ) 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  → 1, the approximate solution of ( 14) is obtained with We rewrite (15) with the Taylor expansion of () at  0 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  0 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 firstorder approximation, as follows, for a given prior  0 , 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  0 , 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,   = 1/ max (     ), where   is the Jacobian matrix of () calculated at the th iteration   and  max (     ) 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  obs we employ the discrepancy principle as a stopping rule for HPIM, which determines the stopping index  * =  * (,  obs ) by for some sufficiently large  > 0. In fact, this stopping rule renders this method as a regularization method.

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.
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.
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.(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.

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 noisecontaminated 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  = 20, 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.

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.

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

4. 1 .Figure 3 :
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.

Figures 3 (
Figures 3(a), 3(b), 4(a), 4(b), 5(a), and 5(b), using(27), are calculated and listed in Table1.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

Figure 4 :
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.

Figure 5 :
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.

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