Friday, January 21, 2011

Gaussian integration

Chapter 8 in the 4th edition of Intermediate Physics for Medicine and Biology covers Biomagnetism: the measurement of the magnetic field produced by electrical currents in nerve and muscle. One issue that arises during biomagnetic recordings is that the magnetic field is not measured at a point, but is averaged over a pickup coil. Therefore, when comparing theoretical calculations to experimental data, you need to integrate the calculated magnetic field over the coil.

One way to do this is Gaussian quadrature, which approximates the integral by a weighted sum. Homework problem 40 in Chapter 8 shows a three-point Gaussian quadrature formula for integrating over a circular coil. At the end of the problem Russ Hobbie and I write “Higher-order formulas for averaging the magnetic field can be found in Roth and Sato (1992).” The reference is to “Roth, B. J. and S. Sato (1992) Accurate and efficient formulas for averaging the magnetic field over a circular coil. In M. Hoke, S. N. Erne, T. C. Okada, and G. L. Romani, eds. Biomagnetism: Clinical Aspects. Amsterdam, Elsevier.” This book is actually the proceedings of the 8th International Conference on Biomagnetism, held in Munster, Germany on August 19-24, 1991. I didn’t attend that meeting, but my colleague Susumu Sato did. Sato is a senior scientist in the Epilepsy Research Branch of the National Institute of Neurological Disorders and Stroke, part of the National Institutes of Health in Bethesda, Maryland. When I worked with him he had an active research program in magnetoencephalography (MEG), including a large and expensive shielded room and a multi-channel SQUID magnetometer.

The introduction of our paper states
“The MEG is measured by detecting the magnetic flux through a pickup coil, usually circular, that is coupled to a SQUID magnetometer. Often the source of the MEG is modeled as a current dipole, whose position, orientation and strength are determined iteratively by fitting the MEG data to a dipolar magnetic field pattern. To obtain an accurate result, this dipole field must be integrated over the pickup coil area to obtain the magnetic flux. Since this integration is repeated for each dipole considered in the iteration, the numerical algorithm used to estimate this integral should be efficient. In this note, several integration formulas are presented that allow the magnetic field to be integrated over the coil area quickly with little error. These formulas are examples of a general technique of approximating multiple integrals described by Stroud [1].”
Reference [1] is to: Stroud AH (1971) Approximate Calculation of Multiple Integrals, Prentice-Hall, Englewood Cliffs, New Jersey, pp. 277-289.

I remember deriving several of these formulas independently before discovering Stroud's textbook (it is always deflating to find you’ve been scooped). The derivation requires solving a system of nonlinear equations (which I rather enjoyed). Each formula requires evaluating the magnetic field at N points, and the integral is accurate to mth order. We presented a 1-point formula accurate to first order, a 3-point formula accurate to second order (this was the formula examined in the homework problem), a 4-point formula accurate to third order, a 6-point formula accurate to fourth order, a 7-point formula accurate to fifth order, and a 12-point formula accurate to seventh order.

The general formulation of Gaussian quadrature was developed by Carl Friedrich Gauss (1777-1855), one of the greatest mathematicians of all time. Gauss's name appears often in the 4th edition of Intermediate Physics for Medicine and Biology, including the Gaussian function (Chapter 4), Gauss's law (Chapter 6), the cgs unit for the magnetic field of a gauss (Chapter 8), the Fast Fourier transform (FFT, Chapter 11) about which we write "the grouping used in the FFT dates back to Gauss in the early nineteenth century," and the Gaussian Probability Distribution (Appendix I).

1 comment:

  1. Do you have any experience feeding anatomically correct (as recorded by something like MRI) data into a bidomain solver?

    ReplyDelete