Sparse reconstruction of log-conductivity in current density impedance tomography

TL;DR

This paper introduces a novel non-linear optimization method for sparse reconstruction of log-conductivities in current density impedance tomography, combining $L^1$ regularization and anisotropic diffusion to improve image resolution.

Contribution

It presents a new framework that enhances image quality in impedance tomography by integrating sparsity promotion and edge enhancement techniques.

Findings

Demonstrates superior reconstruction quality over existing methods.

Effectively captures high-contrast conductivity patterns.

Numerical experiments validate the approach's robustness.

Abstract

A new non-linear optimization approach is proposed for the sparse reconstruction of log-conductivities in current density impedance imaging. This framework comprises of minimizing an objective functional involving a least squares fit of the interior electric field data corresponding to two boundary voltage measurements, where the conductivity and the electric potential are related through an elliptic PDE arising in electrical impedance tomography. Further, the objective functional consists of a $L^1$ regularization term that promotes sparsity patterns in the conductivity and a Perona-Malik anisotropic diffusion term that enhances the edges to facilitate high contrast and resolution. This framework is motivated by a similar recent approach to solve an inverse problem in acousto-electric tomography. Several numerical experiments and comparison with an existing method demonstrate the effectiveness of the proposed method for superior image reconstructions of a wide-variety of log-conductivity patterns.

Tables

Table 1. Table 1: Convergence properties of the VIP algorithm for various values v 𝑣 v of θ , c 1 , c 2 𝜃 subscript 𝑐 1 subscript 𝑐 2 \theta,c_{1},c_{2}

\hlineB3 Values	$\theta$	$c_1$	$c_2$
\hlineB3 $v\ge 2$	No convergence	No convergence	Slow Convergence
	No convergence	Good convergence	Slow Convergence
	No convergence	Convergence	Slow Convergence
	Convergence	Convergence	Slow Convergence
$v\ll 1$	Convergence	Slow Convergence	Good Convergence
\hlineB3

Equations

$$-\nabla \cdot (\sigma (x)\nabla u(x))$$

$$\sigma (x)\frac{\partial u}{\partial \nu }(x)$$

$$H(\sigma (x))=\sigma (x)|\nabla u(x)|,x\in \Omega .$$

$$\displaystyle\nabla\cdot\Bigg{(}\dfrac{H}{|\nabla u|}\nabla{u}\Bigg{)}=0$$

$$\frac{H}{|\nabla u|}\frac{\partial u}{\partial \nu }=f$$

$$-$$

$$u(x,y){|}_{\Gamma }=f_D(x,y),$$

$$\mathcal{L}(u,\sigma ,f_D)=0,$$

$$H(\sigma )=e^{\sigma }|\nabla u|$$

$$J(\sigma ,u_1,u_2)=$$

$$+\gamma \Vert \sigma {\Vert }_{L^1(\Omega )}+\frac{\delta }{2}{\int }_{\Omega }\log (1+|\nabla \sigma (x,y){|}^2)dxdy$$

$$\underset{\sigma }{\min }$$

$$\mathcal{L}(u_1,\sigma ,f_D^1)=0,$$

$$\mathcal{L}(u_2,\sigma ,f_D^2)=0.$$

$$J_{PM}(\sigma )={\int }_{\Omega }\log (1+|\nabla \sigma (x,y){|}^2)dxdy.$$

$$\hat{J}(\sigma )=J(\sigma ,u_1(\sigma ),u_2(\sigma )),$$

$$\underset{\sigma \in L_{ad}}{\min }\hat{J}(\sigma ).$$

$$\hat{J}={\hat{J}}_1+{\hat{J}}_2,J_i:L_{ad}\to {\mathbb{R}}^{+},i=1,2,$$

$${\hat{J}}_1(\sigma )=\frac{{\alpha }_1}{2}\Vert e^{\sigma }|\nabla u_1|-g_1^{\delta }{\Vert }_{L^2(\Omega )}^2+\frac{{\alpha }_2}{2}\Vert e^{\sigma }|\nabla u_2|-g_2^{\delta }{\Vert }_{L^2(\Omega )}^2+\frac{\beta }{2}\Vert \sigma {\Vert }_{L^2(\Omega )}^2$$

$$+\frac{\delta }{2}{\int }_{\Omega }\ln (1+|\nabla \sigma (x,y){|}^2)dxdy,$$

$${\hat{J}}_2(\sigma )=\gamma \Vert \sigma {\Vert }_{L^1(\Omega )}.$$

$$\partial\hat{J}(\bar{\sigma}):=\Bigg{\{}{\phi\in L^{*}_{ad}:\liminf_{\sigma\rightarrow\bar{\sigma}}\dfrac{\hat{J}(\sigma)-\hat{J}(\bar{\sigma})-\langle\phi,\sigma-\bar{\sigma}\rangle}{\|\bar{\sigma}-\sigma\|_{2}}}\geq 0\Bigg{\}},$$

$$\partial \hat{J}(\sigma )=\nabla {\hat{J}}_1(\sigma )+\partial {\hat{J}}_2(\sigma ),$$

$$\partial (\alpha \hat{J})=\alpha \partial \hat{J}.$$

$$0\in \partial \hat{J}({\sigma }^{\ast }),$$

$$-\nabla {\hat{J}}_1({\sigma }^{\ast })\in \partial {\hat{J}}_2({\sigma }^{\ast }).$$

$$⟨\nabla {\hat{J}}_1({\sigma }^{\ast })+\lambda ,\sigma -{\sigma }^{\ast }⟩\ge 0,\forall \sigma \in L_{ad}.$$

$$\lambda \in {\Lambda }_{ad}:=\{\lambda \in L^2(\Omega ):0\le \lambda \le \gamma ,\text{mbox}a.e.in\Omega \}.$$

$${\nabla }_{\sigma }{\hat{J}}_1({\sigma }^{\ast })+{\lambda }^{\ast }+{\lambda }_{{\sigma }_u}^{\ast }-{\lambda }_{{\sigma }_l}^{\ast }=0,$$

$${\lambda }_{{\sigma }_u}^{\ast }\ge 0,{\sigma }_u-{\sigma }^{\ast }\ge 0,⟨{\lambda }_{{\sigma }_u}^{\ast },{\sigma }_u-{\sigma }^{\ast }⟩=0,$$

$${\lambda }_{{\sigma }_l}^{\ast }\ge 0,{\sigma }^{\ast }-{\sigma }_l\ge 0,⟨{\lambda }_{{\sigma }_l^{\ast }},{\sigma }^{\ast }-{\sigma }_l⟩=0,$$

$${\lambda }^{\ast }=\gamma \text{mbox}a.e.on\{x\in \Omega :{\sigma }^{\ast }(x)>0\},$$

$$0\le {\lambda }^{\ast }\le \gamma \text{mbox}a.e.on\{x\in \Omega :{\sigma }^{\ast }(x)=0\}.$$

$${\nabla }_{\sigma }{\hat{J}}_1({\sigma }^{\ast })=$$

