ARTIFICIAL BOUNDARY CONDITIONS AND DOMAIN TRUNCATION IN ELECTRICAL IMPEDANCE TOMOGRAPHY. PART THEORY AND PRELIMINARY RESULTS

. Artiﬁcial boundary conditions have long been an active research topic in numerical approximation of scattering waves: The truncation of the computational domain and the assignment of the conditions along the ﬁctitious boundary must be done so that no spurious reﬂections occur. In inverse boundary value problems, a similar problem appears when the estimation of the unknowns is restricted to a domain that represents the whole domain of the solutions of a partial diﬀerential equation with unknown coeﬃcient. This 2010 Mathematics Subject Classiﬁcation. Primary: 35J20, 35R30, 65M32, 65M55. problem is signiﬁcantly more challenging than general scattering problems, be- cause the coeﬃcients representing the unknown material parameter of interest are not known in the truncated portion and assigning suitable condition on the ﬁctitious boundary is part of the problem also. The problem is addressed by deﬁning a Dirichlet-to-Neumann map, or Steklov-Poincar´e map, on the boundary of the domain truncation. In this paper we describe the procedure, provide a theoretical justiﬁcation and illustrate with computed examples the limitations of imposing ﬁxed boundary condition. Extensions of the proposed approach will be presented in a sequel article.

