|Home | About | Journals | Submit | Contact Us | Français|
We present an adaptive algorithm for solving the inverse problem in electrical impedance tomography. To strike a balance between the accuracy of the reconstructed images and the computational efficiency of the forward and inverse solvers, we propose to combine an adaptive mesh refinement technique with the adaptive Kaczmarz method. The iterative algorithm adaptively generates the optimal current patterns and a locally-refined mesh given the conductivity estimate and solves for the unknown conductivity distribution with the block Kaczmarz update step. Simulation and experimental results and numerical analysis demonstrate the accuracy and the efficiency of the proposed algorithm.
Electrical impedance tomography (EIT) reconstructs the internal conductivity distribution from electrical measurements made on an object's surface. It provides a workable real-time modality with low cost, non- ionization and good temporal resolution (Saulnier 2001). Therefore, EIT is of great interest for clinical applications such as monitoring lung and heart functions and detecting breast cancer (Holder 2005). Because of the ill-posed characteristics of EIT, appropriate reconstruction algorithms need to be chosen to improve the quality of the reconstruction. However, the computational cost grows quickly using conventional forward and inverse solvers. One of the challenging topics in the EIT field is to balance the trade-off between the computational efficiency of the reconstruction algorithms and the accuracy of the reconstructed images.
The computational challenges should be appropriately addressed because many realistic problems require the solution of large scale 2-dimensional (2D) or 3-dimensional (3D) problems up to a certain accuracy (Lionheart 2004, Vauhkonen 2004). One significant contribution to the computational cost is solving the forward problem using the finite element model (FEM) (Boyle et al 2012). The FEM is widely used to obtain the numerical solution of the forward problem in EIT by incorporating domain geometry, boundary conditions and prior information about the conductivity distribution. The accuracy of the numerical approximation is largely dependent on the construction of the finite element meshes. If the mesh discretization is too coarse, it may lead to inaccurate solution of the forward model and low spatial resolution of the reconstructed images. If the mesh discretization is ad-hoc and too fine, it requires a significant amount of memory storage and computational power to generate such a mesh and even more to solve the inverse problem, due to the increase in the number of unknowns that need to be solved for. Another major part of the computational cost comes from computing the sensitivity matrix and solving the inverse problem (Mueller et al 2012). Conventional iterative methods, e.g. Gauss Newton, are cost expensive especially in realistic problems where thousands of conductivity voxels need to be determined (Horesh et al 2005). In order to avoid the large-dimensional matrices during the update step, the Kaczmarz method should be considered as a strong candidate for solving the inverse problem.
The objective of the present work is to strike a balance between the accuracy and the efficiency of the reconstruction by developing an adaptive algorithm. The adaptive algorithm can locally refine the mesh and generate optimal current patterns to improve the accuracy of the reconstructed images with the help of the prior information obtained from previous iterations. The problem of constructing optimal meshes and adaptively determined meshes in EIT has a long history. Borcea et al in 2008 constructed an optimal mesh based on a resistor network model. Molinari et al in 2001 showed that the idea of refining the region where there are large conductivity gradients can effectively improve the efficiency of the forward solver and achieve good spatial resolution with a smaller number of mesh elements. The cost efficiency of the inverse solver can be improved using the adaptive Kaczmarz method, which solves the inverse problem by adaptively applying the optimal current patterns for distinguishing the actual conductivity from the conductivity estimate between each iteration of the block Kaczmarz algorithm. The algorithm proposed here is an extension of our previous work on the adaptive Kaczmarz method (Li et al 2013). It naturally combines the adaptive mesh refinement idea with the adaptive Kaczmarz method and optimizes the reconstruction process with efficient meshing and optimal current patterns. Furthermore, the optimal current patterns are ordered by the corresponding singular values of the Neumann to Dirichlet map in a descending order in the algorithm. The block structure of the adaptive Kaczmarz method provides the flexibility to reconstruct using measurements made using the dominating current patterns, i.e. the current patterns with the largest singular values. Reconstruction with these partial measurements speeds up the data acquisition, improving the temporal resolution of the images, and further improving the efficiency of the reconstruction algorithm.
This paper is organized as follows: first we discuss the idea of adaptive mesh refinement and the Kaczmarz method. Next, the adaptive algorithm which combines the adaptive meshing and the adaptive Kaczmarz method is introduced. Simulation results are shown and reconstructed images are evaluated quantitatively using several image quality metrics. A reconstruction from experimental data is also shown. We show that the proposed adaptive algorithm is able to produce stable and accurate solutions cost-efficiently.
In this paper the inverse problem of EIT is viewed as a least-squares problem. The objective functional to be minimized in the least-square sense is:
where ‖ · ‖2 denotes the L2 norm, ρ is the piecewise constant conductivity distribution, U(ρ) are the voltages that are predicted by a forward model given that the conductivity distribution in the body is ρ, and V are the measured voltages on the electrodes.
The motivation for the adaptive mesh refinement stems from the fact that uniformly refined meshing may generate a large number of elements that are not necessarily needed, resulting in a significant waste of computational power. These elements are not necessary in the sense that: (1) in particular regions additional refinement might not improve the accuracy of the final solution; (2) the refinement step introduces more unknowns needed to be determined in the inverse problem, which makes the already ill-posed problem even harder to solve. It should be noted that a much finer forward mesh is used in the forward computation and a different coarse mesh is used in the inverse computation to avoid the inverse crime. Our focus in the paper is the adaptive refinement of the inverse mesh.
Adaptive refinement techniques automatically make adjustments to the mesh to achieve the desired accuracy with respect to a specific error measure and to do so in an optimal fashion. These techniques have been extensively used in many finite element method applications (Babuška et al 1986, Aiffa 1997). Typically, in the EIT framework the algorithm should start with a reasonably uniform coarse mesh and solve for the conductivity distribution based on this coarse mesh (Molinari 2002). Conductivity gradients around the ith element are computed :
where lij is the length of the edge shared by elements i and j and the sum is over all edges of element i. The gradient is used as guidance in determining the need for local refinement for the next iteration. The coarse elements with higher conductivity gradients will be refined for a smoother reconstruction and better spatial resolution since these elements are generally along the boundaries within the region being imaged.
Triangular elements are used in both the forward and inverse mesh. The h-refinement technique is used to divide one coarse element into 4 elements by approximately connecting the midpoints of all edges of the coarse element (Babuška et al 1983). An example of such refinement is given in Figure 1. During the implementation, the inverse mesh is appropriately chosen so that the elements in the inverse mesh are essentially the superset of the elements in the forward mesh. This particular step will make the conductivity mapping from the coarse mesh to the fine mesh much easier and avoid irregular node issues in the forward solver. From the simulation results which will be shown in the next section, it is observed that 2 or 3 refinement steps are sufficient to reduce the image error and provide good resolution in the reconstructed images.
The adaptive Kaczmarz method, which was introduced and described in more detail in (Li et al 2013), combines the block Kaczmarz method and optimal current pattern generation. The block Kaczmarz method is an iterative algorithm to solve the least-squares problem (1) by recursively solving each row submatrix (a set of multiple current excitation patterns in the EIT case) of the system (Natterer et al 2001, Vauhkonen 2004). An iteration is defined to be a single execution to pass through all applied current patterns and the corresponding voltage measurements. The algorithm starts with an initial estimate and a proper partition scheme to divide the Jacobian matrix J into n subsets. At the kth iteration of the ith subset, the solution is updated as follows:
where i = 1, …, n is the subset index, ak,i is the step length at iteration k for subset i, Jk,i is the row submatrix of the Jacobian obtained at iteration k corresponding to the ith subset, Vi and Ui(ρk,i−1) are the measured and predicted voltages corresponding to the ith subset given the conductivity distribution is ρk,i−1.
The complete set of current patterns and measurements should be appropriately partitioned into subsets to optimize the block Kaczmarz method. With the help of optimal current pattern generation, an ordered subset scheme is used to partition and order the subsets of current patterns and measurements by their ability to distinguish the actual conductivity distribution from a conductivity estimate (Isaacson 1986). The optimal current patterns are the currents I that satisfy max , where V and U(ρ(0)) are the measured voltages and the predicted voltages given a conductivity estimate ρ(0) (Isaacson et al 1996).
The flowchart detailing how the adaptive Kaczmarz method works is shown in Figure 2. At each iteration, the adaptive Kaczmarz first generates the optimal current patterns Iopt to best distinguish the actual conductivity distribution from a conductivity estimate ρint. The corresponding Jacobian matrix is partitioned using the ordered subset scheme and the subsets are used in one-iteration of the Kaczmarz method to obtain the updated conductivity estimate ρupd. The more accurate conductivity estimate ρupd is used for the next iteration to find the optimal current patterns Iopt and proceed to find a new conductivity estimate. Therefore, a new conductivity estimate and a set of optimal current patterns are obtained adaptively at each iteration. The algorithm keeps iterating until an accurate conductivity estimate is found, or in the best case, the conductivity estimate and the actual conductivity are not distinguishable by the voltage difference.
The goal of the proposed algorithm is to balance the trade-off between the accuracy and the efficiency of the EIT reconstruction process. Mesh generation and the inverse solver are two computationally expensive steps in EIT reconstruction. With the help of adaptive mesh refinement, it is possible to avoid unnecessary mesh elements and improve the efficiency and the conditioning of the reconstruction. The adaptive Kaczmarz method inverts , which has a much smaller dimension than JT J used in the Gauss-Newton method, at each iteration and adaptively generates optimal current patterns to improve the distinguishability of the system. It solves the inverse problem efficiently without degrading the quality of the reconstructed images. Here, a novel adaptive algorithm is proposed that combines the adaptive Kaczmarz method with the adaptive mesh refinement.
The algorithm starts with a relatively coarse mesh, an initial conductivity estimate and a set of orthonormal current patterns. Then it proceeds to compute the updated distribution using the adaptive Kaczmarz method on the coarse mesh. Conductivity gradients are then computed and the elements with high gradients are refined. The updated mesh, conductivity estimate, and the corresponding current patterns are then used for the next iteration. The flowchart of the adaptive algorithm is shown in Figure 3.
In this section, reconstructed images using the proposed adaptive algorithm are presented and evaluated. Simulations are conducted in the EIDORS environment (Adler and Lionheart 2006). A normalized 2D circular tank with unit radius is modeled. 32 equally-spaced square electrodes are used with the assumptions of the complete electrode model (Cheng et al 1989). These square electrodes are modeled to have height and width of and the contact impedance of 0.0001 Ω.
Two error metrics are used to quantify the performance of the algorithms. The first is the normalized norm difference between the measured voltages and the predicted voltages given a conductivity estimate ρ. This metric quantifies the accuracy of the reconstructed conductivity estimate in terms of the measurement difference:
The other metric is the mean squared error between the reconstruction and the known “ground truth” (the true values of the simulated conductivity distribution), which is widely used in image reconstruction. In this case, the cumulative squared error between the reconstructed image and the original image is computed. The reconstructed conductivity distributions are projected on a square grid with size 256×256, which makes it possible to compare reconstructed images with different inverse meshes. The mean squared error is defined as:
where I is the reconstructed image after projection on the squared grid, and I0 is the ground truth image after projection on the squared grid, m = n = 256 in our case.
In this simulation, two circular targets with radius r = 0.15 are placed inside the unit radius tank. The conductor with conductivity 2 S/m is at (0.4, −0.35), and the insulator with conductivity 0 S/m is at (−0.3, 0.2). The background conductivity is 1 S/m. Trigonometric patterns with the maximum amplitude of 1 mA are used as the initial current patterns for the adaptive algorithm. Reconstructed images of the first 4 iterations are shown in Figure 4. The reconstruction meshes are shown in Figure 5.
Since we choose to refine the elements with high gradients, most of the refined elements are along the edges of two targets or within the targets. The spatial resolution of the reconstructed images is greatly improved with the help of the adaptive mesh refinement. The same conclusion can be drawn from the error metrics shown in Table 1: the normalized norm voltage difference is reduced monotonically from 0.0278 at iteration 1 to 0.0005 at iteration 4, and the mean squared error is reduced from 0.0390 at iteration 1 to 0.0031 at iteration 4.
We also compare the results with the reconstruction using Gauss-Newton Method with a fixed uniform mesh in Table 2. The uniformly refined mesh with 1830 elements is obtained by gradually increasing the number of the elements until it achieve the similar error performance (E1=0.0005 and E2=0.0037) compared to the adaptive refined mesh algorithm. These results illustrate the efficiency of the adaptive algorithm, which increases the mesh density locally with prior information of the solution. Moreover, the size of the matrix to be inverted is 1830 × 1830 in the Gauss Newton case, compared to 192 × 192 which is the size of the matrix to be inverted in the adaptive algorithm at each iteration given the block size is 6. In terms of the memory requirements, the Gauss Newton with 4 iterations requires 45.6 Megabytes, while the adaptive algorithm requires 15.5 Megabytes. The operation count is O(kn3) in Gauss-Newton and O(km2n) in the adaptive Kaczmarz, where k is the number of iterations, m is the number of measurements, and n is the number of elements in the forward computation. Practically, n is much bigger than m and the Kaczmarz update step is more efficient. Reconstructions of the adaptive Kaczmarz using partial measurements will be shown in the later section and it should be noted the computation cost is further reduced by using a subset of the measurements.
We further test the adaptive algorithm on a more complicated lung and heart simulated phantom. A lung and heart phantom is represented by two ellipses and a circle. The background conductivity is 0.424 S/m, the conductivity of the lungs is 0.240 S/m, and the heart conductivity is 0.750 S/m. Trigonometric patterns with the maximum amplitude of 1 mA are used as the initial current patterns for the adaptive algorithm. Reconstructed images from the first 4 iterations are shown in Figure 6. Even with the more complicated targets, the adaptive algorithm can produce accurate reconstructions in 4 iterations. Similar to the two-target simulation results, the improvement of the spatial resolution from iteration to iteration is easily observed. Furthermore, both the norm voltage difference and the mean squared error have been reduced substantially over 4 iterations, which is shown in Table 3.
The block structure of the Kaczmarz method allows us to reconstruct using only a subset of the measurements instead of the full set of measurements. The optimal current pattern generation step in the algorithm orders the current patterns and the corresponding measurements by their ability to best characterize the unknown conductivity distribution, or in other words, by sorting the corresponding singular values of the Neumann to Dirichlet map in descending order. In fact, by examining the singular values associated with the current patterns in the lung and heart simulation, it is found that the singular values associated with the first 8 patterns contain 99.5% of the energy contained in all 31 singular values, the singular values associated with the first 16 patterns contain 99.8% of the total energy in all 31 singular values.
Reconstructing using partial measurements can speed up the data acquisition and improve the temporal resolution of the images. It will also reduce the number of subsets needed in the block Kaczmarz method and accelerate the reconstruction. It is natural to think about whether we can use only partial measurements without degrading the reconstruction quality. In Figure 7, only the first half of the measurements, i.e. those associated with the first 16 current patterns, are used in the reconstruction. Compared to the reconstruction made using the full set of measurements, the image quality is only slightly affected in terms of the mean squared error performance, which is shown in Table 4. The shape and the location of lungs and heart are accurately recovered with few artifacts after 4 iterations. The norm voltage difference is bigger compared to the case with full measurements due to the reduced number of current patterns.
The technique can be extended further by varying the number of measurements for each iteration. The reconstruction quality of the first 2 iterations is governed mostly by the coarse mesh grid, since there are a limited number of elements to represent the conductivity distribution. While the reconstruction quality of the last 2 iterations is governed mostly by the quality of the measurements and the reconstruction algorithm after a sufficient number of elements are generated through the refinement procedure. Therefore, we can start with a small number of measurements for the coarse mesh and gradually increase the measurements as the mesh is refined during iterations.
In the simulation shown in Figure 8, only 6.25% of the measurements are used for reconstruction on the coarse mesh at iteration 1 and 12.5% of the measurements are used for reconstruction at iteration 2. The ratio gradually increases to 25% for iteration 3 and 50% for iteration 4 as the inverse mesh becomes finer. From the simulation results, it can be observed that for the first 2 iterations, even around 10% of the measurements can provide a valid initial estimate for later iterations since the image quality is mostly limited by the coarseness of the mesh. However, as the mesh gets finer, the image quality is governed by the measurement quality and the reconstruction algorithm, therefore, more measurements and current patterns provide additional useful information for the last 2 iterations. The error metrics in Table 5 show that the reduction in the number of measurements does not affect the image quality, since the norm voltage difference and the mean squared error are only slightly increased compared to the case with full measurements.
In this section, we use the adaptive algorithm to reconstruct experimental data collected on a phantom using the ACT3 system (Edic et al 1995) at Rensselaer Polytechnic Institute using the full set of measurements. The phantom consisted of agar heart and lungs in a saline bath in a 2D circular tank with radius of 15 cm. 32 square stainless steel electrodes of size 2.5 cm high and 2.5 cm wide are used. The conductivity of the saline is 0.424 S/m, the conductivity of the agar lungs is 0.240 S/m, and the conductivity of the agar heart is 0.750 S/m. Trigonometric current patterns are applied with a maximum current amplitude of 0.2 mA. The synthesized voltage measurements are computed based on the optimal current patters generated at each iteration. Reconstructed images for the first 3 iterations are shown in Figure 9 and the improvement in the spatial resolution over iterations can be easily observed. It should be noted that the full benefit of the adaptive Kaczmarz is not available here since we did not apply the optimal current patterns to the tank and take actual voltage measurements at each iteration. Instead the synthesized voltage measurements are computed and used at each iteration. Therefore a full implementation in which the optimal currents are applied would be expected to provide even better results.
We have presented an adaptive algorithm which incorporates an adaptive mesh refinement procedure in the adaptive Kaczmarz method. The adaptive algorithm is able to refine the mesh locally and generate the optimal current patterns at each iteration in an adaptive way. The adaptive mesh refinement step avoids the generation of unnecessary elements and improves the spatial resolution efficiently. The adaptive Kaczmarz method solves the inverse problem cost-efficiently and accurately by generating optimal current patterns at each iteration. Simulation results are shown to demonstrate the reconstruction performance of the adaptive algorithm using different targets in a 2D tank. Moreover, the optimal current pattern generation and the block structure of the Kaczmarz methods allow us to use only a subset of measurements to produce accurate reconstructions. This approach further improves the efficiency of the algorithm and the temporal resolution of the reconstruction, making it favorable in real-time applications.
Future work includes studies of alternative mesh refinement techniques, e.g. Delaunay criterion, to obtain better quality meshes. Local triangulation used in the Delaunay technique can ensure the regularity of the elements and potentially increases the accuracy of the forward solver. We also plan to study a good way to pick regularization parameters over iterations. A general principle to follow is that smaller regularization is needed for coarser meshes and larger regularization is needed for finer meshes, depending on the number of unknowns needed to be solved for different meshes. However, a systematic way should be developed to quantify the appropriate regularization over iterations.
This publication was made possible by Grant Number 1R01HL 109854-01 from the National Institutes of Health. Its contents are solely the responsibility of the authors and do not necessarily represent the official views of the National Institutes of Health.