An image reconstruction algorithm of EIT based on pulmonary prior information

Abstract

Using a CT scan of the pulmonary tissue, a human pulmonary model is established combined with the structure property of the human lung tissue using the software COMSOL. Combined with the conductivity contribution information of the human tissue and organ, an image reconstruction method of electrical impedance tomography based on pulmonary prior information is proposed using the conjugate gradient method. Simulation results show that the uniformity index of sensitivity distribution of the pulmonary model is 15.568, which is significantly reduced compared with 34.218 based on the round field. The proposed algorithm improves the uniformity of the sensing field, the image resolution of the conductivity distribution of pulmonary tissue and the quality of the reconstruction image based on pulmonary prior information.

Keywords

electrical impedance tomography (EIT) / prior information / pulmonary model of human / image reconstruction / COMSOL

Introduction

As one of the most important research fields in bio-medical engineering, electrical impedance tomography (EIT) is a non-invasive, functional imaging technique. EIT is a technique developed in the last decade after the morphology and structure imaging methods. Medical research shows that not only does different impedance exists in different tissues (organs) of human bodies, but also impedance would be changed because of pathological or physiological reasons. Consequently, plentiful pathological and physiological information can be obtained from impedance of tissues. For the EIT technique, conducting electrodes are attached to the skin of the subject and small alternating currents (less than 5 mA for safety) are applied to some or all of the electrodes. Then, the resulting electrical potentials are measured. After repeating the process for other configurations of applied current, data of potentials, which are relevant with pathological or physiological status, are delivered to compute for image reconstruction. A 2D/3D image of the conductivity or permittivity of part of the body is calculated by the image reconstruction algorithm [1,2].

The image reconstruction method is a crucial technique in the EIT system. From a mathematical point of view, this reconstruction problem is an ill-posed nonlinear inverse problem, due to the soft-field property of EIT and limitations of measured data. To improve the quality of the reconstructed image, one feasible method is to reduce the property of ill-posedness by combining with some prior information. For the EIT system used for pulmonary examination, prior information mainly contains the conductivity contribution information of human tissue and organ and the structure of pulmonary tissue. After doing basic research on the structure and conductivity of pulmonary tissues and organs, this article proposed a new reconstruction image method based on pulmonary prior information.

Reconstruction model

The general Laplace equation reflects the relationship of current and voltage in the physical models of EIT. Considering both the shunting effect of the electrodes and the contact impedances between the electrodes and tissue, the complete electrode model was used in this article. The physical model and boundary conditions can be written in the form [3-5]:

(1)\begin{array}{ll}\nabla \cdot (\sigma \nabla u)=\mathrm{0},& x\in \mathrm{\Omega },\end{array}

(2)\begin{array}{ll}u+z_l\sigma {\displaystyle \frac{\partial u}{\partial n}}=U_{\mathrm{0}},& x\in e_l,l=\mathrm{1},\mathrm{2},\mathrm{...},L,\end{array}

(3)\begin{array}{ll}{\displaystyle \underset{e_l}{\int }\sigma {\displaystyle \frac{\partial u}{\partial n}}}\mathrm{d}S=I_l,& x\in e_l,l=\mathrm{1},\mathrm{2},\mathrm{...},L,\end{array}

(4)\begin{array}{ll}\sigma {\displaystyle \frac{\partial u}{\partial n}}=\mathrm{0},& x\in \partial \mathrm{\Omega }\backslash {\displaystyle \underset{l=\mathrm{1}}{\overset{L}{\cup }}{e}_l,}\end{array}

where \sigma is conductivity distribution of the medium, u is the potential in the domain \mathrm{\Omega }, U_{\mathrm{0}} is the measured potentials of the boundary \partial \mathrm{\Omega }, I_l is the injected current of the lth electrode, z_l is the contact impedance between the lth electrode and inner medium.

Normally, we choose a round field for the forward modeling, as shown in Fig. 1. However, in the pulmonary EIT system, the round field does not reflect and simulate the structure property of human lungs. Accordingly, the forward modeling based on the CT scan picture of the human chest was necessary, as shown in Fig. 2. The human pulmonary model was built up according to the chest structure by the finite element method (FEM) software COMSOL, as shown in Fig. 3. After determining the conductivity of every tissue, the forward problem was calculated by dividing the whole domain into FEM elements (small triangles). The proposed model reflected the structure property of the human chest and the inhomogeneous distribution of conductivity precisely.

Prior information of the structure and conductivity of tissues and organs of the human chest were integrated in the model. The conductivity of different tissues and organs are listed in Table 1 [6].

Forward problem and Jacobian matrix based on prior information

After building up the human pulmonary model, 30354 FEM elements were obtained by the automatic triangulation function of COMSOL, as shown in Fig. 4. Large numbers of elements improved the precision of the forward calculation, but were not suitable for the inverse problem.