problem is significantly more challenging than general scattering problems, because the coefficients representing the unknown material parameter of interest are not known in the truncated portion and assigning suitable condition on the fictitious boundary is part of the problem also. The problem is addressed by defining a Dirichlet-to-Neumann map, or Steklov-Poincaré map, on the boundary of the domain truncation. In this paper we describe the procedure, provide a theoretical justification and illustrate with computed examples the limitations of imposing fixed boundary condition. Extensions of the proposed approach will be presented in a sequel article.
1. Introduction. Electrical impedance tomography (EIT) is an imaging modality in which the unknown electric conductivity distribution or, more generally, admittivity distribution in an inaccessible region inside a body is estimated based on electric current and voltage measurements on the surface of the body. In the idealized formulation of the problem, the data consist of the complete set of Cauchy data at the boundary, which is tantamount to the knowledge of either the Dirichlet-to-Neumann map, or its inverse, over the boundary. A more realistic model assumes that the current is injected through a finite set of contact electrodes while outside the electrodes the current density vanishes. The corresponding voltages at the same electrodes are measured, and the data consist of a complete set of pairs of injected currents and resulting voltages.
It is not uncommon that the domain in which the conductivity is to be estimated is only a portion of the conducting body. Thus to improve the computational efficiency of the procedure, it would be natural to truncate the computational domain so as to comprise the part of interest only. The domain truncation introduces an artificial boundary, which we will refer to as the truncation surface in the following, along which suitable conditions need to be assigned. These boundary conditions are usually unknown, in particular when the conductivity distribution in the truncated portion of the domain is not available. In this article, we propose a partial Dirichlet-to-Neumann operator along the truncation surface as an appropriate boundary condition. This operator may be known, approximately known or unknown: in the latter cases its estimation becomes part of the EIT inverse problem itself.
To put the proposed approach in the context of the existing literature, we recall two important areas where the local Dirichlet-to-Neumann maps have been suggested as an artificial boundary conditions. In scattering theory, artificial domain truncation has obtained considerable attention, and different non-reflecting boundary conditions at the fictitious outer surface of the computational domain have been proposed. In [5,6,8,9,14] the authors proposed the Dirichlet-to-Neumann map as a perfect non-reflecting boundary condition and showed the approach to be computationally viable.
Unlike in scattering problems, where the outside domain is usually known to be free space, and the Dirichlet-to-Neumann map for outgoing solutions is readily available, in domain decomposition methods the computational domain is partitioned in subdomains that need to be connected so as not to create spurious reflections from the fictitious boundaries. It is in that context that the idea of using the Dirichlet-to-Neumann map, which in the domain decomposition community is usually referred to as Steklov-Poincaré map, was advocated by numerous authors; see [19] and references therein.
The use of the Dirichlet-to-Neumann map advocated in the present article, is indeed closely connected to both of the aforementioned applications. The use of non-local boundary condition based on a partial Dirichlet-to-Neumann map on truncation surfaces proposed in [10,11], is closely related to the proposed approach. In the present paper, the domain truncation problem is carefully analyzed for the EIT problem based on the complete electrode model for a few different geometries, within the Bayesian framework for inverse problems, and the results are validated numerically with phantom data.
2. Domain truncation. In this section we briefly review the forward model of the EIT problem described by the complete electrode model. A general reference for regularity properties of elliptic boundary value problems is [7].
For the sake of definiteness, let Ω ⊂ R n , n = 2, 3, be a bounded, connected domain, whose boundary is a piecewise smooth curvilinear polygon. For consitency regarding the numerical implementations, we assume that σ : Ω → R is a piecewise continuous function, 0 < c ≤ σ(x) ≤ C for some constants c and C; the theoretical considerations can be extended to more general L ∞ conductivities. The contact electrodes, attached on the boundary ∂Ω, are modeled by disjoint curve segments (n = 2) or closed surface patches (n = 3) e ⊂ ∂Ω, 1 ≤ ≤ L. On each electrode, we define a contact impedance z > 0. The complete electrode model (CEM) forward problem can be stated formally as follows.
Problem. Given the conductivity σ and the contact impedances z , find the voltage potential u ∈ H 1 (Ω) and electrode voltages U ∈ R, 1 ≤ ≤ L satisfying in Ω, where the applied electric currents satisfy the conservation of charge condition The CEM boundary value problem is well posed and it has a unique solution u, {U } L =1 ∈ H 1 (Ω) × R L , provided that a ground condition is given, e.g., by setting [20] for details. The proof is based on a coercivity argument.
Let S be a smooth hypersurface S intersecting Ω, and for simplicity, let Γ = S ∩Ω consist of a single component so that S splits Ω into two disjoint subdomains, denoted by Ω 1 and Ω 2 . We shall denote by S j = Ω j ∩ ∂Ω, j = 1, 2, the two disjoint components of the boundary ∂Ω. Suppose that all the electrodes are on the surface of the subdomain Ω 1 , that is, e ⊂ S 1 for = 1, 2, . . . , L, as illustrated schematically in Figure 1.
Introducing the notations σ j = σ Ωj and u j = u Ωj , j = 1, 2, where u is the solution of Problem 2 and denoting by ν j , j = 1, 2, the normal vectors along Γ pointing out of the domain Ω j (see Figure 1), it readily follows that Assuming, for the time being, that the boundary value u Γ = φ ∈ H 1/2 (Γ) is known, we arrive at the following formulation. Find a solution u 1 satisfying where the Neumann data on the interface Γ is determined by solving the mixed boundary value problem in Ω 2 , In order to obtain a solution of the original problem, u 1 must satisfy the continuity condition u 1 Γ = φ. To guarantee that this is the case, we define the local Dirichletto-Neumann map, or Steklov-Poincaré map of the system (6), Λ : where H 1/2 (Γ) * is the dual of H 1/2 (Γ), which can be identified with H −1/2 (Γ) the Sobolev space of distributions on the hypersurface S that are supported in Γ [21]. If the mapping Λ is given, we can replace (5) and (6), with the single autonomous system governing the truncated forward problem can be defined as follows.
Problem. Given the conductivity σ 1 , the contact impedances z , and the mapping Λ on Γ, find the voltage potential u 1 ∈ H 1 (Ω 1 ) and electrode voltages U ∈ R, 1 ≤ ≤ L which satisfies where the applied currents I satisfy the conservation of charge condition (2).
In the following subsection, we prove that this problem is well posed, has a unique solution, and that the solution indeed coincides with the restriction of the solution of Problem 2 restricted to the subdomain Ω 1 .

