Bo Gong , Benjamin Schullcke , Sabine Krueger-Ziolek and Knut Moeller

A clustering based dual model framework for EIT imaging: first experimental results

De Gruyter | Published online: September 12, 2015

Abstract

Electrical Impedance tomography (EIT) imaging suffers greatly from the illposedness of the corresponding inverse problem. This is mainly caused by the high degree of freedom and the relatively large noise. One attempt to circumvent these difficulties is to use dual models. This article introduces a clustering based non-uniform dual model construction. With this framework, finite elements are grouped to reduce the complexity in inverse computations. The simulation and experiment results indicated that the k-means clustering method did not only preserve the sharp variations over conductivity mediums but also greatly filtered out artefacts found in the standard approach.

1 Introduction

Electrical Impedance Tomography is a radiation free imaging method that can monitor e.g. dynamics of volume distribution in the lungs [1, 2]. It attempts to reveal an inner conductivity distribution inside the human body by electrical data obtained via electrodes attached on the skin. By injecting an electrical current through some electrodes, the induced voltages can be recorded, allowing to reconstruct a conductivity image via the relation between current and voltage information at the boundary. In a common routine of lung EIT experiments, electrodes are located around a horizontal 2D plane across the human thorax. For each EIT frame, currents are successively injected into the human body through adjacent electrodes [3]. A conductivity image is then reconstructed by collecting the voltage measurements recorded from the remaining electrodes. We denote the number of voltage measurements in a frame by M.

A common reconstruction strategy is to seek the solution of a nonlinear optimization problem using Tikhonov regularization [4]. The mathematical formulation of such optimization problem is

(1) s = argmin s V m F ( s ) 2 2 + α 2 s s r e f 2 2 ,

where s denotes a conductivity distribution in the domain, Vm represents the voltage measurement on the boundary, F(s) is a nonlinear function of s that represents the induced boundary voltage with respect to s and the injected currents, α is a regularization parameter and sref is a reference conductivity distribution.

A typical attempt to solve this nonlinear optimization problem is to use a finite element method under the assumption that the conductivity in the domain is piecewise constant [4]. Within this setup we could represent Vm as a M × 1vector and s as a N × 1 vector, where N denotes the number of elements contained in the finite element mesh. The image reconstruction procedure consists of two interactive parts, namely the forward model and the inverse solver. In the forward model, a boundary voltage function F (s) can be calculated from the given parameters of the finite elements, while the inverse solver provides an estimate how to update these parameters. From here on, all discussions of this study are based on the finite element framework.

2 Methods

2.1 Gauss-Newton iterative method

To solve the nonlinear optimization problem (see Eq. 1), we could consider a stepwise linear solution known as iterative Gauss-Newton method (GN for short). According to GN, the estimated conductivity distribution ŝ is approximated by the conductivity distributions sn, and F(sn+1) is replaced by a linear approximation

F ( s n + 1 ) = F ( s n ) + J n ( s n + 1 s n )

where Jn is the Jacobian matrix evaluated at sn. The solution of Eq. 1 with respect to sn+1 is signed to be sn+1 for the next iteration. Explicitly, the estimation is updated by

(2) s n + 1 = s n + ( J n * J n + α 2 I ) 1 ( J n * ( V m F ( s n ) ) + α 2 ( s r e f s n ) )

where I is the N × N identity matrix [4]. Usually, due to the sparsity of Jn, the matrix J*nJn has a very large condition number. For the purpose of stabilizing the calculations to get a unique solution, a regularization penalty with parameter is included in Eq. 1, hence also in Eq. 2.

2.2 Dual model

Although a finer mesh can lead to a more precise calculation of F (s) in the forward model, its huge number of unknown parameters the same mesh may result in an even more severely ill-posed inverse problem. To circumvent this difficulty, a dual model framework is used [5], whose workflow is depicted in Figure 1.

Figure 1 Dual model framework: The procedure starts with solving a forward model under a fine mesh. The coarser mesh in the bottom right corner is used to solve the inverse problem. The solution is mapped back to the fine mesh for the next iteration.

Figure 1

Dual model framework: The procedure starts with solving a forward model under a fine mesh. The coarser mesh in the bottom right corner is used to solve the inverse problem. The solution is mapped back to the fine mesh for the next iteration.

2.3 An integrated algorithm

In this subsection, an approach for constructing a dual mesh is demonstrated. Given a domain and the corresponding fine mesh, a first guess of the conductivity distribution can be achieved by applying the Gauss-Newton method with this fine mesh. In this step, it is enough to quickly obtain an initial image in order to get a rough estimation. Based on the initial image, the k-means clustering method is applied to roughly detect the underlying geometry. Here the number of clusters is adjusted according to the complexity of the underlying geometry shown in the initial image as well as the number of elements in the fine mesh. The distance function employed by the k-means clustering method can be a combination of the conductivity difference in the initial image and the Euclidian distance between the finite elements. Eventually, an irregular coarser mesh is generated; each of its elements corresponds to one cluster. It is essential to modify the Jacobian matrix for the coarse mesh. According to theoretical calculations [6], the new Jacobian matrix not only has smaller size but also relatively larger entry values.

The algorithm (CGN, Clustering integrated GN) that integrates k-means clustering into the dual model framework consists of the following steps:

  1. Fix a finite element mesh that is sufficiently fine.

  2. Calculate the initial reconstruction image using this fine mesh with a few iterations of Gauss-Newton method with a stopping threshold very higher than the optimum, or an advance fixed small number of iterations.

  3. Applying k–means clustering, segment the finite elements of the fine mesh into Nc clusters and construct a coarser mesh with these Nc elements (dual model).

  4. Calculate the corresponding Jacobian matrix with size Nc × M for the coarser mesh from the Jacobian matrix for the fine mesh.

  5. Run the Gauss-Newton method iteratively with the dual model to get the final image.

  6. Optionally, further using the fine mesh, applying standard algorithms locally on the area of focus by forcing the remained area as background could provide a better image (lCGN).

The workflow of this integrated algorithm is represented in Figure 2.

It is worth to mention that the clustering based dual mesh construction could be integrated with other inverse problem solvers according to the same scheme. As a special example, instead of using iterative GN method, an alternative is to employ one step GN solver in the above framework (see Subsection 3.2).

Figure 2 The flowcharts of the clustering aided dual model framework (CGN). Optional extension with subsequent locally applied GN (lCGN).

Figure 2

The flowcharts of the clustering aided dual model framework (CGN). Optional extension with subsequent locally applied GN (lCGN).

3 Simulation and experimental results

3.1 Simulations

To evaluate the performance of this integrated algorithm, simulations were executed based on a circular finite element domain with radius 1 (see Figure 3(a)). On the boundary, 16 electrodes were placed equidistantly to inject currents and to get the voltage measurements. The stimulating currents were injected in the adjacent pattern to drive the simulation. With the predefined background conductivity (s0 = 1.0 S·m−1) a measurement of the voltage denoted by Vh could be calculated. After embedding some small contrasts into the background, another voltage measurement denoted by Vih was simulated. Moreover, a 1% white noise was added to the simulated measurements Vih. All the simulated data were calculated through an independent finite element model (mesh in Figure 3(a)) to avoid the so called “inverse crime” problem [7]. While in reconstruction a fixed fine mesh was employed for forward model calculations (mesh in Figure 3(b)). According to Section 2.3 the coarser mesh is built on this fine mesh.

The performance of the integrated algorithm was compared to the standard Gauss-Newton iterative algorithm. The difference images between reconstructed conductivity distributions and the background are shown in Figure 3. In the subsequent simulations, the standard algorithms from the EIDORS toolbox [8] were applied, where the regularization parameter in Eq. 1 was heuristically specified to be 0.01. In simulation, two iterations of the Gauss-Newton algorithms were employed to get the initial image; then fif-teen clusters were obtained through the k-means method with three inner iterations, and finally the iterations based on the dual model framework. For comparison, solely using the fine mesh, the standard Gauss-Newton method was iterated ten times to achieve its optimal solution. Viewing the conductivity as a vector corresponds to the elements of the fine mesh, the reconstruction errors of fifty independent runs were summarized by l1 and l2 norm in Table 1.

Table 1

Mean reconstruction errors of fifty independent runs.

Method GN lCGN
l1(S·m−1) 153.919 47.283
l2(S·m−1) 5.543 3.865
Figure 3 Simulation with simple contrasts. Figure 3(a) shows the ground truth with contrasts of conductivities 2.0 Sm−1 and 0.5 Sm−1. Figure 3(b) is the image obtained via ten iterations of standard GN method, the computation time is 24.832 s. Figure 3(c) demonstrates the piecewise constant conductivity map by applying three iterations of GN on fifteen clusters (CGN). Figure 3(d) is the reconstruction image by applying seven iterations of GN locally based on Figure 3(c), the total computation time of the lCGN is 28.428 s.

Figure 3

Simulation with simple contrasts. Figure 3(a) shows the ground truth with contrasts of conductivities 2.0 Sm−1 and 0.5 Sm−1. Figure 3(b) is the image obtained via ten iterations of standard GN method, the computation time is 24.832 s. Figure 3(c) demonstrates the piecewise constant conductivity map by applying three iterations of GN on fifteen clusters (CGN). Figure 3(d) is the reconstruction image by applying seven iterations of GN locally based on Figure 3(c), the total computation time of the lCGN is 28.428 s.

