An Edge-Preserving D-bar Method for EIT

# A Data-Driven Edge-Preserving D-bar Method for Electrical Impedance Tomography

###### Abstract.

In Electrical Impedance Tomography (EIT), the internal conductivity of a body is recovered via current and voltage measurements taken at its surface. The reconstruction task is a highly ill-posed nonlinear inverse problem, which is very sensitive to noise, and requires the use of regularized solution methods, of which D-bar is the only proven method. The resulting EIT images have low spatial resolution due to smoothing caused by low-pass filtered regularization. In many applications, such as medical imaging, it is known a priori that the target contains sharp features such as organ boundaries, as well as approximate ranges for realistic conductivity values. In this paper, we use this information in a new edge-preserving EIT algorithm, based on the original D-bar method coupled with a deblurring flow stopped at a minimal data discrepancy. The method makes heavy use of a novel data fidelity term based on the so-called CGO sinogram. This nonlinear data step provides superior robustness over traditional EIT data formats such as current-to-voltage matrices or Dirichlet-to-Neumann operators, for commonly used current patterns.

###### Key words and phrases:
Inverse conductivity problem, electrical impedance tomography, image segmentation, CGO solutions, contrast enhancement
###### 1991 Mathematics Subject Classification:
Primary: 65N21, 94A08; Secondary: 42B37.

Sarah Jane Hamilton

Department of Mathematics, Statistics, and Computer Science

Marquette University

Milwaukee, WI 53233, USA

Andreas Hauptmann and Samuli Siltanen

Department of Mathematics and Statistics

University of Helsinki

Helsinki, 00014, Finland

## 1. Introduction

Noise-robust solutions of ill-posed inverse problems are based on regularization strategies. For Electrical Impedance Tomography (EIT), the only proven regularization strategy is the low-pass filtered D-bar Method, which sets high scattering frequencies to zero therefore resulting in smoothed images where sharp features important to applications such as medical imaging are often absent. In this paper, we propose to recover edges in the smoothed EIT reconstructions by applying a deblurring flow stopped at a minimal data discrepancy (Figure 1), guided by a novel data fidelity term based on the so-called CGO sinogram, which provides superior robustness.

Electrical Impedance Tomography (EIT) is an imaging modality where an unknown physical body is probed with electricity using electrodes attached to the surface of the body. The goal is to recover the internal conductivity distribution of the body typically based on current-to-voltage boundary measurements. EIT has applications in medical imaging, underground contaminant detection and industrial process monitoring. See [10] and [30, Chapter 12] for more details and applications of EIT.

We formulate the inverse conductivity problem [8] for two-dimensional EIT in terms of voltage-to-current measurements. Let be the unit disc. We model the conductivity by a bounded measurable function satisfying for almost every . For a prescribed boundary voltage , the voltage potential satisfies the conductivity equation

 (1.1) ∇⋅σ∇u=0, in Ω,u|∂Ω=ϕ, on ∂Ω.

Infinite-precision voltage and current measurements are modeled by the Dirichlet-to-Neumann (DN), or voltage-to-current density, map

where denotes the outward facing unit normal to . The goal of EIT is to recover the conductivity distribution for , approximately, from the knowledge of a practical data matrix satisfying with known noise amplitude and an appropriate data space , see [27] on details for such a space.

EIT is a severely ill-posed inverse problem, in fact it is only stable. By this we mean that small changes in boundary measurements can correspond to large changes in the internal conductivity distribution, and furthermore that noise in the data is amplified exponentially. Therefore, regularization is needed for the noise-robust recovery of from . The forward map is too nonlinear to be covered by the presently available theory of iterative (Tikhonov-type) regularization. So far the only methodology providing proven regularization properties is the so-called D-bar method in dimension two [38, 29, 27]. Regularization for the 3D case is in progress, based on [11, 6, 14].111Personal communication with Kim Knudsen and his team.

There exists a nonlinear Fourier transform that is intimately connected to EIT. Namely, Nachman showed in [33] that one can use infinite-precision EIT data to completely determine , and then apply the inverse transform, via solving a D-bar equation in the scattering variable, to recover the conductivity. However, in practice one never has such infinite-precision noise-free data; instead one works with noisy data . The basic structure of the regularized D-bar method is shown in Figure 2. Practical data only allow stable computation of the nonlinear Fourier transform in a disc centered at the origin in the frequency domain. One can then use the good part of the transform in the inversion, corresponding to a nonlinear low-pass filtering. The cut-off frequency of the nonlinear low-pass filter depends logarithmically on the noise amplitude, tending to infinity asymptotically in the zero-noise limit. Analogously to linear low-pass filtering, the resulting image is smoothed and appears blurred.