$$\displaystyle-\delta\nabla\cdot\Bigg{(}\dfrac{\nabla\sigma^{*}}{1+|\nabla\sigma^{*}|^{2}}\Bigg{)}$$

$$-\nabla \cdot (e^{{\sigma }^{\ast }}\nabla v_1)$$

$$v_1{|}_{\Gamma }$$

Full text

Keywords: Inverse problems, PDE-constrained optimization, proximal methods, edge-enhancement, sparsity patterns, current density impedance imaging.

MSC: 35R30, 49J20, 49K20, 65M08, 82C31

1 Introduction

Electrical impedance tomography (EIT) is an imaging modality, where one attempts to recover the conductivity of a body from the boundary measurement of current and voltage [14]. The underlying inverse problem is highly ill-posed and non-linear yet very important due to its wide range applications in the fields such as medical imaging [54] and engineering [27, 53]. The mathematical formulation of the EIT inverse problem is given by the following conductivity equation

[TABLE]

where $\Omega\subset\mathbb{R}^{n}$ is a convex and bounded domain with Lipschitz boundary and $\Gamma$ is the boundary of $\Omega$. In this model, $\sigma$ is the electrical conductivity, $u$ represents the electric potential and $f$ is the current applied to the boundary.

The reconstructions obtained through the EIT setup usually have high contrast but limited spatial resolution [48]. On the other hand, reconstructions obtained through ultrasound imaging has very high resolution but limited contrast [3, 42]. In recent years, attempts have been made to combine multiple imaging modalities to obtain image reconstructions with both high contrast and high resolution. This led to the emergence of hybrid imaging methods that belong to class of coupled-physics imaging modalities to generate images of superior quality. One of such imaging methods, known as current density impedance imaging (CDII) combines the classical EIT setup with magnetic resonance (MR) scanning [19, 36]. It is alternatively known as magnetic resonance EIT (MREIT). Current or voltage is applied through the electrodes, which give rise to an interior electric field and the corresponding generated magnetic field, represented as $B=(B_{x},B_{y},B_{z})$, is measured by the MR scanner. The corresponding inverse problem is to solve for the conductivity $\sigma$ from $B_{z}$ using the well-known iterative Harmonic $B_{z}$-algorithm [17, 30]. Convergence of the harmonic $B_{z}$ algorithm has been well-studied [17, 30, 31]. In particular, it has been shown that for small contrast values of the target conductivity, the harmonic $B_{z}$-algorithm is stable and convergent, provided we have a good initial guess [30]. Thus, it is not clear that one can recover good quality images for high contrast objects through Harmonic $B_{z}$-algorithm.

An alternate approach to solve the CDII inverse problem is to use the knowledge of interior electric field, which is obtained from the magnetic field. Correspondingly, the magnitude of the interior electric field is also determined [47, 49], which is given by

[TABLE]

The formulation of reconstruction problem is as follows: Given the boundary data $f$ for, possibly, several choices of boundary patterns and the corresponding interior measurement data $H$, find the conductivity distribution $\sigma$. In this framework, we use the internal function $H(\sigma)$ to replace $\sigma$ in the EIT equation (1) to get the following nonlinear equation

[TABLE]

For the CDII inverse problem, the solution to the boundary value problem (3) is crucial but it is difficult to use it in practice because of its highly nonlinear behaviour and also because the data represented by the measured values of $H$ enter as a coefficient of the differential model [48]. Even with the additional measurements, analysis and application of the 1-Laplacian relies on an iterative localized algorithm, wherein one considers an approximation of the CDII problem. This subsequently led to several computational approaches in solving the CDII inverse problem. In [23], it was proved that the linearized problem is elliptic and hence solvable, if there are at least $n$ set of measurements $\{H_{i}(\sigma)\}_{i=1}^{n}$ and corresponding to $n$ boundary data $\{f_{i}\}_{i=1}^{n}$ such that $\nabla u_{i}$ and $\nabla u_{j}$ are nowhere collinear for $i\neq j$. It has been shown in [26] that the solution of the above 1-Laplacian equation with the Neumann boundary condition is non-existent unless additional measurements with different boundary current patterns are used. Recovery of isotropic conductivity in regions where the magnetic field is transversal using two internal current distributions was done using an explicit local formula [28]. Moreover, using the information of two internal current distributions, the authors in [26] uniquely determine the singular support of the conductivity function. In [34], the authors showed that the conductivity in the planer domain can be recovered from a single voltage-current on a part of boundary and the magnitude of one interior current density. In the same article, they also provide sufficient conditions on Dirichlet boundary data to guarantee unique recovery of conductivity. In [35], the recovery of Hölder continuous conductivities have been establised for domains with connected boundary from the interior measurement of the magnitude of one current density. Determination of isotropic conductivity variations from measurements of two current density vector fields was studied in [19]. In [51], authors showed the recovery of planar conductivities by solving the 1-Laplace equation with partial boundary data. The authors in [32, 8, 9] present explicit reconstruction formulae for recovering the conductivity distribution from multiple interior measurements in two and higher dimensions.

The well-known numerical reconstruction algorithm using the internal current distribution is an iterative $J$-substitution algorithm which was first introduced by [24] and subsequently considered in other works, see for e.g., in [21, 25, 35, 36]. It has been shown that the $J$-substitution algorithm is able to reconstruct the conductivity with high resolution [25, 26]. Another numerical reconstruction iterative method is the regularized D-bar method [22] that provides images with high resolution. In [33], the authors use an alternating split Bregman algorithm for solving a minimization problem related to the energy functional corresponding to the 1-Laplacian equation (3). Also, in [20], Picard and Newton type algorithms are implemented to solve the 1-Laplacian problem. But there is not enough evidence to suggest that these existing algorithms (linearized or localized iterative methods) can provide high contrast images, specially for objects with holes or inclusions, which are inherent to CDII reconstructions.

The CDII inversion problem can be viewed as an inverse problem of estimating the conductivity parameter from the conductivity partial differential equation (PDE) (1). A robust way of solving this inverse problem is to formulate it as an optimal control problem, where the condiuctivity parameter is the control variable that drives the interior electric field close to the measured value, with dynamics governed by the conductivity PDE. Such optimal control methods has been used previously in the context of ultrasonically-induced Lorentz force electrical impedance tomography [5], magnetoacoustic tomography [6, 7, 38] and acousto-electric tomography (AET) [1, 41]. In [5], the authors use an optimal control framework to recover the conductivity distribution from the measurements of current induced by static magnetic field through the Lorentz force. They solve the optimal control problem using an orthogonal field method. In [6, 7, 38], the authors use an optimal control approach to reconstruct conductivity distribution of biological tissue from measurements of the Lorentz force induced tissue vibration. In [1, 41], the authors reconstruct log-conductivity in acousto-electric tomography (AET) using an optimal control formulation based on the theory of anisotropic diffusion to potentially obtain reconstructions with high resolution and contrast. The results obtained in these papers demonstrated that such optimal control frameworks are robust and accurate for imaging modalities arising through a partial differential equation (PDE).

