Hybrid topological derivative and gradient-based methods for electrical impedance tomography

The use of electric signals to obtain information about the environment is an old survival strategy in nature. For example, nocturnal fishes found in fresh water (the gymnotiformes, commonly known as South American knifefishes) produce electric fields for navigation. Technology exploiting this idea has been developed based on impedance imaging. The impedance imaging problem consists in producing an image of the electromagnetic properties of a medium by applying electric currents to its exterior surface and measuring voltages on it. Information about the internal electrical properties of a body can be used for nondestructive materials testing [17], locating mineral deposits [34], imaging multiphase fluid flow [39], tracing contaminants [35] or for medical imaging [27]. The range of medical applications is wide, because different tissues have different electromagnetic properties [10]. For example, we can think of monitoring for lung problems (embolies, clots, accumulation of fluids) or blood flow (internal bleeding, heart function), screening for breast cancer, determining the boundary between dead and living cells, detecting temperature changes in hyperthermia treatments and so on.

The electromagnetic properties of interest are typically the electric conductivity and the permittivity. The reconstruction problem consists in finding an approximation of the electric conductivity and the permittivity in the interior of the body from boundary measurements. This problem is nonlinear and ill posed. Given arbitrary data solutions may not exist. When they exist, small changes in the measured data can lead to large changes in the interior conductivity or permittivity. If the boundary potentials are known for all current patterns in reasonable function spaces, then there is no second smooth conductivity distribution which fits these data [31, 32]. It has been shown that positive conductivities can be uniquely recovered from the Dirichlet to Neumann map [38]. In practice, we only have information from a finite number of electrodes. For a revision of theoretical results, see [4] and [30].

A wide variety of strategies to solve this type of inverse problem have been proposed; see the reviews [4, 12] and references therein. When the conductivity is almost constant, we may resort to linear approximations such as backpropagation methods [2], Calderón's approach [11], moment techniques [1] or one-step Newton methods [10, 19]. Between the noniterative methods addressing the full nonlinear problem, we may cite the layer-stripping method [37], schemes based on Nachman's uniqueness proofs [30], techniques to find discontinuities [5], factorization and MUSIC algorithms to locate cavities [23]. Iterative [15, 26] and adaptive [33] schemes are other possibilities. Recently, level set techniques have been extended to characterize interfaces between regions of different conductivity [13, 28]. In [24] and [25], level set techniques are initialized by means of a topological sensitivity analysis. Available methods typically impose technical restrictions on the admittivity or possible objects to proceed, as the assumption about almost constant conductivity in linear approximations. Most techniques focus on the inverse conductivity problem. Methods tracking discontinuities or objects with sharp interfaces usually consider piecewise-constant conductivity [5, 13, 24, 26, 28]. Some of them reconstruct the support of inclusions or just small cavities without being able to make precise the value of the conductivity [5, 23]. The iterative methods presented in [13, 24, 25, 28] show reconstruction tests involving several hundred or thousand iterations. Noniterative methods may suffer from apparent geometrical restrictions [37].

In this paper, we solve the impedance imaging problem by a combination of topological derivative-based schemes and gradient techniques. We distinguish two cases. Initially, we consider the presence of discontinuities corresponding to inclusions of a different material. The geometry of the set of inclusions and the values of the electromagnetic parameters in them have to be identified. There are no restrictions on their size or distribution, and their parameters may be space dependent. The admittivity of the background is assumed to be known and may also be space dependent. In the second framework, we do not consider the presence of sharp interfaces. The conductivities and permittivities of the whole medium are considered to be unknown and vary smoothly. We propose a technique to reproduce their spatial variations.

Both imaging problems can be reformulated as inverse scattering problems and then recast as constrained optimization problems. In the first, we have a matrix with known electrical properties containing several domains with unknown properties. Computing the topological derivative of the shape functional to be minimized we find the first approximation of the number, size and location of the domains. A gradient method produces the first guess of the values of the parameters. This procedure can be iterated to improve the approximation. In the second, we just minimize the functional with respect to the unknown variable parameters by a gradient method. For an introduction to topological derivative techniques the reader may consult [6, 7, 20, 21] and references therein. Time-harmonic problems with spatially dependent coefficients are considered in [7] and general time-dependent problems in [8].

The paper is organized as follows. In section 2, we formulate the impedance imaging problem and give a constrained optimization reformulation. Section 3 considers domains with inclusions. Both the exterior and the interior admittivities may be constant or smooth-space-dependent functions. We suggest a reconstruction technique combining topological derivatives of the shape functional for fixed parameters and a gradient-based method for the conductivities and permittivities for fixed domains. Section 4 proposes a method to reconstruct smoothly varying conductivities and permittivities. (In this case, we do not consider the presence of permittivities.) In section 5, we show some imaging tests using the numerical algorithms introduced in sections 3 and 4. Section 6 summarizes our conclusions and comments on extensions to more realistic settings.

The electrical properties of interest are typically the electric conductivity and permittivity. Inside the body Ω, the electric potential $u(\mathbf {x})$ satisfies

$$\begin{eqnarray} &&\nabla \cdot \gamma (\mathbf {x},\omega ) \nabla u =0, \end{eqnarray} \tag{ 1 }$$

where γ = σ + ıωε, with σ being the electric conductivity, ε the electric permittivity and ω the angular frequency of the applied current. A derivation from Maxwell's equations is given in [12]. When ω = 0, we are left with an inverse conductivity problem.

Currents are applied through electrodes placed on the boundary ∂Ω of the body. Different boundary conditions have been proposed to account for this effect. The simplest model imposes a 'continuous' Neumann condition

$$\begin{eqnarray} &&\gamma {\partial _\mathbf {n}u } = j\quad \mbox{on $\partial \Omega $}, \end{eqnarray} \tag{ 2 }$$

and measures the voltage on the boundary

$$\begin{eqnarray} &&u= V_{{\rm meas}}\quad \mbox{on $\partial \Omega $}. \end{eqnarray} \tag{ 3 }$$

Here, $\mathbf {n}$ denotes the outward unit normal and j stands for the outward-pointing normal component of the current density generated by the electrodes on the surface. The Neumann problem (1)–(2) admits solutions when the current satisfies the compatibility condition ∫_∂Ωj dl = 0, which is precisely the law of conservation of charge. To select a reference potential u, we impose ∫_∂Ωu dl = 0. The measured voltage also satisfies ∫_∂ΩV_meas dl = 0.

In real experiments, only the values of the current and the voltages at a 'discrete' set of electrodes is known. In this paper, we focus on the simplified model (1)–(2). Section 6 comments on possible extensions to more realistic settings.

We want to reconstruct the admittivity γ inside Ω from measurements on the boundary. The impedance tomography problem can be reformulated as an optimization problem: minimize

$$\begin{eqnarray} &&J(\gamma ) = {1\over 2} \int _{\partial \Omega } |u - V_{{\rm meas}} |^{2} \,{\rm d}l \end{eqnarray} \tag{ 4 }$$

when u solves (1)–(2).

If we assume that Ω contains a number of inclusions Ω_{i, j}, the admittivity γ is a piecewise function in Ω with discontinuities at the boundaries of the inclusions. We set Ω_i = ∪^d_{j = 1}Ω_{i, j}, Ω_{i, j} being open connected bounded sets satisfying $\overline{\Omega }_{i,l}\ \cap\ \overline{\Omega }_{i,j}=\emptyset$ for l ≠ j; see figure 1. The admittivity in the matrix $\Omega _{e}=\Omega \setminus \overline{\Omega }_{i}$ is γ_e. We define γ_i in Ω_i as γ_i = γ_{i, j} in Ω_{i, j}. To simplify, we assume γ_e to be known. In this case, (4) is replaced by

$$\begin{eqnarray} &&J(\Omega \setminus \overline{\Omega }_{i},\gamma _{i}) = {1\over 2} \int _{\partial \Omega } |u - V_{{\rm meas}} |^{2} \,{\rm d}l, \end{eqnarray} \tag{ 5 }$$