Therefore, a square mesh (see Fig. 5) was necessary for solving the Jacobian matrix and the calculation of the inverse problem.

Procedures for the Jacobian calculation are as follows [7].

From Eq. (1), select a continuous scalar function \omega for any \omega,

(5){\displaystyle \underset{\mathrm{\Omega }}{\int }\sigma \nabla u\cdot \nabla \omega \mathrm{d}V={\displaystyle \underset{\partial \mathrm{\Omega }}{\int }\omega \sigma {\displaystyle \frac{\partial u}{\partial n}}}}\mathrm{d}S\mathrm{.}

Here, dVand \mathrm{d}Sare the differentiation of volume and surface area. In particular, for \omega ＝u, we have the power conservation formula

(6){\displaystyle \underset{\mathrm{\Omega }}{\int }\sigma {\left|\nabla u\right|}^{\mathrm{2}}}\mathrm{d}V={\displaystyle \underset{\partial \mathrm{\Omega }}{\int }u\sigma {\displaystyle \frac{\partial u}{\partial n}}}\mathrm{d}S.

Substituting Eq. (2) into the right-hand side of Eq. (6), we have

(7){\displaystyle \underset{\mathrm{\Omega }}{\int }\sigma {\left|\nabla u\right|}^{\mathrm{2}}}\mathrm{d}V+{\displaystyle \sum _l{\displaystyle \underset{e_l}{\int }{z}_l{\left(\sigma {\displaystyle \frac{\partial u}{\partial n}}\right)}^{\mathrm{2}}\mathrm{d}S}}={\displaystyle \sum _l{V}_l{I}_l}.

Equation (7) indicates that the power input is dissipated either in the domain \mathrm{\Omega } or the contact impedance layer under the electrodes. Perturbations are taken into Eq. (7): \sigma \to \sigma +\delta \sigma ,u\to u+\delta u,V_l\to V_l+\delta V_l, and the current in each electrode I_l stays constant. Ignoring second-order terms, we have

(8){\displaystyle \underset{\mathrm{\Omega }}{\int }\delta \sigma {\left|\nabla u\right|}^{\mathrm{2}}}\mathrm{d}V+\mathrm{2}{\displaystyle \underset{\mathrm{\Omega }}{\int }\sigma \nabla u\cdot \nabla \delta u}\mathrm{d}V+\mathrm{2}{\displaystyle \sum _l{\displaystyle \underset{e_l}{\int }{z}_l\left(\sigma {\displaystyle \frac{\partial u}{\partial n}}\right)\delta \left(\sigma {\displaystyle \frac{\partial u}{\partial n}}\right)\mathrm{d}S={\displaystyle \sum _l{I}_l\delta V_l}}}.

From Eq. (2), on the lth electrode e_l, we have

(9)\delta \left(\sigma {\displaystyle \frac{\partial u}{\partial n}}\right)={\displaystyle \frac{\mathrm{1}}{z_l}}(\delta V_l-\delta u).

From Eqs. (5) and (9), letting \omega ＝\delta u, we get

(10){\displaystyle \underset{\mathrm{\Omega }}{\int }\delta \sigma {\left|\nabla u\right|}^{\mathrm{2}}\mathrm{d}V+\mathrm{2}{\displaystyle \underset{\partial \mathrm{\Omega }}{\int }\delta u\sigma {\displaystyle \frac{\partial u}{\partial n}}\mathrm{d}S-\mathrm{2}{\displaystyle \sum _l{z}_l{\displaystyle \underset{e_l}{\int }{\displaystyle \frac{\delta u}{z_l}}\sigma {\displaystyle \frac{\partial u}{\partial n}}\mathrm{d}S}}}}+\mathrm{2}{\displaystyle \sum _l\delta V_l{\displaystyle \underset{e_l}{\int }\sigma {\displaystyle \frac{\partial u}{\partial n}}\mathrm{d}S={\displaystyle \sum _l{I}_l\delta V_l}}}.

The second and third terms on the left-hand side of Eq. (10) can be canceled, and the fourth is simplified as \mathrm{2}{\displaystyle \sum _l{I}_l\delta V_l}, therefore,

(11){\displaystyle \sum _l{I}_l\delta V_l=-{\displaystyle \underset{\mathrm{\Omega }}{\int }\delta \sigma {\left|\nabla u\right|}^{\mathrm{2}}\mathrm{d}V}}.

For the current injection pattern, electrode currents can be simplified as a vector \mathit{I}=(I_{\mathrm{1}},I_{\mathrm{2}},\mathrm{...},I_L). To emphasize that the potential distribution depends on the current excitation pattern of the different position, it is denoted by u(\mathit{I}). For instance, the potential distributions for the mth and the dth excitation patterns are represented by u({\mathit{I}}^m)and u({\mathit{I}}^d), respectively. We assume that the injected current is unit current. From Eq. (11), we can obtain the potential difference for two different excitation patterns:

