Reconstruction of a piecewise constant conductivity on a polygonal partition via shape optimization in EIT
Elena Beretta Stefano Micheletti Simona Perotto Matteo Santacesaria
math.AP math.NA 35R30, 65N21 (Primary), 49Q10 (Secondary)



Reconstruction via shape optimization in EIT]Reconstruction of a piecewise constant conductivity on a polygonal partition via shape optimization in EIT


In this paper, we develop a shape optimization-based algorithm for the electrical impedance tomography (EIT) problem of determining a piecewise constant conductivity on a polygonal partition from boundary measurements. The key tool is to use a distributed shape derivative of a suitable cost functional with respect to movements of the partition. Numerical simulations showing the robustness and accuracy of the method are presented for simulated test cases in two dimensions. 




Primary 35R30, 65N21; Secondary 49Q10
[
    M. I. Krivoruchenko^1,2,3
    December 30, 2023
=============================





§ INTRODUCTION


Electrical Impedance Tomography (EIT) is a noninvasive technique, which aims to detect the conductivity inside a body from voltage and current boundary measurements.
The mathematical problem arising from EIT, known as the inverse conductivity problem, was introduced for the first time by A.-P. Calderón in the early 80's <cit.>. Even though it was first motivated by an application in geophysical prospecting <cit.>, EIT has been having big impact also in medical imaging and nondestructive testing of materials <cit.>.

The conductivity problem can be stated mathematically as follows.
Consider a bounded domain Ω⊂^d, with d=2,3, equipped with an electrical conductivity σ∈ L^∞(Ω) such that σ(x) ≥λ >0.  The corresponding Neumann-to-Dirichlet (ND) or current-to-voltage map is the operator _σ: H_0^-1/2(∂Ω) → H_0^1/2(∂Ω), defined by

    _σ(g) = u|_∂Ω,

where H^s_0(∂Ω) = { f ∈ H^s(∂Ω) : ∫_∂Ω f  ds = 0}, g ∈ H_0^-1/2(∂Ω) and u is the unique H^1(Ω)-weak solution of the Neumann problem for the conductivity equation

    {[     -∇· (σ∇ u) = 0,     in Ω,;        σ∂ u/∂ν =g,     on ∂Ω, ].

where ν is the unit outward normal to ∂Ω,
satisfying the normalization condition

    ∫_∂Ω u  ds = 0.
 
The following inverse boundary value problem arises from this framework.




Inverse conductivity problem. Given _σ, find σ in Ω.




Since the seminal paper by A.-P. Calderón, much interesting mathematics has been developed in order to address the issues of uniqueness, stability and reconstruction for this problem. Concerning uniqueness and reconstruction, we mention the breakthrough results in <cit.>.

The conductivity problem is severely ill-posed as was noted by G. Alessandrini in <cit.>. Despite a-priori smoothness assumptions on the unknown conductivity, a logarithmic-type continuous dependence of the conductivity on the data is the best possible one <cit.>. This fact makes crucial the analysis of the instability and of suitable regularization strategies in order to 
obtain successful computational reconstructions.  
Several recovery methods and procedures have been developed in the last decades.
Without being exhaustive, the possible approaches to reconstruction can be divided into two main streams:


  * iterative methods, based on ad-hoc regularization strategies; 

  * direct methods, where an explicit reconstruction formula of the solution is used.

The first group includes variational-type methods, which reduce the inverse problem to a minimization problem for a least-squares constrained type functional with a suitable regularization. A pioneering paper in this direction is represented by <cit.>, which applies one step of a Newton method with a constant conductivity as an initial guess. In <cit.> the authors introduce a Mumford-Shah type functional, in <cit.> a  level set representation and a total variation regularization is introduced, while in <cit.> an augmented Lagrangian method is proposed. All these methods  are particularly suited to recover piecewise constant conductivities.

Concerning direct methods, we would like to mention the factorization method <cit.>, the D-bar method <cit.>, the enclosure method <cit.>, and the monotonicity method <cit.>.

Finally, statistical inversion has shed interesting insights into EIT reconstruction as well <cit.>.
 
Despite the impressive progress, there remains a big interest in developing new algorithms that take advantage of a-priori information arising from applications. 
An a-priori assumption physically relevant in many applications is to assume the conductivity to be of the form 

    σ=∑_j=1^Nσ_jχ_P_j,

where 𝒫={P_j}_j=1^N is a polygonal (polyhedral) partition of the background body Ω, 
χ_P denoting the characteristic function associated with the generic region P⊂Ω.
Such an assumption arises, for example, in geophysics, medical imaging, and nondestructive testing of materials, where the body under investigation contains regions, represented by the subdomains {P_j}_j=1^N, with different electrical properties.   


Additionally, this kind of a-priori information restores the well-posedness of the inverse problem; in particular, in the case of a given known partition, Lipschitz dependence estimates of the coefficients from the data can be shown <cit.>.
From the results obtained for the Helmholtz equation in <cit.> and <cit.>, where Lipschitz stability holds, we expect that a similar result should be true also when the partition is unknown.
As shown in <cit.>, a crucial role to prove Lipschitz stability  is played by the differentiability properties of the Neumann-to-Dirichlet map with respect to motions of the partition and by the derivation of an explicit formula for the derivative. This is also a crucial step towards reconstruction if we use an optimization approach to solve the inverse conductivity problem.

In this paper, we consider a cost functional, J(σ), representing the L^2-norm of the difference between the potential due to the applied current and the measured potential on the boundary. This functional will be minimized in the class of conductivities σ of the form (<ref>). 
To solve this minimization problem, we introduce an iterative gradient type method which requires the computation of the shape derivative of the functional J(σ). This derivative has been obtained rigorously by several authors in the case of a single sufficiently smooth inclusion ω⊂Ω, i.e. for σ = σ_1 χ_ω+σ_2χ_Ω\ω (see, for example,  <cit.>).

In our case, to implement the optimization procedure, we need to differentiate J(σ) with respect to variations of a partition. 
In <cit.>, the authors derived, for the first time, a rigorous formula for the shape derivative of the functional J(σ) in the planar case and for conductivities of the form

    σ=σ_1χ_P+σ_2χ_Ω\ P

where P is a polygon strictly contained in Ω. The shape derivative is expressed in terms of an integral over the boundary of P and has, surprisingly, exactly the same form as the one derived in <cit.> in the case of a smooth interface, despite the presence of singularities of the gradient of the solutions to the conductivity equation at the vertices of P. The extension of this boundary formula seems not to be possible in the case of an arbitrary partition since the singularity of the gradient might become too strong at the vertices <cit.>.

In this paper, we follow the idea suggested in <cit.> and in <cit.> of using a more general distributed shape derivative of the functional expressed in terms of  an integral  over Ω. The advantage of this formula is twofold: on the one hand, it allows to consider piecewise constant conductivities on very general partitions; on the other hand, it is numerically more accurate (we refer to <cit.> for a thorough comparison of the two formulas from a numerical point of view). In <cit.>, the authors establish the distributed shape derivative in the case of a single measurable conductivity inclusion strictly contained in Ω by using a Lagrangian approach. Here, we establish the formula computing directly and rigorously the derivative of J in terms of the material derivative of the solution to a certain boundary value problem (see Lemma 2.1). The formula is valid for any partition 𝒫 in dimension d = 2,3. Successively, we take advantage from this result to implement our reconstruction procedure in the two dimensional case, based on a gradient type method.

The reconstruction algorithm we present is very similar to the one introduced in <cit.>. Nevertheless, there are two major differences in the implementation, which greatly affect the numerical results. The first one is a regularization step, applied at each iteration, where the number of the sides of each polygon in the partition 𝒫 changes in order to preserve a uniform length. This considerably reduces the artifacts that typically appear in EIT reconstructions. The second major difference lies in the choice of the descent direction for the shape of the partition. While this was done by solving an additional variational problem, we propose a more direct computation which exploits the assumptions made on the conductivity.

The plan of the paper is the following one. In Section 2, we state the problem with the main assumptions, derive the shape derivative of the cost functional and show the equivalence with the boundary shape derivative established in <cit.> for suitable two-dimensional partitions. In Section 3, we describe the reconstruction algorithm. In Section 4 we present some numerical examples which corroborate the reliability and the accuracy of the proposed approach. 




§ MATHEMATICAL FRAMEWORK



 §.§ Main Assumptions
 Let Ω⊂^d, d = 2,3, be a bounded domain with Lipschitz boundary. Let = {P_j}_j=1^N be a polytopal partition of Ω, i.e., P_j are open bounded polytopes (i.e., polygons in 2D, polyhedra in 3D) such that:

    ⋃_j = 1^ NP_j = Ω,      P_j ∩ P_k = ∅   for  j ≠ k.


Let N ∈ be an integer with N>1, and λ >0 be a positive real number. We define the space L^∞(Ω,N,λ) as the collection of conductivities σ∈ L^∞(Ω) such that σ(x) ≥λ > 0 for all x ∈Ω, and such that there exists a polytopal partition = {P_j}_j=1^N' with N' ≤ N, such that σ can be written as

    σ = ∑_j=1^N'σ_j χ_P_j,   with σ_l ≠σ_m  if  P_l  is  adjacent to  P_m.


Now let σ̂, σ∈ L^∞(Ω,N,λ). We denote by σ̂ the unknown conductivity and by σ a (generally) different one, which will be used in the reconstruction scheme.











Let M be the number of measurements. For 1 ≤ j ≤ M, let g_j ∈ H_0^-1/2(∂Ω) be a given function representing the applied current density on ∂Ω, and f_j ∈ H^1/2(∂Ω) the corresponding measurement of the voltage on ∂Ω. 
More precisely, f_j =  û_j|_∂Ω, where û_j is a solution of

    {[     -∇· (σ̂∇û_j)   = 0,     in Ω,;        σ̂∂û_j/∂ν =g_j,     on ∂Ω. ].

In order to recover σ̂, we minimize the following Dirichlet least-squares fitting cost functional:

    J(σ) = 1/2∑_j=1^M∫_∂Ω|u_j -f_j|^2  ds,

for σ∈ L^∞(Ω,N,λ), where the state function u_j solves

    {[     -∇· (σ∇ u_j)   = 0,     in Ω,;        σ∂ u_j/∂ν =g_j,     on ∂Ω, ].

with the normalization condition

    ∫_∂Ω u_j   ds= ∫_∂Ω f_j  ds.


Notice that the functional J(σ) depends on the values {σ_j}_j=1^N and the partition = {P_j}_j=1^N, where σ = ∑_j=1^Nσ_j χ_P_j. The next subsections are devoted to the computation of the gradient of J(σ) with respect to the above variables.



 §.§ Gradient of J with respect to  – the shape derivative

The gradient of J with respect to the partition  is actually a shape derivative. We thus want to compute the shape derivative ⟨∇_ J, U ⟩ of the functional J at the partition  in the direction of a vector field U = (U_1,…,U_d). We assume that U ∈ W^1,∞(^d) and U = 0 in a neighborhood of ∂Ω. This derivation has already been carried out in <cit.> in the case of smooth inclusions, and in <cit.> for a single polygonal inclusion. We present here a more general formula that is valid for any finite partition of a domain.



In order to compute the derivative, consider the transformation Φ_t( x) = x + t U( x): ^d →^d, as smooth as U. We assume that t ≤ 1/(2U_W^1,∞) so that Φ_t^-1 exists globally. For instance, Φ_t is a piecewise affine function that moves the nodes of the partition  (i.e., the vertices of the P_j's). Note that Φ_t |_∂Ω = Id, Id denoting the identity mapping. This assumption is not a real restriction since in EIT the boundary ∂Ω is always assumed to be known and fixed. Nevertheless, it is crucial in the derivation of the main results of this section. 

The following matrix-valued functions will be useful in the following:

    A(t)    = (DΦ_t^-1) (DΦ_t^-1)^T(DΦ_t),
       = . d A/dt|_t=0 = (U)I - (DU+DU^T),

where D Φ_t^-1 and DU are the Jacobian matrices of Φ_t^-1 and U, respectively.
  
Consider the deformed partition _t = Φ_t() and σ_t = σ∘Φ_t^-1 the corresponding conductivity. Let G(t) be defined as

    G(t) = J(σ_t).


Let now u be any of the solutions u_j, with Neumann data g, where g is any of the currents g_j, j=1,…,M, and let f be a boundary measurement corresponding to g. First, we need to study the material derivative of the solution u.


The solution u to problem (<ref>) has a material derivative u̇∈ H^1(Ω) that solves

    ∫_Ωσ∇u̇·∇ w dx = - ∫_Ωσ∇ u ·∇ w dx,     ∀ w ∈ H^1(Ω),

with the normalization condition ∫_∂Ωu̇ ds = 0.

 We follow Step 1 and 2 of the proof of <cit.>.
Let u_t be the solution of

    {[ -∇· (σ_t ∇ u_t)  = 0,     in Ω,;    σ_t ∂ u_t/∂ν = g,     on ∂Ω, ].

with the normalization condition ∫_∂Ω u_t ds = ∫_∂Ω f ds.
Then the transported solution ũ_t = u_t ∘Φ_t solves the variational equation:

    ∫_Ωσ A(t) ∇ũ_t ·∇ w dx -∫_∂Ωg w ds= 0,     ∀ w ∈ H^1(Ω).

Subtracting to (<ref>) the variational equation solved by u and dividing by t, we find

    ∫_Ωσ A(t) ∇ũ_t -∇ u/t·∇ w dx = ∫_ΩσI- A(t)/t∇ u ·∇ w dx,   ∀ w ∈ H^1(Ω).

Using ũ_t- u as test function we obtain

    1/2min_x ∈Ωσ(x) ∇ũ_t - ∇ u/t_L^2(Ω)≤A(t) - I/t_∞∇ u _L^2(Ω).

Thus, we have found that (ũ_t - u)/t is bounded in H^1(Ω). Therefore, the sequence is weakly convergent in H^1(Ω) and its weak limit is the material derivative u̇ of u. Passing to the limit in (<ref>), we find that u̇ solves the desired variational formulation (<ref>). 

By the Lax-Milgram lemma, since the right-hand side is an H^-1(Ω) function (because of our assumptions on U), every solution to the variational problem (<ref>) lies in H^1(Ω). In particular, the trace on ∂Ω is well defined.

Actually, we have strong convergence. Plugging w = (ũ_t - u)/t into (<ref>), we obtain

    ∫_Ωσ A(t) ∇ w ·∇ w  dx  = ∫_ΩσI- A(t)/t∇  u ·∇ w  dx = B_1,t+B_2,t,

where 

    B_1,t = ∫_Ωσ ( A(t)-I) ∇  w ·∇ w  dx   and   B_2,t= ∫_ΩσI- A(t)/t∇ũ_t ·∇ w  dx.

Thanks to the weak convergence of (ũ_t - u)/t, we obtain

    B_1,t→ 0   and   B_2,t→ - ∫_Ωσ∇ u ·∇u̇ dx   as  t → 0.

Now, the variational formulation (<ref>) yields B_2,t→∫_Ωσ∇u̇·∇u̇ dx.
So we have found that ∇ũ_t -∇ u/t converges strongly to ∇u̇ in L^2(Ω). Using the normalization conditions for u̇, u and ũ_t (which coincides with u_t on ∂Ω), we get the strong convergence of ũ_t - u/t to u̇ in H^1(Ω), via the Poincaré inequality. 




Now we can derive the formula for the shape derivative of the functional J.

We have

    ⟨∇_ J , U ⟩ = ∑_j=1^M ∫_Ωσ∇ u_j ·∇ z_j dx,

where z_j solves 

    {[          -∇· (σ∇ z_j)        = 0,     in Ω,;             σ∂ z_j/∂ν = f_j-u_j,     on ∂Ω, ].

with the normalization ∫_∂Ωz_j ds = ∫_∂Ω f_j ds.


We need to compute the derivative d G/dt|_t=0, where G is defined in (<ref>). Let u_j,t, ũ_j,t be defined as u_t, ũ_t in the proof of Lemma <ref>, with Neumann data g_j. By the assumption that Φ_t|_∂Ω = Id, we have

    G(t) = 1/2∑_j=1^M ∫_∂Ω(u_j,t-f_j)^2 ds =   1/2∑_j=1^M ∫_∂Ω(ũ_j,t -f_j)^2 ds.

Note that

    G(t)-G(0)/t=∑_j=1^M ∫_∂Ω( ũ_j,t-u_j/t)( ũ_j,t+u_j/2 -f_j) ds.

Using the strong convergence ũ_j,t-u_j/t→u̇_j in H^1(Ω) and the convergence ũ_j,t→ u_j in L^2(∂Ω) (which follows from the fact that ũ_j,t = u_j,t on ∂Ω and u_j,t→ u_j in H^1(Ω)), we can pass to the limit as t → 0 in (<ref>), thus obtaining 

    ⟨∇_ J , U ⟩  = dG/dt|_t=0 = ∑_j=1^M ∫_∂Ω (u_j-f_j)u̇_j ds,

where u̇_j is the material derivative of u_j.
Using the variational formulation of the adjoint z_j, 

    ∫_Ωσ∇ z_j ·∇ w dx + ∫_∂Ω (u_j-f_j) w dσ = 0,     ∀ w ∈ H^1(Ω),

and of the material derivative, (<ref>), we find

    ⟨∇_ J , U ⟩  = - ∑_j=1^M ∫_Ωσ∇ z_j ·∇u̇_j dx = ∑_j=1^M ∫_Ωσ∇ u_j ·∇ z_j dx,

which is the desired formula.


We now recall also the gradient of J with respect to {σ_j }_j=1^N, that will be used in the reconstruction algorithm:


    d J/d σ_j = ∑_k=1^M ∫_P_j∇ u_k ·∇ z_k dx,    j = 1,…,N,

where z_k solves (<ref>) with a normalization (e.g., ∫_∂Ωz_k ds = ∫_∂Ω f_k ds).



 §.§ Equivalence of the distributed and of the boundary integral formulas in the case of a single polygonal inclusion

We emphasize that, in the case of a single poly­gonal inclusion, the shape derivative of the functional J(σ) has been computed rigorously in <cit.>, and expressed as an integral on the boundary of such an inclusion.  It is unclear if this boundary representation of the shape derivative is still valid in the case of an arbitrary partition.

Let P be a polygon strictly contained in Ω and let

    u^+=u|_Ω\P̅,     z^+=z|_Ω\P̅.

and 

    u^-=u|_ P,     z^-=z|_P.



Assume Ω=P∪(Ω\ P), with P a polygon strictly contained in Ω, and let σ|_P = k and σ = 1 outside P.  Then, we have the following equivalent formulas for the shape derivative of the functional J at P in the direction given by U:

    .d G/dt|_t=0   =∑_j=1^M ∫_Ωσ( x) ∇ u_j( x) ·∇ z_j(x) dx 
        = (k-1) ∑_j=1^M∫_∂ P(1/k∂ u_j^+/∂ν∂ z^+_j/∂ν+∇_τu_j ·∇_τz_j)U_ν ds,

where ν is the unit outward normal vector to ∂ P, ∇_τ is the tangential gradient, U_ν = U·ν, functions u_j,z_j satisfy (<ref>) and (<ref>), respectively.



We can assume, without loss of generality, M=1 and denote by u and z  the solutions corresponding to the datum f. Let B_ε be the union of balls of radius ε centered at each vertex of P,
and let us consider the splitting 

    ∫_Ωσ∇ u ·∇ z dx=∫_Ω\ B_εσ∇ u ·∇ z dx+∫_B_εσ∇ u ·∇ z dx.

Denote by 

    Ω_ε^+=Ω\P∪ B_ε,   Ω_ε^-=P\ B_ε,

and observe that, by standard regularity results (see, for example, <cit.>), it can be shown that u^+,z^+∈ H^2(Ω_ε^+), and  u^-,z^-∈ H^2(Ω_ε^-). Then,  applying Green's formula in Ω_ε^±, observing that U has compact support in Ω, and using the following identity

    ∇ u ·∇ z    = -(b) + (U·∇ u) Δ z +(U·∇ z) Δ u
       =-(b)     in Ω_ε^+∪Ω_ε^- =Ω\B̅_ε,

where 

    b = (U·∇ u) ∇ z + (U·∇ z) ∇ u- (∇ u ·∇ z) U,
 
we easily derive

    ∫_Ω\B̅_εσ∇ u ·∇ z dx=∫_∂ P\B̅_ε[σ b]·ν ds+∫_∂ B_εσ b·ν ds,

where we use notation [f] = f^+-f^- to denote the jump off across ∂ P. Using the transmission conditions satisfied by u and z across ∂ P, we end up with the following relation

    ∫_Ω\B̅_ε   σ∇ u ·∇ z dx = 
       (k-1) ∫_∂ P(1/k∂ u^+/∂ν∂ z^+/∂ν+∇_τu ·∇_τz)U_ν ds
       -(k-1) ∫_∂ P∩B̅_ε(1/k∂ u^+/∂ν∂ z^+/∂ν+∇_τu ·∇_τz)U_ν ds+∫_∂ B_εσ b·ν ds,

where U_ν=U·ν.
Hence,

    ∫_Ω   σ∇ u ·∇ z dx = 
       (k-1) ∫_∂ P(1/k∂ u^+/∂ν∂ z^+/∂ν+∇_τu ·∇_τz)U_ν ds
    
       -(k-1) ∫_∂ P∩B̅_ε(1/k∂ u^+/∂ν∂ z^+/∂ν+∇_τu ·∇_τz)U_ν ds_R_1
       +∫_∂ B_εσ b·ν ds_R_2+∫_B_εσ∇ u ·∇ z dx_R_3
    
       =(k-1) ∫_∂ P(1/k∂ u^+/∂ν∂ z^+/∂ν+∇_τu ·∇_τz)U_ν ds+R_1+R_2+R_3.


From <cit.>, we have the following upper bounds of the gradients of u and z in a neighbourhood of the vertices:

    |∇ u (x)|≤ C |x-x̃|^α-1,    |∇ z (x)| ≤ C |x-x̃|^α-1,

for some constant C>0, α > 1/2 and x sufficiently close to a vertex x̃, with x≠x̃.  This, jointly with the regularity assumptions on the vector field U, implies that 

    ∫_Ω   σ∇ u ·∇ z dx = 
    (k-1) ∫_∂ P(1/k∂ u^+/∂ν∂ z^+/∂ν+∇_τu ·∇_τz)U_ν ds+O(ϵ^2α-1),

and, since α>1/2, letting ϵ→ 0 in the last equation, we finally obtain

    ∫_Ωσ∇ u ·∇ z dx=(k-1) ∫_∂ P(1/k∂ u^+/∂ν∂ z^+/∂ν+∇_τu ·∇_τz)U_ν ds,

which ends the proof.



It is straightforward to check that the boundary formula derived in <cit.> (and hence also the last proposition) extends to the case of a finite number of well separated, polygonal inclusions at a positive distance to the boundary of Ω.



The formula for the shape derivative of the functional J can also be obtained following the Lagrangian approach as in <cit.>. 
In fact, defining the following Lagrangian 

    Ł̃(u,z,t,U)    = 1/2∫_∂Ω|u -f|^2 dσ+ ∫_Ωσ A(t)  ∇ u ·∇ z dx - ∫_∂Ω g z dσ,

where A(t) is defined as in (<ref>), while functions u and z satisfy (<ref>) and (<ref>), respectively corresponding to a given Neumann datum g, 
it is possible to show that 

    dG/dt|_t=0= ∂/∂ tŁ̃(u,z,t,U)|_t=0= ∫_Ωσ(x) ∇ u(x) ·∇ z(x) dx.





§ RECONSTRUCTION ALGORITHM

We use the descent gradient method to solve the minimization problem involving the functional J in (<ref>). Although the algorithm holds in dimension 2 and 3, we focus on the planar case, since computational experiments are in 2D.	

The algorithm adopted in this paper differs from reconstruction algorithms available in the literature (<cit.>)
because of two major features, which are separately addressed in the next sections.



 §.§ Regularization
 The first novelty adopted in the reconstruction algorithm consists in 
increasing or decreasing the number of vertices of the polygons in the partition in order to impose some regularity on the reconstruction. This simple trick considerably improves the reconstruction quality, smoothing out artifacts and irregularities that often characterize EIT reconstructions (see Section 4 below for practical examples).



 §.§ Construction of the descent direction
 The second distinguishing feature of the adopted reconstruction algorithm 
concerns the computation of the descent direction associated with the partition vertices. A standard approach consists in solving a discretized version of the equation

    ⟨θ^k, δθ⟩ + ⟨∇_ J(σ^k), δθ⟩ = 0     for every δθ∈ X,

where θ^k ∈ X is the descent direction at iteration k, σ^k is the conductivity obtained at the iteration k (see Section <ref> for a precise definition), while X ⊂ W^1,∞(Ω,^2). There are several approaches to discretize and solve equation (<ref>). We refer to <cit.> for a thorough discussion on this subject. Here we propose a different, more straightforward, derivation. Since the partition is defined by its vertices, we compute the descent direction for each vertex individually, and then we update the partition. The formula we use at iteration k to define the descent direction, θ_l^k ∈^2, for a given vertex, V_l^k, l=1,…,N^k_V, with N^k_V the number of the partition vertices at the iteration k, is:

    θ_l^k = -(⟨∇_ J(σ^k), U^k_l,1⟩, ⟨∇_ J(σ^k), U^k_l,2⟩),

where the vector fields U^k_l,1,U^k_l,2 are chosen as follows:


  (1) U^k_l,1,U^k_l,2∈ W^1,∞(Ω,^2) with  support strictly contained in Ω

  (2) U^k_l,1,U^k_l,2 piecewise linear on the edges of the partition and such that

    U^k_l,1(V_j^k) = (δ_jl,0),      U^k_l,2(V_j^k) = (0,δ_jl),


where δ_jl indicates the Kronecker delta. 
With this choice, the vector 

    (⟨∇_ J(σ^k), U^k_l,1⟩, ⟨∇_ J(σ^k), U^k_l,2⟩)
 
represents, at the iteration k, the gradient of the functional J with respect to the position of a single vertex, V_l^k. 
In practice, in view of a finite element discretization, we choose the vector fields 
U^k_l,1 and U^k_l,2 as 

    U^k_l,1 = (φ^k_l,0),      U^k_l,2 = (0,φ^k_l),

where φ^k_l is the hat function associated with the node V^k_l of a coarse mesh that contains both the vertices and the edges of the partition  but with no additional nodes on the sides of the polygons. More precisely, ϕ^k_l (V^k_j) = δ_jl, and it is piecewise linear. Examples of coarse meshes adapted to a partition can be seen in every plot in Section 4. Thus, any other hat function corresponding to such a coarse mesh would produce the same results, since the shape derivative is supported on the edges of the partition.

We remark that we need to assume the knowledge of the number N of polygons in the partition, since the algorithm is unable to make changes in the topology. A possible approach to overcome this limitation, using topological derivatives and a level set formulation, has been studied in <cit.>.



 §.§ The reconstruction algorithm

Let N^k_V be the number of vertices of the partition ^k = { P^k_j}_j=1^N at iteration k, and { V^k_l }_l =1^N^k_V be the corresponding set of vertices. Let N^k_j be the number of vertices of the polygon P^k_j, for j=1,…,N. At each iteration, we consider the conductivity σ^k = ∑_j=1^N σ_j^k χ_P^k_j. In the algorithm we introduce the parameters δ_1, δ_2 >0 as threshold for the distances between consecutive vertices in the regularization step and a tolerance tol >0 is chosen to control the size of the shape gradient as a stopping criterion. Finally, the step sizes α^j_k > 0 and β_k >0 can be fixed or obtained by line search. However, in all numerical experiments  below, we keep them fixed.



The stopping criteria for the algorithm may be improved following the error analysis in <cit.>. For the sake of simplicity, in the present work, we decided to focus on the new shape updates in the algorithm, neglecting other issues already investigated in earlier works. The 
inclusion of an advanced stopping criterion as well as of a robust line search algorithm for the step sizes constitute possible topics for a future investigation.



§ NUMERICAL TESTS

In this section, we show the effectiveness of the proposed reconstruction scheme. A series of numerical experiments have been carried out in order to assess numerically that the algorithm enable us to recover simultaneously the partition and the values of the conductivity for some significant configurations.

For the sake of simplicity, we consider Ω = (0,1)^2. For each test, we create a mesh adapted to the unknown conductivity that is used to provide the boundary data. For the reconstruction, a coarse mesh adapted to the approximate partition is generated at each iteration. A refined mesh is also constructed in order to compute the state and the adjoint problem.



The boundary data are generated in the following manner. We set σ∂ u/∂ν = 1 on one of the four sides of Ω, -1 on another side and 0 elsewhere. In this way, we obtain 6 independent current patterns, corresponding to 4 electrodes, one for each side. Then, we divide each side of Ω in half. We set σ∂ u/∂ν = 1 on one half, -1 on another half side and 0 elsewhere. This gives 28 independent current patterns, corresponding to 8 electrodes. We iterate this procedure one more time. In this way, we construct sets of 6, 28 and 120 boundary data, corresponding to 4, 8 and 16 electrodes, respectively. See Figure <ref> for a scheme of the electrode position corresponding to the data sets.

We also add a uniform noise to the data. More precisely, given a noiseless boundary measurement f_j ∈ H^1/2(∂Ω), j=1, …,M, the noisy data f̃_j is obtained by adding to f_j a uniform noise in the following way:

    f̃_j (x) = f_j (x) +  εf_L^2(∂Ω),

where x ∈∂Ω is a boundary vertex of the mesh that generated f_j and ε is a uniform random real in (-γ,γ), where γ> 0 is chosen according to the noise level.
To measure the noise level, we use the relative error on the boundary in the L^2-norm, that is the following quantity:

    √(∑_j=1^M f̃_j - f_j^2_L^2(∂Ω))/√(∑_j=1^M  f_j^2_L^2(∂Ω)).


The regularization parameters, δ_1 and δ_2, are chosen experimentally. Since the initial guess is always a regular polygon (or a collection of regular polygons), we choose δ_1 = α_1 δ, δ_2 = α_2 δ, where δ is the length of the side of the initial guess, with α_1 < 1, α_2 > 1.5. This is done in order to have, at each iteration, a partition with edges of similar length.

All the computations are performed using FreeFem++ <cit.>.



 §.§ Shape reconstruction








In this first set of examples, we assume to know the values of the conductivity, and we only reconstruct the shape of the partition. 

In Figure <ref>, we consider a piecewise conductivity which is equal to 10 on a convex pentagon and to 1 in the background. The reconstruction is carried out using 6 boundary measurements. The initial guess is a regular polygon of 14 sides (top-left). We present a noiseless reconstruction, with (top-right) and skipping (bottom-left) the regularization step, and a reconstruction with regularization when 3% of noise is added to the data (bottom-right). The regularization parameters are set to δ_1 = 0.7 δ and δ_2 = 1.8 δ, where δ is the length of the side of the initial guess. Figure <ref> shows that the three reconstructions well identify the shape with a comparable precision.

A non-convex polygon is considered in Figure <ref>, where the target conductivity is 10 inside and 1 in the background. Here, we employ 28 boundary measurements. The initial guess is a regular polygon of 24 sides (top-left). For this phantom, we present a noiseless reconstruction with (top-right) and without (bottom-left) the regularization step, and a reconstruction post regularization in the presence of a 3% of noise added to the data. The regularization parameters are chosen as δ_1 = 0.85 δ and δ_2 = 1.8 δ, where δ is the length of the side of the initial guess. We notice here that the regularization step helps to reconstruct more precisely the non-convex part of the unknown. It is also interesting to observe that the shape is well identified also in the noisy case.

In Figure <ref>, we consider the so-called heart and lung phantom <cit.>, with background conductivity 1, two ellipses with conductivity 0.5 and a disk with conductivity 2 (both ellipses and the disk are approximated with 16-sided polygons). The initial guess coincides with three identical regular polygons of 16 sides each (top-left). The regularization parameters are δ_1 = 0.9 δ and δ_2 = 1.8 δ, where δ is the length of the side of one of the initial guess polygons. The reconstruction is done using 6 boundary measurements. Here, the difference between the regularized (top-right) and the non-regularized (bottom-left) reconstruction is significative. Moreover, the reconstruction is very robust to the 3% noise added to the data (bottom-right).


The three phantoms just presented have different contrast. It is therefore natural to study the dependence of the reconstruction on the contrast. We observed that the algorithm converges faster in case of higher contrast, yet the reconstruction quality is the same. For this reason we decided to not include reconstructions of the same phantom with different contrast values.




 §.§ Reconstruction from a misplaced initial guess


With the example in Figure <ref>, we show the robustness of the algorithm in recovering an unknown shape from a misplaced initial guess. The target is a piecewise conductivity with value 10 inside a square and 1 in the background. The initial guess is an octagon, far from the exact square
(top-left). The regularization parameters are set to δ_1 = 0.8 δ and δ_2 = 1.7 δ, being δ the length of the side of the octagon. The algorithm is able to recover the square (bottom-right). The striking result strongly depends on the choice of the regularization parameters. Indeed, after the first few iterations, the octagon becomes a triangle (top-right). This reduction in the degrees of freedom avoids potential degeneracy of the shape. The triangle eventually changes and becomes a square when approaching the target (bottom-left).



 §.§ Simultaneous reconstruction



In this section, we present reconstructions obtained with the full algorithm, i.e., the conductivity values are also unknown and updated at each iteration.

We focus on the heart and lung phantom, this configuration being characterized by more interesting features compared with the other ones. The background is assumed to be known and is equal to 1, while the coefficients to be recovered are the lungs, with value 0.5, and the heart, with value 2. The reconstructions presented in Figure <ref> are obtained from 28 boundary data. The initial guess (top-left) identifies the shapes the algorithm acts on, while three different values are adopted for the initial conductivity during the reconstruction procedure. The regularization parameters are set to δ_1 = 0.9 δ and δ_2 = 1.8 δ independently of the run, with δ the length of the side of one of the initial guess polygons.

The first reconstruction deals with noiseless data. The values of the initial guess are 0.55 in the lungs and 2.05 in the heart. The reconstructed values are 0.49 and 2.05, respectively, and the shape is well reconstructed for the three inclusions (top-right). The second reconstruction starts from the same initial guess, by adding a 5% noise to the data. The reconstructed conductivity values are 0.37 in the lungs and 2.06 in the heart. The shape of the lungs is better recovered than the heart (bottom-left).

The last reconstruction is essentially blind, where we assume that the initial guess is provided by a constant conductivity equal to 1 (the background). Moreover, data contains 1% of additive noise. Despite the challenging setting, some features are correctly recovered. The reconstructed values are 0.44 in the lungs and 0.94 in the heart. The shape and the values of the lungs are well reconstructed, whereas this is not the case for the heart  (bottom-right), due to the lack of a priori information and the central position of the middle polygon in the initial guess.



 §.§ Sensitivity to the number of measurements



Using the heart and lung phantom (always with value 0.5 in the lungs, 2 in the heart, and 1 in the background), we check now how the reconstructions change using a different set of measurements. The regularization parameters are always chosen as δ_1 = 0.9 δ and δ_2 = 1.8 δ, with δ the length of the side of one of the initial guess polygons.

In Figure <ref>, we present shape reconstructions (conductivity values are known) starting from noiseless data sets of 6, 28 and 120 boundary measurements. While there is an evident improvement passing from 6 to 28 data, the reconstruction quality obtained from 28 and 120 measurements looks very similar in this example. This is due to the ill-posedness of the problem that limits the resolution.



 §.§ Sensitivity to the noise level



In this section, we show how the algorithm is stable to noise in the measurements. In Figure <ref>, we present reconstructions using the full algorithm (unknown shape and values) for the heart and lung phantom. The regularization parameters are always chosen as δ_1 = 0.9 δ and δ_2 = 1.8 δ, being δ the length of the side of one of the initial guess polygons.

The initial guess has values 0.7 for the lungs and 1.5 for the heart, far from the exact values 0.5 and 2, respectively. Data contain 0.5% (top-right), 5% (bottom-left) and 20% (bottom-right) of additive noise, respectively. The tolerance in the stopping criterion is chosen experimentally, according to the noise level: 0.004 for both 0.5% and 5% of noise, 0.02 for 20% of noise. 
We observe that the shape of the lungs is well recovered even with very noisy data, whereas the heart is poorly recovered, due to the values of the initial guess and the central position of the middle disk. The reconstructed values exhibit a low sensitivity to the noise level. They are 0.4 in the lungs and 1.51 in the heart for both 0.5% and 5% of noise, while we obtain 0.37 in the lungs and 1.49 in the heart when considering 20% of noise.









§ CONCLUSIONS


We have presented a new shape optimization approach that is able to well identify a piecewise constant conductivity on a polytopal partition in electrical impedance tomography. Despite the ill-posedness of EIT, the a-priori assumption on the conductivity to be piecewise constant on a polygonal partition regularizes the problem. This fact, coupled with the effectiveness of the algorithm, allows us to recover an unknown partition (or at least a subset of it) even in the case of noisy data and a wrong initial guess. This is possible due to the use of a distributed shape derivative, to a new regularization step and a new computation of the descent direction.  The algorithm performs better when a good approximation of the values of the conductivity is available. When the coefficients are not known, we noticed that it takes more iterations to approximate their values than to identify the shapes. These values could be obtained, for instance, from a one-step reconstruction <cit.> and used as a first guess to recover the interfaces more accurately (see also <cit.> for a hybrid one-step method). The algorithm can be accelerated using more advanced minimization methods and more efficient forward solvers. Possible approaches to accelerate the minimization include Newton-type methods and also reduced order models to speed up the forward and adjoint problems. This would be crucial with a view to practical applications.



§ ACKNOWLEDGMENTS

E. Beretta and M. Santacesaria thank the New York University Abu Dhabi for its kind hospitality that permitted a further development of the present research. The work of E. Beretta and M. Santacesaria was partially supported by GNAMPA (Gruppo Nazionale per l'Analisi Matematica, la Probabilità e le loro Applicazioni).

elsarticle-num

10
url<#>1urlprefixURL href#1#2#2 #1#1

Calder'on1980
A.-P. Calderón, On an inverse boundary value problem, in: Seminar on
  Numerical Analysis and its Applications to Continuum Physics (Rio
  de Janeiro, 1980), Soc. Brasil. Mat., Rio de Janeiro, 1980, pp. 65–73.

Chambers2006
J. E. Chambers, O. Kuras, P. I. Meldrum, R. D. Ogilvy, J. Hollands, Electrical
  resistivity tomography applied to geologic, hydrogeologic, and engineering
  investigations at a former waste-disposal site, Geophysics 71 (6) (2006)
  B231–B239.

Karhunen2010
K. Karhunen, A. Seppänen, A. Lehikoinen, P. J. M. Monteiro, J. P. Kaipio,
  Electrical resistance tomography imaging of concrete, Cement and Concrete
  Research 40 (2010) 137–145.

Sylvester1987
J. Sylvester, G. Uhlmann, A global uniqueness theorem for an inverse boundary
  value problem, Annals of Mathematics 125 (1987) 153–169.

Novikov1988
R. G. Novikov, A multidimensional inverse spectral problem for the equation
  -Δψ+(v(x)-Eu(x))ψ = 0, Functional Analysis and Its
  Applications 22 (4) (1988) 263–272.

Nachman1996
A. I. Nachman, Global uniqueness for a two-dimensional inverse boundary value
  problem, Annals of Mathematics 143 (1996) 71–96.

Alessandrini1988
G. Alessandrini, Stable determination of conductivity by boundary measurements,
  Applicable Analysis 27 (1988) 153–172.

mandache2001
N. Mandache, Exponential instability in an inverse problem for the
  Schrödinger equation, Inverse Problems 17 (5) (2001) 1435.

Cheney1990a
M. Cheney, D. Isaacson, J. Newell, S. Simske, J. Goble, NOSER: An algorithm
  for solving the inverse conductivity problem, International Journal of
  Imaging Systems and Technology 2 (2) (1990) 66–75.

Rondi2001
L. Rondi, F. Santosa, http://dx.doi.org/10.1051/cocv:2001121Enhanced
  electrical impedance tomography via the Mumford-Shah functional, ESAIM
  Control Optim. Calc. Var. 6 (2001) 517–538.
http://dx.doi.org/10.1051/cocv:2001121
  doi:10.1051/cocv:2001121.
<http://dx.doi.org/10.1051/cocv:2001121>

Chung2005electrical
E. T. Chung, T. F. Chan, X.-C. Tai, Electrical impedance tomography using level
  set representation and total variational regularization, Journal of
  Computational Physics 205 (1) (2005) 357–372.

Chan2004
T. F. Chan, X.-C. Tai, Level set and total variation regularization for
  elliptic inverse problems with discontinuous coefficients, Journal of
  Computational Physics 193 (1) (2004) 40–66.

Chen1999
Z. Chen, J. Zou, An augmented Lagrangian method for identifying discontinuous
  parameters in elliptic systems, SIAM Journal on Control and Optimization
  37 (3) (1999) 892–910.

Bruhl2000
M. Brühl, M. Hanke, Numerical implementation of two non-iterative methods for
  locating inclusions by impedance tomography, Inverse Problems 16 (2000)
  1029–1042.

Kirsch2008
A. Kirsch, N. Grinberg, The Factorization Method for Inverse Problems, Oxford
  University Press, USA, 2008.

Siltanen2000
S. Siltanen, J. Mueller, D. Isaacson, An implementation of the reconstruction
  algorithm of A. Nachman for the 2-D inverse conductivity problem,
  Inverse Problems 16 (2000) 681–699.

Ikehata2000a
M. Ikehata, S. Siltanen, Numerical method for finding the convex hull of an
  inclusion in conductivity from boundary measurements, Inverse Problems 16 (4)
  (2000) 1043–1052.

Tamburrino2002
A. Tamburrino, G. Rubinacci, A new non-iterative inversion method for
  electrical resistance tomography, Inverse Problems 18 (6) (2002) 1809.

Kaipio2004a
J. Kaipio, E. Somersalo, Statistical and Computational Inverse Problems, Vol.
  160 of Applied Mathematical Sciences, Springer Verlag, 2004.

Alessandrini2005a
G. Alessandrini, S. Vessella, Lipschitz stability for the inverse conductivity
  problem, Advances in Applied Mathematics 35 (2) (2005) 207–241.

beretta2015stable
E. Beretta, M. V. de Hoop, E. Francini, S. Vessella, Stable determination of
  polyhedral interfaces from boundary data for the Helmholtz equation,
  Communications in Partial Differential Equations 40 (7) (2015) 1365–1392.

beretta2016stable
E. Beretta, M. V. de Hoop, F. Faucher, O. Scherzer, Inverse boundary value
  problem for the Helmholtz equation: quantitative conditional Lipschitz
  stability estimates., SIAM Journal on Mathematical Analysis 48 (6) (2016)
  3962–3983.

hettlich1998determination
F. Hettlich, W. Rundell, The determination of a discontinuity in a conductivity
  from a single boundary measurement, Inverse Problems 14 (1) (1998) 67.

afraites2007shape
L. Afraites, M. Dambrine, D. Kateb, Shape methods for the transmission problem
  with a single measurement, Numerical Functional Analysis and Optimization
  28 (5-6) (2007) 519–551.

pantz2005sensibilite
O. Pantz, Sensibilité de l'équation de la chaleur aux sauts de
  conductivité, Comptes Rendus Mathematique 341 (5) (2005) 333–337.

beretta2017differentiability
E. Beretta, E. Francini, S. Vessella, Differentiability of the Dirichlet to
  Neumann map under movements of polygonal inclusions with an application to
  shape optimization, SIAM Journal on Mathematical Analysis 49 (2) (2017)
  756–776.

piccinini1972
L. C. Piccinini, S. Spagnolo, On the Hölder continuity of solutions of
  second order elliptic equations in two variables, Annali della Scuola Normale
  Superiore di Pisa-Classe di Scienze 26 (2) (1972) 391–402.

bruce1974
R. Bruce Kellogg, On the Poisson equation with intersecting interfaces,
  Applicable Analysis 4 (2) (1974) 101–129.

Laurain2016
A. Laurain, K. Sturm, Distributed shape derivative via averaged adjoint method
  and applications, ESAIM: Mathematical Modelling and Numerical Analysis 50 (4)
  (2016) 1241–1267.

Giacomini2017
M. Giacomini, O. Pantz, K. Trabelsi, Certified Descent Algorithm for shape
  optimization driven by fully-computable a posteriori error estimators, ESAIM:
  Control, Optimisation and Calculus of Variations 23 (3) (2017) 977–1001.

Hiptmair2015
R. Hiptmair, A. Paganini, S. Sargheini, Comparison of approximate shape
  gradients, BIT Numerical Mathematics 55 (2) (2015) 459–485.

Giacomini2015
M. Giacomini, O. Pantz, K. Trabelsi, An a posteriori error estimator for shape
  optimization: application to EIT, in: Journal of Physics: Conference
  Series, Vol. 657, IOP Publishing, 2015, p. 012004.

Ito2008
K. Ito, K. Kunisch,
  http://epubs.siam.org/doi/abs/10.1137/1.9780898718614Lagrange
  Multiplier Approach to Variational Problems and Applications, Society for
  Industrial and Applied Mathematics, 2008.
http://dx.doi.org/10.1137/1.9780898718614
  doi:10.1137/1.9780898718614.
<http://epubs.siam.org/doi/abs/10.1137/1.9780898718614>

bellout1992stability
H. Bellout, A. Friedman, V. Isakov, Stability for an inverse problem in
  potential theory, Transactions of the American Mathematical Society 332 (1)
  (1992) 271–296.

Dogan2007
G. Dogan, P. Morin, R. H. Nochetto, M. Verani, Discrete gradient flows for
  shape optimization and applications, Computer Methods in Applied Mechanics
  and Engineering 196 (37) (2007) 3898–3914.

Hintermuller2008
M. Hintermüller, A. Laurain, Electrical impedance tomography: from topology
  to shape., Control & Cybernetics 37 (4).

Hecht2012
F. Hecht, New development in FreeFem++, Journal of Numerical Mathematics
  20 (3-4) (2012) 251–266.

Mueller2002
J. Mueller, S. Siltanen, D. Isaacson, A direct reconstruction algorithm for
  electrical impedance tomography, IEEE Transactions on Medical Imaging 21 (6)
  (2002) 555–559.

Hamilton2016
S. J. Hamilton, J. M. Reyes, S. Siltanen, X. Zhang, A hybrid segmentation and
  D-bar method for electrical impedance tomography, SIAM Journal on Imaging
  Sciences 9 (2) (2016) 770–793.