where u is a solution of

$$\begin{eqnarray} \left\lbrace \begin{array}{@{}llll}\nabla \cdot \gamma _e \nabla u =0 & \mbox{ in $\Omega _{e}$}, & \nabla \cdot \gamma _i \nabla u =0 &\mbox{ in $\Omega _{i}$},\\ u^{-} - u^{+} =0 & \mbox{on $\partial \Omega _{i}$}, & \gamma _i \partial _{\mathbf {n}} u^{-} -\gamma _e \partial _{\mathbf {n}} u^{+} =0 &\mbox{ on $\partial \Omega _{i}$}, \\ \gamma _{e} \partial _{\mathbf {n}} u = j &\mbox{on $\partial \Omega $}. \end{array} \right. \end{eqnarray} \tag{ 6 }$$

The unit normal $\mathbf {n}$ points outside Ω_e but inside Ω_i, and u⁻ and u⁺ denote the limit values of u on ∂Ω_i from outside and inside Ω_i, respectively (see figure 1).

Zoom In Zoom Out Reset image size

In the case where the medium Ω_e contains inclusions Ω_{i, j} of a different material, discontinuities in the admittivity are expected at the boundaries of the inclusions. To determine its spatial variation we must optimize (5) with respect to the interior domains and their admittivity. Even if the admittivity is unknown, we may locate the regions where it undergoes notable changes by computing the topological derivative of the shape functional. We refer the reader to [6] and references therein for an introduction to topological derivative methods in inverse scattering and shape optimization. Once the first guess for Ω_i = ∪^d_{j = 1}Ω_{i, j} is found, we may approximate the admittivity γ_i optimizing $J(\Omega \setminus \overline{\Omega }_i,\gamma _i)$ with respect to γ_i for Ω_i fixed.

Below we will describe the resulting procedure in three steps. First, we compute the required topological derivatives and propose the first approximation for Ω_i. Then, we use a gradient technique to optimize with respect to γ_i. Finally, we explain how to iterate both procedures to improve the guesses of Ω_i and γ_i. For simplicity, we take Ω to be a bounded set of $\mathbb {R}^2.$ The method applies with straightforward modifications in three dimensions.

3.1. Topological derivative with respect to the domain

The first guess for the inclusions is found evaluating topological derivatives of the cost functional with respect to the domain.

The topological derivative of a shape functional $\mathcal J(R)$, $\mathcal{R} \subset \Omega$ measures its sensitivity when an infinitesimal hole is removed from $\mathcal{R}$ [36]. Let us consider a small ball $B_\varepsilon (\mathbf {x})=B(\mathbf {x},\varepsilon ),\,\mathbf {x}\in \mathcal{R}$, and the domain $\mathcal{R}_\varepsilon :=\mathcal{R} \setminus \overline{B}_\varepsilon (\mathbf {x})$. The topological derivative of $\mathcal{J}(\mathcal{R})$ is defined as

$$\begin{eqnarray} &&D_T(\mathbf {x},\mathcal{R}):=\lim _{\varepsilon \rightarrow 0}\frac{\mathcal{J}(\mathcal{R}_\varepsilon )-\mathcal{J}(\mathcal{R})}{\mathcal{V}(\varepsilon )}, \end{eqnarray} \tag{ 7 }$$

where $\mathcal{V}(\varepsilon )$ is a positive and decreasing function satisfying $\lim _{\varepsilon \rightarrow 0}\mathcal{V}(\varepsilon )=0$, for which this limit is finite and nonzero. The expansion

$$\begin{eqnarray} &&\mathcal{J}(\mathcal{R}_\varepsilon ) = \mathcal{J}(\mathcal{R}) + \mathcal{V}(\varepsilon ) D_T(\mathbf {x},\mathcal{R}) + o(\mathcal{V}(\varepsilon )) \end{eqnarray} \tag{ 8 }$$

for ε small implies that the functional decreases if we remove from $\mathcal{R}$ small balls centered at points $\mathbf {x}$ at which $D_T(\mathbf {x},\mathcal{R}) <0$. This suggests a minimizing strategy. We start by computing the topological derivative of the cost functional with $\mathcal{R}=\Omega$. Then, we select as the first approximation of the inclusions a region Ω_ap⊂Ω at which the topological derivative takes large negative values, so that $\mathcal{J}(\Omega \setminus \overline{\Omega }_{ap}) < \mathcal{J}(\Omega )$. We may update the approximation Ω_ap by setting $\mathcal{R}= \Omega \setminus \overline{\Omega }_{ap}$, computing the topological derivative of $\mathcal{J}(\Omega \setminus \overline{\Omega }_{ap})$ and adding to Ω_ap points at which the new topological derivative is negative and large. This strategy only allows us to add points to the approximations Ω_ap. To be able to add and remove points from our approximated domain at each iteration, we need a value for the topological derivative inside Ω_ap. Formula (7) gives a value for the topological derivative at points $\mathbf {x}\in \Omega \setminus \overline{\Omega }_{ap}.$ An extension to all $\mathbf {x}\in \Omega$ can be defined for our particular type of functional following [7]. The extended topological derivative is defined as

$$\begin{eqnarray} &&D_T(\mathbf {x},\Omega \setminus \overline{\Omega }_{ap}):= \lim _{\varepsilon \rightarrow 0}\frac{\mathcal{J}_{\varepsilon } (\Omega \setminus \overline{\Omega }_{ap})- \mathcal{J}(\Omega \setminus \overline{\Omega }_{ap})}{\mathcal{V}(\varepsilon )},\qquad \mathbf {x}\in \Omega , \end{eqnarray} \tag{ 9 }$$

where $\mathcal{J}_{\varepsilon }(\Omega \setminus \overline{\Omega }_{ap}):= {1\over 2} \int _{\partial \Omega } |u_{\varepsilon } - V_{{\rm meas}} |^{2} {\rm d}l$ and u_ε solves (6) for the domains

$$\begin{equation*} \Omega _{e}:= \Omega \setminus \overline{\Omega }_{i} \qquad {\rm where}\ \ \Omega _i:=\left| \begin{array}{@{}ll@{}}\Omega _{ap} \setminus \overline{B}_\varepsilon (\mathbf {x}),&\quad \mathbf {x}\in \Omega _{ap},\\ [1.5ex] \Omega _{ap} \cup B_\varepsilon (\mathbf {x}),&\quad \mathbf {x}\in \Omega \setminus \overline{\Omega }_{ap}. \end{array} \right. \end{equation*}$$

For the electrical impedance tomography problem in 2D, $\mathcal{V}(\varepsilon )$ can be taken as the measure of a two-dimensional ball: $\mathcal{V}(\varepsilon ) = \pi \varepsilon ^2.$ When $\mathbf {x}\in \Omega \setminus \overline{\Omega }_{ap},$ formula (9) agrees with (7). An expansion of the form (8) still holds, now for $\mathbf {x}\in \Omega$. By adding to Ω_ap points where the topological derivative is large and negative (as they decrease the shape functional to first order) and removing points at which the topological derivative is large and positive (since they increase the shape functional) we expect to generate new approximations for which the cost functional takes smaller values.

An explicit expression for the topological derivative of the cost functional (5) in terms of forward and adjoint fields is given in the next result. The proof is sketched in the appendix. Expressions for the topological derivative using definition (7) with Ω_ap = ∅ in real-valued inverse conductivity problems with piecewise constant conductivity have been obtained in [24, 25] expanding asymptotically a modified functional.

Theorem 3.1. The topological derivative (9) of the cost functional (5) is given by

