Electrical Impedance Tomography: A Fair Comparative Study on Deep Learning and Analytic-based Approaches

This section introduces the mathematical model of the EIT problem and the discrepancy functional used for reconstructing the conductivity $\sigma$. Let $\Omega$ be an open-bounded domain in $\mathbb{R}^{d}\ (d\geq 2)$ with a Lipschitz boundary $\partial\Omega$, and let $\Lambda_{\sigma,N}$ denote the NtD map of problem (1). We employ the usual Sobolev space for the Neumann boundary data $\sigma\frac{\partial u}{\partial n}=j\in\tilde{H}^{-\frac{1}{2}}(\partial\Omega)$, respectively Dirichlet boundary condition $u=\phi\in\tilde{H}^{\frac{1}{2}}(\partial\Omega)$ on $\partial\Omega$. Throughout, we make use of the space $\tilde{H}^{1}(\Omega)$, which is a subspace of the Sobolev space $H^{1}(\Omega)$ with vanishing mean on $\partial\Omega$, i.e., $\tilde{H}^{1}(\Omega)=\{v\in H^{1}(\Omega):\int_{\partial\Omega}v{\rm d}s=0\}$. The spaces $\tilde{H}^{\frac{1}{2}}(\partial\Omega)$ and $\tilde{H}^{-\frac{1}{2}}(\partial\Omega)$ are defined similarly. These spaces are equipped with the usual norms. We normalise the solution of the Neumann problem by enforcing $\int_{\partial\Omega}u{\rm d}s=0$, so that there exists a unique solution $u\in\tilde{H}^{1}(\Omega)$. We denote the Dirichlet-to-Neumann (DtN) map by $\Lambda_{\sigma,D}$. Then we have $\Lambda_{\sigma,N}=\Lambda_{\sigma,D}^{-1}$, i.e., DtN and NtD maps are inverse to each other. In usual regularised reconstruction, we employ the NtD map $\Lambda_{\sigma,N}$, whereas in the D-bar method, we employ the DtN map $\Lambda_{\sigma,D}$.

An EIT experiment consists of applying a current $j$ and measuring the resulting potential $\phi$ on $\partial\Omega$, and it is equivalent to solving a Neumann forward problem with the physical conductivity $\sigma^{\dagger}$, i.e. $\phi=\Lambda_{\sigma,N}j$, on $\partial\Omega$. In practice, the boundary potential measurements are collected experimentally, and thus $\phi$ is only an element of the space $L^{2}(\Gamma)$. see e.g. [22]. Note that the continuum model is mostly academic. A more realistic model is the so-called complete electrode model (CEM) for EIT [100, 53], which models contact impedances and localised electrode geometries. The CEM is finite-dimensional by construction, leading to different mathematical challenges and reconstruction methods.

The solvability, uniqueness and smoothness of the continuum model with respect to $L^{p}$ norms can be derived using Meyers’ gradient estimate [80], as in [92].

In Theorem 2.1, the boundary condition for the problem can be general. Its effect enters the $W^{1,q}$-estimate through the term $\|u\|_{H^{1}(\Omega)}$. In addition, no regularity has been assumed on $\sigma$. Generally, a precise estimate of the constant $Q(\lambda,d)$ is missing, but in the 2D case, a fairly sharp estimate of $Q(\lambda,d)$ was derived in [6].

EIT suffers from a high degree of non-linearity and severe ill-posedness, as typical of many PDE inverse problems with boundary data. However, its potential applications have sparked much interest in designing effective numerical techniques for its efficient solution. Numerous numerical methods have been proposed in the literature; see [14, Section 7] for an overview (up to 2002). These methods can roughly be divided into two groups: regularised reconstruction and direct methods. Below, we give a brief categorisation of conventional reconstruction schemes.

The methods in the first group are of variational type, i.e., based on minimising a certain discrepancy functional. Commonly the discrepancy $J$ is the standard least-squares fitting, i.e., the squared $L^{2}(\partial\Omega)$ norm of the difference between the electrical potential due to the applied current $j$ and the measured potential $\phi$:

