Boundary element method and simulated annealing algorithm applied to electrical impedance tomography image reconstruction

Physics has played a fundamental role in medicine sciences, specially in imaging diagnostic. Currently, image reconstruction techniques are already taught in Physics courses and there is a growing interest in new potential applications. The aim of this paper is to introduce to students the electrical impedance tomography, a promising technique in medical imaging. We consider a numerical example which consists in finding the position and size of a non-conductive region inside a conductive wire. We review the electric impedance tomography inverse problem modeled by the minimization of an error functional. To solve the boundary value problem that arises in the direct problem, we use the boundary element method. The simulated annealing algorithm is chosen as the optimization method. Numerical tests show the technique is accurate to retrieve the non-conductive inclusion.


INTRODUCTION
Electrical impedance tomography (EIT) is a technique to obtain the inner image of an object exploring its electrical properties. It consists to apply an electric voltage (or electric current) stimulus through electrodes positioned on the object external surface. The corresponding response (electric current or voltage) is measured by the same electrodes. The obtained data are supplied to a computer, which has an algorithm to reconstruct the conductivity distribution inside the object. This distri- * olavohmenin.ifsp@gmail.com bution can be interpreted as the interior object image. A scheme of a medical EIT is depicted in Fig. 1.
The first occurrence of imaging the human body in terms of its electrical properties occurred in 1978 by Henderson and Webster [1]. After this study, several research groups have emerged providing significant developments in EIT technique. Current applications of EIT in medical area are the monitoring of pulmonary functions, breast cancer detection, blood flow monitoring inside the heart and imaging of human brain functions. Also, in engineering, EIT applications includes the monitoring of twophase flow in a pipe, detection of corrosion or cracks inside metallic objects and retrieval a pipeline route in the subsoil [2][3][4].
Although the EIT technique presents attractive fea- tures, such as low cost, portability and robustness, it is not widely used yet due to the difficulty to obtain a good image resolution. This difficulty arises because the EIT image reconstruction is an inverse-problem and therefore, it is intrinsically ill-posed [5]. The Hadamard criteria, (i) existence, (ii) uniqueness and (iii) continuous data dependence on the solution are not simultaneously guaranteed. The first criterion is satisfied, since the physical body being imaged certainly has a actual conductivity distribution. With respect to the second one, it has been shown for the Dirichlet-to-Neumann map there is a unique solution for the conductivity distribution inside the domain [6]. However, the third criterion has been the hardest one to be overcome. To deal with it, the researchers have tried to find a good and suitable regularization function or alternative approaches. For example, it is possible to restrain the search-space to speed up the convergence to the actual solution and avoid nonphysical solutions [7]. A typical way to solve the reconstruction EIT problem is to treat it as an optimization problem. In this case, one must minimize an error functional which is a value expressing the discrepancy between two different models of the same problem, one experimentally obtained and the other one numerically calculated. The conductivity distribution that yields the error function global minimum corresponds to the sought image.
Due to the ill-conditioned nature of the inverse EIT problem, the optimization surface presents irregularities such as multiple local minima and almost flat regions. Hence, one must adopt an efficient optimization procedure to handle such topological features and we have chosen the Simulated Annealing (SA) [8,9] . This is a stochastic method in which the main feature is the possibility of eventually accept a solution that increases the objective function. This allows the method to detrap from local minima basins or flat regions and reach the global minimum [10,11].
To carry out the iterative searching, one must solve the direct problem with different conductivity distributions. It corresponds to a boundary value problem (BVP), which models the stimulus/response process. Frequently, this BVP does not have analytical solution, requiring the use of numerical methods, such as the Finite Ele-ment Method (FEM) and the Finite Difference Method (FDM). Here, the Boundary Element Method (BEM) was chosen. This is well suited to EIT problem, since it requires only the discretization of the domain boundary and the solution on the boundary is calculated firstly without calculating the values of the electric potential inside the domain [12,13]. Moreover, BEM offers important advantages such as great flexibility for arbitrary geometries and boundary shape and easiness of implementation [14,15].
Our goal in this paper is to present to students the EIT technique using BEM and SA algorithm. This study was motivated by the difficulty in teaching the physical and mathematical methods/processes related to image reconstruction. Some studies have been carried in experimental approach, such as the assembly proposed by Mylott et al. [16] to study the computed tomography. Considering the students' growing interest in computational physics in the last years, we have proposed a numerical example to introduce the EIT technique. It consists in retrieving a non-conductive region inside a conductive wire. In engineering, this non-conductive region may represent, for example, a manufacturing defect or a normal wear in the wire which could compromise its use in an apparatus. In medicine, it may model an air bubble inside an artery hampering the blood flow.
Our numerical example considers a conductive cylindrical wire of length 10.0 and diameter 1.0 (arbitrary units). The non-conductive region is a sphere of radius 0.3 with center located along the cylinder axis at 7.0 units far from one of its extremities, as shown in Fig. 2. The challenge is summarized as a two variables optimization problem, the center x-coordinate and radius R of the non-conductive region. Due to the domain symmetry, we have taken a longitudinal section of the wire as a two-dimensional domain and placed the left bottom corner at the origin of a Cartesian system x0y. Also, the stimulus/response process consists in applying an electric voltage V ab = 12 V between the wire extremities a and b (V a = 12 V and V b = 0 V) and measuring the current flux on the same extremities. Moreover, the lateral area of the wire is considered electrically insulated, such that the current flux vanishes. In Sec. 2, we review the EIT mathematical modeling based on the functional approach, the Boundary Ele-ment Method and the Simulated Annealing algorithm. In Sec. 3, we show the numerical tests and the results. Finally, in Sec. 4, we discuss the results and conclude.