$$\begin{eqnarray} D_T(\mathbf {x},\Omega \setminus \overline{\Omega }_{ap},\gamma _i)&=& \mbox{ Re} \bigg [ { 2 \gamma _e({\bf x})(\gamma _{i}({\bf x})-\gamma _{e}({\bf x})) \over \gamma _e({\bf x})+ {\gamma _{i}({\bf x})}} \nabla u(\mathbf {x}) \nabla \overline{p}(\mathbf {x}) \bigg ], \quad \mathbf {x}\in \Omega \setminus \overline{\Omega }_{ap}, \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} D_T(\mathbf {x},\Omega \setminus \overline{\Omega }_{ap} ,\gamma _i)&=& \mbox{ Re} \bigg [ { 2 \gamma _i({\bf x})(\gamma _{i}({\bf x})-\gamma _{e}({\bf x})) \over \gamma _i({\bf x}) + \gamma _{e}({\bf x}) } \nabla u(\mathbf {x}) \nabla \overline{p}(\mathbf {x}) \bigg ], \quad \mathbf {x}\in \Omega _{ap}, \end{eqnarray} \tag{ 11 }$$

where u is the solution of the forward problem (6) and p is the solution of the adjoint problem

$$\begin{equation} \fl \left\lbrace \begin{array}{@{}llll}\nabla \cdot \overline{\gamma }_e \nabla p =0 & \quad \mbox{ in $\Omega _{e}$}, & \nabla \cdot \overline{\gamma }_i \nabla p =0 & \quad \mbox{ in $\Omega _{i}$}, \\ [1.3ex] p^{-} - p^{+} =0 &\quad \mbox{ on $\partial \Omega _i$}, & \overline{\gamma }_i \partial _{\mathbf {n}} p^{-} - \overline{\gamma }_e \partial _{\mathbf {n}} p^{+} =0 & \quad \mbox{ on $\partial \Omega _{i}$}, \\ [1.3ex] \overline{\gamma }_{e} \partial _{\mathbf {n}} p = V_{{\rm meas}}-u &\quad \mbox{ on $\partial \Omega $}, \end{array} \right. \end{equation} \tag{ 12 }$$

with Ω_i = Ω_ap.

The admittivity γ_i enters formulas (10) and (11), but, in principle, it is unknown. To initialize our procedure, we select the first approximation of γ_i by perturbing γ_e: γ⁰_i = γ_e + and set Ω_ap = ∅. Then, we evaluate D_T(x, Ω, γ_i) setting γ_i = γ⁰_i. The first approximation to Ω_i is given by the set of points of Ω in which D_T(x, Ω, γ⁰_i) takes large negative values:

$$\begin{equation*} \Omega _{i}^1= \big\lbrace \mathbf {x}\in \Omega | D_{T}\big(\mathbf {x},\Omega ,\gamma _{i}^0\big) < - C_{1} \big\rbrace \end{equation*}$$

for some value C₁ > 0. We refer the reader to [7] for some guidelines about the selection of this constant C₁.

3.2. Variation with the admittivity

The topological derivative $D_{T}(\mathbf {x},\Omega ,\gamma _i^0)$ of the cost functional (5) provides the first approximation of the regions where the admittivity differs from γ_e. Our goal now is to correct the value of γ_i in Ω¹_i. We seek to minimize

$$\begin{eqnarray} &&J(\gamma _i) = {1\over 2} \ \int _{\partial \Omega } |u - V_{{\rm meas}} |^{2} \,{\rm d}l \end{eqnarray} \tag{ 13 }$$

when u is a solution of (6) with Ω_i = Ω¹_i and γ_i = γ⁰_i. To do so, we resort to a gradient method. We update γ⁰_i setting γ¹_i = γ_i⁰ + ηϕ, choosing η and ϕ in such a way that ${{\rm d} J(\gamma _{i}^{0}+\eta \phi ) \over {\rm d} \eta }<0$. Finding an explicit expression for this derivative, we obtain a formula for the corrector profile ϕ in each iteration. Related differential calculus and formulas (Frechet derivatives with respect to the conductivity, for instance) are developed in [3, 4, 15, 22]. A broad review of techniques to reconstruct admittivities and conductivities can be found in [4, 12]. The gradient iteration we propose is easy to combine with topological derivatives and produces reasonable reconstructions in simple geometries (see section 5). The specific derivatives and corrector functions we need are given below.

Theorem 3.2. The derivative with respect to η of the function J(η) := J(γ⁰_i + ηϕ) is

$$\begin{eqnarray} &&{{\rm d} J\big(\gamma _{i}^0+\eta \phi \big) \over {\rm d} \eta } \big |_{\eta =0}= \mbox{ Re}\left[\int _{\Omega _{i}^{1}} \phi \, \nabla u \nabla \overline{p} \,{\rm d}\mathbf {z}\right], \end{eqnarray} \tag{ 14 }$$

where u and p are the solutions of the forward (6) and adjoint (12) problems with Ωi = Ω1i = ∪d1j = 1Ω1i, j and γi = γ0i. Choosing

$$\begin{equation} \phi _{1} = - \nabla \overline{u} \nabla p, \quad \mathbf {x}\in \Omega _{i}^1, \end{equation} \tag{ 15 }$$

and η > 0, we ensure J(γ¹_i) < J(γ_i⁰) for γ¹_i = γ_i⁰ + ηϕ₁, if η is small enough.

Proof. The derivative with respect to η is

$$\begin{equation} \frac{{\rm d}}{{\rm d}\eta }\,J(\eta )\bigg |_{\eta =0} = \mbox{Re}\left[ \int _{\partial \Omega } (\overline{u-V_{{\rm meas}}}) \dot{u} \, {\rm d}l\right], \end{equation} \tag{ 16 }$$

where $\dot{u} =\frac{{\rm d}}{{\rm d}\eta } u_{\eta }\big |_{\eta =0}$ and u_η is the solution of (6) for γ_i = γ⁰_i + ηϕ. We replace (6) by its variational formulation:

$$\begin{equation} u_{\eta }\in H^1_{m}(\Omega ), \quad | \quad \displaystyle b(\eta ; u_{\eta }, v)= \ell (v),\qquad \forall v\in H^1_{m}(\Omega ), \end{equation} \tag{ 17 }$$

where H¹_m(Ω) is the subspace of H¹(Ω) formed by functions of zero mean,

$$\begin{eqnarray*} &&b(\eta ;u,v) := \displaystyle \int _{\Omega _e} \gamma _{e} \nabla u\nabla {\overline{v}} \,{\rm d}\mathbf {z}+\int _{\Omega _i} \big(\gamma _{i}^0+\eta \phi \big) \nabla u\nabla {\overline{v}} \,{\rm d}\mathbf {z}\quad \forall u, v\in H^1_{m}(\Omega ),\\ &&\ell (v) := \displaystyle \int _{\partial \Omega } j \overline{v} \,{\rm d}l, \qquad \forall v\in H^1_{m}(\Omega ). \end{eqnarray*}$$

We avoid computing $\dot{u}$ by introducing an adjoint problem. Let us set

$$\begin{equation*} \mathcal{L}(\eta ;u,p):=J(\eta ) +\mbox{Re}[b(\eta ;u,p)-\ell (p)],\qquad \forall u,p\in H^{1}_{m}(\Omega ). \end{equation*}$$

Choosing u = u_η,

$$\begin{equation*} J(\eta )=\mathcal{L}(\eta ;u_{\eta },p),\qquad \forall p\in H^{1}_{m}(\Omega ), \end{equation*}$$

and

$$\begin{eqnarray} &&\fl \frac{{\rm d}}{{\rm d}\eta }\,J(\eta )\big |_{\eta =0} = \frac{{\rm d}}{{\rm d}\eta }\mathcal{L}(\eta ; u_\eta ,p) \big |_{\eta =0}= \mbox{Re}\left[\frac{{\rm d}}{{\rm d}\eta }b(\eta ; u,p)\big |_{\eta =0} \right] + \mbox{Re} [b(0;\dot{u},p)] \nonumber \\ &&+ \;\mbox{Re}\left[ \int _{\partial \Omega } (\overline{u -V_{{\rm meas}}}) \dot{u} \, {\rm d}l \right]. \end{eqnarray} \tag{ 18 }$$

We select an adjoint field p satisfying