for one single measurement $(j,\phi^{\delta})$. One early approach of this type is given in [28], which applies one step of a Newton method with a constant conductivity as the initial guess. Due to the severe ill-posedness of the problem, regularisation is beneficial for obtaining reconstructions with improved resolution [35, 94, 56]. Commonly used penalties include Sobolev smoothness [78, 58] for a smooth conductivity distribution, total variation [50], Mumford-Shah functional [92], level set method [30] for recovering piecewise constant conductivity, sparsity [40, 60, 59] for recovering small inclusions (relative to the background). The combined functional is given by

where $R(\sigma)$ denotes the penalty, and $\alpha>0$ is the penalty weight. The functional $\Psi(\sigma)$ is then minimised over the admissible set

for some $\lambda\in(0,1)$. The set $\mathcal{A}$ is usually equipped with an $L^{p}(\Omega)$ norm $(1\leq p\leq\infty)$. One may also employ data fitting other than the standard $L^{2}(\partial\Omega)$-norm. The most noteworthy one is the Kohn-Vogelius approach, which lifts the boundary data to the domain $\Omega$ and makes the fitting in $\Omega$ [107, 69, 15]; see also [67] for a variant of the Kohn-Vogelius functional. In practice, the regularized formulations have to be properly discretized, commonly done by means of finite element methods [39, 91, 62, 61], due to the spatially variable conductivity and irregular domain geometry. Newton-type methods have also been applied to EIT [71, 72]. Probabilistic formulations of these deterministic approaches are also possible [64, 38, 34, 12], which can provide uncertainty estimates on the reconstruction.

The methods in the second group are of a more direct nature, aiming at extracting relevant information from the given data directly, without going through the expensive iterative process. Bruhl et al [18, 19] developed the factorisation method for EIT, which provides a criterion for determining whether a point lies inside or outside the set of inclusions by carefully analysing the spectral properties of certain operators. Thus, the inclusions can be reconstructed directly by testing every point in the computational domain. The D-bar method of Siltanen, Mueller and Isaacson [98, 82] is based on Nachman’s uniqueness proof [83] and utilises the complex geometric solutions and nonphysical scattering transform for direct image reconstruction. Chow, Ito and Zou [29] proposed the direct sampling method when there are only very few Cauchy data pairs. The method employs dipole potential as the probing function and constructs an indicator function for imaging the inclusions in EIT, and it is easy to implement and computationally cheap. Other notable methods in the group include monotonicity method [47], enclosure method [54], Calderón’s method [11, 96], and MUSIC [3, 2, 70] among others. Generally, direct methods are faster than those based on variational regularisation, but the reconstructions are often inferior in terms of resolution and can suffer from severe blurring.

These represent the most common model-based inversion techniques for EIT reconstruction. Despite these important progress and developments, the quality of images produced by EIT remains modest when compared with other imaging modalities. In particular, at present, EIT reconstruction algorithms are still unable to extract sufficiently useful information from data to be an established routine procedure in many medical applications. Moreover, the iterative schemes are generally time-consuming, especially for 3D problems. One possible way of improving the quality of information is to develop an increased focus on identifying useful information and fully exploiting a priori knowledge. This idea has been applied many times, and the recent advent of deep learning significantly expanded its horizon from hand-crafted regularisers to more complex and realistic learned schemes. Indeed, recently, deep learning-based approaches have been developed to address these challenges by drawing on knowledge encoded in the dataset or structural preference of the neural network architecture.

We describe the sparsity approach and D-bar method next, and deep learning approaches in Section 3.

The sparsity concept is very useful for modelling conductivity distributions with “simple” descriptions away from the known background $\sigma_{0}$, e.g. when $\sigma$ consists of an uninteresting background plus some small inclusions. Let $\delta\sigma^{\dagger}=\sigma^{\dagger}-\sigma_{0}$. A “simple” description means that $\delta\sigma$ has a sparse representation with respect to a certain basis/frame/dictionary $\{\psi_{k}\}$, i.e., there are only a few non-zero expansion coefficients. The $\ell^{1}$ norm $\delta\sigma$ can promote the sparsity of $\delta\sigma$ [33]

