A Method for Normal Direction Judge Applied in Electrocardiographic Problem

Chengcheng Tang, Jian Wu, Riqing Chen

Department of Biomedical Engineering, Tsinghua University, Beijing, China.

DOI: 10.4236/jbm.2018.61001                

Abstract

Boundary Element Method (BEM) is widely used in electrocardiographic (ECG) problem. Formulations of these problems based on mathematical and numerical approximations of the known source in heart and the volume conductor that can transfer voltages on the surface of the body. To analyze the electric potentials on body surface or epicardial surface, a set of discrete equations derived from a boundary integral equations need to be solved. Solving these equations means to get the potential distribution eventually. In the process of solving, transfer matrix of discrete equations has received considerable attention, how to get an appropriate transfer matrix is an important issue. This paper found that the direction of normal vector could affect the results when calculating the transfer matrix and presents a method analogous to Mesh Current Method to deal with this direction problem. Several simulations have been carried out to verify the accurate results with the correct direction of normal vector using new method within a torso model given simultaneous epicardial and body surface potential recordings.

Keywords

ECG Problem, BEM, Transfer Matrix, Normal Vector, Direction

1. Introduction

Electrocardiographic problem, namely given the potential distribution at epicardial surface to calculate the distribution at body surface (Forward Problem), or vice versa (Inverse Problem), draws lots of attentions in biomedical area these years, as electrocardiography can be a powerful tool for diagnosis in clinic [1] [2] [3]. The potential distribution at epicardial surface and body surface are linked by a transfer matrix determined by the shape of heart and torso. Thus, both forward and inverse problem need to calculate the transfer matrix elements determined by the geometry properties.

The boundary element method (BEM) is often an appropriate way to compute the transfer matrix, despite the different inner physical properties of the conductor medium when the medium is isotropic [4]. Barr et al. gives an overall theory of dealing with forward problem and inverse problem based on BEM [5]. Horácek et al. calculate transfer coefficients based on triangle-to-triangle discretization [6] [7], rather than node-to-node discretization used by Barr et al. to construct appropriate matrices, these matrices are consisted of coefficients that relate to the potentials and their gradients on the epicardial surface of the conductor and the potential on the body surface [8] [9]. The corresponding boundary integral equation should be solved through discretizing the surface into triangle elements [10] [11]. Munck et al. uses analytically integrated elements to discretize the boundary integral equation linearly [12] [13].

The transfer matrix plays an important role in both forward problem and inverse problem. In the calculation of the matrix, the direction of normal vector of the triangle segments should be paid attention to carefully. Some methods are designed to deal with this issue by the commercial programs. Take one method as an example, first, it finds a point that is definitely outside the mesh, and calculates the centers of the triangles. Then, it computes the distances between the centers of the triangles and the outer point which was found previously. After that, it finds the triangle closest to the point that is for sure outside the mesh based on the distances that have been calculated before. Next, the normal vector $\bar{n}$ of this triangle is computed out based on right hand’s rule according to index sequence of the three vertexes, and also the vector from the center of the nearest triangle to the outer point ${\bar{P}}_{out}$ . If ${\bar{P}}_{out}\cdot \bar{n}$ is negative, the order of the vertex index of each triangle would be adjusted to invert $\bar{n}$ . However, this approach lets mistakes occur sometimes, it can’t make all the directions turn into right in some cases, which results in incorrect potential distribution eventually.

Therefore, this paper presents a new method to figure out the correct triangle normal direction to get the accurate potential distribution on surface.

2. Basic Theory andMethod

2.1. Derivation of Transfer Matrix

First, deriving the integral equations. The potential ${\varphi }_o$ at point O can be expressed as Equation (1) using Green’s second law, and the point O is assumed as an observation point which inside the volume but very close to surface,

${\varphi }_o=-\frac{1}{4\pi }{\displaystyle {\int }_{S_H}^{}{\varphi }_H\frac{\bar{r}\cdot \bar{n}}{r^2}d{S}_H}-\frac{1}{4\pi }{\displaystyle {\int }_{S_H}^{}\frac{\nabla {\varphi }_H\cdot \bar{n}}{r}d{S}_H}+\frac{1}{4\pi }{\displaystyle {\int }_{S_B}^{}{\varphi }_B}\frac{\bar{r}\cdot \bar{n}}{r^2}d{S}_B$ . (1)

where $S_B$ is the body surface, $S_H$ is the epicardial surface. r is the distance that extends from the observer to elements of integration dS. ${\varphi }_B$ is the potential distribution on body surface, ${\varphi }_H$ is the potential distribution on epicardial surface. $\bar{r}$ is an unit vector in the direction of r. $\bar{n}$ is an outward pointing vector of unit magnitude normal to surface element dS.

