# The Factorization method for three dimensional Electrical Impedance Tomography

###### Abstract

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.

## 1 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

### 2.1 The continuum forward model

Let be an open bounded set with Lipschitz boundary . Then, according to the continuum model, the electric potential inside produced by the injection of a current on the boundary solves

(1) |

where is the unit normal to directed outward to and the conductivity is a real valued function such that there exists for which

Let us denote . It is well known that whenever , problem (1) has a unique solution . Then, we can define the so- called Neumann-to-Dirichlet (NtD) map corresponding to the conductivity by

where is the unique solution to (1).

### 2.2 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 of electrodes. Let us denote by these electrodes, for each , 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 when . We also introduce the gap between the electrodes. Following the presentation in [13], let us introduce the subspace of of functions that are constant on each electrode and that vanish on the gaps between the electrodes

where stands for the characteristic function of some domain . For simplicity, in the following we will not make the distinction between an element of and the associated vector of . Let us denote by the injected current which is constant and equal to on electrode . The voltage potential then solves

(2) |

where is the so-called contact impedance and is the unknown measured voltage potential on the electrodes. We assume that the contact impedance satisfies for almost all .

For all there exists a unique 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 . Similarly to the continuum model, we define the (finite dimensional) Neumann-to-Dirichlet map associated with the conductivity by

where is the unique solution to (2).

## 3 The inverse problem

### 3.1 Statement of the inverse problem

Let us denote by the conductivity in the presence of an inclusion. We assume that for almost all and that there exists a domain such that

where is either a positive or a negative function on and is the known conductivity of the background and is such that for almost all . We assume moreover that and that is connected with Lipschitz boundary. The inverse problem we treat consists in finding the indicator function from the knowledge of the finite dimensional maps and .

As stated in the next section (see Theorem 3.1), the Factorization method provides an explicit formula of the indicator function from the knowledge of and . Nevertheless, as we explain later on, and actually give a good approximation of the characteristic function of provided that the number of electrodes is sufficiently large and that they cover most of .

### 3.2 The Factorization method to find the support of an inclusion

For each point in , let us define the Green’s function for the background problem with Neumann boundary conditions which is the solution to

(3) |

where is the Dirac distribution at point . Then,

(4) |

is the electric potential created by a dipole located at point of direction such that .

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 .

###### Theorem 3.1.

Take such that and let us denote the eigenvalues and eigenvectors of the self adjoint and compact operator . Then the characteristic function of is given by

(5) |

###### Remark 3.2.

In [16] the authors state a convergence result that justifies the use of the finite dimensional map instead of to obtain an approximation of the characteristic function of . We recall here briefly the main result they obtain. Let us introduce the orthogonal projector defined by

and let us denote by the difference data operator. Then, if the number of electrodes goes to infinity and if the gap between the electrodes decreases sufficiently fast, we deduce from [16, Theorem 8.2] that for every there exists a truncation index such that for a given the sequence

converges when goes to infinity if and only if is in . In this expression, are the eigenvalues and eigenfunctions of the finite dimensional map . In practice, it is not easy to determine the truncation index 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 . This result tells us that the function

should be an approximation of the indicator function of in the sense that should be greater inside than outside .

###### Remark 3.3.

As in the continuous setting (see Remark 3.2), for each point the characteristic function can be defined by using the inverse of the squared norm of the solution to the linear equation

## 4 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.

### 4.1 Data sets and numerical implementation of the indicator function

In practise, one does not know the full and noiseless map but only the map where is the set of currents that one injects in the tested object and for all we have

where denotes the noise in the data. As a consequence, the data consist of the matrix where is the dimension of . 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 by using finite elements implemented with the software FreeFem++ (see [10]) to solve equations (2) for each current in the admissible set of currents . We refer to [21] for more details on the variational formulation we use to solve (2). Then, we build by adding artificial noise to that is

where is a matrix of size whose element is a random number generated with a normal distribution and is a real number chosen such that

(6) |

is a given level of noise. The denotes the term by term multiplication between two matrices and denotes the Frobenius norm.

Let us introduce the singular values and singular vectors of that satisfy

for all . Then, the function introduced in the previous section corresponds to the inverse of