Under certain regularity conditions on $\{\psi_{k}\}$, the problem of minimising $\Psi$ over the set $\mathcal{A}$ is well-posed [58].

Optimisation problems with the $\ell^{1}$ penalty have attracted intensive interest [33, 13, 17, 108]. The challenge lies in the non-smoothness of the $\ell^{1}$-penalty and high-degree nonlinearity of the discrepancy $J(\sigma)$. The basic algorithm for updating the increment $\delta\sigma_{i}$ and $\sigma_{i}=\sigma_{0}+\delta\sigma_{i}$ by minimising $\Psi$ formally reads

where $s>0$ is the step size, $\Lambda_{\sigma_{i},N}^{\prime}$ denotes the Gâteaux derivative of the NtD map $\Lambda_{\sigma_{i},N}$ in $\sigma$, and $\mathcal{S}_{\lambda}(t)=\mbox{sign}(t)\max(|t|-\lambda,0)$ is the soft shrinkage operator. However, a direct application of the algorithm does not yield accurate results. We adopt the procedure in Algorithm 1. The key tasks include computing the gradient $J^{\prime}$ (Steps 4-5) and selecting the step size (Step 6).

The above description follows closely the work [59], where the sparsity algorithm was first developed. There are alternative sparse reconstruction techniques, notably based on total variation [92, 39, 16, 112]. For example, [16] presented an experimental (in-vivo) evaluation of the total variation approach using a linearized model, and the resulting optimisation problem solved by the primal-dual interior point method; and the work [112] compared different optimisers. Due to the non-smoothness of the total variation, one may relax the formulation with the Modica-Mortola function in the sense of Gamma convergence [92, 61].

The D-bar method of Siltanen, Mueller and Isaacson [98] is a direct reconstruction algorithm based on the uniqueness proof due to Nachman [83]; see also Novikov [87]. That is, a reconstruction is directly obtained from the DtN map $\Lambda_{\sigma,D}$, without going through an iterative process. Note that the DtN map $\Lambda_{\sigma,D}$ can be computed as the inverse of the measured NtD map $\Lambda_{\sigma,N}$ when full boundary data is available. Below we briefly overview the classic D-bar algorithm assuming $\sigma\in C^{2}(\Omega)$, with a positive lower bound (i.e., $\sigma\geq c>0$ in $\Omega$), and $\sigma\equiv 1$ in a neighbourhood of the boundary $\partial\Omega$. In this part, we consider an embedding of $\mathbb{R}^{2}$ in the complex plane, and hence we will identify planar points $x=(x_{1},x_{2})$ with the corresponding complex number $x_{1}+{\rm i}x_{2}$, and the product $kx$ denotes complex multiplication. For more detailed discussions, we refer interested readers to the survey [82].

First, we transform the conductivity equation (1) into a Schrödinger-type equation by substituting $\widetilde{u}=\sqrt{\sigma}u$ and setting $q=\Delta\sqrt{\sigma}/\sqrt{\sigma}$ and extending $\sigma\equiv 1$ outside $\Omega$. Then we obtain

Next we introduce a class of special solutions of equation (4) due to Faddeev [36], the so-called complex geometrical optics (CGO) solutions $\psi(x,k)$, depending on a complex parameter $k\in\mathbb{C}\setminus\{0\}$ and $x\in\mathbb{R}^{2}$. These exponentially behaving functions are key to the reconstruction. Specifically, given $q\in L^{p}(\mathbb{R}^{2}),\,1<p<2$, the CGO solutions $\psi(x,k)$ are defined as solutions to

satisfying the asymptotic condition $e^{-{\rm i}kx}\psi(x,k)-1\in W^{1,\tilde{p}}(\mathbb{R}^{2})$ with $2<\tilde{p}<\infty$. These solutions are unique for $k\in\mathbb{C}\setminus\{0\}$ as shown in [83, Theorem 1.1]. Then D-bar algorithm recovers the conductivity $\sigma$ from the knowledge of the CGO solutions $\mu(x,k)=e^{-{\rm i}kx}\psi(x,k)$ at the limit $k\to 0$ [83, Section 3]

