3D Reconstruction for Partial Data Electrical Impedance Tomography Using a Sparsity Prior

In electrical impedance tomography the electrical conductivity inside a physical body is computed from electro-static boundary measurements. The focus of this paper is to extend recent result for the 2D problem to 3D. Prior information about the sparsity and spatial distribution of the conductivity is used to improve reconstructions for the partial data problem with Cauchy data measured only on a subset of the boundary. A sparsity prior is enforced using the $\ell_1$ norm in the penalty term of a Tikhonov functional, and spatial prior information is incorporated by applying a spatially distributed regularization parameter. The optimization problem is solved numerically using a generalized conditional gradient method with soft thresholding. Numerical examples show the effectiveness of the suggested method even for the partial data problem with measurements affected by noise.


1.
Introduction. Sparse reconstruction for electrical impedance tomography (EIT) with full boundary data has been utilized in [9,14,15] and are based on algorithms from [3,4]. A similar approach was used for the 2D partial data problem in [8] by applying a spatially varying regularization parameter; this paper extends the algorithm to the 3D partial data problem. The main contributions are in deriving the Fréchet derivative for the algorithm and in the numerical results in 3D.
The inverse problem in EIT consists of reconstructing an electrical conductivity distribution in the interior of an object from electro-static boundary measurements on the surface of the object. The underlying mathematical problem is known as the Calderón problem in recognition of Calderón's seminal paper [6]. While the Calderón problem can also be considered in two dimensions, physical electric fields are intrinsically three dimensional, and thus the reconstruction problem in EIT should ideally use a 3D reconstruction algorithm to reduce modelling errors in the reconstruction.
Consider a bounded domain Ω ⊂ R 3 with smooth boundary ∂Ω. In order to consider partial boundary measurements we introduce the subsets Γ N , Γ D ⊆ ∂Ω for the Neumann and Dirichlet data respectively. Let σ ∈ L ∞ (Ω) with 0 < c ≤ σ a.e. denote the conductivity distribution in Ω. Applying a boundary current flux g (Neumann condition) through Γ N ⊆ ∂Ω gives rise to the interior electric potential u characterized as the solution to ∇ · (σ∇u) = 0 in Ω, σ ∂u ∂ν = g on ∂Ω, where ν is an outward unit normal to ∂Ω. The latter condition in (1) is a grounding of the total electric potential along the subset Γ D ⊆ ∂Ω. To be precise we define the spaces L 2 (∂Ω) ≡ g ∈ L 2 (∂Ω) ∂Ω g ds = 0 , consisting of boundary functions with mean zero (here ·, · denotes the dual pairing), and the spaces consisting of functions with mean zero on Γ D . Using standard elliptic theory it follows that (1) has a unique solution u ∈ H 1 Γ D (Ω) for any g ∈ H −1/2 (∂Ω). This defines the Neumann- Recently the partial data Calderón problem has been studied intensively. In 3D uniqueness has been proved under certain conditions on Γ D and Γ N [5,13,16,18]. Also stability estimates of log-log type have been obtained for the partial data problem [12]; this suggests that the partial data problem is even more ill-posed and hence requires more regularization than the full data problem which has log type estimates [2].
The data considered here consist of K pairs of Cauchy data taken on the subsets Γ D and Γ N , i.e.
We assume that the unknown conductivity is given as σ = σ 0 + δσ, where σ 0 is a known background conductivity. For some fixed c ∈ (0, 1) and σ 0 ∈ H 1 (Ω) where c ≤ σ 0 ≤ c −1 , define the closed and convex subset Similarly define The inverse problem is then to approximate δσ ∈ A 0 given the data (2). Let {ψ j } ∞ j=1 denote a chosen orthonormal basis for H 1 0 (Ω). For sparsity regularization we approximate δσ by argmin δγ∈A0 Ψ(δγ) using the following Tikhonov functional with the discrepancy terms J k and penalty term P given by for c j ≡ δγ, ψ j H 1 (Ω) . The regularization parameter α j for the sparsity-promoting 1 penalty term P is distributed such that each basis coefficient can be regularized differently; we will return to this in Section 3. It should be noted how easy and natural the use of partial data is introduced in this way, simply by only minimizing the discrepancy on Γ D where the Dirichlet data is known and ignoring the rest of the boundary.
Remark 1. The non-linearity of σ → R σ leads to a non-convex discrepancy term, i.e. Ψ is non-convex. When applying a gradient based optimization method, the best we can hope is to find a local minimum.
This paper is organised as follows: in Section 2 we derive the Fréchet derivative of J k and reformulate the optimization problem using the generalized conditional gradient method as a sequence of linearized optimization problems. In Section 3 we explain the idea of the spatially dependent regularization parameter designed for the use of prior information. Finally, in Section 4 we show the feasibility of the algorithm by numerical examples.
2. Sparse Reconstruction. In this section the sparse reconstruction of δσ based on the optimization problem (4) is investigated for a bounded domain Ω ⊂ R 3 with smooth boundary. The penalty term emphasizes that δσ should only be expanded by few basis functions in the given orthonormal basis. The partial data problem comes into play in the discrepancy term, in which we only fit the data on part of the boundary. Ultimately, this leads to Algorithm 1 at the end of this section.
For fixed g let u be the unique solution to (1). Define the solution operator F g : σ → u and further its trace F g : σ → u| ∂Ω (note that R σ g = F g (σ)). In order to compute the derivative of F g , let γ ∈ A and g ∈ L p (∂Ω) ∩ H −1/2 (∂Ω) for p ≥ 8 5 . Then following the proofs of Theorem 2.2 and Corollary 2.1 in [15] whilst applying the partial boundary Γ D we have The linear map (F g ) γ maps η to w| ∂Ω , where w is the unique solution to − ∇ · (γ∇w) = ∇ · (η∇F g (γ)) in Ω, γ ∂w ∂ν = 0 on ∂Ω, Note that (F g ) γ resembles a Fréchet derivative of F g evaluated at γ due to (5), however A is not a linear vector space, thus the requirement γ, γ + η ∈ A. The first step in minimizing Ψ using a gradient descent type iterative algorithm is to determine a derivative to the discrepancy terms J k . For this purpose the following proposition is applied, and is a special case of [15,Theorem 3.1].
, there is the following estimate with C only depending on c, Ω and q: Now we can formulate the Fréchet derivative of J k . 8 5 , and χ Γ D be a characteristic function on Γ D . Then there exists c ∈ (0, 1) as the bound in A 0 sufficiently close to 1, such and the Fréchet derivative (J k ) δγ of J k on H 1 0 (Ω) evaluated at δγ is given by Proof. For the proof the index k is suppressed. First it is proved that E ∈ L 6/5 (Ω). Write forq ∈ (2, Q(c)) ∩ [ 3 2 , 12 5 ]. Choosing c sufficiently close to 1 leads to Q(c) > 12 5 . By (10) and (11) then |∇F h (γ)|, |∇F g (γ)| ∈ L 12/5 (Ω), and Hölder's generalized inequality entails that E ∈ L r (Ω) with 1 r = 5 12 + 5 12 , i.e. r = 6 5 , The Sobolev embedding theorem [1] implies the embedding H 1 (Ω) → L 6 (Ω) as Ω ⊂ R 3 . Thus E ∈ L 6/5 (Ω) = (L 6 (Ω)) ⊂ (H 1 (Ω)) ⊂ (H 1 0 (Ω)) = H −1 (Ω). Next we prove (9). J δγ η is by the chain rule (utilizing that R γ g = F g (γ)) given as where χ Γ D is enforcing that the integral is over Γ D . The weak formulations of (1), with Neumann data χ Γ D (R γ g − f ), and (6) are Now by letting v ≡ w in (13) and (14), we obtain using the definition w| ∂Ω = (F g ) γ η that We seek to find a direction η for which the discrepancy decreases. As J δγ ∈ H −1 (Ω) it is known from Riesz' representation theorem that there exists a unique function in H 1 0 (Ω), denoted by G(δγ), such that Now η ≡ −G(δγ) points in the direction of steepest descend among the viable directions. Furthermore, since G(δγ)| ∂Ω = 0 the boundary condition δσ| ∂Ω = 0 for the approximation will automatically be fulfilled. Note that G(δγ) is the unique solution to for which (15) is the weak formulation. In each iteration step we need to determine a step size s i for an algorithm resembling a steepest descent δγ i+1 = δγ i − s i G(δγ i ). As in [8] a Barzilai-Borwein step size rule is applied A maximum step size s max is enforced to avoid problems in the situation where With inspiration from [21], s i will be initialized by (16), after which it is thresholded to lie in [s min , s max ] for two chosen positive constants s min and s max . It is noted in [21] that Barzilai-Borwein type step rules lead to faster convergence if we do not restrict Ψ to decrease in every iteration. Therefore, one makes sure that the following so-called weak monotonicity is satisfied, which compares Ψ(δγ i+1 ) with the most recent M steps. Let τ ∈ (0, 1) and M ∈ N, then s i is said to satisfy the weak monotonicity with respect to M and τ if the following is satisfied If (17) is not satisfied, the step size s i is reduced until this is the case.
Here {ψ j } ∞ j=1 is an orthonormal basis for H 1 0 (Ω) in the H 1 -metric, and P A0 is a projection of H 1 0 (Ω) onto A 0 to ensure that (1) is solvable (note that H 1 0 (Ω) does not embed into L ∞ (Ω), i.e. ζ i+1 may be unbounded). By use of the map S β : R → R defined below, known as the soft shrinkage/thresholding map with threshold β > 0, the solution to (18) is easy to find directly (see also [7, Section 1.5]) where d j ≡ δγ i − s i G(δγ i ), ψ j H 1 (Ω) are the basis coefficients for δγ i − s i G(δγ i ). The projection P A0 : where T c is the following truncation that depends on the constant c ∈ (0, 1) in (3) where v < c a.e., c −1 where v > c −1 a.e., v else.
Since σ 0 ∈ H 1 (Ω) and c ≤ σ 0 ≤ c −1 , it follows directly from [20, Lemma 1.2] that T c and P A0 are well-defined, and it is easy to see that P A0 is a projection. It should also be noted that 0 ∈ A 0 since c ≤ σ 0 ≤ c −1 , thus we may choose δγ 0 ≡ 0 as the initial guess in the algorithm, which is appropriate as we expect the solution to be sparse. The algorithm is summarized in Algorithm 1. In the numerical experiments in Section 4 the stopping criterion is when the step size s i gets below a threshold s stop .

