# A nonlinear weighted anisotropic total variation regularization for
electrical impedance tomography

Yizhuang Song, Yanying Wang and Dong Liu This work was supported by Shandong
Provincial Outstanding Youth Fund (Grand No. ZR2018JL002), National Natural
Science Foundation of China (Grand No. 11501336 and No. 61871356) and the
China Postdoctoral Science Foundation (2019T120604, 2018M630795).
(Corresponding author: Dong Liu.)Yizhuang Song and Yanying Wang are with
School of Mathematics and Statistics, Shandong Normal University, Jinan,
Shandong, 250014, P. R. China. Yizhuang Song is also with Center for Post-
doctoral Studies of Management Science and Engineering, Shandong Normal
University, Jinan, Shandong, 250014, P. R. China (e-mail:
ysong@sdnu.edu.cn).Dong Liu is with CAS Key Laboratory of Microscale Magnetic
Resonance and School of Physical Sciences, University of Science and
Technology of China, Hefei 230026, China, also with School of Biomedical
Engineering and Suzhou Institute for Advanced Research, University of Science
and Technology of China, Suzhou, China, and also with CAS Center for
Excellence in Quantum Information and Quantum Physics, University of Science
and Technology of China, Hefei 230026, China. (e-mail: dong.liu@outlook.com).

###### Abstract

This paper proposes a nonlinear weighted anisotropic total variation (NWATV)
regularization technique for electrical impedance tomography (EIT). The key
idea is to incorporate the internal inhomogeneity information (e.g., edges of
the detected objects) into the EIT reconstruction process, aiming to preserve
the conductivity profiles (to be detected). We study the NWATV image
reconstruction by employing a novel soft thresholding based reformulation
included in the alternating direction method of multipliers (ADMM). To
evaluate the proposed approach, 2D and 3D numerical experiments and human EIT
lung imaging are carried out. It is demonstrated that the properties of the
internal inhomogeneity are well preserved and improved with the proposed
regularization approach, in comparison to traditional total variation (TV) and
recently proposed fidelity embedded regularization approaches. Owing to the
simplicity of the proposed method, the computational cost is significantly
decreased compared with the well established primal-dual algorithm. Meanwhile,
it was found that the proposed regularization method is quite robust to the
measurement noise, which is one of the main uncertainties in EIT.

###### Index Terms:

Electrical impedance tomography, anisotropic total variation, regularization,
lung imaging.

##  I Introduction

Electrical Impedance Tomography (EIT) aims to reconstruct the (change of)
conductivity distribution inside objects by injecting a current and measuring
the voltage responses through pairs of surface electrodes mounted on the
object. EIT has the advantages of being noninvasive, portable, low cost,
capable of high temporal resolution, long duration and continuously
monitoring, and much more. These advantages make EIT useful for bedside
medical apparatus in clinical applications. For this, EIT was commercialized
and introduced in medical applications since 1980s [1]. However, EIT has not
yet been widely used in routine clinical applications due to the fact that it
is a diffusive modality.

The EIT reconstruction process is commonly recasted into a (least square
based) data-fitting inverse problem between the boundary measurement and the
computational data. To deal with the ill-posedness, regularization techniques
are widely added to the data-fitting to attract the solution satisfying
applications driven constraints. Depending on the form of different
constraints, the regularization methods can be roughly classified into three
categories: projection based regularization, a priori conductivity based
penalization and learning based regularization.

For projection based regularization, typical examples are truncated singular
value decomposition (tSVD) method [3, 2] and principle component analysis
(PCA) method [4]. Even though these methods are capable of providing stable
reconstructions, they generally produce ringing artifacts due to the dropout
of a certain frequency components [5].

In the case of penalty based regularization, typical examples include Tikhonov
regularization [6, 4] and its multiplicative form [7], monotonicity based
regularization[8], factorization based regularization [5], sparsity based
regularization [9, 10, 11]. These methods are able to provide stable
reconstructions at the cost of blurring the edges of internal inhomogeneities.
Total variation (TV) based regularization [12, 13] and its variants [14, 15,
16] have the advantage of preserving the discontinuities of the internal
structures especially for dealing with the cases of piecewise constant
conductivity distributions. Since TV regularization is non-differentiable, the
so called primal-dual algorithm [13] is usually used to deal with the non-
differentiability. However, this method needs to handle two optimization
problems and hence, it is time consuming [17] and could produce pseudo-edges
and lead staircase effects in the reconstruction. Meanwhile, the Split Bregman
method and the first-order TV regularization [18] are also used at the cost of
decreasing the effect of edge preserving [19]. Lee et al proposed a so-called
fidelity embedded regularization (FER) technique [20]. Using this method, high
quality geometries of the internal inhomogeneities can be obtained. However,
since the regularization in the method does not depends on the internal
structures the accuracy of values of the reconstructed conductivity, which
could be a useful information for clinical use, can not be guaranteed. In
actuality, numerical simulations show that when the regularization parameter
was set to be infinity, the estimated conductivity is usually far away from
the true value (see Section V for details). Recently, regularization based on
manifold learning [21] has been published using the results in [20] as the
training data. Since EIT has not been widely used in clinical applications, it
is difficult to obtain enough training data to improve the performance
continuously.

In this work, we proposed a nonlinear weighted anisotropic TV regularization
approach for EIT reconstruction problem. In comparison to the well established
isotropic TV regularization methods, the nonlinear weighted anisotropic
regularization takes use of a nonlinear weight related to the unknown
conductivity to eliminate the pseudo edges of the internal structures.
Moreover, a soft thresholding formula is derived in the ADMM [22] type
algorithm to accelerate the reconstruction. Comparing with the anisotropic TV
regularization method [13, 18], the proposed regularization approach utilizes
nonlinear weight to pull back the internal edges from a possible distortions
along the coordinate axes [23] and provides more accurate EIT reconstructions.
To contextualize the proposed regularization among contemporaries, Table I
provides a comparison of the pros and cons associated with the existing EIT
regularization methods in the field of EIT. To validate the proposed method,
we conduct 2D and 3D numerical simulations and human EIT lung imaging
experiment. These experiments are carried out to illustrate the main
advantages of the proposed NWATV based method with respect to the existing TV,
first-order TV and FER regularization methods.

The remaining sections of this paper are organized as follows. In Section II
we review the forward and inverse problems in EIT. Following we describe the
proposed nonlinear weighted anisotropic TV regularization approach in Section
III and provide the reconstruction algorithm in Section IV. Next we provide
the numerical and human lung experiments in Section V. Finally, we provide a
discussion in Section VI and conclude the paper in Section VII.

TABLE I: Comparisons of the pros and cons of different regularization schemes
in EIT