Numerically, one can substitute the limit by $k=0$ and evaluate $\mu(x,0)$. The reconstruction of $\sigma$ relies on the use of an intermediate object called non-physical scattering transform $\mathbf{t}$, defined by

with $e_{k}(x):=\exp({\rm i}(kx+\bar{k}\bar{x}))$, where over-bar denotes complex conjugate. Since $\mu$ is asymptotically close to one, $\mathbf{t}(k)$ is similar to the Fourier transform of $q(x)$. Meanwhile, we can obtain $\mu$ by solving the name-giving D-bar equation

where $\bar{\partial}_{k}=\frac{1}{2}(\frac{\partial}{\partial k_{1}}+{\rm i}\frac{\partial}{\partial k_{2}})$ is known as the D-bar operator. To solve the above equation, scattering transform $\mathbf{t}(k)$ is required, which we can not measure directly from the experiment, but $\mathbf{t}(k)$ can be represented using the DtN map. Indeed, using Alessandrini’s identity [1], we get the boundary integral

Note that $\Lambda_{1,D}$ can be analytically computed, and only $\Lambda_{\sigma,D}$ needs to be obtained from the measurements. Here, we will employ a Born approximation using $\psi\approx e^{{\rm i}kx}$, leading to the linearised approximation

This linearised D-bar algorithm can be efficiently implemented. First, one computes the $\mathbf{t}^{\exp}(k)$ from the measured DtN map $\Lambda_{\sigma,D}$, and then one solves the D-bar equation (5). Note that the solutions of (5) are independent for each $x\in\Omega$ and one can efficiently parallelise over $x$. This leads to real-time implementations and is especially relevant for time-critical applications, e.g., monitoring purposes. The fully nonlinear D-bar algorithm would require first computing $\psi$ by solving a boundary integral equation and then computing the scattering transform $\mathbf{t}(k)$.

The above algorithm assumes infinite precision and noise-free data. When the data is noise corrupted with finite measurements, the measured DtN map $\Lambda_{\sigma,D}$ is not accurate, and then the computation of $\mathbf{t}(k)$ becomes exponentially unstable for $|k|>R$. Thus, for practical data, we need to restrict the computations to a certain frequency range so as to stably compute $\mathbf{t}(k)$. Below we choose $R=5$ for noise-free data and $R=4.5,\,R=4$ for 1% and 5% noisy measurements, respectively. This strategy of reducing the cut-off radius for noisy measurements is shown to be a regularisation strategy [68]. The final algorithm can be summarised as outlined below in Algorithm 2.

Besides the D-bar method, there are other analytic and direct reconstruction methods available, e.g., enclosure method [54], monotonicity method [47], direct sampling method [29], and Calderón’s method [11, 96]. The common advantage of these approaches is their computational efficiency, but unfortunately, also the directly inherited exponential instability to noise. While there are strategies to deal with noise, e.g., reducing the cut-off radius, the reconstruction quality does suffer: the reconstructions tend to be overall smooth. Additionally, there may be theoretical limitations to the reconstructions that can be obtained. For example, for the classic D-bar algorithm, it is $C^{2}$ conductivities, and for the enclosure methods, we can only find the convex hull of all inclusions. Thus, it is very interesting to discuss how deep learning can help overcome these limitations.

In practice, reconstructions obtained with the D-bar method suffer from a smoothing effect due to truncation in the scattering transform, which is necessary for finite and noisy data but leaves out all high-frequency information in the data. Thus, we cannot reconstruct sharp edges, and subsequent processing is beneficial. An early approach to overcome the smoothing is to use a nonlinear diffusion process to sharpen edges [46]. In recent years, deep learning has been highly successful for post-processing insufficient noise or artefact-corrupted reconstruction [63].