$$\begin{equation} p\in H^1_{m}(\Omega ), \quad | \quad \displaystyle b(0;v,p)=- \int _{\partial \Omega } (\overline{u -V_{{\rm meas}}})\, v \, {\rm d}l,\quad \forall v\in H^1_{m}(\Omega ). \end{equation} \tag{ 19 }$$

The derivative of J(η) with respect to η, given by (18), reduces then to

$$\begin{equation*} \frac{{\rm d}}{{\rm d}\eta }J(\eta )\bigg |_{\eta =0}= \mbox{Re}\left[ \frac{{\rm d}}{{\rm d}\eta }b(\eta ; u,p) \bigg |_{\eta =0} \right] = \mbox{Re}\left[ \int _{\Omega _{i}^1} \phi \nabla u \nabla \overline{p}\,{\rm d}\mathbf {z}\right]. \end{equation*}$$

 □

We have defined ϕ₁ only in Ω¹_i, but formula (15) holds for all x ∈ Ω. When γ_i is piecewise constant, we expect the resulting γ¹_i to be almost piecewise constant. Another possible choice for γ¹_i is obtained replacing the function ϕ₁ defined in (15) by piecewise constant approximations

$$\begin{equation} \phi _{1} = - \int _{\Omega _{i,j}^1} \nabla \overline{u} \nabla p \,{\rm d}\mathbf {z}, \quad \mathbf {x}\in \Omega _{i,j}^1, \; j=1,\ldots ,d_{1}. \end{equation} \tag{ 20 }$$

This procedure may be iterated to generate a sequence γ^{k + 1}_i = γ_i^k + η_kϕ_k, with ϕ_k given by (15) or (20), where now u and p solve forward and adjoint problems with Ω_i = Ω¹_i.

3.3. Iterative algorithm

In the previous section, we have proposed a mixed topological derivative—gradient-based method to reconstruct admittivities when sharp interfaces are present. It is worthwhile to compare the performance of this hybrid strategy with global methods to reconstruct admittivities without assuming any a priori information on the background. Let us assume now that Ω is made up of a single material and the admittivity γ inside Ω is unknown. We want to reconstruct its spatial variation inside Ω from measurements at the boundary. The gradient technique described in section 3.2 can be adapted to reproduce smooth variations of the admittivity, when applied to the functional J(γ) defined in (4) without any previous knowledge of the background.

Theorem 4.1. The derivative with respect to η of the function J(η) := J(γ^k + ηϕ) is

$$\begin{eqnarray} &&{{\rm d} J(\gamma ^k+\eta \phi ) \over {\rm d} \eta } \big |_{\eta =0}= \mbox{ Re}\left[\int _{\Omega } \phi \, \nabla u_k \nabla \overline{p}_k \,{\rm d}\mathbf {z}\right], \end{eqnarray} \tag{ 23 }$$

where u_k is a solution of (1)–(2) with γ = γ^k, and p_k is a solution of

$$\begin{eqnarray} &&\nabla \cdot \gamma \nabla p =0, \quad \mbox{in $\Omega $}, \qquad \gamma {\partial _\mathbf {n}p } = V_{{\rm meas}}-u, \quad \mbox{on $\partial \Omega $}, \end{eqnarray} \tag{ 24 }$$

with $\gamma =\overline{\gamma }^{k}$, u = u_k. Choosing

$$\begin{equation} \phi _{k} = - \nabla \overline{u}_k \nabla p_k , \quad \mathbf {x}\in \Omega , \end{equation} \tag{ 25 }$$

we ensure J(γ^{k + 1}) < J(γ^k) for γ^{k + 1} = γ^k + ηϕ_k, provided η > 0 is small.

Proof. The proof is similar to that of theorem 3.2, setting now

$$\begin{equation*} \displaystyle b(\eta ;u,v) = \displaystyle \int _{\Omega } (\gamma ^k+\eta \phi ) \nabla u\nabla {\overline{v}} \,{\rm d}\mathbf {z}, \qquad \forall u, v\in H^1_{m}(\Omega ), \end{equation*}$$

$$\begin{equation*} \ell (v) = \displaystyle \int _{\partial \Omega } j \overline{v} \,{\rm d}l, \quad \forall v\in H^1_{m}(\Omega ). \end{equation*}$$

 □

Related differential calculus and formulas are developed in [3, 15]; see also the review [4].

Note that ϕ_k depends on $\mathbf {x}.$ This theorem may be used to generate a sequence γ^{k + 1} = γ^k + η_kϕ_k with ϕ_k given by (25) yielding information on the spatial variation of γ. The resulting algorithm is as follows.

• Choose a value for γ⁰ (it may be generated from the known value at the boundary).
• Set γ^{k + 1} = γ^k + η_kϕ_k with ϕ_k given by (25).
• Stop when the variations in γ^k from one step to the next are small; the functional falls below a threshold or the discrepancy principle ||u^ε_meas − u_k|| ⩽ τε is satisfied.

This algorithm has been tested in section 5, producing reasonable reconstructions in simple geometries, though the hybrid scheme seems to behave better in problems with interfaces. Other choices would be possible, see [3, 14, 15, 40] for instance.

In this section, we test the algorithms proposed in the previous sections. The first examples consider the inverse conductivity problem, obtained when the frequency ω = 0. The permittivity disappears from the problem, which becomes real valued. We focus on piecewise constant conductivities initially. Next, we consider complex-valued examples, and a problem with space-dependent conductivity, used to compare the two approaches described in sections 3 and 4. Prior to the numerical tests we nondimensionalize the problem.

In all our numerical experiments Ω is the unit ball. Synthetic data are generated by solving problem (6) for the following ten current patterns:

$$\begin{equation*} j_\ell (t)=\cos (\ell t),\quad j_{5+\ell }(t)=\sin (\ell t),\qquad \ell =1,\dots ,5,\quad t\in [0,2\pi ). \end{equation*}$$

When σ is piecewise constant, that is,

$$\begin{equation*} \sigma (\mathbf {x})=\left\lbrace \begin{array}{@{}ll}\sigma _e,&\quad \mathbf {x}\in \Omega \setminus \cup _{j=1}^d\Omega _{i,j},\\ \sigma _{i,j},&\quad \mathbf {x}\in \Omega _{i,j},\end{array} \right. \end{equation*}$$

we can use boundary integral methods to solve the direct and adjoint problems (6) and (12). In this case, at the kth step of the iterative method described in section 3.3 we need to define a smooth contour for Ω^k_i. We track the interfaces assuming that each boundary of Ω^k_i can be parameterized as

$$\begin{equation} (c_x,c_y)+r(t)(\cos (2\pi t),\sin (2\pi t)),\quad t\in [0,2\pi ], \end{equation} \tag{ 26 }$$

and solve a least-squares problem to find an approximation of r(t) of the following form:

$$\begin{equation} r(t)\approx a_0+\sum _{m=1}^Ma_m\cos (2\pi m t)+\sum _{m=1}^Mb_m\cos (2\pi m t), \end{equation} \tag{ 27 }$$

see [7] for details. When σ varies smoothly with $\mathbf {x}$, we resort to finite element codes.

For simplicity, in our numerical experiments, objects do not have inner boundaries, like in annular objects, but our method can deal with those kinds of objects. The interested reader may find a reconstruction of an annular region using topological derivatives for a scattering problem in [7]. In our current implementation, using boundary elements we assume that the boundaries are $\mathcal{C}^2$. The approximation (26)–(27) can be substituted by other methods to track the interfaces, for instance by spline curves, easily. Note also that the number of components does not have to be known from the beginning and can change from one iteration to the next (see figure 3). The method can also deal with components that separate or merge without problems (see [6] for an example in a scattering context).

In the first simulations (figure 2–5), we reconstruct the geometry of regions with different real constant conductivities. In all of them the exterior conductivity is σ_e = 1. The measured data were generated by solving the equivalent nonlinear system of integral equations proposed in [16], using a finer mesh and adding a 1% relative error at each observation point. We took 300 grid points on each curve to create synthetic data and 32 for solving the direct and adjoint problems at each step. In all our experiments, we have taken 30 observation points, that is, 30 electrodes, uniformly distributed on the boundary of the unit ball (they are represented by crosses in all figures).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For our first example, we consider the kite-shaped object Ω_i given by the parameterization

