Improving Image Quality in Electrical Impedance Tomography (EIT) Using Projection Error Propagation-Based Regularization (PEPR) Technique: A Simulation Study

Electrical Impedance Tomography is a non linear inverse problem [21] in which the electrical conductivity distribution of a closed domain (Ω) in an object under test called phantom, [28, 29] is reconstructed from the surface potential at the boundary (∂Ω) developed by injecting a sinusoidal current signal [30]. A low frequency (10kHz-1MHz) and low magnitude constant current is injected through an array of 16 surface electrodes [31] surrounding the domain to be imaged using different electrode switching protocols [32] and the boundary potentials are measured [33] although, the boundary currents can also be measured for an applied voltage signal in Applied Potential Tomography [34]. The voltage data collected by the data acquisition system is then processed by an image reconstruction algorithm in a PC.

To calculate the nodal potential for a known conductivity a relationship can be established between the electrical conductivity (σ) and spatial potential (Φ). Electro-dynamics of EIT is governed by a nonlinear partial differential equation, called the Governing Equation [2, 3, 4] of Electrical Impedance Tomography, which is given by,

A relation can be obtained between the voltage measurements made on the boundary (∂Ω) and the domain conductivity can be found [35, 36] as,

Where σ is elemental conductivity values, Ф is the vector of nodal potential and K is the transformation matrix constructed from the elemental conductivities and nodal coordinates.

If K and I are known, Eq.-2 can be solved numerically using the finite element method (FEM) [37] to calculate the nodal potentials of the domain for the known conductivity (σ). It is known as the “forward problem”. Using the Gauss-Newton method [38] applied on EIT, the update vector of σ [35, 36] can be expressed as:

Where, Q is a function of the Jacobian matrix (J) [35, 36] and regularization parameters [35, 36]. Φ_d is the mismatch vector between calculated boundary potential (V_c) and measured boundary potential (V_m). That means if the matrix Q and the surface potentials Φ are known then the elemental conductivity (σ_e) can be mapped. This is known as the “inverse problem” which is discussed in the next section. Using the Modified Newton Raphson (MNR) iterative technique [35, 36], a suitably assumed conductivity vector (initial guess), [σ_o], is modified to [σ_o + Δσ ] for achieving a specified error limit in the calculated and measured voltage ([Δσ] denotes the conductivity update).

Electrical conductivity imaging is a highly nonlinear and ill-posed inverse problem. The response matrix [J^TJ] is a singular matrix. Hence in EIT, a minimization algorithm [35, 36, 38, 39] is used to obtain the approximate solution of the ill-posed inverse problem. In the minimization algorithm, the objective function formed by the difference between the experimental measurement data (V_m) and the computationally predicted data (V_c) is minimized. Generally in inverse problems a least square solution [35, 36] of a minimized object function (s) [38, 39] obtained from the calculated voltage data and the measured voltage data is searched by a Gauss-Newton method based numerical approximation algorithm (explained in the next sub-section) called the inverse solver [40].

If V_m is the measured voltage matrix and f is a function mapping an E-dimensional (E is the number of element in the FEM mesh) impedance distribution into a set of M (number of the measured data available) approximate measured voltage, then the Gauss-Newton algorithm [35, 36, 38, 39] tries to find a least square solution of the minimized object function s defined as [38, 39, 40, 41, 42, 43]:

Where, s_r is the constrained least-square error of the reconstructions, G is the regularization operator and λ is a positive scalar and called as the regularization coefficient.

Now, differentiating Eq.-5 w.r.t. σ, it reduces to:

Where the term f′ = J is known as the Jacobian matrix of size g × h and which is defined by [36, 39]:

Where,

e = 1, 2 … E [E = number of elements in the FEM mesh], h = 1, 2…M [M = (number of data measured per current projections (d)) × (number of current projections (p))]

By the inherent ill-posed nature of EIT, the [f ^′]^Tmatrix in Eq.-6 is always ill-conditioned [22], and hence small measurement errors will make the solution of Eq.-6 change greatly which is made well posed by the regularization term incorporated. Differentiating Eq.-6 w.r.t. σ, it reduces to:

By the Gauss Newton (GN) method,

Neglecting higher terms, the update vector reduces to:

Replacing G^TG by I (Identity matrix) Eq.-10 reduces to:

In general, for k^th iteration (k is a positive integer), the conductivity update vector of Eq.-11 is reduced to:

Where J_k and ( ΔV )_k are the Jacobian and voltage difference matrix respectively at the k^th iteration.

Thus the Gauss-Newton method based inverse solver algorithm gives a regularized solution of the conductivity distribution for the k^th iteration as:

In this paper a PEPR method is proposed to improve the reconstructed image quality in EIT and the resistivity reconstruction is studied in EIDORS. Generally, the choice of an appropriate regularization can be evaluated or found empirically [25] and the regularization parameter is related to an objective function [38, 39]. In our study, the projection error is utilized as the objective function; hence the regularization parameter would be expected to be related to the projection error. Furthermore, the value of regularization is required to vary with projection error. That is to say, a greater regularization value is needed for a larger projection error during iterations. If the projection error is very low, only a small regularization value is needed to regulate the ill-posed process [25]. Based on these considerations, we define the adaptive regularization parameter λ as follows:

Where, Ψ is a const ( taken as 0.01), ΔΦ is the projection error which is defined as the L-2 norm of the difference between the calculated data and the measured data. Hence the regularization parameter λ in the PEPR method is given by:

In formula (15), the regularization parameter λvaries in the range from 0.01/3 to 0.01/2. It includes two primary considerations: (1) Due to the ill-posed characteristic of the inverse problem, the regularization should not be reduced with iterations to a too small value. (2) The regularization value should be lower than the maximum value of diagonal elements in the matrix JJ^T (calculated by the MATLAB code: max(max(JJ^T))). In the PEPR method the regularization parameter is set as a function of the projection error which is produced by the mismatch between the calculated and measured data. In the first iteration, the regularization parameter is calculated with the projection error developed by the boundary data estimated by the forward solver for an initial guessed resistivity. According to the resistivity update vector calculated, the regularization parameter is also modified in all the other iterations. Hence at each iteration, the projection error varies according to the misfit between the model predicted data and experimental data. The variation of mismatch data is integrated with the response matrix and the reconstruction is carried out.

The J^TJ matrix is formed by the first derivative (J) of the forward solution (FS) and I is the approximation of the Hessian [39] which is, actually, the second derivative of the FS. Generally, the order of magnitude of Hessian is less than the J^TJ due to the higher order derivative.

Hence, in Eq.-13, I acts as the quadratic term which is formed with the maximum value of J^TJ. The magnitude of the diagonal elements of λI should be less than the maximum value of J^TJ. The identity matrix I is formed with the diagonal values equal to the maximum value of the matrix J^TJ. In the LMR method the λ is taken as 10^-1 in the first iteration and it is then modified as λ/k in the modified EIDORS, where k is the number of iterations. For further analysis the resistivity reconstruction is also conducted with a range of regularization parameters.

It is observed that the projection error in LMR and PEPR becomes minimum at λ = 0.1 and Ψ = 0.01 respectively. That is why in the first iteration in LMR λ is taken as 0 .1 and in PEPR λ is calculated with Ψ = 0 .01. Furthermore, for better understanding the regularization effects of both the methods, λ is calculated in PEPR technique (using the Eq.-14) with Ψ = 0.01. The result is compared with the reconstruction obtained in LMR in which the iteration starts with λ = 0.01 and then it is decreased by a factor of 10$$\sqrt{10}$$for the other iterations [25].

Resistivity images are reconstructed using simulated boundary data in EIDORS with the PEPR and LMR method and the results are compared. Resistivity images are reconstructed with simulated data for different inhomogeneity geometries. Circular objects (diameter = 60 mm, resistivity = 33 Ωm) with a homogeneous background medium (diameter (D) = 150 mm, resistivity = 2.5 Ωm) are simulated for boundary data generation. Resistivity images are also reconstructed from the boundary data combined with random noise of different percentages. Noisy data are used for reconstruction with the PEPR and LMR methods and the images are compared.

To analyze the proposed method, the normalized projection error (error due to the voltage mismatch), E_V, is calculated in each iteration as:

The normalized solution error norm (E_ρ) is also calculated in each iteration as:

A contrast to noise ratio (CNR) [44, 45] is calculated for the reconstructed images in this work to evaluate the reconstructed images with different regularization techniques. CNR is defined as the ratio of the difference between the average inhomogeneity resistivity (IR_Mean) and the average background resistivity (BR_Mean) divided by the weighted average of the standard deviations in the IR (SD_IR) and BR (SD_BR):

Where ω_I is the fraction of the area of the region of interest with respect to the area of the whole image; ω_B is defined as ω_B=1-ω_I. IR_Mean and BR_Mean are the mean values of the inhomogeneity and the background regions in the reconstructed images.

Percentage of contrast recovery (PCR) [46] is calculated for the images reconstructed by all the regularization technique to compare the reconstruction accuracy. PCR in EIT is defined as the difference between the averaged resistivity within the reconstructed image (IR_Mean) and the reconstructed background (BR_Mean) divided by the difference between the original resistivity of the inhomogeneity (IR_Original) and the background (BR_Original). Hence mathematically the PCR is obtained by the equation:

Coefficient of Contrast (COC) in EIT is defined as the ratio of the mean inhomogeneity resistivity (IR_mean) to mean background resistivity (BR_mean) though in some literature this ratio is termed as the contrast recovery [47]:

CNR, PCR, COC, IR_Mean and IR_Max (maximum values of the reconstructed inhomogeneity resistivity) are calculated for all the images and are compared to assess the reconstructed image quality with PEPR. To evaluate the reconstructed images the elemental resistivity along the phantom diameter (D), connecting the centre of the reconstructed object and the centre of the phantom, is plotted against the length of the phantom diameter. This resistivity plot is termed as the diametric resistivity plot (DRP). The normalized projection error (error due to the voltage mismatch) and the normalized solution error norm are studied for different iterations.

The reconstruction of the simulated phantom with a circular object near electrode 1 (Fig.-2a) shows that the CNR of the resistivity image reconstructed with LMR (Fig.-2b) is 3.08 whereas the CNR of the resistivity image with PEPR (Fig.-2c) method (Table-1) is 3.55. For the same simulated phantom, the PCR of the resistivity image with LMR technique is 32.52 whereas the PCR with PEPR is 46.57 (Table-1). It is also noticed that the COC of the resistivity image with LMR technique is 2.05 whereas the COC with PEPR is 2.46 (Table-1). It is observed that the DRP of the resistivity image with PEPR (Fig.-2d) is more similar to the DRP of the original object.

Fig.2

For the object near electrode 3 (Fig.-3a), the reconstructed images (Fig.-3b-3c) show that the CNR of the resistivity image with LMR technique is 2.99. On the other hand, CNR with PEPR is 3.51 (Table-2). The PCR of the resistivity image with LMR technique is 31.33 whereas the PCR with PEPR is 46.26 (Table-2). It is also noticed that the COC of the resistivity image with LMR is 2.01 but with PEPR, it is 2.45 (Table-2). Result show that the DRP of the resistivity image with PEPR (Fig.-3d) is more similar to the DRP of the original object.

Fig.3

It is observed that for the simulated phantom with a circular object near electrode 5 (Fig.-4a), the CNR of the image reconstructed with LMR (Fig.-4b) is 2.92 whereas the CNR of the resistivity image with PEPR (Fig.-4c) is 3.46 (Table-3). It is also noticed that the PCR of the image with LMR is 30.32 whereas the PCR with PEPR is 44.74. Results also show that the COC with LMR is 1.98 (Table-3) whereas with PEPR technique it is 2.40. It is observed that the DRP of the resistivity image with PEPR (Fig.-4d) is more similar to the DRP of the original object.

Fig.-4

Resistivity images are reconstructed from the boundary data with 25 % random noise for the object near electrode 1, 3 and 5. Images reconstructed from the boundary data with added noise for the object near electrode 1 show that the CNR of the resistivity image with LMR technique (Fig.-5a) is 1.81. On the other hand, CNR of the image with PEPR (Fig.-5b) is 3.17 (Table-4). The PCR of the resistivity image with LMR technique is 27.29 whereas the PCR with PEPR is 63.79 (Table-4). It is also noticed that the COC of the resistivity image with LMR is 1.85 but with PEPR it is 2.90 (Table-4). Result shows that the DRP of the resistivity image with PEPR (Fig.-5c) is more similar to the DRP of the original object.

Fig.-5

Resistivity images reconstructed from the noisy boundary data of the phantom with object near electrode 3 show that the CNR of the image reconstructed with LMR technique (Fig.-6a) is 1.60. On the other hand, CNR with PEPR (Fig.-6b) is 1.89 (Table-5). The PCR of the resistivity image with LMR technique is 23.61 whereas the PCR with PEPR is 35.68 (Table-5). It is also noticed that the COC of the reconstructed image with LMR is 1.73 but with PEPR it is 2.12 (Table-5). It is observed that the DRP of the resistivity image with PEPR (Fig.-6c) is more similar to the DRP of the original object.

Fig.-6