Categories |  Name of regularizer |  Priors |  Pros |  Cons
---|---|---|---|---
Projection method |  tSVD method [3, 2] |  Singular value of $\mathbf{S}>$a threshold |  Easy to implement. |  There is ringing artifacts abandoning the basis corresponding to the small singular value [5].
Principle component analysis (PCA) [4] |  $\sigma$ lies in the solution space $S_{\omega}$ |  Apply geometrical information from MRI and CT etc. |  Need a learning set obtained from the other modalities such as CT, MRI or the reference conductivity value which may not be available. Instability could occur when the priors are not from the patient.
_a priori_ conductivity based penalty |  Tikhonov [6, 4] |  $\|\sigma\|_{L^{2}}$ |  The optimization problem is convex and differentiable. |  Reconstruction results depends on an artificially chosen regularization parameter. Inclusions is blurred.
Multiplicative Tikholnov regularization [7] |  $\left\|\frac{|\nabla\sigma|^{2}+\delta_{n-1}^{2}}{|\nabla\sigma^{n-1}|^{2}+\delta_{n-1}^{2}}\right\|_{L^{2}}$ for measurement data and mesh dependent $\delta_{n-1}$ |  Reconstruction results are independent of the regularization parameter choice. |  A weighted $L^{2}$ regularizer could still blur the image. In each step, numerical integral and differentiation are needed to calculate. Moreover, Gauss-Newton method with step size using line search method is used to solve a nonquadratic cost functional which is time consuming.
Total variational (TV) regularization [13] |  $\int_{\Omega}|\nabla\sigma|$ |  Can preserve the edge of internal inhomogeneities |  Staircase effect and pseudo edge exists. Reconstructing is slow and resolution is high when the prime-dual algorithm is used while resolution is low and reconstructions is fast when the Split-Bregman method is used to minimized the TV constrained optimal problem [19].
TGV [11, 24] |  $TGV_{G}^{\alpha}(\sigma)$ |  Efficiency in removing the staircase effect when using TV for piecewise linear conductivity distribution. |  In some imaging cases such as human lung imaging, piecewise linear is not a reasonable assumption. Balancing efficiency and quality is difficult as that in TV regularization.
First order TV regularization [18, 23] |  $\|\nabla\sigma\|_{L^{1}(\Omega)}$ |  The imaging speed is high |  The capacity of edge preserving is lower than TV. Distort the inclusion boundaries along the coordinate axes.
Monotonicity based method [8] |  The sign of $\frac{\partial\sigma}{\partial t}$ coincides with the breathing process. |  No confusion between breathing and expirtion |  Still a kind of $L^{2}$ regularization which could blur the internal edges.
Factorization based method [5] |  Highlight the right singular vector associate with high singular values. |  Alleviated Gibbs ringing artifacts |  Heuristical argument without strict mathematical theory
Wavelet frame TGV [16] |  $TGV+\|{\bm{\lambda}}\cdot W\sigma_{recon}\|_{0}$ |  Using the $L^{0}$ norm of the wavelet frame to sharpen the TGV reconstructed conductivity |  Primal-dual method is used for the TGV minimization, reconstruction is time consuming. Postprocessing the TGV reconstructions needs more time.
Fidelity-imbedded regularization [20] |  $\sqrt{\sum\limits_{l}|\langle S_{k},S_{l}\rangle|}$ |  Reconstruction results are independent of the regularization parameter choice and fast imaging due to only a direct algebraic inversion is needed to compute. |  The accuracy of the reconstruction can not be guaranteed due to the lack of mathematical background.
Learning method |  Manifold learning method [21] |  $\sigma\in$ a manifold $\mathcal{M}$. |  Imaging quality is higher than the model based regularizer. |  To obtain experimental labelled data $(\\{\mathbf{V}_{n},\sigma_{n})\\}$ is difficult before designing a high quality model based regularization method.
|  Proposed method |  $\left\|\frac{\nabla\sigma}{|\nabla\sigma|^{2}}\right\|_{L^{1}}$ |  Improved preservation of the internal inhomogeneity edges; the computational time is significantly decreased; robust to the additional Gaussian random noise. |  Need to select proper regularization parameters manually.

##  II Forward and inverse problems in EIT

Let $\Omega\subset\mathbf{R}^{n}$ ($n=2,3$) be the imaged object with a smooth
boundary. Let $\gamma=\sigma+i\omega\epsilon$ represent the admittivity
distribution of the region $\Omega$. Here, the conductivity $\sigma$ and
permittivity $\epsilon$ depend on position $\mathbf{r}=(x,y,z)$, angular
frequency $\omega$.

For an $E$-channel EIT system, $E$ electrodes
$\mathcal{E}_{1},\mathcal{E}_{2},\cdots,\mathcal{E}_{E}$ are attached on
$\partial\Omega$, the boundary of $\Omega$. We inject a series of time
harmonic currents with magnitude $I$ mA following e.g. the neighboring
protocol. Under such a protocol, we sequentially inject several currents using
pairs of electrodes $({\mathcal{E}}_{j},{\mathcal{E}}_{j+1})$ for $j=1,...,E$,
where we assign $E+1$ to be $1$. The induced electric potential $u^{j}$ is
governed by the following elliptic partial differential equation (PDE) with
mixed boundary conditions:

