Spatially-Adaptive Variational Reconstructions for Linear Inverse Electrical Impedance Tomography

The inverse electrical impedance tomography (EIT) problem involves collecting electrical measurements on the smooth boundary of a region to determine the spatially varying electrical conductivity distribution within the bounded region. Effective applications of EIT technology emerged in different areas of engineering, technology, and applied sciences. However, the mathematical formulation of EIT is well known to suffer from a high degree of nonlinearity and severe ill-posedness. Therefore, regularization is required to produce reasonable electrical impedance images. Using difference imaging, we propose a spatially-variant variational method which couples sparsity regularization and smoothness regularization for improved EIT linear reconstructions. The EIT variational model can benefit from structural prior information in the form of an edge detection map coming either from an auxiliary image of the same object being reconstructed or automatically detected. We propose an efficient algorithm for minimizing the (non-convex) function based on the alternating direction method of multipliers. Experiments are presented which strongly indicate that using non-convex versus convex variational EIT models holds the potential for more accurate reconstructions.


Introduction
Electrical impedance tomography is an imaging technique that aims to reconstruct the inner conductivity distribution of a medium starting from a set of measured voltages registered by a series of electrodes that are positioned on the surface of the medium. EIT is therefore a nondestructive testing technique, meaning that it allows to analyse the property of a material or structure without causing damage. It can be considered a tomographic modality due to the fact that it generates images of the internal features of a body. However, if compared to other tomography techniques, EIT provides lower spatial resolution outputs, but data can be B S. Morigi acquired relatively fast (EIT temporal resolution is estimated in the order of millisecond) and the apparatus is more manageable.
The first applications of EIT techniques date back to 1930 on geology [35], and to the mid 80s with the first clinical use by Brown and Barber. Since then EIT technology has had a huge development in various application fields. Applications of EIT in biomedical imaging range from clinical imaging for organ monitoring to cell monitoring in tissue engineering [27,28]. The promising advantages of this imaging procedure over X-Ray, CT and MRI are the fact that the device can be brought to the patient and no exposition to radiation is required for data collection. Recently, feasibility studies on biomaterial-engineering applications of EIT have focused on non-invasive approaches to monitoring and tissue reconstruction [9].
In addition to the broad spectrum of possible application fields, EIT has aroused considerable interest also from a mathematical point of view, both in the forward problem and in the inverse problem, since the pioneering work of Calderón [4]. The forward problem corresponds to the prediction of the boundary voltage distribution when the conductivity distribution and the injected current are known. The governing model for the forward problem, known as Complete Electrode Model, is based on an elliptic partial differential equation subject to a set of constraints and Neumann boundary conditions that account for the fitting with real data, discreteness of the electrodes and the extra parameters defined by the measuring devices.
The inverse conductivity problem, as formulated by Calderón [4], is to infer the distributed conductivity based on boundary potential measurements. If the forward problem is characterized by a nonlinear mapping, called Forward Operator, from the conductivity distribution to the measured voltages, the inverse EIT consists of inverting the Forward Operator. However, due to the ill-posedness nature of the problem, some sort of regularization is required to stabilize the reconstruction. Regularization methods make use of some a priori structural information about the unknown conductivity, which proves crucial in the case of limited and noisy data, as is often in practical EIT situations.
For instance, the method commonly known as total variation regularization, aims to enforce sparsity in the magnitude conductivity gradient distribution, and thus it is suitable for reconstructing piecewise constant conductivities, while the well-known Tikhonov regularization is widely used to promote smoothness in the solution. However, in many practical situations, a conductivity distribution can be naturally characterized by homogeneous and constant regions as well as inhomogeneous and smooth parts. This space-variant aspect of the distribution can undermine the choice of the proper regularizer to be applied to the entire distribution.
A spatially-adaptive regularization allows for a locally different regularization effect, so that local, structural properties of the sought-for conductivity can be potentially addressed. The usefulness of this great flexibility is however conditioned to the existence of effective procedures for the automatic estimation of a space-variant map of regularization parameters based on a priori knowledge on the solution. However, unlike in the image restoration scenario, where the input data is a degraded version of the sought-for solution and can therefore provide some information about it, in EIT reconstruction problems, the data on the boundary does not give any useful insights on the solution.
In the present paper, we propose a new spatially-adaptive variational model for the resolution of the linear EIT inverse problem in time difference imaging reconstructions. The functional to be minimized includes the linearized version of the forward operator and two nonlinear regularization terms, driven by a space-variant regularization map that incorporates structural prior information and sparsity. In particular, we consider a smooth convex generalized Tikhonov-type regularizer and a nonsmooth non-convex TV-like regularizer to enforce sparsity in the conductivity reconstruction. The effect of these two regularizers is controlled by a space-variant function which locally defines whether to favour smoothing or sparsification. This variational model is non-convex and nonsmooth so it demands for challenging numerical optimization solutions. We propose an ADMM-based iterative algorithm for the minimization of the proposed variational model. In the field of image and signal processing, the ADMM has been one of the most powerful and successful methods for solving various convex or non-convex optimization problems, as it resorts to the variable splitting which reduces to the alternating resolution of simpler optimization problems. Numerical examples of EIT conductivity reconstructions show that the proposed approach is particularly effective and well suited for a wide range of spatially variant conductivity distributions.
The paper is organized as follows. In Sect. 1.1 we briefly review the mathematical approaches to the inverse EIT problem. In Sect. 2 the forward model is introduced and Sect. 3 presents two widely known variational approaches for the regularization of the reconstruction EIT problem. The proposed variational space-variant model is described in Sect. 4 together with the procedure proposed for the automatic estimation of the space-variant regularization map, Sect. 4.1. The ADMM-based minimization algorithm is illustrated in Sect. 5 and numerical results are reported and commented in Sect. 6. Conclusions are drawn in Sect. 7.

