$\frac{1}{4\pi}{\displaystyle {\int}_{{S}_{B}}^{}{\varphi}_{B}}\frac{\stackrel{\xaf}{r}\cdot \stackrel{\xaf}{n}}{{r}^{2}}d{S}_{B}-{\varphi}_{H}-\frac{1}{4\pi}{\displaystyle {\int}_{{S}_{H}}^{}{\varphi}_{H}\frac{\stackrel{\xaf}{r}\cdot \stackrel{\xaf}{n}}{{r}^{2}}d{S}_{H}}-\frac{1}{4\pi}{\displaystyle {\int}_{{S}_{H}}^{}\frac{\nabla {\varphi}_{H}\cdot \stackrel{\xaf}{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{\stackrel{\xaf}{r}\cdot \stackrel{\xaf}{n}}{{r}^{2}}d{S}_{f}$ . (4)

As it can be seen here, if the direction of $\stackrel{\xaf}{n}$ is wrong, the sign of $\stackrel{\xaf}{r}\cdot \stackrel{\xaf}{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 \stackrel{\xaf}{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 $\stackrel{\xaf}{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 $\stackrel{\xaf}{n}$ must be confirmed before any further calculation.

A new method is proposed to ensure the direction of $\stackrel{\xaf}{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 ${\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 $\stackrel{\xaf}{n}$ in turn. In other words, if applying the theory of Mesh Current method to determine the direction of $\stackrel{\xaf}{n}$ , some principles should be followed: the direction of $\stackrel{\xaf}{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 $\stackrel{\xaf}{n}$ of them are at the same side. This can be seen in Figure 1(b).

In that case, if the direction of $\stackrel{\xaf}{n}$ of a triangle could be determined, the direction of $\stackrel{\xaf}{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 $\stackrel{\xaf}{n}$ is ensured correct and is referred as reference triangle (RT). Check the $\stackrel{\xaf}{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 $\stackrel{\xaf}{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.

Cite this paper

Tang, C. , Wu, J. and Chen, R. (2018) A Method for Normal Direction Judge Applied in Electrocardiographic Problem.*Journal of Biosciences and Medicines*, **6**, 1-8. doi: 10.4236/jbm.2018.61001.

Tang, C. , Wu, J. and Chen, R. (2018) A Method for Normal Direction Judge Applied in Electrocardiographic Problem.

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.