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

Authors: Jing Wang; Bo Han
Journal: Journal of Applied Mathematics
 (2013)

---

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


---

## Body

## 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 Section2, 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 thatL 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 L 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:
(1)∇·(σ∇u)=0,inΩ,σ∂u∂n=0,on∂Ω⋃l=1L‍el,∫elσ∂u∂ndS=Il,forl=1,2,…,L,u+zlσ∂u∂n=Ul,onel,l=1,2,…,L,
where σ is the electrical conductivity, u is the electrical potential, el denotes the lth electrode, zl is the effective contact impedance between the lth electrode and tissue, and n is the unit outward normal to ∂Ω and ∂u/∂n=∇u·n. Furthermore, Ul and Il are the electrical potential and electrical current on the lth 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:
(2)∑l=1LIl=0,∑l=1LUl=0.
We denote vectors consisting of the electrical potentials and currents on the electrodes by U and I; that is, U=(U1,U2,…,UL)∈ℝL and I=(I1,I2,…,IL)∈ℝL, respectively.


### 2.1. The Forward Problem and Finite Element Approximation

For known conductivity distributionσ, injecting current I, and contact impedance zl, the forward problem is to find the internal electrical potentials u in Ω and the voltages U 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 H1(Ω) we denote the L2(Ω)-based Sobolev space with smoothness index 1, and H~1(Ω)={u∈H1(Ω):∫∂Ω‍udS=0}. For convenience, let us denote the space of length L by ℝ⋄L={V∈ℝL:∑l=1L‍Vl=0}. The weak formulation of the forward problem can be stated as follows: given a current vector I∈ℝ⋄L, a conductivity distribution σ, and positive contact impedance zl, to find the solution (u,U)∈H1(Ω)⊕ℝ⋄L such that
(3)b((u,U),(v,V))=∑l=1LIlVl,∀(v,V)∈H1(Ω)⊕ℝ⋄L,
where the strictly elliptic bilinear form b is defined by
(4)b((u,U),(v,V))=∫Ωσ∇u·∇vdx+∑l=1L1zl∫el(u-Ul)(v-Vl)dS.
Indeed, the existence and uniqueness of a solution (u,U)∈H1(Ω)⊕ℝ⋄L have been shown using the Lax-Milgram lemma in [24].Now we consider the piecewise linear function{ϕi} as the finite element basis function, and the solution domain Ω will be divided into small triangles elements. We approximate the potential distribution u within Ω as
(5)uh=∑i=1Nαiϕi
and potentials U on the electrodes as
(6)Uh=∑l=1L-1βjnj,
where N is the number of nodes in finite element mesh, n1=[1,-1,0,…,0]T, n2=[1,0,-1,0,…,0]T∈ℝL, and so forth and αi and βj are the coefficients to be determined. For notational reasons we identify uh and Uh with their coordinates and write uh=(α1,…,αN), Uh=(β1,…,βL-1). By substituting these two approximations into variational equation (3) and by choosing v=ϕi and V=nj, we obtain a discrete linear system
(7)A(σ)b=f,
where the data vector f=(0,I^)T with 0∈ℝN and I^=(I1-I2,I1-I3,…,I1-IL), the determined coefficient vector b=(uh,Uh)T, and the stiffness matrix
(8)A(σ)=(BCCTD),
with
(9)Bij=∫Ωσ∇ϕi·∇ϕjdx+∑l=1L1zl∫elϕiϕjdS,i,j=1,2,…,N,Cil=-1zl∫elϕidS,i=1,2,…,N,l=1,2,…,L-1,Dll=1zl∫eldS,l=1,2,…,L-1.
After solving the determined coefficients b=A(σ)-1f, we can use (6) to find the approximate potentials Uh 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 I and resulting voltage data U 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σ:I∈ℝ⋄L→U∈ℝ⋄L, which maps the applied electrode current data I to the computed electrode voltage data U or denoted by U(σ) to emphasize the dependence on the conductivity parameter as follows:
(10)U(σ)=RσI
and RσI=ΥDu∈ℝ⋄L, where ΥD:H~1(Ω)→ℝ⋄L denotes the Dirichlet trace operator, which projects the solution u of the forward problem to the electrode voltages U.Now we denote the solution operator byF(σ;I), such that
(11)U(σ)=F(σ;I),
where F relates the conductivity distribution σ with the corresponding virtual measurements U, often referred to as the parameter-to-observation map, also called the forward operator. Notably, F depends linearly on I but nonlinearly on σ. For a fixed current vector I, we can view F(σ;I) as a function of σ only. It is important to distinguish between the PDE operator A(σ) and the forward operator F(σ). Though both are associated with the model problem (1), their meaning is not the same. Indeed, A(σ) describes the forward problem, while F(σ) actually represents its solution.However, in practical applications, actual observations for EIT are typically potential differences between certain electrodes, denoted byUobs, 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:
(12)Uobs=F(σ)+e,
where e is the noise in the measurements.Thusthe inverse problem for EIT becomes that of finding the approximation conductivity distribution σ from the observed partial noisy knowledge Uobs, also called image reconstruction.