So, two equations can be obtained from Equation (1) for observation positions located on body and epicardial surface, respectively,

$-{\varphi }_B+\frac{1}{4\pi }{\displaystyle {\int }_{S_B}^{}{\varphi }_B}\frac{\bar{r}\cdot \bar{n}}{r^2}d{S}_B-\frac{1}{4\pi }{\displaystyle {\int }_{S_H}^{}{\varphi }_H\frac{\bar{r}\cdot \bar{n}}{r^2}d{S}_H}-\frac{1}{4\pi }{\displaystyle {\int }_{S_H}^{}\frac{\nabla {\varphi }_H\cdot \bar{n}}{r}d{S}_H}=0$ , (2)

$\frac{1}{4\pi }{\displaystyle {\int }_{S_B}^{}{\varphi }_B}\frac{\bar{r}\cdot \bar{n}}{r^2}d{S}_B-{\varphi }_H-\frac{1}{4\pi }{\displaystyle {\int }_{S_H}^{}{\varphi }_H\frac{\bar{r}\cdot \bar{n}}{r^2}d{S}_H}-\frac{1}{4\pi }{\displaystyle {\int }_{S_H}^{}\frac{\nabla {\varphi }_H\cdot \bar{n}}{r}d{S}_H}=0$ . (3)

To solve the Equation (2) and Equation (3) numerically, these equations should convert to simultaneous linear equations. Therefore, the surface S, consist of $S_B$ and $S_H$ , need to subdivide in small triangles, which means the integration area will be partitioned into triangle elements. And Equations (2) and Equation (3) can be simplified by the element of solid angle $d{\Omega }_{ef}^i$ which is subtended at an observation point of the i-th location on surface e by an area element on surface f:

$d{\Omega }_{ef}^i=\frac{\bar{r}\cdot \bar{n}}{r^2}d{S}_f$ . (4)

As it can be seen here, if the direction of $\bar{n}$ is wrong, the sign of $\bar{r}\cdot \bar{n}$ wouldn’t be correct, therefore, it will have influence on the solid angle, which means the result of potential distribution may not be right.

Substituting Equation (4) in Equation (2) and Equation (3), and calculating the surface solid angle integration approximately by discretized triangle solid angle summation as the way shown in Equation (5). The coefficients of the components of gradients term are noted in Equation (6).

$\frac{1}{4\pi }{\displaystyle {\int }_{S_Y}^{}{\varphi }_{Y}d{\Omega }_{XY}}={\displaystyle \underset{1}{\overset{N_Y}{\sum }}{P}_{XY}{\Phi }_Y}$ , (5)

$\frac{1}{4\pi }{\displaystyle {\int }_{S_Y}^{}\frac{\nabla {\varphi }_Y\cdot \bar{n}}{r}d{S}_Y}={\displaystyle \underset{1}{\overset{N_Y}{\sum }}{G}_{XY}{\Gamma }_Y}$ . (6)

The different rows of matrices $P_{XY}$ and $G_{XY}$ corresponds to different locations of the observer on surface X. As the theory showing, each of the P’s and G’s is a matrix of coefficients depending entirely on the geometry. In all cases the number of columns corresponds to the number of locations on the surface Y of integration [1]. Based on the derivation above, the equation about the potential distributions on epicardial surface and body surface can be derived as ${\varphi }_B=Z_{BH}{\varphi }_H$ , where $Z_{BH}$ is the transfer matrix described in Equation (7).

$Z_{BH}={\left(P_{BB}-G_{BH}{G}_{HH}^{-1}{P}_{HB}\right)}^{-1}\left(G_{BH}{G}_{HH}^{-1}{P}_{HH}-P_{BH}\right)$ . (7)

2.2. Normal Direction Judge Methods

As it has been shown in the theory, the direction of $\bar{n}$ is very important, the method of Barnard et al. requires the additional computation of the normal to the plane triangle to establish the correct sign. So the direction of $\bar{n}$ must be confirmed before any further calculation.

A new method is proposed to ensure the direction of $\bar{n}$ of each triangle element to be outward. The basic idea of this method comes from the application of Kirchhoff’s laws. There are two applications of Kirchhoff’s laws and Ohm’s law, the Branch Current method and the Mesh Current method. The Mesh Current method can solve a circuit with less variables and less equations by simply determine unknown currents of each mesh in a network. Like its name, mesh current, means a current that loops around the essential mesh. It encompasses all the component of the mesh. First, there are some loops need to be identified within the circuit. As can be seen in Figure 1(a), one loop formed by $V_{DC}$ , $R_1$ and $R_2$ while the other loop formed by $R_2$ , $R_3$ . These current meshes are connected with each other by loops. $I_1$ and $I_2$ is the current of each mesh respectively. The direction of each current is chosen arbitrarily, and the intersecting component will be “went through” by two currents based on the assumption of the method. Second, get all the voltage drops though each resistor according to the direction of mesh currents. Then use the Kirchhoff’s Voltage laws, get equations state as ${\displaystyle {\sum }_{k=1}^n{V}_k}=0$ . Solve a set of equations can get respective mesh current.

