A Real-time D-bar Algorithm for 2-D Electrical Impedance Tomography Data

The aim of this paper is to show the feasibility of the D-bar method for real-time 2-D EIT reconstructions. A fast implementation of the D-bar method for reconstructing conductivity changes on a 2-D chest-shaped domain is described. Cross-sectional difference images from the chest of a healthy human subject are presented, demonstrating what can be achieved in real time. The images constitute the first D-bar images from EIT data on a human subject collected on a pairwise current injection system.


Introduction
The generalized Laplace equation serves as a model for the electric field u(x) propagating in living tissue from a low frequency, low amplitude applied current density on the surface of a body. Denoting the conductivity by σ(x) and a bounded 2-D domain by Ω, the governing equation is ∇ · (σ(x)∇u(x)) = 0, x ∈ Ω. (1) The ideal data is the Dirichlet-to-Neumann (DN) map, or voltage-to-current density map Λ σ , defined by where u is a solution to (1) and ν is the unit normal vector to the surface. Equation (1) serves as the governing equation for the electric field in electrical impedance tomography (EIT), and has a rich mathematical history dating back to the problems posed by Calderón [9]: (1) when does the inverse problem of determining σ from knowledge Λ σ have a unique solution, and (2) how can it be determined? Historical reviews of the answers to these questions can be found in [5,35], and the reader will find that most of the uniqueness proofs have utilized complex geometrical optics (CGO) solutions. Some have also been formulated as constructive proofs [36,7,2], and most of these include PDEs known as D-bar, or∂, equations for the CGO solutions. These works have given rise to a new class of direct (noniterative) reconstruction algorithms known as D-bar methods. D-bar equations are of the form∂w = f, where f may depend on w, and the∂ operator is defined by∂ z =

Background
We begin with an overview of the D-bar method implemented here, both for the reader's convenience and to place the fast implementation in its mathematical context. For further details, see [36,35]. The method begins with a transformation of the generalized Laplace equation with conductivity σ ∈ W 2,p (Ω), for some p > 1, to the Schrödinger equation through the change of variables q(z) = ∆ √ σ(z)/ √ σ(z) andũ(z) = √ σ(z)u(z), where the point z = (x, y) lies in Ω, a bounded, simply connected Lipschitz domain in R 2 . This results in − ∆ũ + q(z)ũ = 0, z ∈ Ω. ( Under the assumption that σ is constant in a neighborhood of the boundary of Ω, one can extend (3) to the whole plane, taking q = 0 outside Ω. Without loss of generality, we will assume σ ≡ 1 in a neighborhood of the boundary. The existence of CGO solutions to (3) in the plane was established by Faddeev [13] in the context of quantum physics and shown by Nachman [36] to always exist for q of the form q(z) = ∆ √ σ(z)/ √ σ(z). Introducing the complex parameter k = k 1 + ik 2 , and identifying the spatial variable z = x + iy with the corresponding point in the complex plane, the CGO solution ψ(z, k) satisfies [36] − ∆ψ(z, k) The CGO solution closely related to ψ is µ(z, k), defined by µ(z, k) ≡ e −ikz ψ(z, k). The conductivity can be obtained directly from µ or ψ through the formula [36] σ(z) = µ 2 (z, 0), z ∈ Ω, and so the key steps in the method are to link the data Λ σ to the function µ(z, k). These links are provided by the scattering transform t(k) and the D-bar equation for µ. The scattering transform is defined by and can be regarded as a nonlinear Fourier transform of q in light of the asymptotic behavior of ψ. The D-bar equation satisfied by µ is where e z (k) ≡ e i (kz+kz) . The scattering transform is related to the DN data through an equation that requires knowledge of ψ on the boundary of Ω: Here Λ 1 denotes the DN map corresponding to the homogeneous conductivity distribution σ ≡ 1. For the fast implementation, we utilize a linearized approximation to the scattering transform, denoted by t exp , which is defined by replacing ψ| ∂Ω by its asymptotic behavior: This approximation was first introduced in [37] and was later studied in [24] where it was shown that the D-bar equation (8) with t(k) replaced by t exp truncated to a disk of radius R in the k-plane has a unique solution which is smooth with respect to z, and the reconstruction is smooth and stable. Further, it was shown that no systematic artifacts are introduced when the method with t exp is applied to piecewise continuous conductivities. In summary, the mathematical steps of the algorithm implemented here are (1) compute t exp (k) on a disk of radius R in the k-plane, (2) solve the D-bar equation (8) for each z in the region of interest on the disk |k| ≤ R, and (3) compute σ from (6).