It is also observed that for the resistivity images reconstructed from the noisy data for the phantom with object near electrode 5, the CNR of the image reconstructed with LMR technique (Fig.-7a) is 1.24. On the other hand, CNR with PEPR (Fig.-7b) is 1.48 (Table-6). The PCR of the resistivity image with LMR technique is 17.59 whereas the PCR with PEPR is 25.91 (Table-6). It is also noticed that the COC of the resistivity image with LMR is 1.54 but with PEPR, it is 1.78 (Table-6). Result shows the DRP of the resistivity image with PEPR (Fig.-7c) is more similar to the DRP of the original object.

Fig.-7

In the above study, for the LMR method the λ is taken as 10^-1 in the first iteration and it is then modified as λ/k in the modified EIDORS, where k is the number of iterations. For further analysis the resistivity reconstruction is also conducted with a range of regularization parameters. Regularization effect, normalized projection error (error due to the voltage mismatch) and the normalized solution error norm are studied for different iterations. Result show that the resistivity image with PEPR (Ψ = 0.01) is found better compared to the image obtained for LMR with a λ = 0.01 at first iteration and then it is decreased by a factor of 10$$\sqrt{10}$$for the other iterations [25]. It is observed that CNR of the image reconstructed with LMR technique (Fig.-8a) is 2.95. On the other hand, CNR with PEPR (Fig.-8c) is 3.26 (Table-7). Due to the over-estimate of the inhomogeneity resistivity (Fig.-8b) the PCR of the resistivity image with LMR technique is very high compared to PEPR. For the similar reason COC is also found higher than the PEPR method. But IR_Max and IR_Mean are both suitable in PEPR method whereas, in LMR technique, they are very absurd compared to the original value (Table-7).

Fig.-8

Resistivity imaging with LMR ( λ = 0 .01) and PEPR ( Ψ = 0.01) techniques is studied for a number of iterations (result for first 12 iterations is presented). It is observed that resistivity images in LMR method (Fig.-9), the reconstruction diverges gradually after second iteration and become unstable with a continuously over-estimated resistivity (Fig.-10). On the other hand, the resistivity reconstruction rapidly converges in the PEPR method (Fig.-11) and gets stable after few iterations with a proper resistivity reconstruction (Fig.-12).

Fig.-9

Fig.-10

Fig.-11

Fig.-12

It is observed that the IR_Max of the resistivity images with PEPR technique is more stable and closer to the original value, except for the first iteration (Fig.-13a). On the other hand, in LMR, IR_Max is larger than the original value (Fig.-13a) except in first iteration (IR_Max < IR_Original) and the second iteration (IR_Max ~ IR_Original). Result show that, for PEPR, the standard deviation of the IR_Max during first twelve iterations is 1.79 whereas it is 17.59 in LMR method (Fig.-13a).

Fig.-13

Result show that IR_Mean of the resistivity images with the PEPR method is more stable whereas it is comparatively largely variable in the LMR method (Fig.-13b). Result show that, for PEPR, the standard deviation of the IR_Mean during the first twelve iterations is 1.79 Ωm whereas it is 8.21 Ωm in LMR method (Fig.-13b).

It is observed that the projection error (E_V) in LMR method is comparatively large and it is gradually increasing after third iteration (Fig.-14a). On the other hand, in PEPR, E_V is comparatively low and almost constant after second iteration (Fig.-114a). Result show that, in the PEPR method, the normalized solution error norm (E_ρ) is less and varies from 0.74 to 0.77 (Fig.-14b). But, in the LMR method, E_ρ is comparatively large and varies from 0.73 to 1.31 (Fig.-14d).

Fig.-14

The projection errors (E_V) are calculated for different values of λ in LMR (Fig.-15a) and different values of Ψ in the PEPR method (Fig.-15b). It is observed that, for the LMR method, the E_V becomes minimum at the third iteration for λ = 1.0 and 0.1 whereas for λ = 0.01, 0.001 and 0.0001, the E_V becomes minimum at the fourth iteration (Fig.-15a). On the other hand, result show that (Fig.-15b), for the PEPR method, the E_V becomes minimum at the third iteration for Ψ = 1.0 and 0.1, 0.01 and 0.001 whereas for Ψ = 0.0001, the E_V becomes minimum at the fourth iteration (Fig.-15a). Hence it is observed that the optimum reconstruction (i.e. the E_V is minimum) occurs in PEPR with Ψ = 0.01 where as the optimum reconstruction occurs in LMR with λ = 0.1 (Table-8). It is also observed that the projection errors are comparatively less (Fig.-15b) in the PEPR method for all the values of Ψ. Result show that the projection errors are comparatively more stable in PEPR and they also become almost constant (Fig.-15b) after 3_rd iteration except for very low Ψ (Ψ= 0.0001).

Fig.-15