In general, noise-robust solutions of ill-posed inverse problems are based on complementing indirect and unstable measurement data with a priori information. The transform domain of the D-bar method organizes the measurement information neatly into a stable part () and unknown/unstable part (). Roughly speaking, information about the high frequencies of the unknown conductivity are missing from EIT data.

What kind of a priori information do we have? In many applications of EIT it is reasonable to assume that the conductivity is piecewise constant with crisp boundaries between the homogeneous regions. Those boundaries have significant high-frequency content which is not stably represented in the data. Modern imaging science provides several options for sharpening blurred images, which may help to recover the edges in the true conductivity. The most straightforward edge-sharpening approach, the Perona-Malik anisotropic diffusion approach [34] (see [37, 40] for its generalizations) proved insufficient in our work due to either too smooth edges or instability at high gradients. We therefore used the technique proposed by Ambrosio and Tortorelli [2] as an approximation to the classical Mumford-Shah image segmentation functional [31]. With this method, areas separate nicely and develop constant values, while the gradient can be controlled by the auxiliary variable of the model. This functional has been widely used in image processing, see for instance [9, 15, 25]. On the other hand the approach is rather uncommon for inverse problems, there are a few examples in signal restoration [22], X-ray tomography [35], and particular for EIT imaging [36]. A recent paper utilizes a Perona-Malik prior for Bayesian inversion in EIT imaging [21].

It is well-known that minimizing the Ambrosio-Tortorelli (AT) functional can be used to sharpen and segment blurry images [2, 37]. The image flows in time according to a nonlinearly deforming diffusivity. High gradients in the initial image will develop sharp edges, while slowly varying regions will tend to constant-value areas in the final image. See Section 3 for more information about AT. Simply minimizing the AT functional for the blurry D-bar reconstruction will introduce edges, but we need to ensure that the change is for the better. In particular in medical imaging applications we need to avoid introducing artifacts. We propose controlling the AT flow via the measured EIT data.

But how should one check the compatibility of the evolving conductivity to the measured data? One option would be to simulate the EIT measurement matrix at each time step for the evolving AT image and ensure that the distance to the measured matrix decreases. However, more robust control is provided by a novel concept called the CGO sinogram. It consists of the traces of the complex geometric optics (CGO) solutions (of the D-bar method) at the boundary, corresponding to low-frequency spectral parameters only. The CGO sinogram is stable to compute from EIT data . Furthermore, it appears to contain geometric information about the conductivity in a far more explicit form than the DN map (see Figure 4), at least for some of the most traditionally used current patterns, e.g., trigonometric current patterns.

While the AT flow will introduce edges, as desired, it also lowers the contrast in the image. To overcome this problem, we introduce a contrast-enhancement step. See Figure 3 for the result. We remark that our method provides a novel connection between the PDE-based inverse problems community and the PDE-based image processing community having a potentially strong impact on both.

The remainder of the paper is organized as follows. Sections 2-4 contain the theory behind the key pieces in the new edge preserving reconstruction method guided by the data-driven contrast enhancement of Section 5. In Section 2, the D-bar method is reviewed. The Ambrosio-Tortorelli functional used for recovering edges is described in Section 3. In Section 4, the novel CGO sinogram is introduced and the data-driven contrast enhancement is presented in Section 5. For the reader’s convenience, Section 6 is dedicated to an explicit description of the algorithm. In Section 7, the proposed algorithm is demonstrated on simulated noisy EIT measurements. A discussion of the results is given in Section 8 and the take home message and conclusions of the paper are given in Section 9.

## 2. A Brief Review of the D-bar Method

By the D-bar method we refer to the EIT algorithm based on the theory introduced in [33], first implemented in [38] and equipped with an explicit regularization step in [27]. Alternative D-bar methods have since emerged which handle complex coefficients [18, 19, 20], less regular conductivities [7, 26, 28] and merely bounded conductivities [4, 5, 3]. However, for the purpose of this article, we proceed with the well-established setting of [27], the only 2D approach with a proven regularization strategy.

The core idea in the D-bar method is to use a nonlinear Fourier transform tailor-made for EIT. To define the transform we need modified exponential functions, also called Complex Geometrical Optics (CGO) solutions.

Assume (for now) that and that in a neighborhood of the boundary. The constant non-unitary condition near can be dealt with by scaling as discussed below. The conductivity equation (1.1) can then be transformed, using the change of variables , to the Schrödinger equation

 [−Δ+q]ψ=0.

Here we define the potential by extending from to all of by setting for and writing

 q(z)=Δ√σ(z)√σ(z).

Note that has compact support in .

We introduce an auxiliary variable and look for CGO solutions satisfying

 [−Δ+q(⋅)]ψ(⋅,k)=0

with the asymptotic condition

 e−ikzψ(z,k)−1∈W1,p(R2).

We associate with by the mapping , so that denotes complex multiplication. For later use we introduce the related CGO solutions .

