ECG dipole source localization in 3D human model

Chia sẻ: _ _ | Ngày: | Loại File: PDF | Số trang:12

Thêm vào BST

Báo xấu

13
lượt xem 1
download

Download Vui lòng tải xuống để xem tài liệu đầy đủ

In this paper, source is considered as a single equivalent dipole and genetic algorithm, an efficient and robust optimization method, is proposed to estimate parameters of source. Some improvements are introduced to enhance performance of conventional genetic algorithm.

Chủ đề:

Bình luận(0) Đăng nhập để gửi bình luận!

Lưu

Nội dung Text: ECG dipole source localization in 3D human model

JOURNAL OF SCIENCE OF HNUE Natural Sci., 2010, Vol. 55, No. 6, pp. 101-112 ECG DIPOLE SOURCE LOCALIZATION IN 3D HUMAN MODEL Dang Thanh Trung(∗) and Nguyen The Loc Hanoi National University of Education Nguyen Thi Thanh Van and Pham Thi Mai Bao College of Technology, Vietnam National University Said Elnaffar College of Information Technology, UAE University, UAE (∗) E-mail: trungdt@hnue.edu.vn Abstract. Computing electrical source in the human heart from potentials on the body surface is an ill-posed problem of electrocardiogram (ECG). It is very difficult to get an exact solution for this problem. In this paper, we would like to introduce a novel method to estimate the solution. For details, we use Finite Element method for the forthcoming problem in con- juction with optimization search algorithm for inverse problems. There are many optimization searches such as: iterative technique, Newton or Down- hill simplex etc., which can be applied to find the best source whose poten- tials best match to the observed potentials in a least square sense. Such methods, however, often converge to a local minimum and their results are affected by initial parameters. In this paper, source is considered as a single equivalent dipole and genetic algorithm, an efficient and robust optimization method, is proposed to estimate parameters of source. Some improvements are introduced to enhance performance of conventional genetic algorithm. In addition, a 3D volume conductor model of the human body, which must be a realistic one, namely accurately reflecting the internal geometrical body structure, is constructed based on an anatomic atlas for numerical tests. A comparison between our approach and one using Downhill Simplex method is implemented. The results show that our approach is stable and may provide a good scheme for solving the ECG inverse problem. Keywords: Bioelectric, Electrocardiography (ECG), Inverse problem, Fi- nite Element Method, Genetic Algorithm. 1. Introduction The ECG inverse problem is the bioelectric inverse problem in which the po- tential values at a limited number of measurement points on boundaries are known and we have to reconstruct the unknown sources generating these data. Solving the 101
Dang Thanh Trung, Nguyen The Loc, Nguyen Thi Thanh Van, Pham Thi Mai Bao and Said Elnaffar ECG inverse problem is made difficult, because it is often ill-posed in the Hadamard sense [1]. The solution does not depend continuously on the data, small errors in the measurement of the torso potentials or thorax geometry can yield unbounded errors in the solution. Because of its importance and usefulness, it has attracted numerous researchers to devote themselves to it in the past decades (Rudy, Mesinger and Rapport 1988, Huiskamp and van Oosterom 1989, Furukawa et al. 1989, T. Musha et al. 1998, Jaakko Malmivuo and Robert Plonsey 1995, Gulrajani et al. 1988, etc.). And several methods are developed for this problem such as solving Gabor-Nelson equations di- rectly (Nelson 1981), combination of Brody shift equation and Levenberg Marquardt method (Gulrajani 1985). Recently, a new approach which yields a more precise so- lution, has been studied to get the solution of inverse problem by combination of a numerical analysis for forward problem, such as Finite Element Method (FEM), Boundary Element Method (BEM), Finite Difference Method (FDM), etc. with an iteration technique such as downhill simplex iteration (T. Musha 1999), simulated annealing (Gerson 1994), Newton-Raphson or Levenberg Marquardt (Xanthis 2006) etc. [1, 3]. These methods, however, may not perform consistently and are also sensitive to the noise in the data. This study proposed using genetic algorithm, an efficient and robust method, to estimate the parameters of electrical sources in the heart. To solve this problem, the best we can do is to build a parametric model and try to adjust the unknown parameters based on the available observations. In this paper, the electrical source is considered as a single moving dipole (e.g 3 parameters for position, 3 parameters for orientations and magnitude) lying entirely within a finite, with inhomogeneous volume conductor to simulate the electrical activity of the heart. Simulation models of a volume conductor (model of thorax) are built from a classical anatomic atlas [11] and stereo matching technique [7]. First, we use a numerical method to solve the forward problem. Because our discussion is limited to an inverse problem in inhomogeneous volume conductor, finite element method [5] is exclusively used in this problem. By this approach, the volume conductor (solution domain) is discretized into a number of finite elements that are connected via nodes. At each node, the governing differential equation is approximated by an algebraic expression, called interpolation function. These interpolation functions are then substituted into the integral equation, integrated and combined with the results from the solution domain to obtain a system of equations. Finally, the system is solved for unknown variables. Next, Genetic Algorithm (GA) [4] is applied to compute the parameters of location, orientation and magnitude of a set of dipoles whose potential fields best match to measurement potentials in least squares sense. Mathematically, it is a very difficult problem because the geometry of the model is inhomogeneous and its 102
ECG dipole source localization in 3D Human Model objective function is very complex. GA is a relatively effective approach for this problem. Although GA might not find the best solution, it can find a near perfect solution with acceptable time. In GA, each solution, called chromosome, is encoded as a bit string (binary form). In the first step of GA, an initial population of random chromosomes is gen- erated. A new population is formed by GA operators such as: selection, evaluation, crossover and mutation, and replaces the current population in the next generation. This process repeats in the specified number of iterations, and chromosomes will be developed generation by generation under control of objective functions to reach the best solution. Because the conventional GA has a very poor local performance, some improvements have been added to enhance its performance. A 3D simulation is constructed from a set of classical anatomic atlas for numer- ical testing. To estimate the performance of our approach, our results are compared with a traditional approach using downhill simplex method [9] (T. Musha et al., 1999). The comparison shows that our results were more accurate and stable than those recovered by the T. Musha method. 2. Modelling 2.1. Dipole modelling In this paper, we encode the solutions as bit strings (Figure 1). The length of string depends on the required precision. The source is considered as a single moving dipole which has the following independent variables: position of dipole r(xc , yc , zc ), orientation (α, β) and magnitude m (Figure 2). For sake of simplicity, in our calculation, we assumed that length of dipole, d, is a constant (in our ex- periment, d=0.1cm). A chromosome which is composed of the parameters of dipole (xc , yc , zc , α, β, m) is encoded as the Figure 3. 103
Dang Thanh Trung, Nguyen The Loc, Nguyen Thi Thanh Van, Pham Thi Mai Bao and Said Elnaffar 2.2. Torso Model A three dimensional volume conductor modeling is required, which must also be a realistic one, namely accurately reflecting the internal geometrical body struc- ture, in order to get the desirable agreement with the measurement. In our experi- ment, a model is constructed from 31 layers (starts from 300th cross-section to 600th cross-section) which was obtained from the visible human viewer project [11] and stereo matching technique [7]. The actual distance between successive cross sections is 1cm. In order to use it in the inverse problem it is imperative to minimize the element number. Since the complexity is mostly due to the curved boundaries be- tween the tissues, the number of elements can be drastically reduced when arbitrary shaped elements are used. Figure 3. Three-dimension mesh model The volume includes regions with different conductivities [10] the conductivi- ties of the body tissues in this study are shown in the Table 1. Table 1. Conductivities of body tissues Element Conductivity (S/m) Thorax 0.21 Lung 0.04 Heart 0.10 For the application of the FEM, the volume conductor is divided into tetra- hedral elements by Deform 3D software. Using this software, a list of nodes or co-ordinates and triangular surfaces for human body structure is created. From this data, Deform 3D software will generate a 3D-mesh model with precision depending on users requests. In our experiment, the model is divided into M=27072 tetrahedral elements and N=5626 nodes. Because the electric field around the source becomes very strong, the mesh points required for FEM are concentrated near source to increase the accuracy in calculation. 104
ECG dipole source localization in 3D Human Model 3. Experimental method 3.1. Problem formulation and Solution Mathematically, most bioelectric field problems can be formulated in terms of Poissons equation as following [1] ∇.σ∇Φ = IV in Ω. (3.1) Where, Φ is the potential field, σ is the electrical conductivity tensor, IV is the current per unit volume within the volume conductor Ω and ∇ is the divergence operator. The forward problem would be solved Equation (3.1) for Φ with known given description of IV , and the Neumann boundary condition (Equation (3.2)) and the Dirichlet condition (Equation (3.3)) [1]: σ∇Φ.n = 0 on ΓT (3.2) Φ = Φ0 on Σ ⊆ ΓT (3.3) Because of the complex geometries and inhomogeneity of volume conductors, it is very difficult to find an exact solution of this problem. Most numerical techniques for solving this kind of problem require that the continuous domain (volume con- ductor) be broken up into discrete elements, called mesh or grid, and then solutions of the volume conductor problem can only be obtained by employing a particular numerical approximation method. In our study, we apply the FEM to solve this problem. The primary idea of this method is discretizing the solution domain into elements which is characterized by interpolation function. The governing differen- tial equation is approximated by a system of equations. The FEM is started by Poissons Equation (3.1) with boundary conditions (3.2), (3.3). Such problem can be formulated in terms of minimizing the potential energy function Π [2]: " 2 2 2 # σ ∂φ ∂φ ∂φ Z Z Π= + + dΩ − IV φdΩ (3.4) 2 ∂x ∂y ∂z Ω Ω The first step in FEM is discretization. For the sake of simplicity, we use tetrahedron elements and the linear approximation function of potential within a typical tetrahedron element, represented by the following equation [5] φ(x, y, z) = a + bx + cy + dz (3.5) 105
Dang Thanh Trung, Nguyen The Loc, Nguyen Thi Thanh Van, Pham Thi Mai Bao and Said Elnaffar Figure 4. A tetrahedron element The total potential energy can be rewritten in the following form by contri- bution of element vectors and matrix [8]: KΦ = H (3.6) The matrix K contains all geometry and conductivity information of the model. Because the basic function is non zero at a few intervals and the number of nodes in mesh is great, the matrix K is sparse and large. Solving linear equation systems (6) with the specific boundary conditions, we can calculate the potential at every node within the volume conductor. Then, the solution of forward problem can be calculated by formulation [8]: Φ = K −1 H (3.7) 3.2. Solution of inverse problem The solution of ECG inverse problem can be reached by using FEM for the cur- rent problem in conjunction with GA in estimation the parameters of dipole sources. Genetic Algorithms (GA) are capable of searching for global functions which cause difficulty for gradient based methods [4]. Principal advantages of GA are domain independence, non-linearity and robustness. Because of these characteristics, GA is suitable for this problem. 3.2.1. Evaluation function The evaluation function (or alternatively called objective functions) is formed to represent the difference between the observed potential fields and the calculated potential fields at the limited nodes and has the following equation: v uN u P0 u (Φi − Ψi )2 u i=1 f =uu N0 (3.8) 2 P Ψi t i=1 106
ECG dipole source localization in 3D Human Model where N0 is the limited number of selected nodes on boundary, Φi , Ψi are correl- atively calculated and measured potentials at these nodes. Our problem is finding the dipole whose potential fields minimize this evaluation function. 3.2.2. Genetic algorithm The genetic algorithm to define the source of the ECG inverse problem may be summarized as follows: - Step 1: Generating an initial population of random chromosomes. To increase diversity of population, the chromosomes are chosen so that they are different from each other. - Step 2: For each individual chromosome in the current population, apply- ing the problem to compute the potential fields within the volume conductor by Equation (3.7) and evaluation function is computed at the selected nodes on the boundary by Equation (3.8). - Step 3: Using the GA operators to create a new population for the next generation: selection, crossover and mutation. • Selection operator: Using the tournament selection. Only one individual from each subgroup whose evaluation function is lowest is chosen to reproduce for candidate population. • Crossover operator: 2-point crossover is used with probability, pc, for each chromosome. Also, a selection algorithm is applied to reduce identical chro- mosomes. • Mutation operator: The purpose of this operator is to allow the GA to avoid the local minima. Each chromosome involves a probability, pm , that an arbitrary bit will be changed from its original state. The new population generated with these operators replaces the old popula- tion. The Step 2 and Step 3 are repeated with this population until a termination condition is satisfied. - Step 4: The generational process is repeated until a termination condition has been reached. In our programme, the process will be terminated after a finite number of iterations. - Step 5: Solution of inverse problem is a chromosome which minimizes the evaluation function in the current population of the last generation. However, conventional GA has a very poor local performance because of the random search. To get a good solution, great computational costs are inevitable. Some improvements are necessary to enhance performance of GA. In this section, we introduce some important modifications that can dramatically improve the per- formance of the conventional GA [4]. 107
Dang Thanh Trung, Nguyen The Loc, Nguyen Thi Thanh Van, Pham Thi Mai Bao and Said Elnaffar • Elitist strategy: Two sets of solutions are stored in our algorithm: a current population and an elite set. The elite set keeps the best solutions at each gen- eration. After the genetic operators are implemented, the elite set is updated by choosing the best individuals of new population. By preserving the good solutions, we can avoid losing some excellent solutions. • Hybrid Algorithm: In many cases, the association of GA with other heuris- tic algorithms can significantly improve the performance of conventional GA. In this study, we propose a hybrid algorithm which combines GA with a local search procedure (hill climbing). To reduce the computational cost, the local search procedure is applied only for the elite set. 4. Simulation Parameters In this simulation, the evaluation function is estimated at N0 =111 electrodes placed surrounding the body. These electrodes were enough to locate the single dipole source. Some parameters used in GA are shown in Table 2. Table 2. Parameters of GA Parameter Attribute Precision 10−5 Bit length 15 Crossover probability (pc ) 0.95 Mutation probability (pm ) 0.015 Population Size 3000 Iteration 1000 Elite Size 2 Group Size 2 We assume the actual trajectory of dipole to be a circle, which has a radius of r = 2cm, lies on a plane with normal vector, ~n = (0.5, 0.7, 1.1), and centre point, (x0 , y0, z0 ) = (35, 33, 25) (Figure 5). The formulations of position, orientation and magnitude of dipole in respect to time, t, are represented by equations: x(t) = r ∗ cos(2πt) + x0 y(t) = r ∗ sin(2πt) + y0 nx ∗ (x − x0 ) + ny ∗ (y − y0 ) + nz ∗ (z − z0 ) = 0 m(t) = 0.1 ∗ cos2 (2πt) + 1 α(t) = 2πt; β = 2πt2 The dipole moments are also changed by time and calculated by the following equa- 108
ECG dipole source localization in 3D Human Model tions: Px (t) = d ∗ m ∗ sin(α) ∗ cos(β) Py (t) = d ∗ m ∗ sin(α) ∗ sin(β) Pz (t) = d ∗ m ∗ cos(α) In the experiment, we chose 50 (Np = 50) dipoles at time points t = 0, 0.02, 0.04, , 0.98. The length of dipole is a constant (d = 0.1cm). With the parameters of GA as in Table 2, the execution time (in elapsed CPU seconds) is approximately 60 seconds per dipole. In Mushas approach, we used randomly 1000 different initial simplexes, and the maximum number of iterations is 5000. In this case, the execution time to determine one dipole is approximately 10s. The programme was run at least 3 times with each approach to confirm the stability of the programme. Figure 6 shows the values of objective function at selected time points in our approach and the T. Mushas approach. Our approach always converged to the global minimum regard- less of the starting parameter estimates while simplex often converged to different solutions for different starting parameter estimates Figure 7 shows some other measures of parameters of dipole. The assumed results are represented with blue lines (solid line) and calculated results are repre- sented with red lines (dot line). The parameters calculated by our approach are shown in Figure 7a. Figure 7b shows the results of the T. Mushas approach. Obvi- ously, our parameters are closer to the parameters of assumption dipoles than those gained from T. Mushas approach. 5. Results and Discusion In order to investigate the solution method, we define Err function as a mea- sure of difference between assumed and calculated solutions in the simulation: v u Np c uP u (fi (t) − fia (t))2 Errf = i=1 t (5.1) Np 109
Dang Thanh Trung, Nguyen The Loc, Nguyen Thi Thanh Van, Pham Thi Mai Bao and Said Elnaffar Figure 7. Component parameters of dipole where, superscripts a, c denote assumption and calculated values, respectively. Np is the number of dipoles. The error of solution in the simulation models can be judged by the average of distance between assumed and calculated dipole location as follows D E= (5.2) Np where D is the summed Euclidean distance of calculated and assumed dipoles and has the equation as following Np X p D= (xci − xai )2 + (yic − yia )2 + (zic − zia )2 (5.3) i=1 Table 3 shows some errors computed by Equations (5.2) and (5.3). 110
ECG dipole source localization in 3D Human Model Table 3. Some errors of dipole components Parameter Our apporach Mushas approach Errx (cm) 0.22453 0.72628 Erry (cm) 0.15717 0.6038315 Errz (cm) 0.14650 0.66305 Er (cm) 0.21412 1.06150 Errp x(mA.cm) 0.03426 0.70728 Errp y(mA.cm) 0.38662 0.27002 Errp z(mA.cm) 0.02256 0.88154 Ep (mA.cm) 0.30087 0.63321 We used 50 dipoles with different orientations, magnitudes and positions as above. These dipoles are enough to estimate the effectiveness of our approach. The source localization was carried out three times to confirm the reliability of the obtained solution. The results in Table 3 show that the errors of dipole locations are acceptable. The average of position errors is approximate 0.22cm, and the average of dipole moment errors is about 0.3mA.cm. Also, the execution time is reasonable, one dipole running on one computer consumed approximately 60 seconds for calculation. Although, our approach execution time is longer than the Mushas one, its still acceptable to get a good solution. These analysis results point out that this approach is a promising way to get stable and highly accurate solutions. 6. Conclusion We represented a novel method to solve the ECG inverse problem in the bio- electric field by finite element method in conjunction with genetic algorithm. We have first introduced finite element methods to solve the problem. The solution of the inverse problem is solved by genetic algorithm. We also introduce some ad- ditional features to overcome the disadvantages of conventional genetic algorithm. The performance of our approach was evaluated based on a 3D realistic torso model and compared with a traditional approach using the downhill simplex method pro- posed by T. Musha. Computer simulation results show that our approach is very effective and stable. The results are as good as what we expected. Further work along these lines is ongoing. Acknowledgment We would like to express our gratitude to the Applied Information Technology Centre for giving us permission to commence this study, to do the necessary research work, to use the equipment and departmental data. We also wish to thank all those in the Faculty of Information Technology (HNUE) for all their help, support, interest and valuable guidelines. 111
Dang Thanh Trung, Nguyen The Loc, Nguyen Thi Thanh Van, Pham Thi Mai Bao and Said Elnaffar REFERENCES [1] Jaakko Malmivuo and Robert Plonsey, 1995. Bioelectromagnetism. Oxfoxd Uni- versity Press, p. 133. [2] R. M. Gulrajani, 1998. The Forward and Inverse Problem of Electrocardiography. IEEE Engineering in Medicine and Biology, 17 (5) pp. 84-101. [3] Robert S. MacLeod, Dana H. Brooks, 1998. Recent Progress in Inverse Problem in Electrocardiology. IEEE Engineering in Medicine and Biology, 17 pp. 73-83. [4] D. T. Pham and D. Karaboga, 2000. Intelligent optimization techniques: genetic algorithm, tabu search, simulated annealing and neural network. Springer, p. 51. [5] P. P. Silvester, R. L. Ferrari, 1996. Finite Elements for Electrical Engineers. New York, Cambridge University Press, p. 28. [6] Jonathan Richard Shewchuk, 2002. Delaunay Refinement Algorithms for Tri- angular Mesh Generation. Computational Geometry: Theory and Applications, 22(1-3), pp. 21-74. [7] Di Stefano, M. Marchinonni, S.Mattoccia, G. Neri, 2004. A fast Area-Based Stereo Matching Algorithm. Image and Vision Computing, Vol 22, pp. 983-1005. [8] P. T. M. Bao, D. T. Trung, N. T. Loc, Said Elnaffar, 2009. ECG Dipole Source Localization by Genetic Algorithm. The 4rd International Conference on Knowl- edge, Information and Creativity Support Systems, Seoul, Korea. [9] J. Hara, T.Musha, W.R.Shankle, 2002. Approximating dipoles from human EEG activity: the effect of dipole source configuration on dipolarity using single dipole models. IEEE Engineering in Medicine and Biology Society, Volume 46, Issue 2, pp. 125-129. [10] Stanley Rush, J.A.Abildskov, Richard McFee, 1963. Resistivity of body tissue at low frequencies. Circ. Res, 22(1), pp. 40-50. [11] Y. Chang, P. Coddington and K. Hutchens, 1999. The NPAC/OLDA visible hu- man viewer Computer Science Department (Adelaide University, Adelaide, Aus- tralia). http://www.dhpc.adelaide.edu.au/projects/vishuman2/ 112