In the context of the deep D-bar method, we are given an initial analytic reconstruction operator $\mathcal{R}_{\text{d-bar}}$ that maps the measurements (i.e., the DtN map $\Lambda_{\sigma,D}$ for EIT) to an initial image, which suffers from various artefacts, primarily over-smoothing. Then a U-Net $\mathcal{G}_{\theta}$ [93] is trained to improve the reconstruction quality of the initial reconstructions, and we refer to the original publication [45] for details on the architecture. Thus, we could write this process as $\sigma\approx\mathcal{G}_{\theta}(\mathcal{R}_{\text{d-bar}}(\Lambda_{\sigma,D}))$, where the network $\mathcal{G}_{\theta}$ is trained by minimising the $\ell^{2}$-loss of D-bar reconstructions to ground-truth images. Specifically, given a collection of $N$ paired training data $\{(\sigma_{i}^{\dagger},\Lambda_{\sigma_{i}^{\dagger},D}^{\delta})\}_{i=1}^{N}$ (i.e., ground-truth conductivity $\sigma_{i}^{\dagger}$ and the corresponding noisy measurement data $\Lambda_{\sigma_{i}^{\dagger},D}^{\delta}$), we train a DNN $\mathcal{G}_{\theta}$ by minimising the following empirical loss

This can be viewed as a specialised denoising scheme to remove the artefacts in the initial reconstruction $\mathcal{R}_{\text{d-bar}}(\Lambda_{\sigma_{i}^{\dagger},D}^{\delta})$ by the D-bar reconstructor $\mathcal{R}_{\text{d-bar}}$. The loss $\mathcal{L}(\theta)$ is then minimised with respect to the DNN parameters $\theta$, typically by the Adam algorithm [66], a very popular variant of stochastic gradient descent. Once a minimiser $\theta^{*}$ of the loss $\mathcal{L}(\theta)$ is found, given a new test measurement $\Lambda_{\sigma,D}^{\delta}$, we can obtain the reconstruction $\mathcal{G}_{\theta^{*}}(\mathcal{R}_{\text{d-bar}}(\Lambda_{\sigma,D}^{\delta}))$. Thus at the testing stage, the method requires only additional feeding of the initial reconstruction $\mathcal{R}_{\text{d-bar}}(\Lambda_{\sigma,D}^{\delta})$ through the network $\mathcal{G}_{\theta^{*}}$, which is computationally very efficient. This presents one distinct advantage of a supervisedly learned map.

Several extensions have been proposed. Firstly, the need to model boundary shapes in the training data can be eliminated by using the Beltrami approach [7] instead of the classic D-bar method. This allows for domain-independent training [44]. A similar motivation is given by replacing the classic U-net that operates on rectangular pixel domains with a graph convolutional version; this way learned filters are domain and shape-independent [49, 48]. Similarly, the reconstruction from Calderón’s method [11, 96] can be post-processed using U-net, leading to the deep Calderón’s method [21]. Distinctly, the deep Calderón’s method is capable of directly recovering complex valued conductivity distributions. Finally, even the enclosure method can be improved by predicting the convex hull from values of the involved indicator function [97].

The deep sampling method (DDSM) [42] is based on the direct sampling method (DSM) due to Chow, Ito and Zou [29]. Using only one single Cauchy data pair on the boundary $\partial\Omega$, The DSM constructs a family of probing functions $\{\eta_{x,d_{x}}\}_{x\in\Omega,d_{x}\in\mathbb{R}^{n}}\subset H^{2\gamma}(\partial\Omega)$ such that the index function defined by

takes large values for points near the inclusions and relatively small values for points far away from the inclusions, where $|\cdot|_{Y}$ denotes the $H^{2\gamma}(\partial\Omega)$ seminorm in and the duality product $\langle f,g\rangle_{\gamma,\partial\Omega}$ is defined by

where $-\Delta_{\partial\Omega}$ denotes the Laplace-Beltrami operator, and $(-\Delta_{\partial\Omega})^{\gamma}$ its fractional power via spectral calculus. Let the Cauchy difference function $\varphi$ be defined by

Then the index function $\mathcal{I}(x,d_{x})$ can be equivalently rewritten as