PROBLEM STATEMENT AND METHODS
Here, we present the mathematical and numerical basis of EIT problem, specially focused on our numerical example.

I. EIT mathematical formulation
The mathematical modeling of the EIT problem is obtained from the Maxwell's equations. Since the stimulus processing is made at low frequency, the inductive effects can be ignored. Also, considering the domain fully ohmic, we can neglect the capacitive effects. The electric field is E = − ∇φ, where φ is the electrical potential and ∇ is the gradient operator. The current density is J = −σ E = −σ ∇φ, being σ the electrical conductivity. Considering no internal current sources, ∇ · J = 0 and the partial differential equation (PDE), governs the electrical potential φ inside the domain Ω.
To simulate the stimulus/response process, we have two boundary conditions given by φ = φ and −σ ∂φ where φ is the voltage profile and J is the normal component of the electric current flux, both on the boundary ∂Ω.
The numerical example considers a two-phase medium problem. One of them (the spherical non-conductive region) has null-conductivity and the other one (the rest of the wire) has constant and uniform non-vanishing conductivity. In this case, the interior of the non-conductive inclusion is interpreted as the exterior of the domain and the conductivity distribution σ corresponds to the position and size of the non-conductive region. Also, Eq. (1) is simplified to the Laplace equation Considering unitary the conductivity of the non null phases (σ = 1), the boundary conditions become − ∂φ(x, y) ∂n = J(x, y), for (x, y) ∈ ∂Ω. (4) The stimulus-response can be mathematically modeled through the Dirichlet-to-Neumann map, in which a known voltage profile φ is imposed on the boundary and the electric current fluxes J are calculated on the same boundary. Otherwise, if a known electric current flux profile J is imposed and the voltage φ is calculated, the model is called Neumann-to-Dirichlet map. Although in practice the EIT problem is formulated as a Dirichletto-Neumann map, mathematically, it represents a mixed problem since it is imposed null current flux for the internal boundary, i.e., Eq. (4) vanishes.
The direct problem consists in solving the BVP given by Eqs. (2) and (3), considering a knowing conductivity distribution σ to obtain the current flux from Eq. (4). However, in EIT image reconstruction, σ is unknown and we have an inverse problem. One approach to solve it and to obtain σ is to construct and minimize an error functional (objective function) that compares the measurements obtained from two different models, the actual model and the numerical model.
For the Dirichlet-Neumann map, the actual model corresponds to a collection of measurements, J actual , which are the current fluxes on the electrodes obtained from experimental assembly. This model contains the unknown actual conductivity distribution, σ actual , which must be found. The numerical model considers a known prospective conductivity distribution, σ prosp , to solve the direct problem to obtain J num , which corresponds to the numerical current fluxes on the same electrodes. The error functional must be defined as a discrepancy measure between J actual and J num . Here, we have defined the error functional e(σ prosp ) as the mean square error between the experimental and numerical current fluxes, where J (i) actual and J (i) num are, respectively, the actual and numerical fluxes in the electrode i and n is the number of electrodes.
The idea of the optimization process is to perform an iterative search to find the prospective conductivity distribution that minimizes Eq. (5). In each iteration, a prospective solution, σ prosp , is generated and supplied to Eqs. (2) and (3) to yield a numerical flux J num . This flux is compared with the actual one J actual to calculate the error functional, from Eq. (5). This process is repeated so that e (σ prosp ) ≈ 0 and, consequently, σ prosp ≈ σ actual within an acceptable error level.

II. Boundary element method
The essence of BEM is to convert the field equation to integral equations on the domain boundary through the reciprocal relation, where φ 1 and φ 2 are two general solutions of Eq. (2), and through the fundamental solution, which is a particular solution of Eq. (2). From Eqs. (6) and (7), it is possible to obtain two integral equations, one for the points (ξ, η) inside the domain Ω and another for the points (ξ, η) on the boundary ∂Ω. In the discretized form, these integral equations are for (ξ, η) ∈ Ω, and Φ(x, y; ξ, η)ds(x, y), The external boundary is discretized in N ext straight elements, ∂Ω (k) , k = 1, 2, ..., N ext , and the internal one in N int elements, ∂Ω (k) , k = N ext + 1, N ext + 2, ..., N T . The coordinates (x (k) , y (k) ) of the elements extremities are supplied to the program and, by convention, the elements are numbered following the counter clockwise direction for the external boundary and clockwise for the internal one. This guarantees the normal unitary vector n always points to the outside of the domain. Also, for each boundary (external and internal), the final extremity of the last element is made to coincide with the initial extremity of the first one. In this case, the domain and the inclusion are approximated by polygons. Equation (9) produces a linear system of N T equations with N T unknowns, either φ or J. Solving this system provides the pair (φ, J) for all boundary elements. Hence, the direct problem solution of EIT is completed. If the solution φ at some internal point is required, it can be calculated through Eq. (8).

III. Simulated annealing algorithm
The simulated annealing algorithm is a stochastic optimization method that has its origins in the Metropolis acceptance criterion when two system configurations are compared [17]. It is based on the probability to find the system with energy E, given by the Boltzmann weight, e −E/kB T , where T is the temperature and k B is the Boltzmann constant.
In 1983, Kirkpatrick et al. [18] solved the salesman problem using the Metropolis criterion adding an important differential. They adopted a cooling schedule for the temperature to control the search stochasticity. In the thermodynamic framework, if the temperature of a liquid material is slowly cooled down, the atoms arrange themselves to form a structure (a perfect crystal), which has the lowest internal energy state. However, if the cooling process is not sufficiently slow, the final structure is not perfect and the internal energy is not the lowest one.
The analogy with an optimization process is made considering the objective function f ( x) as the energy E of the system, the solutions x i as the system configurations and the temperature T becomes a control parameter of the process. The optimization process is done iteratively starting from an initial solution x 0 and temperature T 0 . At each iteration, k searches are performed. Each new solution, x i+1 (i = 0, 1, 2, ..., k), is generated through a pre-defined visitation distribution and compared with the current one, x i , to be, or not, accepted according to the where ∆f = f ( x i+1 ) − f ( x i ) and k B = 1. After the last solution x k is evaluated, the temperature is decreased through a cooling schedule and the search process restarts. The last solution accepted at the previous temperature is taken as the initial solution to the current one.
The most common ways to generate a new solution x i+1 from the previous one x i are through uniform or Gaussian deviates. Here, we have chosen the Gaussian one in the following form where T 0 is the initial temperature, T t is the temperature at the iteration t, λ is a parameter that depends on the search space range and η is a random number with normal distribution N (0, 1). According to Eqs. (12) and (13), at high temperatures, solutions can be generated relatively far from the previous one, characterizing a global search. Also, there is high probability to accept a new solution that increases the objective function. As the temperature decreases, the search becomes essentially local and the acceptance of an uphill error solution becomes less likely to occur. Hence, the temperature T plays a fundamental role. The choice of a suitable cooling schedule is essential to ensure a good performance of the optimization process, i.e., to reach the global minimum of the objective function with as few steps as possible. The most commons cooling schedules in the literature are the logarithmic and the geometric ones. We have chosen the last one in which the temperature T decreases with the iteration t according to where T 0 is the initial temperature and α ∈ (0, 1) is the cooling rate.

NUMERICAL TESTS AND RESULTS
To solve the numerical example proposed in Sec. 1, we developed in MATLAB c language a program that includes the BEM and the SA algorithm routines.
To construct the error functional, Eq. (5), one must confront two models of the same problem, the actual and the numerical ones. The actual model is simulated numerically, since it is not the goal of this paper to explore the experimental part of the problem. Hence, the actual current flux, J actual , is obtained solving numerically the Eq. (2) through BEM with the boundary conditions, φ (0, y) = 12 V and φ (10, y) = 0 V , for y ∈ [0, 1], and J (x, 0) = J (x, 1) = 0, for x ∈ [0, 10]. Also, it was adopted J (x, y) = 0 for the points (x, y) on the boundary of the non-conductive spherical region. The numerical model is defined as the actual one. However, its nonconductive region, called prospective inclusion, has center x-coordinate (x p ) and radius (R p ) variables in order to carry out the searching process.
For both models, the external boundary discretization was made using 220 elements, being 100 elements for each side, top and bottom, and 10 elements for each side, left and right. The internal inclusion was discretized in 80 elements, totalizing 300 boundary elements.

I. Test 1 -Error functional behavior
The first test explores the behavior of the error functional surface. The variables x p and R p were systematically modified over the search space

II. Test 2 -Applying the Simulated annealing
In this test, we evaluated the performance of the developed program to retrieve the position and radius of the non-conductive region. Since SA algorithm is stochastic, its convergence to the actual solution is not always guaranteed. Therefore, it is necessary to run the program many times to assess its performance. Also, to make a more careful analysis, we have divided the Test 2 in three parts.
In the first part, Test 2.a, the intent was to find only the center x-coordinate of the non-conductive region. We set the prospective inclusion radius equal to the actual one (R p = 0.3). The initial solution of x p was generated randomly in the range [1.0, 9.   Finally, the third part, Test 2.c, considers the search for both, position and radius. In this case, we joined the two programs used in the previous parts in a single one. More specifically, in each iteration (at each temperature) the program carried out two separated searching, one for the position and another for the radius. Due the difficulty to carry out the optimization searching with two variables, we adopted 2000 iterations, with k = 1, and α = 0.97. As the previous parts, we set T 0 = 1000, λ = 0.8, to generate a new x p , and λ = 0.04, to generate a new R p . Again, the initial solutions of x p and R p were generated randomly in the same ranges of the previous parts. For each part of the Test 2, the program was executed 50 times and the final solution x p and/or R p found in each one was stored. Table I shows the minimum, maximum, mean and standard deviation of the final solutions considering all of the runnings for each part. For the best running of the Test 2.c, the solution found was x p = 7.0005 and R p = 0.3000, with error functional of 2.6 × 10 −17 . For the best running, the evolution of the error functional, the prospective center x-coordinate and the prospective radius during the optimization process are shown in Fig. (4).

DISCUSSION AND CONCLUSION
We have presented to students the electrical impedance tomography technique and the application of retrieving the size and position of an non-conductive inclusion inside of a conductive wire using the boundary element method and simulated annealing algorithm. Numerical simulations assessed either the behavior of the error functional in the search-space and the performance of the implemented program. The results show that the problem has a difficult error surface to be optimized. Despite this, the implemented program has presented a satisfactory accuracy to find the non-conductive region. Finally, the presented algorithms and methods have a wide scope of application and should be better appreciated in topics of undergraduate and graduate physics programs.