In the theory of Mesh Current method, if all mesh currents are chosen in a same direction, either clockwise or counterclockwise, the direction of two mesh currents that go through the intersecting side of two meshes must be opposite. So, that kind of law can be used in determining the direction of $\bar{n}$ in turn. In other words, if applying the theory of Mesh Current method to determine the direction of $\bar{n}$ , some principles should be followed: the direction of $\bar{n}$ gives the direction of rotation according to the right-hand rule, the direction of rotation of each triangle element shows a direction of each side of the triangle. If two triangles have contrary direction on their common side, it means the directions of $\bar{n}$ of them are at the same side. This can be seen in Figure 1(b).

In that case, if the direction of $\bar{n}$ of a triangle could be determined, the direction of $\bar{n}$ of the neighborhood of this triangle also could be determined. Therefore, in this new method, it needs to figure out a reference normal vector (RNV) at first, then determine the others’ direction according to the direction of

(a) (b)

Figure 1. (a) A circuit for Mesh Current method; (b) Geometry for two adjacent triangle elements. It contains two triangles with four nodes. Triangle t1 with n1, n2 and n3, triangle t2 with n2, n3 and n4. If the order of t1’s nodes is n1, n2, n3 and the order of t2’s nodes is n3, n2, n4, the direction of normal vector of these two triangles are the same.

the RNV. The triangle element of which $\bar{n}$ is ensured correct and is referred as reference triangle (RT). Check the $\bar{n}$ of triangle elements in the neighborhood, which share one of RT’s sides, and then check the neighborhood of these three triangles, and so on. Eventually, we can finish examination on all triangle elements of the model.

3. Simulation Results

The result of simulation can be carried out when using the one auto-solid angle approximation [3]. Data of epicardial surface potential distributions include a number of different sequences of cardiac excitation and repolarization, along with the geometric location (x, y, z coordinates) of each body and heart electrode. For the purpose of contrast, this paper introduces the result of node-to- node discretization as a reference distribution.

The simple model this paper chooses to simulate has 602 epicardial points and 771 body surface points, with a cylinder-shaped heart, which is shown in Figure 2(a) and Figure 2(c). The surface $S_B$ of this model is represented as 1538 triangle elements. The surface $S_H$ of this model is represented as 1200 triangle elements. The comparatively complex model we apply our method for has 337 epicardial points and 771 body surface points, with a relatively real-shaped heart, which is shown in Figure 2(b) and Figure 2(c). The surface $S_B$ of this model is represented as 1538 triangle elements. The surface $S_H$ of this model is represented as 976 triangle elements.

3.1. Simple Model

The simulation on simple model shows in Figure 3. Figure 3(a) gives the potential distribution on epicardial surface which has been set as a known distribution. Figure 3(b) shows the potential distribution on body surface, which is a reference result here. And the potential distribution on body surface calculated by a commercial program and the program with new method is shown in Figure 3(c) and Figure 3(d), respectively.

As it can be seen in Figure 3(a), the potential on heart surface is set to increase from the bottom to the top monotonously. Consequently, Figure 3(b) and Figure 3(d) show a reasonable monotonous potential distribution on the body surface. However, the bottom part of Figure 3(c) have significant differences with Figure 3(b) and Figure 3(d), it shows an unreasonable reverse below the middle area, which means it is inaccurate.

3.2. Comparatively Complex Model

To verify the feasibility of the new method, a more genuine heart model is used for compares. The simulation results on comparatively complex model are shown in Figure 4. Figure 4(a) gives the potential distribution on epicardial surface which has been set as a known distribution. Figure 4(b) shows the potential distribution on body surface, which is a reference result here. And the

(a) (b) (c)

Figure 2. (a) Simple heart model; (b) Comparatively complex heart model; and (c) Body model.

(a) (b) (c) (d)

Figure 3. (a) Potential distribution on epicardial surface; (b) Reference potential distribution on body surface; (c) Potential distribution on body surface calculated by a commercial program; (d) Potential distribution on body surface calculated by the method this paper presented.

potential distribution of body surface calculated by a commercial program and program of the new method are show in Figure 4(c) and Figure 4(d) respectively.