Motivated by the relation between the index function $\mathcal{I}(x,d_{x})$ and the Cauchy difference function $\varphi$ and to fully make use of multiple pairs of measured Cauchy data, Guo and Jiang [42] proposed the DDSM, employing DNNs to learn the relationship between the Cauchy difference functions $\varphi$ and the true inclusion distribution. That is, DSSM construct and train a DNN $\mathcal{G}_{\theta}$ such that

where $\{\varphi_{i}\}_{i=1}^{N}$ correspond to $N$ pairs of Cauchy data $\{g_{\ell},\Lambda_{\sigma,N}g_{\ell}\}_{\ell=1}^{N}$. Guo and Jiang [42] employed a CNN-based U-Net network for DDSM, and later [41] designed a U-integral transformer architecture (including comparison with state-of-the-art DNN architectures, e.g., Fourier neural operator, and U-Net). In our numerical experiments, we choose the U-Net as the network architecture for DDSM as we observe that U-Net can achieve better results than the U-integral transformer for resolution $64\times 64$. For higher resolution cases, the U-integral transformer seems to be a better choice due to its more robust ability to capture long-distance information. The following result [42, Theorem 4.1] provides some mathematical foundation of DDSM.

The idea of DDSM was extended to diffusive optical tomography in [43]. Ning et al. [86] employ the index functions obtained from the DSM as the input of the DNN for solving inverse obstacle scattering problems.

Li et al. [74] proposed using CNN to directly learn the map from the measured data and the conductivity distribution. The employed network architecture is based on LeNet and refined by applying dropout layer and moving average. The CNN architecture used in the numerical experiments below is shown in Fig. 1. Since the number of injected currents and the discretisation size differ from that in [74], we modify the input size, network depth, kernel size, etc. The input size is 32 × 64. The kernel size is 5 × 5 with zero-padding max pooling rather than average pooling is adopted to gain better performance. The sigmoid activation function used in LeNet causes a serious saturation phenomenon, which can lead to vanishing gradients. So, ReLU is chosen as the activation function below. A dropout layer is added to improve the generalisation ability of this model. One-half of the neurons before the first fully connected layers are randomly discarded from the network during the training process. It can reduce the complex co-adaptation among neurons so that the network can learn more robust features. In addition, a dropout layer has been proven to be very effective in training large datasets.

Chen et al. [26] proposed a novel deep learning architecture by adding a fully connected layer before the U-Net structure. The input of the network is given by the difference voltage $u_{\sigma}^{\delta}-u_{\sigma_{0}}$. Inspired by a linearized approximation of the EIT problem for a small perturbation of conductivity distribution $\sigma-\sigma_{0}$:

where $\mathbf{J}$ donates the sensitivity matrix, the method first generates an initial guess of the conductivity distribution $\sigma$ from the linear fully connected(FC) layer followed by a ReLU layer and then feeds it to a denoising U-Net model to learn the nonlinear relationship further. Thus we could write this process as $\sigma\approx\mathcal{G}_{\theta}(u_{\sigma}^{\delta}-u_{\sigma_{0}})=\mathcal{G}_{\theta_{2}}(\mathcal{G}_{\theta_{1}}(u_{\sigma}^{\delta}-u_{\sigma_{0}}))$ with $\mathcal{G}_{\theta_{1}}=\text{FC+ReLU}$ and $\mathcal{G}_{\theta_{2}}=\text{U-Net}$. The authors also proposed an initialisation strategy to further help obtain the initial guess, i.e., the weights $\theta_{1}$ of the fully connected layer are initialised with the least-squares solution using training data. The weights $\theta_{2}$ for the U-Net are initialised randomly as usual. Then, all weights $\theta=\theta_{1}\cup\theta_{2}$ are updated during the training process. According to the numerical results shown in [26], this special weight initialization strategy can reduce the training time and improve the reconstruction quality. With a trained network, different from the deep D-bar and DDSM methods, the methods FC-UNet and CNN based on LeNet only involve a forward pass of the trained network for each testing example.

Based on our numerical experience, dropping the ReLU layer following the fully connected layer can provide better reconstruction results, at least for the examples in section 4. Thus, for the numerical experiments, we employ the FC-UNet network as shown in Fig. 2, in which only a linear fully connected layer is employed before the U-Net.

