3D CGO-Based Experimental EIT]Fast Absolute 3D CGO-Based Electrical Impedance Tomography on Experimental Tank Data

Hamilton et al.]S. J. Hamilton, P. A. Muller, D. Isaacson,  V. Kolehmainen,  J. Newell,  O. Rajabi Shishvan,  G. Saulnier, and  J. Toivanen

S. J. Hamilton is with the Department of Mathematical and Statistical Sciences; Marquette University, Milwaukee, WI 53233 USA,  email:
D. Isaacson is with the Department of Mathematics, Rensselaer Polytechnic Institute, Troy, NY 12180, USA
V. Kolehmainen and J. Toivanen are with the Department of Applied Physics, University of Eastern Finland, FI-70210 Kuopio, Finland
P. A. Muller is with the Department of Mathematics & Statistics; Villanova University, Villanova, PA 19085 USA
J. Newell is with the Department of Biomedical Engineering, Rensselaer Polytechnic Institute, Troy, NY 12180, USA
O. Rajabi Shishvan and G. Saulnier are with the Department of Electrical and Computer Engineering, University at Albany - SUNY, Albany, NY 12222, USA
J. Toivanen is with the Department of Applied Physics, University of Eastern Finland, FI-70210 Kuopio, Finland

Objective: To present the first 3D CGO-based absolute EIT reconstructions from experimental tank data. Approach: CGO-based methods for absolute EIT imaging are compared to traditional TV regularized non-linear least squares reconstruction methods.  Additional robustness testing is performed by considering incorrect modeling of domain shape.  Main Results: The CGO-based methods are fast, and show strong robustness to incorrect domain modeling comparable to classic difference EIT imaging and fewer boundary artefacts than the TV regularized non-linear least squares reference reconstructions. Significance:  This work is the first to demonstrate fully 3D CGO-based absolute EIT reconstruction on experimental data and also compares to TV-regularized absolute reconstruction. The speed (1-5 seconds) and quality of the reconstructions is encouraging for future work in absolute EIT.

§ INTRODUCTION

The main objective of this paper is to demonstrate the feasibility, speed, and robustness of producing 3D absolute (static) images of the electrical conductivity inside a tank, from experimental Electrical Impedance Tomography (EIT) data measured on an array of electrodes on the tank’s surface by the ACT5 <cit.> adaptive current tomography system, using CGO-based reconstruction algorithms.  The primary contributions of this work are that it presents the first use of Complex Geometrical Optics (CGO) based methods to produce absolute (static) images of the conductivity from experimentally measured EIT data in 3D, and studies their robustness under modeling errors.

“Absolute", or “static,” conductivity images are images of the conductivity inside a body made from one set of experimental measurements made on the surface of the body at one time <cit.>. Alternatively, “dynamic," or “time-difference," imaging uses two sets of data measured at two different times to make an image of the change in the conductivity that took place between the two times the measurements were made. We present both types of images made from experimental data by CGO-based algorithms in 3D in this paper.

Applications of EIT began with geophysical exploration over a century ago <cit.> and have since expanded into several fields including the transport of fluids and gases <cit.>, and biomedical imaging <cit.>.  Systems for monitoring lung function in real-time are now commercially available and clinical trials are in progress to determine the extent to which some of these systems might be used to guide mechanical ventilation <cit.>. These systems typically make and display two dimensional images of the differences in conductivity between two states, e.g., lungs filled with air and lungs depleted of air, in order to take advantage of dynamic (difference) imaging methods and algorithms <cit.>.  The desire to improve the diagnosis of cancer and stroke has motivated the development of systems and methods capable of imaging the absolute or static internal conductivity and permittivity in 3D <cit.>.  The interested reader is referred to <cit.> for further review of applications of EIT.

Most absolute (static) EIT reconstruction focuses on solving a simplified linearized problem or iteratively solving an optimization-based method which requires repeated solutions to the forward problem which can become costly for highly dense meshes.  CGO-based methods are direct methods in that they do not require iteration. They have the capability to solve the full nonlinear mathematical inverse problem, and do not require repeated solutions to the forward problem, e.g. via the finite element method (FEM).  The D-bar inversion algorithm, which is a specific type of CGO-based inversion algorithm in 2D, has been used to make both absolute and difference images in 2D in real-time from experimental data measured on tanks, and for difference imaging on human subjects in laboratory and clinical settings.  See <cit.> for a recent review of the 2D D-bar method and its applications, and <cit.> for the theoretical foundation.

In 3D, existence and uniqueness of solutions <cit.> can be shown for a 3D D-bar type equation but the constructive proof, upon which the reconstruction algorithms are built, bypasses the 3D D-bar equation instead using a high (non-physical) frequency limit connecting the nonlinear scattering data and the linear Fourier data.  Advances in the numerical implementation of 3D CGO-based methods are more recent. The first numerical implementation of 3D CGO-based methods on simulated electrode data using current injection on the surface of the sphere was presented in <cit.>.  A 3D CGO-based inversion algorithm, regularization scheme, rigorous proof of stability under certain hypotheses, and examples reconstructing the conductivity inside a sphere from numerically simulated Dirichlet data on the entire surface of the sphere (without electrodes) were given in <cit.>. Until now it has been an open question if 3D CGO-based algorithms could be used with experimental data. This paper answers this question affirmatively by showing that the conductivity can be recovered rapidly and robustly from experimental data measured on 32 electrodes on the surface of a rectangular prism building on the work of <cit.>.

It was demonstrated in <cit.> that it is possible to make electrical impedance tomography (EIT) Systems that produce both 3D static and dynamic images of the interior of the chest showing heart and lung regions, as well as changes in those regions due to ventilation and perfusion. These systems used linearized and iterative methods. The former are fast but less accurate and the latter are slower. Recent 3D CGO inversion algorithms and their analysis, when applied to synthetic or simulated data, suggested that they have the potential to be fast and more accurate  <cit.>.  Here we test their capabilities on experimental EIT data.

In this paper we compare the following methods for making static 3D images from experimental tank data, collected by the ACT5 imaging system:

* A CGO linearization method based on Calderón’s original proposal introducing CGO ideas into the subject <cit.>. We call this algorithm Calderón’s method.

* Two CGO-based methods for solving the full non-linear inversion problem based on Sylvester and Uhlmann’s original uniqueness proof and Nachman’s constructive uniqueness proof in <cit.>.  The first is the $\texp$ algorithm, and the second is called $\tzero$, <cit.> both of which are simplifications of the original constructive proof.

* A more traditional iterative inversion method based on optimization using a Total Variation (TV) regularization.

The EIT data were measured on 32 electrodes attached to the six surfaces of a rectangular tank using the ACT5 system. The system can adaptively determine patterns of currents to apply to all 32 electrodes that result in voltages on those electrodes that are proportional to the applied currents. These “eigen currents" form a discrete orthogonal set and provide improved voltage data for reconstructing the conductivity inside the tank from a limited amount of current or power that can be applied to the tank. The theory of adaptive current tomography systems is given in <cit.>.

The remainder of the paper is organized as follows.  Section <ref> provides a brief review of the mathematics of EIT and encompasses the methods used in the work: the CGO and reference reconstruction algorithms, experimental setup for the ACT5 data collection, robustness tests that will be explored, and metrics used to evaluate the results.  The results section, Section <ref>, presents slice, 3D, and isosurface renderings of the recovered conductivities for correct and incorrect domain modeling.  Section <ref> contains a discussion of the results and conclusions are drawn in Section <ref>.

§ METHODS