Related Works
Impedance measurements are characterized by high temporal resolution and low-cost instrumentation. This fueled intense research in EIT image reconstruction in diverse fields over the last few decades.
Early image reconstruction attempts by Santosa and Vogelius [33] relied on the backprojection approach. This is derived from the Radon transform and is commonly applied for the analysis of computed tomography data. This approach was characterized by important limitations, such as low resolution and severe artifacts that were partly addressed by the Landweber method [37], obtained by modifying the generalized inverse matrix method. However, the strictly bi-dimensional approach and the lack of prior information on the domain to be reconstructed makes the Landweber method useful only for very simplified experimental setups.
Variational approaches are currently the most successful and widely used EIT reconstruction methods. A notable example is Tikhonov reconstruction, which is a regularized least squares approach that contains a penalty term to favour smooth conductivity reconstructions [40]. Generalized Tikhonov approaches allow for a more general smoothing penalty term which includes approximations of differential operators [23] or structural priors [36]. Another common regularization strategy relies on the 1 norm to promote the reconstruction of piece-wise constant conductivity regions since it maintains the discontinuities in the reconstructed profiles. This regularization method is generally referred to as Total Variation [3]. This feature makes it particularly suited for medical and industrial applications (e.g. definition of inter-organ and inter-phase boundaries).
The possibility of adding further prior information about the domain to be reconstructed has been a crucial advancement for EIT reconstruction algorithms, that has allowed for reallife problems to be tackled more effectively. Among these methods, sparse Bayesian learning has recently been used to improve the spatial resolution of the reconstructed conductivity distributions [24] and reduce systematic artifacts in 3D volumes [25]. Alternatively, subdomain methods and similar approaches rely on prior information to limit the reconstruction to specific regions of the domain [18]. While these methods are generally application-specific, they are associated with increased accuracy and a reduced reconstruction time [17,26,32]. Several branches of artificial intelligence, such as neural networks, swarm intelligence, and evolutionary methods have been applied to the EIT image reconstruction problem. In particular, particle swarm optimization has been used in the past to improve the reconstruction quality while reducing the number of iterations [30] and to optimize the current injection pattern [38]. More recently, neural networks have been extensively applied to increase image quality and reduce reconstruction time [10,11]. The D-Bar method has lately gained traction for EIT imaging, thanks to a quick non-iterative approach relying on the Fourier and inverse transform [13] which has also been combined with deep learning techniques to improve reconstruction quality [12].
In this work, we focus on variational methods for EIT image reconstruction, due to the proven wide applicability of these approaches in real-life scenarios [19].

Formulation of the Forward EIT Problem
Let ⊂ R d , d = 2, 3 be a bounded simply connected C ∞ domain, which represents the measurement domain, and σ : → R a strictly bounded measurable function (σ (x) ≥ c > 0) which denotes the electric conductivity. In the absence of interior current sources and in the quasi-electrostatic case, Maxwell's system describing electromagnetic fields inside the domain reduces to the following second-order elliptic partial differential equation which models the electric potential u ∈ H 1 ( ) inside the body where U (x) is the electrical potential distribution on the boundary. Let σ : H ∂n , and is defined as with n the outer normal at x ∈ ∂ . Hence problem (1) can be expressed with Neumann boundary conditions The boundary value problem (1)-(4), for a known function σ (x) in and data I (x), given for all x ∈ ∂ , is referred to as Continuum, forward mathematical model for EIT, which implies complete knowledge of all boundary data.
More recently, in order to provide a more physically realistic mathematical model, some electrode models have been introduced which consider a priori knowledge on the experimental setup: the Gap Model, which allows for regular gaps in the boundary data, the Shunt Electrode Model (SEM), that assumes a finite number of electrodes with associated constant potential, and the Complete Electrode Model (CEM), introduced in [7] and [34], which adds a complex impedance for each electrode in order to model the metal electrode, conductive gel and chemical interaction at the skin-electrode interface.
Let L be the number of electrodes, which is known. Following the CEM model, the boundary of can be modelled as a family of subsets of ∂ : • Boundary with electrodes • Boundary without electrodes˜ where |E l | is the measure (length for d = 2 or area for d = 3) of electrode E l . This leads to replace the Neumann boundary conditions (4) by two weaker conditions where I l is the known current sent to the lth electrode E l , and the current density is zero on the boundaries between the electrodes (on˜ ). The CEM model takes into account the extra resistance (due to an electrochemical effect) between the electrode and the tank, id est the formation of a thin, highly resistive layer. A parameter z l , called effective contact impedance, defines the contact impedance between the object and the electrode E l . Boundary conditions (2) on E l , l = 1, . . . , L, are therefore replaced in the CEM model by the following Robin boundary conditions where U l denotes the unknown voltage to be measured on the lth electrode and predicted by the model. Summarizing, the accurate CEM model for the forward EIT problem can be formulated as follows: where x ∈ R d is the position vector. The forward EIT problem consists in determining the potentials u in , and U l on the electrodes when the applied currents I l , conductivity σ and contact impedances z l are known; its solution amounts to solving the boundary problem (8). In [34] the authors proved that the CEM model for the EIT forward problem has a unique solution and the solution depends continuously on the input current I . Specifically, the existence of the solution is obtained by assuming the charge conservation law L l=1 I l = 0, and the uniqueness is obtained by choosing a ground level for the voltages, for example by L l=1 U l = 0. Numerically, we will consider approximate solutions of (8) obtained by applying the finite element method [22].
The forward problem can be restricted to a relation between the inner conductivity σ and the boundary voltages U (x) and can be modelled by a mapping which is referred to as Forward OperatorF where S = {σ ∈ L ∞ ( ) | σ ∇u = 0}.
The discrete forward operator, denoted by F(σ, I ), obtained with the aid of FEM discretization of (8), relates, for given σ , the applied electrode current data I with the computed electrode voltage data U , namely U = F(σ, I ), and depends linearly on I but nonlinearly on σ . For a fixed current vector I , we can view F(σ ; I ) as a function of σ only, which will be denoted simply by F(σ ).
An EIT experimental setup follows a given stimulation pattern, which establishes as the n M measurements are collected, each obtained by injecting current from an electrodes pair and then measuring the corresponding induced voltage U m on another pair of electrodes. Let the domain be discretized into n T subdomains {τ j } n T j=1 and let σ be considered constant on each of them. The forward operator F can be seen as a nonlinear vector map from R n T → R n M .
As the measured data is unavoidably noisy, we can adopt an additive Gaussian noise model for the measurement error, assuming the following noisy forward observation model where U m ∈ R n M represents the vector of all the measured electrode potentials whose dimension n M depends on the choice of a measurement protocol, andn ∈ R n M is a zeromean Gaussian distributed measurement noise vector.

