The Factorization method for three dimensional Electrical Impedance Tomography

The use of the Factorization method for Electrical Impedance Tomography has been proved to be very promising for applications in the case where one wants to find inhomogeneous inclusions in a known background. In many situations, the inspected domain is three dimensional and is made of various materials. In this case, the main challenge in applying the Factorization method consists in computing the Neumann Green's function of the background medium. We explain how we solve this difficulty and demonstrate the capability of the Factorization method to locate inclusions in realistic inhomogeneous three dimensional background media from simulated data obtained by solving the so-called complete electrode model. We also perform a numerical study of the stability of the Factorization method with respect to various modelling errors.


Introduction
Electrical Impedance Tomography (EIT) is an imaging technique that allows retrieval of the conductivity distribution inside a body by the application of a current to its boundary and measurement of the resulting voltage. The portability and the low cost of electronic devices capable of producing such data makes it an ideal tool for non destructive testing or medical imaging. In the last few years, several imaging devices that use EIT have produced interesting results in the field of medical imaging (see [11]) as well as for non destructive testing (see [22]).
The EIT inverse problem has been intensively studied in the mathematical literature since Calderón's paper in 1980 (see [5]) that formulates the inverse boundary value problem. Then, many authors have been interested in proving uniqueness for the inverse problem and in providing efficient algorithms to find the conductivity from boundary measurements (see [17] and references therein for a complete review on this subject). One of these algorithms is the so-called Factorization method, introduced by Kirsch in [15] for locating obstacles from acoustic scattering data and then extended by Brühl in [2] to the EIT inverse problem. The main advantage we see in this technique is that it places fewer demands on the data since it only locates an embedded inhomogeneity but does not give the conductivity value inside the inclusion. The Factorization method and its regularisation have been studied by Lechleiter et al. in [16] in the context of the so-called complete electrode model which was shown by Isaacson et al. in [6] to be close to real-life experimental devices.
The main purpose of this paper is to show that the Factorization method can be successfully applied to realistic three dimensional EIT inverse problems. To our knowledge, the only work presenting three dimensional reconstruction by using the Factorization method in EIT is due to Hanke et al. [9] and they focus on the case where the probed domain is the half space. Therefore, we would like to extend this to more complex geometries with inhomogeneous background conductivities. To do so, the main difficulty consists in computing the Neumann Green's function of the studied domain since this function is needed to apply the Factorization method. In two dimensions with homogeneous background, one can use conformal mapping techniques to obtain the Neumann Green's function of many geometries from that of the unit circle (see [3,14] for example). Clearly, this technique is no longer available in three dimensions and we have to use other ideas. We choose to follow a classical idea presented in [1], and more recently used in [8] in the EIT context in two dimensions, because it can be extended to three dimensional problems and allows to compute Neumann Green's functions for non constant conductivities. It consists in splitting the Green's function into a singular part (which is known) plus a regular part and to compute the regular part which is the solution to a well posed boundary value problem with numerical methods such as finite elements or boundary elements.
We show with various numerical experiments and for different geometries that this technique can be successfully applied in three dimensions to obtain EIT images from simulated data. To make our experiments more realistic, the simulated data we used were produced by using the complete electrode model with a limited number of electrodes which covered a part of the inspected domain, as it is the case in many applications. We first compare the results obtained by using the Neumann Green's function of the searched domain and the one of the free space. Then we study the influence of various parameters on the quality of the reconstructions and we conclude our numerical evaluation of the method by studying the influence on the reconstructions of various modelling errors such as errors in electrode positions, in the shape of the probed domain and in the background conductivity.
In section 2 we present briefly the two main models for the direct EIT problem and introduce some notations and definitions. In section 3 we present the inverse problem and state the main theoretical foundation of this paper. Finally, in section 4, we present our numerical implementation of the method and in section 5 we give numerical reconstructions for different domains obtained with noisy simulated data.
2 Forward models in electrical impedance tomography

The continuum forward model
Let Ω ⊂ R 3 be an open bounded set with Lipschitz boundary Γ. Then, according to the continuum model, the electric potential u inside Ω produced by the injection of a current I on the boundary Γ solves where ν is the unit normal to Γ directed outward to Ω and the conductivity σ ∈ L ∞ (Ω) is a real valued function such that there exists c > 0 for which Let us denote L 2 (Γ) := {f ∈ L 2 (Γ) | Γ f ds = 0}. It is well known that whenever I ∈ L 2 (Γ), problem (1) has a unique solution u ∈ H 1 (Ω). Then, we can define the so-called Neumann-to-Dirichlet (NtD) map corresponding to the conductivity σ by where u is the unique solution to (1).

The complete electrode model
A more accurate model, the so-called complete electrode model, takes into account the fact that in practise the current is applied on Γ through a finite number M ∈ N of electrodes. Let us denote by (E j ) j=1,...,M these electrodes, for each j = 1, . . . , M , E j is a connected open subset of Γ with Lipschitz boundary and non zero measure. According to the experimental setups, we assume in addition that the distance between two electrodes is strictly positive, that is We also introduce G M := Γ \ ∪ M i=1 E i the gap between the electrodes. Following the presentation in [13], let us introduce the subspace of L 2 (Ω) of functions that are constant on each electrode and that vanish on the gaps between the electrodes where χ V stands for the characteristic function of some domain V . For simplicity, in the following we will not make the distinction between an element f of T M and the associated vector (f i ) i=1,...,M of R M . Let us denote by I ∈ T M the injected current which is constant and equal to I i on electrode i. The voltage potential u then solves where z ∈ L 2 (Γ) is the so-called contact impedance and U ∈ T M is the unknown measured voltage potential on the electrodes. We assume that the contact impedance satisfies z(x) ≥ c > 0 for almost all x ∈ Γ. For all I ∈ T M there exists a unique (u, U ) ∈ H 1 (Ω) × T M that solves problem (2) (see [20] for more details about the model and the proof of existence and uniqueness) and this solution depends continuously on I. Similarly to the continuum model, we define the (finite dimensional) Neumann-to-Dirichlet map associated with the conductivity σ by where (u, U ) is the unique solution to (2).

Statement of the inverse problem
Let us denote by σ ∈ L ∞ (Ω) the conductivity in the presence of an inclusion. We assume that σ(x) ≥ c > 0 for almost all x ∈ Ω and that there exists a domain D ⊂ Ω such that where γ ∈ L ∞ (D) is either a positive or a negative function on D and σ 0 ∈ C 0,1 (Ω) is the known conductivity of the background and is such that σ 0 (x) ≥ c > 0 for almost all x ∈ Ω. We assume moreover that D ⊂ Ω and that Ω \ D is connected with Lipschitz boundary. The inverse problem we treat consists in finding the indicator function χ D from the knowledge of the finite dimensional maps Σ M σ 0 and Σ M σ . As stated in the next section (see Theorem 3.1), the Factorization method provides an explicit formula of the indicator function χ D from the knowledge of Λ σ 0 and Λ σ . Nevertheless, as we explain later on, Σ M σ 0 and Σ M σ actually give a good approximation of the characteristic function of D provided that the number of electrodes M is sufficiently large and that they cover most of Γ.

The Factorization method to find the support of an inclusion
For each point z in Ω, let us define the Green's function for the background problem with Neumann boundary conditions N (·, z) ∈ L 2 (Ω) which is the solution to where δ z is the Dirac distribution at point z ∈ Ω. Then, is the electric potential created by a dipole located at point z ∈ Ω of direction d ∈ R 3 such that |d| = 1.
We use these dipole test functions in the next Theorem (which is a slightly reformulated version of [2,Proposition 4.4]) in the expression of the characteristic function of the inclusion D.
Theorem 3.1. Take d ∈ R 3 such that |d| = 1 and let us denote (λ i , ψ i ) i=1,...,+∞ the eigenvalues and eigenvectors of the self adjoint and compact operator Λ σ 0 − Λ σ . Then the characteristic function χ D of D is given by Remark 3.2. An equivalent (and maybe more classical) statement of Theorem 3.1 is the following range test for d a given unit vector of R 3 . Then, (5) is given by the inverse L 2 squared norm of g z that solves |Λ σ 0 − Λ σ | 1/2 g z = φ d z and which is finite if and only if this equation has a solution in L 2 (Γ).
In [16] the authors state a convergence result that justifies the use of the finite dimensional map Σ M σ 0 − Σ M σ instead of Λ σ 0 − Λ σ to obtain an approximation of the characteristic function of D. We recall here briefly the main result they obtain. Let us introduce P M : L 2 (Γ) → T M the L 2 (Γ) orthogonal projector defined by Then, if the number of electrodes M goes to infinity and if the gap between the electrodes decreases sufficiently fast, we deduce from [16,Theorem 8.2] that for every M there exists a truncation index 0 < R(M ) < M such that for a given z ∈ Ω the sequence ..,M are the eigenvalues and eigenfunctions of the finite dimensional map Σ M . In practice, it is not easy to determine the truncation index R(M ) from the data (see [4] for a heuristic method) but since we will consider the case where we have few measurements, we will see that no regularisation is actually needed and we will take R(M ) = M − 1. This result tells us that the function should be an approximation of the indicator function of D in the sense that χ M d should be greater inside D than outside D.

Numerical implementation
In this section we give a precise definition of the data we used for inversion and we present our implementation of the Factorization method and of the computation of the dipole test functions.

Data sets and numerical implementation of the indicator function
In practise, one does not know the full and noiseless map Σ M but only the mapΣ M : V → T M where V ⊂ T M is the set of currents that one injects in the tested object and for all I ∈ V we haveΣ M I = Σ M I + where denotes the noise in the data. As a consequence, the data consist of the M × N matrix Σ M where N is the dimension of V . Each column of this matrix is a vector that contains the voltage at each electrode. In the following, we will consider the case of synthetic data. To produce them, we compute a noiseless map Σ M by using finite elements implemented with the software FreeFem++ (see [10]) to solve equations (2) for each current I in the admissible set of currents V . We refer to [21] for more details on the variational formulation we use to solve (2). Then, we buildΣ M by adding artificial noise to Σ M that iŝ where N is a matrix of size M × N whose (i, j) element is a random number generated with a normal distribution and η is a real number chosen such that is a given level of noise. The · denotes the term by term multiplication between two matrices and · denotes the Frobenius norm.
As mentioned before, we take the truncation index R(M ) as large as we can which is equal to the dimension of the range of Σ M . Finally, to limit artefacts, we will use the function

Computation of the dipole potentials
In three dimensions, the Green's function of the Laplace equation in the free-space is given by and for any unit vector d ∈ R 3 and points x = z we define the associated dipole potential We remark that the image principle used in two dimensions to compute the Neumann Green's function for a circle is not valid anymore in three dimensions for a sphere. Nevertheless (see [1,8] for example), for σ 0 ∈ C 0,1 (Ω) and for z, x ∈ Ω we can decompose φ d z as and V d z (x) solves the following conductivity problem: One can actually compute ∇ xφ d z (x, σ 0 (z)) and realise that this function is singular for x = z but in L 2 (Ω). Therefore, by elliptic regularity we deduce that V d z (x) is in H 1 (Ω) and we compute it by using the finite element software FreeFem++ to obtain a good approximation of φ d z (x i ) for each electrode position x i and each point z in Ω. In the case of an homogeneous background, we have verified the accuracy of our approximation by comparing the finite element solution V d z (x) of (7) with a boundary element solution computed with the software BEM++ (see [19]). In the case of a homogeneous background (σ 0 is constant) there is no source term inside Ω in equations (7). Therefore, as it has been observed in [3], one can first compute the Neuman-to-Dirichlet map Λ σ 0 associated with the continuum model and obtain the solutions to equation (7) with a simple matrix vector product.
In the following, we will study the impact of using the zero average dipole of the free space:

Numerical simulations and error analysis
In this section, we show that the Factorization method successfully applies to three dimensional imaging problems in rather complicated geometries and with partial covering of the boundary Γ with electrodes. We also test the sensibility of the method with respect to various experimental errors such as errors in the shape of the domain, errors in the background conductivity and errors in the electrode's placement. In what follows we choose to not show plots of the indicator functions but plots of an iso-surface of the indicator functions Ind and Ind in red on the images. The choice of the iso-surface is a complicated question and to our knowledge there is no systematic way to do it; we arbitrarily choose to show the iso-surface of value 0.9 for the indicator function which we normalise between 0 and 1. We show in Figure 5 how this parameter influences the quality of the reconstruction. To overcome this known difficulty, a solution would be to use the indicator function obtained with the Factorization as an initial guess for a level set approach as it is done in [18] in the context of acoustic scattering or use it for regularisation of a linear inversion method as it is proposed in [7].

Preliminary experiments
In Figures 1 to 4   cases, the inclusion is a sphere (of radius 1 for the cylinder and 10 for the head shape) and the conductivity value in this inclusion is double that of the background. The different positions of the centre of the sphere are reported in Table 1. In order to compare the reconstructed object with the true one, we plot in green the projection of the true object on the different planes delimiting the plotting region. We also choose a quantitative estimate of the quality of the reconstruction by introducing the relative error on the location of the barycenter which is defined by where C true is the barycenter of D, C Est (respectively C Est ) is the barycenter of the region delimited by the closed iso-surface of level 0.9 of the normalised indicator function obtained from Ind (respectively from Ind). Finally, diam(Ω) stands for the diameter of the computational domain Ω. The errors on the location of the barycenter for the different tests presented in this section are reported in the caption of the plots. When D is not simply connected, we compute E c and E c for each simply connected component and we give the mean of the errors obtained for the different objects. In all this experiments, the noise added to the simulated data is taken such that δ = 1% in (6).
First setting: cylinder with one ring of electrodes (Figure 1) We consider a cylindrical domain Ω of radius 10 and height 7 with one ring of electrodes containing 32 electrodes (see Figure 1(a)). The space V of input currents is made of 16 independent vectors of R 32 corresponding to the so-called opposite current pattern (see [17,Chapter 12] for more details on different current patterns). This means that all the electrodes are set to 0 except for two of them that are geometrically opposite to each other. One of these is set to 1, the other one to −1 so that the constraint on the input current is satisfied. In this experiment, the background conductivity is taken constant and equal to 1 while the contact impedance value is z = 5 for all the electrodes (we keep this value for all experiments). In Figure 1(b) we see that with one inclusion, the algorithm with φ d z as test function finds the correct (x,y) location of the inhomogeneity but not the z location whereas the algorithm withφ d z (Figure 1(c)) does not even find the (x, y) location. The reason why Ind gives the correct (x, y) location but not the z one is because we only have one ring of electrodes at the same z position. We also run an experiment with two well separated inclusions and the algorithm with φ d z (Figure 1(d)) seems to give quite accurate results in this case as well while the use ofφ d z (Figure 1(e)) only gives a rough idea of the location of the objects. Therefore, in the following we will mainly use the indicator function Ind.
Second setting: cylinder with two rings of electrodes (Figure 2) This experiment is similar to the first one (the domain Ω is the same) except that we have two rings of 20 electrodes each. The injection protocol is again a pairwise injection which corresponds to the opposite injection pattern for each ring of electrodes. Then, the data consist  of a 40 × 20 matrix. The results are similar to the previous case, except that the resolution in the z direction in much better in this case since we accurately find the (x, y, z) location of the inhomogeneity if we use the correct dipole function (Figures 2(b) and 2(c)) This experiment illustrates that it is not necessary to have electrodes covering the entire domain to have good quality reconstructions since we do not have any electrode on the top and on the bottom of the cylinder and still we find the z location with a good accuracy.  Reconstruction for a cylindrical geometry with two rings of electrodes by adding δ = 1% of noise to the synthetic data (same colours as in Figure 1).

Third setting: human head with homogeneous background (Figures 3)
We use a head shape domain with 31 electrodes that covers a portion of the physically accessible part of the head (see Figure 3(a)). In this setting, we show that the Factorization method still performs well. The protocol injection is very similar to the previous one and the data correspond to a 31 × 20 Neumann-to-Dirichlet matrix. We tried two positions for the inclusion which is a sphere of radius 10 and of conductivity 2σ 0 with a background conductivity of σ 0 = 1. We still use a constant contact impedance which is z = 5 for all the electrodes. In both cases the location of the inclusion is accurately found.
Fourth setting: human head with inhomogeneous background (Figure 4) We use the same head shape geometry as previously, but this time the conductivity σ 0 of the background is not constant. We take a piecewise constant conductivity with values 1.5 × 10 −4 in the yellow part, 2.0 × 10 −5 in the green part and 4.4 × 10 −4 in the red part of the domain (see Figure 4(a)). These values correspond approximately to the conductivity of the skin, the skull and the brain respectively of a human head (see [12] and references therein). The inclusion is still a sphere of radius 10 and of conductivity 2σ 0 = 8.8 × 10 −4 located at the same places as in the previous setting. Let us mention the fact that for the two locations, the inclusion is inside the brain (the red part of the computational domain). The reconstructions with the exact dipole tests functions are very precise but when we use the dipole function of the free space we observe a misplacement of the reconstructed object. This difference is actually significant when the inclusion is close the inhomogeneous layers (compare Figure 4(d) with Figure 4(e)). The result obtained with the free space dipole function will be used for comparison in section 5.3.
(b) One inclusion in the middle using φ d z ; Ec = 0.007.
(c) One inclusion in the back side using φ d z ; Ec = 0.01. Figure 3: Reconstruction for a head shape geometry with homogeneous background by adding δ = 1% of noise to the data (same colours as in Figure 1).
(b) One inclusion in the middle using φ d z ; Ec = 0.01.
(c) One inclusion in the middle using φ d z ; Ec = 0.08.
(d) One inclusion in the bask side using φ d z ; Ec = 0.01.
(e) One inclusion in the bask side using φ d z ; Ec = 0.13. Figure 4: Reconstruction for a head shape geometry with inhomogeneous background by adding δ = 1% of noise to the data (same colours as in Figure 1).   Figure 1).

Influence of the truncation value and of the inclusion's size on the reconstruction
In the next set of experiments we illustrate the influence of the choice of iso-surface ( Figure  5) and of the size of the inclusion ( Figure 6). All these simulations are performed in the same geometry as in Figure 4 and the conductivity value of the background and of the inclusion are also the same. Therefore, one can compare these results to the reference results presented in Figure 4.
Regarding the influence of the value of the iso-surface that defines the boundary of the inclusion, we remark that it greatly affects the size of the reconstructed object but not its location. As it has already been observed in previous work, in its actually state, the Factorization method is probably not the right tool to estimate the size of a defect since this quantity depends too much on the choice of the iso-surface value. Nevertheless, whatever the value of the isosurface, the object we reconstruct has the correct shape and the correct location.
These remarks also apply to our second test ( Figure 6) where we plot the iso-surface of value 0.9 but this time for inclusions of different sizes. Nevertheless, let us mention the fact that even for a small inclusion (Figure 6(a)) the method locates accurately the defect.

Robustness to modelling errors
We conclude our numerical analysis of the Factorization method in the context of EIT for brain imaging by studying the influence of modelling errors and noise on the reconstructions.
First of all, in Figure 7 we repeat the same experiment as in Figure 4(b) and we plot the iso-surface of value 0.9 of the indicator function Ind for different level of noise added to the simulated data. The size of the reconstructed defect is strongly affected by the noise level but even for 10% of noise we still obtain a very accurate estimate of the location of the inclusion (see Figure 7(c)).
We also test the sensitivity of the method with respect to errors on the shape of the domain Ω and on the electrode's placement. Indeed, this is of crucial importance for applications since one usually has only a rough idea of the shape of Ω and of the electrode positions. We consider in this experiment the homogeneous head that has already been used in Figure 3. We choose to not use the inhomogeneous head of Figure 4 since we would like to decouple errors in the background conductivity and errors in the shape of the domain. To simulate such modelling   (c) One inclusion in the back side using φ d z ; Ec = 0.08. Figure 8: Reconstruction for a deformed head shape geometry with homogeneous background by adding δ = 1% of noise to the data (same colours as in Figure 1). errors, we deform the boundary of the domain Ω by expanding it in the direction x and z and by contracting it in the direction y and to obtain a perturbed domain Ω ε . See Figure 8(a) for a superposition of cut views of Ω and Ω ε , in blue we plot a slice of Ω and in orange a slice of Ω ε . We remark that by doing so we also generate a domain with inexact electrode positioning. Let us denote E ε ∈ R 31×3 , respectively E ∈ R 31×3 , the center of the electrodes corresponding to Ω ε , respectively Ω. We then simulate the Neumann-to-Dirichlet matrix data in the domain Ω ε with inexact electrode positions E ε and use the Green's function of Ω evaluated at the exact electrode position E to produce an image. In Figure 8 we report results for Ω ε and E ε being such that the mean relative errors on the shape and the electrode positions are where dist stands for the distance function and E i ∈ R 3 , respectively E i ε ∈ R 3 , contains the position of the i th electrode on ∂Ω, respectively on ∂Ω ε . The main conclusion we can draw from these experiments is that the Factorization method we presented in this paper is stable with respect to errors on the shape of the measurement domain.
The last source of modelling errors that certainly affects the Factorization method is in the background conductivity value. For this last experiment we go back to the inhomogeneous head introduced in Figure 4. We generate data for a noisy background conductivity σ γ 0 given for γ > 0 by (1 + γ)σ 0 in the yellow part, (1 − γ)σ 0 in the green part, (1 + γ)σ 0 in the blue part (8) for σ 0 being as in Figure 4 and the colours refer to the regions in Figure 4(a). Then, we build the indicator function Ind by using the dipole function associated with the background conductivity σ 0 . In Figure 9 we plot the iso-surface of value 0.9 of the indicator function Ind for an inclusion located at (40, 40, 0). In view of Figure 4 we expect that this position is more affected by error in the background conductivity than the middle position since in the case of the middle position (Figures 4(b) and 4(c)) even the indicator function Ind gave accurate reconstructions. The results reported in Figure 9 show that up to 20% of noise on the background conductivity, the Factorization method still performs reasonably well and gives an accurate estimate of the object location. Let us also keep in mind that in situations where one does not have any knowledge of the background conductivity, the free space dipole functions are in practise good candidates to probe the domain (see Figures 4(c)and 4(e)).

Conclusions
The main focus of this paper was to show via numerical experiments that the Factorization method can be used to solve realistic three dimensional imaging problem with EIT data. We presented a possible implementation of the computation of the Neumann Green's function for three dimensional problems and we provided reconstructions obtained with the Factorization algorithm by using the numerically computed Green's function. In the considered test cases the obtained results are rather accurate. The Factorization method actually finds the location of an inclusion with 1% accuracy in the challenging and realistic case of a head shape domain with homogeneous and inhomogeneous background conductivities from noisy simulated data obtained with 31 electrodes that cover only a part of the domain's boundary. We have also show that the Factorization method is robust with respect to noise on the measurements and  Figure 9: Reconstruction of a lateral inclusion by using φ d z with error on the background conductivity (same colours as in Figure 1).
various modelling errors such as errors on electrode positioning, on the shape of the domain and on the background conductivity.
A direct extension of this work would be to perform an experimental study and to incorporate some systematic criterion to choose the truncation level for valuable singular values and to pick an iso-surface of the indicator function. We think for example of the criterions introduced in [4]. Another possible extension, which is probably the main interest of the Factorization method for applications, is to couple the result obtained by the Factorization algorithm with optimisation algorithms as it is proposed for example in [7].