In this paper, we consider a similar optimization framework developed in [41] for reconstructing the log-conductivity in CDII. We formulate a minimization problem, where given interior electric field intensity data, we aim at determining the variation in conductivity from a known background conductivity. We, further, assume that this variation demonstrates a sparsity pattern. Such patterns arise frequently in several tomographic imaging scenarios, for e.g. in blood vessel tomographic reconstructions [29]. This is incorporated in our model through a $L^{2}-L^{1}$ regularization term in our objective functional. To obtain sharp edges and, thus, improve spatial resolution of the reconstructed images, we use a Perona-Malik anisotropic diffusion filtering term in our functional. The resulting optimality system gives rise to an elliptic adjoint equation with a $L^{2}$ source term. Classical cell-nodal finite difference schemes are not applicable for solving such equations. We, thus, use a averaged cell-nodal scheme to solve such equations. Finally, we solve the optimization problem using a variable inertial proximal scheme that efficiently handles the non-differentiable terms in the objective functional. We demonstrate through several examples that our method can be used to obtain superior quality reconstructions for objects with holes and inclusions. In this context we would like to remark that our framework can also be used for obtaining superior reconstruction of anisotropic conductivity distributions, in ultrasonically-induced Lorentz force electrical impedance tomography and magnetoacoustic tomography. Further, our proposed framework can also be used in other hybrid imaging modalities like quantitative photoacoustic tomography [16] and quantitative thermoacoustic tomography [4] to obtain better reconstructions.

The article is organized as follows: In the Section 2, we formulate the minimization problem for the CDII. In the Section 3, we present some theoretical results about our optimization problem. We also characterize the optimality system. The variable inertial proximal scheme and the averaged cell-nodal schemes to solve the optimization problem are discussed in Section 4. In the Section 5, we present simulation results of our CDII framework and compare them with the reconstructions obtained using the Picard scheme proposed in [20], which validate our framework for CDII and demonstrate the effectiveness of our method to reconstruct wide variety of objects with corners, holes and inclusions. A section on conclusions completes our work.

2 A minimization problem

We consider the conductivity equation in ${\mathbb{R}}^{2}$ arising in EIT

[TABLE]

where $\Omega\subset\mathbb{R}^{2}$ is bounded, $\Gamma$ is the boundary of $\Omega$, $e^{\sigma}$ is the conductivity coefficient and $u\in H^{1}_{f_{D}}(\Omega)=\{u\in H^{1}(\Omega):u=f_{D}\mbox{ on }\Gamma\}$ is the electric potential.

We assume that $\sigma$ is a sparse conductivity coefficient which we want to recover, given the fact that the conductivity of the background is 1. The conductivity equation (4) can also be written as

[TABLE]

where $\sigma(x,y)\in L_{ad}=\{\sigma\in H_{0}^{1}(\Omega):\sigma_{l}\leq\sigma(x,y)\leq\sigma_{u},~{}\forall(x,y)\in\Omega\}$, $\sigma_{u}>0$ and $\sigma_{l}<0$ .

We consider an optimization-based approach for reconstructing $\sigma$ given $H_{1}(\sigma),H_{2}(\sigma)$, where

[TABLE]

is the interior electric field corresponding to the voltage potential $u$. We consider the following cost functional

[TABLE]

where $u_{1},u_{2}$ satisfy (4) with boundary data $f_{D}^{1},f_{D}^{2}$. We now consider the following minimization problem

[TABLE]

The term $\gamma\,\|\sigma\|_{L^{1}(\Omega)}$, $\gamma>0$ in the functional, defined in (5), implements a $L^{1}$ regularization of the minimization problem that promotes sparsity patterns in the reconstruction of conductivity. Such a regularization method mirrors the well known compressed-sensing technique; see [11]. In recent past, optimal control with $L^{1}$ cost functionals has become a topic of major interest [50], because one obtains sparse controls through this procedure, which finds numerous applications. The motivation for sparse log-conductivity patterns is based on the assumption that the background conductivity is known to be 1 in a substantial part of the domain $\Omega$ after normalization and varies considerably from this value in correspondence to different kind of objects present within the domain.

The combined $L^{2}$-$L^{1}$ regularization allows for the reconstruction of conductivity, and thus the imaging of, possibly, irregular objects inside $\Omega$. This does not serve the ultimate goal of reconstructing objects like tissues in medical imaging, which are more regular, save for the edges that eventually define them. We infuse this additional aprior knowledge into our model through the last term in our functional (5) that, commonly, appears in the field of anisotropic diffusion. Such a term plays an important role in dampening image noise while keeping significant parts of the image content such as edges and other anatomical details that are of utmost importance in the interpretation of the image. Anisotropic diffusion means non-uniform diffusion in different directions. The regions where $|\nabla\sigma|$ is very small corresponds to noise and thus, the process of smoothening occurs. At the edges or singularities of an object, where the value of $|\nabla\sigma|$ is large, there is a small amount of smoothening and this preserves the edges. A standard technique to implement anisotropic diffusion, in order to obtain a good contrast, is to use a total variation (TV) regularization [13, 43]. But, this regularization method gives rise to a non-differentiable term in the functional (5), thus requiring more sophisticated optimization algorithms. On the other hand, anisotropic diffusion is inherent to Perona-Malik (PM) filtering [37]. It is well-known that the diffusion process governed by the PM equation leads to a decrease in the total variation during its evolution [44]. We, thus, choose the energy functional of the Perona-Malik equation for anisotropic diffusion [37]. One can note that the PM regularization term is differentiable and, thus, easier to handle than the TV regularization term.

Mathematically, one can consider the PM filtering as the gradient flow generated by the non-convex and lower semi-continuous functional given by

[TABLE]

We refer to [12, 44] for a general introduction to anisotropic diffusion and a detailed discussion on the PM functional. Further, in [41], the PM model was used in the reconstruction of log-conductivities in AET and it was observed that the reconstructions obtained demonstrated superior contrast and resolution. Thus, for the current setup in CDII, we use a similar PM anisotropic diffusion filter to facilitate high contrast and high resolution images.

3 Theory of the minimization problem

In this section, we discuss the existence of solutions of the minimization problem (P) and its characterization through a first-order optimality system. We refer to this minimization problem as the CDII sparse reconstruction problem (CDII-SR). Our analysis of this problem begins with the discussion concerning the existence and uniqueness of weak solutions of $\mathcal{L}(u,\sigma,f_{D})=0$, which can be proved by standard arguments of Riesz representation theorem [18, Chapter 8].

Proposition 1**.**

Let $\sigma\in L_{ad}$ and $f_{D}\in H^{1/2}(\partial\Omega)$. Then the problem (4) has a unique solution in $H^{1}_{f_{D}}(\Omega)$.

The solvability of the CDII inversion problem depends on the type of Dirichlet boundary data $f_{D}^{j},~{}j=1,2$. In this, context, we have the following lemma from [2].