(12)\delta V_{dm}=-{\displaystyle \underset{\mathrm{\Omega }}{\int }\delta \sigma \nabla u({\mathit{I}}^d)\cdot \nabla u({\mathit{I}}^m)\mathrm{d}V}.

The four-electrode data collection method (the so-called neighboring method) was adopted in this article, as shown in Fig. 6. In this method, the current was injected through two adjacent electrodes and the resulting voltages were measured from all other pairs of electrodes. This was repeated for all the electrodes, and the resulting voltages from different adjacent electrodes were obtained.

For the constant current I, the resulting voltage difference V_{ij} of the ith and jth is

(13)\delta V_{ij}=-{\displaystyle \frac{\mathrm{1}}{I}}{\displaystyle \underset{\mathrm{\Omega }}{\int }\delta \sigma \nabla u_i\cdot \nabla u_j\mathrm{d}V}.

To calculate the Jacobian matrix, the conductivity should be discretized by making the conductivity piecewise constant on polyhedral domains, such as voxels. Taking \delta \sigma as the characteristic function of the kth voxel, we have the Jacobian matrix for a fixed current pattern, when -\nabla u=E [7,8],

(14)\mathit{J}={\displaystyle \frac{\partial V_{dm}}{\partial {\sigma }_k}}=-{\displaystyle \underset{\mathrm{v}\mathrm{o}\mathrm{x}\mathrm{e}\mathrm{l}k}{\int }\mathit{E}({\mathit{I}}^d)\cdot \mathit{E}({\mathit{I}}^m)\mathrm{d}V}.

Comparison of sensitivity distribution

Based on the above analysis, the proposed model used for forward calculation has combined prior information, the chest structure property and the conductivity contribution of human tissues and organs. Figures 7 and 8 showed four typical sensitivity distributions of the pulmonary model and round model respectively.

In Fig. 7, (a)-(d) were three-dimensional sensitivity distribution images of measurement from electrode pairs of 2–3, 3–4, 5–6, 7–8 respectively when electrodes 1 and 2 were current excitation electrodes in the pulmonary model, while in Fig. 8, it was the same sensitivity distribution but for the round model.

From these figures, compared with the round model, we can conclude that the sensitivity distribution for the pulmonary model was more uniform, especially that the sensing field did not contain an obvious spike between two non-adjacent electrodes, that the influence regions were widened, and non-linearity was improved significantly.

To quantify the non-linearity of the sensing field, the uniformity index of the sensitivity distribution was defined as follows [9,10]:

(15){\mathit{J}}_{i,j}^{\mathrm{a}\mathrm{v}\mathrm{g}}={\displaystyle \frac{\mathrm{1}}{n}}{\displaystyle \sum _{e=\mathrm{1}}^n{\mathit{J}}_{i,j}(e)},

(16){\mathit{J}}_{i,j}^{\mathrm{d}\mathrm{e}\mathrm{v}}={\left\{{\displaystyle \frac{\mathrm{1}}{n-\mathrm{1}}}{\displaystyle \sum _{e=\mathrm{1}}^n[{\mathit{J}}_{i,j}(e)}-{\mathit{J}}_{i,j}^{\mathrm{a}\mathrm{v}\mathrm{g}}{]}^{\mathrm{2}}\right\}}^{\mathrm{1}/\mathrm{2}},

(17)p_{i,j}={\displaystyle \frac{{\mathit{J}}_{i,j}^{\mathrm{d}\mathrm{e}\mathrm{v}}}{{\mathit{J}}_{i,j}^{\mathrm{a}\mathrm{v}\mathrm{g}}}},

(18)P={\displaystyle \sum _{i,j}\left|p_{i,j}\right|},

where BoldItalic is the sensitivity matrix or Jacobian matrix, P is the uniformity index of the sensing field. Obviously, the smaller the value of P, the more uniform the sensing field will be.

The results show that the uniformity index of the sensitivity distribution of the pulmonary model is 15.568, which is significantly reduced compared with 34.218 based on the round field. Apparently, the sensing field of the pulmonary model is more uniform than that of the round model, and the non-linearity is improved considerably.

Image reconstruction

Image reconstruction for EIT is the calculation of inner conductivity distribution \sigma from measured data of the boundary [11]. Conjugate gradient method was used in this article for the inverse problem calculation. After some linearization procedures, a linear approximation of an EIT model takes the following form:

(19)\mathit{A}\mathit{x}=\mathit{b},

where BoldItalic is the sensitivity matrix (or Jacobian matrix) of the linear system; BoldItalic is the boundary voltage measurements; BoldItalic is the change in conductivity.

The solution to the inverse problem is calculated by minimizing the functional f(\mathit{x}) by the conjugate gradient method.