$$\begin{equation} \mathbf {x}(t)=\bigg (-\frac{1}{2},-\frac{1}{5}\bigg )+\frac{1}{1600} (380\cos (t)+247\cos (2\,t), 570\sin (t) ),\quad t\in [0,2\pi ), \end{equation} \tag{ 28 }$$

with constant conductivity σ_i = 6. This object is represented in figure 2 with a white solid line. Note that this object is not in the set of functions defined by (26) and (27).

The values of the topological derivative when Ω_i = ∅ at a sampling grid defined in the unit ball are represented in figure 2(a). We see that the largest negative values attained by the topological derivative (dark blue colors in the plots) are concentrated inside the true object. Figure 2(b) shows our initial guess Ω¹_i (denoted in the plot as Ω₁) and the topological derivative computed when Ω_i = Ω¹_i. Some of the subsequent iterations are shown in figures 2(b)–(g). The last plot shows the decreasing values of the cost functional at each iteration. The iterative procedure was finished by the discrepancy principle at the 16th iteration. At each step, the domain was defined using (22). The alternative definition (21) provides exactly the same results since the topological derivative does not take large positive values during the iterations. We obtained a good reconstruction of the size, location and shape of the object taking into account that no a priori information is used.

Let us consider now a geometrical configuration with three objects of circular shapes and different sizes (see figure 3). In this example, σ_i = 8 in all the objects. Computing the topological derivative when Ω_i = ∅, we observe that only the biggest object is detected (see figure 3(a)). Therefore, our initial guess shown in figure 3(b) has only one defect. The computation of the topological derivative when Ω_i = Ω₁ allows us to distinguish the right-most object. In the next iteration, Ω_i has two components (see figure 3(c)). The smaller object is detected when computing the topological derivative for Ω_i = Ω₂ (see figure 3(c)). In the subsequent iterations, we have caught the correct number of defects. Some of them are represented in figures 3(d)–(h). The values of the cost functional versus the number of iterations are given in figure 3(i). We get a satisfactory reconstruction after just 11 iterations without having any a priori information. We are able to detect the correct number of objects, their position and approximate size. A better reconstruction of the shapes is possible after a few more iterations.

In figures 2 and 3, we generated monotone sequences of objects in the sense that Ω^k_i⊂Ω_i^{k + 1} for all k. That is, we used the definition of Ω^{k + 1}_i given in (22). The approximate domains increase slowly to resemble the true object in figure 2. Instead, in the example with three objects, we observe that figure 3(c) already includes a spurious region. All the subsequent subplots include spurious regions that cannot be removed. This is an artifact of the monotone method, though the value of the constants c_k and C_k affects the appearance of such regions. In spite of this, the final approximation provides a reasonable reconstruction of the shape, size and location of the objects. On the other hand, if we observe the values of the topological derivative at the 11th step (see figure 3(h)), we realize that there are regions where it takes positive values close to the two objects at the bottom. This indicates that spurious regions are included and that we should use now a non-monotone method if we want to correct them. Defining now Ω^{k + 1}_i as in (21), the spurious regions would be removed in a few steps.

In figure 4, we revisit the first example taking a larger value for C₁. We will use now the non-monotone method. We compute the topological derivative inside and outside the objects and represent the true object by a solid line and the current approximation by dashed lines. Our initial guess is bigger and contains spurious regions that do not belong to the true object (compare figure 4(b) with figure 2(b)). When computing the topological derivative for Ω_i = Ω¹_i we only detect new regions where its values are large and negative (see figure 4(b)). Therefore, in the next iteration Ω¹_i⊂Ω_i². This also happens in the second and third steps (see figures 4(c) and (d)). However, at the fourth iteration, we detect a region with large positive values, indicating that spurious regions are included in Ω⁴_i (see figure 4(e)) and in the next iteration $\Omega _i^4\not\subset \Omega _i^5$. As the method proceeds, it happens that $\Omega ^k_i\not\subset \Omega _i^{k+1}$ for k = 5, ..., 10 (see figures 4(f) and (g) where we represent the results for k = 5 and k = 8). The computation of the topological derivative when Ω_i = Ω¹¹_i provides regions with large negative values, but no large positive values appear (see figure 4(h)). Therefore, Ω¹¹_i⊂Ω_i¹². Our final reconstruction at the 13th step is shown in figure 4(i) with a dashed line. The values of the cost functional through the iterative procedure are given in figure 4(j). For comparison, we have added in figure 4(k) the reconstruction that we obtained in figure 2(g). We observe that our final reconstruction has a similar quality and also that the values of the cost functional when the iterations stop are almost the same (compare figures 4(j) and (l)). This indicates that our non-monotone method is not very sensitive to the initial guess.

Zoom In Zoom Out Reset image size

In our next example, we illustrate the performance of our method for the full problem, that is, for the reconstruction of both the objects and the parameters. We consider again the kite-shaped defect given by the parameterization (28) with unknown conductivity σ_i = 7. To initialize the procedure we took σ⁰_i = 3. The topological derivative for this value when Ω_i = ∅ is represented in figure 5(a). As pointed out in [7, 9], it is more efficient to iterate several times with respect to the conductivity before updating the domains (updating the domains requires one to recompute the boundary element method matrices but updating the parameters does not). For this reason, we perform five iterations with the gradient method before updating the domain with the topological derivative method. The final reconstruction of the object is shown in figure 5(b). The values of the conductivity at each iteration are represented in figure 5(c), and the values of the cost functional versus the number of iterations are given in figure 5(d). The combined algorithm stopped after 20 updates of the object, producing a satisfactory reconstruction. Furthermore, the true conductivity σ_i = 7 was approximated by σ_i = 6.51 which is also a reasonable value.

Zoom In Zoom Out Reset image size

Our next two examples deal with complex admittivities, that is, of the form γ = σ + ıωε with ω = 1. We consider again the kite-shaped object defined in (28) with admittivity γ_i = 7 + 2ı. In the exterior medium, γ_e = 1 + ı.

First we assume that both γ_i and γ_e are known. The topological derivative when Ω_i = ∅ in this case is represented in figure 6(a). The iterative algorithm stops at the 21st iteration, yielding the approximate domain that is represented by a dashed line in figure 6(b). The values of the cost functional are given in figure 6(c). We observe that the quality of the reconstruction is similar with real-valued γ (see figure 2(g)) or complex-valued γ (see figure 6(b)).

Zoom In Zoom Out Reset image size

In figure 7, we consider the same problem as in figure 6, but now the interior admittivity γ_i is assumed to be unknown. We started the method with initial guess γ⁰_i = 3 + 4ı and obtained for Ω_i = ∅ the topological derivative shown in figure 7(a). As in the real-valued case, we perform five iterations with the gradient method before updating the domain at each step. The reconstruction of the object at the 17th iteration with respect to the domain is depicted in figure 7(b). The interior admittivity γ_i = 7 + 2ı is approximated by 6.58 + 2.26ı. Both the values of the real part of the admittivity (the conductivity σ) and its imaginary part (the permittivity ε) versus the number of iterations are given in figures 7(c) and (d), respectively. Again, the reconstructions have the same quality as in the real-valued case (compare figure 7(b) with figure 5(b) for the reconstruction of the object and figures 7(c)–(d) with figure 5(c) for the reconstruction of the parameters).

Zoom In Zoom Out Reset image size

Our last two examples are devoted to the reconstruction of smoothly varying conductivities. In this case, we use a finite element method (FEM) with triangular elements and linear interpolation [29] to solve direct and adjoint problems.

We start by illustrating the performance of the gradient method introduced in section 4. Our goal is to recover the conductivity

$$\begin{equation} \sigma (x,y)=1+5\exp (-25((x-0.5)^2/0.8^2+y^2/1.1^2)), \end{equation} \tag{ 29 }$$