Lemma 3.1** (Boundary data).**

Let $\Omega\subset{\mathbb{R}}^{2}$ be a bounded simply connected open set, whose boundary $\Gamma$ is a simple closed curve. Let $f=(f^{1},f^{2})$ be a mapping $\Gamma\to{\mathbb{R}}^{2}$ which is a homeomorphism of $\Gamma$ onto a convex closed curve $C$, and let $D$ denote the bounded convex domain bounded by $C$. Let $\sigma\in L^{\infty}(\Omega)$, and let $U=(u_{1},u_{2})$ be the $e^{\sigma}$-harmonic mapping whose components $u_{1}$ and $u_{2}$ are solutions to the Dirichlet problem (4) with $f_{D}=f_{D}^{1}$ and $f_{D}=f_{D}^{2}$, respectively, and $f_{D}^{J}\in H^{1}(\Omega)\cap C(\bar{\Omega})$, $J=1,2$. Then $U$ is a homeomorphism of $\Omega$ onto $D$. In particular, for all $\omega\subset\subset\Omega$ we have either $\det(\nabla u_{1},\nabla u_{1})>0$ or $\det(\nabla u_{1},\nabla u_{1})<0$ almost everywhere in $\omega$.

In [23], the authors have shown that in 2D, the boundary condition pair $f_{D}^{1}=x$ and $f_{D}^{2}=y$ satisfies the conditions of Lemma 3.1 and, thus, the corresponding solutions to (4) $u_{1}$ and $u_{2}$ have no critical points and $\nabla u_{1},\nabla u_{2}$ are not collinear in $\bar{\Omega}$. We will use these boundary conditions for our numerical experiments in Section 5.

Next, we consider the Fréchet differentiability of the mapping $u(\sigma)$.

Lemma 3.2**.**

The map $u(\sigma)$ defined by (4) is Fréchet differentiable as a mapping from $L_{ad}$ to $H^{1}_{f_{D}}(\Omega)$.

For the proof of this Lemma, we refer to [23]. Using Lemma 3.2, we introduce the reduced cost functional

[TABLE]

where $u_{i}(\sigma)$, $i=1,2$ denotes the unique solution of (4) given $\sigma$ and $f_{D}^{i},i=1,2$. The constrained optimization problem (P) can be formulated as an unconstrained one as follows

[TABLE]

We next investigate the existence of a minimizer to the CDII-SR problem (P). We first consider the case when $\delta=0$, i.e., the Perona-Malik term in the functional $J$ is absent.

Proposition 2**.**

Let $f_{D}^{1},f_{D}^{2}\in H^{1/2}(\Omega)$ such that $|\nabla u_{1}|>0,|\nabla u_{2}|>0$ and let $\delta=0$. Then there exists a triplet $(\sigma^{*},u_{1}^{*},u_{2}^{*})\in L_{ad}\times H^{1}_{f_{D}^{1}}(\Omega)\times H^{1}_{f_{D}^{2}}(\Omega)$ such that $u_{i}^{*},i=1,2$ are solutions to $\mathcal{L}(\sigma,u_{i},f^{i}_{D})=0,i=1,2$ and $\sigma^{*}$ minimizes $\hat{J}$ in $L_{ad}$.

Proof.

Boundedness from below of $\hat{J}$ guarantees the existence of a minimizing sequence $(\sigma^{m})$. Since $L_{ad}$ is reflexive and $\hat{J}$ is sequentially weakly lower semi-continuous, this sequence is bounded. Therefore it contains a weakly convergent subsequence $(\sigma^{m_{l}})$ in $L_{ad}$, $\sigma^{m_{l}}\rightharpoonup\sigma^{*}$. Correspondingly, the sequence $(u_{1}^{m_{l}},u_{2}^{m_{l}})$, where $u_{i}^{m_{l}}=u_{i}(\sigma^{m_{l}})$, is bounded in $H^{1}_{f_{D}^{1}}(\Omega)\times H^{1}_{f_{D}^{2}}(\Omega)$. Therefore the sequence converges weakly to $(u_{1}^{*},u_{2}^{*})$. Now, using the Rellich Kondrachev compactness theorem in $\mathbb{R}^{2}$, we have that $L_{ad}$ is compactly embedded in $L^{2}(\Omega)$. This results in a strong convergence of the subsequence $\sigma^{m_{l}}$ in $L^{2}(\Omega)$ to $\sigma^{*}$. We, now, consider the weak formulation of the solutions of the elliptic problem (4) and, thus, focus on $\langle\nabla\cdot(\sigma^{m_{l}}\nabla u_{i}^{m_{l}}),\psi\rangle_{L^{2}(\Omega)}$ for any $\psi\in H^{1}_{0}(\Omega)$. Using integration by parts, we have $\langle\nabla\cdot(\sigma^{m_{l}}\nabla u_{i}^{m_{l}}),\psi\rangle_{L^{2}(\Omega)}=-\langle\sigma^{m_{l}}\nabla u_{i}^{m_{l}},\nabla\psi\rangle_{L^{2}(\Omega)}$. From the above discussion, the sequence of products $\sigma^{m_{l}}\nabla u_{i}^{m_{l}}$ is weakly convergent in $L^{2}(\Omega)$, that is, $\langle\sigma^{m_{l}}\nabla u_{i}^{m_{l}},\nabla\psi\rangle_{L^{2}(\Omega)}\to\langle\sigma^{*}\nabla u_{i}^{*},\nabla\psi\rangle_{L^{2}(\Omega)}$. With this preparation and using the continuity of the maps $u_{i}(\sigma)$, it follows that $(u_{1}^{*},u_{2}^{*})=(u_{1}(\sigma^{*}),u_{2}(\sigma^{*}))$, and the triplet $(\sigma^{*},u_{1}^{*},u_{2}^{*})$ minimizes the objective $\hat{J}$. ∎

In the case $\delta\neq 0$, we first note that the function $\log(1+z^{2})$ is not convex. Therefore the PM functional, and, hence, the functional $\hat{J}$ in (5) is not weakly lower semi-continuous on $W^{1,p}(\Omega)$ for any $1<p<\infty$. Nevertheless, $\hat{J}$ is a bounded below, lower semi-continuous Lipschitz functional, for which a minimizer exists, provided that $L_{ad}$ is compact.

3.1 Characterization of local minima

To characterize the solution of our optimization problem through first-order optimality conditions, we write the reduced functional $\hat{J}$ as

[TABLE]

where

[TABLE]

Remark 3.1**.**

The functional $\hat{J}_{1}$ is smooth and possibly non-convex, while $\hat{J}_{2}$ is non-smooth and convex.

We next state some properties of the reduced functional $\hat{J}_{1}(\sigma)$ which can be proved using the arguments in [23, Lemma 3.1].

Proposition 3**.**

The reduced functional $\hat{J}_{1}(\sigma)$ is weakly lower semi-continuous, bounded below and Fréchet differentiable.