(20)f(\mathit{x})={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{x}}^{\mathrm{T}}({\mathit{A}}^{\mathrm{T}}\mathit{A})\mathit{x}-{({\mathit{A}}^{\mathrm{T}}\mathit{b})}^{\mathrm{T}}\mathit{x}\mathrm{.}

The minimization of the functional f(\mathit{x}) means that the gradient of f(\mathit{x}) equals zero:

(21)\nabla f(\mathit{x})=\mathit{A}\mathit{x}-\mathit{b}=\mathrm{0.}

Therefore, the result of equation \mathit{A}\mathit{x}=\mathit{b} can be obtained by minimizing functional f(\mathit{x}).

The conjugate gradient method is suitable for calculating the functions of which the coefficient matrix is symmetrical positive definite (SPD). For the EIT system, coefficient matrix (Jacobian matrix) is not SPD matrix generally. Consequently, the original equation should be converted to the following form:

(22){\mathit{A}}^{\mathrm{T}}\mathit{A}\mathit{x}={\mathit{A}}^{\mathrm{T}}\mathit{b}\mathrm{.}

Since BoldItalic is not column full rank matrix, thus {\mathit{A}}^{\mathrm{T}}\mathit{A} is not a symmetrical positive definite matrix. After the regularization processing, Eq. (22) can be converted to the following form:

(23)({\mathit{A}}^{\mathrm{T}}\mathit{A}+\lambda \mathit{I})\mathit{x}={\mathit{A}}^{\mathrm{T}}\mathit{b}.

Till now, the conjugate gradient method can be used to calculate the solution by reason that the coefficient matrix has been translated to an SPD matrix.

Results

By using the pulmonary model, the simulation result of image reconstruction, as shown in Fig. 9, can clearly reveal the structure and conductivity properties of tissues and organs.

After conducting human lung tests using the EIT system (shown in Fig. 10), the measured data were obtained and saved for the comparison of image reconstruction of the pulmonary model and round model. Figure 11 shows the reconstructed image based on the round model, while Fig. 12 displays the result of the pulmonary model using the real data obtained from the human lung test.

From Figs.11 and 12, the reconstructed result was much better for the pulmonary model than for the round model, especially for the reconstructed image resolution.

Conclusions

Using a CT scan of the pulmonary structure, a human pulmonary model is established combined with the structure property of human lung tissue using the software COMSOL. Combined with the conductivity contribution information of human tissues and organs, an image reconstruction method of electrical impedance tomography based on pulmonary prior information is proposed using conjugate gradient method. Simulation results show that the uniformity index of the sensitivity distribution is reduced from 34.218 of the round field to 15.568 of the pulmonary model. The proposed algorithm improves the uniformity of the sensing field, the image resolution of the conductivity distribution of pulmonary tissue and the quality of reconstruction image based on pulmonary prior information.

References

[1] Wang H X, Wang C Y. Pulmonary electrical impedance tomography system. Chinese Medical Equipment Journal, 2006, 27(4): 3-4 (in Chinese)

[2] Wang H, Gao J B, Luo J P. Review of electrical impedance tomography. Beijing Biomedical Engineering, 2006, 25(2): 209-212 (in Chinese)

[3] Bi D X. Electro-magnetic Field Theory. Beijing: Publishing House of Electronics Industry, 1985 (in Chinese)

[4] Polydorides N, McCann H. Electrode configurations for improved spatial resolution in electrical impedance tomography. Measurement Science and Technology, 2002, 13(12): 1862-1870

[5] Feng C Z, Ma X K. An Introduction to Engineering Electromagnetic Fields. Beijing: Higher Education Press, 2000 (in Chinese)

[6] Dehghani H, Barber D C, Basarab-Horwath I. Incorporating a priori anatomical information into image reconstruction in electrical impedance tomography. Physiological Measurement, 1999, 20(1): 87-102

[7] Polydorides N, Lionheart W R B. 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. 2002, 13(12): 1871-1883

[8] Hu L, Wang H X, Zhao B, Yang W Q. A hybrid reconstruction algorithm for electrical impedance tomography. Measurement Science and Technology, 2007, 18(3): 813-818

[9] Wang H X, Zhang L F, Zhu X M. Optimum design of array electrode for ECT system. Journal of Tianjin University, 2003, 36(3): 307-310 (in Chinese)

[10] Xie C G, Plaskowski A, Beck M S. 8-electrode capacitance system for two-component flow identification, Part I: tomographic flow imaging. IEE Proceedings A – Science, Measurement and Technology, 1989, 136(4): 173-183

[11] Kauppinen P, Hyttinen J, Malmivuo J. Sensitivity distribution simulations of impedance tomography electrode combinations. International Journal of Bioelectromagnetism, 2005, 7(1): 344-347

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg