Boundary integral formulations for the forward problem in magnetic induction tomography

In this paper we present two models for the forward problem of magnetic induction tomography. In particular, we describe the eddy current model, and a reduced simplified model. The error between the reduced and the full model is analyzed in dependence of parameters such as the frequency and the conductivity. In the case of a piecewise constant conductivity we derive a boundary integral formulation for the reduced model. Finally, we comment on numerical results for the forward problem and give a comparison of both models. Copyrightcopyright 2011 John Wiley & Sons, Ltd.


Introduction
Magnetic induction tomography is a non-invasive and contactless imaging method to determine the conductivity and permittivity distribution inside a certain object, like the human body, see, e.g. [1--5]. An array of coils is placed in a ring around the body as depicted in Figure 1. In each coil a time-harmonic current is induced, which generates a time-harmonic magnetic field. This so-called primary magnetic field B p induces eddy currents inside the body, which perturb the magnetic field. The reconstruction is then based on the measurement of the perturbed voltage in an array of receiver coils around the body.
To When considering biologic tissues we may assume a constant permeability (x) = 0 . The coil is modeled by prescribing a given impressed current j i , of which we assume that it is not influenced by the magnetic reaction fields caused by the eddy currents inside the conducting domain. Hence, we conclude is the complex conductivity. When eliminating the magnetic field H, we finally obtain In particular when considering the inverse problem of MIT, for example by using a standard parameter identification approach, the solution of the forward problem depends on the boundary mesh only. Hence it seems to be a natural choice to use a boundary integral equation approach for the approximate solution of the forward problem [6,7], for finite element methods, see, e.g. [8--10]. This paper is organized as follows: In Section 2 we describe the mathematical model of the forward problem, i.e. the material parameters and the modeling of the excitation coils, and we discuss the eddy current model and a reduced simplified model. Bounds for the approximation error of certain quantities computed with the reduced and the eddy current model are presented in Section 3. In Section 4, a boundary element formulation for the reduced model is given, and some numerical results in Section 5 conclude the paper.

Eddy current model
In this section, we formulate the forward problem of MIT by using the eddy current formulation. The use of the eddy current model is justified since the wavelength for operating frequencies of MIT systems, which lie typically between 10 and 100 MHz, is in the range of some micrometers, so the wavelength is small compared to the size of the conductor.
Let ⊂ R 3 be a bounded domain, which represents the conducting object, and let c = R 3 \ be the exterior domain, where we assume a non-conducting material, e.g. air. For x ∈ c we therefore assume (x) = 0. The eddy current model of MIT is obtained by neglecting the displacement currents D in the non-conducting material domain c , i.e. we set (x) = 0 for x ∈ c . Hence, we can rewrite the partial differential equation (2) as a system of two coupled equations, and where we introduce the gauging condition In addition we assume radiation conditions at infinity, i.e.
Moreover, we have to include transmission boundary conditions for x ∈ = * , which can be rewritten as Note that the jump is given by For the solution of the transmission problem (3)-(7), we introduce the decomposition where E p denotes the primary field generated by the exciting coil C without the presence of a conducting object. Correspondingly, E s is the secondary or reaction field caused by the presence of some conducting material. In particular, the primary field E p is a solution of the partial differential equations which can be written as Hence, we obtain a particular solution by the Newton potential Hence, instead of (3) and (4) it remains to solve the following system of partial differential equations, i.e. for x ∈ and and with the radiation conditions If we define the variational formulation of the transmission problem (10), (11), and (12) reads to find E s ∈ H(curl; R 3 ) such that is satisfied for all F ∈ H(curl; R 3 ). The null space of the variational problem (13) can be characterized by the function E (x) = 0 for x ∈ and E (x) =∇ (x) for x ∈ c , where is the unique solution of the exterior Dirichlet boundary value problem We assume that the boundary = * has only one connected component. The variational formulation (13) of the transmission problem (10)-(12) is therefore not uniquely solvable. Hence, we introduce the gauging condition which corresponds to the conservation of charge. We define The variational problem (13) then admits a unique solution E s ∈ V, see [11,12].

A reduced model
The solution of the forward problem using the eddy current model as described in the previous section is computationally rather expensive. Since in most solution algorithms for the inverse problem the forward problem has to be solved quite often, we are interested in a simplified model which also allows a more efficient solution of the forward problem. For this we write the transmission problem (10)- (12) in terms of the A-formulation [13]. Since B is divergence-free, we can represent the magnetic flux density B as the curl of a magnetic vector potential A, we conclude the existence of a scalar potential satisfying where is uniquely determined by the Coulomb gauge By using the decomposition (8) we can write the primary field E p as while for the secondary field E s we obtain

Note that
is a solution of the partial differential equations Now we can rewrite the eddy current model (3) When applying the divergence operator to Equation (17), this gives In addition, we rewrite the transmission boundary condition (7) in terms of A and and obtain In the parameter range of MIT, numerical examples [14] indicate that A s is very small compared to ∇. Therefore, we neglect A s in (19) and (20), i.e. we finally conclude the Neumann boundary value problem where now denotes the potential in the reduced model. Since is not uniquely determined by the Neumann boundary value problem (21) and (22), we introduce the scaling condition Moreover, by neglecting A s in (17) we obtain Hence, we conclude The electric field can finally be obtained by This means that the solution of the full eddy current model reduces to the solution of a Neumann boundary value problem for the Laplace equation, and the evaluation of a Newton potential. Both models are summarized in Table I.
Eddy current model It remains to estimate the error when considering the reduced model instead of the eddy current model, see also [15]. In particular we have to consider the differences − and A s − A s , respectively. For this, we first introduce the Newton potential operator In the case of a vector-valued function u we consider the Newton potential N 0 u component-wise.

Lemma 3.1
Assume ⊂ B r (0). The Newton potential operator N 0 : L 2 ( ) → L 2 ( ) is bounded satisfying Proof By using the Hölder inequality we have The assertion then follows from Schmidt's inequality, i.e.
In particular, we have which concludes the proof.
Let A s be the solution of the eddy current model (17)- (18), in particular by using (18) we can rewrite (17) as Hence we can write A s as the Newton potential Correspondingly, we have where ∇ is chosen such that We therefore conclude Let , ∈ H 1 ( ) be the weak solutions of the Neumann-type boundary value problems (19)- (20) and (21)- (22), respectively. Then there holds the error estimate If we assume q<1, then there holds and Proof From (19) and (21) we first conclude that := − is a solution of the partial differential equation with the Neumann boundary condition Hence, for ∈ H 1 ( ) the weak formulation of the above Neumann boundary value problem reads as from which (31) follows, i.e.

Moreover, with (27) and by using Lemma 3.1 we further have
The variational formulation of the Robin type boundary value problem (19) and (20) reads, for ∈ H 1 ( ), from which the estimate which immediately results in the estimate (32) when we assume q<1. Finally, by using (30) and Lemma 3.1 we have due to (31). Now, (33) follows from (32).

Remark 3.1
As an example we may consider a test problem with the following parameters: In this case, we have where the last result was obtained by using some finite element discretization.

Corollary 3.3
In addition we have an estimate for the error in an arbitrary point Proof By using (30) we have, for x ∈ , The quantity, which is measured in MIT is the voltage in the receiver coil C, i.e.
Hence, we need to evaluate By using integration by parts, and by using where ∇ y (y)×n y for y ∈ k denotes the surface curl of the function .