In this work we compare three CGO-based methods (Calderón's Method, the $\texp$ Method, and the $\tzero$ Method) to a more common iterative Total Variation (TV) regularized non-linear least squares method for absolute EIT image reconstruction.  For comparison, we also include time difference EIT images for the CGO-based methods and compare them to a linearized difference imaging method.  We begin with a brief review of the mathematical problem.

§.§ Mathematical Background

The mathematical problem of reconstructing the internal conductivity, when measurements can be made with infinite precision everywhere on the boundary of a body, is currently called “Calderón’s Problem” by much of the mathematical community since A. Calderón formulated this inverse problem as follows <cit.>: In two or more dimensions, can one find the conductivity, $\sigma(x)$, inside a body, $\Omega$, from all possible electrical measurements made on the surface, $\bndry$, of the body? Here the voltage or potential,  $u(x)$, inside the body, due to an applied surface current density,  $j(x)$, is assumed to satisfy the following low frequency approximation to Maxwell’s Equations, which we will refer to as the conductivity equation, with a Neumann boundary condition:

\begin{eqnarray}
\nabla\cdot\sigma(x)\nabla u(x) &=& 0, \qquad x\in\Omega\subset \R^3 \label{eq:cond}\\
\sigma(x)\frac{\partial u(x)}{\partial \nu(x)} &=& j(x), \qquad x\in\bndry.\nonumber
\end{eqnarray}

Here $\nu =\nu(x)$ denotes the outward pointing unit normal to the surface at $x\in\bndry$, and, $\frac{\partial u(x)}{\partial \nu(x)}=\nu(x)\cdot\nabla u(x)$.  We denote the mapping from applied current density to resulting voltage on the surface, called the “Neumann-to-Dirchlet” (ND) map, by $\mathcal{R}_\sigma$, where $\mathcal{R}_\sigma j\equiv u(x)$ for $x\in\bndry$, and $u(x)$ solves the conductivity equation (<ref>) with Neumann data $j(x)$.

The Calderón problem can also be formulated as the mathematical problem.  If one is given the ND map or operator, $\mathcal{R}_\sigma$, or equivalently its inverse, the Dirichlet-to-Neumann (DN) operator,  $\Lambda_\sigma:=\mathcal{R}_\sigma^{-1}$,  can one find $\sigma(x)$? Calderón showed that if $\sigma(x)$ does not differ too much from a constant then one can recover an approximation to it from the ND map. His short paper showed this could be done using Fourier Transforms in a very clever way, by introducing special solutions, $e^{\zeta\cdot x}$, to the conductivity equation with constant conductivity, i.e. the Laplace equation, where $\zeta$ is a complex valued vector with $\zeta\cdot\zeta=0$. This is possible in two or more dimensions when the conductivity is close to a constant and gave birth to Calderón's method as described in this paper, as well as more powerful methods for solving the full non-linear problem by generalizing Calderón’s solutions to what we now call Complex Geometrical Optics (CGO) solutions, introduced by <cit.> in their landmark paper proving that the Calderón problem in three or more dimensions has a unique solution, where the inversion problem is reduced to a Fourier transform in the limit $|\zeta|\to\infty$. The $\texp$ and $\tzero$ methods described below in  <ref>  follow this strategy, as opposed to the 2D D-bar methods described in <cit.>. This strategy is not possible in 2D where the CGO solutions are found by solving a first order linear PDE involving D-bar operators in the complex vector $\zeta$ and taking it to 0. Constructive proofs using these CGO solutions and D-bar ideas from inverse scattering theory, along with applications to scattering and acoustics were given in <cit.>.

CGO methods were used to prove uniqueness in the more difficult case of 2D where constructive methods for reconstructing the conductivity were given in detail in <cit.> for $\sigma\in\mathcal{C}^2(\Omega)$ and later in <cit.> for $\sigma\in{L}^\infty(\Omega)$. Other pioneering work proving uniqueness under a variety of hypotheses on the conductivity include <cit.>. Extensive references to more recent progress in the analysis and numerical analysis of the Calderón problem can be found in <cit.>; <cit.>.

In what follows we will be interested in the problem of reconstructing an approximation to the internal conductivity from finitely many experimental measurements made with finite precision. Unfortunately, this is an ill-posed problem and, unlike the purely mathematical problem, it does not have a unique solution. Nevertheless, it is sometimes possible to reconstruct useful approximations to the internal conductivity with a finite number of degrees of freedom, or voxels, which we will illustrate by making images from experimental data and comparing them to the actual interior conductivity within the tanks.

The EIT problem for a body with $L$ electrodes, $e_\ell$, $\ell=1,2,\ldots,L$, on its surface, is to find an approximation to the internal conductivity from all the possible electrical measurements made on these $L$ electrodes. In particular we will assume that we apply $L-1$ linearly independent patterns of currents, $\vec{I}^{(k)}$, $k=1,2,\ldots, L-1$, to the $L$ electrodes, and measure the resulting $L-1$ voltage patterns, $\vec{V}^{(k)}$, where $\vec{I}^{(k)}_\ell$ and $\vec{V}^{(k)}_\ell$ denote the applied current, and measured voltage from the $k^{\mbox{\tiny th}}$ pattern, on the $\ell^{\mbox{\tiny th}}$ electrode, for $\ell = 1,2,\ldots,L$.  From conservation of charge, and our choice of ground, we assume

\[\sum_{\ell=1}^L \vec{I}^{(k)}_\ell = \sum_{\ell=1}^L \vec{V}^{(k)}_\ell = 0.\]

The $k^{\mbox{\tiny th}}$ voltage or potential, $u^{(k)}(x)$, resulting inside the body is determined by the conductivity equation, $\nabla\cdot\sigma \nabla u^{(k)}=0$, and the “Complete Electrode Model”, <cit.> where:

\begin{align}
\int_{e_\ell} \sigma(x)\frac{\partial u^{(k)}(x)}{\partial \nu(x)} dS(x) &=I^{(k)}_\ell, \qquad x\in e_\ell\nonumber\\
\sigma(x)\frac{\partial u^{(k)}(x)}{\partial \nu(x)}  & = 0, \qquad x\notin \bigcup_{\ell=1}^L e_\ell\label{eq:CEM}\\
u^{(k)}(x) + z_\ell \sigma(x)\frac{\partial u^{(k)}(x)}{\partial \nu(x)} & = V^{(k)}_\ell, \qquad x\in e_\ell.\nonumber
\end{align}

Here $z_\ell$ is the effective contact, or surface, impedance on the $\ell^{\mbox{\tiny th}}$ electrode.  The current patterns used will be eigenvectors of the Current to Voltage map, which is a matrix approximation to the ND map, for the homogeneous saline tank. They are found numerically by simulating a homogeneous tank for static/absolute imaging, and experimentally by adaptive methods for difference/dynamic imaging. The matrix approximations to the ND maps from the conductivity distribution $\sigma$ are denoted by the $L\by L$ matrices, $\mathbf{R}_\sigma$, where $\mathbf{R}_\sigma\vec{I}^{(k)}=\vec{V}^{(k)}$, and we define $\mathbf{R}_\sigma\vec{1}=\vec{0}$.  Here the vectors $\vec{1}$, $\vec{0}$, denote vectors all of whose components are $1$, or $0$, respectively. The discrete analog of the DN map used in the CGO-based methods is given by $\mathbf{L}_\sigma:=\mathbf{R}_\sigma^{-1}$.

§.§ Calderón's Method

Following Calderón's original paper <cit.>, Calderón's method approximates the conductivity, $\sigma(x)$, from its Fourier transform.  Here we present a brief description of the method and refer the reader to <cit.> for further details. Calderón's method assumes the conductivity is a small perturbation, $\delta\sigma(x)$, from a constant background, $\sigma_b$, i.e. $\sigma(x)=\sigma_b+\delta\sigma(x)$. In this paper, we assume that the background conductivity $\sigma_b=1$.  If the background constant is not one, then the problem can be scaled and unscaled as in <cit.>.

The three steps of Calderón's method in 3-D, as described in <cit.>, are:

Step 1: Use the DN maps $\Lambda_\sigma$ and $\Lambda_1$ to approximate the Fourier transform of the small perturbation in conductivity, $\widehat{\delta\sigma}(z)$, by

\begin{equation}\label{eq:CaldFhat}
\widehat{\delta\sigma}(z)\approx\trev{\Fhat}(z):=-\frac{1}{2\pi^2|z|^2}\int_{\partial\Omega}e^{\pi i(z\cdot x)+\pi(a\cdot x)}\left(\Lambda_\sigma-\Lambda_1\right) e^{\pi i(z\cdot x)-\pi(a\cdot x)}dS(x),
\end{equation}

where $z$ and $a$ satisfy

\begin{equation}\label{eq:za_conditions}
z,a\in\R^3, |z|=|a|,\text{ and }z\cdot a=0,
\end{equation}

and $\Lambda_1$ is the DN map for a constant conductivity of $1$.

Step 2: Take the inverse Fourier transform of $\trev{\Fhat}(z)$:

\begin{equation}\label{eq:cal-FhattoGammaCart}
\delta\sigma^{\mbox{\tiny{\textbf{CAL}}}}(x) \approx\mathcal{F}^{-1}\{\trev{\Fhat}(z)\}(x)=\int_{\mathbb{R}^3}{\trev{\Fhat}(z)e^{-2\pi i(x\cdot z)}dz}.
\end{equation}

Step 3: Add the background to the perturbation to recover the approximate conductivity, $\sigma^{\mbox{\tiny{\textbf{CAL}}}}(x)$:

\begin{equation}\label{eq:sigCal_pert}
\sigma^{\mbox{\tiny{\textbf{CAL}}}}(x)=\sigma_b+\delta\sigma^{\mbox{\tiny{\textbf{CAL}}}}(x).
\end{equation}

The definition of the Fourier transform in Calderón's method is different from that used in the $\mathbf{t}^{\mbox{\tiny{\textbf{exp}}}}$ and $\tzero$ methods described below.  However, each method is consistent with its definition and is consistent with the respective literature on that method.

In this paper, we compute $\Fhat(z)$ in spherical Fourier coordinates. As such, we choose
\[z=|z|(\cos\tilde\phi\sin\tilde\theta,\sin\tilde\phi\sin\tilde\theta,\cos\tilde\theta)\text{ and } a=|z|(\cos\tilde\phi\cos\tilde\theta,\sin\tilde\phi\cos\tilde\theta,-\sin\tilde\theta),\]

for $|z|\geq0$, $0\leq\tilde\phi\leq2\pi$ and $0\leq\tilde\theta\leq\pi$ so that $z$ and $a$ satisfy (<ref>). Then, the inverse Fourier transform in (<ref>) becomes

\begin{equation}\label{eq:cal-FhattoGamma}
\delta\sigma^{\mbox{\tiny{\textbf{CAL}}}}(x) =\int_{0}^{\infty}{\int_0^{2\pi}{\int_0^\pi{|z|^2\sin\tilde\theta\trev{\Fhat}(|z|,\tilde\phi,\tilde\theta)e^{-2\pi i(x\cdot z)}d\tilde\theta d\tilde\phi d|z|}}}.
\end{equation}

Additionally, we implement the use of a mollifier, $\hat{\eta}\left(\frac{z}{y}\right)$, as introduced in <cit.> for some parameter $y\in\mathbb{R}$ to reduce Gibbs phenomenon caused by jump discontinuities in $\sigma(x)$ while recovering $\delta\sigma^{\mbox{\tiny{\textbf{CAL}}}}$

\begin{equation}\label{eq-sigcal-mollified}
\delta\sigma^{\mbox{\tiny{\textbf{CAL}}}}(x) =\int_{0}^{\infty}{\int_0^{2\pi}{\int_0^\pi{|z|^2\sin\tilde\theta\trev{\Fhat}(|z|,\tilde\phi,\tilde\theta)\hat{\eta}\left(\frac{z}{y}\right)e^{-2\pi i(x\cdot z)}d\tilde\theta d\tilde\phi d|z|}}}.
\end{equation}

We implement the same mollifier as used in <cit.>,

\begin{equation}\label{eq-Cald-moll}
\hat{\eta}\left(\frac{z}{y}\right)=e^{-\pi t|z|^2},
\end{equation}

where $y=1/\sqrt{t}$ and $t$ acts as a smoothing parameter with larger $t$ values producing smoother reconstructions with smaller jumps at points of discontinuity in $\sigma(x)$.

Since noise causes (<ref>) to blow up at large $|z|$, we use a non-uniform truncation regularization strategy. A similar regularization strategy was proved stable for the 2-D D-bar method in <cit.>, which also noted non-uniform truncation also produces reliable reconstructions.  In our case, we will first compute $\Fhat(z)$ for $|z|$ within an outer radius of $T_{z_2}$. We keep values of $\Fhat(z)$ whose real and imaginary amplitudes are below a threshold determined by the amplitudes of $\Fhat$ within a smaller radius $|z|\leq T_{z_1}$; $\Fhat$ is set to 0 everywhere else. Both radii are chosen empirically. The inner radius $T_{z_1}$ is chosen as a region in $z$-space where noise does not cause $\Fhat$ to blow up and $T_{z_2}$ is chosen to keep as much reasonable information from $\Fhat$ without introducing holes in the non-zero region of $\Fhat$. As such, our truncated $\Fhat$ is computed by

\begin{equation}\label{eq:truncFhat}
\small \Fhat_R(z)=\left\{\begin{array}{ll}
\Fhat(z),& \text{if }|z|\leq T_{z_2} \text{, and }
|\text{Re}(\Fhat(z))|\leq \underset{|z|\leq T_{z_1}}{\max}|\text{Re}(\Fhat(z))|,\,\, |\text{Im}(\Fhat(z))|\leq \underset{|z|\leq T_{z_1}}{\max}|\text{Im}(\Fhat(z))|\\
0,& \mathrm{else}.
\end{array}\right.
\end{equation}

With our truncated $\Fhat$, equation (<ref>) is truncated in the radial variable, leading to the approximation

\begin{equation}\label{eq-cal-FhattoGamma-reg}
\delta\sigma_R^{\mbox{\tiny{\textbf{CAL}}}}(x) =\int_{0}^{T_{z_2}}{\int_0^{2\pi}{\int_0^\pi{|z|^2\sin\tilde\theta\Fhat_R(|z|,\tilde\phi,\tilde\theta)e^{-\pi t|z|^2}e^{-2\pi i(x\cdot z)}d\tilde\theta d\tilde\phi d|z|}}}.
\end{equation}

For difference images shown in this paper, we only perform Steps 1 and 2 of the method and replace $\Lambda_1$ in (<ref>) with a reference DN map, $\Lambda_{\sigma_{\mbox{\tiny ref}}}$ before computing (<ref>). Thus, the flow for absolute images is

\[\left(\Lambda_\sigma,\Lambda_1\right)\overset{1}{\longrightarrow}\trev{\Fhat}(z) \overset{2}{\longrightarrow}\delta\sigma^{\mbox{\tiny{\textbf{CAL}}}}\overset{3}{\longrightarrow}\sigma^{\mbox{\tiny{\textbf{CAL}}}}\]

and the flow for difference images in this paper is

\[\left(\Lambda_\sigma,\Lambda_{\sigma_{\mbox{\tiny ref}}}\right)\overset{1}{\longrightarrow}\trev{\Fhat}(z) \overset{2}{\longrightarrow}\delta\sigma^{\mbox{\tiny{\textbf{CAL}}}}.\]

§.§.§ Numerical Implementation

In this section, we review the implementation details of Calderón's method introduced in <cit.>.

For Step 1, we compute (<ref>) for $|z|\leq T_{z_2}$ by discretizing the boundary integral as follows,

\begin{eqnarray}
\Fhat(z)&=&-\frac{1}{2\pi^2|z|^2}\int_{\partial\Omega}e^{\pi i(z\cdot x)+\pi(a\cdot x)}\left(\Lambda_\sigma-\Lambda_1\right) e^{\pi i(z\cdot x)-\pi(a\cdot x)}dS(x)\nonumber\\
&\approx& -\frac{1}{2\pi^2|z|^2}\left(\frac{|\partial\Omega|}{L}\right)(e^{\pi i(z\cdot \mathbf{x})+\pi(a\cdot \mathbf{x})})^{\mbox{\tiny \bf T}}\mathbf{Q}\left(\mathbf{L}_\sigma-\mathbf{L}_1\right)\mathbf{Q}^{\mbox{\tiny \bf T}}\left[e^{\pi i(z\cdot \mathbf{x})-\pi(a\cdot \mathbf{x})}\right],\label{eq-discrete-Fhat_bdry}
\end{eqnarray}

where $|\partial\Omega|$ is the surface area of the domain; $L$ is the number of  electrodes; $\mathbf{x}\in\R^{(L\by3)}$ denotes the vector of the Cartesian centers of the $L$ electrodes; $^{\mbox{\tiny \bf T}}$ denotes the traditional, non-conjugate, matrix transpose; $\mathbf{L}_\sigma$ and $\mathbf{L}_1$ denote the discrete matrix approximations to the DN maps $\Lambda_\sigma$ and $\Lambda_1$ respectively; and $\mathbf{Q}\in\R^{L\by \mbox{\tiny $N_{li}$}}$ denotes an orthonormal basis created using  $N_{li}$ linearly independent applied currents over $L$ electrodes as was done in <cit.>. The matrix $\mathbf{L}_1$ is based on the FEM solution of the CEM (<ref>).  Problem-specific mesh details are given below in section <ref>.

We then compute $\Fhat_R$ according to equations (<ref>) and (<ref>). For Step 2, on an equally-spaced $16\times16\times16$ rectangular grid in $x$, the conductivity difference $\delta\sigma_R^{\mbox{\tiny{\textbf{CAL}}}}(x)$ is computed via (<ref>) using a 3D Simpson's rule using $N_{|z|}=10, N_{\tilde{\theta}}=10$, and $N_{\tilde{\phi}}=30$ uniformly-spaced nodes on the $|z|, \tilde{\theta}$, and $\tilde{\phi}$ grids, respectively. As the number of nodes in the Fourier domain increases, so does computation time, but some artefact reduction can be achieved. Empirically, the artefact reduction did not seem significant enough to warrant an increased computational time beyond these parameter choices.

Difference images are computed using (<ref>) and (<ref>), replacing $\mathbf{L}_1$ with a discrete reference DN map, $\mathbf{L}_{\sigma_{\mbox{\tiny ref}}}$ in (<ref>). For the phantom tank experiments in this paper, this reference map is from data collected with a tank filled only with saline matching the experiment and no other inclusions.

The absolute reconstructions of $\sigma_R^{\mbox{\tiny{\textbf{CAL}}}}(x)$ in this paper are produced using (<ref>) replacing $\sigma_b$ with $\sigma_{\mbox{\tiny{best}}}$

\begin{equation*}
\sigma_R^{\mbox{\tiny{\textbf{CAL}}}}(x)=\sigma_{\mbox{\tiny{best}}}+\delta\sigma_R^{\mbox{\tiny{\textbf{CAL}}}}(x).
\end{equation*}

The solution is then interpolated to a $64\times64\times64$ rectangular grid.

Following <cit.>, $\sigma_{\mbox{\tiny{best}}}$ is given by

\begin{equation}\label{eq-gam-best}
\sigma_{\mbox{\tiny{best}}}=\frac{\sum_{k=1}^K\sum_{\ell=1}^LU_\ell^kU_\ell^k}{\sum_{k=1}^K\sum_{\ell=1}^LU_\ell^kV_\ell^k},
\end{equation}

where $U_\ell^k$ is the $k^{th}$ simulated voltage pattern measured on electrode $\ell$ with a homogeneous conductivity of $1$ S/m and $V_\ell^k$ is the $k^{th}$ voltage pattern measured on electrode $\ell$ for the inhomogeneous conductivity $\sigma$. The simulated voltages are the same voltages used to compute $\mathbf{L}_1$, as described in <ref>.

§.§ The $\texp$ and $\tzero$ methods

Both the $\texp$ and $\tzero$ methods are derived from the constructive proofs of <cit.> and involve special solutions called Complex Geometrical Optics (CGO) solutions <cit.>.  A brief summary is included here for the reader's convenience.  For further details see <cit.>.

Assuming that the conductivity is a constant $\sigma_c=1$ in a neighborhood of the boundary $\bndry$, the real-valued conductivity equation (<ref>) can be transformed to the Schrödinger equation

\begin{equation}\label{eq:schro}
\left(-\Delta + q(x)\right) \widetilde{u}(x)=0, \qquad x\in\R^3,
\end{equation}

via the change of variables $q(x)=\frac{\Delta\sqrt{\sigma(x)}}{\sqrt{\sigma(x)}}$ and $\widetilde{u}(x)=u(x)\sqrt{\sigma(x)}$, by extending $\sigma(x)\equiv1$ for all $x\in\R^3\setminus\bar{\Omega}$.  For $\zeta(\xi)\in \mathcal{V}_\xi$, unique CGO solutions exist to the transformed problem

\[\left(-\Delta + q(x)\right) \psi(x,\zeta)=0, \qquad x\in\R^3,\]

where $\psi(x,\zeta)\sim e^{ix\cdot\zeta}$ for large $|x|$ or $|\zeta|$, and

\begin{equation}\label{eq:Vxi}
\trev{\mathcal{V}_\xi} = \left\{\zeta\in\mathbb{C}^3\middle|\zeta^2=0, \;\; \left(\xi+\zeta\right)^2=0\right\}, \quad \text{for each $\xi\in\mathbb{R}^3$},
\end{equation}

where $\zeta^2 = \zeta\cdot\zeta$ and $\zeta$ is a purely auxiliary parameter.  The conductivity $\sigma(x)$ can then be recovered from the DN map $\Lambda_\sigma$ as follows.

For each $x\in\bndry$ and $\zeta\in \trev{\mathcal{V}_\xi}$, solve the Fredholm integral equation of the Second Kind,

\begin{equation}\label{eq:psi-bndry}
\psi(x,\zeta) = e^{ix\cdot\zeta} - \int_{\bndry} G_\zeta(x-y)\left(\Lambda_\sigma - \Lambda_1\right)\psi(y,\zeta)\;dS(y),
\end{equation}

\[G_\zeta(x)= \frac{e^{ix\cdot\zeta}}{(2\pi)^3}\int_{\R^3} \frac{e^{ix\cdot k}}{|k|^2 + 2k\cdot\zeta}dk, \quad x\in\R^3\setminus\{0\},\]

denotes the Faddeev Green's function <cit.>.  Then, evaluate the scattering data

\begin{equation}\label{eq:t-scat-tBIE}
\mathbf{t}(\xi,\zeta)=\int_{\partial\Omega} e^{-ix\cdot(\xi+\zeta)}\left(\Lambda_{\sigma}-\Lambda_1\right)\psi(x,\zeta)\;dS(x).
\end{equation}

For $|\zeta|$ large, the Schrödinger potential $q(x)$ can be recovered via the inverse Fourier transform

\begin{equation}\label{eq:t-q}
q(x)\approx \mathcal{F}^{-1}\left\{\mathbf{t}(\xi,\zeta)\right\}(x)=\frac{1}{(2\pi)^3}\int_{\mathbb{R}^3} e^{i x\cdot\xi} \mathbf{t}(\xi,\zeta)\;d\xi, \quad x\in\mathbb{R}^3.
\end{equation}

The conductivity is then recovered by solving the boundary value problem

\begin{equation}\label{eq:QtoSigma}
\left\{\begin{array}{rclcl}
(-\Delta + \trev{q}(x)) \tilde{u}(x) & =& 0 & \qquad & x\in\Omega\subset\mathbb{R}^3\\
\tilde{u}(x) & =& 1 && x\in\partial\Omega,
\end{array}
\right.
\end{equation}

for $\widetilde{u}(x)$ and evaluating $\sigma(x)=\left(\tilde{u}(x)\right)^2$.  This is the full nonlinear reconstruction method.

Replacing the CGO solutions $\psi(x,\zeta)$ by their asymptotic behavior $e^{ix\cdot\zeta}$ in the scattering data (<ref>) via

\begin{equation}\label{eq:texp-scat}
\mathbf{t}^{\mbox{\tiny{\textbf{exp}}}}(\xi,\zeta)=\int_{\partial\Omega} e^{-ix\cdot(\xi+\zeta)}\left(\Lambda_{\sigma}-\Lambda_1\right)e^{ix\cdot\zeta}\;dS(x).
\end{equation}

yields a `Born approximation' typically called the $\texp$ approximation for EIT.  Using this approximate scattering data in place of the fully nonlinear $\mathbf{t}(\xi,\zeta)$, one proceeds with the recovery of an approximate potential $\qexp(x)$ via (<ref>) and conductivity $\sigexp(x)$ via (<ref>).  The flow is:

\[\left(\Lambda_{\sigma},\Lambda_1\right)\overset{1}{\longrightarrow}\texp(\xi,\zeta) \overset{2}{\longrightarrow}\qexp(x) \overset{3}{\longrightarrow}\sigexp(x).\]

An intermediate approximation can be computed by replacing the Faddeev Green's function $G_\zeta(x)$ in the single layer potential, in (<ref>), for the traces of the CGOs, with the standard Green's function $G_0(x) = \frac{1}{4\pi|x|}$ for the Laplacian operator.  Thus, one solves

\begin{equation}\label{eq:psi0-bndry}
\psi^0(x,\zeta) = e^{ix\cdot\zeta} - \int_{\bndry} G_0(x-y)\left(\Lambda_\sigma - \Lambda_1\right)\psi^0(y,\zeta)\;dS(y),
\end{equation}

for the CGOs $\psi^0(x,\zeta)$, avoiding the exponentially growing Faddeev Green's function.  A corresponding approximation to the scattering data is then computed by using $\psi^0(x,\zeta)$ in place of $\psi(x,\zeta)$ in (<ref>) and then continuing to recover $\qzero(x)$ and $\sigma^0(x)$.  The flow is then

\[\left(\Lambda_{\sigma},\Lambda_1\right)\overset{1}{\longrightarrow}\psi^0(x,\zeta)\overset{2}{\longrightarrow}\tzero(\xi,\zeta) \overset{3}{\longrightarrow}\qzero(x) \overset{4}{\longrightarrow}\sigma^0(x).\]

We point out that the methods, as outlined above, assumed that the conductivity was a constant $\sigma_c=1$ near the boundary of the domain.  As mentioned in section <ref>, the problem can be scaled and unscaled as in <cit.> when the constant  $\sigma_c$ is not one.  In practice, we estimate the best-fit constant conductivity fit to the data $\sigbest$ as given by (<ref>). Explicitly, as in <cit.>, we scale the DN map by using $\frac{1}{\sigbest}\Lambda_\sigma$ in place of $\Lambda_\sigma$, and re-scale at the end using $\sigexp(x)=\sigbest\left(\tilde{u}^{\mbox{\tiny exp}}(x)\right)$ and $\sigma^0(x)=\sigbest\left(\tilde{u}^0(x)\right)$.

In this work we consider both the $\texp$ and $\tzero$ methods.  We note that this is the first time that $\tzero$ has been implemented on non-continuum DN data and the first time that either $\texp$ or $\tzero$ have been demonstrated on experimental 3D EIT data, for absolute or time-difference EIT imaging.

§.§.§ Numerical Implementation

Here we provide the numerical details pertinent to the implementation of the $\texp$ and $\tzero$ algorithms outlined above.  As with the Calderón method above, the main idea is to expand functions in the same orthonormal basis $\mathbf{Q}$ as described in section <ref>.

Following <cit.>, for each electrode center $x_\ell$, expand $\psi(x,\zeta)$ and $e^{ix\cdot\zeta}$ as

\begin{equation}\label{eq:basis_expansions}
\psi(x_\ell,\zeta)\approx \sum_{j=1}^{N_{li}} b_j(\zeta)\mathbf{Q}^j_\ell, \qquad% \ell=1, \ldots L \label{eq:psi_expansion}\\
e^{ix_\ell\cdot\zeta}\approx \sum_{j=1}^{N_{li}} c_j(\zeta)\mathbf{Q}^j_\ell, \qquad \ell=1, \ldots L, %\label{eq:exp_expansion}
\end{equation}

where $\mathbf{Q}^j_\ell$ denotes the $(\ell,j)$ entry of $\mathbf{Q}$.  Then, the boundary integral equation (<ref>) can be approximated as follows

\begin{eqnarray}
\psi^0(x_\ell,\zeta)
&=&= e^{ix_\ell\cdot\zeta} - \int_{\bndry} G_0(x_\ell-y)\left(\Lambda_\sigma - \Lambda_1\right)\psi^0(y,\zeta)\;dS(y)\nonumber\\
&\approx& e^{ix_\ell\cdot\zeta} - \sum_{\ell'=1}^L\int_{E_{\ell'}} G_0(x_\ell-y)\left(\Lambda_\sigma - \Lambda_1\right)\psi^0(y,\zeta)\;dS(y)\nonumber\\
&\approx& e^{ix_\ell\cdot\zeta} - \left[\sum_{\ell'=1}^L\int_{E_{\ell'}} G_0(x_\ell-y)dS(y)\right]\left[\left(\mathbf{L}_\sigma - \mathbf{L}_1\right)\psi^0(y_{\ell'},\zeta)\right],\label{eq:psi_bie_expan1}
\end{eqnarray}

where $E_{\ell'}$ denotes the $\ell'^{\mbox{\tiny th}}$ extended electrode where $\bigcup_{\ell'=1}^L E_{\ell'}=\bndry$ and the $E_{\ell'}$ are mutually disjoint <cit.>.  Note the true electrodes need not cover the surface $\bndry$, only the extended (mathematical) electrodes that we will use to discretize the integral.  Then, using the expansions from (<ref>) in (<ref>), we have

\begin{eqnarray}
\sum_{j=1}^{N_{li}} b_j(\zeta)\mathbf{Q}^j_\ell
&\approx& \sum_{j=1}^{N_{li}} c_j(\zeta)\mathbf{Q}^j_\ell - \left[\sum_{\ell'=1}^L\int_{E_{\ell'}} G_0(x_\ell-y)dS(y)\right]\left[\left(\mathbf{L}_\sigma - \mathbf{L}_1\right)\sum_{j=1}^{N_{li}} b_j(\zeta)\mathbf{Q}^j_{\ell'}\right] \nonumber\\
&=& \sum_{j=1}^{N_{li}} c_j(\zeta)\mathbf{Q}^j_\ell - \left[\sum_{\ell'=1}^L\int_{E_{\ell'}} G_0(x_\ell-y)dS(y)\right]\left[\sum_{j=1}^{N_{li}} b_j(\zeta)f_j\left(y_{\ell'}\right)\right], \label{eq:eq:psi_bie_expan2}
\end{eqnarray}

where $f_j\left(y_{\ell'}\right)$ represents the action of $\left(\mathbf{L}_\sigma - \mathbf{L}_1\right)$ on $\mathbf{Q}^j$ evaluated at $y_{\ell'}$, which can be computed as the $(\ell',j)$ entry of $\mathbf{Q}\left(\mathbf{L}_\sigma - \mathbf{L}_1\right)$.  Define

\begin{equation}\label{eq:G0-mat}
\widetilde{\mathbf{G}}_0(\ell,\ell')=\begin{cases} \frac{1}{4\pi\left|x_\ell - y_{\ell'}\right|} & \ell\neq\ell'\\
0 & \ell=\ell',
\end{cases}
\end{equation}

where we have removed the singularities at $x_\ell=y_\ell$.  Assuming $|E_\ell|=\frac{|\bndry|}{L}$ for each $\ell=1,\ldots,L$, and using (<ref>) in (<ref>) we find

\begin{eqnarray}
\sum_{j=1}^{N_{li}} b_j(\zeta)\mathbf{Q}^j_\ell
&\approx& \sum_{j=1}^{N_{li}} c_j(\zeta)\mathbf{Q}^j_\ell - \frac{|\bndry|}{L}\sum_{j=1}^{N_{li}} b_j(\zeta) \sum_{\ell'=1}^L \widetilde{\mathbf{G}}_0(\ell,\ell')f_j\left(y_{\ell'}\right) \nonumber\\
&\approx& \sum_{j=1}^{N_{li}} c_j(\zeta)\mathbf{Q}^j_\ell - \frac{|\bndry|}{L}\sum_{j=1}^{N_{li}} b_j(\zeta) \left[\widetilde{\mathbf{G}}_0\mathbf{Q}\left(\mathbf{L}_\sigma - \mathbf{L}_1\right)\right](\ell,j) \nonumber
\end{eqnarray}

or in matrix form,

\[\mathbf{Q}\vec{b} = \mathbf{Q}\vec{c} - \frac{|\bndry|}{L}\widetilde{\mathbf{G}}_0\mathbf{Q}\left(\mathbf{L}_\sigma - \mathbf{L}_1\right)\vec{b}.\]

The solution to this equation can be found by solving the following system for the unknowns $\vec{b}$,

\begin{equation}\label{eq:psi0-bie-matrixform}
\left(I+A\right)\vec{b} = \vec{c},
\end{equation}

where $A=\frac{|\bndry|}{L}\mathbf{Q}^{\mbox{\tiny \bf T}}\widetilde{\mathbf{G}}_0\mathbf{Q}\left(\mathbf{L}_\sigma - \mathbf{L}_1\right)$, and $I$ is the identity matrix.  If the extended electrodes $E_\ell$ are not uniform in size, then one could compute as weighted sum replacing the uniform weight $|E_\ell|=\frac{|\bndry|}{L}$ as appropriate.

Next, following <cit.>, the scattering data $\tzero(\xi,\zeta(\xi))$ can be computed for all $|\xi|$ less than a chosen truncation radius $T_{\xi}$ via

\begin{equation}\label{eq:tzero-matrix}
\tzero(\xi,\zeta)\approx
\begin{cases}
\frac{|\bndry|}{L} \left[e^{-i\mathbf{x}\cdot(\xi + \zeta)}\right]^{\mbox{\tiny \bf T}}\mathbf{Q}\left(\mathbf{L}_\sigma - \mathbf{L}_1\right)\vec{b} & |\xi|<T_{\xi}\\
0 & else,
\end{cases}
\end{equation}

since $\vec{b}=\mathbf{Q}^{\mbox{\tiny \bf T}}\psi^0(\mathbf{x},\zeta)$, where $\mathbf{x}\in\R^{(L\by3)}$ is the same as in Section <ref>, i.e. a vector storing the centers of the electrodes.

Next, the approximate potential $\qzero$ is recovered by computing the inverse Fourier transform of the truncated scattering data $\tzero$ using a Simpson's rule in 3D,

\[\qzero(x)=\frac{1}{(2\pi)^3}\int_{[-T_\xi,T_\xi]^3} e^{i x\cdot\xi} \tzero(\xi,\zeta(\xi))\;d\xi.\]

Alternatively, one could use an IFFT to achieve additional speedup, taking care with quadrature points and the particular form of the kernel.  Following <cit.>, the conductivity $\sigma^0(x)$ was recovered by first solving the boundary value problem (<ref>), using the PDE toolbox in Matlab using a mesh with approximately $21,000$ 3D elements, then computing $\sigma^0(x)=\sigma{\mbox{\tiny{best}}}\left(\tilde{u}(x)\right)^2$.  The 3D visualizations were obtained by interpolation to a $64\times64\times64$ rectangular grid using Matlab's scatteredInterpolant function.
The $\zeta$ values were computed following the minimal-zeta approach outlined in <cit.>.  We remark that while the non-existence of exceptional points for the solution of the boundary integral equation (<ref>) is proven for large $|\zeta|$ <cit.>, as well as small $|\zeta|$  <cit.>; it is still an open question for the intermediate values required here to perform the computation on a computer.

The $\texp$ reconstructions were obtained in an analogous fashion to those of $\tzero$, this time bypassing the boundary integral equation (<ref>) and directly computing

\[
\texp(\xi,\zeta)\approx
\begin{cases}
\frac{|\bndry|}{L} \left[e^{-i\mathbf{x}\cdot(\xi + \zeta)}\right]^{\mbox{\tiny \bf T}}\mathbf{Q}\left(\mathbf{L}_\sigma - \mathbf{L}_1\right)\mathbf{Q}^{\mbox{\tiny \bf T}}\left[e^{i\mathbf{x}\cdot\zeta}\right] & |\xi|<T_{\xi}\\
0 & else,
\end{cases}\]

where we have replaced the vector of coefficients $\vec{b}$ with the expansion of the asymptotic behavior of $e^{ix\cdot\zeta}$ given by $\mathbf{Q}^{\mbox{\tiny \bf T}}\left[e^{i\mathbf{x}\cdot\zeta}\right]$.

Difference imaging can be performed with the $\tzero$ and $\texp$ methods by replacing the matrix $\mathbf{L}_1$ with $\mathbf{L}_{\sigma_{\mbox{\tiny ref}}}$ and computing $\sigma_{\mbox{\tiny diff}(x)} = \sigma(x) - \sigbest$ in the final step.

§.§ Reference methods

Total Variation regularized non-linear least squares reconstructions will serve as the reference reconstructions for the CGO-based absolute imaging cases considered here.  A classic linear difference imaging scheme, also reviewed below, will serve as the reference for the CGO-based difference images.

§.§.§ Absolute imaging with TV regularization

A widely used numerical approach for absolute EIT is the total variation (TV) regularized non-linear least squares minimization

\begin{equation}
\hat\sigma = {\rm arg} \min_{\sigma>0}\{\Vert L_e \left( V-U(\sigma,z_\ast) \right)\Vert^2
+ \alpha {\rm TV}_\beta (\sigma)\},
\label{lssol}
\end{equation}

where $U(\sigma,z)$ is the finite dimensional forward map, $z_\ast \in \mathbb{R}^L$ are the fixed electrode contact impedances obtained from an initialization step of the minimization, $L_e$ is the Cholesky factor of the noise precision matrix of $e$ so that $L_e^{\rm T}L_e = \Gamma_e^{-1}$, scalar valued $\alpha$ is the regularization parameter and ${\rm TV}_\beta (\sigma)$ is the (smooth) TV regularization functional <cit.>

\begin{equation}\label{stvpsi}
{\rm TV}_\beta (\sigma) = \int_\Omega \sqrt{\| \nabla \sigma \|^2 +
\beta }\ {\rm d} r,
\end{equation}

where $\beta$ is the (fixed) smoothing parameter. The forward model $U(\sigma,z)$ in (<ref>) is based on the finite element (FEM) discretization of the complete electrode model <cit.>. For details of the FEM model, see
<cit.>. In the FEM model, the electric conductivity is approximated as a linear combination of the piecewise linear nodal basis functions in a uniform tetrahedral mesh of $N$ nodes, leading to vector of unknowns $\sigma \in \mathbb{R}^N$, and the electric potential is approximated similarly in a significantly more dense tetrahedral mesh with refinements near the electrodes.
The non-linear optimization in (<ref>) is solved by a lagged Gauss-Newton method equipped with
a line search algorithm. The line search is implemented using bounded minimization such that the non-negativity $\sigma > 0$ is enforced. For more details of the method, see <cit.>.

The fixed contact impedances $z_\ast$ and initial (constant) conductivity estimate for (<ref>) are obtained from the solution of the non-linear least squares problem

\begin{equation}
(\sigma_0, z_\ast) = {\rm arg} \min_{\sigma_c, z >0} \{\Vert L_e \left( V-U(\sigma_c,z) \right)\Vert^2\},
\label{lsinit}
\end{equation}

where the scalar $\sigma_c \in \mathbb{R}$ is the coefficient of a spatially constant conductivity image $\sigma_c \mathbf{1}$ and $z\in \mathbb{R}^L$. The non-linear least squares problem (<ref>) is solved by a Gauss Newton optimization.

§.§.§ Linear difference imaging

In linear difference imaging, see e.g. <cit.>, the objective is to reconstruct the change in conductivity between two measurements (V1; V2). The reference linear difference imaging approach of this paper uses linearized approximations of the observation models

\begin{eqnarray}
V_1 \approx U(\sigma_0,z_\ast) + J_\sigma (\sigma_1-\sigma_0) + e_1\\
V_2 \approx U(\sigma_0,z_\ast) + J_\sigma (\sigma_2-\sigma_0) + e_2,
\end{eqnarray}

where the Jacobian matrix $J_\sigma$ of $U(\sigma,z_\ast)$ is evaluated at $\sigma = \sigma_0$. With these linearized models, the difference in the measurements becomes

\begin{eqnarray}
\nonumber &\delta V &= V_2 - V_1\\
\nonumber &&= \left(U(\sigma_0) + J_\sigma (\sigma_2-\sigma_0) + e_2\right) - \left(U(\sigma_0) + J_\sigma (\sigma_1-\sigma_0) + e_1\right) \\
\label{equ:dV}
&&= J_\sigma \delta\sigma + \delta e,
\end{eqnarray}

where $\delta \sigma = \sigma_2 - \sigma_1$ and $\delta e = e_2 - e_1$.
Now, the inverse problem is to reconstruct $\delta \sigma$ based on the difference data $\delta V$ and the model (<ref>). A widely used formulation for the linear problem is to use the (generalized) Tikhonov regularization with a smoothness promoting regularization functional

\begin{equation}
\label{equ:minLD}
\hat{\delta\sigma} = \arg \min_{\delta \sigma} \left\{ || L_{\delta e}(\delta V - J_\sigma \delta\sigma) ||^2 + % p_{LD}(\delta\sigma)
||L_p \delta \sigma||^2
\right\},
\end{equation}

where $L_{\delta e}$ is the Cholesky factor of the noise precision matrix of $\delta e$ so that $L_{\delta e}^{\rm T}L_{\delta e} = \Gamma_{\delta e}^{-1} = \left(\Gamma_{e_1} + \Gamma_{e_2}\right)^{-1}$.
In this paper, the regularization matrix $L_p$ is constructed by utilization of a distance based covariance function. More specifically, we set $L_p^{\mathrm{T}}L_p = \Gamma_p^{-1}$, where the (prior covariance) matrix $\Gamma_p$ is constructed using the distance based correlation function

\begin{align}
\label{equ:GammaSmooth}
\Gamma_{p}(i,j) = \mathrm{std(}\sigma)^2\exp\left(-\frac{\|x_i - x_j \|^2}{2a^2}\right), \ \ \ i, j = 1,\ldots, N,
\end{align}

where the parameter $a$ controls the correlation length and can be solved by setting the distance $\|x_i - x_j \|$ to a selected value $d$ (e.g. half the radius of the target) and setting $\Gamma_{p}(i,j)$ to the desired covariance for that distance (e.g. 1% of variance).

§.§ Experimental Setup

ACT5 <cit.> is a 32 electrode parallel EIT instrument that uses 32 current sources to apply patterns of current to the target and 32 voltmeters to measure the resulting voltages. The current sources in ACT5 adjust the delivered current to compensate for current lost through shunt impedance, including that introduced by the capacitance of cables that connect the instrument to the electrodes, enabling the desired current patterns to be applied with high precision <cit.>.  ACT5 can apply sinusoidal currents at frequencies in the range of 11 kHz to 1 MHz with a peak amplitude of 0.25 mA. The voltmeters measure voltages  up to 0.5 V peak with a maximum signal-to-noise ratio (SNR) of 96 dB. Because the current source compensates for shunt capacitance of the cables, the system operates with grounded-shield cables <cit.>.

A test tank was built from 3/8-inch-thick Plexiglas with interior dimensions of 17.0cm x 25.5cm x 17.0cm.  Thirty-two electrodes were fabricated of 16-gauge 316 stainless steel, each 80.0 mm square.  The inside surfaces of the tank were milled out in the shape and depth of the electrodes so that the interior surfaces were flat.  The nominal gaps between electrodes were 5.0 mm. The 32 electrodes were placed four on each end and six on each side.  Bolts through the center of each electrode passed through the Plexiglas, allowing connections to the cables from the ACT5 electronics.  Five sides of the tank were glued together, and provision was made to secure the top with threaded rods at all four corners.  Access holes of 0.6 mm at various sites in the top to allow threads to suspend targets, and 25-Ga hypodermic needles to complete the filling of the tank after the lid was in place.  In use, a mixture of saline at room temperature at the desired background conductivity is made and measured using a Oakton CON 6+ handheld conductivity meter.

Targets were made the day before experiments by adding NaCl and a few drops of food dye to distilled water until the desired conductivity was obtained, as measured by an immersible conductivity meter.  The solution was then heated slowly with stirring as 4% agar (Fisher Scientific) was added until the temperature reached $85^{\circ}$C$-90^{\circ}$C.  The mixture was poured into 5.27 cm inside-diameter spherical molds and allowed to harden overnight.  Test cells were also filled at this time so the final conductivities could be verified using the ACT5 instrument.

Data were first collected of the tank filled with saline of the desired conductivity (24 mS/m) with no targets present.  Then, targets were added to the tank supported on fine toothpicks from below.  Photos were captured of the tank, with each target in place, to verify its/their position.  The conductivity values of the spheres used were approximately 290 mS/m. Averaged data, over 100 frames, were used in the subsequent conductivity reconstructions.  The experimental data will be made publicly available on github.[See <https://github.com/sarahjhamilton/open3D_EIT_data>]

Note that even though we are testing on a rectangular prism tank, the mathematical reconstruction algorithms are not limited to this geometry.

Experimental setups. Left: Top view of the one target experiment. Middle: Top view of two target experiment. Right: side view of height of targets above floor of tank. Note that the targets were measured at 290 mS/m in test cells.

§.§ Robustness Tests

In addition to testing the reconstruction methods with the correct domain modeling, we consider a moderate, as well as strong, mis-modeling of the domain.  The incorrect domain modeling comes into play in the following places for the CGO-based methods.  First in the simulation of the $\mathbf{L}_1$ DN data that requires solving the forward EIT problem with a known conductivity $\sigma\equiv 1$ S/m and the applied current patterns.  Next, the incorrect domain modeling presents as incorrect information about the centers of the electrodes $\mathbf{x}$ and surface area for the domain $|\bndry|$; the subsequent reconstructions of $\sigcal$, $\sigmaexp$, and $\sigma^0$ are recovered on the incorrect domain. For the TV regularized method, the incorrect domain modeling is present in the FEM based forward map $U(\sigma,z)$ used in the minimization problem (<ref>) and computational mesh. For the linear difference imaging method, the incorrect domain modeling is present in the Jacobian matrix of the forward map $J_{\delta \sigma}$ in the minimization problem (<ref>) and computational mesh.

Accurate domain modeling can be challenging in practice.  To address this we explored a domain with a moderately incorrect domain shape, $18$cm x $27$cm x $19$cm, as well as a domain with a significantly incorrect domain shape, $20$cm x $35$cm x $25$cm[Due to the smaller ratio of longest side to shortest side, the $16\times16\times16$ $X$-grid used for Calderón's method as mentioned in <ref> does not capture the whole domain, so a $32\times32\times32$ $x$-grid is used in the significantly incorrect domain case].  Electrodes were uniformly distributed as in the original box design, but clearly had different locations and spacing than the truth due to the new box dimensions.  Figure <ref> (center and right) depict the larger box domains. Each reconstruction method was then tested assuming that the actual experimental voltage and current measurements were coming from the boxes shown in Figure <ref>.

(0,115)Correct Model (17x25.5x17)
(180,115)Box 18x27x19
(320,115)Box 20x35x25

3D renderings of simulated boxes used in robustness testing.

Time-difference reconstructions were performed using basal measurements with only 24 mS/m saline in the experimental tank.  This data was then used in the reconstruction methods discussed in <ref> and <ref>.

§.§.§ Forward modeling and Regularization Parameters

For the CGO-based methods, the matrix approximation $\mathbf{L}_1$ to the DN map $\Lambda_1$ was formed using simulated voltage data, using EIDORS <cit.>.  The EIDORS data is produced by solving the forward conductivity problem with the finite element approximation of the Complete Electrode Model (<ref>) with $\sigma\equiv 1$ S/m.  In each domain modeling case, see Figure <ref>, the forward problem was solved on a box of the corresponding dimensions, with 8cm by 8cm electrodes, using approximately $250,000$ nodes and 1.2 million elements with high refinement used near the electrodes, and the default contact impedance in EIDORS. For the reference methods, the FEM based forward problem was solved using the same meshes for the electric potential as with the CGO-based methods and the conductivities were approximated in sparser, uniform meshes with $8,918$; $9,410$; $14,678$ nodes and $43,690$; $45,106$; $72,438$ elements, respectively, for the different domain modeling cases.

The regularization parameters for each of the reconstruction methods were chosen as follows. The regularization parameters $T_{z_1}$ and $T_{z_2}$ for Calderón's method were empirically chosen to be 1.4 and 1.7, respectively, and the mollifying parameter was chosen to be $t=0.1$ for all reconstructions. These values were chosen empirically to provide the best visual reconstructions.
The truncation radius $T_\xi$ for the nonlinear scattering data used in the $\texp$ and $\tzero$ methods was chosen empirically in the range [10.5, 12] for each case shown here.  The localization of the target(s) did not appear to change significantly with $T_\xi$, however the contrast does appear sensitive to this parameter choice.  Further discussion is given below in Section <ref>. For the total variation regularized method, the regularization parameter $\alpha$ and the smoothing parameter $\beta$ were selected by computing a series of reconstructions with varying $\alpha$ and $\beta$ values, in the range [1e-6, 100] for $\alpha$ and [1e-4, 0.1] for $\beta$, and by choosing the parameters that gave the best visual quality of reconstructions in the correct domain model case. The chosen values $\alpha = 0.01$ and $\beta = 0.001$ were used in all the test cases with different domain models.  For the linear difference imaging method, the visually best reconstructions were obtained using a standard deviation of conductivity $\text{std}(\sigma) = 16\sigma_0$ and parameter $a$ calculated using a covariance of 1 $\%$ of the variance at a correlation distance $d = l_{\text{box}}/8$, where $l_{\text{box}}$ is the length of the longest side of the box that is used as the domain model.

§.§ Evaluation Metrics

As the main focus of this work is experimental data reconstruction, comparing to a known `truth' with zero error is not possible. The conductivity values for the targets were measured at approximately 290 mS/m and the targets were placed to be roughly centered in a sub-cube of the overall box.  These data allow us to approximate the localization error (LE) of our reconstructions and maximum conductivity in each reconstructed target.

LE, as in <cit.>, measures the  distance between the centroids of the reconstructed targets and those of the true targets by

\begin{equation}\label{eq:LE}
\text{LE}=\sqrt{(x_1^{\mbox{\tiny recon}}-x_1^{\mbox{\tiny truth}})^2+(x_2^{\mbox{\tiny recon}}-x_2^{\mbox{\tiny truth}})^2+(x_3^{\mbox{\tiny recon}}-x_3^{\mbox{\tiny truth}})^2}.
\end{equation}

A localization error of $0$ is ideal, and signifies the targets are reconstructed in the correct location. In our experiments, the centers of the spherical targets were placed at the centers of the nearest electrode in each direction and in (<ref>),  $(x_1^{\mbox{\tiny truth}},x_2^{\mbox{\tiny truth}},x_3^{\mbox{\tiny truth}})$ is our best estimate of the true location of the target center, acknowledging possible errors within millimeters of the truth. Following <cit.>, to obtain the centers of the reconstructed targets, we use a thresholded segmentation of the reconstructions to identify targets with MATLAB's . The thresholds for segmentation were used to empirically identify regions with more than a certain percentage of the maximum reconstructed conductivity. We note that the volume of the segmented targets changes with the choice of threshold, but as we increase the percentage (i.e. threshold) the locations of the centroids stabilize (up to millimeters), from which we computed the LE.

As a major focus of this work is to study how the algorithms tolerate modeling errors, we report a scaled version of the LE,

\begin{equation}\label{eq:LE_scaled}
\text{Scaled LE}=\frac{\sqrt{(x_1^{\mbox{\tiny recon}}-x_1^{\mbox{\tiny truth}})^2+(x_2^{\mbox{\tiny recon}}-x_2^{\mbox{\tiny truth}})^2+(x_3^{\mbox{\tiny recon}}-x_3^{\mbox{\tiny truth}})^2}}{25.5},
\end{equation}

where we have divided the LE from (<ref>) by 25.5cm, the length of the longest side of the true tank.  This scaling value for the computing the LE metric is kept fixed across incorrect modeling scenarios.

§ RESULTS

We first present absolute EIT reconstructions with correct domain modeling, for the one and two targets, and for the four reconstruction methods (Figure <ref>).  Cross-sections in the center of the bottom row of electrodes of the box are shown on the left; this slice cuts through the center(s) of the physical target(s).  Isosurfaces, shown in the middle column, are extracted from the 3D reconstructions to indicate the locations, size, and number of targets.  The isosurfaces were extracted using MATLAB's  with threshold segmentation in the range of 60% to 85% of the maximum recovered conductivity in the tank.  A uniform threshold across all examples and reconstruction methods was not used as the blurring varied across methods.  Note that since the isosurfaces are dependent on this threshold, they can omit information from the full reconstruction. Therefore, the third column shows a full 3D rendering of conductivity in the domain. This 3D rendering is produced in MATLAB by displaying 100 equally-spaced $x_1$-$x_3$ slices set to $0.1$ transparency of each slice.

One target absolute imaging evaluation metrics across all domain modeling scenarios.

Calderón     $\texp$      $\tzero$        TV
3[0]*Scaled LE           correct domain       0.042         0.094       0.098       0.043

box 18x27x19         0.075         0.131       0.119       0.021
box 20x35x25         0.241         0.250      0.231       n/a
3[0]*$\begin{array}{c}\text{Max. Target}\\ \text{Conductivity (mS/m)}\end{array}$    correct domain       63.77        300.00     308.67     44.26

box 18x27x19         76.94        303.32     284.67     70.14
box 20x35x25         85.51        305.63     260.88     n/a

Two target absolute imaging evaluation metrics across all domain modeling scenarios. Note: Target 1 is the same target as in the one-target case.

1|r|                         Calderón    $\mathbf{t}^{\mbox{\tiny{\textbf{exp}}}}$      1c|$\mathbf{t}^{\mbox{\tiny{\textbf{0}}}}$     1c|TV

6[0]*Scaled LE             3*Target 1     correct domain       0.056         0.098       0.105       0.034

box 18x27x19         0.103         0.135       0.133       0.021
box 20x35x25         0.268         0.267       0.257       n/a
2-72.5                                                 3*Target 2     correct domain       0.062        0.195       0.180       0.031

box 18x27x19         0.120         0.191       0.214       0.057
box 20x35x25         0.277         0.311       0.295       n/a
7[0]*$\begin{array}{c}\text{Max. Target}\\\text{Conductivity (mS/m)}\end{array}$     3*Target 1     correct domain       55.61        310.95     315.50     55.57

box 18x27x19         73.53        216.29     322.39     68.29
box 20x35x25         80.00        334.58     264.72     n/a
2-72.5                                             3*Target 2     correct domain       60.47        234.85     226.35     53.94

box 18x27x19         67.94        153.82     203.84     64.22
box 20x35x25         77.81        285.06     244.80     n/a

Two target difference imaging evaluation metrics across all domain modeling scenarios. Note: Target 1 is the same target as in the one-target case.

1|r|                         Calderón    $\mathbf{t}^{\mbox{\tiny{\textbf{exp}}}}$      1c|$\mathbf{t}^{\mbox{\tiny{\textbf{0}}}}$     1c|Linear

6[0]*Scaled LE             3*Target 1     correct domain       0.034         0.095       0.103       0.049

box 18x27x19         0.022         0.116       0.124       0.065
box 20x35x25         0.113         0.242       0.243       n/a
2-72.5                                                 3*Target 2     correct domain       0.053         0.100       0.094       0.050

box 18x27x19         0.043         0.116       0.115       0.029
box 20x35x25         0.112         0.230      0.233       n/a
6[0]*$\begin{array}{c}\text{Max. Target}\\\text{Conductivity (mS/m)}\end{array}$     3*Target 1     correct domain       57.18        277.83     295.63     484.83

box 18x27x19         61.72        272.80     282.05     227.95
box 20x35x25         56.30        253.46     243.88     227.95
2-72.5                                             3*Target 2     correct domain       53.84        266.91     284.25     490.88

box 18x27x19         59.91        294.17     306.30     208.05
box 20x35x25         55.07        293.61     282.38     n/a

Figure <ref> shows the absolute EIT reconstructions for the moderately mismodelled domain whereas Figure <ref> shows the absolute EIT reconstructions for the significantly mismodelled domain.  Difference images are shown in Figure <ref>.  Metrics of maximum conductivity per target and scaled localization error are shown in Tables <ref> and <ref> for the absolute images and Table <ref> for the difference images. Bolded entries correspond to the best metric value in each row.

(68,400)One Target
(298,400)Two Targets

Absolute image reconstructions comparing the CGO methods to the regularized method with correct domain modeling.  Slices, isosurfaces, and 3D renderings of the conductivity are shown.  Note the truth targets had a measured conductivity of approx 290 mS/m.

(68,400)One Target
(298,400)Two Targets

Absolute image reconstructions comparing the CGO methods to the regularized method with moderately incorrect domain modeling, using a box of size 18cm x 27cm x 19cm.  Note the truth targets had a measured conductivity of approx 290 mS/m.

(68,400)One Target
(298,400)Two Targets

Absolute image reconstructions comparing the CGO methods to the regularized method with largely incorrect domain modeling, using a box of size 20cm x 35cm x 25cm.  Note the truth targets had a measured conductivity of approx 290 mS/m.

(20,395)Correct modeling
(170,395)18x27x19 modeling
(320,395)20x35x25 modeling

Difference image reconstructions comparing the CGO methods to a typical linear method.  Slices and 3D renderings of the conductivity are shown for the correct domain modeling, and increasing levels of error in domain modeling.  Note the truth targets had a measured conductivity difference from the background of approx 266 mS/m.

§ DISCUSSION

Beginning with the correct domain modeling scenario, each method, Calderón, $\texp$, $\tzero$, and TV produced reconstructions that clearly show the target(s) with good localization. The TV images are sharpest, as expected.  The $\texp$ and $\tzero$ methods achieved the best contrast and best approximated the estimated experimental conductivity.  The Calderón and TV methods achieved lower contrast but better scaled localization error.  Notably, as shown in Figure <ref>, in terms of artefacts, the absolute CGO reconstructions are as clean as their corresponding difference images.  The CGO methods required simulated voltages data for the same experimental setup but with a conductivity $\sigma\equiv1$ S/m.  Even though this data was not tuned to the EIT machine (noise and contact impedances) the methods produced high quality absolute reconstructions.

Moving into incorrect domain modeling, we see all methods handled the moderate domain mismodeling quite well and the CGO methods handled the severe domain mismodeling as well.  No significant boundary artefacts are seen in the CGO reconstructions for the 18cm x 27cm x 19cm and only moderate boundary artefacts in the TV reconstructions.  In the more severe mis-modeling case, assuming the data is coming from a much larger box (20cm x 35cm x 25cm) instead of the true domain (17cm x 25.5cm x 17cm), the $\texp$ and $\tzero$ methods still produce good reconstructions showing the correct number of targets and correct region of the tank.  The localization is worst in the $x_3$ direction showing the targets slightly elongated.  While artefacts in the Calderón 20x35x25 reconstructions are present, the method still does detect the true objects as the most conductive.  The TV algorithm failed to identify any of the targets for both cases with the 20x35x25 box.
The artefacts in the TV reconstructions in these cases with incorrect domain model could be expected as
the regularized non-linear least squares minimization-based absolute imaging EIT algorithms are known to be highly sensitive to modeling errors such as errors in modelling of the domain shape, see e.g. <cit.>.

We can see the effect of the underlying assumption for Calderón's method that the conductivity is a small perturbation from a constant.  In the physical experiments, the conductor(s) had roughly twelve times the conductivity of the background saline. This is seen in the underestimation of maximum conductivity from Calderón's method.  However, unsurprisingly, this method does a good job of localizing the targets, in all cases except for the absolute image from the significant domain modeling error.  The artefacts seen in the Calderón reconstructions are consistent with what has been observed for absolute imaging with Calderón in 2D, such as Gibbs phenomenon for all absolute reconstructions and higher reconstructed conductivity at electrode locations when the domains were incorrectly modelled.  Across all modeling scenarios, the regularized TV method struggled with the high contrast targets, giving maximum target conductivities well below the measured values. The same effect was seen in simulated scenarios where, however, maximum conductivities of lower contrast targets were reliably recovered. This is at least partially explained by saturation of contrast distinquishability of the measurements, i.e. the effect of increasing the conductivity of the target inclusion on the measurements gradually diminishes and until certain level the TV regularized method can no longer discern between high and even higher conductivities.

We note that the electrodes used in this experiment were very large and the structure of the domain, a box with corners, may exacerbate some of the modeling and/or hardware challenges.  Nevertheless, the study provides informative results on the feasibility of absolute EIT reconstruction in 3D.

The difference images from Calderón are able to handle the stronger mismodeling of the domain, as are the $\texp$ and $\tzero$ methods.  The strong mismodeling proved too severe for the linear difference imaging reference method, which did not manage to identify the targets, see Figure <ref>.

While the $\texp$ and $\tzero$ CGO methods did reliably recover the contrast and approximate location of the targets across examples studied here, they do appear more sensitive than their 2D D-bar based counterparts in regards to the regularization parameter, $T_\xi$, used in the truncation of the nonlinear scattering data.  Figure <ref> displays the effect that $T_\xi$ plays on the scaled localization error and maximum recovered conductivity value for the single target, correct domain modeling case. Recall that a secondary nonuniform truncation is also enforced where scattering data with magnitudes exceeding 20 for the real or imaginary parts are set to zero.  Adjusting that value will also have an effect on the reconstruction.  The value of 20 was chosen in this work for its overall reliability across examples. The contrast appears more sensitive than the localization error.  As in <cit.>, the minimum $\zeta$ parameterization of $\xi$ was used, as the scattering data is a function of both $\zeta\in\C^3$ and $\xi\in\R^3$ in 3D instead of just $k\in\C$ as in 2D D-bar based methods.  Alternative parameterizations and a more detailed study of the effect of the regularization parameters, while interesting, are outside the scope of this work.

(15,135)Scaled Localization Error
(200,135)Max Target Cond. (mS/m)

Comparison of the effect of the truncation value $T_\xi$ of the scattering radius in the $\texp$ and $\tzero$ CGO methods for Scaled Localization Error (left) and Maximum value of the recovered target.  Max conductivity values (right) for $T_\xi=11,\;11.5$, off the plot, spiked into the 1500-5000 mS/m range.

In terms of speed, when running reconstructions on a MacBook Pro with a 2.3 GHz Dual-Core Intel® Core i5 processor, Calderón reconstructions on a $16\times16\times16$ $x$-grid which are interpolated to a $64\times64\times64$ $x$-grid take 1 to 2 seconds without optimizing for parallelization. This increases to 6-8 seconds per reconstruction when the initial $x$-grid is $32\times32\times32$. When running on a PC with a AMD EPYC 7702P 64-Core Processor 2.00 GHz, the reconstruction times are 0.6-0.7 seconds and 4 seconds, respectively, again without optimizing for parallelization.  On an 2015 iMac with a 4 GHz Quad-Core Intel® Core i7 processor, the $\texp$ and $\tzero$ methods require 3-4 seconds/recon using a $21\by21\by21$ $x$-grid for the potential $q(x)$ or 6-8 seconds/recon when using a $41\by41\by41$ $x$-grid for $q(x)$.  The timings are non-optimized with the highest computational cost coming from computing the inverse Fourier transform and solving the boundary value problem using FEM.  The regularized TV, and linear difference, reconstructions averaged 2-3 hrs, and 3 minutes, respectively, when computed on a server with 256GB of RAM and two 10 core Intel® Xeon® CPU E5-2630 v4 @2.20GHz processors. We remark that the rather long computation times of the TV regularized non-linear least squares approach are caused by the 3D problem with a computationally rather challenging geometry as the sufficient accuracy of the CEM  forward model (<ref>) necessitates significant mesh refinement near the electrodes, leading to the large number of degrees of freedom (approx 250,000 nodes) in the FEM based forward model. The FEM model needs to be solved multiple times in the line search at each iteration of the Gauss-Newton method and with the mesh used each forward solution takes approximately 80s computation time.

§ CONCLUSIONS

In this work, we presented the first 3D absolute EIT reconstructions from CGO-based methods on experimental 3D tank data, and compared them to the current standard, a total variation regularized non-linear least squares approach.  We demonstrated that, with correct domain modeling, quality 3D absolute reconstructions can be obtained by all of the methods, comparable to the quality seen in linear difference imaging. All methods, Calderón, $\texp$, $\tzero$, and TV reasonably handled the moderate domain modeling error within little noticeable change in localization error and target contrast.  For the large modeling error case, the $\texp$ and $\tzero$ methods correctly identified the targets with high contrast, additional artefacts were introduced into the Calderón reconstruction, and the error proved too significant for the TV method.  The computational cost of the CGO reconstruction is trivial compared to TV (non-optimized, less than 1 sec/recon for Calderón, approximately 5 sec/recon for $\texp$ and $\tzero$, compared to 2-3 hours per reconstruction for TV).

§ ACKNOWLEDGMENTS

Research reported in this paper was supported by the National Institute of Biomedical Imaging and Bioengineering of the National Institutes of Health under award numbers R21EB028064 (SH and JN) 1R01EB026710-01A1 (GS, DI, JN, ORS, and the development of the ACT5 device). The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.  JT and VK were supported by the Academy of Finland (Project 336791, Finnish Centre of Excellence in Inverse modeling and Imaging), the Jane and Aatos Erkko Foundation and Neurocenter Finland.