The Linearized Forward Operator and the Sensitivity Matrix
Calderón in [4] showed that the forward mappingF, defined asF(σ ) = σ (U ) with σ being the DtN map (3), is Fréchet differentiable with respect to σ , and gave an explicit formula for the derivativeF which leads to the definition of the best linear approximation ofF at some fixed value of σ . In the following the formula for the derivative and the linear approximation are derived using the weak formulation of the forward problem. An equivalent integral formulation of (1) is obtained by considering the standard variational formulation of the Dirichlet problem (1), using the test function v ∈ H 1 ( ) such that, if u, v are solutions of (1)-(2) with Dirichlet data U ,V respectively, then where for u = v, Q σ (U , U ) represents the total power dissipated in when the distribution voltage on the boundary is U , and the last term in (12) is the power input across the boundary. In particular, Calderón in [4] showed that, considering a change in conductivity from σ 0 to σ 0 + δ, then IfF (σ ) is the Fréchet derivative of the forward mapF(σ ) at conductivity σ then Note that the discrete forward operator F is a nonlinear vector field. SinceF is Fréchet differentiable, F is a matrix, called the Jacobian of F and denoted by J . In order to derive a formula for J we follow [6], using the discretization of the domain introduced in Sect. 2, and assuming an experimental protocol of n M measures. Let δ in (15) be the characteristic function of the jth subdomain, then each element of the symmetric matrix J ∈ R n M ×n T is defined as The row index i corresponds to the ith measurement, associated with the dth driving potential u d and mth measurement potential u m , while the column index j corresponds to the subdomain τ j . Matrix (16) is called the sensitivity matrix since it measure the sensitivity of a boundary voltage to a change in conductivity. Finally, (14) leads to the linear approximation of the discrete operator F where J defined in (16) is calculated at the initial conductivity estimate σ 0 .

Inverse EIT Problem
In the EIT forward problem, the potential distribution u inside satisfies an elliptic equation with Cauchy data, that is with data given on a part of the boundary defined in (6); this is itself an ill-posed problem, according to the Hadamard criteria.
The EIT inverse problem aims to reconstruct the conductivity σ inside the body starting from a finite set of measured voltages U m on the boundary ∂ by inverting the nonlinear Forward Operator F(σ ). However, the inverse map F −1 is typically discontinuous (it is unbounded), the measurements are usually noisy, and thus the EIT inverse is a severely illposed reconstruction problem, as well described in [1]. This practically means that small changes in boundary measurements can correspond to large changes in the internal conductivity distribution. Consequently, regularization techniques have been adopted in order to stabilize the inversion and thus ensuring convergence of reconstruction algorithms. At this aim, we need some additional information about the conductivity distribution to restrict the conductivity to a proper space.
Therefore, assuming the degradation model (11), the natural way to reconstruct the conductivity data is to consider the following non-linear regularized least squares problem which includes a regularization term R(σ ) that constrains the solution to satisfy a priori information about the conductivity distribution, and a regularization parameter λ > 0 that controls the trade-off amidst data fitting and bounding the derivatives of σ . Iterative numerical optimization methods to solve the nonlinear least squares problem (EITNL) are necessarily time-consuming. A conventional approach for a more efficient solution of the EIT inverse problem relies on the linearized model (17) of the nonlinear forward operator F, which is accurate for a sufficiently small change about an initial conductivity distribution σ 0 . The reconstruction of a small conductivity change is thus obtained by solving the following linear regularized least squares problem where δσ = σ − σ 0 and δU m = U m − F(σ 0 ). It is important to highlight the fact that due to the linearization of the problem, the fidelity term in (EITL) leads to an underdetermined linear system of equations to be solved, which further demands for a regularization. The functional R(·) is often assumed to be either linear or nonlinear. A standard choice for the prior R(·) is a linear term within a L 2 norm framework Among the several choices for the linear operator L in literature, the one suggested by the NOSER algorithm in [6] takes L as a positive diagonal matrix given by the diagonal elements of J T J . Another typical simple choice which promotes smoothness in the solution is the well-known Tikhonov regularization, where L is the discrete Laplacian filter to penalize for non-smooth regions in the conductivity. The regularizer R(·) can also be chosen to be nonlinear. The most common choice in this case is the Total Variation (TV) functional, introduced in [31] and proposed in EIT inverse problem by Borsic [3]: Total Variation is a L 1 norm regularization penalty, hence it allows to preserve discontinuities in the reconstructed conductivity. This structural characteristic, as described in [3], is common to several fields where EIT has been applied: in medicine, where the conductivity discontinuities are defined by internal organ boundaries, as each organ has different electric features, or in industrial process tomography where they are defined by the different phase interfaces of a multiphasic fluid. The gain in the accuracy of the reconstruction is however obtained at the expense of the loss of differentiability for the objective function, thus leading to consider non-smooth optimization methods. The different notation σ and δσ aims to underline the fact that the former represents an absolute conductivity distribution while the latter represents a relative change of conductivity distribution at a posterior state with respect to an initial state with conductivity distribution σ 0 . Linearized models are generally involved in the reconstruction of such conductivity change δσ . The idea of computing the difference in conductivity between two states is of interest to many applications, and it is known as difference imaging. However, the linear approximation cannot reconstruct large contrasts or complex geometries so the process must be applied iteratively. At each iteration, the solution of a regularized linear problem (EITL) consists of a small conductivity change δσ n = σ n+1 − σ n , and a new computation for the Jacobian is performed. The difference imaging is to some extent tolerant to approximation errors since absolute errors in the measured voltages are partially canceled when the difference of the measurements is computed. Due to this property the difference imaging has been favoured over absolute imaging, and is considered in the present work. From this point onwards for the simplicity of notation σ will be regarded as conductivity difference and U m as voltage difference.

