The D-bar Method for Electrical Impedance Tomography – Simulated Data

Summary: Electrical Impedance Tomography (EIT) aims to recover the internal electrical conductivity of a physical body from measurements of voltages and currents at the boundary of the body. EIT has applications in medical imaging, underground prospecting, and nondestructive testing. The image reconstruction problem of EIT is a nonlinear and severely ill-posed inverse problem. The D-bar method is a non-iterative regularized reconstruction method based on low-pass filtering a nonlinear Fourier transform. This page contains Matlab routines implementing the D-bar method for simulated EIT data.

Authors: Jennifer Mueller, Samuli Siltanen and Janne Tamminen.

Software: Matlab, including PDE toolbox (probably works also with Octave).

Download package as zip file: DbarEIT_Matlab

Literature: Chapters 12-15 of the book Linear and nonlinear inverse problems with practical applications by Jennifer L Mueller and Samuli Siltanen (SIAM 2012). Below we refer to this book by [MS2012].

How to compute regularized reconstructions from EIT data, given that the inverse problem is both nonlinear and extremely sensitive to modelling errors and measurement noise? One can of course use variational regularization: write a functional consisting of a data discrepancy term and a regularization term, and minimize using an iterative algorithms such as the Gauss-Newton method. While this is a reasonable approach, there are two drawbacks. The iteration may get stuck to a local minimum without you noticing, and it is computationally expensive to solve a direct problem in every iteration.

The D-bar method bypasses these problems: it is a direct method based on computing a nonlinear Fourier transform of the conductivity from EIT data and then inverting the transform. There is no need to solve any direct problem. Regularization is provided by low-pass filtering in the nonlinear frequency domain. The reconstruction is blurry because of the low-pass filter, but this can be helped in post-processing by anisotropic diffusion, frequency-domain filling, or deep learning.

The motivation for this blog post comes from the following medical imaging scenario, where a row of electrodes is attached around the chest of a patient:

We aim to image the heart and lungs of the patient using EIT. We think of maintaining a given voltage potential at each electrode and measuring the resulting electric current through the electrodes. This measurement is repeated for several voltage patterns. Note that practical EIT devices typically feed currents and measure voltages; however, in this blog post we consider voltage-to-current measurements for mathematical convenience. This is not a huge deal since one can always switch between the two data types computationally.

Assuming that there are 32 electrodes in total, we can use at most 31 linearly independent voltage patterns since one of the electrodes is considered as the ground potential to which the other potentials are compared to. We use trigonometric voltage patterns, approximated by continuous sine curves at the boundary:

Furthermore, we approximate the above three-dimensional (3D) situation with a two-dimensional computational
model. Our virtual patient will be modelled by a two-dimensional (2D) disc with varying electrical conductivity
inside. This 2D approximation, while obviously incorrect, works surprisingly well in practice.

Next we show how to simulate the above kind of voltage-to-current data and how to recover the inner
conductivity from the boundary measurements using the D-bar method.

Definition of the conductivity

Our simulated conductivity is defined in the file heartNlungs.m, and you can plot it using the routine DbarEIT01_heartNlungs_plot.m. The result should look something like this:

Here the background conductivity is 1, the heart is filled with blood and has higher conductivity 2, and the lungs are filled with air and have lower conductivity 0.5.

Wonder about the white color at the bottom of the colorbar above? It is for is for creating the white background in Matlab. Check out this page to learn more.

Simulation of EIT data

We simulate by first computing the current-to-voltage map using the Finite Element Method (FEM) and then inverting. In what follows, the voltage-to-current map is called DN map for Dirichlet-to-Neumann, and the current-to-voltage map is called ND map for Neumann-to-Dirichlet.

FEM is an efficient method for solving elliptic partial differential equations. The solution of the conductivity equation is given as a linear combination of basis functions that are piecewise linear in a triangular mesh. The mesh is constructed by the routine DbarEIT02_mesh_comp.m containing the parameter Nrefine. Here are plots of the triangle meshes with Nrefine=0 and Nrefine=2:

In practice it is a good idea to use Nrefine=5 for accurate results. The mesh is saved to a file called data/mesh.mat. Note that the routine DbarEIT02_mesh_comp.m creates a subdirectory called data. If the subdirectory data already exists, Matlab will show a warning which you can safely ignore.

The simulation of the continuum-model ND map involves solving Neumann problems of the form

containing a Fourier basis function defined by

Both the Neumann condition and the conductivity given to Matlab’s PDE toolbox routines as user-defined functions. Let us see how to do that. First of all, the user-defined routines need to define the functions at appropriate points related to the FEM mesh. Here are some important points:

The conductivity is constant inside each triangle, so it is naturally specified at the centers of triangles shown in (a). This is implemented in the routine FEMconductivity.m. The solution produced by FEM is linear inside each triangle, so the gradient is constant inside each triangle. Therefore it is appropriate to specify the Neumann data at the midpoints (c) of boundary segments. This is done in routine BoundaryData.m. Finally, we need to find the trace of the solution for constructing the ND map; this can be done simply by picking the values of the FEM solution at the vertices (b) located at the boundary.

The routine DbarEIT03_ND_comp.m uses FEM to compute a matrix approximation to the ND map and saves it to the file data/ND.mat. See Chapter 13.2 of [MS2012] for mathematical details of this computation. Note: you may see some warnings in Matlab when runningDbarEIT03_ND_comp.m, but you can safely ignore them.

Computation of the scattering transform via the boundary integral equation

We need to construct a grid in the complex k-plane for the evaluation of the complex-valued scattering transform t(k), also known as nonlinear Fourier transform. Run the routine DbarEIT05_Kvec_comp.m. You should see something like this:

The next step is to compute the traces of Complex Geometric Optics (CGO) solutions by solving a boundary integral equation. This is done in routine DbarEIT05_psi_BIE_comp.m. This routine calls the function solveBIE.m. The traces are then fed to a boundary integral formula evaluating the scattering transform. The evaluation is done in routine DbarEIT06_tBIE_comp.m. This is how the scattering transform (also known as nonlinear Fourier transform) looks like:

Reconstruction from scattering transform using the D-bar method

Now we are ready to reconstruct the conductivity. The solution of the D-bar equation is based on generalized Vainikko’s method, which is organized into the files GV_grids.m, GV_project.m, GV_prolong.m and DB_oper.m.

Run DbarEIT08_tBIErecon_comp.m, and you should see this image:

The left image shows the original conductivity distribution, and the right image shows the reconstruction using R=4. The colors (and thus numerical values) are directly comparable. Note that the reconstruction underestimates the value at the “heart.”

You can try a bigger frequency cutoff radius as well. Open the file DbarEIT07_tBIErecon_comp.m and change R=4 to R=6, say. The result looks like this:

This time the reconstruction slightly overestimates the conductivity of the “heart;” this is a typical feature of the truncated nonlinear Fourier transform.

When you use larger cutoff frequencies R, you may need more accuracy in the computational grid. To achieve that, change M = 8 into M=9 or even higher. You can monitor the change in the reconstruction when M grows; when the change from M to M+1 becomes negligible you can stop increasing M. Note that higher M increases memory demands and computation time.