Fast absolute 3D CGO-based electrical impedance tomography on experimental tank data

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óns Problem' by much of the mathematical community since A. Calderón formulated this inverse problem as follows (Calderón 1980): In two or more dimensions, can one find the conductivity, σ(x), inside a body, Ω, from all possible electrical measurements made on the surface, ∂Ω, 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 Maxwells equations, which we will refer to as the conductivity equation, with a Neumann boundary condition:

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\nabla }}\cdot \sigma (x){\rm{\nabla }}u(x) & = & 0,\qquad x\in {\rm{\Omega }}\subset {{\mathbb{R}}}^{3}\\ \sigma (x)\displaystyle \frac{\partial u(x)}{\partial \nu (x)} & = & j(x),\qquad x\in \partial {\rm{\Omega }}.\end{array}\end{eqnarray} \tag{ 2.1 }$$

Here ν = ν(x) denotes the outward pointing unit normal to the surface at x ∈ ∂Ω, and, $\tfrac{\partial u(x)}{\partial \nu (x)}=\nu (x)\cdot {\rm{\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 ∈ ∂Ω, and u(x) solves the conductivity equation (2.1) 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, ${{\rm{\Lambda }}}_{\sigma }:= {{ \mathcal R }}_{\sigma }^{-1}$, can one find σ(x)? Calderón showed that if σ(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^ζ·x , to the conductivity equation with constant conductivity, i.e. the Laplace equation, where ζ is a complex valued vector with ζ · ζ = 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óns solutions to what we now call Complex Geometrical Optics (CGO) solutions, introduced by Sylvester and Uhlmann (1986, 1987) 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 ∣ζ∣ → ∞ . The ${{\bf{t}}}^{{\bf{\exp }}}$ and t⁰ methods described below in section 2.3 follow this strategy, as opposed to the 2D D-bar methods described in Mueller and Siltanen (2020). 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 ζ 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 Nachman (1988), Novikov (1988).

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 Nachman (1996) for $\sigma \in {{ \mathcal C }}^{2}({\rm{\Omega }})$ and later in Astala and Päivärinta (2006) for σ ∈ L^∞(Ω). Other pioneering work proving uniqueness under a variety of hypotheses on the conductivity include Langer (1933), Tikhonov (1949), Druskin (1982), Kohn and Vogelius (1984), Parker (1984), Novikov (1988). Extensive references to more recent progress in the analysis and numerical analysis of the Calderón problem can be found in Borcea (2002), Mueller and Siltanen (2012), Delbary and Knudsen (2014), Uhlmann (2014), Feldman et al (2019), Mueller and Siltanen (2020), Knudsen and Rasmussen (2022).

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_ℓ , ℓ = 1, 2,...,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,...,L − 1, to the L electrodes, and measure the resulting L − 1 voltage patterns, ${\vec{V}}^{(k)}$, where ${\vec{I}}_{{\ell }}^{(k)}$ and ${\vec{V}}_{{\ell }}^{(k)}$ denote the applied current, and measured voltage from the kth pattern, on the ℓth electrode, for ℓ = 1, 2,...,L. From conservation of charge, and our choice of ground, we assume

$$\begin{eqnarray*}&&\displaystyle \sum _{{\ell }=1}^{L}{\vec{I}}_{{\ell }}^{(k)}=\displaystyle \sum _{{\ell }=1}^{L}{\vec{V}}_{{\ell }}^{(k)}=0.\end{eqnarray*}$$

The kth voltage or potential, u^(k)(x), resulting inside the body is determined by the conductivity equation, ∇ · σ∇u^(k) = 0, and the 'Complete Electrode Model', (Cheng et al 1989, Somersalo et al 1992) where:

$$\begin{eqnarray}\begin{array}{rcl}{\int }_{{e}_{{\ell }}}\sigma (x)\displaystyle \frac{\partial {u}^{(k)}(x)}{\partial \nu (x)}{dS}(x) & = & {I}_{{\ell }}^{(k)},\qquad x\in {e}_{{\ell }}\\ \sigma (x)\displaystyle \frac{\partial {u}^{(k)}(x)}{\partial \nu (x)} & = & 0,\qquad x\notin \bigcup _{{\ell }=1}^{L}{e}_{{\ell }}\\ {u}^{(k)}(x)+{z}_{{\ell }}\sigma (x)\displaystyle \frac{\partial {u}^{(k)}(x)}{\partial \nu (x)} & = & {V}_{{\ell }}^{(k)},\qquad x\in {e}_{{\ell }}.\end{array}\end{eqnarray} \tag{ 2.2 }$$

Here z_ℓ is the effective contact, or surface, impedance on the ℓ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 σ are denoted by the L × L matrices, R_σ , where ${{\bf{R}}}_{\sigma }{\vec{I}}^{(k)}={\vec{V}}^{(k)}$, and we define ${{\bf{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 ${{\bf{L}}}_{\sigma }:= {{\bf{R}}}_{\sigma }^{-1}$.

Following Calderón's original paper (Calderón 1980), Calderón's method approximates the conductivity, σ(x), from its Fourier transform. Here we present a brief description of the method and refer the reader to Calderón (1980), Bikowski and Mueller (2008), Muller et al (2017), Hamilton et al (2021) for further details. Calderón's method assumes the conductivity is a small perturbation, δ σ(x), from a constant background, σ_b , i.e. σ(x) = σ_b + δ σ(x). In this paper, we assume that the background conductivity σ_b = 1. If the background constant is not one, then the problem can be scaled and unscaled as in Isaacson et al (2004), Hamilton et al (2021).

The three steps of Calderón's method in 3-D, as described in Hamilton et al (2021), are:

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

$$\begin{eqnarray}&&\widehat{\delta \sigma }(z)\approx \hat{{\rm{F}}}(z):= -\displaystyle \frac{1}{2{\pi }^{2}| z{| }^{2}}{\int }_{\partial {\rm{\Omega }}}{e}^{\pi i(z\cdot x)+\pi (a\cdot x)}\left({{\rm{\Lambda }}}_{\sigma }-{{\rm{\Lambda }}}_{1}\right){e}^{\pi i(z\cdot x)-\pi (a\cdot x)}{dS}(x),\end{eqnarray} \tag{ 2.3 }$$

where z and a satisfy

$$\begin{eqnarray}&&z,a\in {{\mathbb{R}}}^{3},| z| =| a| ,\mathrm{and}\,z\cdot a=0,\end{eqnarray} \tag{ 2.4 }$$

and Λ₁ is the DN map for a constant conductivity of 1.

Step 2: Take the inverse Fourier transform of $\hat{{\rm{F}}}(z)$:

$$\begin{eqnarray}&&\delta {\sigma }^{{\bf{CAL}}}(x)\approx {{ \mathcal F }}^{-1}\{\hat{{\rm{F}}}(z)\}(x)={\int }_{{{\mathbb{R}}}^{3}}\hat{{\rm{F}}}(z){e}^{-2\pi i(x\cdot z)}{dz}.\end{eqnarray} \tag{ 2.5 }$$

Step 3: Add the background to the perturbation to recover the approximate conductivity, σ^CAL (x):

$$\begin{eqnarray}&&{\sigma }^{{\bf{CAL}}}(x)={\sigma }_{b}+\delta {\sigma }^{{\bf{CAL}}}(x).\end{eqnarray} \tag{ 2.6 }$$

The definition of the Fourier transform in Calderón's method is different from that used in the ${{\bf{t}}}^{{\bf{\exp }}}$ and t⁰ 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 $\hat{{\rm{F}}}(z)$ in spherical Fourier coordinates. As such, we choose

$$\begin{eqnarray*}&&z=| z| (\cos \tilde{\phi }\sin \tilde{\theta },\sin \tilde{\phi }\sin \tilde{\theta },\cos \tilde{\theta })\,\mathrm{and}\,a=| z| (\cos \tilde{\phi }\cos \tilde{\theta },\sin \tilde{\phi }\cos \tilde{\theta },-\sin \tilde{\theta }),\end{eqnarray*}$$

for ∣z∣ ≥ 0, $0\leqslant \tilde{\phi }\leqslant 2\pi$ and $0\leqslant \tilde{\theta }\leqslant \pi$ so that z and a satisfy (2.4). Then, the inverse Fourier transform in (2.5) becomes

$$\begin{eqnarray}&&\delta {\sigma }^{{\bf{CAL}}}(x)={\int }_{0}^{\infty }{\int }_{0}^{2\pi }{\int }_{0}^{\pi }| z{| }^{2}\sin \tilde{\theta }\hat{{\rm{F}}}(| z| ,\tilde{\phi },\tilde{\theta }){e}^{-2\pi i(x\cdot z)}d\tilde{\theta }d\tilde{\phi }d| z| .\end{eqnarray} \tag{ 2.7 }$$

Additionally, we implement the use of a mollifier, $\hat{\eta }\left(\tfrac{z}{y}\right)$, as introduced in Calderón (1980) for some parameter $y\in {\mathbb{R}}$ to reduce Gibbs phenomenon caused by jump discontinuities in σ(x) while recovering δ σ^CAL

$$\begin{eqnarray}&&\delta {\sigma }^{{\bf{CAL}}}(x)={\int }_{0}^{\infty }{\int }_{0}^{2\pi }{\int }_{0}^{\pi }| z{| }^{2}\sin \tilde{\theta }\hat{{\rm{F}}}(| z| ,\tilde{\phi },\tilde{\theta })\hat{\eta }\left(\displaystyle \frac{z}{y}\right){e}^{-2\pi i(x\cdot z)}d\tilde{\theta }d\tilde{\phi }d| z| .\end{eqnarray} \tag{ 2.8 }$$

We implement the same mollifier as used in Bikowski and Mueller (2008), Hamilton et al (2021),

$$\begin{eqnarray}&&\hat{\eta }\left(\displaystyle \frac{z}{y}\right)={e}^{-\pi t| z{| }^{2}},\end{eqnarray} \tag{ 2.9 }$$

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 σ(x).

Since noise causes (2.3) to blow up at large ∣z∣, we use a non-uniform truncation regularization strategy. A similar regularization strategy was proved stable for the 2D D-bar method in Knudsen et al (2009), which also noted non-uniform truncation also produces reliable reconstructions. In our case, we will first compute $\hat{{\rm{F}}}(z)$ for ∣z∣ within an outer radius of ${T}_{{z}_{2}}$. We keep values of $\hat{{\rm{F}}}(z)$ whose real and imaginary amplitudes are below a threshold determined by the amplitudes of $\hat{{\rm{F}}}$ within a smaller radius $| z| \leqslant {T}_{{z}_{1}};$ $\hat{{\rm{F}}}$ 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 $\hat{{\rm{F}}}$ to blow up and ${T}_{{z}_{2}}$ is chosen to keep as much reasonable information from $\hat{{\rm{F}}}$ without introducing holes in the non-zero region of $\hat{{\rm{F}}}$. As such, our truncated $\hat{{\rm{F}}}$ is computed by

$$\begin{eqnarray}{\hat{{\rm{F}}}}_{R}(z)=\left\{\begin{array}{ll}\hat{{\rm{F}}}(z), & \mathrm{if}| z| \leqslant {T}_{{z}_{2}},\mathrm{and}\,| \mathrm{Re}(\hat{{\rm{F}}}(z))| \leqslant \mathop{\max }\limits_{| z| \leqslant {T}_{{z}_{1}}}| \mathrm{Re}(\hat{{\rm{F}}}(z))| ,\,| \mathrm{Im}(\hat{{\rm{F}}}(z))| \leqslant \mathop{\max }\limits_{| z| \leqslant {T}_{{z}_{1}}}| \mathrm{Im}(\hat{{\rm{F}}}(z))| \\ 0, & \mathrm{else}.\end{array}\right.\end{eqnarray} \tag{ 2.10 }$$

With our truncated $\hat{{\rm{F}}}$, equation (2.8) is truncated in the radial variable, leading to the approximation

$$\begin{eqnarray}&&\delta {\sigma }_{R}^{{\bf{CAL}}}(x)={\int }_{0}^{{T}_{{z}_{2}}}{\int }_{0}^{2\pi }{\int }_{0}^{\pi }| z{| }^{2}\sin \tilde{\theta }{\hat{{\rm{F}}}}_{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{eqnarray} \tag{ 2.11 }$$

For difference images shown in this paper, we only perform Steps 1 and 2 of the method and replace Λ₁ in (2.3) with a reference DN map, ${{\rm{\Lambda }}}_{{\sigma }_{\mathrm{ref}}}$ before computing (2.11). Thus, the flow for absolute images is

$$\begin{eqnarray*}&&\left({{\rm{\Lambda }}}_{\sigma },{{\rm{\Lambda }}}_{1}\right)\mathop{\longrightarrow }\limits^{1}\hat{{\rm{F}}}(z)\mathop{\longrightarrow }\limits^{2}\delta {\sigma }^{{\bf{CAL}}}\mathop{\longrightarrow }\limits^{3}{\sigma }^{{\bf{CAL}}}\end{eqnarray*}$$

and the flow for difference images in this paper is

$$\begin{eqnarray*}&&\left({{\rm{\Lambda }}}_{\sigma },{{\rm{\Lambda }}}_{{\sigma }_{\mathrm{ref}}}\right)\mathop{\longrightarrow }\limits^{1}\hat{{\rm{F}}}(z)\mathop{\longrightarrow }\limits^{2}\delta {\sigma }^{{\bf{CAL}}}.\end{eqnarray*}$$

2.2.1. Numerical implementation

In this section, we review the implementation details of Calderón's method introduced in Hamilton et al (2021).

For Step 1, we compute (2.3) for $| z| \leqslant {T}_{{z}_{2}}$ by discretizing the boundary integral as follows,

$$\begin{eqnarray}\begin{array}{rcl}\hat{{\rm{F}}}(z) & = & -\displaystyle \frac{1}{2{\pi }^{2}| z{| }^{2}}{\int }_{\partial {\rm{\Omega }}}{e}^{\pi i(z\cdot x)+\pi (a\cdot x)}\left({{\rm{\Lambda }}}_{\sigma }-{{\rm{\Lambda }}}_{1}\right){e}^{\pi i(z\cdot x)-\pi (a\cdot x)}{dS}(x)\\ & \approx & -\displaystyle \frac{1}{2{\pi }^{2}| z{| }^{2}}\left(\displaystyle \frac{| \partial {\rm{\Omega }}| }{L}\right){\left({e}^{\pi i(z\cdot {\bf{x}})+\pi (a\cdot {\bf{x}})}\right)}^{{\bf{T}}}{\bf{Q}}\left({{\bf{L}}}_{\sigma }-{{\bf{L}}}_{1}\right){{\bf{Q}}}^{{\bf{T}}}\left[{e}^{\pi i(z\cdot {\bf{x}})-\pi (a\cdot {\bf{x}})}\right],\end{array}\end{eqnarray} \tag{ 2.12 }$$

where ∣∂Ω∣ is the surface area of the domain; L is the number of electrodes; ${\bf{x}}\in {{\mathbb{R}}}^{(L\times 3)}$ denotes the vector of the Cartesian centers of the L electrodes; ^T denotes the traditional, non-conjugate, matrix transpose; L_σ and L₁ denote the discrete matrix approximations to the DN maps Λ_σ and Λ₁ respectively; and ${\bf{Q}}\in {{\mathbb{R}}}^{L\times {N}_{{li}}}$ denotes an orthonormal basis created using N_li linearly independent applied currents over L electrodes as was done in Hamilton et al (2021). The matrix L₁ is based on the FEM solution of the CEM (2.2). Problem-specific mesh details are given below in section 2.6.

We then compute ${\hat{{\rm{F}}}}_{R}$ according to equations (2.10) and (2.12). For Step 2, on an equally-spaced 16 × 16 × 16 rectangular grid in x, the conductivity difference $\delta {\sigma }_{R}^{{\bf{CAL}}}(x)$ is computed via (2.11) 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 (2.11) and (2.12), replacing L₁ with a discrete reference DN map, ${{\bf{L}}}_{{\sigma }_{\mathrm{ref}}}$ in (2.12). 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}^{{\bf{CAL}}}(x)$ in this paper are produced using (2.6) replacing σ_b with σ_best

$$\begin{eqnarray*}&&{\sigma }_{R}^{{\bf{CAL}}}(x)={\sigma }_{\mathrm{best}}+\delta {\sigma }_{R}^{{\bf{CAL}}}(x).\end{eqnarray*}$$

The solution is then interpolated to a 64 × 64 × 64 rectangular grid.

Following Isaacson et al (2004), σ_best is given by

$$\begin{eqnarray}&&{\sigma }_{\mathrm{best}}=\displaystyle \frac{{\sum }_{k=1}^{K}{\sum }_{{\ell }=1}^{L}{U}_{{\ell }}^{k}{U}_{{\ell }}^{k}}{{\sum }_{k=1}^{K}{\sum }_{{\ell }=1}^{L}{U}_{{\ell }}^{k}{V}_{{\ell }}^{k}},\end{eqnarray} \tag{ 2.13 }$$

where ${U}_{{\ell }}^{k}$ is the k^th simulated voltage pattern measured on electrode ℓwith a homogeneous conductivity of 1 S m⁻¹ and ${V}_{{\ell }}^{k}$ is the k^th voltage pattern measured on electrode ℓfor the inhomogeneous conductivity σ. The simulated voltages are the same voltages used to compute L₁, as described in section 2.6.1.

Both the ${{\bf{t}}}^{{\bf{\exp }}}$ and t⁰ methods are derived from the constructive proofs of Nachman (1988), Novikov and Khenkin (1987), Novikov (1988) and involve special solutions called Complex Geometrical Optics (CGO) solutions (Sylvester and Uhlmann 1987). A brief summary is included here for the reader's convenience. For further details see Delbary et al (2012), Hamilton et al (2021), Delbary and Knudsen (2014), Nachman (1988).

Assuming that the conductivity is a constant σ_c = 1 in a neighborhood of the boundary ∂Ω, the real-valued conductivity equation (2.1) can be transformed to the Schrödinger equation

$$\begin{eqnarray}&&\left(-{\rm{\Delta }}+q(x)\right)\widetilde{u}(x)=0,\qquad x\in {{\mathbb{R}}}^{3},\end{eqnarray} \tag{ 2.14 }$$

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

$$\begin{eqnarray*}&&\left(-{\rm{\Delta }}+q(x)\right)\psi (x,\zeta )=0,\qquad x\in {{\mathbb{R}}}^{3},\end{eqnarray*}$$

where ψ(x, ζ) ∼ e^ix·ζ for large ∣x∣ or ∣ζ∣, and

$$\begin{eqnarray}&&{{ \mathcal V }}_{\xi }=\left\{\zeta \in {{\mathbb{C}}}^{3}| {\zeta }^{2}=0,\,{\left(\xi +\zeta \right)}^{2}=0\right\},\quad \mathrm{for}\,\mathrm{each}\,\xi \in {{\mathbb{R}}}^{3},\end{eqnarray} \tag{ 2.15 }$$

where ζ² = ζ · ζ and ζ is a purely auxiliary parameter. The conductivity σ(x) can then be recovered from the DN map Λ_σ as follows.

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

$$\begin{eqnarray}&&\psi (x,\zeta )={e}^{{ix}\cdot \zeta }-{\int }_{\partial {\rm{\Omega }}}{G}_{\zeta }(x-y)\left({{\rm{\Lambda }}}_{\sigma }-{{\rm{\Lambda }}}_{1}\right)\psi (y,\zeta ){dS}(y),\end{eqnarray} \tag{ 2.16 }$$

where

$$\begin{eqnarray*}&&{G}_{\zeta }(x)=\displaystyle \frac{{e}^{{ix}\cdot \zeta }}{{\left(2\pi \right)}^{3}}{\int }_{{{\mathbb{R}}}^{3}}\displaystyle \frac{{e}^{{ix}\cdot k}}{| k{| }^{2}+2k\cdot \zeta }{dk},\quad x\in {{\mathbb{R}}}^{3}\setminus \left\{0\right\},\end{eqnarray*}$$

denotes the Faddeev Green's function (Faddeev 1966). Then, evaluate the scattering data

$$\begin{eqnarray}&&{\bf{t}}(\xi ,\zeta )={\int }_{\partial {\rm{\Omega }}}{e}^{-{ix}\cdot (\xi +\zeta )}\left({{\rm{\Lambda }}}_{\sigma }-{{\rm{\Lambda }}}_{1}\right)\psi (x,\zeta ){dS}(x).\end{eqnarray} \tag{ 2.17 }$$

For ∣ζ∣ large, the Schrödinger potential q(x) can be recovered via the inverse Fourier transform

$$\begin{eqnarray}&&q(x)\approx {{ \mathcal F }}^{-1}\left\{{\bf{t}}(\xi ,\zeta )\right\}(x)=\displaystyle \frac{1}{{\left(2\pi \right)}^{3}}{\int }_{{{\mathbb{R}}}^{3}}{e}^{{ix}\cdot \xi }{\bf{t}}(\xi ,\zeta )d\xi ,\quad x\in {{\mathbb{R}}}^{3}.\end{eqnarray} \tag{ 2.18 }$$

The conductivity is then recovered by solving the boundary value problem

$$\begin{eqnarray}\left\{\begin{array}{ll}(-{\rm{\Delta }}+q(x))\tilde{u}(x)=0 & x\in {\rm{\Omega }}\subset {{\mathbb{R}}}^{3}\\ \tilde{u}(x)=1 & x\in \partial {\rm{\Omega }},\end{array}\right.\end{eqnarray} \tag{ 2.19 }$$

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 ψ(x, ζ) by their asymptotic behavior e^ix·ζ in the scattering data (2.17) via

$$\begin{eqnarray}&&{{\bf{t}}}^{{\bf{\exp }}}(\xi ,\zeta )={\int }_{\partial {\rm{\Omega }}}{e}^{-{ix}\cdot (\xi +\zeta )}\left({{\rm{\Lambda }}}_{\sigma }-{{\rm{\Lambda }}}_{1}\right){e}^{{ix}\cdot \zeta }\,{dS}(x).\end{eqnarray} \tag{ 2.20 }$$

yields a 'Born approximation' typically called the ${{\bf{t}}}^{{\bf{\exp }}}$ approximation for EIT. Using this approximate scattering data in place of the fully nonlinear t(ξ, ζ), one proceeds with the recovery of an approximate potential ${{\bf{q}}}^{{\bf{\exp }}}(x)$ via (2.18) and conductivity ${\sigma }^{{\bf{\exp }}}(x)$ via (2.19). The flow is:

$$\begin{eqnarray*}&&\left({{\rm{\Lambda }}}_{\sigma },{{\rm{\Lambda }}}_{1}\right)\mathop{\longrightarrow }\limits^{1}{{\bf{t}}}^{{\bf{\exp }}}(\xi ,\zeta )\mathop{\longrightarrow }\limits^{2}{{\bf{q}}}^{{\bf{\exp }}}(x)\mathop{\longrightarrow }\limits^{3}{\sigma }^{{\bf{\exp }}}(x).\end{eqnarray*}$$

An intermediate approximation can be computed by replacing the Faddeev Green's function G_ζ (x) in the single layer potential, in (2.16), for the traces of the CGOs, with the standard Green's function ${G}_{0}(x)=\tfrac{1}{4\pi | x| }$ for the Laplacian operator. Thus, one solves

$$\begin{eqnarray}&&{\psi }^{0}(x,\zeta )={e}^{{ix}\cdot \zeta }-{\int }_{\partial {\rm{\Omega }}}{G}_{0}(x-y)\left({{\rm{\Lambda }}}_{\sigma }-{{\rm{\Lambda }}}_{1}\right){\psi }^{0}(y,\zeta ){dS}(y),\end{eqnarray} \tag{ 2.21 }$$

for the CGOs ψ⁰(x, ζ), avoiding the exponentially growing Faddeev Green's function. A corresponding approximation to the scattering data is then computed by using ψ⁰(x, ζ) in place of ψ(x, ζ) in (2.17) and then continuing to recover q⁰(x) and σ⁰(x). The flow is then

$$\begin{eqnarray*}&&\left({{\rm{\Lambda }}}_{\sigma },{{\rm{\Lambda }}}_{1}\right)\mathop{\longrightarrow }\limits^{1}{\psi }^{0}(x,\zeta )\mathop{\longrightarrow }\limits^{2}{{\bf{t}}}^{0}(\xi ,\zeta )\mathop{\longrightarrow }\limits^{3}{{\bf{q}}}^{0}(x)\mathop{\longrightarrow }\limits^{4}{\sigma }^{0}(x).\end{eqnarray*}$$

We point out that the methods, as outlined above, assumed that the conductivity was a constant σ_c = 1 near the boundary of the domain. As mentioned in section 2.2, the problem can be scaled and unscaled as in Isaacson et al (2004), Hamilton et al (2021) when the constant σ_c is not one. In practice, we estimate the best-fit constant conductivity fit to the data σ_best as given by (2.13). Explicitly, as in (Hamilton et al 2021), we scale the DN map by using $\tfrac{1}{{\sigma }_{\mathrm{best}}}{{\rm{\Lambda }}}_{\sigma }$ in place of Λ_σ , and re-scale at the end using ${\sigma }^{{\bf{\exp }}}(x)={\sigma }_{\mathrm{best}}\left({\tilde{u}}^{\exp }(x)\right)$ and ${\sigma }^{0}(x)={\sigma }_{\mathrm{best}}\left({\tilde{u}}^{0}(x)\right)$.

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

2.3.1. Numerical implementation

Here we provide the numerical details pertinent to the implementation of the ${{\bf{t}}}^{{\bf{\exp }}}$ and t⁰ algorithms outlined above. As with the Calderón method above, the main idea is to expand functions in the same orthonormal basis Q as described in section 2.2.

Following DeAngelo and Mueller (2010), for each electrode center x_ℓ , expand ψ(x, ζ) and e^ix·ζ as

$$\begin{eqnarray}&&\psi ({x}_{{\ell }},\zeta )\approx \displaystyle \sum _{j=1}^{{N}_{{li}}}{b}_{j}(\zeta ){{\bf{Q}}}_{{\ell }}^{j},\qquad {e}^{{{ix}}_{{\ell }}\cdot \zeta }\approx \displaystyle \sum _{j=1}^{{N}_{{li}}}{c}_{j}(\zeta ){{\bf{Q}}}_{{\ell }}^{j},\qquad {\ell }=1,\ldots L,\end{eqnarray} \tag{ 2.22 }$$

where ${{\bf{Q}}}_{{\ell }}^{j}$ denotes the (ℓ, j) entry of Q. Then, the boundary integral equation (2.21) can be approximated as follows

$$\begin{eqnarray}\begin{array}{rcl}{\psi }^{0}({x}_{{\ell }},\zeta ) & = & ={e}^{{{ix}}_{{\ell }}\cdot \zeta }-{\int }_{\partial {\rm{\Omega }}}{G}_{0}({x}_{{\ell }}-y)\left({{\rm{\Lambda }}}_{\sigma }-{{\rm{\Lambda }}}_{1}\right){\psi }^{0}(y,\zeta ){dS}(y)\\ & \approx & {e}^{{{ix}}_{{\ell }}\cdot \zeta }-\displaystyle \sum _{{\ell }^{\prime} =1}^{L}{\int }_{{E}_{{\ell }^{\prime} }}{G}_{0}({x}_{{\ell }}-y)\left({{\rm{\Lambda }}}_{\sigma }-{{\rm{\Lambda }}}_{1}\right){\psi }^{0}(y,\zeta ){dS}(y)\\ & \approx & {e}^{{{ix}}_{{\ell }}\cdot \zeta }-\left[\displaystyle \sum _{{\ell }^{\prime} =1}^{L}{\int }_{{E}_{{\ell }^{\prime} }}{G}_{0}({x}_{{\ell }}-y){dS}(y)\right]\left[\left({{\bf{L}}}_{\sigma }-{{\bf{L}}}_{1}\right){\psi }^{0}({y}_{{\ell }^{\prime} },\zeta )\right],\end{array}\end{eqnarray} \tag{ 2.23 }$$

where ${E}_{{\ell }^{\prime} }$ denotes the ${\ell }^{\prime} \mathrm{th}$ extended electrode where ${\bigcup }_{{\ell }^{\prime} =1}^{L}{E}_{{\ell }^{\prime} }=\partial {\rm{\Omega }}$ and the ${E}_{{\ell }^{\prime} }$ are mutually disjoint (Hyvönen 2009). Note the true electrodes need not cover the surface ∂Ω, only the extended (mathematical) electrodes that we will use to discretize the integral. Then, using the expansions from (2.22) in (2.23), we have

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \sum _{j=1}^{{N}_{{li}}}{b}_{j}(\zeta ){{\bf{Q}}}_{{\ell }}^{j} & \approx & \displaystyle \sum _{j=1}^{{N}_{{li}}}{c}_{j}(\zeta ){{\bf{Q}}}_{{\ell }}^{j}-\left[\displaystyle \sum _{{\ell }^{\prime} =1}^{L}{\int }_{{E}_{{\ell }^{\prime} }}{G}_{0}({x}_{{\ell }}-y){dS}(y)\right]\left[\left({{\bf{L}}}_{\sigma }-{{\bf{L}}}_{1}\right)\displaystyle \sum _{j=1}^{{N}_{{li}}}{b}_{j}(\zeta ){{\bf{Q}}}_{{\ell }^{\prime} }^{j}\right]\\ & = & \displaystyle \sum _{j=1}^{{N}_{{li}}}{c}_{j}(\zeta ){{\bf{Q}}}_{{\ell }}^{j}-\left[\displaystyle \sum _{{\ell }^{\prime} =1}^{L}{\int }_{{E}_{{\ell }^{\prime} }}{G}_{0}({x}_{{\ell }}-y){dS}(y)\right]\left[\displaystyle \sum _{j=1}^{{N}_{{li}}}{b}_{j}(\zeta ){f}_{j}\left({y}_{{\ell }^{\prime} }\right)\right],\end{array}\end{eqnarray} \tag{ 2.24 }$$

where ${f}_{j}\left({y}_{{\ell }^{\prime} }\right)$ represents the action of $\left({{\bf{L}}}_{\sigma }-{{\bf{L}}}_{1}\right)$ on Q^j evaluated at ${y}_{{\ell }^{\prime} }$, which can be computed as the $({\ell }^{\prime} ,j)$ entry of ${\bf{Q}}\left({{\bf{L}}}_{\sigma }-{{\bf{L}}}_{1}\right)$. Define

$$\begin{eqnarray}{\widetilde{{\bf{G}}}}_{0}({\ell },{\ell }^{\prime} )=\left\{\begin{array}{ll}\tfrac{1}{4\pi \left|{x}_{{\ell }}-{y}_{{\ell }^{\prime} }\right|} & {\ell }\ne {\ell }^{\prime} \\ 0 & {\ell }={\ell }^{\prime} ,\end{array}\right.\end{eqnarray} \tag{ 2.25 }$$

where we have removed the singularities at x_ℓ = y_ℓ . Assuming $| {E}_{{\ell }}| =\tfrac{| \partial {\rm{\Omega }}| }{L}$ for each ℓ = 1,...,L, and using (2.25) in (2.24) we find

$$\begin{eqnarray*}\begin{array}{rcl}\displaystyle \sum _{j=1}^{{N}_{{li}}}{b}_{j}(\zeta ){{\bf{Q}}}_{{\ell }}^{j} & \approx & \displaystyle \sum _{j=1}^{{N}_{{li}}}{c}_{j}(\zeta ){{\bf{Q}}}_{{\ell }}^{j}-\displaystyle \frac{| \partial {\rm{\Omega }}| }{L}\displaystyle \sum _{j=1}^{{N}_{{li}}}{b}_{j}(\zeta )\displaystyle \sum _{{\ell }^{\prime} =1}^{L}{\widetilde{{\bf{G}}}}_{0}({\ell },{\ell }^{\prime} ){f}_{j}\left({y}_{{\ell }^{\prime} }\right)\\ & \approx & \displaystyle \sum _{j=1}^{{N}_{{li}}}{c}_{j}(\zeta ){{\bf{Q}}}_{{\ell }}^{j}-\displaystyle \frac{| \partial {\rm{\Omega }}| }{L}\displaystyle \sum _{j=1}^{{N}_{{li}}}{b}_{j}(\zeta )\left[{\widetilde{{\bf{G}}}}_{0}{\bf{Q}}\left({{\bf{L}}}_{\sigma }-{{\bf{L}}}_{1}\right)\right]({\ell },j)\end{array}\end{eqnarray*}$$

or in matrix form,

$$\begin{eqnarray*}&&{\bf{Q}}\vec{b}={\bf{Q}}\vec{c}-\displaystyle \frac{| \partial {\rm{\Omega }}| }{L}{\widetilde{{\bf{G}}}}_{0}{\bf{Q}}\left({{\bf{L}}}_{\sigma }-{{\bf{L}}}_{1}\right)\vec{b}.\end{eqnarray*}$$

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

$$\begin{eqnarray}&&\left(I+A\right)\vec{b}=\vec{c},\end{eqnarray} \tag{ 2.26 }$$

where $A=\tfrac{| \partial {\rm{\Omega }}| }{L}{{\bf{Q}}}^{{\bf{T}}}{\widetilde{{\bf{G}}}}_{0}{\bf{Q}}\left({{\bf{L}}}_{\sigma }-{{\bf{L}}}_{1}\right)$, and I is the identity matrix. If the extended electrodes E_ℓ are not uniform in size, then one could compute as weighted sum replacing the uniform weight $| {E}_{{\ell }}| =\tfrac{| \partial {\rm{\Omega }}| }{L}$ as appropriate.

Next, following Hamilton et al (2021), the scattering data t⁰(ξ, ζ(ξ)) can be computed for all ∣ξ∣ less than a chosen truncation radius T_ξ via

$$\begin{eqnarray}{{\bf{t}}}^{0}(\xi ,\zeta )\approx \left\{\begin{array}{ll}\tfrac{| \partial {\rm{\Omega }}| }{L}{\left[{e}^{-i{\bf{x}}\cdot (\xi +\zeta )}\right]}^{{\bf{T}}}{\bf{Q}}\left({{\bf{L}}}_{\sigma }-{{\bf{L}}}_{1}\right)\vec{b} & | \xi | \lt {T}_{\xi }\\ 0 & {else},\end{array}\right.\end{eqnarray} \tag{ 2.27 }$$

since $\vec{b}={{\bf{Q}}}^{{\bf{T}}}{\psi }^{0}({\bf{x}},\zeta )$, where ${\bf{x}}\in {{\mathbb{R}}}^{(L\times 3)}$ is the same as in section 2.2.1, i.e. a vector storing the centers of the electrodes.

Next, the approximate potential q⁰ is recovered by computing the inverse Fourier transform of the truncated scattering data t⁰ using a Simpson's rule in 3D,

$$\begin{eqnarray*}&&{{\bf{q}}}^{0}(x)=\displaystyle \frac{1}{{\left(2\pi \right)}^{3}}{\int }_{{\left[-{T}_{\xi },{T}_{\xi }\right]}^{3}}{e}^{{ix}\cdot \xi }{{\bf{t}}}^{0}(\xi ,\zeta (\xi ))d\xi .\end{eqnarray*}$$

Alternatively, one could use an IFFT to achieve additional speedup, taking care with quadrature points and the particular form of the kernel. Following Hamilton et al (2021), the conductivity σ⁰(x) was recovered by first solving the boundary value problem (2.19), using the PDE toolbox in Matlab using a mesh with approximately 21, 000 3D elements, then computing ${\sigma }^{0}(x)=\sigma \mathrm{best}{\left(\tilde{u}(x)\right)}^{2}$. The 3D visualizations were obtained by interpolation to a 64 × 64 × 64 rectangular grid using Matlab's scatteredInterpolant function. The ζ values were computed following the minimal-zeta approach outlined in Delbary and Knudsen (2014). We remark that while the non-existence of exceptional points for the solution of the boundary integral equation (2.16) is proven for large ∣ζ∣ (Sylvester and Uhlmann 1987, Nachman 1988), as well as small ∣ζ∣ (Cornean et al 2006); it is still an open question for the intermediate values required here to perform the computation on a computer.

The ${{\bf{t}}}^{{\bf{\exp }}}$ reconstructions were obtained in an analogous fashion to those of t⁰, this time bypassing the boundary integral equation (2.17) and directly computing

$$\begin{eqnarray*}{{\bf{t}}}^{{\bf{\exp }}}(\xi ,\zeta )\approx \left\{\begin{array}{ll}\displaystyle \frac{| \partial {\rm{\Omega }}| }{L}{\left[{e}^{-i{\bf{x}}\cdot (\xi +\zeta )}\right]}^{{\bf{T}}}{\bf{Q}}\left({{\bf{L}}}_{\sigma }-{{\bf{L}}}_{1}\right){{\bf{Q}}}^{{\bf{T}}}\left[{e}^{i{\bf{x}}\cdot \zeta }\right] & | \xi | \lt {T}_{\xi }\\ 0 & {else},\end{array}\right.\end{eqnarray*}$$

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

Difference imaging can be performed with the t⁰ and ${{\bf{t}}}^{{\bf{\exp }}}$ methods by replacing the matrix L₁ with ${{\bf{L}}}_{{\sigma }_{\mathrm{ref}}}$ and computing σ_diff(x)=σ(x) − σ_best in the final step.

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.

2.4.1. 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{eqnarray}&&\hat{\sigma }=\arg \mathop{\min }\limits_{\sigma \gt 0}\left\{{\parallel {L}_{e}\left(V-U(\sigma ,{z}_{* })\right)\parallel }^{2}+\alpha {\mathrm{TV}}_{\beta }(\sigma )\right\},\end{eqnarray} \tag{ 2.28 }$$

where U(σ, z) is the finite dimensional forward map, ${z}_{* }\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}={{\rm{\Gamma }}}_{e}^{-1}$, scalar valued α is the regularization parameter and TV_β (σ) is the (smooth) TV regularization functional (Rudin et al 1992)

$$\begin{eqnarray}&&{\mathrm{TV}}_{\beta }(\sigma )={\int }_{{\rm{\Omega }}}\sqrt{\parallel {\rm{\nabla }}\sigma {\parallel }^{2}+\beta }\ {\rm{d}}r,\end{eqnarray} \tag{ 2.29 }$$

where β is the (fixed) smoothing parameter. The forward model U(σ, z) in (2.28) is based on the finite element (FEM) discretization of the complete electrode model (Somersalo et al 1992). For details of the FEM model, see Kaipio et al (2000), Vauhkonen et al (1998), Vauhkonen (1997). 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 (2.28) 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 σ > 0 is enforced. For more details of the method, see Toivanen et al (2021).

The fixed contact impedances z_* and initial (constant) conductivity estimate for (2.28) are obtained from the solution of the non-linear least squares problem

$$\begin{eqnarray}&&({\sigma }_{0},{z}_{* })=\arg \mathop{\min }\limits_{{\sigma }_{c},z\gt 0}\{\parallel {L}_{e}\left(V-U({\sigma }_{c},z)\right){\parallel }^{2}\},\end{eqnarray} \tag{ 2.30 }$$

where the scalar ${\sigma }_{c}\in {\mathbb{R}}$ is the coefficient of a spatially constant conductivity image σ_c 1 and $z\in {{\mathbb{R}}}^{L}$. The non-linear least squares problem (2.30) is solved by a Gauss Newton optimization.

2.4.2. Linear difference imaging

In linear difference imaging, see e.g. Bagshaw et al (2003), Barber and Seagar (1987), 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}_{* })+{J}_{\sigma }({\sigma }_{1}-{\sigma }_{0})+{e}_{1}\end{eqnarray} \tag{ 2.31 }$$

$$\begin{eqnarray}&&{V}_{2}\approx U({\sigma }_{0},{z}_{* })+{J}_{\sigma }({\sigma }_{2}-{\sigma }_{0})+{e}_{2},\end{eqnarray} \tag{ 2.32 }$$

where the Jacobian matrix J_σ of U(σ, z_*) is evaluated at σ = σ₀. With these linearized models, the difference in the measurements becomes

$$\begin{eqnarray}\begin{array}{l}\delta V=\,{V}_{2}-{V}_{1} & \,\\ =\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)\\ \,={J}_{\sigma }\delta \sigma +\delta e,\end{array}\end{eqnarray} \tag{ 2.33 }$$

where δ σ = σ₂ − σ₁ and δ e = e₂ − e₁. Now, the inverse problem is to reconstruct δ σ based on the difference data δ V and the model (2.33). A widely used formulation for the linear problem is to use the (generalized) Tikhonov regularization with a smoothness promoting regularization functional

$$\begin{eqnarray}&&\hat{\delta \sigma }=\arg \mathop{\min }\limits_{\delta \sigma }\left\{\parallel {L}_{\delta e}(\delta V-{J}_{\sigma }\delta \sigma ){\parallel }^{2}+\parallel {L}_{p}\delta \sigma {\parallel }^{2}\right\},\end{eqnarray} \tag{ 2.34 }$$

where L_{δ e} is the Cholesky factor of the noise precision matrix of δ e so that ${L}_{\delta e}^{{\rm{T}}}{L}_{\delta e}={{\rm{\Gamma }}}_{\delta e}^{-1}={\left({{\rm{\Gamma }}}_{{e}_{1}}+{{\rm{\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}^{{\rm{T}}}{L}_{p}={{\rm{\Gamma }}}_{p}^{-1}$, where the (prior covariance) matrix Γ_p is constructed using the distance based correlation function (Lieberman et al 2010)

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{p}(i,j)=\mathrm{std}{\left(\sigma \right)}^{2}\exp \left(-\displaystyle \frac{\parallel {x}_{i}-{x}_{j}{\parallel }^{2}}{2{a}^{2}}\right),\quad i,j=1,\ldots ,N,\end{eqnarray} \tag{ 2.35 }$$

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 Γ_p (i, j) to the desired covariance for that distance (e.g. 1% of variance).

ACT5 (Rajabi Shishvan et al 2021) 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 (Saulnier et al 2020). 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 (Abdelwahab et al 2020).

A test tank was built from 3/8-inch-thick Plexiglas with interior dimensions of 17.0 cm × 25.5 cm × 17.0 cm. 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 °C–90 °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 (Figure 1). 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. ⁷

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

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 L₁ DN data that requires solving the forward EIT problem with a known conductivity σ ≡ 1 S m⁻¹ and the applied current patterns. Next, the incorrect domain modeling presents as incorrect information about the centers of the electrodes x and surface area for the domain ∣∂Ω∣; the subsequent reconstructions of σ^CAL , ${\sigma }^{{\bf{\exp }}}$, and σ⁰ are recovered on the incorrect domain. For the TV regularized method, the incorrect domain modeling is present in the FEM based forward map U(σ, z) used in the minimization problem (2.28) and computational mesh. For the linear difference imaging method, the incorrect domain modeling is present in the Jacobian matrix of the forward map J_{δ σ} in the minimization problem (2.34) and computational mesh.

Zoom In Zoom Out Reset image size

Accurate domain modeling can be challenging in practice. To address this we explored a domain with a moderately incorrect domain shape, 18cm x 27cm x 19 cm, as well as a domain with a significantly incorrect domain shape, 20 cm × 35 cm × 25 cm ⁸ . 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 2 (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 2.

Zoom In Zoom Out Reset image size

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 section 2.3.1 and section 2.4.2.

2.6.1. Forward modeling and regularization parameters

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 Hamilton et al (2021), measures the distance between the centroids of the reconstructed targets and those of the true targets by

$$\begin{eqnarray}&&\mathrm{LE}=\sqrt{{\left({x}_{1}^{\mathrm{recon}}-{x}_{1}^{\mathrm{truth}}\right)}^{2}+{\left({x}_{2}^{\mathrm{recon}}-{x}_{2}^{\mathrm{truth}}\right)}^{2}+{\left({x}_{3}^{\mathrm{recon}}-{x}_{3}^{\mathrm{truth}}\right)}^{2}}.\end{eqnarray} \tag{ 2.36 }$$

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 (2.36), $({x}_{1}^{\mathrm{truth}},{x}_{2}^{\mathrm{truth}},{x}_{3}^{\mathrm{truth}})$ is our best estimate of the true location of the target center, acknowledging possible errors within millimeters of the truth. Following Hamilton et al (2021), to obtain the centers of the reconstructed targets, we use a thresholded segmentation of the reconstructions to identify targets with MATLAB's regionprops3. 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{eqnarray}&&\mathrm{ScaledLE}=\displaystyle \frac{\sqrt{{\left({x}_{1}^{\mathrm{recon}}-{x}_{1}^{\mathrm{truth}}\right)}^{2}+{\left({x}_{2}^{\mathrm{recon}}-{x}_{2}^{\mathrm{truth}}\right)}^{2}+{\left({x}_{3}^{\mathrm{recon}}-{x}_{3}^{\mathrm{truth}}\right)}^{2}}}{25.5},\end{eqnarray} \tag{ 2.37 }$$

where we have divided the LE from (2.36) by 25.5 cm, 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.