CGO solutions were originally introduced by Faddeev in [16] and later reinvented in the context of inverse problems in [39]. By [33] we know that the solutions exist and are unique for any for the potentials we consider here. (Other kinds of potentials may have exceptional points, or values with no unique . See [32] for more details.)

The regularized D-bar reconstruction algorithm is comprised of the following steps:

 Λδσ1⟶tR(k)2⟶σR(z).
• From boundary measurements to scattering data .
For each fixed , solve the Fredholm integral equation of the second kind for

 (2.1) ψδ(z,k)=eikz−∫∂ΩGk(z−ζ)[Λδσ−Λ1]ψδ(ζ,k)dS(ζ),

where is the DN map for the constant unit conductivity, and is the Faddeev Green’s function [16], with asymptotics matching , defined by

 Gk(z):=eikzgk(z),gk(z):=1(2π)2∫R2eiz⋅ξ|ξ|2+2k(ξ1+iξ2)dξ.

Evaluate the scattering transform for the cut-off frequency using the boundary traces from (2.1)

 (2.2) tR(k):={∫∂Ωei¯k¯z[Λδσ−Λ1]ψδ(z,k)dS(z),|k|
• From scattering data to conductivity .
Fix and solve the D-bar equation

 (2.3) ¯¯¯∂kμR(z,k)=14π¯ktR(k)e(z,−k)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯μR(z,k),

where , saving . The regularized conductivity is recovered by

 (2.4) σR(z)=(μR(z,0))2.

If near , but is instead a constant , the entire problem can be scaled as follows. Let denote a scaled conductivity which then has a value of 1 near . The corresponding scaled DN map is then computed by

 Λ~σ=σ0Λσ.

After recovering from Step 2 of the D-bar algorithm above, undo the scaling by multiplying by yielding the correct .

The algorithm described above has been used successfully on both simulated and experimental data [38, 29, 23, 27]. In many applications it is well known that high frequency features such as jump discontinuities and clear edges are present, e.g. organ boundaries. However, the low-frequency truncation needed in the regularized D-bar algorithm results in smoothed/blurred reconstructions where these high-frequency features are often absent. This calls for post-processing of the D-bar reconstruction to reintroduce the missing features. We propose minimizing a functional in a process known as Ambrosio-Tortorelli image segmentation.

## 3. Diffusive Image Segmentation

Consider the following image processing problem. We begin with a smooth image , which is the result of a blurring process applied to a clean, piecewise constant image . How can we recover from ?

In 1985, Mumford and Shah introduced the following functional for the purpose of detecting boundaries in general images [31]:

 (3.1) EMS(u,K)=∫Ω∖K|∇u|2dx+β∫Ω(u−˜u)2dx+α|K|,

where denotes a curve segmenting , the length of , and the two parameters are used for weighting the terms. The idea is to find the minimum of over images and curves ; the minimizing image is then considered an edge-preserving reconstruction of from .

As is unknown and singular, numerically minimizing (3.1) is a challenging task; in particular formulating a gradient-descent method with respect to is not straightforward. Therefore, Ambrosio and Tortorelli [2] proposed an elliptic approximation to (3.1) by introducing an edge-strength function for controlling the gradient of . The Ambrosio-Tortorelli (AT) functional is defined by

 (3.2) EAT(u,v)=∫Ωβ(u−˜u)2+v2|∇u|2+α(ρ|∇v|2+(1−v)24ρ)dx.

The additional parameter specifies, roughly speaking, the edge width of . Then -converges to as [2], which can be understood as solutions of (3.2) converge to solutions of (3.1) with the parameter , for further explanations see [9].

The advantage of (3.2) over (3.1) is that the minimizer can be obtained by an artificial time evolution formulated via a coupled PDE as the gradient-descent equations with an imposed homogeneous Neumann boundary condition:

 (3.3)

The equations in (3.3) are referred to as the AT flow. With these gradient descent equations, implementing the numerical minimization algorithm is now a straightforward task (use finite differences or a parabolic finite element solver). Numerous modifications of the functionals (3.1) and (3.2) have been proposed, e.g. substituting the squared norm with in the spirit of Total Variation [1], or adjusting the auxiliary function as in [15].

The effect of (3.2) can be visualized as follows. First, keeping fixed, one sees that the first equation (3.3) minimizes the functional

 ∫Ωβ(u−˜u)2+v2|∇u|2dx=β∥u−˜u∥2L2(Ω)+∥v|∇u|∥2L2(Ω).

The representation on the right hand side is a typical form for stating regularization problems. The first fidelity term controls the distance of the obtained solution from the initial reconstruction, whereas the regularization term weights with respect to the positive (here fixed) edge-strength function . The balance between the fidelity and regularization terms can be controlled by the regularization parameter .

In the same manner, keeping fixed, from the second equation in (3.3), we have a minimizer for a functional that can be written as

 ∫Ω4ρ|∇v|2+1+4ρα|∇u|2ρ⎛⎝11+4ρα|∇u|2−v⎞⎠2dx.

Here we clearly see the auxiliary variable can be interpreted as a smooth approximation to the Perona-Malik filter [34], given by

 (3.4) g(|∇u|2)=11+4ρα|∇u|2.

To state a minimizing algorithm on the equations (3.3), one needs to specify an initial guess of . This is where the interpretation (3.4) comes in handy and we can set the first approximation for the auxiliary variable as .

Perona and Malik based their edge-aware smoothing on the diffusion equation

 (3.5) ∂tu=div(g(|∇u|2)∇u).

Historically, the model in (3.5) has a strong effect on noise removal and the smoothing of images, by keeping edges stable for a long time in the process. However, one downside of this particular diffusion concept is that the limiting function is a constant and hence a stopping criterion during the iteration is needed [40]. In contrast, the AT functional in (3.2) converges to the Mumford-Shah segmentation functional (3.1), for which minimizers are known to be piecewise constant with respect to the discontinuity set [13]. By using the AT functional instead, the problem of defining a proper stopping criterion is shifted to choosing the correct parameters and , and making it possible to adapt the minimization problem to the reconstruction task at hand.

## 4. The CGO Sinogram

The central idea of this study is to complement the regularized D-bar method by applying an edge-introducing image processing algorithm to a blurred EIT reconstruction. However, when manipulating the EIT image we need to ensure that we are improving the image. This is accomplished by monitoring the resulting error in a novel way (see Figure 3).

The obvious approach is to monitor the data discrepancy , where is the measured data and denotes the enhanced reconstruction, but the basic observation behind our new data fidelity term is that, for the commonly used trigonometric current patterns, the DN map encodes geometric information nonlinearly in a very complicated, non-intuitive, and often unstable way. For example, take a look at the three conductivities shown in the left column of Figure 4. They each have one circular inclusion of conductivity two embedded into a homogeneous background of unit conductivity. The middle column shows the matrix approximations to the DN map. Can you deduce the location of the inclusion from the DN map? We didn’t think so.

We recall from Section 2 the related CGO solutions with asymptotic behaviour

 μ(z,k)−1∈W1,p(R2).

Now take a look at the right column in Figure 4. There we show the absolute value of the difference of the CGO solution and its limiting value 1, i.e. , as a function of the spectral frequency angles (vertical axis) and physical angles (horizontal axis) where both angles range from to . The locations of the inclusions are immediately discernible!

Clearly, Figure 4 is just one very simple example. See Figure 5 for a more complicated situation, where we add a circular inclusion to a heart-and-lungs phantom. The location of the inclusion is clearly indicated in the difference of the two CGO sinograms.

We believe that the numerical evidence presented in Figures 4 and 5 reflects a more general fact: calculating the CGO sinogram

 (4.1) Sσ(θ,φ,r):=μ(eiθ,reiφ)−1=exp(−irei(φ+θ))ψ(eiθ,reiφ)−1

from the EIT measurements of is stable (for inside the region of proven stability) and encodes the geometry of the conductivity in a useful and more transparent way.

## 5. Contrast Enhancement

The Ambrosio-Tortorelli (AT) segmentation flow, discussed in Section 3 above, can transform a blurry image into a sharper image. However, the AT flow comes with a reduction in contrast, which in turn is a key benefit of EIT imaging. To overcome this obstacle, we propose using a data-driven contrast enhancement technique kept in check by the CGO sinogram.

We search for a contrast enhanced (corrected) , such that the values are stretched or damped according to the resulting error in CGO sinogram. We utilize a two parameter model and consider values greater and lower 1 independently, where 1 is the value near the boundary. For that let

 f(z):=σ\tiny AT(z)−1,

denote the difference between the conductivity after AT flow and the constant 1. Recall the a priori constants and known to bound the conductivity . In practice, such ballpark bounds are readily available. For example, in chest imaging, using an applied current with frequency 100 kHz, internal conductivities range from around Siemens/meter (e.g. fat) to Siemens/meter (e.g. heart tissue) [12].

Denote by and the following minimum and maximum values

 m:=minz∈Ωf(z),M:=maxz∈Ωf(z),

assuming these values are nonzero. Define the scaling function for the scaling parameter by

 fs,t(z)=⎧⎨⎩sf(z)m(c−1)for z satisfying f(z)<0% tf(z)M(C−1)for z satisfying f(z)≥0

and set

 (5.1) σs,t(z):=1+fs,t(z).

The optimal contrast enhanced conductivity within the bounds and is then determined by minimizing the CGO sinogram error over the scaling parameter , i.e.

 (5.2) (s0,t0):=argmin(s,t)∈[0,1]2{∥Sδσ(θ,φ,r)−Sσs,t(θ,φ,r)∥2L2(T2)},

where is the CGO sinogram corresponding to the noisy measurement of , and is computed from . Then (5.1) is used to define

 σ\tiny CE(z):=σs0,t0(z).

Various algorithms are available to find an approximation to the optimal and in (5.2). In this introductory work we apply a derivative free pattern search algorithm which is known to be stable in the presence of noise: the DIRECT (DIviding RECTangles) pattern search for global optimization [24, 17].

In Figure 6 we see clear evidence that minimizing with respect to the CGO sinogram is advantageous to the DN map for preserving important image features. The figure shows the heart and lungs test problem with conductivity reconstructed by the D-bar algorithm applied to noisy EIT measurements (left) and the optimal contrast-enhanced images chosen by minimizing the error of CGO sinogram (middle) and DN map (right) to the measured data. Clearly, the image chosen by the CGO sinogram minimization more accurately portrays the original than the one chosen by the DN map. Even though both solutions were able to minimize the -error to the true phantom, from the introductory example in Figure 1, the features preserved differ immensely (i.e. the heart is barely visible in the DN guided minimization).

## 6. A Data-Driven Edge-Preserving D-bar Algorithm

The aim of this paper is to combine the strengths of each the methods described above in Sections 2-4, with the data-driven contrast enhancement of Section 5. We compute the reconstruction of the conductivity with the regularized D-bar method and reintroduce edges by a data-driven post-processing of the image. The post-processing is monitored by the CGO sinogram error which incorporates geometric information about the reconstruction. We thus propose the following Data-Driven Edge-Preserving D-bar Algorithm:

• Fix and compute the regularized D-bar reconstruction using the D-bar algorithm described in Section 2.

• Fix a radius located in the stable disc: .

1. Compute the CGO sinogram from the original noisy data by solving the boundary integral equation (2.1) for and setting

 Sδσ(θ,φ,r)=μδ(eiθ,reiϕ)−1.
2. Calculate the CGO sinogram for the D-bar image.

3. Record the relative CGO sinogram error

 E0:=∥Sδσ(θ,φ,R)−SσR(θ,φ,R)∥L2(T2)∥Sδσ(θ,φ,R)∥L2(T2).
• Reintroduce edges to via minimization of the Ambrosio-Tortorelli functional by solving the gradient descent equations (3.3).

1. Initialize the constants and .

2. Calculate the initial approximation for the auxiliary function as defined by (3.4).

3. Choose a time step and begin solving (3.3) iteratively.

• Check the flow every , e.g. , time steps as follows. For the -th check:

1. Denote the image to be checked .

2. Determine , the optimal contrast enhanced version of using the CGO sinogram optimization method described in Section 5.

3. Record the CGO sinogram error for .

4. If , return to the AT flow with (the non-contrast enhanced version) and repeat steps (i)-(iii). If not, or if a maximum number of iterations is reached, set

 σNEW(z):=σ\tiny CEj(z)

and the algorithm is complete.

## 7. Computational results

We tested the algorithm on simulated noisy EIT measurement data for test cases of potential interest for the medical and industrial communities.

To simulate the EIT measurements, we solved the Neumann problem corresponding to the conductivity equation

 ∇⋅σ∇u=0,z∈Ω⊂R2σ∂u∂ν=ϕj,z∈∂Ω,

for , representing 32 linearly independent current patterns. In this paper is the unit disc, and we used the Fourier basis functions

 ϕj(eiθ)=1√2eijθ,z=eiθ∈∂Ω.

The matrix approximation to the DN maps and (the DN map corresponding to a uniform unit conductivity) were computed for each test problem using the standard methods described in [30]. The numerical implementation of the regularized D-bar method was first fully described in [27] and later explained in more detail in [30, Section 13.2], including freely available Matlab code.

A discussion of a numerical implementation with finite differences of the AT flow can be found in [15]. In this study, we used the Matlab integrated PDE solver for the equations (3.3) and for each test problem below the choice of parameters is given.

Additional zero-mean random Gaussian noise was added to the DN matrix so that

 (7.1) ∥Λσ−Λδσ∥H1/2(∂Ω)→H−1/2(∂Ω)≤δ

using the methods described in [27]. While types of noise and their magnitudes vary among EIT devices, a benchmark of around 0.01% has been obtained [10].

The resulting CGO sinogram and conductivity errors stated below are given as relative errors. The relative error of CGO sinograms corresponding to the original conductivity and to some reconstruction was measured by

 (7.2) ∥Sδσ(θ,φ,R)−Sσd(θ,φ,R)∥L2(T2)∥Sδσ(θ,φ,R)∥L2(T2).

Similarly we measured the relative error of the reconstructed conductivities to the known original via

 (7.3) ∥σ−σd∥Lp(Ω)∥σ∥Lp(Ω)

with or . We note that (7.2) can be measured for real data cases, contrary to (7.3) for which the knowledge of the correct conductivity is needed.

### 7.1. A Heart-and-Lungs Phantom

The leftmost image in Figure 7 shows our piecewise constant phantom. The lungs have conductivity 0.5, the heart has conductivity 2, and the background has unit conductivity.

The D-bar reconstruction was obtained using a truncation radius of in the scattering data. The D-bar method allows computing the reconstruction at arbitrary points in the -plane, and we chose to reconstruct directly at the points of the FEM mesh (of 33025 elements) to be used in the AT flow. We have added additional noise of amplitude to the DN map satisfying (7.1), and corresponding to noise far surpassing the benchmark for measurement noise. The D-bar reconstruction is shown in the middle image in Figure 7; it has a relative CGO sinogram error of 20.3%.

The AT flow (Step 3) was computed with parameters

 α=200,β=0.1,ρ=0.1,

and the contrast enhanced solution was computed every fifth iteration using the CGO sinogram with , well within the observed reliable region for the regularized D-bar method. The boundary constants were roughly chosen as and . A summary of the obtained solutions is illustrated in Figure 8.

The optimal solution was obtained at iteration 45 with a relative error in CGO sinogram of 17.28%. Figure 9 displays the error in CGO sinogram as well as the error of reconstruction to the true conductivity throughout the evolution of the AT flow.

### 7.2. An Industrial Pipe Phantom

The second test case is an example from the oil industry. It roughly models a pipe with oil (top layer, conductivity 1.2), water (middle layer, conductivity 2.0), and sand (bottom layer, conductivity 0.3) similar to the test problem of [3].

The initial D-bar reconstruction is computed on the same mesh with a truncation radius of and additional noise of in the measured DN map. The parameters for the AT flow and contrast enhancement were chosen the same as in the previous example. The CGO sinograms were computed with a small radius of far within the reliable region.

The minimization of the CGO sinogram error with the proposed algorithm did not converge and hence the algorithm stopped at the maximum number of iterations chosen as .

The error in the CGO sinogram has been halved from in the original D-bar reconstruction, to in the improved reconstruction obtained with the new algorithm (left plot of Figure 11). The right plot shows the errors for contrast enhanced versions of the conductivities after the AT flow. Their error increased, after an initial decrease, from of the D-bar reconstruction to of the reconstruction obtained by the AT flow after 200 iterations and of the corresponding contrast enhanced solution. However, we point out that in practice such a comparison (right plot of Figure 11) is not possible as the true conductivity is unknown, and only the left plot is possible, where in fact we see a clear decrease in the CGO error.

## 8. Discussion

The heart-and-lungs phantom of Section 7.1 is of particular interest, since it represents, an admittedly simplified version of, a major application of electrical impedance tomography in the medical field: monitoring the blood and air flow in a patient’s heart and lungs. With the proposed algorithm we were able to clearly distinguish the left and right lung as well as introduce clear edges, see Figure 7. As discussed above, a level of noise is reasonably high for EIT measurements and hence gives a good impression of how the algorithm behaves with noise corrupted data. The original D-bar recovered conductivity (see the middle image in Figure 7) has good contrast but lacks sharpness, as is typical for D-bar reconstructions.

A summary of the AT flow is illustrated in Figure 8. Notice that the AT flow gradually divides the two lungs, and that the separated areas converge to constant values. As seen in the right plot of Figure 9, the AT flow (blue line) clearly minimizes the -error of the evolved conductivity to the true conductivity. Furthermore, the behavior of its corresponding CGO sinogram error (blue plus sign markers in the left plot), initially decreasing but followed by a sharp increase, reinforces the need for a contrast enhancement step. Note that the CE solutions have a similar behavior in both CGO sinogram error and the reconstruction error of the conductivity, and that the CE evolved conductivities have smaller relative errors to the true conductivity than their non CE counterparts.

The second test case (Section 7.2) is an example from the oil industry. When measuring a pipe with oil (top layer), water (middle layer), and sand (bottom layer) one wants to know how much oil is transported in the pipe, making it important to distinguish the levels of each substance clearly, i.e. their edges. The new reconstruction has clear and sharp edges dividing the different substances and one can tell easily how much of each substance is present in the pipe. As one can see in Figure 10, the structure of the pipe has been reintroduced to the reconstruction, delivering a realistic view of the imaged target.

The right plot in Figure 11 shows that the error in reconstructed conductivities did not decrease as nicely as in the Heart-and-Lungs phantom. In fact, the relative conductivity reconstruction error increased from of the D-bar reconstruction to of the contrast enhanced solution of the last iteration 200 in the AT flow. A reason for the higher error can be seen in Figure 10, the contrast enhancement increases the conductivity of the top layer as well as the middle layer, which produces a higher error in the middle layer. This suggests that more sensitive models for the contrast enhancement may be needed. Nevertheless, the information contained in the contrast enhanced solution is far more useful for evaluation of the target and suggests that the CGO sinogram contains more information about the reconstruction’s geometry.

## 9. Conclusions

A novel edge-preserving D-bar method with a data-driven contrast enhancement was introduced and tested on simulated EIT measurement data. The algorithm works as advertised by both sharpening and enhancing the contrast of the reconstruction, even in the presence of additional noise added to the Dirichlet-to-Neumann boundary measurements.

Key to the approach is the invention of the CGO sinogram, a more reliable and geometrically transparent quantity than the DN map. The CGO sinogram provides an important breakthrough towards new uses of this nonlinear data, and based on our findings we are excited what further theoretical analysis on the stability of the CGO sinogram will reveal.

Our results also have implications outside the realm of D-bar methods. The traditional regularization approach for EIT reconstructions is to find the minimizer of a functional of the form

 (9.1) ∥Λδσ−Λσ′∥Y+α∥σ′∥X,

where the norm corresponds, for instance, to Tikhonov or Total Variation regularization, and is a regularization parameter. Replacing the traditional data fidelity term in (9.1) by an analogous term based on the CGO sinogram leads to the functional

 (9.2) ∥Sσ(θ,φ,r)−Sσ′(θ,φ,r)∥2L2(T2)+α∥σ′∥X,

where denotes the two-dimensional torus. Based on the evidence seen in Figures 4 and 5, we strongly suspect that using (9.2) would lead to superior reconstructions compared to (9.1). A similar comment applies to the likelihood distributions used in Bayesian inversion approaches for EIT.

## Acknowledgments

The study was supported by the SalWe Research Program for Mind and Body (Tekes - the Finnish Funding Agency for Technology and Innovation grant 1104/10) and by the Academy of Finland (Finnish Centre of Excellence in Inverse Problems Research 2012–2017, decision number 250215).

## References

• [1] R. Alicandro, A. Braides, and J. Shah, Approximation of non-convex functionals in GBV, 1998.
• [2] L. Ambrosio and V. M. Tortorelli, Approximation of functionals depending on jumps by elliptic functionals via -convergence, Communications on Pure and Applied Mathematics, 43 (1990), pp. 999–1036.
• [3] K. Astala, J. Mueller, L. Päivärinta, A. Perämäki, and S. Siltanen, Direct electrical impedance tomography for nonsmooth conductivities, Inverse Problems and Imaging, 5 (2011), pp. 531–549.
• [4] K. Astala and L. Päivärinta, A boundary integral equation for Calderón’s inverse conductivity problem, in Proc. 7th Internat. Conference on Harmonic Analysis, Collectanea Mathematica, 2006.
• [5] K. Astala and L. Päivärinta, Calderón’s inverse conductivity problem in the plane, Annals of Mathematics, 163 (2006), pp. 265–299.
• [6] J. Bikowski, K. Knudsen, and J. L. Mueller, Direct numerical reconstruction of conductivities in three dimensions using scattering transforms, Inverse Problems, 27 (2011), p. 015002.
• [7] R. M. Brown, Global uniqueness in the impedance imaging problem for less regular conductivities, SIAM Journal on Mathematical Analysis, 27 (1996), pp. 1049–1056.
• [8] A.-P. Calderón, On an inverse boundary value problem, in Seminar on Numerical Analysis and its Applications to Continuum Physics (Rio de Janeiro, 1980), Soc. Brasil. Mat., Rio de Janeiro, 1980, pp. 65–73.
• [9] A. Chambolle, Image segmentation by variational methods: Mumford and shah functional and the discrete approximations, SIAM Journal on Applied Mathematics, 55 (1995), pp. 827–863.
• [10] M. Cheney, D. Isaacson, and J. C. Newell, Electrical impedance tomography, SIAM Review, 41 (1999), pp. 85–101.
• [11] H. Cornean, K. Knudsen, and S. Siltanen, Towards a -bar reconstruction method for three-dimensional eit , Journal of Inverse and Ill-Posed Problems, 14 (2006), pp. 111–134.
• [12] I. N. R. Council, Dielectric properties of body tissues, 2013.
• [13] E. De Giorgi, M. Carriero, and A. Leaci, Existence theorem for a minimum problem with free discontinuity set, Archive for Rational Mechanics and Analysis, 108 (1989), pp. 195–218.
• [14] F. Delbary, P. C. Hansen, and K. Knudsen, Electrical impedance tomography: 3D reconstructions using scattering transforms, Applicable Analysis, 0 (0), pp. 1–19.
• [15] E. Erdem and S. Tari, Mumford-shah regularizer with contextual feedback, Journal of Mathematical Imaging and Vision, 33 (2009), pp. 67–84.
• [16] L. D. Faddeev, Increasing solutions of the Schrödinger equation, Soviet Physics Doklady, 10 (1966), pp. 1033–1035.
• [17] D. E. Finkel, Direct optimization algorithm user guide, tech. rep., Center for Research in Scientific Computation, North Carolina State University, 2003.
• [18] E. Francini, Recovering a complex coefficient in a planar domain from Dirichlet-to-Neumann map, Inverse Problems, 16 (2000), pp. 107–119.
• [19] S. Hamilton, C. Herrera, J. L. Mueller, and A. VonHerrmann, A direct D-bar reconstruction algorithm for recovering a complex conductivity in 2-D, Inverse Problems, 28 (2012), p. 095005.
• [20] S. J. Hamilton and J. L. Mueller, Direct EIT reconstructions of complex admittivities on a chest-shaped domain in 2-D, IEEE Transactions on Medical Imaging, 32 (2013), pp. 757–769.
• [21] L. Harhanen, N. Hyvönen, H. Majander, and S. Staboulis, Edge-enhancing reconstruction algorithm for three-dimensional electrical impedance tomography, ArXiv e-prints, arXiv:1406.1279 (2014).
• [22] T. Helin and M. Lassas, Hierarchical models in statistical inverse problems and the mumford–shah functional, Inverse problems, 27 (2011), p. 015008.
• [23] D. Isaacson, J. Mueller, J. Newell, and S. Siltanen, Imaging cardiac activity by the D-bar method for electrical impedance tomography, Physiological Measurement, 27 (2006), pp. S43–S50.
• [24] D. R. Jones, C. D. Perttunen, and B. E. Stuckman, Lipschitzian optimization without the lipschitz constant, Journal of Optimization Theory and Applications, 79 (1993), pp. 157–181.
• [25] M. Jung, X. Bresson, T. F. Chan, and L. A. Vese, Nonlocal Mumford-Shah regularizers for color image restoration, IEEE Trans. Image Process., 20 (2011), pp. 1583–1598.
• [26] K. Knudsen, A new direct method for reconstructing isotropic conductivities in the plane, Physiological Measurement, 24 (2003), pp. 391–403.
• [27] K. Knudsen, M. Lassas, J. Mueller, and S. Siltanen, Regularized D-bar method for the inverse conductivity problem, Inverse Problems and Imaging, 3 (2009), pp. 599–624.
• [28] K. Knudsen and A. Tamasan, Reconstruction of less regular conductivities in the plane, Communications in Partial Differential Equations, 29 (2004), pp. 361–381.
• [29] J. Mueller and S. Siltanen, Direct reconstructions of conductivities from boundary measurements, SIAM Journal on Scientific Computing, 24 (2003), pp. 1232–1266.
• [30] J. Mueller and S. Siltanen, Linear and Nonlinear Inverse Problems with Practical Applications, vol. 10 of Computational Science and Engineering, SIAM, 2012.
• [31] D. Mumford and J. Shah, Boundary detection by minimizing functionals, in IEEE Conference on Computer Vision and Pattern Recognition, 1985.
• [32] M. Music, P. Perry, and S. Siltanen, Exceptional circles of radial potentials, Inverse Problems, 29 (2013), p. 045004.
• [33] A. I. Nachman, Global uniqueness for a two-dimensional inverse boundary value problem, Annals of Mathematics, 143 (1996), pp. 71–96.
• [34] P. Perona and J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Transactions on Pattern Analysis and Machine Intelligence, 12 (1990), pp. 629–639.
• [35] R. Ramlau and W. Ring, A mumford–shah level-set approach for the inversion and segmentation of x-ray tomography data, Journal of Computational Physics, 221 (2007), pp. 539–557.
• [36] L. Rondi and F. Santosa, Enhanced electrical impedance tomography via the Mumford-Shah functional, ESAIM Control Optim. Calc. Var., 6 (2001), pp. 517–538.
• [37] J. Shah, A common framework for curve evolution, segmentation and anisotropic diffusion, in IEEE Conference on Computer Vision and Pattern Recognition, IEEE, 1996, pp. 136–142.
• [38] S. Siltanen, J. Mueller, and D. Isaacson, An implementation of the reconstruction algorithm of a. nachman for the 2-d inverse conductivity problem, Inverse Problems, 16 (2000), pp. 681–699.
• [39] J. Sylvester and G. Uhlmann, A global uniqueness theorem for an inverse boundary value problem, Annals of Mathematics, 125 (1987), pp. 153–169.
• [40] J. Weickert, Anisotropic Diffusion in Image Processing, Teubner Stuttgart, 1998.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters