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

Método dos elementos de contorno e algoritmo de recozimento simulado aplicados na reconstrução de imagem da tomografia de impedância elétrica

ARTIGOS GERAIS

Olavo H. Menin^{I, 1}; Vanessa Rolnik^II; Alexandre S. Martinez^II

^IInstituto Federal de Educação, Ciência e Tecnologia de São Paulo, Sertãozinho, SP, Brasil

^IIFaculdade de Filosofia Ciências e Letras de Ribeirão Preto, Universidade de São Paulo, Ribeirão Preto, SP, Brasil

ABSTRACT

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 electrical 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.

Keywords: electrical impedance tomography, boundary element method, simulated annealing algorithm, inverse problem, optimization.

RESUMO

A física tem tido um papel fundamental nas ciências médicas, especialmente em diagnósticos por imagem. Atualmente, as técnicas de reconstrução de imagem já são ensinadas nos curso de Física e existe um crescente interesse em possíveis novas aplicações. O objetivo deste trabalho é apresentar aos alunos a tomografia de impedância elétrica, uma promissora técnica de imageamento em medicina. Para isso, consideramos um exemplo numérico que consiste em encontrar a posição e o tamanho de uma região não condutora no interior de um fio condutor. Nós revisamos o problema inverso da tomografia de impedância elétrica modelado pela minimização de um funcional de erro. Para resolver o problema de valor de contorno que surge no problema direto, nós usamos o método dos elementos de contorno. O algoritmo de recozimento simulado foi escolhido como método de otimização. Testes numéricos mostram que a técnica é precisa para encontrar a inclusão não condutora.

Palavras-chave: tomografia de impedância elétrica, método dos elementos de contorno, algoritmo de recozimento simulado, problema inverso, otimização.

1 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 distribution can be interpreted as the interior object image. A scheme of a medical EIT is depicted in Fig. 1.

Figure 1
- Layout of the electrical impedance tomography. An electrodes belt is attached on the body surface. A voltage source applies a stimulus on the electrodes and the corresponding re- sponse is measured and supplied to a computer. A specific software reconstructs the interior body image.

The first occurrence of imaging the human body in terms of its electrical properties occurred in 1978 by Henderson and Webster ¹. 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 two-phase flow in a pipe, detection of corrosion or cracks inside metallic objects and retrieval a pipeline route in the subsoil ^[2-4].

Although the EIT technique presents attractive features, 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 Element 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 . Also, the stimulus/response process consists in applying an electric voltage between the wire extremities b and 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.

Figure 2 - Scheme of the numerical example: cylindrical conductive wire with length 10:0 and diameter 1:0 with a non-conductive region of radius 0:3 located 7:0 units far from extremity a.

In Sec. 2, we review the EIT mathematical modeling based on the functional approach, the Boundary Element 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.

2 Problem statement and methods

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

2.2 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 , where is the electrical potential and is the gradient operator. The current density is , being the electrical conductivity. Considering no internal current sources, 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 , where is the voltage profile and 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 ,, the boundary conditions become

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 are calculated on the same boundary. Otherwise, if a known electric current flux profile 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 Dirichlet-to-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 _prosp as the mean square error between the experimental and numerical current fluxes,

where 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_actual. This flux is compared with the actual one J_num to calculate the error functional, from Eq. (5). This process is repeated so that _prosp 0 and, consequently, within an acceptable error level.

2.2 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 and 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

for , where

The external boundary is discretized in N_ext straight elements, , and the internal one in N_int elements, . The coordinates 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. Solving this system provides the pair 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).

2.3 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 . It is based on the probability to find the system with energy E, given by the Boltzmann weight, , where T is the temperature and k_b is the Boltzmann constant.

In 1983, Kirkpatrick et al.^[28] 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 xi 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 E and temperature x0. At each iteration, k searches are performed. Each new solution,, is generated through a pre-defined visitation distribution and compared with the current one, x_i, to be, or not, accepted according to the Metropolis criterion. If , the new solution is accepted and if , the new solution is accepted with probability

where . 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 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₀ is the initial temperature, T_t is the temperature at the iteration t, is a parameter that depends on the search space range and n 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₀ is the initial temperature and is the cooling rate.

3 Numerical tests and results

To solve the numerical example proposed in Sec. , we developed in MATLAB 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, , for , and , for . Also, it was adopted 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 non-conductive region, called prospective inclusion, has center (x)-coordinate Rp and radius (Rp) 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.

3.1 Test 1 - Error functional behavior

The first test explores the behavior of the error functional surface. The variables x_p and RP were systematically modified over the search space [2.0,8.0]×[0.1,0.4] in a regular mesh of 61×61=3721 points. The error functional is calculated for each position xp and radius RP generating the error functional surface, as shown in the plot of Fig.3 with error functional axis in logarithm scale. It is possible to see almost flat regions and a narrow channel that contains a prominent global minimum. These characteristics require a powerful optimization method to reach the global minimum.

Figure 3 - Error functional surface in logarithm scale. It presents almost at regions and a channel that contains a prominent global minimum.

3.2 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 (Rp=0.3). The initial solution of was generated randomly in the range [1.0,9,0] and, to generate a new solution xp from the previous one by Eq. (13), it was chosen . In the second one, Test 2.b, the program searched only the non-conductive region radius. In this case, we set the prospective inclusion center x-coordinate equal to the actual one (xp=7.0). The initial solution of Rp was generated randomly, but in the range [0.05,0.45]. The new solution from the previous one by Eq. (13) was generated considering For both, Test 2.a and Test 2.b, the program was executed considering 1000 iterations, with k =1, and we set =0.95 and T₀ = 1000.

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 to generate a new x_p, and , 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 1 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).

Table 1 - Minimum, maximum, mean and standard deviation of the nal solutions x^p and/or R^p for the three parts of the Test 2, considering 50 runnings.

Figure 4
- Evolution with the iteration for the best running of Test 2.c of the (a) error functional (b) prospective center x-coordinate and (c) prospective radiu

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.

Acknowledgments

ASM acknowledges Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) (305738/2010-0). VR acknowledges Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) (2008/01284-4).

Recebido em 24/7/2012

Aceito em 2/2/2013

Publicado em 24/4/2013

• [1] R.P. Henderson and J.G. Webster , IEEE Transactions on Biomedical Engineering 25, 250 (year 1978).
• [2] W. Lionheart , N. Polydorides and A. Borsic , Electrical Impedance Tomography: Methods, History, and Applications, (Institute of Physics Publishing, Bristol, year 2005), Chap. chapter The reconstruction problem.
• [3] J.F. Abascal , S.R. Arridge , D. Atkinson , R. Horesh , L. Fabrizi , M.D. Luca , L. Horesh , R.H. Bayford and D.S. Holder , journal 43, 258 (year 2008).
• [4] A. Osella , G. Chao and F. Sanchez , American Journal of Physics 69, (year 2001).
• [5] L. Borcea, Inverse Problems 18, R99 (2002).
• [6] A.I. Nachman, The Annals of Mathematics, Second Series 143, 71 (1996).
• [7] O.H. Menin and V. Rolnik, International Journal of Modern Physics C 22, 825 (2011).
• [8] E. Aarts and J. Korst, Simulated Annealing and Boltzmann Machines: A Stochastic Approach to Combinatorial Optimization and Neural Computing (J. Wiley & Sons, 1988).
• [9] D. Goldberg, Genetic Algorithms in Search, Optimization and Machine Learning (Addison-Wesley Longman Publish, Boston, 1989).
• [10]T.C. Martins, E.D.L.B. Camargo, R.G. Lima, M.B.P. Amato and M.S.G. Tsuzuki, in Engineering in Medicine and Biology Society, EMBC, 2011 Annual International Conference of the IEEE (2011).
• [11] H. Kim, C. Boo and Y. Lee, International Journal of Control, Automation and System 2, 211 (2005).
• [12] V. Rolnik, O.H. Menin, G.L.C. Carosio and P.J. Seleghim, in ABCM Symposium Series in Mechatronics, Vol. 4 (2010) pp. 835{843.
• [13] R. Duraiswami, K. Sarkar and G.L. Chanine, Engineering Analysis with Boundary Elements, 13(1998).
• [14] C.A. Brebbia and J. Dominguez, Boundary Elements Methods An Introductory Course, 2nd ed. (Computational Mechanics Publications Southampton Boston,1992).
• [15] W.T. Ang, A Beginner's Course in Boundary Element Methods (Universal Publishers, Boca Raton, Florida, 2007).
• [16]E. Mylott, R. Klepetka, J.C. Dunlap and R. Widenhorn, European Journal of Physics 32, 1227 (2011).
• [17]N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller and E. Teller, Journal of Chemistry Physics 21, 1087 (1953).
• [18] S. Kirkpatrick, C. Gellat and M. Vecchi, Science 220, 671 (1983).

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

Método dos elementos de contorno e algoritmo de recozimento simulado aplicados na reconstrução de imagem da tomografia de impedância elétrica

Publication Dates

• Publication in this collection
  05 July 2013
• Date of issue
  June 2013

History

• Received
  24 July 2012
• Accepted
  02 Feb 2013