We now define the subdifferential of a non-smooth functional.

Definition 3.1** (Subdifferential).**

If $\hat{J}$ is finite at a point $\sigma$, the Fréchet subdifferential of $\hat{J}$ at $\sigma$ is defined as follows [15]

[TABLE]

where $L^{*}_{ad}$ is the dual space of $L_{ad}$. An element $\phi\in\partial\hat{J}(\sigma)$ is called a subdifferential of $\hat{J}$ at $\sigma$.

In our setting, we have the following

[TABLE]

since $\hat{J}_{1}$ is Fréchet differentiable by Prop. 3. Moreover, for each $\alpha>0$, it holds that

[TABLE]

The following proposition gives a necessary condition for a local minimum of $\hat{J}$ (see [41]).

Proposition 4** (Necessary condition).**

If $\hat{J}=\hat{J}_{1}+\hat{J}_{2}$, with $\hat{J}_{1},\hat{J}_{2}$ given by (8), attains a local minimum at $\sigma^{*}\in L_{ad}$, then

[TABLE]

or equivalently

[TABLE]

The following variational inequality holds for each $\lambda\in\partial\hat{J}_{2}(\sigma^{*})$ (see [50]).

[TABLE]

Using the definition of $\hat{J}_{2}$ in (8) and the fact that $L_{ad}$ is reflexive, the inclusion $\lambda\in\partial\hat{J}_{2}(\sigma^{*})$ gives the following characterization of space of $\lambda$

[TABLE]

A pointwise analysis of the variational inequality (10) leads to the existence of a non-negative functions $\lambda_{\sigma_{l}}^{*},\lambda_{\sigma_{u}}^{*}\in L^{2}(\Omega)$ that correspond to Lagrange multipliers for the inequality constraints in $L_{ad}$. We, thus, have the following first-order optimality system.

Proposition 5** (First-order necessary conditions).**

The optimal solution of the minimization problem (7) can be characterized by the existence of $(\lambda^{*},\lambda_{\sigma_{l}}^{*},\lambda_{\sigma_{u}}^{*})\in\Lambda_{ad}\times L^{2}(\Omega)\times L^{2}(\Omega)$ such that

[TABLE]

The conditions (12)-(15) are known as the complementarity conditions for $(\sigma^{*},\lambda^{*})$.

To determine the gradient $\nabla_{\sigma}\hat{J}_{1}$, we use the adjoint approach (see for e.g., [39, 40]). This gives the following reduced gradient of $\hat{J}_{1}$

[TABLE]

where $u_{1},u_{2}$ satisfy the forward equations $\mathcal{L}(u_{1},\sigma^{*},f_{D}^{1})=0,~{}\mathcal{L}(u_{2},\sigma^{*},f_{D}^{2})=0$, respectively, and $v_{1},v_{2}$ satisfy the adjoint equations

[TABLE]

[TABLE]

The complementarity conditions (12)-(15) can be rewritten in a compact form as follows. Define

[TABLE]

Then the triplet $(\lambda^{*},\lambda_{\sigma_{l}}^{*},\lambda_{\sigma_{u}}^{*})$ is obtained by solving the following equations

[TABLE]

(see [50]). For each $k\in{\mathbb{R}}^{+}$, define the following quantity

[TABLE]

The following lemma determines the complementarity conditions (12)-(15) in terms of $E$ (see [50, Lemma 2.2]).

Lemma 3.3**.**

The complementarity conditions (12)-(15) are equivalent to the following

[TABLE]

where $\mu$ is defined in (19).

Using the gradients in (16) and Lemma 3.3, the optimality conditions (11)-(15) for the CDII-SR problem can be rewritten as follows

Proposition 6**.**

A local minimizer $(u_{1},u_{2},\sigma^{*})$ of the problem (P) can be characterized by the existence of $(v_{1},v_{2},\mu^{*})\in H^{1}_{0}(\Omega)\times H^{1}_{0}(\Omega)\times L_{ad}$, such that the following system is satisfied

[TABLE]

4 Numerical solution of the CDII-SR problem

In this section, we discuss numerical optimization and approximation schemes to solve the CDII-SR problem. We also describe the Picard-type iterative scheme used by the authors in [20] to compare the results of our scheme with the Picard scheme.

4.1 Variable inertial proximal method for CDII-SR

In this context, we first discuss proximal methods that consists of identifying a smooth and a non-smooth part in the reduced objective $\hat{J}(\sigma)$. Thus, we consider the following optimization problem

[TABLE]

We assume that $\nabla_{\sigma}\hat{J}_{1}(\sigma)$, given in (16) is Lipschitz continuous and the upper bound for the Lipschitz constant is obtained using a backtracking search scheme, which will be discussed later. Also, from (8), we have that $\hat{J}_{2}(\sigma)$ is a continuous, convex, and nondifferentiable functional. The formulation of proximal methods depends, essentially, on the following lemma [41]

Lemma 4.1**.**

Let $\hat{J}_{1}$ be differentiable with a Lipschitz continuous gradient with Lipschitz constant $L(\hat{J}_{1})$. Then the following holds

[TABLE]

for all $L\geq L(\hat{J}_{1})>0$.

We note that $L:=L(\hat{J}_{1})$ represents the smallest value of $L$ such that (24) holds true.

In a proximal scheme, one usually minimizes an upper bound of the objective functional at each iteration, instead of minimizing the functional directly. From Lemma 4.1, we obtain the following

[TABLE]

where equality holds if $\sigma=\widetilde{\sigma}$. Furthermore, we have the following equation

[TABLE]

Using the definition of $\hat{J}_{2}(\sigma)=\gamma\|\sigma\|_{L^{1}(\Omega)}$, we have the following lemma from [45] that helps in characterizing the solution of (25).

Lemma 4.2**.**

The following equation holds

[TABLE]

where the left-hand side represents the proximal function and the projected soft thresholding function on the right-hand side is defined as follows

[TABLE]

Using this lemma, the solution to (25) is given by

[TABLE]

This gives rise to the following iterative scheme

[TABLE]

starting from a given $\sigma_{0}$ and is known as the iterative shrinkage-thresholding algorithm (ISTA) scheme [45]. We note that the argument of $\mathbb{S}^{L_{ad}}_{\frac{\gamma}{L}}$ represents a gradient update in a steepest descent scheme with a fixed step size $s=1/L$ in conjunction with a regularized PM filter [41]. Further, to accelerate the ISTA scheme described above, one can consider a sequence $\{t_{k},v_{k}\}$ [45, 46] such that

[TABLE]

and

[TABLE]

This gives us the following update for the optimization variable $\sigma_{k}$

[TABLE]

Replacing $v_{k}$ in (29) with (28), and assuming that ${\nabla}\hat{J}_{1}(\sigma_{k})\approx{\nabla}\hat{J}_{1}(v_{k})$, we obtain the following iterative scheme [46]

[TABLE]

where $\sigma_{-1}=\sigma_{0}$.