Variational formulation.
To show the unique solvability of Problem 2, we consider first the mixed boundary value problem (6). Multiplying the differential equation by a test function v and integrating by parts we arrive at the weak form which we write as Let Φ ∈ H 1 (Ω 2 ) be an arbitrary but fixed function such that Φ| Γ = φ, and denote by H 1 (Ω 2 , Γ 0 ) the closure of the set {u ∈ C 2 (Ω 2 ) ∩ C 1 (Ω 2 ) | u Γ = 0} with respect to the Sobolev norm. Then the weak form of the problem (6) is: The existence and uniqueness of the solution u with fixed Φ follows by the Lax-Milgram lemma from the coercivity of B 2 thanks to the Poincaré inequality. The identity (9) allows us to define the Steklov-Poincaré map (7). Similarly, we can write the weak form of the boundary value problem (8). Following the derivation in [20], multiplication by a test function w integration by parts leads to the identity On the other hand, Therefore, equations (10) and (11) together imply that which we write as where U ∈ R L and V ∈ R L are vectors with components U and V , respectively. We prove now the following result.
Theorem 2.1. There exists a unique pair (u 1 , U ) ∈ H 1 (Ω 1 ) × R L satisfying the equation (12) for all (w, V ) ∈ H 1 (Ω 1 )×R L and the ground condition (3). Moreover, is the solution of the weak form corresponding to Problem 2.
Proof. As in [20], we define the space It was shown in [20] that the quadratic form is coercive. Observe that B 1 corresponds to a weak form of the boundary value problem with vanishing Neumann data along the interface Γ. We prove the existence of the solution of the weak form equation (12), showing that the quadratic form is bounded, symmetric and positive semi-definite. The boundedness follows from the mapping properties (7) of the operator Λ. To show the positive definiteness, assume that w ∈ H 1 (Ω 1 ), and let u 2 ∈ H 1 (Ω 2 ) be the weak solution of the mixed boundary value problem (6) with the Dirichlet data Then, as pointed out above, u 2 satisfies the identity (9) or, equivalently, using the definition of the local Dirichlet-to-Neumann map, Observe that if w = 0 in the quotient space H, then u 2 = constant, and B 2 (u 2 , u 2 ) = 0. Therefore the bilinear form (12) is well defined on the quotient space H, that is, its value does not depend on the choice of a representative in the equivalent class.
To show the symmetry, let w and w be two functions in H 1 (Ω 1 ), and u 2 and u 2 the corresponding functions in H 1 (Ω 2 ). It follows from the symmetry of the quadratic form B 2 that which is the desired property. The Lax-Milgram theorem assures the existence of the solution in H (Ω 1 ), unique up to an additive constant determined by the ground condition (3).
It remains to show that u 1 is the restriction to Ω 1 of the weak solution u in Ω. Given that u 1 is a solution of (12), let u 2 ∈ H 1 (Ω 2 ) be defined as above. For each test function v ∈ H 1 (Ω), define v j = v Ωj ∈ H 1 (Ω j ), j = 1, 2. Letting w = u 1 and combining (12) with w = v 1 and (13), it follows that This implies that the function u defined piecewise via its restrictions u Ωj = u j , j = 1, 2, is indeed in H 1 (Ω) and (u, U ) satisfies the weak form of the complete electrode model, for all test functions (v, V ) ∈ H (Ω), completing the proof.
3. Inverse problem. We are now ready to state two versions of the EIT inverse problem with domain truncation, both of which will be studied in the light of computed examples. In the following, we refer to a set of L − 1 linearly independent current patterns satisfying the conservation of charge condition (2) as a current frame, and we denote it by The corresponding set of voltages which satisfy the ground condition (3) is denoted by Problem.
Let Ω be a bounded domain divided into two non-overlapping subdomains Ω 1 and Ω 2 , with known contact impedances. We consider the following pair of inverse problems: Given the current-voltage data (I , U ), (a) estimate the conductivity σ 1 = σ Ω1 , assuming that an approximation Λ 0 of the local Dirichlet-to-Neumann operator Λ is given; (b) estimate both the conductivity σ 1 = σ Ω1 and the local Dirichlet-to-Neumann operator Λ from the data. 3.1.1. A rectangular geometry. The geometric setting of our first example is analogous to that used in [16,17], and can be thought of as a two-dimensional model for mammography guided EIT, see [12,13]. The domain is a rectangle, with electrodes attached on the boundary segments along x 1 = 0 and x 1 = H, while the lateral segments along x 2 = ±W have no electrodes, simulating breast tissue between two compression plates which are equipped with contact electrodes, see Figure 2. In the computed examples, we restrict the measurement to one half of the domain, Ω 1 , while the other half, Ω 2 , is truncated, that is, The Dirichlet-to-Neumann (DtN) map is defined along the interface Γ = {(x, 0) | 0 < x < H}. Assuming that the conductivity is constant in Ω 2 , σ 2 = σ 0 = constant, we have an explicit formula for the DtN mapping: Observe that the Fourier coefficients λ n grow linearly with n, since the Dirichletto-Neumann map is a anti-smoothing operator, and the infinite sum converges in the weak sense as a quadratic form. 3.1.2. The half-space geometry. The second example is a two-dimensional model for the EIT problem in geophysics [15], which is also relevant in some medical imaging configurations, see [22]. The domain Ω is the lower half plane, and the electrodes are attached on a finite interval −R < x 1 < R of the boundary, see Figure 3. In this case the domain of interest is a semi-disc, while the rest of the half space is truncated, Figure 3.

