IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38, NO. I , JANUARY 1991
82
Moving Dipole Inverse Solutions Using Realistic Torso Models Christopher J. Purcell and Gerhard Stroink
Abstract-A noniterative numerical solution for the potentials on the surfaces of a piecewise homogeneous volume conductor due to a current dipole is described. This forward solution has been used in electric and magnetic single moving dipole (SMD) inverse solutions that employ a torso volume conductor model whose boundaries are specified numerically. Thus, the volume conductor model used by the inverse solutions need not be limited to simple geometric shapes; torso models of realistic shape can be used.
I. INTRODUCTION
M
OVING dipole inverse solutions can be used to locate a source of electric current in the human body from noninvasive surface recordings of electric potential or magnetic field. Gulrajani et al. [ 11 have reviewed inverse solution methods and applications, which have seen widest use in EEG and MEG studies of evoked responses and epilepsy, and in ECG and MCG studies of WPW syndrome and myocardial infarction. The most popular fitting algorithm used for inverse solutions is that of Marquardt [2]. The algorithm minimizes the differences between the measured surface distributions of electric potential or magnetic field, and distributions calculated using a forward solution. The forward solution requires a source model, a volume conductor model, and a description of the measurement process. The computer time and memory requirements for the inverse solution are dominated by the forward solution since it must be repeated many times. In particular, it is dominated by the volume conductor treatment; so much so, that most inverse solutions described to date have been limited to using volume conductor models with a high degree of symmetry such as concentric spheres or an infinite half space. This is because only in these simple cases are infinite series or closed form expressions for the surface potentials and fields known. He et al. [3] have used a homogeneous head model with realistic outer boundary in an inverse solution to interpret EEG data, but did not describe the forward portion of their inverse algorithm. Our previous modeling studies [4] showed there are circumstances where the use of a semi-infinite volume conductor model in electric and magnetic inverse solutions can lead to errors in the location and orientation of cardiac dipoles. This motivated the development of an inverse solution that uses a volume conductor model with realistically shaped torso boundaries. The key to making this practical in terms of computer resources was the use of the noniterative forward solution described below. Manuscript received February 17, 1989. This work was supported by the Medical Research Council of Canada, the Natural Sciences and Engineering Research Council of Canada, and the Nova Scotia Heart and Stroke Foundation. The authors are with the Departments of Physics, Physiology, and Biophysics, Dalhousie University, Halifax, N.S. B3H 355, Canada. IEEE Log Number 9040626.
NONITERATIVEFORWARD SOLUTION The discretized integral equation for the potentials on the boundaries of a piecewise homogeneous conductor containing a current source can be written as [5]
4
=
4 m + wq5
(1)
where the boundaries have been tessellated with N triangular area elements, q5 is an N dimensional vector of the potentials on all area elements, and q5m is a vector of the mean infinite medium potentials at the same elements. w is an N X N matrix which specifies the conductivity interfaces using the solid angles subtended by every surface area element from every other element. This equation can be solved iteratively using GaussSeidel or Jacobi iteration, after “deflating” the matrix w following the methods of Lynn and Timlake [6]. The magnetic fields can then be calculated by an additional numerical integration [7]. The iterative solution is costly in terms of computer time since it involves more than N 2 floating point multiplications per iteration times m iterations. E.g., we find using N = 320 requires m 10 for satisfactory convergence. Rewriting (1) as
-
(I -
U)-’&
= q5
where I is the identity matrix, does not help since ( I - w ) has no inverse. However, if this same form is used but with the deflated solid angle matrix W [6], we arrive at a useful expression :
(I -
a)-’+-
=
$.
(3)
Here ( I - W ) can be inverted and the solution vector $, which will later have to be “inflated” [6] to recover q5, can then be obtained with N 2 multiplications, once the inverse ( I - 5 ) - ’ is known. Thus, this form yields a forward solution for the potentials that will run m times faster than the iterative solution. If the torso geometry must be changed frequently, the overhead of computing the inverse of ( I - W ) must be considered. This inverse requires of the order of N 3 floating point multiplications [8]. The break even point can be estimated as occurring when N 3 nN2 = nmN2 (4)
-
-
+
where the left-hand side of (4) represents an order of magnitude estimate for the number of multiplications required by the noniterative forward solution, the right-hand side represents an estimate for the number required by the iterative forward solution, and n is the number of forward solutions that are to be found using the same torso geometry. The break even point comes after n = N / m - 1 forward solutions are to be calculated, e.g., by the forward section of an inverse solution.
00 18-9294191/O 100-0082$01.OO 0 1991 IEEE
83
PURCELL A N D STROINK: M O V I N G DIPOLE INVERSE SOLUTIONS
When an electric or magnetic SMD inverse solution based on the Marquardt algorithm is performed using a volume conductor model shaped like a human torso every iteration of the Marquardt algorithm will require -six forward solutions for the potentials. These are needed to calculate derivatives of the norm of the residuals with respect to the six parameters of a single dipole. Derivatives calculated from finite differences must be used for dealing with volume conductors whose boundaries have a complicated shape. This is because closed form or series expressions for the derivatives are not available except for the simple geometries noted above. In our cardiac studies we have found that the Marquardt al15 iterations for convergence using a gorithm often needs 320 triangles. Thus, torso shaped volume conductor with N m 90 forward solutions may be required per inverse solution. These estimates indicate that the time required to perform the matrix inversion will have been recouped after only a single inverse solution.
-
-
-
METHODS To evaluate the noniterative method, the U matrix of a uniform 320 triangle tessellation of a homogeneous sphere with radius 1.O m and conductivity 1.O mho/m was constructed, and the inverse of ( I - W ) calculated. The matrix inversion was done by LU decomposition [8]. The potentials at the triangle centroids due to sets of three orthogonal current dipoles of strength 1 A-m situated at five different radii in the sphere were computed using the iterative forward solution and the non-iterative solution. These were compared with potentials obtained using the “exact” closed form expressions of Wilson and Bayley [9] for the potential at any point in a sphere, evaluated at the triangle centroids.
RESULTS Both methods gave equivalent accuracy, as measured by calculating the average deviations [8] of the differences between the iterated (and noniterated) potentials 4i, ( i = 1, N ), and the “exact” potentials ai. The average deviation (normalized to the peak potential and expressed as a %) is given by 100 ADev = Na,,,,, i =
I
I(&
-
ai)- ( 4 - @)I
(5)
where
(4 - a)
c
l N (4; N ~ = I
= -
ai)and a,,,,,
=
max
{ai}.
The average deviations are shown in Table I for the five radial locations and the three orthogonal directions studied. The average deviations of the iterated and noniterated potentials agreed to four significant figures at all five locations and three orientations studied. Consequently only a single set of entries corresponding to both iterated and non-iterated results are shown in Table I. The average deviations were found to depend on the dipole orientation, hence Table I contains entries for the three orthogonal dipole directions studied. This dependency probably arises from the symmetries of the tessellation of the sphere we employed. The average deviations tend to decrease as the source is moved radially outwards. This is at first glance counter-intuitive, especially when compared with the column of Table I containing the maximum fractional error. The entries in this
TABLE I AVERAGE DEVIATIONS ( ADEv) OF DIFFERENCES OF POTENTIALS AT 320 TRIANGLES CALCULATED BY ITERATIVE(AND NONITERATIVE METHODS) WITH RESPECTTO POTENTIALS CALCULATED USINGEXACT EXPRESSIONS FOR DIPOLES ORIENTED ALONGX, Y , Z DIRECTIONS AND SITUATED ON THE Z AXISOF A HOMOGENEOUS SPHERE OF RADIUS1 .O M AT VARIOUS DISTANCES ( r ) FROM THE CENTER. ALSOGIVENIS THE MAXIMUM FRACTIONAL ERROR (WITH RESPECTTO PEAKPOTENTIAL amax) FOUNDON ANYOF THE 320 TRIANGLES AND FOR ALL THREE DIPOLE ORIENTATIONS.
ADev ( % )
r (m)
X
Y
Z
0.0 0.2 0.4
1.3
0.61 0.60 0.68 0.68 0.55
0.91 0.63 0.42 0.27
0.6 0.8
1.2 1.1 0.91 0.58
1.3
Max Fractional Error ( % ) 3.5 4.5 7.6 15.0
41.0
column are the largest errors (normalized by the peak potential) found on all 320 triangles, and for all three dipole orientations, at each radial position. The maximum fractional errors increase, as the dipoles are moved radially outwards. This is consistent with the expectation that as the source is moved towards the surface, the surface potential gradient will increase, and the approximation implicit in (1), that the potential is constant over each triangle, will begin to break down. The behavior of both these error measures with radial coordinate of the sources can be reconciled. While at a relatively few triangles near the source the calculated potentials may develop substantial error as the source approaches the surface, on average the source is getting further away from the majority of the surface triangles. Hence it is reasonable for the average deviation to decrease. DISCUSSION The iterative and noniterative methods give equivalent accuracy for these five source locations. The noniterativc forward solution, however, runs on average 43 times faster than the iterative solution (not counting the time required to do the matrix inversion). In particular, it takes 78 s of CPU time on a VAX 11/780 to compute the solid angle matrix (required by both methods), 69 s to compute the potentials (at 320 triangles) due to a dipole with the iterative solution, 385 s to do the matrix inversion of ( I - W ) and 1.6 s to compute the potentials using the noniterative solution. These times verify the earlier rough estimates: that it will pay to use the noniterative solution even for a single inverse solution in a 320 triangle torso taking into account our earlier comment that -90 forward solutions will be required per inverse solution. The noniterative forward solution described above has been used in SMD electric and magnetic inverse solutions to fit 117 lead body surface potential and 56 lead magnetocardiographic data [ 101 measured from normal subjects and myocardial infarct patients. These inverse solutions use a homogeneous torso model with a realistic outer boundary. The model consists of 320 triangles, a reduced tessellation of the outer boundary of Horacek’s torso model [ 111. These algorithms can solve the ECG and MCG SMD inverse problems in mean times of 50 and 130 s, respectively, on a VAX 11/780, which makes it practical to perform a substantial number of inverse solutions per subject.
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38. NO. I. JANUARY 1991
84
In conclusion, in situations where the quality of the data justifies using a realistically shaped torso model, the noniterative forward solution is well suited for use in electric and magnetic inverse solutions to locate current dipole sources. ACKNOWLEDGMENT We greatly appreciate the advice and encouragement of Dr. B. Milan Horacek of the Department of Physiology and Biophysics, Dalhousie University. REFERENCES [ I ] R. M. Gulrajani, F. A. Roberge, and P. Savard, IEEE Trans. Biomed. Eng., vol. BME-31, pp. 903-910, Dec. 1984. [2] D. W. Marquardt, J . Soc. Indust. Appl. Math., vol. 11, pp. 431441, 1963. [3] B. He, T. Musha, Y. Okamoto, S. Homma, Y. Nakajima, and T. Sato, IEEE Trans. Biomed. Eng., vol. BME-34, pp. 406-414, June 1987. [4] C. J. Purcell, G. Stroink, and B. M. Horacek, IEEE Trans. Biomed. Eng., vol. BME-35, pp. 671-678, Sept. 1988. [ 5 ] R. C. Barr, T. C. Pilkington, J. P. Boineau, and M. S. Spach, IEEE Trans. Biomed. Eng., vol. BME-13, pp. 88-92, Feb. 1966. [6] M. S. Lynn and W. P. Timlake, SIAM J. Numer. Anal., vol. 5, no. 2, pp. 303-322, 1967. [7] B. M. Horacek, IEEE Trans. Magn., vol. Mag-9, pp. 440-444, 1973. [8] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes: The Art of Scientific Computing. New York: Cambridge Univ., 1987.
[9] F. N. Wilson and R. H. Bayley, Circulation, vol. 1, pp. 84-92, 1950. [IO] G. Stroink, C. Purcell, R. Lamothe, R. Merritt, B. M. Horacek, and B. J. ten Voorde, in Biomagnet. '87, 6th Internat. Conf. Biomagnet., K. Atsumi et a l . , Eds. Tokyo Denki Univ. Press, Tokyo, 1988, pp. 74-81. [ I l l B. M. Horacek, Adv. Cardiol., vol. 10, pp. 51-57, 1974.
Christopher J. Purcell was born in 1951 in Halifax, N.S., Canada. He received the doctorate degree in physics from Dalhousie University, Halifax, in 1988. His thesis work was a comparison of electric and magnetic body surface mapping of cardiac fields. He did postdoctoral research on electroencephalography and magnetoencephalography at the Francis Bitter National Magnet Laboratory, Massachusetts Institute of Technology during 1988-1990, and was an invited researcher of the NTT Basic Research Laboratory, Musashino-Shi, Tokyo, during 1989. He is presently with the Defence Research Establishment Atlantic, Dartmouth, N.S., doing research on finite element modeling of acoustic transducers.
Gerhard Stroink was born in 1942. He received the engineering physics degree from the Technical University of Delft, the Netherlands, and the Ph.D. degree in physics from McGill University, Montreal, P.Q., Canada. He joined the Department of Physics, Dalhousie University, Halifax, in 1973 and is presently a Professor in that department. He is also an Associate Professor in the Department of Physiology and Biophysics. His research interests include electro- and magnetocardiology, and the magnetic properties of amorphous alloys.
82
Moving Dipole Inverse Solutions Using Realistic Torso Models Christopher J. Purcell and Gerhard Stroink
Abstract-A noniterative numerical solution for the potentials on the surfaces of a piecewise homogeneous volume conductor due to a current dipole is described. This forward solution has been used in electric and magnetic single moving dipole (SMD) inverse solutions that employ a torso volume conductor model whose boundaries are specified numerically. Thus, the volume conductor model used by the inverse solutions need not be limited to simple geometric shapes; torso models of realistic shape can be used.
I. INTRODUCTION
M
OVING dipole inverse solutions can be used to locate a source of electric current in the human body from noninvasive surface recordings of electric potential or magnetic field. Gulrajani et al. [ 11 have reviewed inverse solution methods and applications, which have seen widest use in EEG and MEG studies of evoked responses and epilepsy, and in ECG and MCG studies of WPW syndrome and myocardial infarction. The most popular fitting algorithm used for inverse solutions is that of Marquardt [2]. The algorithm minimizes the differences between the measured surface distributions of electric potential or magnetic field, and distributions calculated using a forward solution. The forward solution requires a source model, a volume conductor model, and a description of the measurement process. The computer time and memory requirements for the inverse solution are dominated by the forward solution since it must be repeated many times. In particular, it is dominated by the volume conductor treatment; so much so, that most inverse solutions described to date have been limited to using volume conductor models with a high degree of symmetry such as concentric spheres or an infinite half space. This is because only in these simple cases are infinite series or closed form expressions for the surface potentials and fields known. He et al. [3] have used a homogeneous head model with realistic outer boundary in an inverse solution to interpret EEG data, but did not describe the forward portion of their inverse algorithm. Our previous modeling studies [4] showed there are circumstances where the use of a semi-infinite volume conductor model in electric and magnetic inverse solutions can lead to errors in the location and orientation of cardiac dipoles. This motivated the development of an inverse solution that uses a volume conductor model with realistically shaped torso boundaries. The key to making this practical in terms of computer resources was the use of the noniterative forward solution described below. Manuscript received February 17, 1989. This work was supported by the Medical Research Council of Canada, the Natural Sciences and Engineering Research Council of Canada, and the Nova Scotia Heart and Stroke Foundation. The authors are with the Departments of Physics, Physiology, and Biophysics, Dalhousie University, Halifax, N.S. B3H 355, Canada. IEEE Log Number 9040626.
NONITERATIVEFORWARD SOLUTION The discretized integral equation for the potentials on the boundaries of a piecewise homogeneous conductor containing a current source can be written as [5]
4
=
4 m + wq5
(1)
where the boundaries have been tessellated with N triangular area elements, q5 is an N dimensional vector of the potentials on all area elements, and q5m is a vector of the mean infinite medium potentials at the same elements. w is an N X N matrix which specifies the conductivity interfaces using the solid angles subtended by every surface area element from every other element. This equation can be solved iteratively using GaussSeidel or Jacobi iteration, after “deflating” the matrix w following the methods of Lynn and Timlake [6]. The magnetic fields can then be calculated by an additional numerical integration [7]. The iterative solution is costly in terms of computer time since it involves more than N 2 floating point multiplications per iteration times m iterations. E.g., we find using N = 320 requires m 10 for satisfactory convergence. Rewriting (1) as
-
(I -
U)-’&
= q5
where I is the identity matrix, does not help since ( I - w ) has no inverse. However, if this same form is used but with the deflated solid angle matrix W [6], we arrive at a useful expression :
(I -
a)-’+-
=
$.
(3)
Here ( I - W ) can be inverted and the solution vector $, which will later have to be “inflated” [6] to recover q5, can then be obtained with N 2 multiplications, once the inverse ( I - 5 ) - ’ is known. Thus, this form yields a forward solution for the potentials that will run m times faster than the iterative solution. If the torso geometry must be changed frequently, the overhead of computing the inverse of ( I - W ) must be considered. This inverse requires of the order of N 3 floating point multiplications [8]. The break even point can be estimated as occurring when N 3 nN2 = nmN2 (4)
-
-
+
where the left-hand side of (4) represents an order of magnitude estimate for the number of multiplications required by the noniterative forward solution, the right-hand side represents an estimate for the number required by the iterative forward solution, and n is the number of forward solutions that are to be found using the same torso geometry. The break even point comes after n = N / m - 1 forward solutions are to be calculated, e.g., by the forward section of an inverse solution.
00 18-9294191/O 100-0082$01.OO 0 1991 IEEE
83
PURCELL A N D STROINK: M O V I N G DIPOLE INVERSE SOLUTIONS
When an electric or magnetic SMD inverse solution based on the Marquardt algorithm is performed using a volume conductor model shaped like a human torso every iteration of the Marquardt algorithm will require -six forward solutions for the potentials. These are needed to calculate derivatives of the norm of the residuals with respect to the six parameters of a single dipole. Derivatives calculated from finite differences must be used for dealing with volume conductors whose boundaries have a complicated shape. This is because closed form or series expressions for the derivatives are not available except for the simple geometries noted above. In our cardiac studies we have found that the Marquardt al15 iterations for convergence using a gorithm often needs 320 triangles. Thus, torso shaped volume conductor with N m 90 forward solutions may be required per inverse solution. These estimates indicate that the time required to perform the matrix inversion will have been recouped after only a single inverse solution.
-
-
-
METHODS To evaluate the noniterative method, the U matrix of a uniform 320 triangle tessellation of a homogeneous sphere with radius 1.O m and conductivity 1.O mho/m was constructed, and the inverse of ( I - W ) calculated. The matrix inversion was done by LU decomposition [8]. The potentials at the triangle centroids due to sets of three orthogonal current dipoles of strength 1 A-m situated at five different radii in the sphere were computed using the iterative forward solution and the non-iterative solution. These were compared with potentials obtained using the “exact” closed form expressions of Wilson and Bayley [9] for the potential at any point in a sphere, evaluated at the triangle centroids.
RESULTS Both methods gave equivalent accuracy, as measured by calculating the average deviations [8] of the differences between the iterated (and noniterated) potentials 4i, ( i = 1, N ), and the “exact” potentials ai. The average deviation (normalized to the peak potential and expressed as a %) is given by 100 ADev = Na,,,,, i =
I
I(&
-
ai)- ( 4 - @)I
(5)
where
(4 - a)
c
l N (4; N ~ = I
= -
ai)and a,,,,,
=
max
{ai}.
The average deviations are shown in Table I for the five radial locations and the three orthogonal directions studied. The average deviations of the iterated and noniterated potentials agreed to four significant figures at all five locations and three orientations studied. Consequently only a single set of entries corresponding to both iterated and non-iterated results are shown in Table I. The average deviations were found to depend on the dipole orientation, hence Table I contains entries for the three orthogonal dipole directions studied. This dependency probably arises from the symmetries of the tessellation of the sphere we employed. The average deviations tend to decrease as the source is moved radially outwards. This is at first glance counter-intuitive, especially when compared with the column of Table I containing the maximum fractional error. The entries in this
TABLE I AVERAGE DEVIATIONS ( ADEv) OF DIFFERENCES OF POTENTIALS AT 320 TRIANGLES CALCULATED BY ITERATIVE(AND NONITERATIVE METHODS) WITH RESPECTTO POTENTIALS CALCULATED USINGEXACT EXPRESSIONS FOR DIPOLES ORIENTED ALONGX, Y , Z DIRECTIONS AND SITUATED ON THE Z AXISOF A HOMOGENEOUS SPHERE OF RADIUS1 .O M AT VARIOUS DISTANCES ( r ) FROM THE CENTER. ALSOGIVENIS THE MAXIMUM FRACTIONAL ERROR (WITH RESPECTTO PEAKPOTENTIAL amax) FOUNDON ANYOF THE 320 TRIANGLES AND FOR ALL THREE DIPOLE ORIENTATIONS.
ADev ( % )
r (m)
X
Y
Z
0.0 0.2 0.4
1.3
0.61 0.60 0.68 0.68 0.55
0.91 0.63 0.42 0.27
0.6 0.8
1.2 1.1 0.91 0.58
1.3
Max Fractional Error ( % ) 3.5 4.5 7.6 15.0
41.0
column are the largest errors (normalized by the peak potential) found on all 320 triangles, and for all three dipole orientations, at each radial position. The maximum fractional errors increase, as the dipoles are moved radially outwards. This is consistent with the expectation that as the source is moved towards the surface, the surface potential gradient will increase, and the approximation implicit in (1), that the potential is constant over each triangle, will begin to break down. The behavior of both these error measures with radial coordinate of the sources can be reconciled. While at a relatively few triangles near the source the calculated potentials may develop substantial error as the source approaches the surface, on average the source is getting further away from the majority of the surface triangles. Hence it is reasonable for the average deviation to decrease. DISCUSSION The iterative and noniterative methods give equivalent accuracy for these five source locations. The noniterativc forward solution, however, runs on average 43 times faster than the iterative solution (not counting the time required to do the matrix inversion). In particular, it takes 78 s of CPU time on a VAX 11/780 to compute the solid angle matrix (required by both methods), 69 s to compute the potentials (at 320 triangles) due to a dipole with the iterative solution, 385 s to do the matrix inversion of ( I - W ) and 1.6 s to compute the potentials using the noniterative solution. These times verify the earlier rough estimates: that it will pay to use the noniterative solution even for a single inverse solution in a 320 triangle torso taking into account our earlier comment that -90 forward solutions will be required per inverse solution. The noniterative forward solution described above has been used in SMD electric and magnetic inverse solutions to fit 117 lead body surface potential and 56 lead magnetocardiographic data [ 101 measured from normal subjects and myocardial infarct patients. These inverse solutions use a homogeneous torso model with a realistic outer boundary. The model consists of 320 triangles, a reduced tessellation of the outer boundary of Horacek’s torso model [ 111. These algorithms can solve the ECG and MCG SMD inverse problems in mean times of 50 and 130 s, respectively, on a VAX 11/780, which makes it practical to perform a substantial number of inverse solutions per subject.
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38. NO. I. JANUARY 1991
84
In conclusion, in situations where the quality of the data justifies using a realistically shaped torso model, the noniterative forward solution is well suited for use in electric and magnetic inverse solutions to locate current dipole sources. ACKNOWLEDGMENT We greatly appreciate the advice and encouragement of Dr. B. Milan Horacek of the Department of Physiology and Biophysics, Dalhousie University. REFERENCES [ I ] R. M. Gulrajani, F. A. Roberge, and P. Savard, IEEE Trans. Biomed. Eng., vol. BME-31, pp. 903-910, Dec. 1984. [2] D. W. Marquardt, J . Soc. Indust. Appl. Math., vol. 11, pp. 431441, 1963. [3] B. He, T. Musha, Y. Okamoto, S. Homma, Y. Nakajima, and T. Sato, IEEE Trans. Biomed. Eng., vol. BME-34, pp. 406-414, June 1987. [4] C. J. Purcell, G. Stroink, and B. M. Horacek, IEEE Trans. Biomed. Eng., vol. BME-35, pp. 671-678, Sept. 1988. [ 5 ] R. C. Barr, T. C. Pilkington, J. P. Boineau, and M. S. Spach, IEEE Trans. Biomed. Eng., vol. BME-13, pp. 88-92, Feb. 1966. [6] M. S. Lynn and W. P. Timlake, SIAM J. Numer. Anal., vol. 5, no. 2, pp. 303-322, 1967. [7] B. M. Horacek, IEEE Trans. Magn., vol. Mag-9, pp. 440-444, 1973. [8] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes: The Art of Scientific Computing. New York: Cambridge Univ., 1987.
[9] F. N. Wilson and R. H. Bayley, Circulation, vol. 1, pp. 84-92, 1950. [IO] G. Stroink, C. Purcell, R. Lamothe, R. Merritt, B. M. Horacek, and B. J. ten Voorde, in Biomagnet. '87, 6th Internat. Conf. Biomagnet., K. Atsumi et a l . , Eds. Tokyo Denki Univ. Press, Tokyo, 1988, pp. 74-81. [ I l l B. M. Horacek, Adv. Cardiol., vol. 10, pp. 51-57, 1974.
Christopher J. Purcell was born in 1951 in Halifax, N.S., Canada. He received the doctorate degree in physics from Dalhousie University, Halifax, in 1988. His thesis work was a comparison of electric and magnetic body surface mapping of cardiac fields. He did postdoctoral research on electroencephalography and magnetoencephalography at the Francis Bitter National Magnet Laboratory, Massachusetts Institute of Technology during 1988-1990, and was an invited researcher of the NTT Basic Research Laboratory, Musashino-Shi, Tokyo, during 1989. He is presently with the Defence Research Establishment Atlantic, Dartmouth, N.S., doing research on finite element modeling of acoustic transducers.
Gerhard Stroink was born in 1942. He received the engineering physics degree from the Technical University of Delft, the Netherlands, and the Ph.D. degree in physics from McGill University, Montreal, P.Q., Canada. He joined the Department of Physics, Dalhousie University, Halifax, in 1973 and is presently a Professor in that department. He is also an Associate Professor in the Department of Physiology and Biophysics. His research interests include electro- and magnetocardiology, and the magnetic properties of amorphous alloys.