Fast implementation
A fast implementation in Matlab on a 12 core Mac Pro with two 2.66 GHz 6 core Intel Xeon processors and Matlab's Parallel Computing Toolbox is capable of computing reconstructions at less than the data acquisition rate of 16 frames/s, or 0.0625 s/frame, of the ACE 1 pairwise current injection EIT system at CSU [32]. This demonstrates the feasibility of CGO methods for real-time reconstructions.
In fact, we consider two options for the parallel computations. Ideally, in real-time reconstructions, data is collected, demodulated, and fed directly to the reconstruction algorithm, one frame at a time. In this configuration, only the loop over the z-values in the solution of the D-bar equation is trivially parallelizable. This accounts for over 95% of the computation time and can be used to obtain real-time reconstructions at a rate of 0.0621 s/frame on 7 cores on a coarse mesh of 562 elements, as reported in Section 4. This approach is structured as shown in Algorithm 1.
If a time delay of approximately one second is acceptable to the user, the algorithm can be parallelized over the frames. In this configuration, the data is collected, demodulated, and sent to a buffer from which multiple frames are sent to the reconstruction algorithm in batch. This approach results in an algorithm that is over 99% parallelizable, and the frame rate is even faster. As shown in Section 4, reconstructions were computed at a rate of 0.0215 s/frame on 64 cores or 0.0578 s/frame on 12 cores on a mesh of 1931 elements. The computational method for this approach is shown in Algorithm 2.
We now discuss the computational steps in detail, including numerical approximations of various functions and operations, and numerical solution of the D-bar equation. The computation of the matrix approximation to the DN map Λ σ can be accomplished efficiently with inner products, as explained in [35]. For bipolar current patterns, as applied by ACE 1, the current pattern matrix is first transformed to an orthonormal basis {J m l }, m = 1, . . . , N and l = 1, . . . , L, where N is the number of linearly independent current patterns, and L is the number of electrodes. The number of linearly independent current patterns depends on the particular choice of current pattern. For a bipolar current pattern that skips α electrodes between injection electrodes, N = L − α − 1. The voltages are transformed accordingly through the change-of-basis formula to the set denoted by {V m l }. The discrete ND map R σ d for a given data set σ d is approximated by 1: Setup Phase. Define parameters and compute functions independent of the dynamic DN data, including: -Physical parameters -Boundary parameterization and arclength function -Current pattern matrix J -Demodulate the reference data set and form the DN matrix for the reference data set L σ re f -Computational grids in k-plane and z-plane Load a single frame of the measured data and demodulate 3: Compute matrix approximation to DN map, L σ d 4: Compute t exp simultaneously for all z using vector operations 5: Compute the pointwise multiplication operator T R simultaneously for all z using vector operations 6: parfor all z in domain do 7: Solve D-bar equation for µ exp R (z, k) 8: where A is the area of an electrode, s l is the arclength function s(θ) for the boundary evaluated at the angle θ l corresponding to the center of the l th electrode, and ∆θ l = θ l − θ l−1 (where θ −1 = θ L ). Since the voltages sum to zero, we can then compute the matrix approximation The scattering transform is computed for |k| ≤ R, where R is chosen empirically. While better reconstructions can often be obtained by considering a non-uniform truncation radius R, we consider a disk for simplicity. Non-uniform truncation would not result in appreciable loss of computational speed since over 95% of the computation time is spent in the solution of the D-bar equation. Since difference images from a reference frame are being reconstructed here, we implement the scattering transform t exp dif , introduced in [23], in which the DN map for the conductivity σ = 1 is replaced by that of a reference frame Λ σ re f . Then The functions e ikz and e ikz are expanded in the orthonormalized current pattern basis. The coefficients of the expansions of the functions e ikz | ∂Ω and e ikz | ∂Ω in the orthonormalized current pattern basis vectors are computed in the setup phase of the algorithm, and are denoted by Then, denoting the discrete inner product over two vectors u and v of length L by (u, v) L , The fast evaluation of this formula is accomplished using inner products and vector operations. A single grid and multigrid (2-grid) method were introduced in [29] for the fast computation of Lippmann-Schwinger type equations that arise in D-bar methods for EIT, closely based on the method of Vainikko [40]. The convolutions are computed as FFTs, and the solution of the resulting linear system by a matrix-free method, such as GMRES. For difference images, and in particular the data sets considered here, the frame-to-frame change in the data is sufficiently small that the method nearly always converges in one inner and one outer iteration of GMRES. In such cases, the 1-grid method is significantly faster than the 2-grid method described in [29]. Compute matrix approximation to DN map, L σ d