since

solves (see Remark 3.3). As mentioned before, we take the truncation index as large as we can which is equal to the dimension of the range of . Finally, to limit artefacts, we will use the function

as an indicator function of where is a set of unit vectors of . The choice of the number of dipole directions is based on experimental observations.

### 4.2 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 and points 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 and for we can decompose as

where

and solves the following conductivity problem:

(7) |

One can actually compute and realise that this function is singular for but in . Therefore, by elliptic regularity we deduce that is in and we compute it by using the finite element software FreeFem++ to obtain a good approximation of for each electrode position and each point in . In the case of an homogeneous background, we have verified the accuracy of our approximation by comparing the finite element solution of (7) with a boundary element solution computed with the software BEM++ (see [19]). In the case of a homogeneous background ( 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 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: on the quality of the reconstructions. To this end, we introduce the new indicator function

where corresponds to with replaced by .

## 5 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 and 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 for the indicator function which we normalise between and . 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].

### 5.1 Preliminary experiments

In Figures 1 to 4 we show reconstructions for various geometries, for constant or piecewise constant background conductivities and for inclusions located at different locations. In all 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 is the barycenter of , (respectively ) is the barycenter of the region delimited by the closed iso-surface of level of the normalised indicator function obtained from (respectively from ). Finally, 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 is not simply connected, we compute and 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 in (6).

Figures | Position |
---|---|

Figures 1(b), 1(c) and 2(b) | (0,5,2) |

Figures 1(d), 1(e) and 2(c) | (0,5,2) and (5,-2,2) |

Figures 3(b), 4(b) and 4(c) | (0,0,0) |

Figures 3(c), 4(d) and 4(e) | (40,40,0) |

#### First setting: cylinder with one ring of electrodes (Figure 1)

We consider a cylindrical domain of radius and height with one ring of electrodes containing electrodes (see Figure 1(a)). The space V of input currents is made of independent vectors of 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 except for two of them that are geometrically opposite to each other. One of these is set to , the other one to so that the constraint on the input current is satisfied. In this experiment, the background conductivity is taken constant and equal to while the contact impedance value is for all the electrodes (we keep this value for all experiments). In Figure 1(b) we see that with one inclusion, the algorithm with as test function finds the correct (x,y) location of the inhomogeneity but not the location whereas the algorithm with (Figure 1(c)) does not even find the location. The reason why Ind gives the correct location but not the one is because we only have one ring of electrodes at the same position. We also run an experiment with two well separated inclusions and the algorithm with (Figure 1(d)) seems to give quite accurate results in this case as well while the use of (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 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 matrix. The results are similar to the previous case, except that the resolution in the direction in much better in this case since we accurately find the 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 location with a good accuracy.

#### Third setting: human head with homogeneous background (Figures 3)

We use a head shape domain with 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 Neumann-to-Dirichlet matrix. We tried two positions for the inclusion which is a sphere of radius and of conductivity with a background conductivity of . We still use a constant contact impedance which is 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 of the background is not constant. We take a piecewise constant conductivity with values in the yellow part, in the green part and 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 and of conductivity 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.

### 5.2 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 iso-surface, 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.

### 5.3 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 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 errors, we deform the boundary of the domain by expanding it in the direction and and by contracting it in the direction 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 , respectively , the center of the electrodes corresponding to , respectively . We then simulate the Neumann-to-Dirichlet matrix data in the domain with inexact electrode positions and use the Green’s function of evaluated at the exact electrode position to produce an image. In Figure 8 we report results for and being such that the mean relative errors on the shape and the electrode positions are

where stands for the distance function and , respectively , contains the position of the 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 given for by

(8) |

for 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 . In Figure 9 we plot the iso-surface of value 0.9 of the indicator function Ind for an inclusion located at . 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 gave accurate reconstructions. The results reported in Figure 9 show that up to 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)).

## 6 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 accuracy in the challenging and realistic case of a head shape domain with homogeneous and inhomogeneous background conductivities from noisy simulated data obtained with 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 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].

## Acknowledgement

The research of N. C. is supported by the Medical Research Council Grant MR/K00767X/1. The authors are grateful to the EIT group in the Department of Medical Physics and Bioengineering of University College London for providing the meshes of the head. The authors would like to thank the anonymous referees for their valuable comments and suggestions.

## References

- [1] K. A. Awada, D. R. Jackson, J. T. Williams, D. R. Wilton, S. B. Baumann, and A. C. Papanicolaou. Computational aspects of finite element modeling in EEG source localization. IEEE Transactions on Biomedical Engineering, 44(8), 1997.
- [2] M. Brühl. Explicit characterization of inclusions in electrical impedance tomography. SIAM J. Math. Anal., 32(6):1327–1341, 2001.
- [3] M. Brühl and M. Hanke. Numerical implementation of two noniterative methods for locating inclusions by impedance tomography. Inverse Problems, 16:1029–1042, 2000.
- [4] M. Brühl and M. Hanke. Recent progress in electrical impedance tomography. Inverse Problems, 19:865–890, 2003.
- [5] A. P. Calderón. On an inverse boundary value problem. Mat. apl. comput., 25(2–3):133–138, 2006.
- [6] K.-S. Cheng, D. Isaacson, J. C. Newell, and D. G. Gisser. Electrode models for current computed tomography. IEEE Transactions on Biomedical Engineering, 36(9):918–924, 1989.
- [7] M. K. Choi, B. Harrach, and J. K. Seo. Regularizing a linearized eit reconstruction method using a sensitivity based factorization method. available at http://www.mathematik.uni-stuttgart.de/harrach/publications/FM_regularization.pdf.
- [8] H. Haddar and G. Migliorati. Numerical analysis of the Factorization method for EIT with piecewise constant uncertain background. Inverse Problems, 29(6):065009, 2013.
- [9] M. Hanke and B. Schappel. The factorization method for electrical impedance tomography in the half-space. SIAM J. Appl. Math., 68(4):907–924, 2008.
- [10] F. Hetch. FreeFem++. www.freefem.org/ff++.
- [11] D. S. Holder, editor. Electrical impedance tomography: methods, history and applications. Institute of Physics, 2005.
- [12] L. Horesh. Novel approaches in modelling and image reconstruction for multi-frequency electrical impedance tomography. PhD thesis, University College London, 2006.
- [13] N. Hyvönen. Complete electrode model of electrical impedance tomography: approximation properties and characterization of inclusions. SIAM J. Appl. Math., 64(3):902–932, 2004.
- [14] N. Hyvönen, H. Hakula, and S. Pursiainen. Numerical implementation of the factorization method within the complete electrode model of electrical impedance tomography. Inverse Problems and Imaging, 1:299–317, 2007.
- [15] A. Kirsch. Characterization of the shape of a scattering obstacle using the spectral data of the far field operator. Inverse Problems, 14(6), 1998.
- [16] A. Lechleiter, N. Hyvönen, and H. Hakula. The factorization method applied to the complete electrode model of impedance tomography. SIAM J. Appl. Math., 68(4):1097–1121, 2008.
- [17] J. Mueller and S. Siltanen. Linear and nonlinear inverse problems with practical applications. SIAM, Computational Science and Engineering, 2012.
- [18] D. Nicolas. Couplage de méthodes d’échantillonnage et de méthodes d’optimisation de formes pour des problèmes de diffraction inverse. PhD thesis, Ecole Doctorale de l’Ecole Polytechnique, 2013.
- [19] W. Śmigaj, S. Arridge, T. Betcke, J. Phillips, and M. Schweiger. Solving boundary integral problems with BEM++. ACM Trans. Math. Software, (submitted).
- [20] E. Sommersalo, M. Cheney, and D. Isaacson. Existence and uniqueness for electrode models for electric current computed tomography. SIAM J. Appl. Math., 52(4):1023–1040, 1992.
- [21] P. J. Vauhkonen, M. Vauhkonen, T. Savolainen, and J. P. Kaipio. Three-dimensional electrical impedance tomography based on the complete electrode model. IEEE Transactions on Biomedical Engineering, 46(9):1150–1160, 1999.
- [22] T. York. Status of electrical impedance tomography in industrial applications. Journal of Electronic Imaging, 10:608–619, 2001.