Algorithm 1 Sparse Reconstruction for Partial Data EIT
Set δγ 0 := 0. While stopping criteria not reached . Compute step length s i by (16), and decrease it till (17) is satisfied.
3. Prior Information. Prior information is intrinsically linked to the penalty term P for Tikhonov-like functionals, and the regularization parameter determines how much this prior information is enforced. In the case of sparsity regularization this implies knowledge of how sparse we expect the solution is in general. Instead of applying the same prior information for each basis function, a distributed parameter is applied. Let where α is a usual regularization parameter, corresponding to the case where no prior information is considered about specific basis functions. The µ j ∈ (0, 1] will be used to weight the penalty depending on whether a specific basis function should be included in the expansion of δσ. The µ j are chosen as µ j = 1, no prior on c j , ∼ 0, prior that c j = 0, i.e. if we know that a coefficient in the expansion of δσ should be non-zero, we can choose to penalize that coefficient less.
3.1. Applying the FEM Basis. In order to improve the sparsity solution for finding small inclusions, it seems appropriate to include prior information about the support of the inclusions. There are different methods available for obtaining such information assuming piecewise constant conductivity [11,17] or real analytic conductivity [10]. The idea is to be able to apply such information in the sparsity algorithm in order to get good contrast in the reconstruction while maintaining the correct support, even for the partial data problem.
Suppose that as a basis we consider a finite element method (FEM) basis {ψ j } N j=1 for the subspace V h ⊆ H 1 0 (Ω) of piecewise affine functions on each element. Let δγ ∈ V h with mesh nodes {x j } N j=1 , then δγ(x) = N j=1 δγ(x j )ψ j (x) and ψ j (x k ) = δ j,k , i.e. for each node there is a basis function for which the coefficient contains local information about the expanded function; this is convenient when applying prior information about the support of an inclusion.
When applying the FEM basis for mesh nodes {x j } N j=1 , the corresponding functional is It is evident that the penalty corresponds to determining inclusions with small support, and prior information on the sparsity corresponds to prior information on the support of δσ.
We cannot directly utilize (20) due to the FEM basis not being an orthonormal basis for H 1 0 (Ω), and instead we suggest the following iteration step as in [8]: Note that the regularization parameter will depend quite heavily on the discretization of the mesh, i.e. for the same domain a good regularization parameter α will be much larger on a coarse mesh than on a fine mesh. Instead we can weight the regularization parameter according to the mesh cells, by having α j ≡ αβ j µ j . This leads to a discretization of a weighted L 1 -norm penalty term: where f µ : Ω → (0, 1] is continuous and f µ (x j ) = µ j . The weights β j consists of the node volume computed in 3D as 1/4 of the volume of supp(ψ j ) (if using a mesh of tetrahedrons). This corresponds to splitting each cell's volume evenly amongst the nodes, and it will not lead to instability on a regular mesh. This will make the choice of α almost independent of the mesh, and will be used in the numerical examples in the following section.
Remark 2. The corresponding algorithm with the FEM basis is the same as Algorithm 1, except that the update is applied via (21).