A boundary element method for the reduced model
In this section we derive a boundary element formulation for the reduced model as described in the previous section. The variational formulation of the Neumann boundary value problem (21), (22) and (23) is to find ∈ H 1 ( ) such that for all ∈ H 1 ( ). For a piecewise constant conductivity (x) we consider a non-overlapping domain decomposition Instead of the global Neumann boundary value problem (21) and (22) we now consider the local boundary value problems, by using (16) together with the transmission boundary conditions, see (7), For the solution of the local partial differential equation in (38) we use the local Dirichlet to Neumann map is the associated Steklov-Poincaré operator [16]. Let H 1/2 ( S ):= H 1 ( ) | S be the skeleton trace space of H 1 ( ). We then have to solve a variational problem to find ∈ H 1/2 ( S ) such that is satisfied for all ∈ H 1/2 ( ). Since the bilinear form in the variational formulation (39) is bounded and H 1/2 ( S )-elliptic, see, e.g. [17], unique solvability of the variational formulation (39) follows.
To describe the application of the local Steklov-Poincaré operators which are involved in the variational formulation (39) we use the symmetric boundary integral operator representations where is the single layer integral operator, is the double layer integral operator, and is the so-called hypersingular boundary integral operator. The mapping properties of all boundary integral operators and therefore of the local Steklov-Poincaré operators S k are well known, see, e.g. [18--21]. For a symmetric boundary element discretization of the variational formulation (39) we introduce a sequence of admissible boundary element meshes S,h with a globally quasi-uniform mesh size h. Let S 1 h ( S ) be the associated boundary element space of piecewise linear continuous basis functions i . By S 1 h ( k ) = S 1 h ( S ) | k we denote the localized boundary element space of local basis functions k,i , and by k = A k we describe the localization of the global degrees of freedom. The symmetric boundary element approximation of the variational problem (39) results in the linear system, see, e.g. [16], are local boundary element matrices, and S 0 h ( k ) are local boundary element spaces of, e.g. piece-wise constant basis functions k, . Moreover, the right-hand side in (41) is given locally as The stability and error analysis of the symmetric boundary element discretization of the variational problem (39) is well established, see, e.g. [16], and the references given therein.

Numerical results
As conducting domain we first consider the cylinder where the transmitting coil is modeled as a current loop of radius 0.04 which is centered at (−0.14, 0, 0.1) , see Figure 2. The vector normal to the current loop points into the direction of the x 1 -axis, i.e., n x = (1, 0, 0) . Inside the cylinder we place a ball with radius r = 0.02, whose center lies in the point (−0.06, 0, 0.1) . The background conductivity of the cylinder is = 0.1, and the conductivity of the inscribed ball is inc . Figure 2 shows the magnitude of the electric field |n×(E| ×n)| for inc = 0.1. In Figure 3 we give a comparison of the reduced model with the full     Based on the above results we conclude that the reduced model describes an appropriate approximation of the full eddy current model as used in MIT models.
In a second example we consider the model of a human thorax with two lungs, see Figure 5. The volume mesh consists of 83 514 volume elements and 15 641 volume nodes, while the boundary element mesh consists of 13 076 boundary elements ‡ http://www.hpfem.jku.at/ngsolve/. muscle = 0.3618 S / m, while the conductivity of the lungs is lung = 0.2716 S / m, see [23]. The center of the transmitting coil was placed in the point (0, −0.2, 0) , the normal vector of the coil is given by (0, 1, 0) , and its radius is 0.05. In Figure 5 we plot the magnitude of the tangential trace of the electric field, i.e. |n×(E| ×n)|. The field lines of the primary magnetic field B of the secondary magnetic field B s are given in Figure 6.

Conclusions and outlook
In this paper we derived two models which describe the forward problem of MIT, an eddy current problem and a reduced model. We proved estimates for the error between the reduced and the eddy current model, and we formulated a boundary element method for the reduced model. Numerical examples show that the reduced model is a good approximation for the eddy current model in the parameter range of MIT.
To be able to reconstruct the complex conductivity distribution inside the human body with MIT, an inverse problem has to be solved. For the reconstruction of the location and of the shape of organs, a shape reconstruction approach in combination with a level set method can be used. For such an approach, the application of boundary element methods seems to be advantageous, since in every step of the level set method only the boundary has to be remeshed. The solution of the inverse problem demands a very fast solution of the forward problem on meshes with a high number of degrees of freedom. To establish a fast solver for the forward problem, fast boundary element methods may be employed such as the fast multipole method or hierarchical matrices.