3.2 Experimental results

In this subsection experiments were implemented to achieve the changes of conductivities between two ventilation states, namely the start and end point of one inspiratory phase. The breathing data were acquired from a healthy volunteer by the commercial EIT device (PulmoVista500, Draeger, Luebeck, Germany). Two groups of comparisons were performed off line to test the algorithms (see Figure 4), the first group was between one step GN solver and one step CGN, the second group was between conjugate gradient solver (CG) and Clustering integrated CG (CCG).

Figure 4 Image reconstruction of conductivity change using experimental data. Images use the same colour scheme. Upper row: One step GN and one step CGN with the regularization parameter equals 0.01. Lower row: Iterative CG solver and CCG with the error threshold equals 0.001. The clustering integrated algorithms (CGN, CCG) employed ten clusters achieved from k-means method.

Figure 4

Image reconstruction of conductivity change using experimental data. Images use the same colour scheme. Upper row: One step GN and one step CGN with the regularization parameter equals 0.01. Lower row: Iterative CG solver and CCG with the error threshold equals 0.001. The clustering integrated algorithms (CGN, CCG) employed ten clusters achieved from k-means method.

4 Discussion and conclusion

The reconstructed images (Figure 3, Figure 4) demonstrate that in spite of reducing the degree of freedom in the inverse problem, the k-means method has the potential to roughly preserve the geometry around the area of interest through clustering. As indicated in the simulations and experiments, the integrated method also filtered out the undesired artefacts on the background. In addition, the simulation result (Figure 3(d)) suggests that performing standard algorithms locally on the area of interest could improve the final image.

A challenge for this integrated method is to achieve an optimal segmentation of the initial image by k-means clustering. Selection of parameters for clustering should regard to the complexity of the contrasts, the number of elements in the fine mesh and the level of noise.

Further investigation might consider adaptively determining the parameters, such as the number of clusters or the distance function used in the k-means clustering method, to get better and more controllable segmentation. In addition, individual CT information, e.g. incited by [9], could provide a helpful source for clustering.

Acknowledgment

This work was supported by the BMBF grant no. 03FH038I3 (MOSES).

Author's Statement

    Conflict of interest: Authors state no conflict of interest. Material and Methods: Informed consent: Informed consent has been obtained from all individuals included in this study. Ethical approval: The research related to human use has been complied with all the relevant national regulations, institutional policies and in accordance the tenets of the Helsinki Declaration, and has been approved by the authors’ institutional review board or equivalent committee.

References

[1] Z. Zhao, et al., “The EIT-based global inhomogeneity index is highly correlated with regional lung opening in patients with acute respiratory distress syndrome,” BMC Res Notes, 2014, vol. 7, p. 82. Search in Google Scholar

[2] Z. Zhao, et al., “Regional airway obstruction in cystic fibrosis determined by electrical impedance tomography in comparison with high resolution CT,” Physiol Meas, Nov 2013, vol. 34, pp. N107-14. Search in Google Scholar

[3] S. Leonhardt and B. Lachmann, “Electrical impedance tomography: the holy grail of ventilation and perfusion monitoring?,” Intensive Care Med, Dec 2012, vol. 38, pp. 1917-29. Search in Google Scholar

[4] D. S. Holder, “Electrical Impedance Tomography: Methods, History and Applications, Institute of Physics ” Institute of Physics Publishing, 2004, pp. 3-64. Search in Google Scholar

[5] K. D. Paulsen, et al., “A dual mesh scheme for finite element based reconstruction algorithms,” IEEE Trans Med Imaging, 1995, vol. 14, pp. 504-14. Search in Google Scholar

[6] N. Polydorides and W. R. B. Lionheart, “A Matlab toolkit for three-dimensional electrical impedance tomography: a contribution to the electrical impedance and diffuse optical reconstruction software project,” Measure. Sci. Tech., 2002, vol. 13 pp. 1871–1883. Search in Google Scholar

[7] W. R. Lionheart, ”EIT reconstruction algorithms: pitfalls, challenges and recent developments,” Physiol Meas, Feb 2004, vol. 25, pp. 125-42. Search in Google Scholar

[8] A. Adler and W. R. Lionheart, “Uses and abuses of EIDORS: an extensible software base for EIT,” Physiol Meas, May 2006, vol. 27, pp. S25-42. Search in Google Scholar

[9] Z. Zhao, et al., “Individual thorax geometry reduces position and size differences in reconstructed images of electrical impedance tomography,” J Xray Sci Technol, 2014, vol. 22, pp. 797-807. Search in Google Scholar

Published Online: 2015-9-12
Published in Print: 2015-9-1

© 2015 by Walter de Gruyter GmbH, Berlin/Boston

This article is distributed under the terms of the Creative Commons Attribution Non-Commercial License, which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.