Enhancing D-bar reconstructions for electrical impedance tomography with conformal maps

TL;DR

This paper introduces methods to improve 2D electrical impedance tomography reconstructions by applying conformal maps, enabling more efficient and focused imaging in specific regions of interest.

Contribution

It proposes novel uses of conformal maps, including domain transformations and magnification, combined with specialized current patterns for enhanced D-bar reconstructions.

Findings

Transforming domains to the unit disk improves computational efficiency.

Magnification of regions of interest enhances local reconstruction accuracy.

Numerical tests demonstrate improved imaging in simulated data scenarios.

Abstract

We present a few ways of using conformal maps in the reconstruction of two-dimensional conductivities in electrical impedance tomography. First, by utilizing the Riemann mapping theorem, we can transform any simply connected domain of interest to the unit disk where the D-bar method can be implemented most efficiently. In particular, this applies to the open upper half-plane. Second, in the unit disk we may choose a region of interest that is magnified using a suitable M\"obius transform. To facilitate the efficient use of conformal maps, we introduce input current patterns that are named conformally transformed truncated Fourier basis; in practice, their use corresponds to positioning the available electrodes close to the region of interest. These ideas are numerically tested using simulated continuum data in bounded domains and simulated point electrode data in the half-plane. The connections to practical electrode measurements are also discussed.

Equations

$$\sigma \in L_{+}^{\infty }(\Omega ):=\{\kappa \in L^{\infty }(\Omega )|\mathrm{ess}\inf \kappa >0\}$$