## 2.1. The Forward Problem and Finite Element Approximation

For known conductivity distributionσ, injecting current I, and contact impedance zl, the forward problem is to find the internal electrical potentials u in Ω and the voltages U 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 H1(Ω) we denote the L2(Ω)-based Sobolev space with smoothness index 1, and H~1(Ω)={u∈H1(Ω):∫∂Ω‍udS=0}. For convenience, let us denote the space of length L by ℝ⋄L={V∈ℝL:∑l=1L‍Vl=0}. The weak formulation of the forward problem can be stated as follows: given a current vector I∈ℝ⋄L, a conductivity distribution σ, and positive contact impedance zl, to find the solution (u,U)∈H1(Ω)⊕ℝ⋄L such that
(3)b((u,U),(v,V))=∑l=1LIlVl,∀(v,V)∈H1(Ω)⊕ℝ⋄L,
where the strictly elliptic bilinear form b is defined by
(4)b((u,U),(v,V))=∫Ωσ∇u·∇vdx+∑l=1L1zl∫el(u-Ul)(v-Vl)dS.
Indeed, the existence and uniqueness of a solution (u,U)∈H1(Ω)⊕ℝ⋄L have been shown using the Lax-Milgram lemma in [24].Now we consider the piecewise linear function{ϕi} as the finite element basis function, and the solution domain Ω will be divided into small triangles elements. We approximate the potential distribution u within Ω as
(5)uh=∑i=1Nαiϕi
and potentials U on the electrodes as
(6)Uh=∑l=1L-1βjnj,
where N is the number of nodes in finite element mesh, n1=[1,-1,0,…,0]T, n2=[1,0,-1,0,…,0]T∈ℝL, and so forth and αi and βj are the coefficients to be determined. For notational reasons we identify uh and Uh with their coordinates and write uh=(α1,…,αN), Uh=(β1,…,βL-1). By substituting these two approximations into variational equation (3) and by choosing v=ϕi and V=nj, we obtain a discrete linear system
(7)A(σ)b=f,
where the data vector f=(0,I^)T with 0∈ℝN and I^=(I1-I2,I1-I3,…,I1-IL), the determined coefficient vector b=(uh,Uh)T, and the stiffness matrix
(8)A(σ)=(BCCTD),
with
(9)Bij=∫Ωσ∇ϕi·∇ϕjdx+∑l=1L1zl∫elϕiϕjdS,i,j=1,2,…,N,Cil=-1zl∫elϕidS,i=1,2,…,N,l=1,2,…,L-1,Dll=1zl∫eldS,l=1,2,…,L-1.
After solving the determined coefficients b=A(σ)-1f, we can use (6) to find the approximate potentials Uh 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 I and resulting voltage data U 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σ:I∈ℝ⋄L→U∈ℝ⋄L, which maps the applied electrode current data I to the computed electrode voltage data U or denoted by U(σ) to emphasize the dependence on the conductivity parameter as follows:
(10)U(σ)=RσI
and RσI=ΥDu∈ℝ⋄L, where ΥD:H~1(Ω)→ℝ⋄L denotes the Dirichlet trace operator, which projects the solution u of the forward problem to the electrode voltages U.Now we denote the solution operator byF(σ;I), such that
(11)U(σ)=F(σ;I),
where F relates the conductivity distribution σ with the corresponding virtual measurements U, often referred to as the parameter-to-observation map, also called the forward operator. Notably, F depends linearly on I but nonlinearly on σ. For a fixed current vector I, we can view F(σ;I) as a function of σ only. It is important to distinguish between the PDE operator A(σ) and the forward operator F(σ). Though both are associated with the model problem (1), their meaning is not the same. Indeed, A(σ) describes the forward problem, while F(σ) actually represents its solution.However, in practical applications, actual observations for EIT are typically potential differences between certain electrodes, denoted byUobs, 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:
(12)Uobs=F(σ)+e,
where e is the noise in the measurements.Thusthe inverse problem for EIT becomes that of finding the approximation conductivity distribution σ from the observed partial noisy knowledge Uobs, also called image reconstruction.


## 3. Homotopy Perturbation Inversion Method