Numerical Examples.
In this section we illustrate, through a few examples, the numerical algorithm implemented by use of the finite element library FEniCS [19]. First we consider the full data case Γ D = Γ N = ∂Ω both without and with prior information, and then we do the same for the partial data case. For the following examples Ω is the unit ball in R 3 . The numerical phantom consists of a background conductivity with value 1, a smaller ball inclusion with value 2 centred at (−0.09, −0.55, 0) and with radius 0.35, and two large ellipsoid inclusions with value 0.5. One ellipsoid is centred at (−0.55 sin( 5 12 π), 0.55 cos( 5 12 π), 0) and with semi-axes of length (0.6, 0.3, 0.3). The other ellipsoid is centred at (0.45 sin( 5 12 π), 0.45 cos( 5 12 π), 0) and with semi-axes of length (0.7, 0.35, 0.35). The two ellipsoids are rotated respectively 5 12 π and − 5 12 π about the axis parallel to the Z-axis and through the centre of the ellipsoids; see Figure 1. In this paper we do not consider choice rules for α; it is chosen manually by trial and error. The parameters are chosen as σ 0 ≡ 1, M = 5, τ = 10 −5 , s min = 1, s max = 1000, and the stopping criteria is when the step size is reduced below s stop = 10 −3 . Let Y m n denote Laplace's spherical harmonics of degree n and order m, with real form The Neumann data consists ofỸ m n for −n ≤ m ≤ n and n = 1, 2, . . . , 5, i.e. a total of K = 35 current patterns. For the partial data examples a half-sphere is used for local data Γ = Γ N = Γ D , and the corresponding Neumann data are scaled to have the same number of periods as the full data examples.
When applying prior information, the coefficients µ j are chosen as 10 −2 where the support of δσ is assumed, and 1 elsewhere. The assumed support is a 10% dilation of the true support, to show that this inaccuracy in the prior information still leads to improved reconstructions.
For the simulated Dirichlet data, the forward problem is solved on a very fine mesh, and afterwards interpolated onto a different much coarser mesh in order to avoid inverse crimes. White Gaussian noise has been added to the Dirichlet data {f k } K k=1 on the discrete nodes on the boundary of the mesh. The standard deviation of the noise is chosen as max k max xj ∈Γ D |f k (x j )| as in [8], where = 10 −2 corresponding to 1% noise.  Figure 2 shows 2D slices of reconstructions from full boundary data. It is seen that the reconstructions attain the correct contrast, and close to the boundary gives good approximations to the correct support for the inclusions. Using the overestimated support as prior information gives vastly improved reconstruction further away from the boundary. This holds for the entire 3D reconstruction as seen in the bottom part of Figure 2, and makes it possible to get a reasonable separation of the inclusions.
From Figure 3 2D slices of partial data reconstructions are shown, and it is evident that far from the measured boundary the reconstructions suffer severely. Reconstructing with data on the lower part of the sphere gives a reasonable reconstruction with correct contrast for the ball inclusion, however the larger inclusions are hardly reconstructed at all.
With data on the top half of the sphere yields a reconstruction with no clear separation of the ellipsoid inclusions, which is much improved by use of the overestimated support. There is however an artefact in one of the reconstructed inclusions that could correspond to data from the ball inclusion, which is not detected in the reconstruction even when the additional prior information is used.
The reconstructions shown here are consistent with what was observed in [8] for the 2D problem, and it is possible to reconstruct the correct contrast even in the partial data case, and also get decent local reconstruction close to the measured boundary. However, the partial data reconstructions seems to be slightly worse in 3D when no prior information about the support is applied.