The results shown in Figure 4 are analogous to the results in Figure 3. The potential on heart surface also increases from the bottom to the top monotonously. Figure 4(b) and Figure 4(d) still keep this property. But Figure 4(c) appears some abnormal distribution on the shoulder area, which exhibits its inaccuracy.

4. Conclusion

The direction of the normal vector of discretized triangle element is quite

Figure 4. (a) Potential distribution on epicardial surface; (b) Reference potential distribution on body surface; (c) Potential distribution on body surface calculated by a commercial program; (d) Potential distribution on body surface calculated by the method this paper presented.

important for transfer matrix calculation in the electrocardiographic problem. It could be wrong sometimes if we do not deal with it appropriately. This paper presents a method analogous to mesh current method to judge the normal direction of triangles. Principle and process of this method are analyzed and explained in detail. The new method shows its accuracy in a simple and a more realistic heart shape model in the electrocardiographic problem. This method also has an inconvenient part, namely, an initial $\bar{n}$ of a triangle element needs to be determined first. Sometimes it should be checked manually if needed. What’s more, though the method does not require that either surface has any particular shape such as that of a sphere, it can’t work correctly when geometry is extremely complex. However, the new method this paper presented can deal with most situations with high accuracy.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Coll-Font, J., et al. (2014) New Additions to the Toolkit for Forward/Inverse Problems in Electrocardiography within the SCIRun Problem Solving Environment. Computing in Cardiology Conference, Cambridge, MA, 2014, 213-216.
[2]	Bear, L.R., Cheng, L.K., LeGrice, I.J., Sands, G.B., Lever, N.A., Paterson, D.J. and Smaill, B.H. (2015) Forward Problem of Electrocardiography. Arrhythmia and Electrophysiology, 8, 677-684. https://doi.org/10.1161/CIRCEP.114.001573
[3]	MacLeod, R. and Buist, M. (2010) The Forward Problem of Electrocardiography. In: Macfarlane, P.W., van Oosterom, A., Pahlm, O., Kligfield, P., Janse, M. and Camm, J., Eds., Comprehensive Electrocardiology, Springer Science & Business Media, 247-298.
[4]	Bear, L., Dubois, R. and Zemzemi, N. (2016) Optimization of Organ Conductivity for the Forward Problem of Electrocardiography. Computing in Cardiology Conference, IEEE, Vancouver, BC, 2016, 385-388.
[5]	Barr, R.C. and Spach, M.S. (1977) Relating Epicardial to Body Surface Potential Distributions by Means of Transfer Coefficients Based on Geometry Measurements. Biomedical Engineering, IEEE Transactions on BME, 24, 1-11. https://doi.org/10.1109/TBME.1977.326201
[6]	Horácek, B.M. and Clements, J.C. (1997) The Inverse Problem of Electrocardiography: A Solution in Terms of Single- and Double-Layer Sources on the Epicardial Surface. Mathematical Biosciences, 144, 119-154. https://doi.org/10.1016/S0025-5564(97)00024-2
[7]	Van Oosterom, A. and Strackee, J. (1983) The Solid Angle of a Plane Triangle. IEEE Trans Biomed Eng., BME, 30, 125-126. https://doi.org/10.1109/TBME.1983.325207
[8]	Tikhonov, A.N. and Arsenin, V.Y. (1977) Solutions of Ill-Posed Problems. V.H. Winston & Sons, Washington DC.
[9]	Hansen, P.C. and O’Leary, D.P. (1993) The Use of the l-Curve in the Regularization of Discrete Ill-Posed Problems. Siam Journal on Scientific Computing, 14, 1487-1503. https://doi.org/10.1137/0914086
[10]	Assecondi, S., Hallez, H., D’Asseler, Y. and Lemahieu, I. (2007) Comparison of Different Auto-Solid Angle Approximations in BEM for EEG Dipole Source Localization. Joint Meeting of the International Symposium on Noninvasive Functional Source Imaging of the Brain and Heart and the International Conference on Functional Biomedical Imaging, Nfsi-Icfbi IEEE, 2007, 70-73.
[11]	Yao, B., Pei, S. and Yang, H. (2016) Mesh Resolution Impacts the accuracy of Inverse and Forward ECG Problems. 38th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), Orlando, FL, 4047-4050.
[12]	de Munck, J.C. (1992) A Linear Discretization of the Volume Conductor Boundary Integral Equation Using Analytically Integrated Elements. IEEE Transactions on Biomedical Engineering, 39, 986-990. https://doi.org/10.1109/10.256433
[13]	Friston, K.J. (2006) Statistical Parametric Mapping: The Analysis of Functional Brain Images. Science Press, Beijing.