The above discussion is valid for any $L\geq L(\hat{J}_{1})$. However, since the quantity $s=1/L$ represents the step size in a gradient update, we use a backtracking line search algorithm to determine the optimal step size in each iteration. This leads to the computation of an upper bound $L_{k}$ that satisfies $L_{k}\geq L(\hat{J}_{1})$ at each iteration step. Thus, we define our variable step size as $s_{k}=1/L_{k}$ and substitute $\tau$ in (26) with $\gamma s_{k}$. The variable step size causes the factor $\frac{(t_{k-1}-1)}{t_{k}}$ in (28) to be non-optimal and we replace it by the fixed inertial parameter $\theta$. This leads to a variable inertial proximal (VIP) scheme, which is described in Algorithm 4.1.

With our VIP scheme, we aim at determining an optimal $\sigma\in L_{ad}\subset H^{1}_{0}(\Omega)$. But in the update step of the algorithm, we have the argument of the thresholding function $\mathbb{S}^{L_{ad}}_{\frac{\gamma}{L}}$ as $\sigma_{k}-s_{k}\nabla_{\sigma}\hat{J}_{1}(\sigma_{k})$. The term $\nabla_{\sigma}\hat{J}_{1}(\sigma_{k})$ is only in $L^{2}(\Omega)$ and the resulting update gives us the argument of $\mathbb{S}^{L_{ad}}_{\frac{\gamma}{L}}$ in $L^{2}(\Omega)$, which is not desired. We, thus, use the $H^{1}$ gradient instead of the $L^{2}$ gradient, which are related by the equation $((\nabla_{\sigma}\hat{J}_{1})_{H^{1}},v)_{H^{1}(\Omega)}=(\nabla_{\sigma}\hat{J}_{1},v)_{L^{2}(\Omega)}$ for all $v\in H^{1}(\Omega)$. But such a $H^{1}$ gradient results in a highly diffused $\sigma$ with blurred edges. We, instead, consider a weighted $H^{1}$ product that represents a suitable denoising of the $\nabla_{\sigma}\hat{J}_{1}(\sigma)$. We apply the denoising operator $R(c)=(I-c\,\Delta)^{-1}$ with a small denoising parameter $c$ and define $(\nabla_{\sigma}\hat{J}_{1})_{H^{1}}=R(c)\,\nabla_{\sigma}\hat{J}_{1}$. Note that a higher value of $c$ results in a greater blurring of the edges along with noise removal. On the other hand, since the PM term in the functional $J$ promotes better resolution with edge-enhancement, we choose the value of $c$ in proportion to the weight $\delta$ of the PM functional term.

We summarize the variable inertial proximal (VIP) scheme for our CDII-SR setup in Algorithm 4.1 below

Algorithm 4.1** (Variable inertial proximal (VIP) method).**

*Input: $\beta$, $\hat{J}_{1}$, $\sigma_{0}=\sigma_{-1}$, $L_{ad}$, $TOL$, $n>1$, *$L_{0}>0$

Initialize:* $E_{0}=1$, $k=0$, choose $\theta\in(0,1)$ and $c_{1}<2$ and $c_{2}>0$;* 2. 2.

While$\|E_{k-1}\|>TOL$* do* 3. 3.

Compute $\nabla_{\sigma}\hat{J}_{1}(\sigma_{k})$ 4. 4.

*Backtracking: Find the smallest nonnegative integer $i$ such that with *

*    $\tilde{L}=n^{i}L_{k-1}$*

[TABLE]

where $\tilde{\sigma}=\mathbb{S}^{L_{ad}}_{\gamma\,s}\left(\sigma_{k}-s\,(\nabla_{\sigma}\hat{J}_{1})_{H^{1}}(\sigma_{k})+\theta(\sigma_{k}-\sigma_{k-1})\right)$, $s=c_{1}(1-\theta)/(\tilde{L}+2c_{2})$, 5. 5.

Set $L_{k}=\tilde{L}$ and $s_{k}=c_{1}(1-\theta)/(L_{k}+2c_{2})$. 6. 6.

$\sigma_{k+1}=\mathbb{S}^{L_{ad}}_{\gamma\,s_{k}}\left(\sigma_{k}-s_{k}\,(\nabla_{\sigma}\hat{J}_{1})_{H^{1}}(\sigma_{k})+\theta(\sigma_{k}-\sigma_{k-1})\right)$** 7. 7.

$\mu_{k}=-\alpha\sigma_{k}-(\nabla_{\sigma}\hat{J}_{1})_{H^{1}}(\sigma_{k})$** 8. 8.

$E_{k}=E(\sigma_{k},\mu_{k})$** 9. 9.

$k=k+1$** 10. 10.

end

In the VIP algorithm, we need to compute the reduced gradient $\nabla_{\sigma}\hat{J}_{1}$. This, in turn, requires an accurate numerical solution of the forward and the corresponding adjoint EIT problems as given in Proposition 6. For the forward EIT equation (4), we use the cell-nodal finite-difference approximation. We consider a sequence of uniform grids $\{\Omega_{h}\}_{h>0}$ given by

[TABLE]

where $N$ represents the number of cells in each direction and $h=\dfrac{(b-a)}{N}$ is the mesh size. The corresponding cell-nodal scheme for (4), at the grid point $(x_{i},y_{j})$, is given as follows

[TABLE]

where $\sigma_{i\pm 1,j}=\sigma(x_{i}\pm h,y_{j})$, $\sigma_{i,j\pm 1}=\sigma(x_{i},y_{j}\pm h)$. The required intermediate values of $\sigma$ are computed as follows

[TABLE]

The Dirichlet boundary data $f_{D}$ is included in the usual way in the right-hand side of the algebraic equation.

For the adjoint equations (17) and (18), we first note that the cell nodal finite difference scheme is not applicable to the right-hand side term in both the equations as they are of the form $G=\nabla\cdot F$, where $F$ is in $L^{2}(\Omega)$. We modify the cell nodal scheme by replacing the nodal value of $G$ at $(x_{i},y_{j})$ with a cell average of $G$ given as follows

[TABLE]

where the cell $C_{ij}$ is defined by

[TABLE]

Since $G=\nabla\cdot F$, using the divergence theorem we have

[TABLE]

The above integral can be approximated with a midpoint quadrature rule along each edge of $C_{ij}$. This results in the following approximation

[TABLE]

where $F=(F^{1},F^{2})$.

4.2 Picard-type algorithm

To solve the CDII inverse problem, the authors in [20] use a Picard-type iterative algorithm. This algorithm was the fundamental approach that was used to solve the CDII inverse problem expressed as the 1-Laplacian equation given in (3) [35, 36] was known as the J-substitution algorithm. To start the Picard algorithm, we use an initial guess for $\sigma$ as $\sigma_{0}$. We then solve (4) for $u_{1}$ in the first iteration and then obtain $\sigma_{1}$ from the interior data (2) corresponding to this value of $u_{1}$. In the next iteration, we solve for $u_{2}$ using $\sigma_{1}$ and again obtain the next iterate $\sigma_{2}$ from the interior data (2) corresponding to this value of $u_{2}$. We apply this method repeatedly, going back and forth between the two sets of interior data, till the difference in the $L^{2}$ norm of two successive iterates of $\sigma$ are small. The algorithm is summarized below.