that is plotted in figure 8(a). To do that, we start with the value σ⁰ = 1. The gradient method stopped at the 17th iteration. The reconstructed function is shown in figure 8(g). Some of the intermediate approximations are presented in figures 8(b)–(f). We do not recover the function σ(x, y) sharply, but we obtain a reasonable reconstruction. The true function σ is a radially increasing function that attains its maximum at point (0.5,0). Our reconstruction provides an almost radially increasing function attaining its maximum value close to the point (0.55,0). However, the region where σ ≠ 1 is bigger for the reconstructed function than for the true one, and the values of σ are underestimated in the region close to the point (0.55,0) and overestimated far from that region (note that the scale is the same in all the plots in figure 8).

Zoom In Zoom Out Reset image size

In our last example, we consider the problem of recovering an obstacle with unknown smoothly varying conductivity. The obstacle is the ellipse centered at (0.5,0) and semi-axes 0.3 and 0.4, represented in figures 9(a), (b) and (d), are represented by a solid line. The unknown true interior conductivity is given by the function σ_i defined in (29), and the exterior σ_e = 1 is assumed to be known. Note that the function σ(x, y) defined as σ(x, y) = σ_i(x, y) inside the ellipse and σ(x, y) = 1 outside it (represented in figure 9(f)), has a discontinuity at the boundary of the ellipse but is close to the smooth function σ(x, y) considered in the previous example (see figure 8(a)).

Zoom In Zoom Out Reset image size

To recover both, the geometry of the object and the function σ_i, we use the method described in section 3.3. We take as the initial value σ⁰_i = 1.5 and compute the topological derivative when Ω_i = ∅. The largest negative values of the topological derivative are concentrated inside the object, as seen in figure 9(a). Figure 9(b) shows our initial guess Ω₁. We perform now ten iterations with the gradient method to update the interior conductivity, recovering the function represented in figure 9(c). We fix this function and compute the topological derivative again (see figure 9(b)). After 18 iterations with respect to the domain, our reconstruction of the ellipse is represented in figure 9(d). This approximation is quite satisfactory. The reconstructed conductivity is presented in figure 9(e). Our reconstruction is not sharp, but we get reasonable information about the behavior of σ_i, recognizing a radially increasing function. The small discontinuities that are observed inside the ellipse in figure 9(e) correspond to updates in the domains at the iterative procedure (recall that the gradient method is applied to reconstruct σ_i only inside the current approximate domain).

Comparing with our previous example, the combination of topological derivatives with the gradient method improves the reconstruction of the conductivity (compare figures 8(g) and 9(e)) since the center of the ellipse is better reconstructed and the values of σ are closer to the real ones. However, the second strategy uses information on the background, whereas the first one assumes no a priori knowledge of the structure of the medium. It is worth noting that although the overall reconstruction is satisfactory, the maximum error in the region where σ takes the largest values is 2.19 in the first case and 1.64 in the second one. To compare both methods, we have used the data for the conductivity of figure 9 in the gradient method (although it was proposed to recover smooth functions σ). The reconstructed conductivity is represented in figure 10(c). Note that the conductivities in figures 8(g) and 10(c) are almost indistinguishable at first sight. Figure 10 also includes the true conductivity and its reconstruction using the hybrid method for a better comparison.

Zoom In Zoom Out Reset image size

Although more research on the influence of the threshold choices and the way to combine iterations on the quality of the reconstruction remains to be done, the approximations of the variable conductivities could be poor in practice. Nevertheless, they can be of interest to give a good initial guess for a more sophisticated technique. Topological derivatives provided a good reconstruction of the object even when the reconstruction of σ_i was not sharp. The gradient method could be substituted by a more refined technique to recover the conductivity such as those described in [3, 14, 40], which combined with topological derivatives could provide a better algorithm for recovering both objects and their space-dependent conductivities. Further research on the behavior of the schemes in more complex geometries remains to be done.

We have presented a numerical reconstruction technique for electrical impedance tomography problems that allows us to handle both smoothly varying admittivities, and admittivities with sharp discontinuities due to the presence of different objects inside the medium.

The inverse conductivity and impedance problems are the mathematical problems that must be solved for electrical impedance tomography systems to be able to generate images. These systems apply currents to the boundary of a body, measure the induced voltages at this boundary and use this information to reconstruct the electromagnetic properties inside the body.

In real experiments, only the values of the current and the voltages at a 'discrete' set of electrodes are known. Thus, a 'discrete' Neumann condition replaces (2):

$$\begin{eqnarray} \gamma {\partial _\mathbf {n}u } =0& &\qquad \mbox{ between the electrodes,} \end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray} \int _{e_{k}} \gamma {\partial _\mathbf {n}u } \, {\rm d}l = I_{k}, & &\qquad k=1,\ldots ,N, \end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray} &&\sum _{k=1}^N I_{k} = 0, \end{eqnarray} \tag{ 32 }$$

where I_k is the current sent to the kth electrode and e_k denotes the part of ∂Ω that corresponds to the kth electrode.

At the electrodes, we measure the voltages V_k. Dirichlet conditions u = V_k at e_k have been shown not to be realistic [12]. Instead, the following constraint is used:

$$\begin{eqnarray} &&u + z_{k} \gamma {\partial _\mathbf {n}u }= V_{k} \quad \mbox{ \rm at the electrodes $e_{k},$ $k=1,\ldots ,N$,} \\ && \sum _{k=1}^N V_{k}=0, \nonumber \end{eqnarray} \tag{ 33 }$$

where z_k is the effective contact impedance or surface impedance.

Our methods can be extended to handle this more realistic setting. However, the numerical codes are much more complex. Testing our methods for this model will be the subject of forthcoming work.

This research has been supported by MEC grants FIS2011-28 838-C02-02, TRA2010–18 054 and PR2009-0014. AC thanks J B Keller for suggesting the problem and for nice discussions. The authors also thank the referees for useful comments which have helped to improve the paper.

This appendix is devoted to the proof of theorem 3.1. We adapt the strategy used in [6, 7] for inverse scattering problems.

The idea is the following. First, we give a variational formulation of the boundary value problem. Next, we deform the domain Ω_i along the vector field $\bf V$ and compute the shape derivative in the deformed domains. To do so, we differentiate the transformed functionals with respect to a control parameter τ. The resulting expression involves the derivative of the solutions u_τ in the deformed domains. The computation of this derivative is avoided establishing a relationship with the derivatives of modified functionals, in which a free state can be selected to eliminate ${{\rm d} u_{\tau } \over {\rm d} \tau }.$ This adjoint state solves the so-called adjoint problem. Integrating by parts the resulting expression, we find the desired formula for the shape derivative. The topological derivative is then computed using the relationship with shape derivatives established by [18]. The vector field V representing the direction of change of the domain is selected in such a way that ${\bf V}(\mathbf {z})=-\,{\bf n}(\mathbf {z}),$ $\mathbf {z}\in \partial B_\varepsilon (\mathbf {x})$. The normal vector n points inside the ball. This field $\bf V$ is extended to $\mathbb {R}^2$ in such a way that it vanishes away from a narrow neighborhood of ∂B_ε. Selecting the deformations

$$\begin{equation*} \varphi _\tau (\mathbf {x}):=\mathbf {x}+\tau {\bf V}(\mathbf {x}), \quad \mathbf {x}\in \mathbb {R}^2, \end{equation*}$$

it follows that the topological derivative of a shape functional $\mathcal J(R)$ is given by

$$\begin{equation} D_T(\mathbf {x},\mathcal{R})=\lim _{\varepsilon \rightarrow 0}\frac{1}{\mathcal{V}^{\prime }(\varepsilon )}\,\left.\frac{{\rm d}}{{\rm d}\tau }\mathcal{J}(\varphi _\tau (\mathcal{R}\setminus B_\varepsilon ({\bf x})))\right|_{\tau =0}, \qquad \mathbf {x}\in \mathcal{R}, \end{equation} \tag{ A.1 }$$