A Variational Space-Variant Model for the Inverse EIT Problem
The proposed variational model aims at reconstructing the difference conductivity, from now on named σ , on the planar domain , and consists of the following minimization problem: where λ > 0 is the regularization parameter, ∇σ is the conductivity gradient defined on , and η : ⊂ R 2 → [0, 1] is a space variant function that works as a trade-off between the smooth convex quadratic regularization term, which penalizes large gradients of the conductivity, and hence favours smooth reconstructions, and the non-convex nonsmooth regularization term, which promotes sparsity in the conductivity distribution. The first term of J in (18), so-called fidelity term, represents the consistency to the measurements given on the boundary of and follows the linearization (17).
The regularization is driven by the a priori information on the sought for conductivity distribution, here represented by η(x). In many applications of EIT it is reasonable to assume that the conductivity is piecewise constant with sharp boundaries between the homogeneous regions. The sparsity-promoting penalty term of J in (18) relies on the parameterized, nonconvex function φ( · ; a) : [ 0, +∞) → R, with parameter a ≥ 0, which controls the degree of non-convexity and will be referred to as the concavity parameter. As function φ, a rescaled and reparameterized version of the minmax concave penalty function [41] was chosen, yielding the following piecewise quadratic function: The attractiveness of such non-convex penalty resides in its ability to promote sparsity more strongly than the 1 -norm TV penalty [20].

Spatially Varying Á Function
A space variant weighting function η(x) has been introduced in (18) to get a good balance between the two regularization terms. The goal of η(x) is to distinguish the local regions where to prioritize the action of the smoothing term and the ones where to promote sparsity instead, so that coexistent smoothly inhomogeneous and piecewise constant conductivity distributions can be accurately reconstructed. At this aim, a certain prior knowledge about the solution is then required. However, unlike what happens in image restoration problems where the input data is a degradated version of the sought for solution and can therefore provide some information about it, in EIT reconstruction problems, solution and data are defined on different spaces. Two scenarios can thus be identified: (a) the configuration inside the domain (tank) is known or at least visible, which allows to preliminary construct an ad hoc discretization of the domain on which to assign appropriate η values in the different regions; (b) the setup is blind, meaning no information is provided about the content of the domain (tank).
The procedure followed in case (b) in our experiments has been to consider the solution σ * 0 obtained by a preliminary step of the proposed method with η(x) ≡ 0, ∀x ∈ . This generates a piece-wise constant reconstruction σ * 0 where sharp variations are detected. Inspired by the edge-map functions used in image processing, the η-map is then obtained by applying the bounded non-negative continuous and monotonically descending function where t = ∇σ * 0 2 , and κ is a positive constant.

Space and Operator Discretization
Let the planar computational domain ⊂ R 2 being discretized by a 2D triangulated domain the set of n T triangles and E = {e j } n E j=1 the set of n E edges. We approximate the conductivity distribution σ : T → R by a piecewise-constant function over the elements of M, thus, the conductivity function attains a constant value σ j , j = 1, . . . , n T over triangle τ j ∈ T .
The gradient operator applied to a piecewise-constant function σ vanishes to zero everywhere but the mesh edges, thus ∇σ : E → R. Therefore, adapting a forward finite difference scheme, the discrete gradient operator is represented by a matrix D ∈ R n E ×n T having two nonzero elements per row defined as where the nonzero entries correspond to position of adjacent triangles τ j , τ k sharing the edge e i . Applying homogeneous Neumann boundary conditions, ∇σ | ∂ h = 0, the rows of D corresponding to the boundary edges can be eliminated from D, thus, n E will refer to the number of inner mesh edges. The Jacobian operator J in (16) is discretized on h following the implementation in the EIDORS software [29].
The weak formulation of the fidelity term in (18) vanishes to zero everywhere but on the electrodes through which the boundary voltage U m is measured, that is on ∂ h . Hence the fidelity term simplifies as Considering any injection-measurement protocol consisting of n M measurements U m ∈ R n M , each characterized by a pair of measuring electrodes m i is the portion of the boundary subset of , then we get where the simplification in the last row assumed that all the electrodes in the model are of the same length |E|/2.
The two regularizers in (18) are expressed in terms of the gradient ∇σ , which is nonvanishing only on the boundaries of each triangle of M, thus the weak formulation can be consequently rewritten in the element-wise form as a sum over all the inner mesh edges. In particular, the first regularization term in Eq. (18) can be rewritten as and, analogously, the second regularization term in Eq. (18) takes the discrete form where we denoted by η j = η(x)| e j , j = 1, . . . , n E the sampling of the η(x) function restricted at each edge, and we considered that, since the integrands are constants, the integrals over the edges reduced to the length of the edges e j , denoted as l j .

Some Issues on the Optimization Problem
We conclude by highlighting some properties of functional J in (26).
is lower semicontinuous and coercive. In particular, the mapping h attains its infimum inf h ∈ [−∞, +∞] at some point in R n .

Remark 1
Let J ∈ R n M ×n T , n M < n T , be the matrix defined in (16), with rank equal to n M , and D ∈ R n E ×n T , n E > n T , be defined in (21) with rank equal to n T − 1 and nullD = span{ (1, 1, . . . , 1) T }. It is reasonable to assume that property i) in Proposition 1 holds for the matrices J and D, which means that a constant conductivity changeσ does not belong to the null space of J . Assumeσ is in the null space of J , so that Jσ = 0, then according to the linearized model (17), F(σ ) − F(σ 0 ) =Ū m = 0, which would lead to the "absurd" conclusion that the corresponding boundary voltage does not change with respect to a constant change from an initial conductivity distribution σ 0 =σ . This suggests that any linear combination of the generating vectors of the null space of J cannot represent the constant vector. This has been experimentally verified in the conducted experiments.
On the other hand, property i) in Proposition 1 for D and J matrices is required when the EIT inverse problem is solved by generalized Tikhonov regularization approaches which have been commonly used in electrical impedance tomography in the past [15,36]. Proposition 2 Let φ( · ; a) : R + → R be a non-convex penalty function as the one defined in (19) and let the parameters (λ, η, a) satisfy conditions λ, a > 0, 0 ≤ η ≤ 1. Then, the functional J ( · ; λ, η, a) in (26) is proper, continuous, bounded from below by zero, and coercive, hence J ( · ; λ, η, a)

admits global minimizers.
Proof From the definition of the MC penalty in (19) and the form of J in (26), it follows easily that J is proper, continuous, and bounded from below by zero.
Due to the non-convex MC penalty function φ, the functional J is non-smooth and could be convex or non-convex depending on the parameter a. Since matrix J ∈ R n M ×n T , with n M < n T , cannot be full column rank, and the penalty is additive separable, we cannot resort to the convex non-convex strategy to derive conditions on the parameter a so that the associated functional J in (26) is convex, see [20].
The non-convex MC penalty, the last in (26), is continuous and bounded from above by (1 − η)l/2 and from below by zero, hence does not affect coerciveness of J .
Let f 1,2 := 1 2 ·−c 2 2 where c is a constant. They satisfy the assumption ii) in Proposition 1, with J and D satisfying assumption i) in Proposition 1 according to Remark 1, then the sum of the first two terms of J in (26) can be written as f 1 (J v) + f 2 (Dv), and it is lower semicontinuous and coercive. This implies that the total function J is also coercive and, thus, admits minimizers over R n T , equivalently, the variational model (26) admits solutions.
We finally notice that the special case of constant η = 0, requires that the penalty function φ is coercive, which is not the case of the MC penalty (19). However, many other non-convex penalty functions can be equivalently used, for example the φ log (t; a) := log(1 + at)/a, see [21].

ADMM-Based Numerical Optimization
In order to numerically solve the non-convex non-smooth minimization problem (26) over the triangulated domain h , an efficient ADMM-based iterative minimization algorithm is proposed, where we provide an additional condition to guarantee the uniqueness of the solution for one of the two ADMM sub-problems.
Towards this aim, we introduce the auxiliary variable t ∈ R n E , and we resort to the variable splitting technique, such that problem (26) can be reformulated in the following equivalent discrete form: subject to : t = Dσ .
To solve the constrained optimization problem (29)-(30) we define the augmented Lagrangian functional where β > 0 is a scalar penalty parameter and ρ ∈ R n E is the vector of Lagrange multipliers associated with the linear constraint t = Dσ in (30). The following saddle-point problem is then considered: Given vectors σ (k) and ρ (k) computed at the kth iteration (or initialized if k = 0), the (k + 1)th iteration of the ADMM-based iterative scheme applied to the solution of the saddle-point problem (31)-(32) is split into the following three sub-problems: We remark that the Jacobian J in (31) is computed with respect to an initial conductivity distribution σ 0 and is not updated throughout the iterative scheme.
In the following sections the methods to tackle the two minimization sub-problems (33) and (34) for the primal variables will be described in detail.

Solving the Sub-problem for t
Recalling the separability property of φ(·; a), defined in (19), the minimization sub-problem for t in (33) can be rewritten in the following element-wise (edge-per-edge) form: where constant terms have been omitted. Therefore, solving (36) is equivalent to solving n E one-dimensional problems for t j , j = 1, . . . , After some algebraic manipulation, we can rewrite (37) in the following form where Based on results in [16,Proposition 7.1], we can ensure strong convexity of the subproblems for the primal variables t j , j = 1, . . . , n E , if and only if the following condition holds: In case that (40) is satisfied, following [ [16], Proposition7.1], the unique solutions of the strongly convex problems in (38) can be obtained in closed-form by applying the softthresholding operator:

Solving the Sub-problem for
The minimization sub-problem for σ in (34) can be rewritten as follows: The first-order optimality conditions of the quadratic minimization problem (42) lead to the following linear system: The solution of (43) is unique when the coefficient matrix of the equivalent linear least square problem has full rank, that is which holds since, as stated in Remark 1, the null spaces of J and D intersect trivially: The n T × n T coefficient matrix of the linear system (43) is symmetric positive definite and remains constant along the iterations, thus the solution of (43) can be efficiently obtained by conjugate gradients method or by computing the Cholesky factorization once and for all, depending on the size of the problem.
The proposed ADMM-based algorithm, applied to solve the inverse EIT problem, relies on the key role of the space-adaptive η(x) function. In the following section we validate the algorithm, named Spatial Adaptive Electrical Impedance Tomography η (SAEiT η ), on synthetic data.
The special case η ≡ 1 (SAEiT 1 ) reduces the model (18) to the well-known, convex L 2 -Tikhonov model. In this case the first-optimality conditions for (18) involve the solution of the following system of linear equations for σ where L e = diag(l 1 , . . . , l n E ) is the diagonal matrix of the edge lengths. Summarizing, we can apply the proposed algorithm for the solution of the minimization problem (18) either using the space-variant η-map, SAEiT η , which locally balances the two regularization terms, or using a constant value for η over the entire h domain. In particular, in the following section we will denote by SAEiT c , c ∈ [0, 1] the choice of η ≡ c. Special cases are represented by c = 1, denoted as SAEiT 1 , which corresponds to the L 2 -Tikhonov model, and c = 0, denoted as SAEiT 0 , which corresponds to the non-convex TV-like-L 2 model.

Remark 2
The minimization sub-problems are all strictly convex and admit a unique solution under proper conditions [16]. However, this is not sufficient to guarantee the convergence of the overall ADMM algorithm. It is in fact well-known that, for non-convex problems, the convergence of the ADMM has not been well assessed, despite the fact it has been widely used in applications. We will further investigate this topic in a future work.

Numerical Experiments
In this section we evaluate the proposed SAEiT η method on a set of synthetic 2D experiments. All examples simulate a circular tank slice of unitary radius with a boundary ring of 16 electrodes; drive current value is set to 0.1 mA and conductivity of the background liquid is set to 1 ( * m) −1 . Measurements are simulated through opposite injection-adjacent measurement protocol via EIDORS software using a generic forward mesh of n T = 39488 triangles and n V = 19937 vertices. A coarse backward mesh was employed for inverse EIT solutions, without overlapping elements with the ones of the fine forward mesh so that to avoid the so called inverse crime in EIT [39]. The performance is assessed both qualitatively and quantitatively. The quantitative analysis is performed via the following slice metric: which measures how well the original conductivity distribution is reconstructed in case a ground truth (GT) conductivity distribution σ GT is known. In order to allow for a mesh independent comparison, the conductivity distributions are evaluated on an image structure of dimension 576 × 576 pixels. In all the experiments, the regularization parameter λ in (18) as well as the parameters of all the compared algorithms have been hand-tuned so as to achieve the lowest possible s norm error. In all the examples but Example 5, the setup is considered blind, that is no a priori information about the conductivity distribution is given. In this case, defined as case (b) in Sect. 4.1, the construction of the η-map consists of the following two steps: (i) SAEiT 0 : the spatially adaptive method is applied with η ≡ 0 to generate a preliminary 'observed image' solution σ * 0 ; (ii) η-map is generated via the edge-detection thresholding function (20) applied on ∇σ * 0 . The parameter κ in (20) for these examples is set as κ = 1.5 max |∇σ * 0 |.

Example 1: Benefits of the Á-map in SAEiT Á
In this example we illustrate the full potential of our method with the space-variant map η(x) ∈ [0, 1) in the reconstruction of several conductivity distributions, illustrated in Fig. 1, synthetically built to validate the proposed method. Moreover, we compare its performance with some state-of-the-art regularization algorithms for inverse EIT, implemented via EIDORS, a specific software for EIT problems. Specifically, for the class of nonlinear EIT methods (EIT-NL), we considered the Gauss-Newtons algorithm with linear regularization terms (IGNL) and non-linear regularization terms (Iterative Reweighted Least Squares, IRLS, [14]); while for the methods belonging to the class of linearized approaches (EIT-L), we compared with the One-step Gauss-Newton with Tikhonov priors, the (OGNT), the Primal-Dual Interior Point algorithm for Total Variation regularization (TVPIM, [5]), and the Newton One-Step Algorithm (NOSER [6]), which performs only the first step of the Newtons method for the solution of the (EIT-NL) problem. In Fig. 1(first row) the Empty Measurement Chamber setup is reported (EMPTY), illustrated on a mesh of n V = 750 vertices and n T = 1402 triangles. This setup, sampled on the generic forward mesh, has been used as a baseline for the time difference reconstruction. In Fig. 1 different kinds of inclusions inside the scaffold are illustrated and considered in order to test the different performances of the methods: piece-wise constant (PC-test), smooth (SMtest), and mixed piece-wise constant and smooth (PCSM-test). The colormap on the right represents the intensity of the conductivity which is imposed over the background mesh.
The η maps are generated by the blind procedure described in the beginning of Sect. 6. By the way of illustration, in Fig. 2 we illustrate the two-steps procedure for the PCSMtest setup. The preliminary conductivity distribution σ * 0 , shown in Fig. 2b, is obtained by applying SAEiT 0 to the measurements U m = F(σ ), where σ is the PCSM-test distribution illustrated in Fig. 2a. Finally, by the edge-detection procedure on ∇σ * 0 , the η-map is detected and illustrated in Fig. 2c.
The obtained space-variant η-map is then used in the proposed SAEiT η algorithm; the reconstructed piece-wise smooth conductivity distribution σ * is shown in Fig. 3 (second row) and compared with the conductivity reconstructed by TVPIM algorithm, shown in the first row of Fig. 3. Our spatially adaptive proposal, as expected, provides more accurate reconstructions for the piecewise-constant region as well as for the smooth area, which is free from staircase effects. In Fig. 3 (third row) the cross-sections along the reconstructed conductivity distributions are shown, which confirm the improvements.
From the several numerical tests reported in Table 1 in terms of error s values, the TVPIM algorithm turned out to perform better than the other competitors, hence in the remainder of this numerical section the comparisons will be limited to this method.

Example 2: SAEiT 0 Versus TV-L 2
In order to validate the performance due to the non-convex φ(·; a) penalty in our model, we set η ≡ 0, which reduces the variational model (18) to the non-convex TV-like-L 2 model (SAEiT 0 ), and compare its performance with a TV-L 2 model solved by the primal-dual internal-point method (TVPIM) proposed in [5].
The results obtained by applying TVPIM and SAEiT 0 methods to the reconstruction of the PC-test setup are reported in Fig. 4, first and second row, respectively. The qualitative benefits obtained by the sparsity-inducing regularizer in SAEiT 0 w.r.t. the TV regularizer in TVPIM range from reducing contrast loss to preservation of corners and finer details. The  slice error metric s , indicated in brackets, clearly confirms these advantages and suggest that SAEiT 0 holds the potential for promoting sparsity of the vector Dσ more effectively than the TV regularizer.
Both the compared models SAEiT 0 and TV-L 2 belong to the (EITL) class which considers a linear fidelity term, resorting on the linearized model (17) of the nonlinear forward operator F(σ ). To evaluate the effect of this model simplification on the reconstruction accuracy, we regenerated the input measurements U m by using the linear forward model, that is U m = J σ , and rerun the two algorithms. The reconstructed conductivities reported in Fig. 5, show a visible improvement especially in the square shaped reconstructions while the circular shape profile is accurately reconstructed only by imposing the sparsity penalty in SAEiT 0 . The plots of a cross-section along the conductivity distributions together with the ground truth conductivity are reported both in case the measurements U m are achieved by the non-linear forward operator F(σ ) (Fig. 4), and by its linearization (Fig. 5).
The error metric U , defined as , determines how well the reconstructed conductivity σ * matches the input voltages U m obtained with respect to a given ground truth distribution σ GT . For both algorithms, this  error reduces as expected from 1 × 10 −2 , for measurements U m = F(σ GT ) generated by the nonlinear forward model, to 1 × 10 −18 for U m = J σ GT generated by the linearized model, thus confirming the weakness of the linearized fitting term.

Example 3: Robustness to Noisy Measurements
In order to evaluate the robustness of the proposed method to noisy measurements we corrupted the voltage measurements generated from the solution of the forward model by additive white Gaussian noise. In particular, considering the PCSM-test setup, we generated the noisy measurement vector U m following the perturbed acquisition model (11), by adding to the difference of voltage a vectorn ∼ N (0, s 2 ) of Gaussian noise characterized by standard deviation s and zero-mean. We assume as data quality measure the Signal-to-Noise Ratio in dB, defined by the following formula SNR(U m ) = 10 log 10 which represents the ratio between the standard deviation of the data with respect to the noise variation.
To a given set of SNR values S N R(U m ) = {30, 40, 60}d B is corresponding a set of noise standard deviations s, which can be expressed in terms of percentage value of the standard variation of the measurements std(U m ), that is s = {3.16, 1, 0.1}% of std(U m ). Figure 6 reports the reconstructions from noisy measurements U m for decreasing noise perturbations, left to right, applying TVPIM (first row) and SAEiT η (second row). From a visual inspection we can see that SAEiT η algorithm outperforms the TV-L 2 method, even for severe noise corruptions on the input data U m . This is partially motivated by the fact that the space-variant η-map, even if generated from noisy measurements U m , is well recovered, as illustrated in Fig. 6, third row. This is confirmed as well by the relative errors s reported in Fig. 6, which improve, especially for TV-L 2 , as the noise level decreases.

Example 4: Generic Reconstruction Mesh
In the previous examples the η-map was preliminary exploited to define the backward mesh so that to take advantage of the detected objects for mesh edges alignment to the object boundaries. When, instead, a generic backward mesh is used, then the reconstruction quality inevitably worsens. However, thanks to the shrinkage property that characterizes both TV regularizer as well as the non-convex penalty proposed in SAEiT η , the performance remains of high level. In order to validate this performance, we consider the reconstruction of the PCSM-test setup, on a generic backward mesh with n V = 2625 vertices and n T = 5056 triangles, illustrated in Fig. 7 (left), where the ground truth conductivity is overimposed in false colors. As expected, the straight boundaries of the rectangular inclusion do not overlay exactly on the mesh edges, thus, the perimeter of the object is longer than the original profile. The effects of how this generic backward mesh influences the η-map generation are illustrated in Fig. 8, where σ * 0 obtained by applying SAEiT 0 (Fig. 8, left column) and the η-map detected from it (Fig. 8, middle column) are reported for L = 8, L = 16 and L = 32 electrodes in the first, second and third row, respectively. The measurements U m have been linearly obtained in order to validate the actual improvement due to the increasing number of electrodes, avoiding the inaccuracy due to the linearized fidelity term.
In Fig. 8 (right) the 3D view of the reconstructed conductivities σ * obtained by SAEiT η are shown. Comparing the reconstruction results in Fig. 8 by rows, we can notice how the quality of the η-map as well as of the reconstructed conductivity σ * increases with number of electrodes used to generate the measurements. This is confirmed by the reconstruction error s reported in the first row of Table 2 which decreases almost by the same factor as the increase of the number of electrodes. In the second row of Table 2 we considered the same reconstruction problem with measurements U m given by the nonlinear operator U m = F(σ ). The improvement in s follows the same behavior when the number of electrodes increases from L = 8 to L = 16. However, increasing the number of electrodes to L = 32 there is no longer an improvement that justifies the higher cost of the new actual setup.
The comparison between the two rows in Table 2 highlights that the linearization error, of which the results in the second row are affected, eliminates the reconstruction improvement obtained by increasing the number of electrodes. Nevertheless, this limitation stays in line with other works [2,14], where the authors claim that the improvement in the reconstruction for L = 32 electrodes is only marginal with respect to the reconstruction with measurement obtained from L = 16 electrodes.
Finally, the reconstructed conductivities σ * obtained by the TVPIM and SAEiT η methods are compared and shown in Fig. 9, together with the attained error values, in case L = 16 and U m = F(σ ). A cross-section comparison along a vertical line through the middle of the reconstructed conductivity is shown in Fig. 9 on the right. The shrinkage of the boundary length, which is a property that mostly characterizes the TV regularizer, and tends to round the corners of the sharp inclusion, is well visible in the TVPIM reconstruction in Fig. 9 (left) and thus further influences the accuracy of the reconstruction results.

Example 5: Á-map Construction by Auxiliary Image
In the previous examples the numerical experiments have been conducted in a blind scenario [case (b) in Sect. 4.1], thus the previous η-maps have been constructed by the procedure outlined at the beginning of Sect. 6.
In contrast, this last example illustrates how to embed a priori information on the domain content, represented by a photographic image, to construct an ad-hoc backward mesh and to design a proper η-map. This is identified as case (a) in Sect. 4.1. Several applications can be framed in this scenario. For example, in tissue engineering, an interesting goal is to measure mineral deposition by differentiated stem cells cultured in a hydrogel scaffold. The position of the object of interest inside the domain is known, while the useful information is hidden and encoded in the conductivity distribution.
Given a photograph of the experimental setup as shown in Fig. 10, the η-map and the backward mesh are constructed by the following steps: • Region segmentation of the interior of the domain (tank) to detect and label the background and the contained objects; • Edge detection of the inner N detected regions and boundaries extractions (closed contour curves) C = {C s } N s=1 ; • Backward mesh construction by scaling of the interior of the tank onto a unitary circle; • Building of the η-map by mapping the set C onto the backward mesh. Figure 10a shows a tank with a sensor array of 16 electrodes and two inclusions, a hydrogel scaffold with smooth conductivity distribution and a constant conductivity triangular shaped object, immersed in a buffer for cell cultures that keeps the pH of the experiment under control. In Fig. 10b the inner inclusions have been segmented and labelled; Fig. 10c shows the boundary curves generated by the edge detection procedure on the segmented image. Finally, Fig. 10d shows the detected space variant η-map imposed on the backward mesh which contains n V = 1301 vertices and n T = 2504 triangles.
The simulation was conducted considering an EIT setup composed of 16 electrodes and an injected current with intensity 0.1 mA; the background conductivity was set to The GT conductivity distribution which aims to be reconstructed is illustrated in Fig. 11, imposed on an ad hoc backward mesh generated from the photograph of the experimental setup in Fig. 10.
Reconstructions of the conductivity variation distributions with respect to an empty tank are shown in Fig. 12. A visual inspection reveals that TV-L 2 generates some artifacts around the rubber triangle and completely flattens the variation inside the scaffold, while SAEiT η well reconstructs both the sharp step variation and the inhomogeneous smooth distribution, which is confirmed by the lower s value for SAEiT η .

Conclusion
In this paper, a spatially adaptive non-convex variational model is proposed for dealing with the reconstruction of conductivity distributions in EIT problems that present sharp as well as inhomogeneous variations, as encountered in situations of bioengineering/medical interest. The reconstructions are enhanced by incorporating prior structural information into the regularization through an η-map which is either automatically detected, in case that no a priori knowledge is given, or derived by a photographic image of the acquisition setup. The η-map weights differently two regularization terms. For the reconstruction of sharp variations a non-convex penalty is introduced, while a convex generalized Tikhonov-type regularizer is applied in smoothly varying regions. Experiments demonstrate how the non-convex penalty can promote sparsity better than the convex 1 norm penalty (TV). However, since the matrix J does not have full-column rank, a convex non-convex formulation to choose the parameter a in such a way that the total cost functional J in (18) is convex is not possible if φ is separable as in (19), see [20]. Future directions will consider non-convex non-separable penalty functions as proposed in [20], that do maintain strong convexity of the cost functional J for any matrix J . This would ensure the uniqueness of the regularized solution and the convergence of the ADMM optimization algorithm. However, it will involve a strategy to set the matrix of free parameters suitably derived to face the specific EIT problem.
Finally, from the conducted experiments emerged the crucial role of the backward mesh and the benefit to adapt it to the structural information of the sought for distribution. Therefore, another interesting future direction could be to incorporate a space-adaptive refinement of the backward mesh, driven by a dynamic η-map. Including prior information about the boundary location is feasible in tissue engineering applications, where a picture of the culture well clearly shows the boundaries of the scaffold but misses the crucial information embedded in the conductivity distribution. In a future work, the authors plan to exploit the presented method to measure and map the mineral deposition by alive differentiated stem cells in real time.
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.