5:
Compute t exp simultaneously for all z using vector operations 6: Compute the pointwise multiplication operator T R simultaneously for all z using vector operations 7: for all z in domain do 8: Solve D-bar equation for µ exp R (z, k) 9: σ(z) ← (µ exp R ) 2 (z, 0) 10:

end for 11: end parfor
To apply this method, the D-bar equation (8) is formulated as an integral equation To construct the computational k-grid, we define the square [−s, s] 2 where s ≥ R, and we choose M = 2 n for some positive integer n. The step-size for the k-grid is defined to be h = 2s/(M − 1), and the final size of the grid is then M × M. We further choose a computational z-grid of domain points, which need not be equally spaced. Equation (13) can be written compactly as the linear system where T R is the pointwise multiplication operator defined by and the action of the operator A R is given by In our fast implementation, we compute T R simultaneously for all z-values in the computational grid using vector operations, which in Matlab is more efficient than computing each T R separately inside the z-loop shown in step 6 of Algorithm 1 and step 7 of Algorithm 2. The action of the operator A R can be approximated by where β(k) = (πk) −1 is the Green's function for the∂ operator, and · denotes element-wise multiplication. In Matlab, this operation can be performed efficiently using IFFTN and FFTN. Once the solution µ exp R (z, k) to the D-bar equation has been found, the conductivity is given by σ(z) = (µ exp R ) 2 (z, 0). Note that the D-bar equation requires the solution of the linear system (14) for all k values in the computational grid, even though we are ultimately only interested in k = 0. It is also to be noted that this method allows for the reconstruction of σ pointwise in Ω, independent from any other z value, and so it is trivially parallelizable over the values in the computational z-grid.
For maximum computational speed, we invoked Matlab with the flag -singleCompThread, which limits Matlab to a single computational thread. This choice is compatible with the Parallel Computing Toolbox, and will restrict each parallel computation to a single core, which proved to be advantageous in the runtimes. Several choices of numerical solvers for the linear system (14) were compared for the fast implementation. Here, we use a customized GMRES, as opposed to the built-in Matlab routine. This avoids significant overhead in calling external .m files and performing unnecessary error-checking.

Results and discussion
To demonstrate the feasibility of clinically useful real-time reconstructions using the D-bar method, we present reconstructions from data collected on the ACE 1 (Active Complex Electrode) EIT system at CSU. The ACE 1 system is a bipolar current injection system with 32 active electrodes operating at a user-specified frequency up to 125 kHz. Details of the hardware can be found in [32]. For the results presented here, difference images of perfusion in a crosssection of the chest of a healthy male subject sitting upright and holding his breath are presented. 360 frames of data were collected at 16 frames/s at 125 kHz and current amplitude 0.823 mA. One frame was chosen as a reference data set and 359 difference images were computed using the fast implementation in Section 3 on a uniform k-grid of size 16 by 16 (256 elements) and a truncation radius of R = 3.8. All programming was in Matlab and utilized the Parallel Computing Toolbox for the parallel solution of the D-bar equation.
In Tables 1a and 2a we compare the performance for each version of the algorithm with various numbers of cores in parallel on a 12 core Mac Pro with two 2.66 GHz 6 core Intel Xeon processors using three different spatial grids consisting of a coarse mesh with 562 elements, a medium mesh with 1931 elements, and a fine mesh with 5,916 elements. In Tables 1b and 2b we compare the performance on a 64 core Linux system with four 2.3 GHz 16 core processors and 512 GB of RAM on the same three spatial grids.
Comparing all four tables, it is immediately clear that all runtimes utilizing the same number of cores were faster on the Mac Pro than on the Linux system, which is likely due in part to the difference in processor speed. It is also evident that while adding increasing numbers of parallel cores continues to improve runtimes for Algorithm 2 in Table 2, the runtimes for Algorithm 1 in Table 1 reach maximum efficiency with a smaller number of cores, after which adding additional parallel cores actually slows the computation. This optimum number of cores increases as the z-mesh size increases.
It is well-known in parallel computing that the efficiency gained by adding additional parallel processors follows a "law of diminishing returns," embodied by Amdahl's law [1], which states that the maximum theoretical speedup s obtainable when using n processors in parallel is given by where p is the proportion of the program that is parallelizable. One can see that as we increase n to ∞, the maximum theoretical speedup goes to 1/ (1 − p). To compute the actual speedup obtained using n processors, we use s actual (n) = T (1)/T (n), where T (n) is the runtime with n cores in parallel. In Figure 1, the actual speedups obtained on the Mac Pro using both versions of the parallel algorithm are shown along with the theoretical maximum speedups predicted by (18). One can see that the results obtained using Algorithm 2 are much closer to Amdahl's ideal values than the results obtained using Algorithm 1. We can also see that the larger the size of the z-mesh, the more efficient the parallelization becomes; this is predicted by (18) since increasing the number of elements in the z-mesh also increases p, the parallelizable portion of the program. In Figure 2 the same plots are included for the 64 core Linux system. It is evident that the additional cores do not provide speedup for Algorithm 1, and the divergence from the theoretically predicted speedup by Amdahl's Law increases with the number of cores, but Amdahl's Law also predicts the leveling off of speedup at around 20 to 25 cores; in fact, this happens at around 10 to 12 cores for parallelization over mesh values (Algorithm 1). However, for Algorithm 2, the speedup continues for all 64 cores, although it is clearly starting to level out at around 60 cores. Figure 3 contains four frames in the reconstruction of the human chest data displayed in the three z-meshes to illustrate the resolution provided by these mesh choices. The figure depicts changes due to perfusion. The heart is at       the top, and red represents high conductivity and blue low conductivity. The images are displayed on the same scale. While the topmost mesh is quite coarse, the lungs and heart are clearly visible, and changes are evident. The medium mesh is significantly better, and probably provides the best compromise for real-time imaging, but this implementation comes with the price of an approximately 1 second delay. The very fine mesh is included to illustrate highly desirable spatial resolution and the associated runtimes, given in Tables 1 and 2. showing changes due to perfusion in the chest of a healthy human subject. The heart is at the top, and red represents high conductivity and blue low conductivity with respect to the reference frame. The images are displayed on the same scale.
From the clinical perspective, difference images are sufficient for some applications, such as real-time detection of a pneumothorax or atelectasis. For other applications, such as distinguishing between blood, water, and mucus in the lung, absolute images may be required, which may require longer computation times. Some applications may also require finer spatial resolution than that presented here. However, the data is preserved for off-line or delayed reconstruction, and improvements such as the use of additional cores in parallel, faster FFTs, and faster processors are likely to yield improved runtimes in the future.

Conclusions
The results presented here show for the first time the D-bar method applied to human chest data collected on a pairwise current injection system. Conductivity changes due to perfusion are clearly visible in the images. The fast implementation demonstrates the clinical potential of the D-bar algorithm as a reconstruction algorithm for real-time bedside imaging.

Acknowledgments
The project described was supported by Award Number 1R21EB016869-01A1 from the National Institute Of Biomedical Imaging And Bioengineering. The content is solely the responsibility of the authors and does not necessarily represent the official view of the National Institute Of Biomedical Imaging And Bioengineering or the National Institutes of Health.