where $\mathcal{V}^{\prime }(\varepsilon )$ is the derivative of the function $\mathcal{V}(\varepsilon )$. Note that $\varphi _\tau (\mathcal{R}\setminus B_\varepsilon ({\bf x}))= \mathcal{R}\setminus \varphi _\tau ( B_\varepsilon ({\bf x})).$ Let us recall that $\mathcal{V}(\varepsilon )=\pi \varepsilon ^2$.

The proof is organized in seven steps. Steps 1–6 assume Ω_ap = ∅ and $\mathcal{R}=\Omega$. Step 7 indicates the changes when Ω_ap ≠ ∅.

Step 1: variational formulation of the forward problem. We replace (6) by its variational formulation:

$$\begin{equation} u\in H^1_{m}(\Omega ), \quad | \quad \displaystyle b(\Omega _i; u, v)= \ell (v),\quad \forall v\in H^1_{m}(\Omega ), \end{equation} \tag{ A.2 }$$

where

$$\begin{eqnarray*} &&\displaystyle b(\Omega _i;u,v) := \displaystyle \int _{\Omega _e} \gamma _{e} \nabla u\nabla {\overline{v}} \,{\rm d}\mathbf {z}+\int _{\Omega _i} \gamma _{i}\nabla u\nabla {\overline{v}} \,{\rm d}\mathbf {z},\qquad \forall u, v\in H^1_{m}(\Omega ),\\ &&\ell (v) := \displaystyle \int _{\partial \Omega } j \overline{v} \,{\rm d}l, \qquad \forall v\in H^1_{m}(\Omega ), \end{eqnarray*}$$

and H¹_m(Ω) is the subspace of H¹(Ω) formed by functions of zero mean.

Step 2: transformed functionals. By definition, the field $\bf V$ decays fast enough away from ∂Ω_i = Γ_ε to ensure φ_τ(∂Ω) = ∂Ω. Let us set Ω_{i, τ} := φ_τ(Ω_i) and $\Omega _{e,\tau }^{\prime }:=\varphi _\tau (\Omega _e)=\Omega \setminus \overline{\Omega }_{i,\tau }$. The cost functional in the transformed domains is

$$\begin{equation} J(\Omega \setminus \overline{\Omega }_{i,\tau },\gamma _i)= \frac{1}{2}\, \int _{\partial \Omega } |u_\tau -V_{{\rm meas}}|^2 \,{\rm d}l, \end{equation} \tag{ A.3 }$$

where u_τ is a solution of

$$\begin{equation} u_\tau \in H^{1}_{m}(\Omega ), \quad | \quad \displaystyle b(\Omega _{i,\tau }; u_\tau , v)= \ell (v),\qquad \forall v\in H^1_{m}(\Omega ), \end{equation} \tag{ A.4 }$$

with

$$\begin{eqnarray*} b(\Omega _{i,\tau };u,v)&:=&\!\int _{\Omega _{e,\tau }} \!\!\!\gamma _{e} \nabla _{\mathbf {z}_{\tau }} u\nabla _{\mathbf {z}_{\tau }} {\overline{v}} \,{\rm d}\mathbf {z}_{\tau } +\int _{\Omega _{i,\tau }}\!\!\! \gamma _{i}\nabla _{\mathbf {z}_{\tau }} u\nabla _{\mathbf {z}_{\tau }} {\overline{v}} \,{\rm d}\mathbf {z}_{\tau }, \ \ \forall u, v\in H^1_{m}(\Omega ). \end{eqnarray*}$$

For τ = 0, we recover (5) and (A.2).

Step 3: adjoint states. We introduce the modified functional

$$\begin{equation} \fl \mathcal{L}(\Omega _{i,\tau };u,p):=J(\Omega \setminus \overline{\Omega }_{i,\tau },\gamma _i) +\mbox{Re}[b(\Omega _{i,\tau };u,p)-\ell (p)], \;\; \forall u, p\in H^1_{m}(\Omega ). \end{equation} \tag{ A.5 }$$

The derivative of J in the direction $\bf V$ is then

$$\begin{eqnarray} \fl {{\rm d} \over \,{\rm d}\tau } J(\Omega \setminus \overline{\Omega }_{i,\tau },\gamma _i) \bigg |_{\tau =0} &=& \frac{{\rm d}}{{\rm d}\tau }\mathcal{L}(\Omega _{i,\tau };u_\tau ,p) \bigg |_{\tau =0} = \mbox{Re}\left[\frac{{\rm d}}{{\rm d}\tau }b(\Omega _{i,\tau };u,p) \bigg |_{\tau =0} \right] \nonumber \\ [1.5ex] &&+ \mbox{Re} [b(\Omega _i;\dot{u},p) ] + \mbox{Re}\left[\int _{\partial \Omega } (\overline{u -V_{{\rm meas}}}) \dot{u} \, {\rm d}l \right]. \end{eqnarray} \tag{ A.6 }$$

If p satisfies

$$\begin{equation} \fl p\in H^1_{m}(\Omega ), \quad | \quad \displaystyle b(\Omega _i;v,p)= - \int _{\partial \Omega } (\overline{u -V_{{\rm meas}}})\, v \, {\rm d}l,\quad \forall v\in H^1_{m}(\Omega ), \end{equation} \tag{ A.7 }$$

then the derivative is given by the first term in (A.6). The adjoint problem (A.7) is equivalent to (12).

Step 4: derivative of the transformed functionals. Let us compute the first term in (A.6). For any u, p ∈ H¹_m(Ω),

$$\begin{eqnarray} &&\fl \frac{{\rm d}}{{\rm d}\tau }b(\Omega _{i,\tau };u,p)\bigg |_{\tau =0} = \int _{\Omega _e} \gamma _{e} (\nabla u\cdot \nabla \overline{p}) \, \nabla \cdot {\bf V} \,{\rm d}\mathbf {z}-\int _{\Omega _e} \gamma _{e}((\nabla {\bf V}+\nabla {\bf V}^\top ) \nabla u)\cdot \nabla \overline{p} \, {\rm d}\mathbf {z}\, \nonumber \\ &&+ \int _{\Omega _i} \gamma _{i} (\nabla u\cdot \nabla \overline{p}) \,\nabla \cdot {\bf V} \, {\rm d}\mathbf {z}-\int _{\Omega _i} \gamma _{i} ((\nabla {\bf V}+\nabla {\bf V}^\top ) \nabla u)\cdot \nabla \overline{p} \, {\rm d}\mathbf {z}. \end{eqnarray} \tag{ A.8 }$$

Integrating by parts, using the equations satisfied by u and p, and the fact that V vanishes on ∂Ω, the volume integrals on Ω_i and Ω_e and integrals on ∂Ω vanish. Only the contribution to (A.8) from the interface Γ_i = ∂Ω_i remains:

$$\begin{eqnarray} \fl \frac{{\rm d}}{{\rm d}\tau }J(\Omega \setminus \overline{\Omega }_{i,\tau },\gamma _i)\bigg |_{\tau =0} \!\!&=&\mbox{Re}\bigg [\frac{{\rm d}}{{\rm d}\tau }b(\Omega _{i,\tau };u,p)\bigg |_{\tau =0}\bigg ]\nonumber \\ [1.ex] &=&\mbox{Re} \bigg [\int _{\Gamma } \!\!\bigg ( \gamma _{i} \left({\gamma _{i}\over \gamma _{e}}-1\right) \partial _{\mathbf {n}} u^- \partial _{\mathbf {n}} \overline{p}^-\! + (\gamma _{i}-\gamma _{e}) \partial _{\bf t} u^- \partial _{\bf t} \overline{p}^- \! \bigg ) \,\! {\rm d}l \bigg ] \end{eqnarray} \tag{ A.9 }$$