$\left\\{\begin{split}&\nabla\cdot(\sigma\nabla u^{j})=0,\quad\mbox{ in
}\Omega\\\ &\sigma\nabla u^{j}\cdot\mathbf{n}=0,\quad\mbox{ on
}{\partial\Omega\setminus\mathop{\cup}\limits_{i=1}^{E}{\mathcal{E}}_{i}}\\\
&\int_{\mathcal{E}_{j}}\sigma\frac{\partial u^{j}}{\partial{\bf
n}}dS=I=-\int_{\mathcal{E}_{j+1}}\sigma\frac{\partial u^{j}}{\partial{\bf
n}}dS\\\ &\left.u^{j}+z_{i}\sigma\frac{\partial u_{j}}{\partial{\bf
n}}\right|_{{\mathcal{E}}_{i}}=U^{j}_{i},\qquad i,j=1,2,\cdots,E\\\
&\int_{\mathcal{E}_{k}}\sigma\frac{\partial u^{j}}{\partial{\bf n}}dS=0,\quad
k\in\\{1,2,\cdots,E\\}\setminus\\{j,j+1\\}.\end{split}\right.$ (1)

Here, $z_{i}$ is the contact impedance between the electrode and
$\partial\Omega$ and $U_{i}^{j}$ is the voltage potential on
${\mathcal{E}}_{i}$ caused by the injection currents using the electrode pair
$({\mathcal{E}}_{j},{\mathcal{E}}_{j+1})$. Using an EIT measurement device,
the quantity $V_{i}^{j}[\sigma]=U_{i}^{j}[\sigma]-U_{i+1}^{j}[\sigma]$ is
measurable. Note that $z_{i}$ is unknown, however, for
$i\in\\{1,2,\cdots,E\\}\setminus\\{j-1,j,j+1\\}$, the contact impedance can be
neglected [20] and hence $V_{i}^{j}\approx
u^{j}|_{{\mathcal{E}}_{i}}-u^{j}|_{{\mathcal{E}}_{i+1}}$. To reconstruct
$\sigma$, EIT uses the following reciprocity principle

$V_{i}^{j}[\sigma]\approx
u^{j}|_{{\mathcal{E}}_{i}}-u^{j}|_{{\mathcal{E}}_{i+1}}=\frac{1}{I}\int_{\Omega}\sigma(\mathbf{r})\nabla
u^{i}(\mathbf{r})\cdot\nabla u^{j}(\mathbf{r})d\mathbf{r}.$ (2)

EIT uses $E\times(E-3)$ data $V_{i}^{j}$ and the relation (2) to reconstruct
the unknown quantity $\sigma$. However, the above problem is nonlinear and
ill-posed. Linearization is usually applied to deal with the non-linearity and
regularization is used to handle the ill-posedness.

To be precise, we assume that $\sigma$ is a perturbation of a constant
$\sigma_{0}$, that is $\sigma(\mathbf{r})=\sigma_{0}+\delta\sigma(\mathbf{r})$
for $\mathbf{r}\in\Omega$. Here, $\delta\sigma$ has compact support in
$\Omega$. Then $\delta\sigma$ approximately satisfies the following integral
equation

$\int_{\Omega}\delta\sigma(\mathbf{r})\nabla u^{i}_{0}(\mathbf{r})\cdot\nabla
u^{j}_{0}(\mathbf{r})d\mathbf{r}=I(V_{i}^{j}[\sigma]-V_{i}^{j}[\sigma_{0}])=:I\delta
V_{i}^{j}[\sigma],$ (3)

where $u^{i}_{0}$ is the solution of the equation (1) with $\sigma$ replaced
by $\sigma_{0}$.

Discretize the domain $\Omega$ as $\Omega=\cup_{k=1}^{N}T_{k}$, where $T_{k}$
is triangular element, then $\sigma|_{T_{k}}$ can be considered as a constant
$\sigma_{k}$. $\sigma$ can be approximated by the vector
$\bm{\sigma}:=(\sigma_{1},\sigma_{2},\cdots,\sigma_{N})^{T}$. Then we convert
the formula (3) to the following system

$\mathbf{S}{\delta\bm{\sigma}}=\delta\mathbf{V}$ (4)

where $\mathbf{S}$ is an $E(E-3)\times N$ _sensitivity matrix_ whose element
$S_{pq}=\frac{1}{I}\int_{T_{q}}\nabla u_{0}^{m}(\mathbf{r})\cdot\nabla
u_{0}^{n}(\mathbf{r})d\mathbf{r}$ , where $m=[\frac{p-1}{N-3}]+1$,
$n={p-12[\frac{p-1}{N-3}]+2}$, $\delta\mathbf{V}\in\mathbf{R}^{E(E-3)}$ whose
$q$-th element is $V_{q}=\delta V_{m}^{n}$ and the operator $[\cdot]$
represents the largest integer less than $\cdot$.

##  III Nonlinear weighted anisotropic TV regularization

To solve the equation (4), we reformulate it to the following least squares
problem

$\delta\bm{\sigma}^{*}=\mathop{{\rm arg\
min}}\limits_{\delta\bm{\sigma}}\|\mathbf{S}\delta\bm{\sigma}-\delta\mathbf{V}\|^{2}_{L^{2}(\Omega)}.$
(5)

Let $\mathbf{s}_{l}$ denotes the $l$-th column vector of $\mathbf{S}$ for
$l=1,2,\cdots,N$. The correlation between the $l$ and $k$-the column vectors
$\langle\frac{\mathbf{s}_{l}}{|\mathbf{s}_{l}|},\frac{\mathbf{s}_{k}}{|\mathbf{s}_{k}|}\rangle$
decreases rapidly as the distance between $T_{l}$ and $T_{k}$ increases [20];
this makes the condition number of the matrix $\mathbf{S}^{T}\mathbf{S}$ to be
approximately infinity. Hence, the minimization problem (5) is unstable
against the measurement errors and noises in $\delta\mathbf{V}$. We use the
following minimization problem to approximate (5),

$\delta\bm{\sigma}^{*}_{\lambda}=\mathop{{\rm arg\
min}}\limits_{\delta\bm{\sigma}}\left\\{\frac{1}{2}\|\mathbf{S}\delta\bm{\sigma}-\delta\mathbf{V}\|^{2}_{L^{2}(\Omega)}+\lambda\mbox{Reg}(\delta\bm{\sigma})\right\\}.$
(6)

Here, instead of minimizing the term
$\|\mathbf{S}\delta\bm{\sigma}-\delta\mathbf{V}\|^{2}_{L^{2}(\Omega)}$, we
also minimize the other regularization term $\mbox{Reg}[\delta{\bm{\sigma}}]$,
where $\lambda$ is a parameter balancing the fidelity term
$\|\mathbf{S}\delta\bm{\sigma}-\delta\mathbf{V}\|^{2}_{L^{2}(\Omega)}$ and the
regularization term. The choice of regularization depends on the a priori
information of the conductivity.

Following, for ease of explanation we write
$\Omega=\cup_{i=1}^{n_{I}}D_{i}\cup(\Omega\setminus\overline{\cup_{i=1}^{n_{I}}D_{i}})$,
where $D_{i}$ ($i=1,2,\cdots,n_{I}$) represent $n_{I}$ internal structures
with $n_{I}$ a positive integer. Notice that
$\zeta(\mathbf{r}):=1/|\nabla\sigma(\mathbf{r})|^{2}$ can be considered as
indicators of edges of conductivity inhomogeneities. To be precise, as ${\bf
r}\to\partial D_{i}$, $\zeta(\mathbf{r})\to 0$ and
$\zeta(\mathbf{r})\approx+\infty$ for $\mathbf{r}$ away from $\partial D_{i}$
($i=1,2,\cdots,n_{I}$). Moreover, $|\nabla\sigma|$ is relatively sparse in
lung imaging due to the piecewise constant structure of $\sigma$ at a given
time. Based on the above observations, we construct the nonlinear weighted
regularization term $\|\nabla\sigma/|\nabla\sigma|^{2}\|_{L^{1}(\Omega)}$.
Hence, (6) is changed to

$\delta\bm{\sigma}^{*}_{\lambda}=\mathop{{\rm arg\
min}}\limits_{\delta\bm{\sigma}}\left\\{\frac{1}{2}\|\mathbf{S}\delta\bm{\sigma}-\delta\bm{V}\|^{2}+\lambda\|\mathbf{p}\cdot\mathbf{D}(\delta\bm{\sigma})\|_{l_{1}}\right\\},$
(7)

Here,
$\mathbf{p}=(\zeta(\delta\bm{\sigma});\zeta(\delta\bm{\sigma}))\in\mathbf{R}^{2N}$,
$\mathbf{D}(\delta\bm{\sigma})=(\mathbf{D}_{x}{\delta\bm{\sigma}};\mathbf{D}_{y}{\delta\bm{\sigma}})\in\mathbf{R}^{2N}$,
where $\mathbf{D}_{x}$, $\mathbf{D}_{y}\in\mathbf{R}^{N\times N}$ are
respectively the first order difference operators along the $x$ and $y$
directions.

Given the fact that the minimization problem (7) is non-differentiable, we
adopt the well established Alternating Direction Method of Multipliers (ADMM)
[22] to solve it. The initial guesses
$(\mathbf{z}_{0},\mathbf{p}_{0},\mathbf{y}_{0})$, $\delta\bm{\sigma}$ are
updated via the following schemes,

$\displaystyle\delta\bm{\sigma}_{n+1}$ $\displaystyle=\mathop{{\rm arg\
min}}\limits_{\delta\bm{\sigma}}\mathcal{L}_{\rho}(\delta\bm{\sigma},\mathbf{z}_{n},\mathbf{p}_{n};\mathbf{y}_{n});$
(8a) $\displaystyle\mathbf{z}_{n+1}$ $\displaystyle=\mathop{{\rm arg\
min}}\limits_{\mathbf{z}}\mathcal{L}_{\rho}(\delta\bm{\sigma}_{n+1},\mathbf{z},\mathbf{p}_{n};\mathbf{y}_{n});$
(8b) $\displaystyle\mathbf{p}_{n+1}$
$\displaystyle=(\zeta(\delta\bm{\sigma}_{n+1});\zeta(\delta\bm{\sigma}_{n+1}));$
(8c) $\displaystyle\mathbf{y}_{n+1}$
$\displaystyle=\mathbf{y}_{n}+\rho(\mathbf{D}\delta\bm{\sigma}_{n+1}-\mathbf{z}_{n+1}).$
(8d)

Here, $\mathbf{z}\in\mathbf{R}^{2N}$ is an auxiliary variable,
$\mathcal{L}_{\rho}$ is the augmented Lagrangian functional defined as

$\begin{split}\mathcal{L}_{\rho}(\delta\bm{\sigma},\mathbf{z},\mathbf{p};\mathbf{y}):=\frac{1}{2}\|\mathbf{S}\delta\bm{\sigma}-\delta\mathbf{V}\|^{2}+\lambda\|\mathbf{p}\cdot\mathbf{z}\|_{{l_{1}}}\\\
+\frac{\rho}{2}\|\mathbf{D}\delta\bm{\sigma}-\mathbf{z}\|_{l_{2}}^{2}+\mathbf{y}^{T}(\mathbf{D}\delta\bm{\sigma}-\mathbf{z}),\end{split}$

where $\mathbf{y}$ is the Lagrangian multiplier and $\rho$ is a scalar penalty
parameter.

Meanwhile, from (8a) and (8b) we can derive

$\delta\bm{\sigma}_{n+1}=(\frac{1}{\rho}\mathbf{S}^{T}\mathbf{S}+\mathbf{D}^{T}\mathbf{D})^{-1}(\frac{1}{\rho}\mathbf{S}^{T}\delta\mathbf{V}+\mathbf{D}^{T}\mathbf{z}_{n}-\mathbf{D}^{T}\frac{\mathbf{y}_{n}}{\rho}),$
(9)

and

$\mathbf{z}_{n+1}[k]=h_{\frac{\lambda\mathbf{p}_{n}[k]}{\rho}}(\mathbf{D}\delta\bm{\sigma}_{n+1}[k]+\frac{\mathbf{y}_{n}[k]}{\rho}).$
(10)

Here, $\mathbf{z}_{n+1}[k]$ is the $kth$ element of $\mathbf{z}_{n+1}$ for
$k=1,2,\cdots,2N$, $h_{g}(\cdot)$ is the soft thresholding formula defined as
[26]

$h_{g}(\cdot)=\left\\{\begin{array}[]{ll}\cdot-g\mbox{sgn}(\cdot),&|\cdot|>g,\\\
0,&\mbox{otherwise},\end{array}\right.$ (11)

where sgn is the sign function. We provide a proof of formula (10) in the
Appendix section.

##  IV Reconstruction algorithm

In this section, we summarize the reconstruction algorithm utilizing the
proposed nonlinear weighted anisotropic TV regularization. Note that in the
human experiment there is unavoidable noise and artifacts in the measured
data, we first pre-process the data as that in [20] to reduce the artifacts
caused by rib cage movement[30].

###  IV-A Data pre-processing and modeling error blocking

Instead of using the direct measurement $\delta\mathbf{V}$, we use the pre-
processed data $\delta\mathbf{V}^{\diamond}$. Here,

$\delta\mathbf{V}^{\diamond}=\delta\mathbf{V}-\mathbf{S}_{bdy}(\mathbf{S}^{T}_{bdy}\mathbf{S}_{bdy}+\lambda_{b}\mathbf{I})^{-1}\mathbf{S}_{bdy}^{T}\delta\mathbf{V},$

where $\mathbf{S}_{bdy}$ is a sub-matrix of $\mathbf{S}$ including the columns
related with boundary elements, $\lambda_{b}$ is a regularization parameter
and $\mathbf{I}$ is the identity matrix.

Due to the reconstruction problem is iteratively solved, the modeling error
could propagate in the forward problem, especially the estimated conductivity
near the electrodes during the iterations. The modeling error propagation
could heavily deteriorate the accuracy of the reconstructed images. To this
end, we assume that the conductivity near the electrodes is invariant with
respect to patients breathing and it only changes in the domain
$\widetilde{\Omega}$, where $\widetilde{\Omega}\subset\subset\Omega$
represents the domain including two lungs. For each iteration, we modify
$\delta{\bm{\sigma}}_{n+1}$ by

$\delta{\bm{\sigma}}_{n+1}\leftarrow\delta{\bm{\sigma}}_{n+1}\chi_{\widetilde{\Omega}},$
(12)

where $\chi_{\widetilde{\Omega}}$ is the characteristic function of the domain
$\widetilde{\Omega}$. Subsequently, we will nominate modified nonlinear
weighted anisotropic TV (MNWATV) as the proposed regularization using the
preprocessed data.

###  IV-B Peudocode for the proposed method

We summarize the above procedures as the following algorithm in the form of
pseudocode.

Algorithm 1 The proposed algorithm

0: Measured data $\mathbf{V}$, the maximum iteration number $M$, the
parameters $\lambda$, $\rho$, tolerance $tol$ and the sensitivity matrix
$\mathbf{S}$

0: conductivity distribution $\bm{\sigma}$

1: Initialize: $\mathbf{y}_{0}=\mathbf{0}$, $\mathbf{z}_{0}=\mathbf{0}$,
$\mathbf{p}_{0}=\mathbf{1}$, $n=1$.

2: while $n<M$ do

3: Update $\bm{\sigma}_{n+1}$ using (9) and set
$\bm{\sigma}_{n+1}\leftarrow\bm{\sigma}_{n+1}\chi_{\widetilde{\Omega}}$.

4: Update $\mathbf{z_{n+1}}$ using (10).

5: Update $\mathbf{p}_{n+1}$ using (8c).

6: Update $\mathbf{y}_{n+1}$ using (8d).

7: if $\|\bm{\sigma}_{n+1}-\bm{\sigma}_{n}\|<tol$ then

8: break

9: end if

10: end while

##  V Numerical and Human Experiments

To test the performance of the proposed algorithm, we carried out numerical
examples and human EIT lung imaging experiment.

###  V-A Experiments Setup

In the experiments, we attach 16 electrodes to the boundary of the subject as
uniform as possible. In all the experiments, we inject currents with an
amplitude of 1.0 mA and measure the voltages between adjacent pairs of
electrodes. To avoid having $0$ in the denominator, in all the experiments we
approximate $\zeta(\delta\bm{\sigma}_{n+1})$ by
$\frac{1}{|\mathbf{D}\delta\bm{\sigma}_{n+1}|^{2}+\delta}$ for a small
$\delta>0$.

The results of the experiments are compared to methods employing TV, first-
order TV and FER regularizers. For the comparison, in the numerical
experiments, we computed the relative error (RE), the peak signal-to-noise
ratio (PSNR) to evaluate the accuracy of reconstructed images. At the $n$-th
step, $RE(n)$ is defined as

$RE(n)=\frac{\|\bm{\sigma}_{n}-\bm{\sigma}^{\ast}\|_{2}}{\|\bm{\sigma}^{\ast}\|_{2}},$

and $PSNR(n)$ is defined as

$PSNR(n)=10\log_{10}\frac{\max(\bm{\sigma}_{n}\odot\bm{\sigma}_{n})}{MSE(n)}.$

Here $\cdot\odot\cdot$ represents the component wise multiplication.
$\bm{\sigma}^{\ast}$ represents the true conductivity and
$MSE(n)=\frac{1}{N^{2}}\sum_{i=1}^{N^{2}}\left[\bm{\sigma}_{n}(i)-\bm{\sigma}^{\ast}(i)\right]^{1/2}$
is the mean square error between $\bm{\sigma}^{\ast}$ and $\bm{\sigma}_{n}$.

Since we do not know the true conductivity in the human lung imaging
experiment, to define the $\widetilde{RE}(n)$ and $\widetilde{PSNR}(n)$,
$\sigma^{*}$ was replaced by $\sigma_{TV}$, which is the reconstructed
conductivity using the TV regularization. Additionally, we also compare the
CPU time in the reconstruction process. The reconstructions are carried out
using a workstation with 2.00 GHz Inter (R) Xeon (R) Gold 6138 CPU, 256 GB
memory, Windows 10 operating system.

To decrease the computational cost, we use the point electrodes at the center
of the true electrodes instead of using the true circular electrodes. [29]
provide a theoretical result towards the reasonability of such an
approximation. We present the reconstruction results in the image with
$256\times 256$ pixels. Moreover, in all the experiments, we set
$tol=10^{-5}$.

###  V-B 2D numerical experiments

We construct a disk with radius of $0.1$ m centered at the original point
$(0,0)$. We put two ellipses inside the disk to simulate the lungs. The center
of the constructed ellipses are respectively $(0.04,-0.01)$ and
$(-0.04,-0.01)$. The size of the ellipses vary to form 10 different models
(see the first row of Fig. I). For the $k$-th model ($k=1,2,\cdots,10$), the
major and minor semi-axes of the left ellipse are respectively
$\mathbf{a}_{l}=(0.012+0.001k,0.024+0.002k)$,
$\mathbf{b}_{l}=(-0.012-0.001k,0.006+0.0005k)$ while the values of the right
ellipse are respectively $\mathbf{a}_{r}=(-0.012-0.001k,0.024+0.002k)$,
$\mathbf{b}_{r}=(0.012+0.001k,0.006+0.0005k)$. We set the conductivity of
background and inclusion to be 1.0 S/m and 1.1 S/m respectively. One of the 10
models is shown in Fig. 1 (a). In the inversion, a finite element mesh with
1024 triangular elements and 545 nodes as shown in Fig. 1 (b) is used.

(a)

(b)

Figure 1: (a) The 2D numerical model and its conductivity distribution. (b)
Nodes and elements used in the inversion process.

We first solve the two-dimensional forward problem (1) to obtain the boundary
voltage datum. Using this datum we reconstruct the conductivity images using
the proposed algorithm and three existing methods. Gaussian random noise with
SNR 50dB is added to the data to test the anti-noise performance. We set
$\lambda=5\times 10^{-13},\rho=1\times 10^{-10},\delta=0.01,M=20$ for the
proposed method. We explain the selected parameters in Section VI. For the FER
method, we only consider the case when the regularization parameter is set to
be $\infty$. The reference conductivity $\sigma_{0}$ is set to be background
conductivity 1.0. The parameters of the other regularization methods are
empirically chosen.The results are shown in Fig. I.

Figure 2: The reconstruction results for 2D model using four different
regularizers. The first row depicts the ground truth conductivity
distribution, the second to the fifth row illustrate the reconstruction
results using the FER, the TV, the first-order TV and the proposed
regularizers.

Fig. 3 (a)-(c) illustrates the profiles along the solid, the dash and the dot
dash lines respectively shown in the 7th model in Fig. I. Fig. 4 (a)-(b) shows
the behavior of $RE(n)$ and $PSNR(n)$ as $n$ increases for each regularizers
except FER since it is a direct reconstruction method. In Table II we compare
the computational time for the reconstructions using the aforementioned
regularizers.

(a)

(b)

(c)

Figure 3: Profiles of the reconstructed conductivity along the solid line (a),
the slash line (b) and the dot dash line (c) shown in Fig. I. The meaning of
each colorful lines in (a)-(c) are shown in each subfigure.

(a)

(b)

Figure 4: Behaviors of $RE(n)$ (a) and $PSNR(n)$ (b) of 2D simulation.
Meanings of each colorful lines are shown in the subfigure.

From the comparisons, the other three regularizers can produce more accurate
images than FER does. Moreover, for the proposed method, the computational
time is significantly reduced while maintaining a satisfactory accuracy. On
the other hand the accuracy of the first-order TV regularizer is obviously
decreased.

###  V-C 3D numerical experiment

We do the numerical simulation using a 3D adult human thorax model in EIDORS
[28]. The geometry of the model is shown in Fig. 5 (a). We attach 16 circular
electrodes along the boundary of the center slice. The geometry is discretized
by 5,379 nodes 26,358 tetrahedrons. To test the algorithm, we set a
conductivity contrast as follows: the conductivity of the background is set to
be 1.0 S/m while the values in the two lung shape anomalies are set to be 0.8
S/m. Fig. 5 (b) illustrates the true conductivity distribution at the center
slice. Fig. 5 (c) provide the discretization information of the reconstruction
which includes 637 nodes and 1,172 triangles.

(a)

(b)

(c)

Figure 5: The 3D adult human thorax model. LABEL:sub@Side_view The geometry
and internal structure. LABEL:sub@conductivity_distribution Ground truth
conductivity distribution at the center slice. LABEL:sub@Imaging_plane Imaging
plane.

We first solve the three-dimensional forward problem (1) to obtain the
boundary voltage datum. Using this datum we reconstruct the conductivity
images using the proposed algorithm and three existing methods. The results
are shown in Fig. 6. In Fig. 6, (a) shows the reconstruction using TV
regularization and primal-dual algorithm at the 10th step, (b) depicts the
results of FER for the regularization parameter set to be $\infty$, (c)
illustrates the reconstructions using the first-order TV regularizer and the
ADMM method for the 10th step and (d) is the results using the proposed method
at the 10th step. The regularization parameters for TV and first-order TV are
selected manually as optimal as possible. The regularization parameters for
the proposed method are set to be $\lambda=1\times 10^{-10},\rho=1\times
10^{-9},\delta=0.06$ and $M=10$. For the comparison, in Fig. 7 (a)-(c) we
respectively plot the profiles along the solid, the dash and the dot dash
lines respectively shown in Fig. 6 (d). Fig. 8 (a)-(b) show the behavior of
$RE(n)$ and $PSNR(n)$ as $n$ increases for the methods excluding FER. In Table
II we compare the computational time for the reconstructions using different
regularizers.

From the comparisons, we obtained similar conclusions as the 2D case does.

(a)

(b)

(c)

(d)

Figure 6: The reconstruction results using four different regularizers. (a) is
the reconstruction using the primal-dual algorithm for TV regularizer at the
10th steps. (b) depicts the result of FER for the regularization parameter set
to be $\infty$. (c) is the reconstructions using the first-order TV
regularizer at the 10th iteration. (d)is the results using the NWATV at the
10th step.

(a)

(b)

(c)

Figure 7: Profiles of the reconstructed conductivity along the solid line (a),
the slash line (b) and the dot dash line (c) shown in Fig. 6 (d). The meaning
of each colorful lines in (a)-(c) are shown in each subfigure.

(a)

(b)

Figure 8: Behaviors of $RE(n)$ (a) and $PSNR(n)$ for the 3D simulation. Meanings of each colorful lines are shown in the subfigure. TABLE II: Computational times of the the 7th model in Fig. I for 2D simulation and the 3D simulations. Dimension | FER | TV | First-order TV | NWFOTV
---|---|---|---|---
2D | 0.029 | 2.311 | 0.644 | 0.629
3D | 0.028 | 1.733 | 0.398 | 0.428

###  V-D EIT lung imaging experiment

The human lung EIT experiment is approved by the ethics committee of the
science and technology division, Shandong Normal University. In this
experiment, we attach 16 electrodes [31] around the exterior of the object’s
torso as shown in Fig. 9 (a). We choose the disposable Ag/AgCl
electrocardiogram (ECG) electrodes with the size of $20\times 27$ mm. We use
the 3D scanner [32] to scan the body with electrodes and obtain the geometry
and the positions of the electrodes as accurate as possible. The geometry
obtained is shown in Fig. 9 (b). To obtain the position of the electrodes, we
use point electrode approximations which lie on the center of the true
electrodes [29]. Then we discretize the imaging slice with point electrodes
into 614 nodes and 1132 triangles. The discretization result is shown in Fig.
9 (c). During the reconstruction, we assign the reference conductivity
$\sigma_{0}=0.0107$ S/m [28] as the conductivity of the lungs at the state of
the end expiration and use this $\sigma_{0}$ to obtain the sensitivity matrix
$\mathbf{S}$ in (4).

(a)

(b)

(c)

Figure 9: (a) Experiment setup, (b) the geometry obtained from 3D scanner, (c)
a 2D imaging slice with point electrodes approximations.

The human EIT data was collected with the Sciospec 16-channel EIT system [33].
The frequency of the injection current is 10 kHz and the speed of the data
acquisition is 30 frames/s. Using the datum obtained from the EIT scanner, we
reconstruct 10 frames of the human lung images. The time step for each frame
is around 0.66 s ($t_{n}-t_{n-1}\approx 0.66$ s for $n=1,2,\cdots,10$). In
Fig. I we compare the reconstructions using TV, FER, first-order TV, the NWATV
and the modified NWATV regularizer (MNWATV). In this experiment the
regularization parameters for the MNWATV regularizer are set to be
$\lambda=1\times 10^{-14}$, $\rho=1.25\times 10^{-5}$, $\delta=7\times
10^{-8}$, $M=4$ and $\lambda_{b}=1\times 10^{-7}$. Again for the regularizers
of TV and the first order TV, the regularization parameters are set to be as
optimal as we can.

Due to the noise in the measured datum, in order to block the error
propagation in the iteration process, in each step we modify the
reconstruction result by formula (12). The results of the reconstructions
using each regularizer are shown in Fig. I. Especially the results of the
modified reconstructions are shown in the last row of Fig. I. Table III
illustrates the behavior of $\widetilde{RE}(n)$ and $\widetilde{PSNR}(n)$ for
the frames $t_{3}$, $t_{6}$ and $t_{10}$. In Table IV we provide the
computational time when using different regularizers.

As we can see through the use of the proposed method, we obtain a similar
reconstruction as observed in TV images and a more accurate image than FER.
Moreover, the computational cost is acceptable in clinical applications since
it produces up to 6 frames per seconds,which is much higher than the typical
respiratory rate.

Figure 10: 10 frames of human lung images using different regularizers. The
first to the third row illustrate the reconstruction results using the FER,
the TV and the first-order TV regularizers while the fourth to the fifth row
illustrate the results using the NWATV and MNWATV regularizers respectively
for $t=t_{1},\cdots,t_{10}$.

TABLE III: The behaviors of $\widetilde{RE}(n)$ and $\widetilde{PSNR}(n)$ for the frames $t_{3}$, $t_{6}$ and $t_{10}$ | FER, $\lambda=\infty$ | First-order TV | MNWATV
---|---|---|---
| $\widetilde{RE}$
$t_{3}$ | 0.0107 | 0.0039 | 0.0021
$t_{6}$ | 0.0089 | 0.0037 | 0.0022
$t_{10}$ | 0.0081 | 0.0034 | 0.0017
| $\widetilde{PSNR}$
$t_{3}$ | 19.13 | 27.83 | 33.70
$t_{6}$ | 20.63 | 28.19 | 35.91
$t_{10}$ | 21.36 | 28.71 | 35.00

TABLE IV: Comparison of the five methods in terms of computational time for human experiments. | Computational time (s)
---|---
| TV | FER | First-order TV | NWATV | MNWATV
$t_{1}$ | 2.055 | 0.020 | 0.116 | 0.109 | 0.153
$t_{2}$ | 2.131 | 0.017 | 0.104 | 0.105 | 0.140
$t_{3}$ | 2.056 | 0.019 | 0.119 | 0.108 | 0.169
$t_{4}$ | 2.071 | 0.018 | 0.109 | 0.113 | 0.163
$t_{5}$ | 2.081 | 0.021 | 0.116 | 0.112 | 0.149
$t_{6}$ | 2.088 | 0.022 | 0.107 | 0.129 | 0.158
$t_{7}$ | 2.140 | 0.020 | 0.105 | 0.122 | 0.199
$t_{8}$ | 2.136 | 0.019 | 0.110 | 0.109 | 0.176
$t_{9}$ | 2.157 | 0.02 | 0.107 | 0.108 | 0.135
$t_{10}$ | 2.161 | 0.019 | 0.112 | 0.1149 | 0.135

##  VI Discussion: Regularization Parameter Selection and Future Work

TV regularization has the advantage of preserving the inter-medium
discontinuities [34] especially for a piecewise constant images and definitely
has promising in EIT clinical applications. [35] introduced the TV
regularization in EIT, [13] introduce a primal dual interior-point-method for
the minimizing process. However, this method is time consuming and
consequently it is difficult to apply in real time applications. The proposed
nonlinear weighted anisotropic TV regularization method reduces the
computational time by using a weighted anisotropic TV regularizer. Moreover,
this method is promising in producing a better result than using only
traditional TV (isotropic TV) or anisotropic TV regularization. This is
because the weight $1/|\nabla\sigma|^{2}$ could avoid possible pseudo-edges in
isotropic TV or distortion along the coordinates axes in anisotropic TV
regularizations [23].

For regularization based reconstructions, the choices of the regularization
parameter(s) are always critical. In this paper, the choices of $\rho$ and
$\delta$ depend on the mesh size used in the process of inversion. $\rho$ has
been chosen in such a way that $\frac{1}{\rho}\mathbf{S}^{T}\mathbf{S}$ and
$\mathbf{D}^{T}\mathbf{D}$ in the similar order so that the formula (9) makes
sense. $\delta$ is chosen in the similar way as $\rho$ such that $\delta$ is
in the similar order as $|\mathbf{D}\bm{\sigma}|$. This explains why the
parameter $\rho$ and $\delta$ is chosen so small. From the formula (10), only
the ratio $\lambda/\rho$ is meaningful. Hence, the parameter $\lambda$ has to
be chosen according to the selection of $\rho$. Fig 11 (a) and (b) illustrate
the change of $RE$ and $PSNR$ with respect to $\delta$ and $\lambda/\rho$,
respectively. As we can see from this figure, for a fixed $\lambda/\rho$, $RE$
and $PSNR$ are almost invariant with respect to $\delta$. Moreover, there is
optimal choice of $\lambda/\rho$ to minimize the $RE$ and maximize $PSNR$.

Future studies should cover a strict mathematical theory on the convergence
analysis of the iteration process. For clinical applications, a potential
consideration is to add the low rank constraints of the sequential images [9]
in the reconstruction.

(a)

(b)

Figure 11: Changes of $RE$ and $PSNR$ with respect to $\frac{\lambda}{\rho}$
and $\delta$ in the 10th image of two-dimensional numerical simulation.

##  VII Conclusion

In this paper, we propose a nonlinear weighted anisotropic TV regularization
method in EIT to reduce the computational time and improve the capability of
edge preservation in comparison with TV regularization based imaging. To
validate the advantages of the proposed regularization method corroborating
with the well established FER, TV, and first-order TV regularizers, we carried
out 2D, 3D simulations and human EIT lung imaging. In the testing campaign, it
was shown that the proposed method reduces the computational time
significantly while providing the reconstruction images similar to or even
better than the traditional TV method.

In this section we provide the proofs of the formula (10).

###  -A Proof of the formula (10)

Direct calculation yields that

$\begin{split}\mathbf{z}_{n+1}&=\mathop{\arg\min}_{\mathbf{z}}\lambda\|\mathbf{p}_{n}\cdot\mathbf{z}\|_{{l_{1}}}\\\
&\qquad\qquad+\frac{\rho}{2}\|\mathbf{D}\delta\bm{\sigma}_{n+1}-\mathbf{z}\|_{l_{2}}^{2}+\mathbf{y}_{n}^{T}(\mathbf{D}\delta\bm{\sigma}_{n+1}-\mathbf{z})\\\
&=\mathop{\arg\min}_{\mathbf{z}}\|\mathbf{D}\delta\bm{\sigma}_{n+1}-\mathbf{z}\|_{l_{2}}^{2}\\\
&\qquad\qquad+\frac{2}{\rho}\mathbf{y}_{n}^{T}(\mathbf{D}\delta\bm{\sigma}_{n+1}-\mathbf{z})+2\frac{\lambda}{\rho}\|\mathbf{p}_{n}\cdot\mathbf{z}\|_{{l_{1}}}\\\
&=\mathop{\arg\min}_{\mathbf{z}}\left\|\mathbf{z}-(\mathbf{D}\delta\bm{\sigma}_{n+1}+\frac{\mathbf{y}_{n}}{\rho})\right\|_{{l_{2}}}^{2}\\\
&\qquad\qquad+\frac{2\lambda}{\rho}\|\mathbf{p}_{n}\cdot\mathbf{z}\|_{l_{1}}-\left\|\frac{\mathbf{y}_{n}}{\rho}\right\|_{{l_{2}}}^{2}\\\
&=\mathop{\arg\min}_{\mathbf{z}}\left\|\mathbf{z}-(\mathbf{D}\delta\bm{\sigma}_{n+1}+\frac{\mathbf{y}_{n}}{\rho})\right\|_{{l_{2}}}^{2}+\frac{2\lambda}{\rho}\|\mathbf{p}_{n}\cdot\mathbf{z}\|_{l_{1}},\end{split}$

where the last equality comes from the fact that
$\|\mathbf{y}_{n}/\rho\|_{{l_{2}}}^{2}$ is independent of $\mathbf{z}$. From
the definition of $l_{p}$ ($p=1,2$) norm we obtain that

$\mathbf{z}_{n+1}=\mathop{\arg\min}\sum_{k=1}^{N}F_{\rho,\lambda}(\mathbf{z}[k]),$

where $F_{\rho,\lambda}[\mathbf{z}[k]]$ is defined as follows

$\begin{split}F_{\rho,\lambda}[\mathbf{z}[k]]&=\left\|\mathbf{z}[k]-\left(\mathbf{D}\delta\bm{\sigma}_{n+1}[k]+\frac{\mathbf{y}_{n}[k]}{\rho}\right)\right\|_{{l_{2}}}^{2}\\\
&\qquad\qquad+\frac{2\lambda}{\rho}|\mathbf{p}_{n}[k]\mathbf{z}[k]|.\end{split}$

Hence

$\begin{split}\mathbf{z}_{n+1}[k]=&\mathop{\arg\min}F_{\rho,\lambda}(\mathbf{z}[k]).\end{split}$
(13)

Next we will calculate the minimizer of $F_{\rho,\lambda}(\mathbf{z}[k])$,
$\mathbf{z}_{n+1}[k]$ explicitly. Indeed, from the first-order optimal
condition we obtain

$\begin{split}0&=\left.\frac{1}{2}\frac{\partial
F_{\rho,\lambda}(\mathbf{z}_{n+1}[k])}{\partial
z}\right|_{z=\mathbf{z}_{n+1}[k]}\\\
&=\mathbf{z}_{n+1}[k]-\left(\mathbf{D}\delta\bm{\sigma}[k]+\frac{\mathbf{y}[k]}{\rho}\right)+\frac{\lambda\mathbf{p}[k]}{\rho}\mathrm{sgn}(\mathbf{z}_{n+1}[k]).\end{split}$

Hence $\mathbf{z}_{n+1}[k]$ satisfies the following identity

$\mathbf{z}_{n+1}[k]=\mathbf{D}\delta\bm{\sigma}_{n+1}[k]+\frac{\mathbf{y}_{n}[k]}{\rho}-\frac{\lambda\mathbf{p}_{n}[k]}{\rho}\mathrm{sgn}(\mathbf{z}_{n+1}[k]).$
(14)

From (14) and the fact that $\lambda\mathbf{p}[k]/\rho\geq 0$ for
$k=1,2,\cdots,N$, it is obvious that

$\mathrm{sgn}(\mathbf{z}_{n+1}[k])=\mathrm{sgn}(\mathbf{D}\delta\bm{\sigma}_{n+1}[k]+\mathbf{y}_{n}[k]/\rho).$
(15)

Noting that from the basics of calculus, the candidates of
$\mathbf{z}_{n+1}[k]$ is either the stationary point or the non-differentiable
point of $F_{\rho,\lambda}$. Hence

  * •

For the case
$|\mathbf{D}\delta\bm{\sigma}_{n+1}[k]+\mathbf{y}_{n}[k]/\rho|>\frac{\lambda\mathbf{p}_{n}[k]}{\rho}$.

From (11), (14) and (15) we obtain

$\mathbf{z}_{n+1}[k]=(I-\lambda\mathbf{p}_{n}[k]/\rho\cdot\mbox{sgn})[\mathbf{D}\delta\bm{\sigma}_{n+1}[k]+\mathbf{y}_{n}[k]/\rho]$

or $\mathbf{z}_{n+1}[k]=0$, where $I$ is the identity map.

Since

$\begin{split}&~{}~{}~{}~{}F_{\rho,\lambda}(\mathbf{D}\delta\bm{\sigma}_{n+1}[k]+\mathbf{y}_{n}[k]/\rho]\mp\lambda\mathbf{p}_{n}[k]/\rho)\\\
&=\frac{\lambda^{2}}{\rho^{2}}(\mathbf{p}_{n}[k])^{2}+\frac{2\lambda}{\rho}\mathbf{p}_{n}[k]\cdot\\\
&\qquad\qquad\qquad\qquad\left[\pm\left(\mathbf{D}\delta\bm{\sigma}_{n+1}[k]+\frac{\mathbf{y}_{n}[k]}{\rho}\right)\mp\frac{\lambda}{\rho}\mathbf{p}_{n}[k]\right]\\\
&=\left\\{\begin{split}&\frac{2\lambda}{\rho}\mathbf{p}_{n}[k]\left(\mathbf{D}\delta\bm{\sigma}_{n+1}[k]+\frac{\mathbf{y}_{n}[k]}{\rho}\right)-\frac{\lambda^{2}}{\rho^{2}}\mathbf{p}_{n}^{2}[k],\\\
&\qquad\qquad\qquad\mbox{if
}\mathbf{D}\delta\bm{\sigma}_{n+1}[k]+\frac{\mathbf{y}_{n}[k]}{\rho}>\frac{\lambda\mathbf{p}_{n}[k]}{\rho}\\\
&-\frac{2\lambda}{\rho}\mathbf{p}_{n}[k]\left(\mathbf{D}\delta\bm{\sigma}_{n+1}[k]+\frac{\mathbf{y}_{n}[k]}{\rho}\right)+3\frac{\lambda^{2}}{\rho^{2}}\mathbf{p}_{n}^{2}[k],\\\
&\qquad\qquad\qquad\mbox{if
}\mathbf{D}\delta\bm{\sigma}_{n+1}[k]+\frac{\mathbf{y}_{n}[k]}{\rho}<-\frac{\lambda\mathbf{p}_{n}[k]}{\rho}\\\
\end{split}\right.\\\
&<\left\|\mathbf{D}\delta\bm{\sigma}_{n+1}[k]+\frac{\mathbf{y}_{n}[k]}{\rho}\right\|_{l_{2}}^{2}=L(0).\end{split}$

Hence in this case, (15) holds.

  * •

For the case
$|\mathbf{D}\delta\bm{\sigma}_{n+1}[k]+\mathbf{y}_{n}[k]/\rho|\leq\frac{\lambda\mathbf{p}_{n}[k]}{\rho}$.

In this case, if $\mathbf{z}_{n+1}[k]\neq 0$ we have the following relation

$\begin{split}&\mathrm{sgn}((I-\lambda\mathbf{p}_{n}[k]/\rho\mbox{sgn})[\mathbf{D}\delta\bm{\sigma}_{n+1}[k]+\mathbf{y}_{n}[k]/\rho])\\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad+\mathrm{sgn}(\mathbf{z}_{n+1}[k])\equiv
0.\end{split}$

This contradicts with the relation (15). Hence the only minimizer of (13) is
$\mathbf{z}_{n+1}[k]=0$.

Finally, we obtain the formula (10).

## References

  * [1] D. C. Barber and B. H. Brown, “Applied potential tomography,” _J. Phys. E: Sci. Instrum._ , vol. 17, no. 9, pp. 723-733, 1984.
  * [2] J. N. Tehrani, A. McEwan, C. Jin and A. van Schaik “L1 regularization method in electrical impedance tomography by using the L1-curve (Pareto frontier curve)” _Appl. Math. Model._ vol. 36, 1095-1105, 2012.
  * [3] P. C. Hansen “The truncated SVD” as a method for regularization” _BIT_ vol. 27, pp. 534-553, 1987.
  * [4] M. Vauhkonen, D. Vadász, P. A. Karjalainen, E. Somersalo, and J. P. Kaipio, “Tikhonov regularization and prior information in electrical impedance tomography,” _IEEE Trans. Med. Imag._ , vol. 17, no. 2, pp. 285-293, 1998.
  * [5] M. K. Choi, B. Harrach, and J. K. Seo, “Regularizing a linearized EIT reconstruction method using a sensitivity-based factorization method,” _Inverse Probl. Sci. Eng._ , vol. 22, no. 7, pp. 1029-1044, 2014.
  * [6] M. Cheney, D. Isaacson and J. C. Newell, “Electrical impedance tomography,” _SIAM Rev._ , vol. 41, no. 1, pp. 85-101, 1999.
  * [7] K. Zhang, M. Li, F. Yang, S. Xu and A. Abubakar, “Three-dimensional electrical impedance tomography with multiplicative regularization,” _IEEE Trans. Biomed. Eng._ , vol. 66, no. 9, pp. 2470-2480, 2019.
  * [8] L. Zhou, B. Harrach, and J. K. Seo, “Monotonicity-based electrical impedance tomography for lung imaging,” _Inv. Probl._ , vol. 34, no. 4, p. 045005, 2018.
  * [9] Q. Wang, F. Li, J. Wang, X. Duan and X. Li “Towards a Combination of Low Rank and Sparsity in EIT Imaging” _IEEE Access_ vol. 7, pp. 156054-156064, 2019.
  * [10] J. Wang “Non-convex $\ell_{p}$ regularization for sparse reconstruction of electrical impedance tomography” _Inverse Probl. Sci. En._ vol. 29, pp. 1-22, 2020.
  * [11] Y. Shi, X. Kong, M. Wang, Y. Wu and L. Yang, “A non-convex $L_{1}$ norm penalty-based total generalized variation model for reconstruction of conductivity distribution” _IEEE Sens. J._ vol. 20, pp. 8137-8146, 2020.
  * [12] L. I. Rudin, S. Osher and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” _Phys. D_ , vol. 60, pp. 259-268, 1992.
  * [13] A. Borsic, B. M. Graham, A. Adler, and W. R. B. Lionheart, “In vivo impedance imaging with total variation regularization,” _IEEE Trans. Med. Imag._ , vol. 29, no. 1, pp. 44-54, 2010.
  * [14] J. Li, S. Yue, M. Ding, Z. Cui and H. Wang “Adaptive $L_{p}$ regularization for electrical impedance tomography” _IEEE Sens. J._ vol. 19, pp. 12297-12305, 2019.
  * [15] Y. Shi, J. Liao, M. Wang, Y. Li, F. Fu and M. Soleimani “Total fractional-order variation regularization based image reconstruction method for capacitively coupled electrical resistance tomography” _Flow. Meas. Instrum._ vol. 82, pp. 1-8, 2021.
  * [16] Y. Shi, Z. Tian, M. Wang, X. Kong, L. Li and F. Fu, “A wavelet frame constrained total variation generalized vational model for image conductivity distribution” _Inverse Probl. Imag._ in presee.
  * [17] C. Wu and X. Tai, “Augmented Lagrangian method: dual methods, and split Bregman iteration for ROF, Vectorial TV, and high order models,” _SIAM J. Imag. Sci._ , vol. 3, no. 3, pp. 300-339, 2010.
  * [18] Y. M. Jung and S. Yun, “Impedance imaging with first-order TV regularization,” _IEEE Trans. Med. Imag._ , vol. 34, no. 1, pp. 193-202, 2015.
  * [19] Z. Zhou, G. S. dos Santos, T. Dowrick, J. Avery, Z. Sun, H. Xu and D. S. Holder, “Comparison of total variation algorithms for electrical impedance tomography,” _Physiol. Meas._ , vol. 36, no. 6 pp. 1193-1209, 2015.
  * [20] K. Lee, E. J. Woo, and J. K. Seo, “A fidelity-embedded regularization method for robust electrical impedance tomography,” _IEEE Trans. Med. Imag._ , vol. 37, no. 9, pp. 1970-1977, 2018.
  * [21] J. K. Seo, K. C. Kim, A. Jargal, K. Lee, and B. Harrach, “A learning-based method for solving ill-posed nonlinear inverse problems: A simulation study of lung EIT,” _SIAM J. Imag. Sci._ , vol. 12, No. 3, pp. 1275-1295, 2019.
  * [22] S. Boyd, N. Parikh, E. Chu, B. Peleato and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” _Found. and Trends in Mach. Learn._ , vol. 3, no. 1, pp. 1-122, 2011.
  * [23] G. González, V. Kolehmainen, A. Seppänen, “Isotropic and anisotropic total variation regularization in electrical impedance tomography” _Comput. Math. with Appl._ vol. 74, pp. 564-576, 2017.
  * [24] B. Gong, B. Schullcke, S. Krueger-Ziolek, F. Zhang, U. Mueller-Lisse, K. Moeller, “Higher order total variation regularization for EIT reconstruction” _Med. Biol. Eng. Comput._ vol. 56, pp. 1367-1378, 2018.
  * [25] B. Harrach and J. K. Seo, “Exact shape-reconstruction by one-step linearization in electrical impedance tomography,” _SIAM J. Math. Anal._ , vol. 42, no. 4, pp. 1505-1518, 2010.
  * [26] I. Daubechies, M. Defrise, and C. De Mol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” _Commun. Pure Appl. Math._ , vol. 57, no. 11, pp. 1413-1457, 2004
  * [27] A. Adler, R. Guardo, and Y. Berthiaume, “Impedance imaging of lung ventilation: Do we need to account for chest expansion?” _IEEE Trans. Biomed. Eng._ , vol. 43, no. 4, pp. 414-420, Apr. 1996.
  * [28] A. Adler and W. R. B. Lionheart, “Uses and abuses of EIDORS: an extensible software base for EIT,” _Physiol. Meas._ , vol.27, no. 5, pp.S25-S42, 2006.
  * [29] M. Hanke, B. Harrach, and N. Hyvönen, “Justification of point electrode models in electrical impedance tomography,” _Math. Models Methods Appl. Sci._ , vol. 21, no. 6, pp. 1395-1413, 2011.
  * [30] D. Liu, V. Kolehmainen, S. Siltanen, A-M. Laukkanen and A. Seppänen, “Nonlinear difference imaging approach to three-dimensional electrical impedance tomography in the presence of geometric modeling errors,” _IEEE Transactions on Biomedical Engineering_ , vol. 63, no. 9, pp. 1956-1965, 2016.
  * [31] Electrode company. _Xi’an friendship medical electronics company_. Accessed: Aug. 24, 2021. [Online]. Available: http://www.xafdec.com
  * [32] Sense 3D Scanner. _3D systems_. Accessed: Aug. 24, 2021. [Online]. Available: https://support.3dsystems.com/s/article/Sense-Scanner
  * [33] Sciospec company. https://www.sciospec.com/.
  * [34] A. Javaherian, M. Soleimani, K. Moeller, A. Movafeghi and R. Faghihi, “An accelerated version of alternating direction method of multipliers for TV minimization in EIT,” _Appl. Math. Model._ , vol, 40, pp. 8985-9000, 2016.
  * [35] D. C. Dobson and F. Santosa, “An image-enhancement technique for electrical impedance tomography,” _Inv. Probl._ , vol. 10, no. 2, pp. 317-334, 1994.