# Linearization-based direct reconstruction for EIT using triangular Zernike
decompositions

Antti Autio  Department of Mathematics and Systems Analysis, Aalto University,
P.O. Box 11100, 00076 Helsinki, Finland<EMAIL_ADDRESS>,  Henrik
Garde  Department of Mathematics, Aarhus University, Ny Munkegade 118, 8000
Aarhus C, Denmark<EMAIL_ADDRESS>,  Markus Hirvensalo  Department of
Mathematics and Systems Analysis, Aalto University, P.O. Box 11100, 00076
Helsinki, Finland<EMAIL_ADDRESS>and  Nuutti Hyvönen  Department
of Mathematics and Systems Analysis, Aalto University, P.O. Box 11100, 00076
Helsinki, Finland<EMAIL_ADDRESS>
###### Abstract.

This work implements and numerically tests the direct reconstruction algorithm
introduced in [Garde & Hyvönen, SIAM J. Math. Anal., 2024] for two-dimensional
linearized electrical impedance tomography. Although the algorithm was
originally designed for a linearized setting, we numerically demonstrate its
functionality when the input data is the corresponding change in the current-
to-voltage boundary operator. Both idealized continuum model and practical
complete electrode model measurements are considered in the numerical studies,
with the examined domain being either the unit disk or a convex polygon.
Special attention is paid to regularizing the algorithm and its connections to
the singular value decomposition of a truncated linearized forward map, as
well as to the explicit triangular structures originating from the properties
of the employed Zernike polynomial basis for the conductivity.

Keywords: electrical impedance tomography, linearization, direct
reconstruction, Zernike polynomials, singular value decomposition

2020 Mathematics Subject Classification: 65N21, 65N20, 35R30, 35R25

##  1\. Introduction

In electrical impedance tomography (EIT) the aim is to gain information on the
conductivity distribution inside a physical body from current and potential
measurements on its boundary; see the review articles [2, 4, 22] and the
references therein for more information on EIT. The reconstruction problem of
EIT can be mathematically formulated as the ill-posed task of inverting the
nonlinear operator that maps the internal conductivity of the examined body to
the boundary measurements.

The main goal of this paper is to implement and numerically test the direct
reconstruction algorithm for two-dimensional EIT introduced in [10, Theorem
1.3], capable of rapidly reconstructing any (even unbounded)
$L^{2}$-perturbation to a homogeneous background conductivity. The article
[10] only considered the idealized continuum model (CM) of EIT that allows
current input and potential measurement everywhere on the object boundary. In
this work, we also test the algorithm on data simulated by resorting to the
complete electrode model (CEM) [5, 20] as well as on real-world data from a
water tank [13]. Moreover, even though the reconstruction algorithm was
originally designed for a linearized version of EIT, our numerical studies are
exclusively based on difference data corresponding to the nonlinear forward
map of either the CM or the CEM. The aim is to demonstrate that an algorithm
based on a linearization can systematically produce good quality
reconstructions in EIT.

The basic version of the algorithm is formulated assuming the examined object
is modeled as the unit disk $D$, but it can also be applied to any smooth
enough simply-connected domain by resorting to the Riemann mapping theorem, as
explicitly outlined in the paragraph following [10, Theorem 1.3]. In the
framework of the linearized CM, one assumes that the available data is the
Fréchet derivative of the forward map, evaluated at the unit conductivity, in
the direction of the to-be-reconstructed conductivity perturbation. This is a
linear operator sending input currents to measured boundary potentials. To be
more precise, the algorithm takes as its input the elements in the infinite
matrix representation of this boundary map with respect to the standard
Fourier basis. As hinted above, in our numerical studies these matrix elements
are replaced by those corresponding to the actual change in the Neumann-to-
Dirichlet (ND) map when the conductivity inside $D$ is perturbed. The
reconstruction algorithm itself is comprised of a set of explicit linear
formulas that map the aforementioned matrix elements to the coefficients,
$c_{j,k}$, in the expansion of the conductivity perturbation in the Zernike
polynomial basis [24]; see also [1] for a similar approach.

The Zernike polynomials are indexed by two indices: $j\in\mathbb{Z}$
controlling the Fourier frequency in the polar angle and $k\in\mathbb{N}_{0}$
controlling (for a fixed $j$) the behavior in the radial variable, with
increasing $k$ indicating more oscillations closer to the origin. As
explicitly proven in [10], reconstructing the coefficients $c_{j,k}$ becomes
more unstable as $k$ increases. On the other hand, it turns out that $c_{j,k}$
with a fixed angular index $j$ only affect the $j$th diagonal in the matrix
representation of the associated Fréchet derivative. In consequence, a certain
angular frequency in the conductivity perturbation can be reconstructed by
only considering a single diagonal in the data matrix, and the whole
reconstruction algorithm can thus be decoupled into independent smaller tasks
indexed by $j$. What is more, the linear dependence of the $j$th diagonal in
the data matrix on $c_{j,k}$ can be represented by an infinite triangular
matrix. The formulas comprising the direct reconstruction algorithm correspond
to explicitly inverting these triangular matrices for all $j\in\mathbb{Z}$.
Note that singular vectors of the linearized forward map must also obey such a
decoupling, that is, a right-hand singular vector may only contain a single
angular frequency and a left-hand singular vector can only have nonzero
elements on a single diagonal of its matrix representation.

The decoupling of the angular frequencies allows, in particular, to separately
regularize the subsystems corresponding to different angular indices $j$. This
leads to increasing level of regularization as a function of $\lvert j\rvert$.
Numerically considering radial indices beyond, say, $k=15$, is impractical due
to the increasing ill-posedness of the subproblems for higher radial indices.
Hence, any reasonable regularization method can be used for the decoupled
subproblems without loosing the instantaneous nature of the reconstruction
algorithm. We consider two simple approaches accompanied by the Morozov
discrepancy principle:

  1. (i)

truncated singular value decomposition (SVD) and

  2. (ii)

choosing the sizes of the decoupled subsystems as a function of $\lvert
j\rvert$ by monitoring the absolute values of the elements on their diagonals.

The former has a better theoretical justification, but it cannot fully exploit
the triangular structure of the decoupled systems as can the latter.

###  1.1. Article structure

Section 2 introduces the linearized CM for the unit disk and extends it for
square-integrable conductivity perturbations following [10]. In Section 3, the
linearized forward map is written out explicitly with respect to the Zernike
polynomial basis for the conductivity perturbation and the Fourier basis for
the boundary measurements. In addition, the triangular structures of the
subsystems corresponding to different angular Zernike indices are studied. The
implementation of regularized reconstruction algorithms are considered in
Section 4, while Section 5 explains how data for the algorithm can be
measured/simulated on polygonal domains or based on the CEM. The numerical
experiments are presented in Section 6.

##  2\. Linearization of the CM of EIT

This section introduces the linearization of the CM with respect to square-
integrable conductivity perturbations following the material in [10]. In
particular, our approach only requires considering linearizations for the CM
in the unit disk $D\subset\mathbb{R}^{2}$. In the numerical experiments with
polygonal domains or the CEM, we rely on methods for replicating or
approximating CM boundary data for the unit disk based on (electrode)
measurements on other simply-connected planar domains; see [8, 17].

Consider an isotropic conductivity

$\gamma\in L^{\infty}_{+}(D)\coloneqq\big{\\{}\kappa\in
L^{\infty}(D):\text{ess}\inf(\mathop{\textup{Re}}\kappa)>0\big{\\}}$

in $D$ and a normal current density

$f\in L^{2}_{\diamond}(\partial D)\coloneqq\big{\\{}g\in L^{2}(\partial
D):\langle g,1\rangle_{L^{2}(\partial D)}=0\big{\\}}$

on $\partial D$. According to the CM, the corresponding electric potential $u$
in $D$ weakly satisfies a Neumann problem for the conductivity equation:

$-\nabla\cdot(\gamma\nabla u)=0\quad\text{in }D,\qquad\nu\cdot(\gamma\nabla
u)=f\quad\text{on }\partial D,$ (2.1)

where $\nu$ is the exterior unit normal of $\partial D$. The variational form
of (2.1) is to find $u\in H^{1}(D)$ that satisfies

$\int_{D}\gamma\nabla u\cdot\overline{\nabla v}\,\mathrm{d}x=\langle
f,v\rangle_{L^{2}(\partial D)}\qquad\text{for all }v\in H^{1}(D),$ (2.2)

for which the Lax–Milgram lemma yields a unique solution

$u^{\gamma}_{f}\in H^{1}_{\diamond}(D)\coloneqq\\{w\in H^{1}(D)\,:\,\langle
w|_{\partial D},1\rangle_{L^{2}(\partial D)}=0\\}.$

Note that the use of the subspace $H^{1}_{\diamond}(D)\subset H^{1}(D)$
corresponds to a systematic way of choosing the ground level of potential.

The associated ND, or current-to-voltage, map

$\Lambda(\gamma):f\mapsto u^{\gamma}_{f}|_{\partial D}$ (2.3)

is a compact self-adjoint operator in $\mathscr{L}(L^{2}_{\diamond}(\partial
D))$. In fact, in the considered two-dimensional setting with a smooth domain,
$\Lambda(\gamma)$ belongs to the space $\mathscr{L}_{\rm
HS}(L^{2}_{\diamond}(\partial D))$ of Hilbert–Schmidt operators. The forward
map of the CM is defined as

$\Lambda:\left\\{\begin{array}[]{l}\gamma\mapsto\Lambda(\gamma),\\\\[2.84526pt]
L^{\infty}_{+}(D)\to\mathscr{L}_{\rm HS}(L^{2}_{\diamond}(\partial
D))\subset\mathscr{L}(L^{2}_{\diamond}(\partial D)).\end{array}\right.$

It is well-known that $\Lambda$ is Fréchet differentiable at any $\gamma\in
L^{\infty}_{+}(D)$, with the corresponding derivative $D\Lambda(\gamma)$ in
$\mathscr{L}(L^{\infty}(D),\mathscr{L}_{\rm HS}(L^{2}_{\diamond}(\partial
D)))\subset\mathscr{L}(L^{\infty}(D),\mathscr{L}(L^{2}_{\diamond}(\partial
D)))$; see, e.g., [9].

In this work, we are interested in the Fréchet derivative of $\Lambda$ at the
unit conductivity, i.e., $F=D\Lambda(1)$. $F$ satisfies

$\langle(F\eta)f,g\rangle_{L^{2}(\partial D)}=-\int_{D}\eta\nabla
u_{f}^{1}\cdot\overline{\nabla u_{g}^{1}}\,\mathrm{d}x$ (2.4)

for any $\eta\in L^{\infty}(D)$ and all $f,g\in L^{2}_{\diamond}(\partial D)$.
In fact, (2.4) uniquely characterizes $F$ as an element of
$\mathscr{L}(L^{\infty}(D),\mathscr{L}(L^{2}_{\diamond}(\partial D)))$, and it
can thus be used as a definition.

According to the elliptic regularity theory, the gradients of the harmonic
functions $u_{f}^{1}$ and $u_{g}^{1}$ belong to $H^{1/2}(D)$, with their
norms, respectively, bounded by the $L^{2}(\partial D)$-norms of $f$ and $g$.
By the continuous embedding $H^{1/2}(D)\hookrightarrow L^{4}(D)$, the formula
(2.4) can thus be used to extend $F$ to an operator in
$\mathscr{L}(L^{2}(D),\mathscr{L}(L^{2}_{\diamond}(\partial D)))$. This is how
we interpret $F$ in what follows; consult [10] for the details of this
extension argument.

##  3\. Explicit representation for $F$ and its inverse

In this section, we continue to assume that the investigated domain is the
unit disk $D$ and give an explicit representation for the bounded linear map

$F:\left\\{\begin{array}[]{l}\eta\mapsto D\Lambda(1)\eta,\\\\[2.84526pt]
L^{2}(D)\to\mathscr{L}(L^{2}_{\diamond}(\partial D))\end{array}\right.$

with respect to the Zernike polynomial basis [24] for $L^{2}(D)$ and the
standard Fourier basis for $L^{2}_{\diamond}(\partial D)$. In fact, we also
present an explicit characterization for the inverse of $F$ on its range. Take
note that the injectivity of $F$ follows from the original uniqueness proof of
Calderón [3] that generalizes for square-integrable conductivity perturbations
in two dimensions; see also [10, Section 6] for an alternative proof. We still
partially follow the material and presentation in [10].

###  3.1. Zernike polynomials and Fourier basis

As in [10], we define the Zernike polynomials in the polar coordinates
$(r,\theta)$ as

$\psi_{j,k}(r,\theta)=\sqrt{\frac{\lvert j\rvert+2k+1}{\pi}}R^{\lvert
j\rvert}_{\lvert j\rvert+2k}(r)\,\mathrm{e}^{\mathrm{i}j\theta},\qquad
j\in\mathbb{Z},\ k\in\mathbb{N}_{0},$ (3.1)

where

$R^{\lvert j\rvert}_{\lvert j\rvert+2k}(r)=\sum^{k}_{i=0}(-1)^{i}\binom{\lvert
j\rvert+2k-i}{i}\binom{\lvert j\rvert+2k-2i}{k-i}\,r^{\lvert j\rvert+2k-2i}$

is a Zernike radial function. As mentioned in Section 1, the two indices
$j\in\mathbb{Z}$ and $k\in\mathbb{N}_{0}$ control the angular and radial
characteristics of the Zernike polynomials, with larger absolute values of $j$
and $k$ indicating faster spatial oscillations.

When scaled as in (3.1), the Zernike polynomials form an orthonormal basis for
$L^{2}(D)$, and they are used for representing the conductivity perturbation
$\eta$ as

$\eta=\sum_{j=-\infty}^{\infty}\eta_{j}=\sum_{j=-\infty}^{\infty}\sum_{k=0}^{\infty}c_{j,k}\psi_{j,k},$
(3.2)

where we call $\eta_{j}$ the _$j$ th angular frequency in $\eta$_, defined by

$\eta_{j}=\sum_{k=0}^{\infty}c_{j,k}\psi_{j,k},\qquad j\in\mathbb{Z}.$

In the following, we implicitly assume such a Zernike polynomial expansion
with square-summable coefficients
$\\{c_{j,k}\\}_{j\in\mathbb{Z},k\in\mathbb{N}_{0}}\subset\mathbb{C}$ for the
to-be-reconstructed conductivity perturbation $\eta$.

The employed orthonormal basis for $L^{2}_{\diamond}(\partial D)$ is a
standard Fourier basis

$f_{m}(\theta)=\frac{1}{\sqrt{2\pi}}\mathrm{e}^{\mathrm{i}m\theta},\qquad
m\in\mathbb{Z}^{\prime},$

where $\mathbb{Z}^{\prime}=\mathbb{Z}\setminus\\{0\\}$ in accordance with all
elements of $L^{2}_{\diamond}(\partial D)$ having zero mean. With these
choices, the mapping $F:L^{2}(D)\to\mathscr{L}(L^{2}_{\diamond}(\partial D))$
is completely characterized by the coefficients

$a^{j,k}_{m,n}=\langle(F\psi_{j,k})f_{m},f_{n}\rangle_{L^{2}(\partial
D)},\qquad j\in\mathbb{Z},\ k\in\mathbb{N}_{0},\ m,n\in\mathbb{Z}^{\prime}.$
(3.3)

Indeed, for a given conductivity perturbation $\eta$ expanded as in (3.2), the
mapping $F\eta\in\mathscr{L}(L^{2}_{\diamond}(\partial\Omega))$ is represented
in the Fourier basis as an infinite matrix given element-wise as

$a_{m,n}(\eta)=\langle(F\eta)f_{m},f_{n}\rangle_{L^{2}(\partial
D)}=\sum_{j=-\infty}^{\infty}a_{m,n}(\eta_{j})=\sum_{j=-\infty}^{\infty}\sum_{k=0}^{\infty}a^{j,k}_{m,n}c_{j,k},\qquad
m,n\in\mathbb{Z}^{\prime},$ (3.4)

where

$a_{m,n}(\eta_{j})=\sum_{k=0}^{\infty}a^{j,k}_{m,n}c_{j,k},\qquad
j\in\mathbb{Z},\ m,n\in\mathbb{Z}^{\prime},$ (3.5)

only depends on the Zernike coefficients defining $\eta_{j}$. In particular,
$\\{a_{m,n}(\eta)\\}_{m,n\in\mathbb{Z}^{\prime}}$ can be considered the data
for the linearized inverse problem of the CM in the unit disk.

###  3.2. Matrix representation of the linearized model

Our reconstruction algorithm is based on the fundamental observation that most
matrix elements in (3.3) vanish. Indeed, [10, Eq. (4.5)] immediately gives

$a^{j,k}_{m,n}=-\frac{1}{\sqrt{\pi}}\,\frac{\sqrt{\lvert
j\rvert+2k+1}}{\min\\{\lvert m\rvert,\lvert n\rvert\\}+\lvert
j\rvert+k}\,\prod^{k}_{i=1}\frac{\min\\{\lvert m\rvert,\lvert
n\rvert\\}-i}{\lvert j\rvert+\min\\{\lvert m\rvert,\lvert n\rvert\\}+k-i}$
(3.6)

if $n=m+j$, $k<\min\\{\lvert m\rvert,\lvert n\rvert\\}$ and $mn>0$, and
$a^{j,k}_{m,n}=0$ for all other combinations of $j\in\mathbb{Z}$,
$k\in\mathbb{N}_{0}$ and $m,n\in\mathbb{Z}^{\prime}$. When $k=0$, the product
in (3.6) is defined to take the value $1$. In particular, if
$\\{a^{j,k}_{m,n}\\}_{n,m\in\mathbb{Z}^{\prime}}$ are interpreted as elements
of an infinite matrix (cf. (3.4)), then for a fixed $j$ but any $k$, all
nonzero elements in the matrix lie on its $j$th diagonal. This means that the
infinite matrix defined by (3.5) also only has nonzero elements on its $j$th
diagonal, so by linearity we have

$a_{m,m+j}(\eta)=a_{m,m+j}(\eta_{j}),\qquad j\in\mathbb{Z},\
m\in\mathbb{Z}^{\prime}.$ (3.7)

In other words, the $j$th angular frequency $\eta_{j}$ in the conductivity
perturbation, i.e., the Zernike coefficients
$\\{c_{j,k}\\}_{k\in\mathbb{N}_{0}}$, only contributes to the $j$th diagonal
in the matrix representation (3.4) of $F\eta$, which is the decoupling in the
angular frequencies we were looking for.

To explain how this decoupling can be exploited in solving for the
coefficients $\\{c_{j,k}\\}_{j\in\mathbb{Z},k\in\mathbb{N}_{0}}$ of $\eta$
from the data $\\{a_{m,n}(\eta)\\}_{m,n\in\mathbb{Z}^{\prime}}$, let us define
two infinite vectors indexed by the angular index $j\in\mathbb{Z}$:

$\mathbf{c}^{j}\coloneqq[c_{j,k-1}]_{k=1}^{\infty},\qquad\mathbf{a}^{j}\coloneqq\left\\{\begin{array}[]{ll}\big{[}a_{m,m+j}(\eta)]_{m=1}^{\infty}&\quad\text{if
}j\geq 0,\\\\[5.69054pt]
\big{[}a_{-m,-m+j}(\eta)]_{m=1}^{\infty}&\quad\text{if
}j<0.\end{array}\right.$ (3.8)

The Zernike coefficients defining $\eta_{j}$ form $\mathbf{c}^{j}$, whereas
$\mathbf{a}^{j}$ is composed of certain elements on the $j$th diagonal of
$\\{a_{m,n}(\eta)\\}_{m,n\in\mathbb{Z}^{\prime}}$. Our aim is to present the
linear relationship between $\mathbf{a}^{j}$ and $\mathbf{c}^{j}$ as a
multiplication by an infinite lower triangular matrix.

Let us start with the case $j\geq 0$. According to (3.7) and (3.5), and since
$a^{j,k}_{m,n}=0$ unless $n=m+j$, $k<\min\\{\lvert m\rvert,\lvert n\rvert\\}$
and $mn>0$,

$\mathbf{a}^{j}_{m}=\sum_{k=0}^{\infty}a^{j,k}_{m,m+j}c_{j,k}=\sum_{k=0}^{m-1}a^{j,k}_{m,m+j}c_{j,k}=\sum_{k=1}^{m}a^{j,k-1}_{m,m+j}\mathbf{c}^{j}_{k},\qquad
m\in\mathbb{N},$ (3.9)

which is the sought-for infinite triangular system for $j\geq 0$. Analogously,
for $j<0$,

$\mathbf{a}^{j}_{m}=\sum_{k=0}^{\infty}a^{j,k}_{-m,-m+j}c_{j,k}=\sum_{k=0}^{m-1}a^{j,k}_{-m,-m+j}c_{j,k}=\sum_{k=1}^{m}a^{j,k-1}_{-m,-m+j}\mathbf{c}^{j}_{k},\qquad
m\in\mathbb{N}.$ (3.10)

Due to obvious symmetries in (3.6), $a^{-\lvert j\rvert,k-1}_{-m,-m-\lvert
j\rvert}=a^{\lvert j\rvert,k-1}_{m,m+\lvert j\rvert}$, and thus the
coefficients of the triangular systems in (3.9) and (3.10) only depend on $j$
via $\lvert j\rvert$. Hence,

$\displaystyle F^{\lvert j\rvert}\mathbf{c}^{j}$
$\displaystyle=\mathbf{a}^{j},\qquad j\in\mathbb{Z},$ (3.11)

where the $\lvert j\rvert$-dependent infinite lower triangular system matrix
is defined by

$F^{\lvert j\rvert}_{m,k}=\left\\{\begin{array}[]{ll}a^{\lvert
j\rvert,k-1}_{m,m+\lvert j\rvert}&\text{if }1\leq k\leq
m\in\mathbb{N},\\\\[2.84526pt] 0&\text{otherwise}.\end{array}\right.$ (3.12)

In particular, note that the diagonal elements $F^{\lvert j\rvert}_{m,m}$ are
nonzero for all $m\in\mathbb{N}$ and $j\in\mathbb{Z}$, which means that one
can, in principle, solve for the components of $\mathbf{c}^{j}$ via forward
substitution.

######  Remark 3.1.

Observe that less than a half of the $j$th diagonal of
$\\{a_{m,n}(\eta)\\}_{n,m\in\mathbb{Z}^{\prime}}$ is included in
$\mathbf{a}^{j}$. As this is already sufficient for uniquely determining
$\eta_{j}$, the data that is left out must be redundant. Indeed, it follows
straightforwardly from (3.4) and (3.6) that the $j$th diagonal of
$\\{a_{m,n}(\eta)\\}_{n,m\in\mathbb{Z}^{\prime}}$ is symmetric with respect to
its midpoint, i.e., the (possibly virtual) element $a_{-j/2,j/2}(\eta)$.
Moreover, the $\lvert j\rvert$ central elements on the $j$th diagonal vanish
since their indices have opposite signs (see the text after (3.6)). In fact,
two quadrants of $\\{a_{m,n}(\eta)\\}_{n,m\in\mathbb{Z}^{\prime}}$ are zeros,
and either of the two remaining ones carry all information contained in the
matrix.

In particular, one could equivalently define the $j$th data vector by
employing the opposite ends of the diagonals:

$\mathbf{a}^{j}=\left\\{\begin{array}[]{ll}\big{[}a_{-m-j,-m}(\eta)]_{m=1}^{\infty}&\quad\text{if
}j\geq 0,\\\\[5.69054pt] \big{[}a_{m-j,m}(\eta)]_{m=1}^{\infty}&\quad\text{if
}j<0.\end{array}\right.$ (3.13)

If the whole matrix $\\{a_{m,n}(\eta)\\}_{n,m\in\mathbb{Z}^{\prime}}$ is
available, this observation makes it possible to slightly reduce the effect of
measurement noise by using the average of the two equivalent versions of
$\mathbf{a}^{j}$ in (3.8) and (3.13) as the data in numerical considerations.

###  3.3. Inversion of the infinite triangular systems

The angularly decoupled triangular systems in (3.11) imply a direct algorithm
for reconstructing the coefficients $c_{j,k}$ one angular frequency at a time,
by resorting to simple forward substitution. Since the resulting formulas have
already been worked out in [10, Theorem 1.3] — although without emphasizing
the underlying triangular structure or identifying the method as forward
substitution —, we settle for writing the resulting reconstruction formulas in
our notation:

$\displaystyle\mathbf{c}^{j}_{k+1}$ $\displaystyle=-\sqrt{\pi(\lvert
j\rvert+2k+1)}\;\binom{\lvert j\rvert+2k}{k}\;\mathbf{a}^{j}_{k+1}$
$\displaystyle\phantom{{}=}-\sum^{k}_{i=1}\mathbf{c}^{j}_{i}\frac{\sqrt{(\lvert
j\rvert+2k+1)(\lvert j\rvert+2i-1)}}{\lvert j\rvert+k+i}\,\binom{\lvert
j\rvert+2k}{k-i+1},\qquad k=0,1,2,\dots\,,$ (3.14)

for any $j\in\mathbb{Z}$. When $k=0$, the sum on the right-hand side of (3.3)
is interpreted to vanish.

By virtue of the definition of $\mathbf{a}^{j}$ via (3.4) and (3.8), one only
needs to apply two boundary currents $f_{\pm 1}$ and measure (all Fourier
frequencies in) the corresponding potentials to reconstruct
$c_{j,0}=\mathbf{c}^{j}_{1}$ for all $j\in\mathbb{Z}$. Moving to a next
$k$-level always ‘costs’ applying two more Fourier currents in increasing
order with respect to the absolute value of their angular frequencies. On the
other hand, if one resorts to the dual definition of $\mathbf{a}^{j}$ in
(3.13), reconstructing $c_{j,0}$ for all $j\in\mathbb{Z}$ requires applying
all Fourier boundary currents but only recording the two lowest Fourier modes
in the resulting boundary potentials. In this case, moving to a next $k$-level
always requires recording two more Fourier frequencies in all measured
boundary potentials.

##  4\. Regularization for truncated systems

This section discusses the implementation of regularized reconstruction
algorithms based on truncations of the triangular systems (3.11). To this end,
we introduce the notation

$\mathbb{Z}_{M}\coloneqq\\{-M,\dots,M\\}\quad\textup{and}\quad\mathbb{Z}^{\prime}_{M}\coloneqq\mathbb{Z}_{M}\setminus\\{0\\}$

and assume that the available data is the truncated matrix
$\\{a_{m,n}\\}_{m,n\in\mathbb{Z}^{\prime}_{M}}$ for some fixed
$M\in\mathbb{N}$. This means that we assume to be able to apply the Fourier
boundary currents indexed by $\mathbb{Z}^{\prime}_{M}$ and record the
corresponding Fourier frequencies in the resulting relative boundary
potentials. As a consequence, the infinite linear systems (3.11) are to be
replaced by finite-dimensional ones due to a shortage of data:

$\displaystyle F^{\lvert j\rvert,M}\mathbf{c}^{j,M}$
$\displaystyle=\mathbf{a}^{j,M},$ (4.1)

where $F^{\lvert j\rvert,M}$ is a lower triangular square matrix and
$\mathbf{a}^{j,M}$ and $\mathbf{c}^{j,M}$ are, respectively, composed of a
certain number of first components in $\mathbf{a}^{j}$ and $\mathbf{c}^{j}$
defined by (3.8).

The length of the vectors $\mathbf{a}^{j,M}$ vary as a function of $j$ as they
correspond to different diagonals of a finite square matrix. Indeed,
$\mathbf{a}^{j,M}\in\mathbb{C}^{M-\lvert j\rvert}$, which also reveals that
the systems (4.1) are nonempty only for $j\in\mathbb{Z}_{M-1}$. To summarize,

$F^{\lvert j\rvert,M}\in\mathbb{R}^{(M-\lvert j\rvert)\times(M-\lvert
j\rvert)},\quad\mathbf{a}^{j,M},\mathbf{c}^{j,M}\in\mathbb{C}^{M-\lvert
j\rvert},\qquad j\in\mathbb{Z}_{M-1},$

defined, respectively, by the top left corner of $F^{\lvert j\rvert}$ in
(3.12) and the first $M-\lvert j\rvert$ components of the vectors in (3.8).
Thanks to the triangular structure of (4.1), the truncated systems (4.1) are
still uniquely solvable. In the framework of the linearized CM, they are also
exact in the sense that the components of the solution vectors
$\mathbf{c}^{j,M}$ are precisely the corresponding Zernike coefficients of the
to-be-reconstructed conductivity perturbation $\eta$. In particular, the
explicit inversion formulas (3.3) are still valid, but only for
$j\in\mathbb{Z}_{M-1}$ and $k=0,\dots,M-\lvert j\rvert-1$.

###  4.1. Truncated SVD

The first considered regularization technique for the $2M-1$ triangular
systems (4.1) is truncated SVD combined with the Morozov discrepancy
principle. Let us denote by
$\mathbf{a}^{M},\mathbf{c}^{M}\in\mathbb{C}^{M^{2}}$, respectively, the
concatenations of the data vectors $\mathbf{a}^{j,M}\in\mathbb{C}^{M-\lvert
j\rvert}$ and Zernike coefficient vectors
$\mathbf{c}^{j,M}\in\mathbb{C}^{M-\lvert j\rvert}$ for $j\in\mathbb{Z}_{M-1}$.
The systems (4.1) can be written as a single matrix equation

$F^{M}\mathbf{c}^{M}=\mathbf{a}^{M},$ (4.2)

where

$F^{M}={\rm
diag}\big{(}F^{M-1,M},F^{M-2,M},\dots,F^{0,M},\dots,F^{M-2,M},F^{M-1,M}\big{)}\in\mathbb{R}^{M^{2}\times
M^{2}}$ (4.3)

is a block diagonal matrix whose structure is visualized in Fig. 4.1 for
$M=10$.

Figure 4.1. The structure of $F^{M}$, for $M=10$, as a block diagonal matrix
with lower triangular blocks.

As a consequence, singular systems for the individual $F^{l,M}$,
$l=0,\dots,M-1$, completely define an SVD for the whole $F^{M}$. Indeed, each
singular value of any $F^{l,M}$, $l=0,\dots,M-1$, is also a singular value for
$F^{M}$, with the corresponding singular vectors of $F^{M}$ being suitably
zero-padded versions of those of $F^{l,M}$. Moreover, excluding the case
$l=0$, each singular value of $F^{l,M}$, $l=1,\dots,M-1$, appears twice in an
SVD for $F^{M}$ since $F^{l,M}$ also appears twice on the diagonal of $F^{M}$.

Computing an SVD for $F^{M}$ thus boils down to computing such for the $M$
small matrices $F^{l,M}$, $l=0,\dots,M-1$, i.e., for slightly more than a half
of the diagonal blocks defining $F^{M}$. This is computationally tractable in
practice. However, to enable adopting truncated SVD as a regularization
method, the singular values of the complete system matrix $F^{M}$ also need to
be ordered. For the sake of symmetry, we always include in a truncated SVD of
$F^{M}$ either both or neither of the two right-left singular vector pairs of
$F^{M}$ associated to the same singular value of a repeated block $F^{l,M}$,
$l=1,\dots,M-1$, in (4.3). This leads to $M(M+1)/2$ potentially different
singular values for $F^{M}$. The ordering of these $M(M+1)/2$ singular values
is then achieved with the help of a mapping

$\tau:\\{1,2,\dots,M(M+1)/2\\}\to\mathbb{N}_{0}^{M}$

that works as follows: $\tau_{l}(p)$, i.e., the $l$th component of $\tau(p)$,
indicates how many of the largest $p$ singular values of $F^{M}$ correspond to
$F^{l,M}$. The map $\tau$ can be formed and tabulated offline if $M$ is known
prior to taking the measurements.

Let $(F^{\lvert j\rvert,M})^{\dagger}_{q}$ denote a truncated Moore–Penrose
pseudoinverse of $F^{\lvert j\rvert,M}$, i.e., the pseudoinverse that has been
formed after all but the largest $q\in\\{0,\dots,M-\lvert j\rvert\\}$ singular
values of $F^{\lvert j\rvert,M}$ have been replaced by zeros in an SVD for
$F^{\lvert j\rvert,M}$. For $p\in\\{1,2,\dots,M(M+1)/2\\}$, acting as a
truncation index/regularization parameter, we define the truncated SVD
solution for (4.1) via

$\mathbf{c}^{j,M,p}=(F^{\lvert j\rvert,M})^{\dagger}_{\tau_{\lvert
j\rvert}(p)}\mathbf{a}^{j,M},\qquad j\in\mathbb{Z}_{M-1}.$

The complete truncated SVD solution $\mathbf{c}^{M,p}\in\mathbb{C}^{M^{2}}$
for (4.2) is then defined as the concatenation of $\mathbf{c}^{j,M,p}$ for
$j\in\mathbb{Z}_{M-1}$.

Assume now that the available data is a noisy version
$\mathbf{a}^{M}_{\delta}$ of $\mathbf{a}^{M}$, approximately satisfying

$\lVert\mathbf{a}^{M}_{\delta}-\mathbf{a}^{M}\rVert_{2}=\delta>0,$ (4.4)

where $\lVert\,\cdot\,\rVert_{2}$ denotes the Euclidean norm.111Employing the
Euclidean norm for the diagonals of the data matrix corresponds to the
Frobenius norm for the matrix itself, which is in line with the
Hilbert–Schmidt topology for the ND map in [10]. Let
$\mathbf{c}^{M,p}_{\delta}$ be a truncated SVD solution corresponding to the
noisy data $\mathbf{a}^{M}_{\delta}$. Truncated SVD solutions are combined
with the Morozov discrepancy principle in the standard manner, i.e., by
choosing the smallest truncation index $p\in\\{1,2,\dots,M(M+1)/2\\}$ that
satisfies

$\lVert
F^{M}\mathbf{c}^{M,p}_{\delta}-\mathbf{a}^{M}_{\delta}\rVert_{2}\leq\omega\delta,$
(4.5)

where $\omega\geq 1$ is a ‘fudge’ factor that aims to account for
linearization and numerical errors. The choice of $\omega$ is considered in
Section 6.

###  4.2. Truncation of the triangular subsystems

A downside of truncated SVD is that it does not fully exploit the triangular
structure of (4.1) and, in particular, the explicit solution formula (3.3).
Our second regularization approach aims to address this shortcoming by
directly resorting to (3.3), but only up to a $\lvert j\rvert$-dependent
threshold for the radial index $k$.

It is straightforward to prove that the diagonal elements in each $F^{j,M}$,
$j=0,\dots,M-1$, are monotonically decreasing in their absolute values.
Running the recursion (3.3) until, say, $k+1=q\leq M-\lvert j\rvert$
corresponds to applying forward substitutions to operate on the first $q$
elements in $\mathbf{a}^{j,M}$ with the inverse of the triangular square
matrix composed of the first $q$ rows and columns in $F^{\lvert j\rvert,M}$.
The most unstable single step in this procedure is thus dividing by the last
involved diagonal element $F^{\lvert j\rvert,M}_{q,q}$. Motivated by this
observation, our second (more heuristic) regularization strategy is based on
satisfying the Morozov discrepancy principle by running (3.3) for all
$j\in\mathbb{Z}_{M-1}$ in such a way that the smallest absolute value of a
diagonal element involved in a forward substitution is maximized. Although
such an idea does not fully account for the accumulation of errors in forward
substitutions, the smallest absolute value of a diagonal element in a lower
triangular matrix is in any case closely related to the stability of forward
substitutions for solving a corresponding linear system; see [14, Chapter 8]
for more information.

To simplify our notation, let $(F^{\lvert j\rvert,M})^{-1}_{q}$ denote the
inverse of the top left triangular $q\times q$ block in $F^{\lvert j\rvert,M}$
and let $\mathbf{a}^{j,M,q}$ be the vector composed of the first $q$ elements
in $\mathbf{a}^{j,M}$. Applying $(F^{\lvert j\rvert,M})^{-1}_{q}$ to
$\mathbf{a}^{j,M,q}$ is equivalent to running the algorithm (3.3) until
$k+1=q$ for the considered $j\in\mathbb{Z}_{M-1}$, with
$q\in\\{0,\dots,M-\lvert j\rvert\\}$. Note that $F^{M}$ contains $M(M+1)/2$
potentially distinct diagonal elements, and let

$\upsilon:\\{1,2,\dots,M(M+1)/2\\}\to\mathbb{N}_{0}^{M}$

order the absolute values of those elements in the same way that $\tau$ orders
singular values: $\upsilon_{l}(p)$ indicates how many of the in-absolute-value
largest $p$ diagonal elements of $F^{M}$ are found on the diagonal of
$F^{l,M}$. Note that $\upsilon$ can be formed and tabulated offline without
any prior knowledge apart from some upper bound on $M$.

Given $p\in\\{1,2,\dots,M(M+1)/2\\}$, the truncated triangular solution to
(4.1) is defined via

$\mathbf{c}^{j,M,p}=\begin{bmatrix}(F^{\lvert
j\rvert,M})^{-1}_{\upsilon_{\lvert
j\rvert}(p)}\mathbf{a}^{j,M,\upsilon_{\lvert j\rvert}(p)}\\\\[5.69054pt]
\mathbf{0}\end{bmatrix},\qquad j\in\mathbb{Z}_{M-1},$ (4.6)

where $\mathbf{0}$ is a zero-vector of length $M-\lvert
j\rvert-\upsilon_{\lvert j\rvert}(p)$. The complete truncated triangular
solution $\mathbf{c}^{M,p}\in\mathbb{C}^{M^{2}}$ for (4.2) is then defined as
the concatenation of $\mathbf{c}^{j,M,p}$ for $j\in\mathbb{Z}_{M-1}$.

Assume then that the available data is a noisy version
$\mathbf{a}^{M}_{\delta}$ of $\mathbf{a}^{M}$, approximately satisfying (4.4),
and let $\mathbf{c}^{M,p}_{\delta}$ be a corresponding truncated triangular
solution. Truncated triangular solutions are combined with the Morozov
discrepancy principle by choosing the smallest truncation index
$p\in\\{1,2,\dots,M(M+1)/2\\}$ that satisfies

$\lVert
F^{M}\mathbf{c}^{M,p}_{\delta}-\mathbf{a}^{M}_{\delta}\rVert_{2}^{2}=\sum_{j\in\mathbb{Z}_{M-1}}\lVert\widetilde{\mathbf{a}}^{M,j,\upsilon_{\lvert
j\rvert}(p)}_{\delta}\rVert_{2}^{2}\leq\omega^{2}\delta^{2},$ (4.7)

where $\omega\geq 1$ is again a ‘fudge’ factor and
$\widetilde{\mathbf{a}}^{M,j,\upsilon_{\lvert j\rvert}(p)}_{\delta}$ denotes
the vector composed of the last $M-\lvert j\rvert-\upsilon_{\lvert
j\rvert}(p)$ elements in $\mathbf{a}^{M,j}_{\delta}$. The equality in (4.7)
follows from the use of the exact inverse matrix for the top left
$\upsilon_{\lvert j\rvert}(p)\times\upsilon_{\lvert j\rvert}(p)$ subsystem in
(4.6).

##  5\. Simulating data on piecewise smooth domains and with the CEM

In this section, we explain how data measured on more general planar domains
or via electrode measurements can be used as input for the regularized
algorithms introduced in the preceding section. Although the algorithms have
been designed for linearized measurements, all data used as their input in
this work are obtained from boundary potential measurements corresponding to
nonlinear forward models of EIT.

###  5.1. Other planar domains with CM

Assume the imaged object is modeled as a simply-connected domain
$\Omega\subset\mathbb{R}^{2}$ with a Lipschitz boundary that is composed of a
finite number of $C^{1,\alpha}$ smooth arcs. Let $\Psi:\Omega\to D$ be a
bijective conformal mapping whose existence is guaranteed by the Riemann
mapping theorem. Denote the inverse of $\Psi$ by $\Phi$ and their complex
derivatives by $\Psi^{\prime}$ and $\Phi^{\prime}$. The mappings $\Psi$ and
$\Phi$ have continuously differentiable extensions to the relative interiors
of the $C^{1,\alpha}$ arcs composing $\partial\Omega$ and their images under
$\Psi$, respectively [18]. The absolute values $\lvert\Psi^{\prime}\rvert$ and
$\lvert\Phi^{\prime}\rvert$ restricted onto $\partial\Omega$ and $\partial D$,
respectively, are the Jacobian determinants of the corresponding boundary
transformations.

On $\partial\Omega$ define

$\widetilde{f}_{m}\coloneqq(f_{m}\circ\Psi)\lvert\Psi^{\prime}\rvert,\qquad
m\in\mathbb{Z}^{\prime},$ (5.1)

with $f_{m}$ still denoting the $m$th Fourier basis function on $\partial D$.
By performing a change of variables on each of the $C^{1,\alpha}$ arcs
composing $\partial\Omega$, it follows that

$\int_{\partial\Omega}\lvert\widetilde{f}_{m}\rvert^{2}\,\mathrm{d}s=\int_{\partial\Omega}\lvert
f_{m}\circ\Psi\rvert^{2}\lvert\Psi^{\prime}\rvert^{2}\,\mathrm{d}s=\int_{\partial
D}\lvert f_{m}\rvert^{2}\lvert\Phi^{\prime}\rvert^{-1}\,\mathrm{d}s<\infty$

since the singularities of $\lvert\Phi^{\prime}\rvert^{-1}$ at the preimages
of the corners on $\partial\Omega$ are integrable [18, Theorem 3.2 & Exercise
3.4.1]. A similar argument also yields

$\int_{\partial\Omega}\widetilde{f}_{m}\,\mathrm{d}s=\int_{\partial
D}f_{m}\,\mathrm{d}s=0.$

To summarize, $\widetilde{f}_{m}\in L^{2}_{\diamond}(\partial\Omega)$ is a
proper boundary current on $\partial\Omega$ for all $m\in\mathbb{Z}^{\prime}$.

Assume that $\Omega$ is characterized by the conductivity
$\widetilde{\gamma}\in L^{\infty}_{+}(\Omega)$ and drive the current density
$\widetilde{f}_{m}$ through $\partial\Omega$. The resulting potential
$\widetilde{u}_{m}^{\widetilde{\gamma}}$ in $\Omega$ is the unique weak
solution of

$-\nabla\cdot(\widetilde{\gamma}\nabla\widetilde{u})=0\quad\text{in
}\Omega,\qquad\nu\cdot(\widetilde{\gamma}\nabla\widetilde{u})=\widetilde{f}_{m}\quad\text{on
}\partial\Omega$

in $H^{1}_{\diamond}(\Omega)$. According to [17, Lemma 4.1 & Remark 4.1], the
conformally mapped potential

$u^{\gamma}_{m}\coloneqq\widetilde{u}^{\widetilde{\gamma}}_{m}\circ\Phi$

is a weak solution in $D$ to (2.1) for $f=f_{m}$ and
$\gamma=\widetilde{\gamma}\circ\Phi$.

These observations can be interpreted as guidelines for simulating (or
measuring) relative boundary data on $\Omega$ with the aim to reconstruct
information on the conductivity perturbation
$\widetilde{\eta}\coloneqq\widetilde{\gamma}-1$:

  1. (i)

Apply the current patterns
$\\{\widetilde{f}_{m}\\}_{m\in\mathbb{Z}_{M}^{\prime}}$ from (5.1) on
$\partial\Omega$ for conductivities $\widetilde{\gamma}$ and $1$, and measure
the resulting relative boundary potentials

$\widetilde{g}_{m}^{\widetilde{\eta}}\coloneqq(\widetilde{u}^{\widetilde{\gamma}}_{m}-\widetilde{u}^{1}_{m})|_{\partial\Omega},\qquad
m\in\mathbb{Z}_{M}^{\prime}.$

  2. (ii)

Form

$g^{\eta}_{m}\coloneqq\widetilde{g}_{m}^{\widetilde{\eta}}\circ\Phi,\qquad
m\in\mathbb{Z}_{M}^{\prime},$

that are relative boundary potentials on $\partial D$ for the Fourier current
basis $\\{f_{m}\\}_{m\in\mathbb{Z}_{M}^{\prime}}$, corresponding to the
conductivity perturbation $\eta=\widetilde{\gamma}\circ\Phi-1=\gamma-1$ in
$D$.

  3. (iii)

Form the nonlinear data matrix

$a_{m,n}^{\rm NL}(\eta)\coloneqq\langle
g^{\eta}_{m},f_{n}\rangle_{L^{2}(\partial D)},\qquad
m,n\in\mathbb{Z}_{M}^{\prime}.$ (5.2)

  4. (iv)

Use the matrix $\\{a_{m,n}^{\rm NL}(\eta)\\}_{m,n\in\mathbb{Z}_{M}^{\prime}}$
in place of $\\{a_{m,n}(\eta)\\}_{m,n\in\mathbb{Z}_{M}^{\prime}}$ as the input
for the algorithms of Section 4 to reconstruct an approximation for $\eta$.

  5. (v)

Form the final reconstruction of $\widetilde{\eta}$ in $\Omega$ by composing
the reconstruction in $D$ with $\Psi$.

All employed CM data for the unit disk, i.e., $\Omega=D$, are simulated by
choosing $\Phi={\rm Id}$ in the above scheme. Hence, also with $D$ as the
examined domain, the data used as the input for the reconstruction algorithms
of Section 4 correspond to a nonlinear forward model in our numerical tests.

###  5.2. CEM measurements

Assume next that the measurements of current and voltage on the boundary of
the simply-connected domain $\Omega$ are performed with the help of
$2L+1\in\mathbb{N}\setminus\\{1,2\\}$ contact electrodes
$\\{E_{l}\\}_{l\in\mathbb{Z}_{L}}$ that are modeled as open nonempty connected
subsets of $\partial\Omega$ with mutually disjoint closures. That is, one can
only control the mean-free net current pattern $I\in\mathbb{C}^{2L+1}$ through
the electrodes and measure the resulting electrode potential vector
$U\in\mathbb{C}^{2L+1}$. For simplicity, let us assume that all electrodes
have the same arclength $d>0$. The thin resistive layers at the electrode
object interfaces are modeled by contact impedances
$\\{z_{l}\\}_{l\in\mathbb{Z}_{L}}\subset\mathbb{C}$ that are assumed to have
strictly positive real parts and to stay the same between measurements. The
electrode measurements are modeled by the CEM [5, 20], which is in good
agreement with practical EIT.

As in Section 5.1, suppose that we aim to reconstruct a conductivity
perturbation $\widetilde{\eta}=\widetilde{\gamma}-1\in L^{\infty}(\Omega)$ by
performing electrode measurements on $\partial\Omega$ with
$\widetilde{\gamma}$ and $1$ acting in turns as the internal conductivity of
$\Omega$. Assuming one is able to use a linearly independent set
$\\{I_{l}\\}_{l\in\mathbb{Z}_{L}^{\prime}}$ of $2L$ mean-free current patterns
and freely choose the positions of the available electrodes, the construction
presented in [8, Section 6.2] produces an approximation
$\widehat{g}_{m}^{\widetilde{\eta}}$ based on the measured electrode
potentials for the relative boundary potential
$\widetilde{g}_{m}^{\widetilde{\eta}}$ considered in Section 5.1. If it is
further assumed that $\partial\Omega$ is smooth and
$\mathop{\textup{supp}}\widetilde{\eta}$ is compactly contained in $\Omega$,
then [8, Corollary 6.2] implies

$\displaystyle\lVert\widetilde{g}_{m}^{\widetilde{\eta}}-\widehat{g}_{m}^{\widetilde{\eta}}\rVert_{L^{2}(\partial\Omega)}$
$\displaystyle\leq
C\bigl{(}\bigl{(}\tfrac{1}{2}+L\bigr{)}^{-s}+d^{2}\bigr{)}\lVert\widetilde{f}_{m}\rVert_{H^{s}(\partial\Omega)}$
$\displaystyle\leq
C\bigl{(}\bigl{(}\tfrac{1}{2}+L\bigr{)}^{-s}+d^{2}\bigr{)}\lVert
f_{m}\rVert_{H^{s}(\partial D)}$ $\displaystyle\leq
C\bigl{(}\bigl{(}\tfrac{1}{2}+L\bigr{)}^{-s}+d^{2}\bigr{)}\lvert m\rvert^{s}$
(5.3)

for any $s>1/2$, with $C>0$ being a generic constant (not depending on $m$)
that may change between occurrences. The penultimate step follows from
$\Psi|_{\partial\Omega}$ being smooth due to the smoothness of
$\partial\Omega$ [18].

When the reconstruction algorithms of Section 4 are applied to electrode
measurements, we replace $\widetilde{g}_{m}^{\widetilde{\eta}}$ by
$\widehat{g}_{m}^{\widetilde{\eta}}$ in the guidelines of Section 5.1.
According to (5.2), the approximation provided by
$\widehat{g}_{m}^{\widetilde{\eta}}$ is good for low Fourier frequencies $m$
and small electrodes. In particular, the parameter $M\in\mathbb{N}$
controlling the size of the Fourier data matrix should be chosen smaller than
$L$ for (5.2) to guarantee a reasonable correspondence between
$\widehat{g}_{m}^{\widetilde{\eta}}$ and
$\widetilde{g}_{m}^{\widetilde{\eta}}$ for all involved Fourier frequencies
$m\in\mathbb{Z}_{M}^{\prime}$.

######  Remark 5.1.

In our numerical experiments with electrode measurements, the imaged object is
the unit disk with equiangular electrodes for simulated data and a water tank
with the shape of a right circular cylinder for real-world measurements [13].
The considered conductivity distributions inside the water tank are
homogeneous in the vertical direction, as are the rectangular electrodes that
are equiangularly attached to lateral surface of the tank and extend from its
bottom to the water surface. Due to certain symmetries, such three-dimensional
water tank measurements can be modeled with the two-dimensional CEM in the
unit disk, cf., e.g., [15].

Hence, one only needs to consider approximating the relative boundary
potential $g_{m}^{\eta}=\widetilde{g}_{m}^{\widetilde{\eta}}$ from
measurements on equiangular electrodes on the boundary of the unit disk. This
is achieved by introducing a mean-free current pattern $I_{m}$ via evaluating
a suitably scaled version of $f_{m}$ at the midpoints of the electrodes,
measuring the resulting relative electrode potentials corresponding to the
conductivity perturbation $\eta$, and finally performing trigonometric
interpolation of the relative potentials with respect to the midpoints of the
electrodes. See [8, Sections 5 and 6.2] for the details.

##  6\. Numerical experiments

We employ three numerical techniques for simulating the relative potentials
$g_{m}^{\eta}$, $m\in\mathbb{Z}_{M}^{\prime}$ on the boundary of the unit disk
$D$ corresponding to a given conformally mapped perturbation
$\widetilde{\eta}\in L^{\infty}(\Omega)$, with $\widetilde{\eta}>-1$, in the
imaged domain $\Omega$. Each of them generates approximate point values of
$g_{m}$ on an equidistant grid with $2N+1$ nodes on $\partial D$.
Subsequently, the nonlinear data matrix $\\{a_{m,n}^{\rm
NL}(\eta)\\}_{m,n\in\mathbb{Z}_{M}^{\prime}}$ is computed via fast Fourier
transform (FFT), which corresponds to approximating the integrals in (5.2) by
the periodic trapezoidal rule that has favorable convergence properties. The
values of $M$ and $N$ vary between different tests.

If $\widetilde{\eta}=\eta$ is a scalar multiple of the characteristic function
of a disk compactly contained in $\Omega=D$, the relative data $g_{m}^{\eta}$
is simulated by exploiting Möbius transformations. This technique is based on
$\\{f_{m}\\}_{m\in\mathbb{Z}^{\prime}}$ being the eigenfunctions with
explicitly known eigenvalues for
$\Lambda(1+\kappa\chi_{D_{\rho}})-\Lambda(1)$, where $\kappa>-1$ is a constant
and $\chi_{D_{\rho}}$ is the characteristic function of the open disk
$D_{\rho}$ of radius $0<\rho<1$ centered at the origin. See [8, Example 7.2]
and [16, Example 4] for more details. This approach leads to highly accurate
data, enabling studying the limits of the reconstruction method with noiseless
data as well as the effect of the linearization error. Alternatively, one
could use exact characterizations such as those in [7, Section 3].

For other types of $\widetilde{\eta}$ and/or domains $\Omega$, the absolute
potentials on $\partial\Omega$ for the current patterns
$\\{\widetilde{f}_{m}\\}_{m\in\mathbb{Z}_{M}^{\prime}}$ corresponding to the
conductivities $1+\widetilde{\eta}$ and $1$ are simulated by a finite element
method (FEM) with piecewise linear basis functions using the Python package
scikit-fem [11]. If $\Omega=D$, the finite element (FE) mesh is constructed so
that the boundary nodes directly form an equidistant grid, which makes
collecting the approximate nodal values of the relative potentials
$g_{m}^{\eta}$ straightforward. If $\Omega\not=D$ is a simply-connected
polygon, the Schwarz–Christoffel toolbox [6] is utilized to numerically
introduce the needed Riemann mapping $\Psi$. The boundary nodes of the FE mesh
in $\Omega$ are then chosen as the images of an equidistant grid on $\partial
D$ under $\Phi=\Psi^{-1}$. The simulated point values of the relative boundary
potentials $\widetilde{g}^{\widetilde{\eta}}_{m}$,
$m\in\mathbb{Z}_{M}^{\prime}$, at the boundary nodes of the FE mesh are then
precisely the values of $g_{m}^{\eta}$, $m\in\mathbb{Z}_{M}^{\prime}$, on the
equidistant grid of $\partial D$, cf. Section 5.1.

When considering electrode measurements on $\Omega=D$, the CEM forward
problems for the conductivities $1+\eta$ and $1$ are numerically solved with
FEM based on piecewise linear basis functions for the potential in $\Omega$
and certain macro elements for the electrode voltages [23]. Constructing an
approximation $\widehat{g}_{m}^{\eta}$ for the relative boundary potential
$\widetilde{g}_{m}^{\widetilde{\eta}}=g_{m}^{\eta}$ is performed as explained
in Remark 5.1, with the midpoints of the electrodes serving as the uniform
grid on $\partial D$, resulting in $N=L$. Note that only accounting for the
point values of $\widehat{g}_{m}^{\eta}$ at the electrode centers when
evaluating the integrals in (5.2) does not lead to a reduction in accuracy,
since the periodic trapezoidal rule is exact when computing (5.2) with
$g_{m}^{\eta}$ replaced by $\widehat{g}_{m}^{\eta}$; see [19] and [8, Eq.
(4.8)].

When considering noisy simulated measurements, we add independent realizations
of Gaussian noise to the real and imaginary parts of all elements in the
nonlinear data matrix $\\{a_{m,n}^{\rm
NL}(\eta)\\}_{m,n\in\mathbb{Z}_{M}^{\prime}}$. We say that the data contains
$100\sigma\%$ of additive noise, with $\sigma\in[0,1]$, if the noise added to
$\mathop{\textup{Re}}(a_{m,n}^{\rm NL})$ is drawn from
$\mathcal{N}(0,\sigma^{2}(\mathop{\textup{Re}}(a_{m,n}))^{2})$ and that added
to $\mathop{\textup{Im}}(a_{m,n}^{\rm NL})$ is drawn from
$\mathcal{N}(0,\sigma^{2}(\mathop{\textup{Im}}(a_{m,n}))^{2})$ for all
$m,n\in\mathbb{Z}_{M}^{\prime}$. For simulated measurements, we assume to know
the statistics of the additive noise and choose $\delta$ in (4.5) and (4.7) to
be the square root of the expected value for the squared Euclidean norm of the
noise vector, i.e., the square root of the sum of the variances of the
independent noise components. On the other hand, the value for the fudge
factor $\omega$ in (4.5) and (4.7) is case-dependent.

###  6.1. Visualization of right singular vectors of $F^{M}$

We start the numerical experiments by computing the right singular vectors of
$F^{M}$ for $M=32$, interpreting them as coefficients in a truncated Zernike
polynomial expansion via (4.2), (4.1), (3.8) and (3.2), and visualizing the
resulting right singular functions over $D$. We do not prove that the singular
functions or values of $F^{M}$ approximate such for the untruncated operator
$F$. We also remark that the order in which singular functions with certain
qualitative properties appear in the SVD depends on $M$.

Due to the block diagonal structure (4.3), each right singular function of
$F^{M}$ is a linear combination of Zernike polynomials with a fixed absolute
angular index $\lvert j\rvert$. In fact, there are $M$ singular functions for
$j=0$ and $2(M-\lvert j\rvert)$ for each $\lvert j\rvert\in\\{1,\dots,M-1\\}$.
The first right singular function is a linear combination of Zernike
polynomials for $j=0$ over different radial indices $k$. The second largest
singular value has multiplicity two, and the corresponding two orthonormal
singular functions can be chosen to be linear combinations of Zernike
polynomials for $j=1$ and $j=-1$, respectively, with the same coefficients for
each $k$ in (3.2). This is the convention we choose for any fixed $\lvert
j\rvert$ in what follows. The next $13$ singular values also have multiplicity
two, and they correspond to $\lvert j\rvert=2,\dots,14$. Tables 6.1 and 6.2,
respectively, show the real and imaginary parts of the first nine right
singular functions that exhibit higher angular frequencies and faster decay
toward the origin as $\lvert j\rvert$ increases.

$j=0$ | $j=\pm 1$ | $j=\pm 2$ | $j=\pm 3$ | $j=\pm 4$
---|---|---|---|---
|  |  |  |
|  |  |  |
Table 6.1. The real parts of the first nine right singular functions of $F^{M}$ for $M=32$. $j=0$ | $j=\pm 1$ | $j=\pm 2$ | $j=\pm 3$ | $j=\pm 4$
---|---|---|---|---
|  |  |  |
|  |  |  |
Table 6.2. The imaginary parts of the first nine right singular functions of
$F^{M}$ for $M=32$.

At a certain point, singular functions that are new linear combinations of
Zernike polynomials over $k$ for an already visited $\lvert j\rvert$ start to
appear. For our choice of $M=32$, the first such singular function is the 30th
that corresponds to the lowest angular frequency $j=0$. The next ones are the
45th and 46th, 53rd and 54th, as well as 59th and 60th singular functions that
are associated to $j=\pm 1$, $j=\pm 2$ and $j=\pm 3$, respectively. Their real
and imaginary parts, illustrated in Tables 6.3 and 6.4, exhibit similar
angular oscillations as the first singular functions for the same $\lvert
j\rvert$ in Tables 6.1 and 6.2. However, these ‘second order’ singular
functions in Tables 6.3 and 6.4 carry an additional radial oscillation close
to the boundary of $D$.

30th | 45th and 46th | 53rd and 54th | 59th and 60th
---|---|---|---
singular function, | singular functions, | singular functions, | singular functions,
$j=0$ | $j=\pm 1$ | $j=\pm 2$ | $j=\pm 3$
|  |  |
|  |  |
Table 6.3. The real parts of the first seven right singular functions of $F^{M}$ that are linear combinations of Zernike polynomials for some already visited angular index $j$ with $M=32$. 30th | 45th and 46th | 53rd and 54th | 59th and 60th
---|---|---|---
singular function, | singular functions, | singular functions, | singular functions,
$j=0$ | $j=\pm 1$ | $j=\pm 2$ | $j=\pm 3$
|  |  |
|  |  |
Table 6.4. The imaginary parts of the first seven right singular functions of
$F^{M}$ that are linear combinations of Zernike polynomials for some already
visited angular index $j$ with $M=32$.

As expected, the first ‘third order’ singular function corresponds to $j=0$,
and it is altogether the 67th one. The next angular frequencies $j=\pm 1$ are
visited for the third time at the 74th and 75th singular functions. As a
general rule, the smaller the singular value corresponding to a given absolute
angular index $\lvert j\rvert$ is, the more oscillations there are in the
radial direction in the associated singular functions. To illustrate this, the
left-hand image of Fig. 6.5 plots the real parts of the first three singular
functions (2nd, 45th, 74th) for $j=1$ along the real axis. Finally, the right-
hand image of Fig. 6.5 presents, for the sake of completeness, the first 100
singular values of $F^{M}$ counted according to their multiplicity.

(a) The real parts of the first three right singular functions for $j=1$ along
the horizontal axis.

(b) The first 100 singular values with the ones corresponding to the singular
functions in the left-hand image highlighted.

Figure 6.5. The radial components of certain right singular functions and the
first hundred singular values of $F^{M}$ for $M=32$.

###  6.2. CM and a discoidal inclusion in $D$

The first numerical tests on the actual reconstruction algorithms consider the
target in Fig. 6(a), corresponding to a discoidal inclusion of radius $0.2$
centered at $(\tfrac{1}{4},\tfrac{\sqrt{3}}{4})$ inside the unit disk. The
conductivity of the inclusion is $1.2$ and, as in all our tests based on
simulated measurements, that of the background is $1$. We visualize all
targets and reconstructions as perturbations of the background conductivity,
which explains the colormap of Fig. 6.6. Highly accurate data $\\{a_{m,n}^{\rm
NL}(\eta)\\}_{m,n\in\mathbb{Z}_{M}^{\prime}}$ is generated for $M=32$ using
Möbius transformations as briefly explained at the beginning of Section 6.

First, the truncated SVD algorithm described in Section 4.1 is run on the
noiseless data with minimal truncation. According to a subjective visual
inspection over many truncation indices, $p=492$ gives a close to optimal
reconstruction that is shown in Fig. 6(b) and demonstrates that the support of
a target inclusion can indeed be faithfully reconstructed despite the employed
linearization if the available data is very accurate, cf. [12]. This
conclusion remains valid even for inclusions with higher contrasts. Although
the reconstruction in Fig. 6(b) also reproduces the conductivity level of the
inclusion rather accurately, this would no longer be the case if the
conductivity of the inhomogeneity differed considerably more from the
background value, due to the linearization error.

Next both the truncated SVD method and the algorithm based on truncation of
triangular subsystems from Section 4.2 are applied to data with $1\%$ of
additive noise. We set $\omega=1$ in (4.5) and (4.7), which means that there
is no attempt to compensate for the employed linearization or numerical errors
via the choice of the fudge factor. The resulting reconstructions are shown in
Fig. 6(c) and Fig. 6(d), and they correspond to the truncation indices $p=174$
and $p=154$ in the respective algorithms. In particular, far fewer Zernike
coefficients are considered in Figs. 6(c) and 6(d) than in Fig. 6(b). The two
regularization methods seem to work about equally well, with the SVD-based
approach better recapturing the unperturbed background but being worse at
estimating the strength of the conductivity perturbation.

(a) The target conductivity perturbation.

(b) Very accurate data and minimal SVD-based regularization with $p=492$.

(c) $1\%$ of noise, SVD-based regularization with $\omega=1$, $p=174$.

(d) $1\%$ of noise, truncation of triangular subsystems with $\omega=1$,
$p=154$.

Figure 6.6. Reconstructions of a discoidal inclusion inside the unit disk from
highly accurate data generated using a Möbius transformation with $M=32$.

The Zernike coefficients $c_{j,k}$ that contribute to the reconstructions in
Figs. 6(b), 6(c) and 6(d), i.e., the coefficients that are accepted/estimated
by the corresponding regularization algorithms, are plotted in Fig. 6.7. For
both algorithms, the number of considered radial indices $k$ decreases as a
function of $|j|$, with this effect being more pronounced for the truncated
SVD method.

(a) The SVD-based reconstruction in Fig 6(b).

(b) The SVD-based reconstruction in Fig 6(c).

(c) Truncation of triangular subsystems in Fig 6(d).

Figure 6.7. The accepted Zernike indices for the reconstructions in Fig 6.6.

###  6.3. CM for other types of perturbations and domains

In this section the ‘nonlinear’ data matrix $\\{a_{m,n}^{\rm
NL}(\eta)\\}_{m,n\in\mathbb{Z}_{M}^{\prime}}$, with $M=16$, is simulated by
FEM using $1024$ boundary nodes. This leads to FE meshes with approximately
130,000 nodes and 260,000 triangles for all considered geometries.

(a) The target conductivity perturbation.

(b) $1\%$ of noise, truncated SVD with $\omega=1$, $p=151$.

(c) $1\%$ of noise, truncated triangular subsystems with $\omega=60$, $p=120$.

(d) $10\%$ of noise, truncated SVD with $\omega=1$, $p=114$.

(e) $10\%$ of noise, truncated triangular subsystems with $\omega=7$, $p=78$.

Figure 6.8. Reconstructions of a wave-like conductivity perturbation in the
unit disk with $M=16$.

Figure 8(a) illustrates the unit disk with a target conductivity perturbation
that behaves as a sine wave with amplitude 0.1 in the vertical direction and
is homogeneous in the horizontal direction. The other four images in Fig. 6.8
show reconstructions produced by the truncated SVD method (Section 4.1) and
truncation of triangular subsystems (Section 4.2) for the noise levels of 1%
and 10%. The SVD-based reconstructions still correspond to the trivial choice
$\omega=1$ in (4.5). However, for the truncation of triangular subsystems, the
fudge factor needs to be manually tuned in (4.7) to enable reasonable
reconstructions, leading to the heuristic choices $\omega=60$ and $\omega=7$,
respectively, for the lower and higher noise levels.

The reconstructions in Fig. 6.8 are able capture the qualitative structure of
the target conductivity for both noise levels, although the periodic pattern
is somewhat distorted in all reconstructions. This effect is most visible
close to the boundary and center of the domain, with reconstructions with the
higher noise level suffering from stronger artifacts in this regard. Be that
as it may, Fig. 6.8 demonstrates that ignoring the linearization error may
still lead to useful reconstructions in EIT even if the supports of the to-be-
reconstructed conductivity perturbations are not restricted to the interior of
the imaged domain.

A possible explanation for the requirement to tune the Morozov principle for
the algorithm based on the truncation of triangular subsystems is its
sensitivity to numerical inaccuracies in those elements of the data matrix
that correspond to high spatial frequencies: in Section 6.2 all simulated data
were highly accurate, whereas the FEM simulations of this section are bound to
suffer from numerical inaccuracies that are presumably the largest for the
measurements corresponding to the highest spatial frequencies in the set of
employed Fourier current patterns
$\\{f_{m}\\}_{m\in\mathbb{Z}_{16}^{\prime}}$. As all our remaining tests
consider either data simulated by FEM or/and CEM measurements, we only employ
the SVD-based regularization method in what follows.

We then test the truncated SVD algorithm on the two polygonal domains shown in
Figs. 6.9 and 6.10. The reason for only considering convex polygons as models
for the imaged object $\Omega$ is the following: The derivative of the
conformal mapping $\Psi^{\prime}|_{\partial\Omega}$, needed in (5.1) for
defining the appropriate current basis on $\partial\Omega$, has a zero at each
corner of $\partial\Omega$ with interior angle less than $\pi$ and a
singularity at each corner with interior angle wider than $\pi$; see, e.g.,
[21, Section 3.2] for more details. Hence, by settling for convex polygons,
one can simulate the required measurements on $\Omega$ without having to solve
Neumann problems for singular boundary current densities.

(a) The target conductivity perturbation.

(b) Reconstruction.

Figure 6.9. An SVD-based reconstruction of a smooth conductivity perturbation
in a square for $M=16$, $\omega=1$ and $1\%$ of additive noise, resulting in
$p=163$.

(a) The target conductivity perturbation.

(b) Reconstruction.

Figure 6.10. An SVD-based reconstruction of a smooth conductivity perturbation
in a polygon for $M=16$, $\omega=1$ and $1\%$ of additive noise, resulting in
$p=145$.

The considered conductivity perturbations are smooth Gaussian humps centered
so far from the boundary that their supports are compactly contained in
$\Omega$ within the numerical precision. The SVD-based algorithm, with
$\omega=1$, is applied to simulated data with 1% of additive noise to produce
the reconstructions on the right in Figs. 6.9 and 6.10. Although the algorithm
finds the positions of the Gaussian humps, the reconstructions are not as
localized as the target perturbations, and they also exhibit wave-like
artifacts in the background.

###  6.4. CEM data for $D$

Consider an electrode configuration with $32$ electrodes that are attached
equiangularly to the unit disk and cover, in total, half of its
boundary.222Although our model for the CEM in Section 5.2 assumes an odd
number of electrodes, the generalization of the approximation in [8] for an
even number of electrodes is straightforward. The contact impedance is set to
$z=1$ on all electrodes. CEM measurements are simulated for $31$ linearly
independent current patterns by a FEM with $33,281$ piecewise linear basis
functions and appropriate refinements of the mesh at the electrodes (cf.
[23]); this is done for both the mere background conductivity and the
background conductivity contaminated by the perturbation composed of two
Gaussian humps with opposite signs depicted in Fig. 11(a). As briefly
explained in Section 5.2, the ‘measured’ relative electrode potentials can
then be used to construct an approximation for the data matrix
$\\{a_{m,n}^{\rm NL}(\eta)\\}_{m,n\in\mathbb{Z}_{M}^{\prime}}$. See [8] for
the details.

We choose $M=16$ as the maximal absolute value of the considered spatial
Fourier frequencies and add $1\%$ of noise to the data. Observe that in
addition to the synthetic measurement noise, there are two other significant
sources of error in the data: the approximation between the CM and the CEM
quantified by (5.2) and the employment of ‘nonlinear’ data for a linear
reconstruction method. We employ the SVD-based algorithm with $\omega=1$ to
form the reconstruction presented in Fig. 11(b). The reconstruction captures
the basic characteristics of the target perturbation reasonably well.

(a) The target conductivity perturbation.

(b) Reconstruction.

Figure 6.11. An SVD-based reconstruction of a smooth conductivity perturbation
from simulated CEM data with 32 electrodes, $M=16$, $\omega=1$ and $1\%$ of
additive noise, resulting in $p=108$.

###  6.5. Water tank measurements

In the final experiment, we consider data from a right cylindrical water tank
with 16 electrodes documented in [13]. Note that all measurements in [13] are
performed with electrodes that are homogeneous in the vertical direction and
extend from the bottom of the tank exactly up to the water surface. The
embedded inclusions are also homogeneous in the vertical direction and break
the water surface, and they are either perfect insulators (plastic) or can be
modeled as ideal conductors (metal). Due to the insulating, i.e., homogeneous
Neumann, boundary conditions at the bottom and top of the water layer, the
measurements can thus be modeled by the two-dimensional CEM in the unit disk.
Moreover, conductivity perturbations reconstructed by resorting to such a
model could be presented in the SI units by appropriately scaling with the
dimensions of the tank and the conductivity of the saline filling the tank.
However, since it is not clear what the precise value for the conductivity of
the saline is in [13], the presented reconstructions are to be understood as
changes to the unit conductivity inside the unit disk without a proper
interpretation in SI units.

We employ the current pattern set 5 for the target cases 2.3, 4.4 and 5.2 as
numbered in [13]. The three targets are shown at top in Fig. 6.12, and the set
5, dubbed “all against 1” in [13], corresponds to feeding current through the
right-most electrode and letting it out in turns through the other ones,
resulting in 15 linearly independent current patterns. Instead of the set 5,
one could as well use any other current patterns in [13], as long as they form
a basis for the space of mean-free vectors of length 16. Note that the
measurements for the empty tank, i.e., the case 1.0 in [13], are also needed
for building an approximation for the data matrix $\\{a_{m,n}^{\rm
NL}(\eta)\\}_{m,n\in\mathbb{Z}_{M}^{\prime}}$, as reviewed in Section 6.4. In
particular, the applicability of the approximation (5.2) is based on an
assumption that the contacts at the electrodes are (approximately) the same
for all considered target cases.

Figure 6.12. Three SVD-based reconstructions with $M=8$ from real-world water
tank difference data measured at 16 equiangular electrodes. For each target,
the ‘optimal’ truncation index $p$ has been chosen according to a visual
inspection. Top: the targets. Bottom: the reconstructions.

We choose $M=8$ and apply the SVD-based reconstruction algorithm to the three
sets of data corresponding to the configurations on the top row of Fig. 6.12.
Since we do not know the noise level (or the noise model) exactly, we decide
to tune the truncation index $p$ via trial and error so that each of the three
reconstructions is as good as possible according to a subjective visual
inspection. The resulting reconstructions are shown on the bottom row of Fig.
6.12. According to our noise model and the Morozov discrepancy principle, the
used regularization parameters correspond to noise levels between $0.1\%$ and
$1\%$. Although the reconstructions in Fig. 6.12 are not based on any
systematic way of choosing the level of regularization, they in any case
demonstrate that the presented algorithm can under optimal circumstances
produce reasonable reconstructions in practical EIT.

##  7\. Concluding remarks

This paper introduced, implemented and numerically tested two regularized
variants of the direct reconstruction method for the linearized CM of two-
dimensional EIT presented in [10]. Both algorithms exploit the possibility to
present the linearized problem in the unit disk with the help of a block
diagonal (infinite) matrix with lower triangular blocks, which enables
tackling the triangular subsystems separately in the inversion process. More
precisely, one of the regularized algorithms is based on a semi-heuristic
truncation of the triangular subsystems and the other one on truncating their
SVDs. According to our numerical tests, the latter variant is more stable and
capable of producing high-quality reconstructions from accurate simulated
data, as well as reasonably good reconstructions from noisy measurements on
convex polygonal domains and even from real-world water tank data. In
particular, ignoring the linearization error does not seem to significantly
impede the practical applicability of the introduced reconstruction method.

### Acknowledgments

This work is supported by the Academy of Finland (decisions 353080, 353081, 359181) and the Aalto Science Institute (AScI). HG is supported by grant 10.46540/3120-00003B from Independent Research Fund Denmark | Natural Sciences.

## References

  * [1] A. Allers and F. Santosa.  Stability and resolution analysis of a linearized problem in electrical impedance tomography.  Inverse Problems, 15:515, 1991.
  * [2] L. Borcea.  Electrical impedance tomography.  Inverse problems, 18:R99–R136, 2002.
  * [3] A. P. Calderón.  On an inverse boundary value problem.  In Seminar on Numerical Analysis and its Applications to Continuum Physics, pages 65–73. Soc. Brasil. Mat., Rio de Janeiro, 1980\.
  * [4] M. Cheney, D. Isaacson, and J.C. Newell.  Electrical impedance tomography.  SIAM Rev., 41:85–101, 1999.
  * [5] K.-S. Cheng, D. Isaacson, J. S. Newell, and D. G. Gisser.  Electrode models for electric current computed tomography.  IEEE Trans. Biomed. Eng., 36:918–924, 1989.
  * [6] T. A. Driscoll.  Algorithm 756; a MATLAB toolbox for Schwarz–Christoffel mapping.  ACM Trans. Math. Soft., 22:168–186, 1996.
  * [7] H. Garde.  Comparison of linear and non-linear monotonicity-based shape reconstruction using exact matrix characterizations.  Inverse Probl. Sci. Eng., 26(1):33–50, 2018.
  * [8] H. Garde and N. Hyvönen.  Mimicking relative continuum measurements by electrode data in two-dimensional electrical impedance tomography.  Numer. Math., 147(3):579–609, 2021.
  * [9] H. Garde and N. Hyvönen.  Series reversion in Calderón’s problem.  Math. Comp., 91:1925–1953, 2022.
  * [10] H. Garde and N. Hyvönen.  Linearised Calderón problem: Reconstruction and Lipschitz stability for infinite-dimensional spaces of unbounded perturbations.  SIAM J. Math. Anal., 2024.  Accepted, arXiv:2204.10164.
  * [11] T. Gustafsson and G. D. McBain.  scikit-fem: A Python package for finite element assembly.  Journal of Open Source Software, 5(52):2369, 2020.
  * [12] B. Harrach and J. K. Seo.  Exact shape-reconstruction by one-step linearization in electrical impedance tomography.  SIAM J. Math. Anal., 42:1505–1518, 2010.
  * [13] A. Hauptmann, V. Kolehmainen, N. M. Mach, T. Savolainen, A. Seppänen, and S. Siltanen.  Open 2D electrical impedance tomography data archive, 2017.  arXiv:1704.01178.
  * [14] N. J. Higham.  Accuracy and stability of numerical algorithms.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, second edition, 2002.
  * [15] N. Hyvönen and L. Mustonen.  Smoothened complete electrode model.  SIAM J. Appl. Math., 77:2250–2271, 2017.
  * [16] N. Hyvönen and L. Mustonen.  Generalized linearization techniques in electrical impedance tomography.  Numer. Math., 140:95–120, 2018.
  * [17] N. Hyvönen, L. Päivärinta, and J. Tamminen.  Enhancing D-bar reconstructions for electrical impedance tomography with conformal maps.  Inverse Probl. Imaging, 12:373–400, 2017.
  * [18] Ch. Pommerenke.  Boundary behaviour of conformal maps.  Springer-Verlag, 1992.
  * [19] J. Saranen and G. Vainikko.  Periodic integral and pseudodifferential equations with numerical approximation.  Springer Monographs in Mathematics. Springer-Verlag, Berlin, 2002.
  * [20] E. Somersalo, M. Cheney, and D. Isaacson.  Existence and uniqueness for electrode models for electric current computed tomography.  SIAM J. Appl. Math., 52:1023–1040, 1992.
  * [21] L. N. Trefethen.  Numerical computation of the Schwarz–Christoffel transformation.  SIAM J. Sci. Stat. Comput., 1:82–102, 1980.
  * [22] G. Uhlmann.  Electrical impedance tomography and Calderón’s problem.  Inverse Problems, 25:123011, 2009.
  * [23] M. Vauhkonen.  Electrical impedance tomography with prior information, volume 62.  Kuopio University Publications C (Dissertation), 1997.
  * [24] F. Zernike.  Beugungstheorie des schneidenverfahrens und seiner verbesserten form, der phasenkontrastmethode.  Physica, 1:689–704, 1934.