3.1.3.
A circular geometry. In the third model geometry that we consider, the electrodes are attached to the boundary of a unit disc, , and the truncation surface is a circle of radius R < 1, so that Figure 4. This geometry is particularly well suited to illustrate how significantly artifacts from domain truncation not handled appropriately, can dominate the solution. The DtN map with constant conductivity σ 2 = σ 0 on Ω 2 is where the harmonic component λ n is similar to that in the half space model, λ n = σ 0 n πR . We discretize the forward problem as described in [16,18], using the finite difference (FD) scheme on a rectangular p × m grid in Ω. It was shown in the cited references that the FD scheme with the complete electrode model, if properly organized, leads to a symmetric positive definite linear system. We modify the matrix so as to account for the DtN map replacing the homogenous Neumann condition on Γ in Ω 1 by the proper boundary condition. To find a computationally efficient discrete approximation for the operator (15) which preserves the symmetry of the approximation facilitating the computation of the Jacobians, we use a Galerkin collocation scheme: The Dirichlet data along Γ is approximated by a piecewise linear function, while the Neumann data is discretized by collocation at the grid points. This leads to the approximation , and ψ k are piecewise continuous rooftop functions with ψ k (x j ) = δ k,j . By truncating the Fourier series approximation at the frequency T , we obtain a matrix approximation of the operator L such that with ω n = nπ/H. The matrix L thus obtained is symmetric and the discretization scheme leads to Fourier components that decay as 1/n. Furthermore, since the Cauchy data along Γ oscillates very slowly, the truncation of the Fourier series is not expected to create problems in the approximation.
To test the approximation numerically, we consider first the case in which the formula (15) is exact: The dimensions of the computational domain are W = 6 cm and H = 3 cm, respectively, conductivity over the whole body, and thus over Ω 2 , is a constant, σ = σ 0 = 0.5 mS/cm. In the simulation, we apply a current pattern in which the electric current is injected through the electrode e 4 and drawn out through electrode e 8 (see Figure 2), that is, I = i 0 [0, 0, 0, 1, 0, 0, 0, −1] T where i 0 = 1.0 mA. Since the current flows through the body close to the truncation surface, it is reasonable to assume that the corresponding voltages are most sensitive to the domain truncation. We compute the voltages in three different ways: (1) Using the entire domain Ω with no truncation, to provide the ground truth; (2) restricting the computation to Ω 1 alone with homogenous Neumann condition on Γ; (3) using the boundary condition provided by the boundary map (15) with σ 0 equal to the true value. To minimize the effects of discretization, we use a fine grid, in which p = 300, m = 600.
It is clear from Figure 5, which shows the computed voltages at the electrodes corresponding to the applied current pattern in the three different ways, that the The conductivity distribution corresponding to a highly conducting inclusion in the truncated domain. The background conductivity is 0.5 mS/cm, and the conductivity of the inclusion is 1.5 mS/cm. voltages corresponding to the full domain and those computed with the proposed method match almost perfectly, while the homogenous Neumann condition produces a visible discrepancy.
In the second computed test, we assume that the truncated domain Ω 2 contains a significant conducting inclusion. The background conductivity in Ω 1 and Ω 2 is as before, σ 0 = 0.5 mS/cm, while the conductivity of the inclusion is 1.5 mS/cm. The inclusion is shown in Figure 6. We perform the same numerical tests as in the previous case, but adjust slightly the constant conductivity appearing in (15).
The approximate DtN boundary condition absorbs very well the truncation error even if the true inclusion on the right is not taken into account. However, even better performance is achieved if the background constant conductivity in the formula (15) is slightly adjusted to account for the energy drop of the current in domain Ω 1 . The computed voltages corresponding to the current feed between the electrodes adjacent to the truncated domain shown in Figure 7 show an almost perfect overlay of the curves that correspond to the entire domain and the truncated domain with the DtN boundary condition. Tables 1-2 show how the residual norms at the electrodes e 4 and e 8 , adjacent to the truncation boundary, change with the background constant used in the computation of the DtN map. The best performance is obtained by increasing the background by a factor of k = 1.2.  Table 1.