Step 5: topological derivative when γ_i and γ_e are constant and Ω_ap = ∅. For any subset B⊂Ω, we set $\mathcal{J}(\Omega \setminus \overline{B})= J(\Omega \setminus \overline{B},\gamma _{i})$. We use identities (A.1) and (A.9) with $\Omega _i=B_\varepsilon (\mathbf {x})$ to obtain

$$\begin{eqnarray} &&\fl D_T(\mathbf {x},\Omega ,\gamma _i)=\lim _{\varepsilon \rightarrow 0}\frac{1}{\mathcal{V}^{\prime }(\varepsilon )} \,\mbox{Re} \bigg [\int _{\Gamma _\varepsilon } \bigg ( \gamma _{i} \left({\gamma _{i}\over \gamma _{e}}-1\right) \partial _{\mathbf {n}} u_{\varepsilon }^- \partial _{\mathbf {n}} \overline{p}_{\varepsilon }^- + (\gamma _{i}-\gamma _{e}) \partial _{\bf t} u_{\varepsilon }^- \partial _{\bf t} \overline{p}_{\varepsilon }^- \bigg ) \, {\rm d}l \bigg ],\nonumber\\ &&{\bf x}\in \Omega . \end{eqnarray} \tag{ A.10 }$$

Here, u_ε and p_ε solve the forward (6) and adjoint (12) problems when $\Omega _{i}=B_\varepsilon (\mathbf {x})$ and $\Gamma = \Gamma _\varepsilon =\partial B_\varepsilon (\mathbf {x})$. The asymptotic behavior of u_ε and p_ε at Γ_ε is obtained expressing them as corrections of u and p:

$$\begin{eqnarray*} &&u_\varepsilon (\mathbf {z}) = u(\mathbf {z}) \chi _{\Omega \setminus \overline{B}_{\varepsilon } }(\mathbf {z}) + v_{\varepsilon } (\mathbf {z}), \qquad p_\varepsilon (\mathbf {z}) = p(\mathbf {z}) \chi _{\Omega \setminus \overline{B}_{\varepsilon } }(\mathbf {z}) + q_{\varepsilon } (\mathbf {z}), \end{eqnarray*}$$

and expanding the remainders in powers of ε. Let us denote by B the unit ball and set $\Omega _{\varepsilon } := (\Omega -\mathbf {x})/ \varepsilon .$ Making the change of variables ${\bf \xi} :=(\mathbf {z}-\mathbf {x})/\varepsilon$, the correction $v_{\varepsilon } ({\bf \xi})$ satisfies

$$\begin{equation} \fl \left\lbrace \begin{array}{@{}ll}\Delta _{{\bf \xi}} v_{\varepsilon } =0& \ \mbox{in $\Omega _{\varepsilon } \setminus \overline{B}$}, \\ [1.5ex] \Delta _{{\bf \xi}} v_{\varepsilon } =0& \ \mbox{in $B$}, \\ [1.5ex] v_{\varepsilon }^- -v_{\varepsilon }^+ = u(\mathbf {z}) = u(\mathbf {x})+ \varepsilon {\bf \xi} \cdot \nabla u(\mathbf {x}) +O(\varepsilon ^2) &\ \mbox{on $\Gamma $}, \\ [1.5ex] \gamma _{i} \mathbf {n}({\bf \xi}) \nabla _{{\bf \xi}} \cdot v_{\varepsilon }^- - \gamma _{e} \mathbf {n}({\bf \xi})\cdot \nabla _{{\bf \xi}} v_{\varepsilon }^+ = \varepsilon \gamma _{e} \mathbf {n}({\bf \xi}) \cdot \nabla u(\mathbf {x}) + O(\varepsilon ^{2}) & \ \mbox{on $\Gamma $}, \\ [1.5ex] \mathbf {n}({\bf \xi})\cdot \nabla _{{\bf \xi}} v_{\varepsilon } = 0 & \ \mbox{on $\partial \Omega _{\varepsilon }$}. \end{array} \right. \end{equation} \tag{ A.11 }$$

Expanding $v_{\varepsilon }({\bf \xi})$ in powers of ε $v_{\varepsilon }({\bf \xi})= v^{(1)}({\bf \xi}) + \varepsilon v^{(2)}({\bf \xi})+ O(\varepsilon ^{2})$, we find $v^{(1)}({\bf \xi}) = u(\mathbf {x}) \chi _{B}({\bf \xi})$ and

$$\begin{equation*} v^{(2)}({\bf \xi})\sim \nabla u(\mathbf {x}) \bigg ({2\gamma _e\over \gamma _e+ \gamma _{i} } {\bf \xi}\chi _B({\bf \xi})+{\gamma _{e}-\gamma _{i} \over \gamma _{e} + \gamma _{i}} {{\bf \xi}\over |{\bf \xi}|^2 } \chi _{\mathbb {R}^{2}\setminus \overline{B}}({\bf \xi})\bigg ) = \nabla u(\mathbf {x}) \cdot g({\bf \xi}), \end{equation*}$$

up to an error term which is O(ε²). Thus,

$$\begin{eqnarray*} u_\varepsilon (\mathbf {z}) &=&u(\mathbf {z}) \chi _{\Omega \setminus \overline{B}_{\varepsilon } }(\mathbf {z}) + u(\mathbf {x}) \chi _{ B_{\varepsilon } }(\mathbf {z}) + \varepsilon \nabla u(\mathbf {x}) \cdot g({\bf \xi}) + O(\varepsilon ^{2}), \quad \mathbf {z}\in \Gamma _{\varepsilon }. \end{eqnarray*}$$

A similar expansion holds for $p_\varepsilon (\mathbf {z})$. This implies

$$\begin{equation*} u_{\varepsilon }(\mathbf {z}) \rightarrow u(\mathbf {x}), \qquad p_{\varepsilon }(\mathbf {z}) \rightarrow p(\mathbf {x}), \qquad \mbox{as $\varepsilon \rightarrow 0$}, \end{equation*}$$

$$\begin{equation*} \nabla u_{\varepsilon }^-(\mathbf {z}) \rightarrow {2\gamma _e\over \gamma _e+ \gamma _{i} } \nabla u(\mathbf {x}), \qquad \nabla p_{\varepsilon }^-(\mathbf {z}) \rightarrow {2\gamma _e\over \gamma _e+ \gamma _{i}} \nabla p(\mathbf {x}),\qquad \mbox{as $\varepsilon \rightarrow 0$}, \end{equation*}$$

uniformly when $|\mathbf {z}-\mathbf {x}|=\varepsilon .$

Taking limits in the two integrals appearing in (A.10) and taking into account that $\mathcal{V}^{\prime }(\varepsilon )=2\pi \varepsilon$, (10) follows.

Step 6: variable parameters. The previous derivation holds when the parameters are constant. When they vary smoothly, we expand them about x and repeat the previous computations. Note that due to the choice of V, only the values of γ_i and γ_e in a narrow band about Γ_ε are relevant for the computation of the shape derivative. For the topological derivative, we set $\gamma _{i}(\mathbf {z})= \gamma _{i}(\mathbf {x})+ \varepsilon \, {\bf \xi}\cdot \nabla \gamma _{i}(\mathbf {x}) + O(\varepsilon ^{2})$ and an analogous expression for γ_e. Then, we check that v⁽¹⁾ and v⁽²⁾ remain unchanged. For piecewise smooth γ_e and γ_i, we use a regularizing approximation; see [7] for details.

Step 7: topological derivative when Ω_ap ≠ ∅. The formula for $\mathbf {x}\in \Omega \setminus \overline{\Omega }_{ap}$ (i.e. (10)) follows in a similar way, replacing Ω_i with $\Omega _{ap}\cup B_{\varepsilon }(\mathbf {x})$ in step 5. The formula for $\mathbf {x}\in \Omega _{ap}$ (i.e. (11)) is obtained following the same steps, taking into account that now $B_\varepsilon (\mathbf {x})$ is part of the exterior domain and that the roles of γ_i and γ_e are interchanged. That is, the interior problem in B has parameter γ_e, whereas γ_i corresponds to the exterior problem. See [7] for further details.