One popular approach to deal withthe inverse problem (12) is the nonlinear least-square method
(13)minσ∥F(σ)-Uobs∥2.The corresponding Euler equation of (13) is given by
(14)F′(σ)T(F(σ)-Uobs)=0.By the homotopy perturbation technique, a fixed point homotopy functionH(v,p):R×[0,1]→R is constructed for (14), which satisfies
(15)H(v,p)=p[F′(v)T(F(v)-Uobs)]+(1-p)(v-σ0)=0,
where p∈[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
(16)H(v,0)=v-σ0=0,H(v,1)=F′(v)T(F(v)-Uobs)=0.
As the parameter p changed continuously from zero to one, the trivial problem H(v,0)=0 is continuously deformed to the original problem H(v,1)=0 that is, v changed from σ0 to the approximation solution σ. In topology, this is called deformation, and σ-σ0 and F′(σ)T(F(σ)-Uobs) are called homotopic.Next, according to the homotopy perturbation technique, we can use the parameterp as a small parameter and assume that the solution of (15) can be expressed as a power series in p,
(17)v=σ0+pσ1+p2σ2+p3σ3+⋯.
As p→1, the approximate solution of (14) is obtained with
(18)σ=limp→1v=σ0+σ1+σ2+σ3+⋯.We rewrite (15) with the Taylor expansion of F(σ) at σ0 and ignore the higher order term, which yields
(19)H(v,p)=p[F′(σ0)T(F(σ0)+F′(σ0)(v-σ0)-Uobs)]+(1-p)(v-σ0)=0.
Substituting (17) into (19), we will have
(20)p[×(pσ1+p2σ2+p3σ3+⋯)-Uobs)F′(σ0)T(×(pσ1+p2σ2+p3σ3+⋯)-Uobs)F(σ0)+F′(σ0)×(pσ1+p2σ2+p3σ3+⋯)-Uobs)×(pσ1+p2σ2+p3σ3+⋯)-Uobs)F′(σ0)T(×(pσ1+p2σ2+p3σ3+⋯)-Uobs)F(σ0)+F′(σ0)]+(1-p)(pσ1+p2σ2+p3σ3+⋯)=0.
Equating the coefficients for the different power of p, we get the following formulations:
(21)p1:F′(σ0)T(F(σ0)-Uobs)+σ1=0,pk+1:F′(σ0)TF′(σ0)σk+σk+1-σk=0,k=1,2,3,….
Then starting with an initial prior estimate σ0 and in turn solving (21), we can obtain
(22)σ1=-F′(σ0)T(F(σ0)-Uobs),σk+1=(I-F′(σ0)TF′(σ0))k(-F′(σ0)T(F(σ0)-Uobs)),k=1,2,3,….Thus by (18) we can get the approximate solution of (14); that is,
(23)σ=limp→1v=σ0+σ1+σ2+σ3+⋯=σ0+(-F′(σ0)T(F(σ0)-Uobs))+(I-F′(σ0)TF′(σ0))(-F′(σ0)T(F(σ0)-Uobs))+(I-F′(σ0)TF′(σ0))2(-F′(σ0)T(F(σ0)-Uobs))+⋯.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 σ0,
(24)σk+1=σk-αkF′(σk)T(F(σk)-Uobs),
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,
(25)σk+1=σk-αk(2I-F′(σk)TF′(σk))×(F′(σk)T(F(σk)-Uobs)),
where αk is an appropriately chosen step length. Since the iteration-varying step length usually shows better performance than the constant value, the step length αk for iterations (24) and (25) was chosen with the same criterion; that is, αk=1/λmax(JkTJk), where Jk is the Jacobian matrix of F(σ) calculated at the kth iteration σk and λmax(JkTJk) denotes the maximum eigenvalue of JkTJk. 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 dataUobs we employ the discrepancy principle as a stopping rule for HPIM, which determines the stopping index k*=k*(δ,Uobs) by
(26)∥Uobs-F(σk*)∥≤τδ<∥Uobs-F(σk)∥,0≤k<k*,
for some sufficiently large τ>0. 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 Figure1, 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.The forward and inverse models. (a) the forward model with 1968 triangulation elements; (b) the inverse model with 492 elements.
(a)
(b)Figure2 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.Three simulation phantoms.
(a)
(b)
(c)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:(27)DE=∥F(σ†)-Uobs∥,RE=∥σ†-σ*∥∥σ*∥,
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. Figures3(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 Figures3(a), 3(b), 4(a), 4(b), 5(a), and 5(b).
Iterations
Figure3(a)
Figure3(b)
Figure4(a)
Figure4(b)
Figure5(a)
Figure5(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.1713Reconstruction 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)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)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.
(a)
(b)
(c)


### 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(28)Unoise=Usyn+ε·n,
where ε is the noise level and n 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 k* and T 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τ=20 at four different noise levels.
ε
HPIM
CLIM
k
*
DE
RE
T (s)
k
*
DE
RE
T (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.6608Reconstruction results by HPIM for three phantoms with different noise levels. (a)ε=5%; (b) ε=3%; (c) ε=1%; (d) ε=0.5%.
(a)
(b)
(c)
(d)


## 4.1. Noise-Free Case

The numerical performances of HPIM and CLIM were first tested by the noise-free data. Figures3(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 Figures3(a), 3(b), 4(a), 4(b), 5(a), and 5(b).
Iterations
Figure3(a)
Figure3(b)
Figure4(a)
Figure4(b)
Figure5(a)
Figure5(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.1713Reconstruction 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)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)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.
(a)
(b)
(c)


## 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(28)Unoise=Usyn+ε·n,
where ε is the noise level and n 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 k* and T 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τ=20 at four different noise levels.
ε
HPIM
CLIM
k
*
DE
RE
T (s)
k
*
DE
RE
T (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.6608Reconstruction results by HPIM for three phantoms with different noise levels. (a)ε=5%; (b) ε=3%; (c) ε=1%; (d) ε=0.5%.
(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.



---

*Source: 454706-2013-12-15.xml*