Algorithm 4.2** (Picard-type algorithm).**

*Input: Initial guess $\sigma_{0}$, maximum iterations = $K$, *$TOL$

Initialize:* $err_{0}=1$, $k=0$.* 2. 2.

While $\|err_{k}\|_{2}>TOL$ and $k<K$ do 3. 3.

$j=(k~{}MOD~{}2)+1$* (To go back and forth between the two interior measurements.* 4. 4.

Solve for $u_{j}^{k}$ from (4), i.e., solve

[TABLE] 5. 5.

Update $\sigma_{k+1}=\dfrac{H_{j}}{|\nabla u_{j}^{k}|}$ 6. 6.

$err_{k+1}=\|\sigma_{k+1}-\sigma_{k}\|_{2}$** 7. 7.

$k=k+1$** 8. 8.

end

The authors demonstrated in their paper that the Picard algorithm was fast and efficient. With 2 data sets and even in the presence of noisy interior data, the Picard method gave good conductivity reconstructions in CDII with a high resolution and less visible artifacts in comparison to other Newton based schemes. Thus, in order to test the performance of our CDII-SR scheme, we will compare the results of our scheme with the results of the Picard scheme in the next section.

5 Numerical experiments

In this section, we validate our CDII-SR framework using different experiments that support the choice of our proposed mathematical formulation and demonstrate its effectiveness in reconstructing a wide variety of objects. We choose the two boundary conditions as $f^{1}_{D}=x,~{}f_{D}^{2}=y$ on $\Gamma$, which is the boundary of $\Omega=(-1,1)\times(-1,1)$. For the experiments, we test our scheme with different values of the regularization parameters $\beta,\gamma,\delta$ in the functional (5), relative to the value of the data-fit parameters $\alpha_{1},\alpha_{2}$. For this purpose, we fix the values of $\alpha_{1}=\alpha_{2}=1.0$. Since the primary objective is to fit the measured interior data, the value of the regularization parameters are chosen to be less than 1. The parameters of our VIP scheme are chosen based on the convergence properties of the scheme given in [46]. We specified the tolerance level as $TOL=10^{-4}$ and maximum number of iterations as 20. But, due to the high non-linearity of the problem, our VIP scheme terminates because of the criteria of maximum number of iterations. For the Picard algorithm, described in Section 4.2, we choose the same stopping criterion. To solve the equation $\mathcal{L}(u_{j}^{k},\sigma_{k},f_{D}^{j})=0$ in Algorithm 4.2, we use the finite difference scheme described in (31).

In all the experiments, the domain $\Omega$ is uniformly discretized into $N=150$ subintervals in both the $x$ and $y$ directions with $h=0.013$. The generation of the synthetic interior electric field data $H^{\delta}$ is done as follows: we first solve for $u$ in (4) with given value of $\sigma$ on a finer mesh with $N=400$ using the finite difference method outlined in Section 4. Then, we compute $\nabla u$ with one-sided finite differences to obtain $H^{\delta}$ on the finer mesh. In the final step, we restrict the obtained $H^{\delta}$ onto the coarser mesh with $N=150$ and choose this as our given data to which we also add noise in some of the experiments.

In Test Case 1, we consider a phantom that is represented by a disk centered at $(0.25,0.25)$ with radius 0.25. The value of $\sigma$ inside the disk is 1 with the background value chosen as 0. The plots of the actual $\sigma$ and the reconstructed $\sigma$ are shown in Figure 1.

Figure 1b shows the reconstruction of $\sigma\in L_{ad}$ without any regularization terms, i.e. $\beta=\gamma=\delta=0$ and no denoising, i.e., $c=0$. Presence of strong artifacts can be observed in this case, which is inherent to the inverse problem and not the algorithm. A study of the pattern of such artifacts are very challenging and is out of the scope of the paper.

Figure 1c shows the result using the CDII-SR scheme without the denoising and the Perona-Malik regularization term, $c=\delta=0$, but with the $L^{2}-L^{1}$ regularization. The $L^{2}-L^{1}$ regularization parameter values were set at $\beta=0.03,~{}\gamma=0.3$. We observe that the artifacts are reduced to some extent, but are still present. Figure 1d shows the reconstruction with the same $L^{2}-L^{1}$ regularization parameter values and the $H^{1}$ denoising parameter value $c=0.001$ but without the PM filter. In this case, we observe that the artifacts diminish by a huge amount, but the edges are more blunt and the value of $\sigma$ is lowered, leading to a loss of resolution and contrast. We correct this loss using the PM regularization term with the parameter value $\delta=0.01$. This can be observed in Figure 1e, where the edges are fairly well seen and the recovered parameter values are very close to the true ones. We also compare our results with the reconstruction in Figure 1f obtained with the Picard algorithm proposed in [20]. We observe a lot of artifacts and a significant loss of contrast in Figure 1f in comparison to the reconstruction shown in Figure 1e, which suggests that our CDII-SR scheme outperforms the Picard scheme.

We also implemented Test Case 1 with different values of the parameters $\theta,c_{1},c_{2}$ in the VIP algorithm along with values of $\beta,\gamma,\delta=0,0.001,0.01,0.1$. The convergence results of the VIP algorithm are listed in Table 1. The variable $v$ represents the value of each one of the parameters $\theta,c_{1},c_{2}$. In the table, each cell in the columns under the individual parameters represent the convergence property of the VIP algorithm for $v$ in the prescribed range, when the other parameters are specified in their range of convergence.

We first remark that varying the parameters in the VIP algorithm results in change of convergence properties of the CDII-SR scheme. But choosing different values of $\theta,c_{1},c_{2}$ in the range of convergence, we did not observe any change in the qualitative or quantitative properties of the reconstructed conductivities. From Table 1, we can see that for values of $\theta\in[0,1),~{}c_{1}\approx 2,~{}c_{2}\ll 1$ we have good convergence of the VIP scheme. For a more detailed discussion of these range of values for convergence, we refer the reader to [46]. Thus, we choose $\theta=0.5,~{}c_{1}=1.9,c_{2}=0.001$ in our subsequent experiments.

In Test Case 2, we consider the heart and lung phantom. It consists of two ellipses representing lungs with the value of $\sigma=1$ and a circular region representing the heart with the value of $\sigma=0.5$. The background value of $\sigma$ is 0. The plots of the actual and reconstructed $\sigma$ for various values of $\beta,\gamma,\delta,c$ are shown in Figure 2.

From Figure 2b, we observe that with a small value of the $L^{1}$ regularization parameter $\gamma=0.01$ in comparison to a higher $L^{2}$ regularization parameter $\beta=0.3$, there are artifacts present in the image. Increasing the value of $\gamma$ to 0.1, keeping $\beta$ fixed results in a decrease of the artifacts as seen in Figure 2c. A further increase in the value of $\gamma$ to 0.3 results in a significant reduction of artifacts as observed in Figure 2d. In Figure 2f, we obtain a reconstruction with high resolution and contrast with a decrease in $\beta$ to 0.03. This is expected because the phantom under consideration in this example is sparse. Thus, with a greater weight for $L^{1}$ regularization term, coupled with the denoising operator controlled by the parameter $c$, described in Sec 4.1, helps in removing artifacts and, thus, promotes high contrast. We further note from Figures 2e and 2g that changing the value of the Perona-Malik regularization term $\delta$ to 0.1 does not visibly change the resolution of the edges. This is because the term $\delta$ is linked to the denoising parameter $c$ that is fixed to be 0.001. For this choice of $c$, choosing $\delta=0.01$ is sufficient to guarantee sharp edges. We further compare our results to the Picard algorithm proposed in [20]. The corresponding reconstruction is shown in Figure 2h. It can be seen that the CDII-SR scheme provides a better contrast image, yet maintaining the same resolution as that of the Picard scheme. Also, there are far more artifacts through the Picard reconstruction method whereas the sparsity assumption and the $H^{1}$ denoising in CDII-SR scheme results in an image with very less artifacts.

Further, to test the robustness of our method, we introduce $10\%$ and $25\%$ multiplicative Gaussian noises in the interior data $H^{\delta}$, which is fed as input to our CDII-SR algorithm. The corresponding reconstructions are shown in Figure 3. We also plot the reconstructions obtained with the Picard algorithm.

We note from Figure 3a that even with $10\%$ noise, we still have a high contrast and high resolution reconstruction, which demonstrates the robustness of our method. When the noise percentage was increased to $25\%$, it led to appearance of some streak artifacts, as seen in Figure 3b. For this purpose, we used a higher value of $L^{1}$ regularization parameter $\gamma=0.5$ and denoising parameter $c=0.01$ to remove artifacts. Due to a higher value of $c$, we increase the value of $\delta=0.1$. We also used the Picard algorithm with the noisy data and the reconstructions showed visible presence of streak artifacts in Figure 3c and 3d.

In Test Case 3, we consider a combination of phantoms, where one is supported on a square annulus $S_{a}=\{(x,y)\in\mathbb{R}^{2}:-0.8<x<-0.7,-0.2<x<-0.1,-0.8<y<-0.7,-0.2<y<-0.1\}$ with $\sigma=3.0$; the other one consists of 2 disks centered at $(0.7,0.7)$ with radius 0.2 and $\sigma=1.0$ and at $(0.55,0.55)$ with radius 0.15 and $\sigma=2.0$. The value of $\sigma$ inside the square annulus has a value -2.0. The plots of the actual and the reconstructed $\sigma$ are shown in Figure 4.

We again observe a similar behavior as in Test Case 2. For the values of $\beta=0.3$ and $\gamma=0.01,0.1$, we obtain reconstructions that have visible artifacts, especially at the corners of the square annulus. This can be seen in Figure 4b and 4c. When $\gamma$ is increased to 0.3 and when $\beta$ is subsequently lowered to 0.03, we have almost no artifacts except for some in front of the smaller disk as seen in Figure 4f and 4g. We also observe that the values inside the annular hole and the inclusion due to the intersection of the two disks are very well reconstructed with nice contrast and high resolution. With the Picard method, we obtain a large number of artifacts and the value inside the hole is also not reconstructed well, thus, leading to loss of contrast as seen in Figure 4h.

We also test the performance of our method in presence of $10\%$ and $25\%$ Gaussian noise in the interior data. The results are shown in Figure 5.

We again observe that our method gives superior quality reconstructions as shown in Figures 5a and 5b. As with Test Case 2, we need greater values of $\beta,\gamma,\delta,c$ to obtain high contrast and high resolution images. From Figures 5c and 5d, we see that, in comparison to the reconstructions with the CDII-SR scheme, the Picard scheme does not provide good reconstructions of objects with holes and inclusions in presence of noise. In particular, the region inside the hole contains a large number of artifacts. This shows the robustness and efficiency of our method in the reocnstructing objects with holes and inclusions.

In this final Test Case 4, we consider another combination of phantoms that result in a lesser sparsity pattern. This experiment is conducted to demonstrate the effect of the $L^{1}$ regularization term in improving contrast. The first phantom is supported on a square annulus $S_{a}=\{(x,y)\in\mathbb{R}^{2}:-0.8<x<-0.7,-0.2<x<-0.1,-0.8<y<-0.7,-0.2<y<-0.1\}$ with $\sigma=3.0$. The value of $\sigma$ inside the square annulus has a value -1.5. The second phantom consists of 4 disks: the first centered at $(0.7,0.7)$ with radius 0.2 and $\sigma=1.0$, the second centered at $(0.55,0.55)$ with radius 0.15 and $\sigma=2.0$, the third centered at $(0,0)$ with radius 0.25 and $\sigma=1.5$ and the last centered at $(0.05,0.6)$ with radius 0.2 and $\sigma=2.5$. The final phantom consists of two heart-shaped objects, represented by cardioids, with the value of $\sigma=4$ in the bigger one and value of $\sigma=3.0$ in the smaller one. The plots of the actual and the reconstructed $\sigma$ are shown in Figure 6.

From Figure 6b, we observe the presence of a large number of artifacts. This is the case when the value of $\gamma=0.01$ is small compared to the value of $\beta$. Increasing the value of $\gamma$ to 0.3 and, subsequently, to 0.5, reduces the number of artifacts as can be seen in Figures 6c and 6d. This shows that the $L^{1}$ regularization indeed leads to a nice contrast in the images, even when the sparsity level is low. Comparing the reconstruction obtained with the Picard algorithm, we see in Figure 6e the presence of a large number of artifacts and loss of contrast, especially in the hole and the inclusions.

6 Conclusions

In this paper, we propose a new framework to facilitate high contrast and high resolution reconstructions in CDII. Our framework is based on formulating the CDII inverse problem as a PDE-constrained optimization problem. In this setup, we minimize an objective functional comprising of least square interior data fitting terms corresponding to two boundary voltage measurements, a $L^{2}-L^{1}$ penalization terms of the log-conductivity that helps promotes sparsity patterns, thus, improving contrast. Additionally we introduce a PM filtering term to sharpen the edges and obtain high resolution images.

We characterized the solution of the optimization problem through an optimality system that was solved using a proximal scheme, coupled with $H^{1}$ denoising to remove artifacts in the reconstructions. We then demonstrated the effectiveness of our proposed scheme through several numerical experiments and compared our results with an existing Picard-type scheme. Our scheme facilitated reconstructions of a wide variety of conductivity patterns with superior contrast and resolution, even in the presence of noise in the data.

7 Acknowledgments

S. Roy was partly supported by the National Cancer Institute, National Institutes of Health, grant number: 1R21CA242933-01.