$$\displaystyle\left\{\begin{array}[]{ll}\nabla\cdot(\sigma\nabla u)=0&\rm{in}\ \Omega,\\[2.84526pt] u=g&\rm{on}\ \partial\Omega.\end{array}\right.$$

$$\Lambda_{\sigma}g:=\sigma\frac{\partial u}{\partial\nu}\big{|}_{\partial\Omega}\in H^{-1/2}_{\diamond}(\partial\Omega),$$

$$H_{\diamond }^s(\partial \Omega ):=\{v\in H_{\diamond }^s(\partial \Omega )|⟨f,1{⟩}_{\partial \Omega }=0\},s\in \mathbb{R},$$

$${\Lambda }_{\sigma }:H^{1/2}(\partial \Omega )/\mathbb{C}\to H_{\diamond }^{-1/2}(\partial \Omega ),$$

$$\displaystyle\left\{\begin{array}[]{ll}\nabla\cdot(\sigma\nabla u)=0&\rm{in}\ \Omega,\\[2.84526pt] {\displaystyle\sigma\frac{\partial u}{\partial\nu}}=f&\rm{on}\ \partial\Omega\end{array}\right.$$

$${\mathcal{R}}_{\sigma }f=u{|}_{\partial \Omega }\in H^{1/2}(\partial \Omega )/\mathbb{C},$$

$${\Lambda }_{\sigma }{\mathcal{R}}_{\sigma }=I\mathrm{and}{\mathcal{R}}_{\sigma }{\Lambda }_{\sigma }=I,$$

$$(-\Delta +q)v=Ev\text{in }{\mathbb{R}}^2,q=\frac{\Delta \sqrt{\sigma }}{\sqrt{\sigma }},E=0,$$

$$\exp (-\mathrm{i}kz)\psi (z,k)-1\in W^{1,\tilde{p}}({\mathbb{R}}^2),1/\tilde{p}=1/p-1/2.$$

$$\mathbf{t}(k)={\int }_{\Omega }{e}^{\mathrm{i}\overline{z}\overline{k}}q(z)\psi (z,k)\mathrm{d}z.$$

$${\partial }_{\overline{k}}\mu (z,k)=\frac{\mathbf{t}(k)e_{-k}(z)}{4\pi \overline{k}}\overline{\mu (z,k)},e_{-k}(z)=\exp (-\mathrm{i}(kz+\overline{k}\overline{z})),k\ne 0,$$

$$\sigma (z)=\mu (z,0{)}^2.$$

$$(S_k\varphi )(z):={\int }_{\partial \Omega }{G}_k(z-y)\varphi (y)\mathrm{d}s(y),$$

$$G_k(z):=e^{\mathrm{i}kz}{g}_k(z),g_k(z):=\frac{1}{(2\pi {)}^2}{\int }_{{\mathbb{R}}^2}\frac{e^{\mathrm{i}z\cdot y}}{|y{|}^2+2k(y_1+\mathrm{i}{y}_2)}\mathrm{d}y.$$

$$\displaystyle\big{(}I+S_{k}(\Lambda_{\sigma}-\Lambda_{1})\big{)}\psi(\,\cdot\,,k)|_{\partial\Omega}=e^{{\rm i}k\,\cdot},$$

$$G_{k}=(G_{0}+\hat{H}_{k})+\Big{(}H_{1}(0)-\frac{1}{2\pi}\log\left|{k}\right|\Big{)}.$$

$$S_k=\frac{1}{2}{\mathcal{R}}_1+{\hat{\mathcal{H}}}_k\mathrm{on}{H}_{\diamond }^{-1/2}(\partial \Omega ),$$

$$I+S_k({\Lambda }_{\sigma }-{\Lambda }_1)$$

$$=\frac{1}{2}(I+{\mathcal{R}}_1{\Lambda }_{\sigma })+{\hat{\mathcal{H}}}_k({\Lambda }_{\sigma }-{\Lambda }_1).$$

$$\mathbf{t}(k)={\int }_{\partial \Omega }{e}^{\mathrm{i}\bar{k}\bar{z}}({\Lambda }_{\sigma }-{\Lambda }_1)\psi (z,k)\mathrm{d}s(z),$$

$$\mu (z,k)=1+\frac{1}{(2\pi {)}^2}{\int }_{{\mathbb{R}}^2}\frac{\mathbf{t}(k^{\prime })}{(k-k^{\prime }){\bar{k}}^{\prime }}{e}^{\mathrm{i}(k^{\prime }z+\overline{k^{\prime }}\overline{z})}\overline{\mu (z,k^{\prime })}\mathrm{d}{k}_1^{\prime }\mathrm{d}{k}_2^{\prime }.$$

$$\partial \Omega =\{z\in \mathbb{C}|z=\gamma (s),s\in \mathbb{R}\}.$$

$$\phi_{n}(s)=\frac{1}{\sqrt{\left|{\partial\Omega}\right|}}\exp\Big{(}\mathrm{i}\frac{2\pi n}{\left|{\partial\Omega}\right|}s\Big{)},\qquad n=-N,...,N,$$

$$A:H^s(\partial \Omega )\to H^r(\partial \Omega ),s,r\in (-1,1),$$

$${\mathbf{A}}_{m,n}=⟨A{\varphi }_n,{\overline{\varphi }}_m{⟩}_{\partial \Omega },m,n=-N,\dots ,N,$$

$$⟨Af,{\overline{\varphi }}_m⟩\approx [\mathbf{A}\mathbf{f}{]}_m,m=-N,\dots ,N.$$

$$\Big{(}\frac{1}{2}({\mathbf{I}}+{\mathbf{R}}_{1}{\mathbf{L}}_{\sigma})+\hat{{\mathbf{H}}}_{k}({\mathbf{L}}_{\sigma}-{\mathbf{L}}_{1})\Big{)}{\mathbf{p}}_{k}={\mathbf{e}}_{k}$$

$${\mathbf{t}}(k)=\int_{\partial\Omega}e^{\mathrm{i}\bar{k}\bar{z}}\mathcal{F}^{-1}\big{(}({\mathbf{L}}_{\sigma}-{\mathbf{L}}_{1}){\mathbf{p}}_{k}\big{)}(z)\,{\rm d}s(z),$$

$${\mathbf{t}}_R(k)=\{\begin{array}{ll}\mathbf{t}(k)& \text{mbox}for|k|<R,\\ 0& \text{otherwise},\end{array}$$

$${\mathbf{t}}_{R,c}(k)=\{\begin{array}{ll}{\mathbf{t}}_c(k)& \text{mbox}for|k|<R,\\ 0& \text{otherwise},\end{array}$$

Full text

The work of NH and JT was supported by the Academy of Finland (decision 267789). The work of LP was supported by Estonian government grant PUT1093.

1. Introduction

Electrical impedance tomography (EIT) is an intensively studied research topic because it is a fruitful mathematical problem and has also applications in medical imaging and non-destructive testing. See [5, 9, 48] for overviews of EIT. To summarize, the goal of EIT is to form an image of the conductivity distribution inside a physical body by using electric measurements on its boundary. This corresponds to a nonlinear, ill-posed inverse problem, which needs to be tackled by some regularized algorithm or by resorting to Bayesian inference. In this work, we use the D-bar method, which is an example of a direct reconstruction algorithm that assumes the knowledge of the idealized boundary measurements, i.e., the Dirichlet-to-Neumann (D-to-N) map. To be more precise, the D-to-N map is linked to the unknown conductivity by equations found by investigating the inverse problem analytically.

The D-to-N map is a theoretical object representing all possible potential-to-current measurements on the boundary of the domain of interest. Alternatively, one may consider the Neumann-to-Dirichlet (N-to-D) map, broadly speaking the inverse of D-to-N map, which sends the applied boundary current to the measured boundary potential. The N-to-D map is arguably the better way to model real measurements: It is smoothening and thus easier to handle numerically, and the approximative connection between real-world electrode measurements and the N-to-D map has also been analyzed [22, 18]. In particular, the boundary condition between the employed electrodes is in practice always of the homogeneous Neumann type, which obviously favors the use of the N-to-D map (cf. [22]). Due to these reasons, we simulate all boundary measurements employed in this paper by computing (suitably truncated versions of) the N-to-D maps for the considered target configurations, and subsequently form the needed (relative) D-to-N maps only when applying the D-bar method to reconstructing conductivities in the unit disk. Such an approach is enabled in two spatial dimensions by the use of conformal mappings as explained in what follows.

In this paper, we consider EIT in two-dimensional simply connected domains; observe that such considerations do have practical relevance since a two-dimensional model can be utilized if the imaged three-dimensional object is cylindrically symmetric with homogeneous Neumann conditions at the top and bottom boundaries. The original imaging problem is transformed to the unit disk with the help of a suitable conformal map — by doing so, we enable the use of our direct reconstruction method in the simplest possible setting. Furthermore, in the unit disk we may use Möbius transforms to ‘magnify’ a region of interest (ROI). Both of these steps are implemented by modifying the applied boundary current densities so that they account for the stretching introduced by the conformal maps and correspond to the standard (truncated) Fourier boundary current basis for the unit disk where the D-bar algorithm is implemented. We call the transformed domain and conductivity the virtual domain and virtual conductivity, respectively, as opposed to the true domain and true conductivity that correspond to the actual measurements. As demonstrated by our numerical experiments, forming the reconstruction in the virtual domain results in an improved outcome inside the ROI.

Magnifying the ROI by a Möbius transform is motivated by practical measurement set-ups where the number of available electrodes is restricted, but their positions can be freely chosen: If the electrodes are placed equiangledly on the boundary of the unit disk, their number, loosely speaking, restricts the number of Fourier modes that can be used as input current densities (see [22, 18] and Appendix A). Similarly, if one wants to get more accurate information about a certain ROI inside the imaged object, the pattern of electrodes should arguably be the densest on the boundary section closest to the ROI, enabling the use of current patterns with higher spatial frequencies there. The magnification by a Möbius transform can thus be thought of as mapping a nonuniform electrode pattern designed for accurate imaging of the ROI onto an equiangled pattern that facilitates a straightforward reconstruction by the D-bar method in the unit disk. We give further evidence for this interpretation by presenting theoretical results on the possibility of approximating a given (continuum) current density by electrode measurements modeled by the point electrode model (PEM) [18]. This would also straightforwardly lead to similar results for the complete electrode model (CEM) [10] with small electrodes.

On a general level, the idea of using conformal mappings in EIT-related considerations is definitely not new: In [39] the Möbius transforms between disks and the half-plane were already presented in the discussion of visibilities and sensitivities of computed tomography. This discussion continued as a question of distinguishability for EIT in [24, 8], using again conformal mappings in [50, 14]. Following the ideas presented in [50], it was recently demonstrated in [23] that the CEM is approximately conformally invariant (with constant contact resistances). Möbius transforms have also been employed in the implementation of the so-called convex source support algorithm for detecting anomalies from a single measurement pair of EIT both in a bounded domain [19] and in a half-plane [20]. This paper uses the same transformation properties of conformal maps as the aforementioned articles, but we also introduce and numerically test the concept of virtual domain with the help of the conformally transformed truncated Fourier basis introduced in [14]. In particular, to the authors’ best knowledge, this work presents the first D-bar reconstructions in a half-space geometry.

Implementing a numerical method for computing conformal maps between planar domains is outside the scope of this manuscript. Therefore, we employ the Schwarz–Christoffel Toolbox [12] in MATLAB to compute numerical conformal maps and their derivatives between the unit disk and arbitrary polygons [47]. On the other hand, Möbius transforms and their derivatives are trivial to implement numerically.

This text is organized as follows. In Section 2, we introduce our mathematical framework. The main ideas of the D-bar method are recalled in Section 3, with special emphasis on its numerical implementation (in the unit disk). Section 4 briefly considers how a transformation induced by a conformal mapping affects the Neumann boundary value problem defining the N-to-D map. The numerical results are presented in Section 5. The treatment of half-space geometry as well as the theoretical results related to the approximation of (relative) N-to-D maps with the help of pointlike electrodes are included in two appendices.

2. The setting

The inverse problem of EIT is mathematically described by the Calderón problem [7]. Let $\Omega\subset{\mathbb{R}}^{2}$ be a simply connected, bounded domain with a Lipschitz boundary and let

[TABLE]

be a strictly positive, real-valued, isotropic conductivity. A boundary voltage distribution $g\in H^{1/2}(\partial\Omega)$ creates an electrostatic interior potential $u\in H^{1}(\Omega)$ that solves the Dirichlet problem

[TABLE]

The resulting current density through the boundary is

[TABLE]

where $\nu\in L^{\infty}(\partial\Omega,{\mathbb{R}}^{2})$ is the outward unit normal of $\partial\Omega$, the linear operator $\Lambda_{\sigma}$ is the D-to-N map and

[TABLE]

is the subspace of $H^{s}(\partial\Omega)$ consisting of those functions/distributions that have vanishing mean. Here and in what follows, $\langle\cdot,\cdot\rangle_{\partial\Omega}:H^{s}(\partial\Omega)\times H^{-s}(\partial\Omega)\to{\mathbb{C}}$ denotes the (bilinear) dual bracket on $\partial\Omega$. Since the boundary potentials $g$ and $g+c$, $c\in{\mathbb{C}}$, for (3) correspond to the interior potentials $u$ and $u+c$, respectively, we may define the D-to-N map as

[TABLE]

where the use of a quotient space reflects the irrelevance of the ground level of potential in modeling the underlying physical phenomenon. It follows from the standard theory of elliptic boundary value problems that $\Lambda_{\sigma}$ is bounded. The Calderón problem is to determine $\sigma$ from the knowledge of $\Lambda_{\sigma}$.

As explained above, in this work we assume to be able to measure (an approximation of) the N-to-D map ${\mathcal{R}}_{\sigma}:H_{\diamond}^{-1/2}(\partial\Omega)\to H^{1/2}(\partial\Omega)/{\mathbb{C}}$ instead of the D-to-N map. If the current density on $\partial\Omega$ is set to $f\in H_{\diamond}^{-1/2}(\partial\Omega)$, then the electrostatic potential $u\in H^{1}(\Omega)/{\mathbb{C}}$ uniquely solves the Neumann problem

[TABLE]

interpreted in the weak sense. Notice that (7) does not fix the ground level of potential and thus the solution is unique only up to an additive constant, i.e., up to the ground level of potential. We then define

[TABLE]

where $u$ is the solution to (7). It is easy to see that ${\mathcal{R}}_{\sigma}$ is the inverse of $\Lambda_{\sigma}$, that is,

[TABLE]

where the identity map $I$ is that of $H_{\diamond}^{-1/2}(\partial\Omega)$ in the first identity and that of $H^{1/2}(\partial\Omega)/{\mathbb{C}}$ in the second one. In particular, ${\mathcal{R}}_{\sigma}$ is bounded. Take note that the boundary quotient space $H^{1/2}(\partial\Omega)/{\mathbb{C}}$ can be replaced by its dual $H^{1/2}_{\diamond}(\partial\Omega)$ at all occurrences in the above definitions if the ground level of potential is systematically chosen so that all electrostatic potentials have vanishing mean over $\partial\Omega$.

The direct reconstruction approach is to find equations that link $\Lambda_{\sigma}$ to $\sigma$ as well as to implement such a connection as a practical algorithm. The D-bar method originates from the work of Ablowitz, Nachman, Beals and Coifman [3, 1] in a related inverse scattering problem. It is based on a non-linear Fourier transform, where exponentially behaving complex geometric optics (CGO) solutions of Faddeev [13] are used. A boundary integral equation deduced by R. G. Novikov [36] and Nachman [34] links $\Lambda_{\sigma}$ to the CGO solutions. The two-dimensional Calderón problem was solved, i.e., its uniqueness was deduced, using the D-bar method by Nachman [35] for conductivities in $W^{2,p}(\Omega)$, $p>1$. The solution technique was later generalized by Brown and Uhlmann [6] for less regular conductivities. Using D-bar techniques and quasiconformal maps, the two-dimensional Calderón problem was finally completely solved by Astala and Päivärinta [2] for $\sigma\in L^{\infty}_{+}(\Omega)$.

The ideas of exploiting conformal transformations presented in this manuscript are arguably as useful for any of the aforementioned variants of the D-bar method and, more generally, for any reconstruction method where the D-to-N or N-to-D map is approximated by a matrix in a truncated Fourier basis. However, we will only numerically test these ideas using Nachman’s original method.

3. The D-bar method for EIT

The origin of the D-bar method lies with the direct solution to a certain inverse scattering problem formulated as the Gel’fand–Calderón (GC) problem [15, 7]: By writing $u=\sigma^{1/2}v$, (3) is transformed into a zero-energy Schrödinger equation

[TABLE]

which necessitates the smoothness assumption $\sigma\in W^{2,p}(\Omega)$. Similarly, the equations governing diffuse optical tomography (DOT) and acoustic tomography (AT) can be recast in the form (10) but with negative and positive energies [11, 46], respectively. We remark what the method of Astala and Päivärinta [2] does not require transforming the conductivity equation into the form (10), and in theory it thus works with more general conductivities $\sigma\in L^{\infty}_{+}(\Omega)$. However, a number of numerical results indicate that Nachman’s method remains functional even for discontinuous conductivities or if the smoothness assumptions are violated in some other way [28, 27, 43].

Complex Geometrical Optics (CGO) solutions $\psi$ are exponentially behaving solutions to the Schrödinger equation (10). They depend on the spatial variable $z\in{\mathbb{C}}$ as well as on a spectral variable $\zeta\in{\mathbb{C}}^{2}$ which satisfies $\zeta\cdot\zeta=E=0$. Such $\zeta$ can be parametrized by a single complex parameter $k\in{\mathbb{C}}$. The exponential behavior is controlled by

[TABLE]

Nachman dubbed the potential $q$ appearing in (10) as of conductivity type and proved that for such potentials the CGO solution $\psi(z,k)$ exists for all $k\in{\mathbb{C}}$. Thus for any $\sigma\in W^{2,p}(\Omega)$ and $k\in{\mathbb{C}}$, one may define the scattering transform, i.e., a certain nonlinear Fourier transform, as

[TABLE]

It can be shown that the function $\mu(z,k)=\psi(z,k)\exp(-\mathrm{i}kz)$ satisfies the following D-bar equation:

[TABLE]

which is uniquely solvable in a certain weighted Sobolev space [35]. The conductivity can then be recovered via

[TABLE]

What remains to be explained is how the CGO solutions as well as the equations (12) and (13) are connected to the boundary measurement $\Lambda_{\sigma}$.

For simplicity, we assume in the following that $\sigma\equiv 1$ in some interior neighborhood of $\partial\Omega$. Let us define a single-layer operator

[TABLE]

where $G_{k}$ is the Faddeev’s Green function

[TABLE]

For all $k\in{\mathbb{C}}\setminus\{0\}$, the trace of the CGO solution $\psi(\,\cdot\,,k)|_{\partial\Omega}\in H^{1/2}(\partial\Omega)$ can be solved from the boundary integral equation

[TABLE]

where $\Lambda_{1}$ is the D-to-N map corresponding to the homogeneous unit conductivity.

When solving (17), the singularity in the kernel of the operator $S_{k}$ can be handled as follows; see [44] for more details. We introduce the harmonic function $H_{k}=G_{k}-G_{0}$, where $G_{0}(z)=-1/(2\pi)\log\left|{z}\right|$ is the fundamental solution for the Laplacian, and note that $H_{k}(z)=H_{1}(kz)-1/(2\pi)\log\left|{k}\right|$. This motivates defining yet another auxiliary function $\hat{H}_{k}(z)=H_{1}(kz)-H_{1}(0)$, for which $\hat{H}_{k}(0)=0$. It follows that the kernel of $S_{k}$ can be restructured as

[TABLE]

Since the latter term is independent of the spatial variable and it is known that $2S_{0}={\mathcal{R}}_{1}$, it is easy see that

[TABLE]

with $\hat{\mathcal{H}}_{k}$ being the integral operator defined by the well-behaving kernel $\hat{H}_{k}$. In particular,

[TABLE]

The operator $\hat{\mathcal{H}}_{k}$ includes the Faddeev’s Green function $G_{k}$ as a part of its kernel.

The scattering transform (12) also admits an alternative form

[TABLE]

which can be evaluated after solving $\psi(z,k)|_{\partial\Omega}$ from (17). Subsequently, the D-bar equation (13) can be written as a Lippman–Schwinger type integral equation, uniquely solvable in $k\in{\mathbb{C}}$ for each fixed $z\in\Omega$:

[TABLE]

Finally, (14) provides the solution to the two-dimensional Calderón problem.

3.1. Standard numerical implementation of the D-bar method

The D-bar method has a standard numerical implementation; see [42, 33, 30] and [25, 26] for its applications to simulated and experimental data, respectively.

Let the Lipschitz continuous, $|\partial\Omega|$-periodic function $\gamma:{\mathbb{R}}\to{\mathbb{C}}$ define a parametrization of $\partial\Omega$ with respect to its arclength, that is, $|\gamma^{\prime}|\equiv 1$ almost everywhere and

[TABLE]

The corresponding $L^{2}(\partial\Omega)$-orthonormal truncated Fourier basis is defined via

[TABLE]

as functions of the arclength parameter. The linear operators ${\mathcal{R}}_{\sigma},{\mathcal{R}}_{1},\Lambda_{\sigma},\Lambda_{1}$ and $\hat{\mathcal{H}}_{k}$, needed in (19) and (20), are approximated by matrices ${\mathbf{R}}_{\sigma},{\mathbf{R}}_{1},{\mathbf{L}}_{\sigma},{\mathbf{L}}_{1}$ and $\hat{{\mathbf{H}}}_{k}$, respectively, with respect to the basis (23). To be more precise, a bounded linear operator (see, e.g., [16, p. 20])

[TABLE]

is replaced in the numerical algorithm by the matrix ${\mathbf{A}}\in{\mathbb{C}}^{(2N+1)\times(2N+1)}$ defined through

[TABLE]

where $m,n=-N$ refers to the first row/column of the matrix. If ${\mathbf{f}}_{n}=\langle f,\overline{\phi}_{n}\rangle_{\partial\Omega}$, $n=-N,\dots,N$, are the Fourier coefficients of an input $f\in H^{s}(\partial\Omega)$, then the Fourier coefficients of $Af$ can be approximated as

[TABLE]

It is straightforward to deduce that the appropriate interpretation of ${\mathbf{A}}$ as a finite-dimensional mapping from $H^{s}(\partial\Omega)$ to $H^{r}(\partial\Omega)$ converges to $A:H^{s}(\partial\Omega)\to H^{r}(\partial\Omega)$ pointwise as $N$ goes to infinity (see, e.g., [38]). Furthermore, the convergence occurs in the operator norm if $A:H^{s}(\partial\Omega)\to H^{r}(\partial\Omega)$ is compact. The case of mean-free spaces, i.e., $A:H^{s}(\partial\Omega)/{\mathbb{C}}\simeq H^{s}_{\diamond}(\partial\Omega)\to H^{r}_{\diamond}(\partial\Omega)\simeq H^{r}(\partial\Omega)/{\mathbb{C}}$ can be treated by simply dropping out the constant basis function corresponding to $n=0$.

To simulate data in a bounded domain,111See Appendix B for the case of the upper half-plane. we compute the necessary direct solutions $u$ of (7) in order to construct the N-to-D matrix ${\mathbf{R}}_{\sigma}$ by a finite element method (FEM) with piecewise linear basis functions. The corresponding D-to-N matrix ${\mathbf{L}}_{\sigma}$ can be obtained by inverting ${\mathbf{R}}_{\sigma}$, though this process becomes unstable for large $N$ since ${\mathcal{R}}_{\sigma}$ is compact and $\Lambda_{\sigma}$ unbounded on $L^{2}_{\diamond}(\partial\Omega)$. To compute the matrix $\hat{{\mathbf{H}}}_{k}$ approximating the operator $\hat{\mathcal{H}}_{k}$ in (19) we need an algorithm for evaluating $G_{k}(z)$ of (16). In the zero energy case this is trivial, as we can use the property $g_{k}(z)=g_{1}(kz)$ and [4, (3.10)] to compute $g_{1}$ using MATLAB’s built-in exponential-integral function expint.

Solving (17) could be skipped by replacing $\psi(z,k)$ with $\exp(\mathrm{i}kz)$ in (20) without considerably affecting the reconstructions (see, e.g., [41]) because the two functions have the same asymptotic behavior as $\left|{z}\right|$ tends to infinity. However, here we solve the full problem by approximating (17) as a matrix equation (cf. (19))

[TABLE]

in the truncated Fourier basis. Subsequently, an approximate scattering transform is computed as

[TABLE]

where $\mathcal{F}^{-1}$ corresponds to evaluating the function with the Fourier coefficients $({\mathbf{L}}_{\sigma}-{\mathbf{L}}_{1}){\mathbf{p}}_{k}$. The integral in (26) is approximated by the trapezoidal rule with respect to a dense enough discretization of the arclength parameter $s\in[0,\left|{\partial\Omega}\right|)$.

In practice, the discretized boundary integral equation (25) is not (stably) solvable for all $k\in{\mathbb{C}}$; typically, problems arise when $|k|$ is large and/or when the data contain noise. However, the D-bar method can be regularized by truncating the scattering transform, that is, introducing

[TABLE]

with $R>0$ suitably defined as a decreasing function of the noise level [29]. In this work, we resort to a somewhat more general truncation

[TABLE]

where ${\mathbf{t}}_{c}(k)=0$ when $\left|{{\mathbf{t}}(k)}\right|>c$ with $c>0$ being a suitable cut-off parameter. We have $\left|{{\mathbf{t}}_{R,c}(k)}\right|\leq\left|{{\mathbf{t}}_{R}(k)}\right|$ for every $k\in{\mathbb{C}}$ so the convergence results of [29] still hold. However, the analysis on why (28) would be a better truncation than (27) is nontrivial and absent to the authors’ best knowledge. In our numerical experiments, the values of $R,c>0$ are chosen after visually investigating the right-hand side of (26) on an extensive grid over $k$ in the complex plane.

The integral equation (21) can be tackled, e.g., by a modified Lippman–Schwinger solver originating from Vainikko’s work [49]: The equation is periodized and solved using GMRES as explained in [30]. The resulting solution is (an approximation of) the function $\mu(z,\,\cdot)$ at the investigated reconstruction point $z\in\Omega$ and a grid of values over the second variable, i.e., $k$. We then set $k=0$ and get the reconstructed conductivity at $z$ via (14).

3.2. The case of unit disk

If the D-bar method is implemented in a domain that is not a disk, the reference N-to-D map $\mathcal{R}_{1}$ (or $\Lambda_{1}$) must also be approximated numerically by, e.g., FEM; cf. (25) and (26). One benefit in transforming the reconstruction problem to the unit disk is that there these auxiliary maps admit trivial diagonal representations, which makes the implementation of the D-bar algorithm more efficient.

Indeed, the reference D-to-N map $\tilde{\Lambda}_{1}$ corresponding to the unit disk $D:=\{z\in{\mathbb{C}}:\left|{z}\right|<1\}$ with homogeneous unit conductivity can be characterized in terms of the Fourier basis (23) as

[TABLE]

Here and in what follows, we denote a Fourier basis function by $\tilde{\phi}_{n}$ if it corresponds to $\Omega=D$, that is, $\tilde{\phi}_{n}(\theta)=1/\sqrt{2\pi}\exp(\mathrm{i}n\theta)$, with $\theta$ being the polar angle. Such ‘tilde notation’ is systematically adopted to indicate which entities are related to $D$. As a matrix with respect to truncated Fourier basis, the relations (29) takes the form

[TABLE]

where the constant basis function has been dropped. This obviously also induces a diagonal representation for the reference N-to-D matrix $\tilde{{\mathbf{R}}}_{1}$.

The following lemma lists other symmetries that streamline the implementation of the D-bar method in the unit disk by reducing the computations required for forming the matrix $\hat{{\mathbf{H}}}_{k}$ appearing in (25). These results are formulated for the matrix ${\mathbf{S}}_{k}$ approximating the single-layer operator (15), but they also hold for $\hat{{\mathbf{H}}}_{k}$ due to the relation

[TABLE]

induced by (18). The proof is omitted since the claims are easy consequences of the symmetries exhibited by the Faddeev Green’s function $G_{k}$.

Lemma 3.1**.**

For all $k\in{\mathbb{C}}\setminus\{0\}$,

[TABLE]

If $\Omega$ is a disk, then also

[TABLE]

where $\odot$ denotes elementwise multiplication of matrices, $\alpha$ is the polar angle of $k\in{\mathbb{C}}$ and $({\mathbf{E}}_{\alpha})_{n,m}=\exp(\mathrm{i}\alpha(m-n))$.

4. Conformal transformation of the domain

Let $D=\{z\in{\mathbb{C}}:\left|{z}\right|<1\}$ still be the open unit disk. According to the Riemann mapping theorem, there is a biholomorphic/conformal map $\Phi:\Omega\to D$. Since $\partial\Omega$ is Lipschitz, $\Phi$ extends to a homeomorphism $\Phi:\overline{\Omega}\to\overline{D}$ [37]. The inverse of $\Phi$ is denoted by $\Psi$. We denote by $J_{\Phi}$ the Jacobian of $\Phi$ (interpreted as a map from ${\mathbb{R}}^{2}$ to itself) and by $|\Phi^{\prime}|=\sqrt{\det J_{\Phi}}$ the absolute value of its (complex) derivative. The following lemma, the basic principle of which is well-known, relates the solution of (7) to that of

[TABLE]

for a certain $\tilde{f}\in H^{-1/2}_{\diamond}(\partial D)$ and $\tilde{\sigma}:=\sigma\circ\Psi$.

Lemma 4.1**.**

Let $u\in H^{1}(\Omega)/{\mathbb{C}}$ be the solution of (7) for some $f\in H^{-1/2}_{\diamond}(\partial\Omega)$. Then, $\tilde{u}:=u\circ\Psi\in H^{1}(D)/{\mathbb{C}}$ is the solution of (33) for $\tilde{f}\in H^{-1/2}_{\diamond}(\partial D)$ defined by

[TABLE]

for all $\tilde{g}\in H^{1/2}(D)/{\mathbb{C}}$.

Proof.

We mimic the proof of [40, Lemma A.3] and start by showing the linear map $L_{\Phi}:\tilde{v}\mapsto\tilde{v}\circ\Phi$ is bounded from $H^{1}(D)/{\mathbb{C}}$ to $H^{1}(\Omega)/{\mathbb{C}}$. Let us denote by $D_{\epsilon}$ the open, origin-centered disk of radius $1-\epsilon$ and by $\chi_{D_{\epsilon}}$ its characteristic function. For any $\tilde{v}\in H^{1}(D)/{\mathbb{C}}$, it holds that

[TABLE]

where we used $J_{\Phi}^{\rm-1}J_{\Phi}^{\rm-T}=(1/\left|{\Phi^{\prime}}\right|^{2})I$, a consequence of the Cauchy–Riemann equations. Since the integrands of the first and the last term of (35) are clearly increasing as $\epsilon\to 0$, the Lebesgue monotone convergence theorem gives

[TABLE]

which shows the boundedness of $L_{\Phi}:H^{1}(D)/{\mathbb{C}}\to H^{1}(\Omega)/{\mathbb{C}}$ by virtue of the Poincaré inequality. (Note that $L_{\Phi}$ is actually a homeomorphism as its inverse is obviously defined by $v\mapsto v\circ\Psi$.)

As the trace map is continuous from $H^{1}(\Omega)/{\mathbb{C}}$ to $H^{1/2}(\partial\Omega)/{\mathbb{C}}$ [21] and $\Phi|_{\partial\Omega}:\partial\Omega\to\partial D$ is a homeomorphism [37], it follows, e.g., by the density of the embedding $C^{\infty}(\overline{D})\hookrightarrow H^{1}(D)$ [16], that any $\tilde{v}\in H^{1}(D)/{\mathbb{C}}$ satisfies

[TABLE]

modulo an additive constant almost everywhere on $\partial\Omega$. Let $\tilde{\gamma}^{-1}:H^{1/2}(\partial D)/{\mathbb{C}}\to H^{1}(D)/{\mathbb{C}}$ be a bounded right-inverse for the (quotient) trace operator in $D$ (cf. [16]). By (36), the boundedness of the quotient trace map and the continuity of $L_{\Phi}$, an arbitrary $\tilde{g}\in H^{1/2}(\partial D)/{\mathbb{C}}$ satisfies

[TABLE]

In particular,

[TABLE]

where the fact $\langle\tilde{f},1\rangle_{\partial D}=0$ was implicitly used. Together with the obvious linearity of $\tilde{f}$, this demonstrates that $\tilde{f}\in H^{-1/2}_{\diamond}(\partial D)$.

A similar computation as (35), Hölder’s inequality and the Lebesgue dominated convergence theorem can be used to show that

[TABLE]

where the second equality corresponds to the variational formulation of (7) with the test function $\tilde{v}\circ\Phi\in H^{1}(\Omega)/{\mathbb{C}}$ [16]. It thus follows from (36) and the definition of $\tilde{f}$ that actually

[TABLE]

which is the variational formulation of (33) and completes the proof. ∎

Remark 4.1**.**

If $f\in L^{2}_{\diamond}(\partial\Omega)$ and the boundary $\partial\Omega$ is regular enough, one can make the interpretation

[TABLE]

in (34). This holds, e.g., if the Lipschitz boundary $\partial\Omega$ consists of a finite number of $C^{1,\alpha}$ smooth arcs, with $0<\alpha\leq 1$, because then $\Psi$ has a continuously differentiable extension onto the open $C^{1,\alpha}$ arcs [37], allowing a change of variables one arc at a time. In other words,

[TABLE]

in the sense of $L^{2}(\partial D)$ and $L^{2}(\partial\Omega)$, respectively.

As we want to use the truncated Fourier current basis $\tilde{f}=\tilde{\phi}_{n}$, $1\leq|n|\leq N$, on the boundary of the virtual reconstruction domain $D$, the transformations (38) indicates that one should drive the (true) current patterns

[TABLE]

through $\partial\Omega$. These are called conformally transformed Fourier basis functions and were mentioned in [14]. See Figure 1 for an example of the boundary current transformation (39) in a rectangular $\Omega$.

In the following, $\Phi$ is often a composition of two conformal maps: one being a Schwarz–Christoffel map produced by the SC-Toolbox [12] sending a given polygon onto the unit disk, the other being a Möbius transform magnifying a certain ROI within the unit disk. (Take note that such a composition is also itself a Schwarz–Christoffel map.)

4.1. Möbius transforms

All conformal maps of $D$ onto itself are given by the Möbius transforms of the form

[TABLE]

where $\alpha\in{\mathbb{R}}$ determines a rotation and $a\in D$ is a complex number that is mapped to the origin. In the following, we systematically choose $\alpha=0$. Such Möbius transforms can be used to magnify a certain part of the unit disk as illustrated by Figure 2 for $a=0.6$. With such a choice, the inhomogeneities inside the ROI close to the right edge of $D$ are magnified.

Our leading idea is that the virtual conductivity on the right in Figure 2 is easier to reconstruct than the true one on the left if employing a finite number of standard Fourier basis functions, i.e., the Fourier basis with respect to the polar angle $\{\tilde{\phi}\}_{1\leq|n|\leq N}$, as the input current densities on the respective boundaries. In other words, if it is a priori known that the ROI lies close to the right-hand edge of the true domain, it is beneficial to employ the conformally transformed Fourier basis functions (39), with $\Phi={\mathcal{M}}_{a}$, as the true input current densities for the true domain on the left, which corresponds to the use of the standard Fourier current basis for the virtual domain on the right. The standard, efficient numerical implementation of the D-bar algorithm in the unit disk described in Section 3.2 can then be applied to the virtual data on the boundary of the virtual domain (cf. Lemma 4.1), and the obtained reconstruction can be subsequently mapped by $\Psi$ back to the true domain. Our numerical experiments in Section 5 demonstrate this procedure results in a more accurate reconstruction in the ROI. (Take note that the true domain does not need to be the unit disk nor the conformal map a Möbius transform when magnifying a ROI.)

In Appendix A, it is shown that relative continuum measurements for the true domain with the conformally transformed Fourier basis currents (39) can be approximated accurately by using a finite number of pointlike electrodes (cf. [18]) if their locations are chosen to be the preimages under $\Phi$ of a uniform grid on the boundary of the virtual domain. In other words, the use of a conformally transformed Fourier basis corresponds in practice to employing a certain nonuniform electrode grid that is densest close to the ROI.

5. Numerical tests

Our numerical experiments are divided into three sections. In the first part, we investigate how reconstructions of a certain conductivity phantom in polygonal domains depend on whether the D-bar method is implemented directly in the true domain by resorting to the standard Fourier current basis (23) or in virtual domain (the unit disk) by using the conformally transformed Fourier current basis (39) on the boundary of the true domain. These tests do not consider the magnification of a ROI, but only aim at demonstrating that both approaches produce approximately as good reconstructions. This suggests that the choice between the two options should be made by comparing their computational expense.

In the second part, we numerically demonstrate that it is possible to produce more accurate reconstructions of a given ROI by magnifying it with a suitable conformal map. We consider two measurement geometries: In the first one, $\Omega=D$, meaning that the magnification is performed by a mere Möbius transform. In the second case, $\Omega$ is a polygon and the needed conformal transform $\Phi$ is thus a Schwarz–Christoffel map [47] (or a composition of a Schwarz–Christoffel map and a Möbius transform as in our numerical implementation).

In the third part, we consider EIT measurements in a half-plane setting; simulation of (relative) EIT data for the half-plane is considered in Appendix B. The idea is to map the upper half-plane onto the unit disk by utilizing a family of Möbius transforms, each of which magnifies different details in the target conductivity. By forming the corresponding reconstructions in the unit disk and mapping them back to the half-plane, one obtains reconstructions that highlight different details in the conductivity of the upper half-plane. In practice, the use of such a family of Möbius transforms can be interpreted as moving an array of electrodes along the horizontal axis; see Appendices A and B.

All measurement data for bounded target domains are simulated by solving the Neumann problem (7) for the employed current patterns by FEM with piecewise linear basis functions. For the half-plane geometry, the (relative) EIT data (for the virtual domain) are simulated by resorting to layer potential techniques, as explained in Appendix B. In both cases, the utilized discretizations are so dense that the (relative) amount of numerical noise is assumed to be insignificant for the employed lowish spatial Fourier frequencies; if not stated otherwise, we use the truncation index $N=16$, for all families of Fourier-like current patterns (cf. (23) and (39)).

For each considered reconstruction domain, the $\tilde{{\mathbf{H}}}_{k}$ matrix needed in (25) is formed on a uniform grid of $128\times 128$ points over the square $-12\leq k_{1},k_{2}\leq 12$, $k=k_{1}+\mathrm{i}k_{2}$. To evaluate (26), we use the trapezoidal rule with 256 uniform nodes. The integral equation (21) is discretized over a uniform $k$-grid of $2^{7}\times 2^{7}$ points. Our choice for the truncation parameters $R$ and $c$ in (28) is somewhat heuristic: We first plot the scattering transform over the full $k$-grid mentioned above, then pick a reasonable $c>0$, and finally choose the smallest $R\in{\mathbb{N}}$ for which the connected area where ${\mathbf{t}}_{R,c}={\mathbf{t}}$ around the origin is as large as possible; see Figure 4 for an example.

Finally, note that all considered conductivities are piecewise constant, which contradicts the $W^{2,p}$ assumption of Nachman’s method [35] but anyway results in reasonable reconstructions (cf. [28, 27, 43]).

5.1. Reconstructions without magnification

We consider three different true domains:

1. (1)

$\Omega_{1}=D$ is the unit disk, 2. (2)

$\Omega_{2}\subset D$ is a rectangle, 3. (3)

$\Omega_{3}\subset D$ is a more complicated polygonal domain.

The target conductivity $\sigma_{1}=\sigma:D\to{\mathbb{R}}_{+}$ for $\Omega_{1}$ is composed of three embedded inclusions in homogeneous unit background; see the top left image of Figure 3. The red inclusion in the middle has constant conductivity 3, whereas the blue and light red ones closer to the boundary are characterized by conductivity levels 0.2 and 2, respectively. The target conductivities for $\Omega_{2}$ and $\Omega_{3}$ are defined as $\sigma_{2}=\sigma|_{\Omega_{2}}$ and $\sigma_{3}=\sigma|_{\Omega_{3}}$, respectively; observe that ${\rm supp}(\sigma-1)\subset\Omega_{j}$, $j=2,3$, so that the target conductivity equals one in some interior neighborhood of the boundary for all three target domains. The bottom row of Figure 3 illustrates the conformally mapped conductivities $\tilde{\sigma}_{j}=\tilde{\sigma}_{j}\circ\Psi_{j}$, $j=2,3$, in the unit disk. Here, $\Psi_{j}$, $j=2,3$, is a Schwarz–Christoffel map sending the unit disk onto $\Omega_{j}$, $j=2,3$, respectively, with the extra requirement $\Psi_{j}(0)=0$ that fixes the map up to rotations of the unit disk.

Figure 4 shows the real parts of the truncated scattering transforms (26) for the five geometries of Figure 3 with $c=10$. The reconstructions produced by the D-bar method are illustrated in Figure 5, with the top row showing the target conductivities for comparison. The reconstructions formed in the original domains using the respective standard Fourier current basis (23) are shown in the middle row. The bottom row presents the reconstructions corresponding to the employment of the conformally transformed Fourier current basis (39) on $\partial\Omega_{j}$, $j=2,3$, subsequently applying the unit disk version of the D-bar method (see Section 3.2) to the resulting boundary data (cf. Lemma (4.1)) to reconstruct the virtual conductivity $\tilde{\sigma}_{j}$, $j=2,3$, in $D$, and finally mapping the reconstruction back to the true domain with the help of $\Psi_{j}$, $j=2,3$. The reconstructions in the second and third rows emphasize different characteristic of the target conductivities. The (single) reconstruction in $\Omega_{1}$ is arguably the best one, but for $\Omega_{2}$ and $\Omega_{3}$ the two reconstruction techniques seem to function approximately as well. Hence, it seems that the choice between the two of them should be made by comparing the respective computational cost, that is, whether the possibility of implementing the D-bar method in the unit disk (see Section 3.2) compensates for the extra cost of forming the conformal map between the true domain and the unit disk. The answer to this question depends on many things such as, e.g., whether the reference N-to-D map corresponding to unit conductivity for the true domain can measured or needs to be computed and whether there exists an efficient technique for evaluating the needed conformal map for the considered measurement geometry.

5.2. Magnification of a ROI

In this section, we test the idea of magnifying a ROI by a conformal map. The two target conductivities and true domains are shown in the top row of Figure 6. The left-hand domain is the unit disk with several embedded circular inhomogeneities; the first quadrant with the three small inclusions is considered as the ROI. The right-hand domain is the polygon $\Omega_{3}$ from the previous section with the same conductivity pattern consisting of three inclusions in a homogeneous background; the bottom antler enclosing the pink inclusion serves as the ROI.

The middle row of Figure 6 shows two virtual conductivities in the unit disk obtained by magnifying the ROI with a certain conformal map; in both cases, the requirement that the black dot in the top row is mapped to the corresponding one in the middle row fixes $\Phi$ (up to rotations of the unit disk). For the left-hand setup, such $\Phi$ is realized by the Möbius transform ${\mathcal{M}}_{a}$ of (40) with the free parameter $a=0.6\exp(\mathrm{i}\theta/4)$. For the right-hand configuration $\Phi$ is a Schwarz–Christoffel map, which we construct by composing a suitable Möbius transform with a preliminary Schwarz–Christoffel map that sends the origin to itself. The highest, i.e. the sixteenth, conformally transformed Fourier current pattern (cf. (39)) for each domain is shown in the bottom row of Figure 6. Driving these current densities through the boundaries of the true domains corresponds to the use of the standard Fourier boundary current $\tilde{\phi}_{16}$ for the virtual configurations in the middle row of Figure 6.

The first row of Figure 7 shows the target conductivities within the ROIs. The conductivity reconstructions produced by the D-bar method implemented directly in the true target domains with the standard Fourier current patterns (23) are presented in the middle row. The bottom row illustrates the reconstructions corresponding to the use of conformally transformed Fourier currents (39) for the true domain and an application of the standard D-bar method in the virtual reconstruction domain (i.e., the unit disk). The reconstructions in Figure 7 demonstrate that the magnification by a suitable conformal map results in a more accurate reconstruction of the conductivity within the ROI. For all reconstructions in Figure 7, we set $c=20$ in (26), while $R>0$ was chosen as explained just before Section 5.1.

5.3. Reconstructions in the half-plane

Our final numerical experiment considers the unbounded case $\Omega=\{z\in{\mathbb{C}}\ |\ {\rm Im}(z)>0\}$; see Appendix B for information on the well posedness of the forward problem of EIT in such a half-space geometry as well as on the employed layer potential techniques for simulating the needed boundary data. The Möbius transforms used for mapping $\Omega$ onto the unit disk are of the form

[TABLE]

where $\xi=1.2\mathrm{i}$ is fixed and the horizontal translation parameter $b\in{\mathbb{R}}$ takes five different values, namely $b=-3,-1.5,0,1.5,3$. In particular, ${\mathcal{M}}_{b}$ maps the point $b+\xi$ to the origin.

The target conductivity consists of three inclusions in homogeneous unit background; it is visualized in the left-hand column of Figure 8 over the subset $\Omega_{0}=\{z\in{\mathbb{C}}\ |\ \left|{\mbox{Re\,}(z)}\right|\leq 10,\ 0\leq\mbox{Im\,}(z)\leq 3\}\subset\Omega$. The right-hand column of Figure 8 shows the virtual conductivities $\tilde{\sigma}=\sigma\circ\Psi$ in the unit disk for the conformal map $\Phi={\mathcal{M}}_{b}$ and our five choices of $b\in{\mathbb{R}}$. The gray area represents ${\mathcal{M}}_{b}(\Omega\setminus\Omega_{0})$. Observe that the parameter $b$ determines the horizontal position that is magnified the most under ${\mathcal{M}}_{b}$ (or, actually, shrunk the least).

For this half-plane setting, we use only ten conformally transformed Fourier currents on $\partial\Omega$, i.e., set $N=5$ in (39). The black dots on the boundary of $\Omega$ in Figure 8 represent a possible configuration of small electrodes for approximating the conformally transformed Fourier currents (39) corresponding to the particular $\Phi={\mathcal{M}}_{b}$; the images of these ‘electrodes’ under ${\mathcal{M}}_{b}$ form an (incomplete) uniform grid on the boundary of the virtual domain as shown in the right-hand column of Figure 8 (cf. Appendix A and B). To be more precise, for each $\Phi={\mathcal{M}}_{b}$ the parts of (the PEM approximation for) the current patterns (39) supported outside the interval $[-2.5+b,2.5+b]$ on the real axis are ignored, i.e., set to zero, and the corresponding potential measurements are also only taken on this interval and set to zero on the rest of $\partial\Omega$. The effect of these truncations is not further analyzed, but increasing the width of the measurement interval would certainly improve the reconstructions to a certain extent — as would increasing the number of employed conformally transformed Fourier currents on $\partial\Omega$.

The reconstructions are once again formed by the D-bar method in the virtual domain (i.e., the unit disk) and mapped back to the half-plane by ${\mathcal{M}}_{b}^{-1}$. The five reconstructions corresponding to our five values for the parameter $b$ are shown in Figure 9, with the vertical gray lines indicating the midpoints of the ‘measurement arrays’, i.e., the five values of $b$. It is obvious that the parameter $b$ pinpoints the area that is reconstructed most accurately: The blue and light red inclusions are clearly visible in the first ($b=-3$) and the last ($b=3$) reconstruction, respectively. In addition, there are arguably some traces of the dark red inclusion, lying the furthest from $\partial\Omega$, in the third reconstruction ($b=0$). One could postprocess the five reconstructions of Figure 9 to obtain a ‘composite reconstruction’, but such image processing is outside the scope of this work.

6. Concluding remarks

By utilizing conformal maps, we have introduced and numerically tested a method for transferring the D-bar method from an arbitrary simply connected planar domain to the unit disk, where the algorithm can be implemented most efficiently. In particular, we have presented the first D-bar reconstructions in a half-plane. By exploiting the degrees of freedom in the choice of a conformal map between a given simply connected planar domain and the unit disk (characterized by the Möbius transforms sending the unit disk onto itself), we have also introduced a method for magnifying a ROI and thus obtaining more accurate reconstructions of its conductivity. The magnification of a ROI is based on the use of so-called conformally transformed Fourier current patterns that exhibit fast spatial oscillations close to the ROI. However, we have proved that idealized continuum boundary data corresponding to such boundary current densities can be accurately approximated by real-world EIT measurements if the available (small) electrodes are positioned in a suitable nonuniform configuration; see Appendix A.

Appendix A Approximations by point electrodes

The aim of this appendix is to explain how the positioning of the available electrodes on the boundary of the true domain is related to the ability to approximate the output of the relative N-to-D map for given virtual boundary current densities on the boundary of the unit disk. To put it short, the electrode pattern should be densest on those boundary sections where the conformally transformed current densities oscillate the most.

For simplicity, we assume in the following that $\partial\Omega$ is of the class $C^{\infty}$, although the considered distributional Neumann boundary value problems could be defined even in domains with piecewise $C^{1,\alpha}$ boundaries [40]. As in Section 4, $\Phi:\Omega\to D$, with $\Psi=\Phi^{-1}$, is a conformal map of the true domain $\Omega$ onto the virtual domain, i.e., the open unit disk. In addition, it is assumed that the conductivity $\sigma$ equals one in some interior neighborhood of $\partial\Omega$. Under these assumptions, the relative N-to-D map

[TABLE]

is bounded for any $r,q\in{\mathbb{R}}$ [32]. As $\tilde{\sigma}=\sigma\circ\Psi$ also equals one in some interior neighborhood of $\partial D$, the analogous result obviously holds for $\tilde{{\mathcal{R}}}_{\tilde{\sigma}}-\tilde{{\mathcal{R}}}_{1}$ as well, with $\tilde{{\mathcal{R}}}_{\tilde{\sigma}},\tilde{{\mathcal{R}}}_{1}$ being N-to-D maps for the unit disk.

The PEM models the electrodes as pointlike boundary current sources. Let $\tilde{z}_{m}\in\partial D$, $m=0,\dots,M-1$, be a set of virtual equiangular electrodes on the boundary of the unit disk and define the corresponding true point electrodes on $\partial\Omega$ to be $z_{m}=\Psi(\tilde{z}_{m})$, $m=0,\dots,M-1$. Hence, the point electrodes for the true domain are the preimages of equidistant electrodes on $\partial D$ under the bijective and smooth map $\Phi:\overline{\Omega}\to\overline{D}$. These fixed electrode positions are used in the following notation: for arbitrary $g\in H^{r}(\partial\Omega)$ and $\tilde{g}\in H^{r}(\partial D)$, $r>1/2$, we write

[TABLE]

In particular, if $g=\tilde{g}\circ\Phi|_{\partial\Omega}$, then $[g]_{M}=[\tilde{g}]_{M}$.

Let ${\mathbb{C}}_{\diamond}^{M}$ be the zero-mean subspace of ${\mathbb{C}}^{M}$ and denote by $\mathcal{I}\in{\mathbb{C}}_{\diamond}^{M}$ a net current pattern through the pointlike electrodes $\{z_{m}\}_{m=0}^{M-1}$ on the true domain boundary $\partial\Omega$. The relative current-to-potential operator $\mathscr{M}:{\mathbb{C}}_{\diamond}^{M}\to{\mathbb{C}}^{M}/{\mathbb{C}}$ of the PEM is defined via [18]

[TABLE]

where $\delta_{z}\in H^{-1/2-\epsilon}(\partial\Omega)$, $\epsilon>0$, is the Dirac delta distribution supported at $z\in\partial\Omega$. In other words, the net electrode currents $\mathcal{I}$ are driven through the pointlike electrodes and the corresponding relative potentials are measured at the same locations. Notice that $\mathscr{M}$ is well-defined due to (41) and the Sobolev embedding theorem. From the practical standpoint, it is important to note that the discrepancy between $\mathscr{M}$ and the relative current-to-potential map of the CEM [10, 45] is of the order $O(d^{2})$, with $d>0$ being the maximal diameter of the electrodes [18].

Suppose a given virtual continuum current pattern $\tilde{f}\in H^{r}_{\diamond}(\partial D)$, $r>1/2$, on the boundary of the unit disk satisfies $[\tilde{f}]_{M}\in{\mathbb{C}}_{\diamond}^{M}$. In particular, this holds if $\tilde{f}$ is an element of the Fourier current basis (23), with $N<M$, on $\partial D$ since $\tilde{z}_{m}$, $m=0,\dots,M-1$, are equiangular. We define the corresponding true electrode net currents as

[TABLE]

The following theorem demonstrates that the virtual relative potential measurement on $\partial D$ corresponding to $\tilde{f}$ can be approximated by the true electrode measurements that correspond to driving the net currents $\mathcal{I}^{\tilde{f}}\in{\mathbb{C}}_{\diamond}^{M}$ through the point electrodes $\{z_{m}\}_{m=0}^{M-1}\subset\partial\Omega$.

Theorem A.1**.**

Assume that $\tilde{f}\in H^{r}_{\diamond}(\partial D)$, $r\in{\mathbb{N}}$ and $[\tilde{f}]_{M}\in{\mathbb{C}}_{\diamond}^{M}$. Then, it holds that

[TABLE]

where $C=C(r,\tilde{\sigma})>0$ is independent of $M\in{\mathbb{N}}\setminus\{1\}$ and $\tilde{f}$.

Proof.

First of all, it follows from the compatibility of conformal maps and Neumann boundary values composed of Dirac deltas that

[TABLE]

which is to be understood modulo an additive constant; see the proof of [17, Theorem 3.2]. To shorten the notation, we denote

[TABLE]

in what follows.

Let $r\in{\mathbb{N}}$ and $v\in H^{r}(\partial D)/{\mathbb{C}}$ be arbitrary, and denote by $v_{0}\in H^{r}_{\diamond}(\partial D)$ the zero-mean element in the equivalence class $v\in H^{r}(\partial D)/{\mathbb{C}}$. In particular, take note that $\|v\|_{H^{r}(\partial D)/{\mathbb{C}}}=\|v_{0}\|_{H^{r}(\partial D)}$. Since the mean-free distribution $\tilde{f}_{M}-\tilde{f}\in H^{-r}_{\diamond}(\partial D)$ does not see an additive constant, we have

[TABLE]

Observe that the right-hand side of (A) corresponds to the quadrature error of the trapezoidal rule for $\tilde{f}v_{0}$ over $\partial D$ on an equidistant grid. Hence, it follows from [31, Theorem 1.1 & Remark 1.2] that

[TABLE]

where $s$ denotes the arclength parameter on $\partial D$ and $C=C(r)>0$. By expanding the derivative on the right-hand side and using the triangle and Schwarz inequalities, one easily sees that

[TABLE]

We now take the supremum over those $v\in H^{r}(\partial D)/{\mathbb{C}}$ that have unit norm in order to conclude that

[TABLE]

for $r\in{\mathbb{N}}$.

Due to the boundedness of $\tilde{{\mathcal{R}}}_{\tilde{\sigma}}-\tilde{{\mathcal{R}}}_{1}:H_{\diamond}^{r}(\partial D)\to H^{q}(\partial D)/{\mathbb{C}}$ for any $q\in{\mathbb{R}}$, it follows from (47) that

[TABLE]

and some $C^{\prime\prime}(q,r,\tilde{\sigma})>0$. The assertion now follows by combining (42), (45) and the Sobolev embedding theorem. ∎

The optimal convergence rate in (44) follows from the equidistant placement of the virtual point electrodes on $\partial D$, that is, from the efficiency of the trapezoidal rule for approximating integrals with periodic integrands (cf., e.g., [31]). This was achieved by choosing the positions of the true electrodes on $\partial\Omega$ to be the conformal preimages of the the virtual equidistant ones. By comparing this observation, e.g., with the Möbius transform (40) used for magnifying a certain ROI in case $\Omega=D$, it becomes obvious that the true electrode pattern on $\partial\Omega$ is the densest close to the ROI.

Appendix B Simulating data for a half-plane

Let $\Omega$ be the (unbounded) open upper half-plane

[TABLE]

Although the general treatment of EIT in unbounded domains would require some extra considerations, here we exclusively restrict our attention to the following setting that is straightforward to handle and suits our purposes [20]:

• I

The boundary current densities are (finite) linear combinations of delta distributions, i.e.,

[TABLE]

where $\mathcal{I}\in{\mathbb{C}}_{\diamond}^{M}$ and $z_{m}$, $m=0,\dots,M-1$, are points on $\partial\Omega$, i.e., on the horizontal axis.

• II

The conductivity $\sigma:\Omega\to{\mathbb{R}}_{+}$ is composed of a finite number of compactly supported $C^{2}$ inclusions with constant conductivity levels in unit background.

Under these assumptions, the Neumann problem (7) has a unique solution $u\in H^{1}_{\rm loc}(\Omega)$ under the decay condition

[TABLE]

In particular, the relative boundary potentials corresponding to current patterns of the form (48), i.e.,

[TABLE]

are smooth and can be simulated using layer potentials. See [20, Appendix] for the details. In consequence, the relative current-to-potential map of the PEM remains well-defined for the half-space geometry.

What is more, the argumentation in the proof of Theorem A.1 remains valid. That is, if (i) $\Psi$ is a conformal map of the unit disk $D$ onto the half-space $\Omega$, i.e., a certain Möbius transform, (ii) the points $\{\tilde{z}_{m}\}_{m=0}^{M-1}$ form an uniform grid on $\partial D$ and (iii) $z_{m}=\Psi(\tilde{z}_{m})$, $m=0,\dots,M-1$, then the estimate (44) is true for $\tilde{f}\in H^{r}_{\diamond}(\partial D)$ satisfying $[\tilde{f}]_{M}\in{\mathbb{C}}_{\diamond}^{M}$. (The only part in the proof of Theorem A.1 that actually refers to $\Omega$ is (45) which can be relatively easily validated for a half-space.)

When generating the data for the half-space geometry in our numerical examples, we first simulate the relative current-to-potential operator of the PEM (42) by employing the boundary element code from [20]. Subsequently, we approximate the grid values of the (virtual) relative potential measurements

[TABLE]

for the Fourier current basis $\tilde{f}=\tilde{\phi}_{n}$, $1\leq|n|\leq N$, on the boundary of the virtual reconstruction domain $D$ with the help of (44). These can then be used for approximating the Fourier basis representation for the (relative) D-to-N map on the boundary of the virtual domain $\partial D$.