In addition, by employing the FC-UNet to extract structure distribution and a standard CNN to extract conductivity values, a structure-aware dual-branch network was designed in [25] to solve EIT problems.

Generating simulated data consists of three main parts, which we describe below. The codes for data generation are available at https://github.com/dericknganyu/EIT_dataset_generation.

In the 2D setting, we generate $N$ circular phantoms $\{P_{i}\_{i=1}^{N}$, all restricted to the unit circle centred at the origin, i.e., $\Omega=\{(x,y):x^{2}+y^{2}\leq 1\}$ in the Cartesian coordinates or $\{(r,\theta):r\leq 1,\theta\in[0,2\pi]\}$ in polar coordinates. The phantoms are generated randomly. Firstly, we decide on the maximum number $M\in\mathbb{N}$ of inclusions. Each phantom would then contain $n$ inclusions, where $n\in\mathcal{U}\{1,\ldots,M\}$, the uniform distribution over the set $[1,\ldots,M]$. To mimic realistic scenarios in medical imaging, the inclusions are elliptical and are sampled such that when $n>1$, the inclusions do not overlap. Since the inclusions are elliptical, each inclusion, $(E_{j})_{j=1}^{n}$ is characterised by a centre $C_{j}=(h_{j},k_{j})$, an angle of rotation $\alpha_{j}$, a major and minor axis $a_{j}$ and $b_{j}$ respectively. The parametric equation of an ellipse $E_{j}$ is thus given by

To mimic realistic scenarios in medical imaging, the inclusions are sampled to avoid contact with the boundary $\partial\Omega$ of the domain $\Omega$. For an inclusion $E_{j}$, we have $x^{2}+y^{2}<0.9$ for any $(x,y)\in E_{j}$. In this way, all phantoms have inclusions contained within $\Omega$. We illustrate this in Algorithm 3.

Each phantom $P_{i},i\in\{1,2,\ldots,N\},$ has $(E_{j})_{j=1}^{n}$ inclusions, with $n\in\mathcal{U}\{1,\ldots,M\}$. For each $E_{j}$, we assign a conductivity $\sigma^{i}_{j}\in\Sigma_{j}:=\mathcal{U}(0.2,0.8)\cup\mathcal{U}(1.2,2.0)$. The background conductivity is set to $1$. In this way, given a point $(x,y)\in P_{i}$ in the domain/phantom, the conductivity $\sigma_{i}(x,y)$ at that point is therefore given by

Fig. 3(b) shows an example of a phantom generated in this way.

Next, for any simulated $\sigma$, we solve the forward problem (1) using the Galerkin finite element method (FEM) [71, 39], for the injected currents $g_{1}$ and $g_{2}$ in (15) around the boundary $\partial\Omega$. The points $(x,y)\in P_{i}$ are thus nodes in the finite element mesh shown in Fig. 3(a)

We use the MATLAB PDE toolbox in the numerical experiment to solve the forward problem.

In real-life situations, the conductivities of the inclusions are rarely constant. Indeed, usually, there are textures on internal organs in medical applications. Motivated by this, we take a step further in generating phantoms, with inclusions having variable conductivities. This introduces a novel challenge to the EIT problem, and we seek to study its impact on different reconstruction algorithms. The procedure to generate simulated data remains unchanged. However, $\sigma^{i}_{j}$ in equation (14) becomes

where $f:\mathbb{R}^{2}\ni(x,y)\mapsto\tfrac{1}{2}\left(\sin{k_{x}x}+\sin{k_{y}y}\right)\in[-1,1]$, $R_{\alpha_{j},C_{j}}$ is the rotation of centre $C_{j}=(h_{j},k_{j})$ and angle $\alpha_{j}$, with respect to the centre and angle of the ellipse $E_{j}$ respectively; and $s$ applies a scaling so that the resulting $\sigma^{i}_{j}$ is either within the range $[0.2,0.8]$ or $[1.2,2.0]$. Fig. 3(c) shows an example phantom.