Current
Computed voltage differences as a function of the background conductivity, which is assumed constant, in the approximate DtN map. The conductivity is k × σ 0 , where σ 0 is the conductivity of the background without the inclusion. The table lists the residual norm V(Ω 1 ∪ Ω 2 ) − V (Ω 1 , DtN) , on electrode e 4 (top, closest to Γ) over seven independent current patterns.    Table 2.
Computed voltage differences as a function of the background conductivity, which is assumed constant, in the approximate DtN map. The conductivity is k × σ 0 , where σ 0 is the conductivity of the background without the inclusion. The table lists the residual norm V(Ω 1 ∪ Ω 2 ) − V(Ω 1 , DtN) , on electrode e 8 (bottom, closest to Γ) over seven independent current patterns.

Inverse problem: Preliminary tests.
To shed light on the artifacts arising from domain truncations, we consider reconstructions of the conductivities corresponding to different boundary conditions. Here, we adhere to the Bayesian framework of inverse problems, in which all unknown are represented as random variables, the randomness representing the uncertainty about their values. More extensive tests with different geometries as well as an extension for poorly known DtN maps will be presented in the accompanying paper [4].
We start by specifying the prior probability density. Assume that Ω 1 is discretized with a rectangular grid with N grid points r j , and σ 1 in Ω 1 is approximated by grid values where σ 0 is a presumably known background value and X j represents the unknown variation around it. We use a Gaussian prior density with a hard positive lower bound constraint, writing the prior for X as where "∝" stands for "equal up to a multiplicative scaling", and π + (X) is a threshold function, where τ > 0 is a small cut-off value. The matrix D ∈ R N ×N is a second order finite difference matrix. In order to guarantee that the prior is proper, we impose boundary conditions that guarantee the invertibility of D. In this problem, however, the design of boundary conditions, in particular on the truncation surface Γ are crucial, since too committal conditions (e.g., Dirichlet) will overrule any boundary effect, leading to a clutter inside Ω 1 which would be difficult to distinguish from other effects. Therefore, we follow the ideas in [1,2,3], and use the Aristotelian boundary conditions, which are built on a hierarchical manner: The smoothness is first imposed in interior points conditional on boundary values, and then a smoothness prior for boundary values is imposed so that the prior variance over the whole grid is approximately constant. For further details, we refer to the cited articles.
To set up the likelihood,we consider a discretized observation model with additive noise, where V ∈ R L(L−1) is the voltage vector containing the measured electrode voltages on L = 8 electrodes, corresponding to a full frame of L − 1 linearly independent current patterns, H is the forward map assigning the computed voltages to an underlying conductivity density σ 1 in Ω 1 discretized as in (16), represents additive noise, which we assume to be zero mean Gaussian with covariance matrix Σ, and independent of X. With these assumptions, the likelihood becomes The posterior distribution of X, given the measured voltages, is found by Bayes' formula, From the posterior density, we compute an approximation of the Maximum a Posteriori (MAP) estimate of the conductivity X using a standard Gauss-Newton iteration with projections on the cone X j + σ 0 > τ . The forward map and Jacobians are computed using the FD scheme as explained in [18].
In the first test, we assume that a single conductive object is in Ω 1 , while Ω 2 has constant conductivity σ 0 that is known. The data is generated using the full The conductivity distribution used to generate the data is shown in the top row. In the middle row, the reconstruction after four steps of Gauss-Newton iterations towards the maximum of the posterior density with a homogenous Neumann boundary condition along the truncation boundary. In the bottom row, the corresponding reconstruction with the DtN boundary condition. Observe that the two reconstructions are plotted in the same color scale, the boundary artifact in the middle row being out of the scale. Right column: In the top row, the conductivity distribution used to generate data is shown. The middle row shows the reconstruction after four steps of Gauss-Newton iterations towards the maximum of the posterior density. The boundary condition along the truncation boundary is the DtN condition corresponding to a constant background of 0.5 mS/cm. In the bottom row, the background conductivity in the DtN operator is scaled by a factor 1.4 to compensate the higher conductivity due to the immersed conductive blob on the right.
domain Ω with a discretization of a 300 × 600 mesh. The conductive object is a smooth Gaussian blob with the maximum conductivity 3.1 mS/cm, the background being at 0.5 mS/cm. Figure 8 (left column) shows the reconstruction of the conductivity with the homogenous Neumann boundary condition as well as with the Dirichlet-to-Neumann boundary condition along the truncation boundary Γ. The truncation index in the approximate Dirichlet-to-Neumann operator is T = 12. We see that the Neumann condition creates a strong conducting artifact along the truncation boundary, which can be explained by using a simple resistor analog, see Figure 9: The elevated apparent conductivity is needed to compensate for the conductance through the truncated part of the body. However, the Dirichlet-to-Neumann boundary condition is able to suppress the artifact completely. Observe that with only four electrodes on both sides, and a smoothness prior, the resolution of the reconstruction is necessarily quite low. In this example, as the emphasis is on the boundary effect, the computational details are not optimized.
In the second test, a conducting object is immersed in the truncated domain Ω 2 , see Figure 8 (right column). The object increases the effective conductivity of the truncated domain Ω 2 and the Dirichlet-to-Neumann map adjusted for the constant background σ 0 underestimates it, leading again to a conductive artifact at the boundary. However, the effect can be compensated almost completely by increasing slightly the background value σ 0 in the non-local boundary condition. In this simulation the constant apparent conductivity σ app was set to σ app = 1.4 σ 0 . The factor could be easily implemented to be part of the optimization problem, akin to [10,11]. In this simulation, no optimization with respect to the factor was performed, the more comprehensive computational work being addressed in the accompanying article [4].
In the preliminary test cases, the truncated domain can be satisfactorily treated by an approximate DtN boundary condition because the conductivity along the truncation boundary is almost constant, while when the conductivity varies along Γ, the constant background approximation is likely to be an insufficient model, which is a good motivation to develop a more general approximations, which is done in [4].

5.
Conclusions. The present article provides a theoretical background for a domain truncation technique in electrical impedance tomography, based on domain decomposition techniques that employ partial Dirichlet-to-Neumann operators, or Steklov-Poincaré operators, on the truncation boundary. Preliminary numerical examples demonstrate that in canonical geometries, e.g., rectangular or circular domains, it is possible to find satisfactory approximate models that under the extra assumption of constant conductivity along the truncation boundary are easy to implement and lead to satisfactory reduction of the boundary artifacts. More general geometries and conductivity configurations, however, require new ideas that are not necessarily based on spectral decompositions of the boundary operator. A viable way is to write the discretized boundary map as a generic unknown matrix, and declare the matrix as part of the inverse problem. The question that arises then is, how to implement prior information about the boundary condition such that the matrix indeed corresponds to a truncated conductivity distribution in some admissible class. A solution to this problem is proposed in the article [4] along with computed examples. Figure 9. A schematic figure explaining the increased apparent conductivity near the truncation boundary. On top, the true currents through the truncated part of the body (red arrows) are neglected when a homogenous Neumann condition is imposed, implying that these currents need to be compensated by higher current density through Ω 1 . The equivalent resistor analogue (bottom row) elucidates the needed decrease in the apparent resistivity for this compensation.