WO1993020749A1 - Eeg spatial enhancement method and system - Google Patents

Abstract

In an electrographic (EEG) system and method the subject's brain waves are recorded using electrodes removably connected to the scalp. The spatial resolution may be improved, in one embodiment, by 'Deblurring', based on digital computer processing measurements of the positions of the electrodes; measurements of the subject's head size, shape, scalp thickness, skull thickness and brain shape; and estimates of the conductivity of the skull, scalp and cerebral spinal fluid. The system includes head support (10), adjustable guides (11), probe (12), and computer (13).

Description

"EEG SPATIAL ENHANCEMENT METHOD AND SYSTEM"

This invention was made with United States Government under Grants R44-MH-27525, R01-MH43324 and R01-NS23550 awarded by the Drug Abuse and Mental Health Administration and the National Institutes of Health. The Government has certain rights in the invention.

Background of the Invention

1. Field of Invention

The invention relates to medical systems and methods and more particularly to an electroencephalograph (EEG) system and method having improved spatial resolution.

2. Related Art

Electroencephalography (EEG) is a noninvasive technique for monitoring brain function based on amplification of brain potentials recorded at the scalp. Although the EEG has been measured for over 60 years, its full potential as a technology for imaging brain function has not yet been realized. Its chief advantages, when compared with other brain imaging technologies, such as positron emission tomography (PET) are: (1) an

unsurpassed millisecond-level temporal resolution necessary for resolving sub-second neurological processes; and (2) simplicity of the sensors (essentially a small piece of metal and an operational amplifier) which results in low cost and ability to record for prolonged periods of time in a naturalistic setting. The chief disadvantage of the traditional EEG method is poor spatial resolution. There are only two obstacles to

substantially improving spacial resolution. The first is that more scalp recording sites are needed. While it is the usual clinical practice to record EEG from 19 sites, it has been demonstrated with recordings from up to 124 sites that additional spatial information is available. U.S. Patent 4,736,751,

incorporated by reference, describes a system using a larger number of electrodes and various digital computer based methods to obtain more information from the brain wave signals. U.S.

patents 4,967,038, 5,038,782 and [patent to issue on Serial No. 07/578,880], incorporated by reference, describe systems for conveniently obtaining EEG recordings from a larger number of electrodes and locating their positions. The second obstacle is that spatial enhancement procedures are need to reduce blur distortion that occurs as electrical potentials generated in the brain at the micro-volt level are volume conducted through brain, cerebral spinal fluid, the low-conductivity skull, and scalp to the recording electrodes at the scalp surface. However, spatial enhancement is possible since the 3 dB point of the point spread function for conductance of potentials from the brain surface to the scalp averages about 2.5 cm and with 128 electrodes spaced evenly on an average adult head, the interelectrode distance is about 2.25 cm. U.S. Patent 4,416,288 to Walter Freeman entitled "Apparatus And Method For Reconstructing Subsurface Electrophysiological Patterns" discusses, in an EEG system, calculation of spatial deconvolution transformation of the detected electrode field potentials. However, no attempt is made to actually measure the shape and thickness of the scalp and skull of the subject, or the actual positions of the electrodes on the scalp, or to correct the potentials for distortion due to transmission through the skull and scalp.

Summary of the Invention

A functional-anatomical brain scanner with a temporal resolution of less than a hundred milliseconds measures the neural substrate of higher cognitive functions, as well as diagnosing seizure disorders. Electrophysiological techniques, such as electroencephalograms (EEGs) , have the requite temporal resolution but their potential spatial resolution has not yet been realized. Progress in increasing the spatial detail of scalp-recorded EEGs and in registering their functional

information with anatomical models of a patient's brain has been inhibited by the lack of a convenient means of rapidly placing many electrodes on the scalp, determining their location,

determining the local thickness and conductivity of the scalp and skull, incorporating this information into a mathematical model, and using the mathematical model to enhance the spatial detail of the EEG signals. In the present invention, the three-dimensional positions of each electrode and the shape of the head are

determined, the local thickness and conductivity of the skull and scalp are determined, and this information is used to reduce blur distortion of scalp-recorded EEGs, in effect mathematically placing the electrode just above the surface of the brain.

In accordance with one aspect of the present invention, a method is presented to improve the spatial resolution of

electroencephalograph (EEG) medical images. The physical location of each of the electrodes in a large set of electrodes placed on the head of a subject is measured. A flexible hat is placed on the subject's head and piezoelectric ribbons, embedded in the hat, are stretched. The ribbons are directed in cross-lines and the extent the ribbons are stretched (tensioned) determines the size of the subject's head in the direction along which the ribbon lies. The head measurements are entered into a pattern recognition computer system, preferably a neural network computer, and the subject's head shape is classified as falling into one class of head shapes. The electrode locations are mathematically derived in a computer system based on the known locations of the corresponding electrodes on a head model in the same head shape class and the size of the subject's head.

In accordance with the present invention, a method is presented to improve the spatial resolution of

electroencephalograph (EEG) medical images using a realistic biophysical model of the passive conducting properties of each subject's head to estimate an electrical potential distribution at the cortical surface. A procedure, called "Deblurring", reduces spatial blur distortion of scalp-recorded EEGs due to transmission through the skull and other tissues. Deblurring estimates potentials at the superficial cerebral cortical surface from EEGs recorded at the scalp using a finite element or surface integral model of each subject's scalp, skull and cortical surface constructed from their Magnetic Resonance Images (MRIs) or from other measurement means such as mechanical measurement of the shape of the outside of the head and ultrasonic measurement of the thickness of the scalp and skull. The Deblurring

procedure requires: (1) measuring the size and shape of the subject's head; (2) measuring the local scalp and skull

thickness; (3) estimating the local conductivity of the skull; and (4) using these measurements in a mathematical model of the head to correct the EEG signals for distortion.

Brief Description of the Drawings

In the drawings:

Figures la-lc are diagrams of a 124-channel scalp montage. Electrode names are based on 10-20 system with additional letter prefixes and numerical suffices. Figure la shows the top view of the head, Figure lb shows the right side of the head, and Figure lc shows the left side of the head. Figure 2 is a perspective view of the apparatus used to measure the electrode coordinates. The research assistant touches each electrode with a magnetic position sensor while the subject is resting his head on a chin rest built into a head support which gently restricts head movement during the

measurements. The 3-D coordinates of each electrode position are transmitted to the data collection computer.

Figure 3 shows the top views of seven representative head shapes.

Figure 4 shows output of automatic edge-detection software of a Midline sagittal Magnetic Resonance Image (MRI) obtained with TR=600 and TR=20 ms. The lines shown correspond to (Right) Output of automatic edge-detection software. The lines shown corresponding to locations of local maxima in the gradient of the image intensity. These contours pass through the points at which the local average intensity is changing most rapidly.

Figure 5 shows automated construction of finite elements within brain, skull and scalp volumes, with terrahedral elements generated in an outer ring in brain.

Figure 6 is a perspective view of the flexible measurement hat having piezo-electric ribbons embedded therein.

Figures 7a and 7b are top plan views of a cross-section of a human head. Figure 7a shows the original and Figure lb shows the Laplacian Derivation (LD) 124-channel evoked potentials (EPs) elicited by 15-Hz stimulation of the middle and right index fingers. The data are mapped onto a computer model of the subject's head derived from his magnetic resonance image (MRI) and show the improvement in spatial detail produced by the LD.

Figures 8a and 8b show a single horizontal magnetic

resonance image (MRI), at Figure 8a, after the outer surface of the scalp, skull and brain have been contoured (small circles). (The rings of circles on the inner surface of the brain do not attempt to outline sulci and gyri.) The resulting finite

elements constructed from these contours are shown at Figure 2b with alternating elements shown for brain, skull and scalp volumes, respectively.

Figures 9a-9d show automated construction of finite elements within brain, skull and scalp volumes, with alternating elements shown for the respective tissue type. Figure 9a is a single ring of scalp. Figure 9b shows the addition of skull and brain elements within one slice bounded by two consecutive horizontal MR images. Figure 9c shows all slices except for the topmost horizontal slice. Figure 9d shows a shaded depiction of all three surfaces superimposed.

Figure 10 shows a block schematic diagram of the Deblurring method. The initial estimate of the cortical potentials,

G (x, y, z) is the starting point for the loop contained inside the large frame. Using the estimated cortical potentials and the volumetric finite elements, the finite element method generates a forward solution of Maxwell's equation and produces an estimate of the scalp potentials. The goodness of fit of the estimated scalp potentials with respect to the actual scalp data is

evaluated. If the goodness of fit is acceptable, the method has converged and ends. If it is not acceptable, the

multidimensional optimization scheme searches for a new estimate of the cortical potentials, and the loop repeats until

convergence is achieved. The round corner boxes represent the procedures and the rectangular boxes represent the intermediate data sets. The boxes outside of the large frame are the fixed data sets for the loop.

Figures Hal, Ha2, llbl and Hb2 are perspective views of a patient showing the top of the brain exposed and showing 15-Hz steady-state somatosensory evoked potentials (EPs) elicited by stimulation of right-hand fingers. The patient was an epileptic patient who was a candidate for surgical treatment of her

seizures. (a) Isopotential map of the original EP on the scalp surface (Fig. Hal) and the deblurred EP computed on the cortical surface (Fig. Ha2). The brain, skull and scalp surfaces were constructed from horizontal magnetic resonance images obtained prior to surgery, one of which is shown in part between the scalp and cortex (Fig. Hal). The small cylinders on the scalp surface represent the recording electrodes. The single large potential maxima in the scalp topography is more clearly defined in the deblurred data, which in addition shows some polarity reversals not visible in the scalp potential distribution. (b) The EP recorded directly on the brain surface is shown (Fig. Hb2) with the small cylinders representing the electrodes of the recording grid. The deblurred EP is shown again (Fig. llbl) for

comparison. The area covered by the grid shows a single large peak similar to that in the deblurred potentials. (c) Figures llcl - Hc3 show the effect of varying the scalp/skull

conductivity ratio (r) on the deblurred data of Figures 5bl. The standard (80), and ratios which are +/- 50% of the standard (40 and 120) are shown. There is not much difference in topography between t=40 and t=80. The topography appears more diffuse and the fine detail changes somewhat for t=120, which models the skull with the smallest conductivity value.

Figures 12al, 12a2, 12a3, 12bl and 12b2 are top plan views of cross-sections of a head. Figures 12al - 12a3 show deblurred steady-state evoked potentials (EPs) elicited by 15-Hz simulation of left middle, right middle and right index fingers of a normal subject. All three cases show the expected contralateral maximum in activity, while the response to right index finger stimulation differs markedly from that for the right middle finger. Figs. 12bl and 12b2 show original (Fig. 12bl) and deblurred (Fig. 12b2) EPs of the same subject in response to 15-Hz stimulation of left middle and right index fingers. The white dots (Fig. 12bl) represent the 122 scalp recording positions. The increase in spatial detail produced by the Deblurring procedure is apparent. The improvement in spatial detail as compared with the Laplacian Derivation of the same data (Figure 7b) is also evident.

Detailed Description

I. EEG RECORDING METHOD

A. Nomenclature for 256 Electrode Placements

A nomenclature is defined for the extra equidistant coronal rows added between the original rows of the conventional

International 10-20 System as well as the extra equidistant electrodes added to fill in the spaces in each row (Figures 1a- 1c). For example, the preferred number of electrodes is 64-512 and most preferred are 256 scalp electrodes. The letter "a" is prefixed to the conventional 10-20 row name to indicate a

position anterior to the existing one_* For example, aO is anterior occipital, aP is anterior parietal, aC is anterior central and aF is anterior frontal. Two additional numbers appear after the location number (1-8) of the 10-20 positions. These two numbers (0-9) indicate the proportional distance to the next anterior and medical 10% electrode position, respectively. For example, an electrode halfway between P3 and aPl would be named P355, while an electrode halfway between P3 and P1 would be named P305. With the original nineteen electrodes of the 10-20 System, the typical distance between electrodes on an average adult male head is about 6 cm, with 124 electrodes, the typical distance is about 2.25 cm; with 256 electrodes, the typical distance is about 1 cm.

B. Smart EEG Hat

In most routine clinical EEGs, electrodes are prepared and positioned on the head individually, while in some labs

commercially available electrode caps with up to 32 built-in electrodes are used. Extending this idea, we developed a more efficient system for recording up to 124 or more EEG channels with little or no preparation of the scalp, as described in U.S. Patent 5,038,782, incorporated by reference herein.

C. Measuring Electrode Positions

Traditionally, individual electrodes are placed according to measurements taken with a tape measure. While this has been sufficient for qualitative examination of strip chart tracings with up to 21 electrodes, greater precision and a method using less time are needed when many more electrodes are used and when one wishes to relate a recording position on the scalp to the underlying brain anatomy. In accordance with the present invention, the position of each electrode on a subject's head is measured with a probe 12 such as the 3-Space Isotrak digitizer by Polhemus Navigation Sciences, that has coils for sensing the three-dimensional position of the probe tip with respect to a magnetic field source attached to the head support 10 (Figure 2). Adjustable guides 11 built into the head support hold the subject's head comfortably in place while the measurements are made. A menu-driven program on a control digital computer 13 is used to select electrodes to be measured and display the digitized position of each electrode on a two-dimensional projection display. Position measurement is accurate to better than 3mm [RMS] . While this procedure provides greatly increased precision over the use of a tape, it can take fifteen minutes to measure the positions of 124 electrodes during which time the patient must keep his head still in the head support device. Such lack of motion for that time period may be difficult or impossible for very young or infirm patients.

An alternative procedure is to attach the transmitter for the Isotrak digitizer to the head with an adjustable elastic band, eliminating the need for the head support device. This alternative is more comfortable for the patient, but does not reduce the time required to measure each electrode.

In accordance with the present invention, in another

alternative procedure, electrode positions are determined by classifying a subject's head shape as one of seven canonical head shapes (Figure 3), determining a scale factor which is the size of the subject's head relative to an average sized head of that canonical shape, and estimating electrode positions as scaled values of the standard electrode positions for that particular canonical head. The patient's head shape is classified by:

(1) measuring the stretch of stretch hat 20 having elastic bands

21 at several points over the head using piezo-electric sensors

22 or ribbons embedded within the heat (see Figure 6);

(2) scaling the vector of stretch measurements, termed the

"stretch vector", to have unit magnitude; and (3) feeding this vector to a neural network or other mathematical computer pattern classifier that has been trained to make the desired

classification The scaling factor is the ratio between the stretch vector magnitudes of the subject and that determined from a canonical head of average size. For each electrode montage provided with the hat and each canonical head shape, accurate measurements of electrode position are made on an average sized model of that head. Appropriately scaling these electrode positions gives a reasonably accurate estimate of actual

electrode position and has the advantage that no additional time is added to the EEG recording session.

Another alternative procedure is to place a flexible grid of electric field sources in close proximity to the head and then measure the potential induced in each electrode when each node in the grid is activated. For example the nodes (electric field sources) may he resonant circuits or switch points which form a grid on a flexible cap. The relative position of each node is calculated by interpolation from the pattern of induced voltages. The grid may be embedded in a hat which fits over the hat containing the recording electrodes or the recording electrodes and field source may be combined This has the advantage of increased precision with no significant increase in recording time.

D. Software System for EEG Analysis

Although commercial EEG computer systems have improved during the past few years, some important limitations still restrict their utility, including: (1) system capacity that is inadequate for recording large numbers of channels (124 or more) and collecting and storing large databases; and (2) spatial sampling and analyses that are inadequate, with too few channels, a lack of effective means to reduce volume conduction distortion, a lack of cross-channel analyses (e.g., crosscovariance and crosspower, correlation and coherence), and a lack of means to investigate the relationship between neuroeiectric data and cortical anatomy and physiology revealed by MRI and PET imaging technologies.

1. Data Collection:

One or two computers are used for stimulus presentation and data collection. One computer is used to present stimuli and gather behavioral response data from the subject, while the other collects physiological data and controls the first computer. When a very fast computer is used, or not many channels are required, a single computer can be used to per- form all operations. A control software system runs the two computers, presents a variety of visual, auditory and somatosensory stimuli to the subject according to flexible task protocols, and digitizes up to 256 channels of evoked brain wave physiological data. Up to 128 EEG traces can be monitored in real time as data are collected. Most parameters of an experiment can be altered via a menu-driven interactive display. A calibration module numerically adjusts the gain of all channels according to the magnitude of a calibration signal which is injected, under computer control, into each electrode. Another module detects gross artifacts and color-codes contaminated data such as eye movements, gross head and body movements and bad electrode contacts. Stimulus, behavioral and physiological data are stored on hard disk according to a self-decoding data description language, and archived to optical disk or magnetic tape. The data are immediately available to researchers at their desks using either remote terminals, or UNIX workstations via Ethernet and Network File System.

2. Data Analysis:

For subsequent data analysis the system also has a number of other functional improvements over its predecessors including: (1) least-squares and 3-D spline Laplacian derivation estimation and spatial deblurring using finite element or integral methods to reduce volume conduction distortion; (2) digital filters with user-specified characteristics; (3) time series analysis inclviding spectral analysis, Wagner Distributions and event-related covariance analysis; ( 4) neural-network-driven pattern recognition to extract optiπal or near-optimal subsets of features for recognizing different experimental clinical categories; and ( 5) anatomical modeling to construct 3-D mathematical models of the brain and head from MRI or CT brain scans. Four on-line, interactive subsystems are used to examine and edit data for residual artifacts (on an individual channel basis if desired), sort data according to stimulus, response or other categories, perform exploratory data analysis, and produce three-dimensional color graphics representations of the brain and head. II. BRAIN MODELING METHODS

Since commercial magnetic resonance image(MRT) analysis packages lack essential features for functional-anatomical integration studies, there is described below algorithms and a software outline to produce 3-D brain models suitable for functional localization studies. Visualization software permits construction of 3-D composites of multiple 2-D image planes, as well as 3-D surface rendering based on surface contours. Since generating surface contours manually is laborious and subject to error, we have largely automated the procedure. We have also automated the alignment of the digitized EEG electrode positions with the scalp surface contours, which is a critical first step in a functional-anatomical analysis. The MRI contour information is used to produce mathematical finite element or surface integral models of the brain and head suitable for equivalent dipole source localization and scalp EEG deblurring procedures.

A. EEG-MRI Alignment Procedure

In order to visualize the brain areas underlying EEG electrodes, a procedure is needed for aligning scalp electrode positions with the MR Images. In our procedure, x, y, z translation and x, y, z axis rotations are computed iteratively to align the digitized positions of the EEG electrodes with the MRI data. This is done for each electrode by finding the distance to the closest point on the scalp surface MRI contours and minimizing the mean distance for all electrodes. With MR images that have a 3 mm inter-slice spacing, the mean error distance is usually better than 2 mm. This is more accurate and less subjective than alignment procedures that use skull landmarks such as the nasion and pre-auricular points located visually in the MR images. Fig. 1a shows the electrodes displayed schematically on a scalp surface reconstructed from horizontal MR scalp contours as described below.

B. 3-D Composite MRI Displays

In order to visualize how the individual MRI slices fit together to form a threedimensional object, a means is needed to juxtapose an arbitrary number of slices from sagital, coronal and horizontal scanning planes into a composite image, and view the jextaposed slices from an arbitrary viewing angle.

C. Contour Extraction and External Surface Rendering

In order to visualize the external three dimensional surface of the brain and head, a means is needed to compute the curved surfaces from the information contained in the set of two-dimensional MRI slices. Scalp and cortical surface contours are extracted from MRIs using two image analysis methods we developed. The first technique uses intensity thresholding, which involves extracting contours along which the image intensity is equal to a defined threshold. This technique is useful for extracting the scalp surface contour which can be discerned easily from the black (approximately zero intensity) background of the image. The technique has six steps: 1) for each pair of adjacent contours, find the point on the second contour closest to the first; 2) calculate the two distances from each of these points to the next point on the adjacent contour; 3) make a triangle using the point with the shortest distance, advancing on that contour; 4) repeat steps two and three until all points are exhausted; 5) repeat steps two to four in the reverse direction along the contours; and 6) piece together the "best" result

The second contour extraction method involves differential intensity analysis. Using this technique, contours that separate image regions with different local average intensities are extracted Resulting contours pass through the points in the image at which the local average image intensity is changing most rapidly. This technique requires no a priori intensity threshold value, and is useful for extracting the cortical surface contour which has a less well-defined image intensity value throughout an image than the scalp surface contour. The first- order and second-order partial derivatives of the image are estimated using 2-D filters, and these derivatives are used to locate the local maxima in the gradient of the image intensity. Highly computationally efficient filtering techniques have been developed for the estimation of the partial derivatives of the image. Figure 6 shows an example of contours corresponding to locations of local maxima in the gradient of the image intensity, which trace sulci.

D. 3-D Cortical Internal Surface Image Model 3-D surface models of the external convexity of the cerebral cortex are not sufficient to visualize and computationally model the cortical surface within fissures and sulci. As a result, we developed an algorithm to model cubic volume elements (voxels). The faces of the voxels lie in the horizontal, coronal, and sagittal planes for which MRI data have been obtained. The MRI data are mapped onto the faces of the voxels to obtain a 3-dimensional display of the image data. The image planes used are averaged horizontal planes that lie halfway between the horizontal planes in which MRI data have been acquired, and coronal and sagittal image planes (that can optionally be synthesized from the acquired horizontal images). The spacing between image planes is 3 mm, which yields voxels with dimensions of 3 mm by 3 mm by 3 mm. The images have pixel dimensions of approximately 1 mm by 1 mm.

The initial set of voxels is the set bounded by those voxels that lie just inside the cortical surface. To view the MRI data at a slightly deeper level, a mathematical morphology erosion operation is used to remove the boundary layer of voxels, thereby exposing the faces of the voxels that lie one layer deeper. By eroding the model iteratively, anatomical structures can be tracked and 3-dimensional models of the structures can be made.

E. 3-D Electrical Brain Model: Finite Element Method

A mathematical 3-D representation of the brain and head are needed to most effectively deblur the EEG recorded at the scalp and to compute equivalent current dipole sources from the scalp EEG. We developed an efficient finite element method for this purpose.

Maxwell's Equation, ∇•( σ∇u ) + V•J = 0 in Ω (1) is frequently used to study the electromagnetic field generated by populations of neurons. Here σ (> 0) is the conductivity tensor at a point R = (x,y,z) in Ω (e.g., a human head), u is the electric potential at R , J is the eiectric current density at R and∇ is the grad vector operator. When Ω represents a human head, the following boundary condition for equation (1) is obtained σ∇u•n = 0 on S(headsurface) (2) since electric current does not flow out of the head in the direction normal to the surfaces. The limited scalp potential measurement specifies another boundary condition. u - U(x, y, z ) on S₁. (3)

Where S₁ is a subsurface of S corresponding to the area covered by the electrodes.

The Finite Element Method (FEM) we developed finds the numerical solution of equation (1) with the boundary condition as specified in equation (2) when the electromagnetic field activities of a head are modeled. Our implementation has the following special features: a) A procedure to generate a list of wedges within a head such, that each wedge has an unique tissue type and the interfaces between two neighboring wedges are perfectly matched to guarantee a correct numerical solution to Eq.(l). b) A proper numbering scheme on a chosen set of nodes and the use of local supporting functions to produce a diagonally condensed block sparse coefficient matrix. c) The use of elementwise local shape functions to represent the local supporting functions. d) The use of a compressed data structure to obviate the need to allocate a very large memory space so that large size problems can be solved on desktop computer workstations. e) Use of a C-language sparse matrix solver incorporating a temporary working array in the process of Cholesky Decomposition to speed up the process of computing the numerical solution.

Advantages of this approach are: 1) it handles the Dirac Delta Function W smoothly by transforming equation (1) into a variational form when J represents a dipole-like kind of source; 2) it allows modeling of the complicated geometires of different tissues within the head by generating finite elements using contours on pairs of adjacent MR images: and 3) it produces a sparse matrix in which many entries are zeroes and thus is convenient and fast to compute on a desktop computer workstation. For example, a brain model with 10,087 FEM nodes can be solved in about 100 minutes on a SUN Sparc-1 workstation (12 MIPS) with 8MB of memory (Figure 5).

II . SPATIAL ENHANCEMENT TECHNIQUES TO REDUCE BLUR DISTORTION Elecincal currents generated by sources in the brain are volume conducted through brain, cerebral spinal fluid, skull and scalp to the recording electrodes. Because of this, potentials due to a localized source are spread over a considerable area of scalp and the po.cnuai measured at a scalp site represents the summation of signals from many sources over much of the brain. In accordance with the current invention, two spatial enhancement methods are described which correct this blur distortion. Both methods use a realistic model of the shape of the head, and neither method requires specification of an arbitrary source model (e.g., current dipoles).

The simpler method computes an accurate csumaie of the surface Laplacian Derivation (LD), which is proportional to local normal current at the scalp. This has the advantage of eliminating the effect of the reference electrode used for recording, and of eliminaung much of the common activity due to either the reference electrode or volume conduction from distant sources. The disadvantages are that the LD does not produce valid values at the outermost ring of electrodes and it does not correct for local differences in skull thickness and conduction propeπies. The most accurate surface LD uses the actual measured electrode positions, and estimates the LD over the actual shape of the head using a 3-D spline algorithm.

According to a method called Deblurring. a further reduction in distonion can be obtained using a finite eiement or boundary integral model of the cortex, cerebral spinal fluid, skull and scalp to estimate the potentials which would actually be recorded on the surface of the brain. Deblurring is a "downward continuation" method in that, without prior knowledge or assumptions about the generating sources, the cortical potential distribution is derived given the scalp potential distribution and a realistic model of the conducting volume between the scalp and cortical surfaces. In the preferred implementation of Deblurring, a transformation matrix is constructed based on the geometry and conductivities of the finite elements which predicts the cortical potentials for any given set of scalp potentials.

For either procedure, a means is needed to align the EEG electrode positions with the outside surface of the scalp. Accordingly, according to the current invention, x, y, z translation and x, y, z axis rotations are computed iteratively to align the measured or computed positions of the EEC electrodes with the external surface of the scalp measured by a magncuc resonance image (MRI). computerized tomographic or ultrasonic scan, or mechanically. This alignment is done for each electrode by finding the distance to the closest point on the scalp surface and minimizing the mean distance for all electrodes. With MRIs that have a 3 mm inter-slice spacing, a mean error distance of better than 2 mm is achieved. This is more accurate and less subjective than alignment procedures that use skull landmarks such as the πasion and pre-auricular points located visually in the MRI. Figures 5a (left) and 6b (left) show the electrodes displayed schcmaucally on a scalp surface reconstructed from hon zontal MRI scalp contours.

THREE-DIMENSIONAL SPUNE LAPLACIAN DERIVATION METHOD

3-D spline Laplacian Derivation

This section describes a new implementation of a common means for achieving spatial enhancement of the EEG signals using a Laplacian operator. Distortion of neuroeiectric signals recorded at the scalp occurs principally because transmission through the low-conductivity skull acts as a spatial low-pass filter. A well-known signal processing technique to enhance higher spatial frequencies is to estimate the rate of change in the signal over a small distance using the Laplacian operator, the discrete approximation of which is referred to as the Laplacian Derivation. The original implementation of the Laplacian Derivation modeled the scalp surface as a plane and calculated a simple second-order finite difference from a rectangular array of equally-spaced electrodes. Spherical and ellipsoidal models have also been used, for which the measured electrode posidons are projected onto this surface, all of the measured potentials are used to calculate a spline-interpolation function, and the Laplacian operator is computed analytically from this function for each electrode position. The advantages of the method described here are construction of a realistic and accurate surface model of the head from the measured electrode positions, use of an accurate three-dimensional local spline interpolation function, and a stable numerical integration method of estimating Laplacian values using Taylor expansions and the least square error minimization technique. Construction of the triangle mesh surface from electrode locations is performed in a step-wise procedure designed for primarily convex surfaces with no folding. In the first step the two closest neighbors of any selected electrode arc found which determines the first triangle in the mesh and the current mesh boundary. The next step finds the closest electrode outside the boundary, adds a new triangle, and revises the mesh boundary. This is repeated until all electrodes have been incorporated into this triangle mesh. Also, after each step, a triangle is added using three consecutive electrodes on the boundary, when necessary, lo make the boundary edge convex. At the end. all adjacent triangles are also checked to see if an obtuse angle might be reduced by substituting the two alternative triangles.

For a given electrode, the information required to estimate the desired local surface laplacian is the local surface geometry associated with this electrode and the local potential distribution. In reality only the polygons defined by the given electrode and the surrounding electrodes, and potential values at these points are used. Local surface construction

For a given electrode, all its surrounding electrodes are found and a sublist of polygons using these electrodes only is formed from the global triangular mesh. If the resulting polygons in the sublist form a complete circle around the given electrode, this point is classified as an internal point and this sublist of polygons is ciassied as the local surface associated with this given electrode. Otherwise, it is classified as an edge point and no Laplacian Derivation operation will take place for this electrode since the information of local potential distribution around this point is not complete.

Local potential distribution

Potential values at the given electrode and its surrounding electrodes, say n, are used to construct the following local spline interpolation function P(x,y,z) to approximate the local potential distribution function.

The coefficients a_j are determined by solving a set of simultaneous linear equations which force the interpolated values at the electrode positions to be equal to the measured potentials, and which satisfies the following:

Function P(x,y,z) defined in Eq.(1) is a spline-like function as it is a linear combination of reciprocals of distances. A nonzero constant w assures that this function is infinitely differentiable. The value of w which has been used is the average inter-electrode distance. The function P (x,y,z) is reference-invariant in the sense that if the measured potentials are re-referenced, the interpolated values change by the same constant.

Laplacian Derivation operation

The Laplacian Derivation operation at electrode (x_o,y_o,z_o) is the following quantity:

where x = f (ξ,η), y = g(ξ,η) and z = h(ξ,η) is the parametric representation of the local surface reconstructed, and x_o =f (ξ_a,η_a ); y_o = g(ξ_a,η_a); zo = h(ξ_a ,η_a ).

The desired quantities and at (ξ_a,η_a) are numerically solved using the three-term

Taylor expansions of twelve surrounding points (ξ_b,η_b) in its parametric space. They are: α

In our programming implementation, (ξ_b,η_b) are twelve equally-spaced points each at a radial distance ε from (ξ_a,η_a) to guarantee numerical stability when equation - 6 is solved by a least-squares method. The function values V(ξ_b,η _b) (b=1,2,...,12) are calculated by the local interpolation function P(f(ξ_b,η_b),g(ξ_b ,η_b),h(ξ_b,η_b) ), and V(ξ_a,η_a) equals P(x_o,y_o,z_o), which in turn is the measured potential value at the given electrode.

This implementation can be made more efficient for Laplacian operations on an EEG time series by viπue of the fact that, when a triangular mesh is made, the local surface, the structure of the local spline function and the Laplacian Derivation operation for any given internal electrode are fixed. Therefore, a Laplacian weight vector associated with an internal electrode is explicitly constructed which yields the desired Laplacian Derivation when multiplied by the measured potentials on this local surface. Figure 7 shows the dramatic improvement in spatial detail obtained with this procedure.

DEBLURRING METHOD

A better improvement in distonion correction is possible by using a finite element or surface intergral method to represent the true geometry of cortex, cerebrospinal fluid, skull and scalp, and estimates of the conductivities of these tissues, to model the potential activity described by Maxwell's equations in the scalp and skull layers. This allows estimation of the potentials which would actually be recorded on the surface of the brain. This can be done without introducing an arbitrary model of the actual number and location of sources because, regardless of where they are generated in the brain, potentials recorded at the scalp must arise from volume conduction from the conical surface through the skull and scalp.

The finite element method (FEM) is used to compute the numerical forward solutions of Maxwell's equation as this method is numerically stable and well-suited for handling the complex geometry and varying conductivities of the head. A multidimensional, bipolar vector space optimization scheme is developed to search for the optimal conical potential distribution whose FEM-based numerical forward solution best fits the actual measured scalp potentials. Since previously, there has been no method of deriving realistic 3-D volumetric finite elements, in accordance with the current invention, we have developed a method for generating finite elements for the entire conducting volume using magnetic resonance images of a subject's head or other measurement means such as mechanical measurement of the shape of the outside of the head and ultrasonic measurment of the thickness of the scalp and skull and outside shape of the brain.

1. Theoretical Formulation

Brain electrical activity may be described by Maxwell's Equation,

∇•( σ(x,y,z)∇u(x,y ,z) ) = F(x,y,z) in Ω, ^{( 7 )} where Ω is the conducting volume, σ(x,y,z) is the conductivity tensor, u(x,yj) is the potential distribution function, and F(x,y,z) is the source function within.

For various electromagnetic problems, Eq. 7 coupled with some boundary conditions specified on the surface of Ω uniquely defines a physical model on Ω. Currently, we choose Ω as the volume bordered by the scalp surface, denoted by S_scalp , and the cortical surface, denoted by S_cortical, extending from the top of the head to a horizontal cut-off surface, denoted by S_cut-off. To warrant a unique solution to the model, the boundary condition on S_scalp is specified as σ∇u(x,y,z)•n = 0 on S_scalp, (8) since electric current does not flow out of the head in the direction normal to the scalp surface

(denoted by n). The boundary condition on S_cortical and S_cut-off is specified as

where G(x,y-z) is the conical potential distribution function. With this kind of choice of Ω, source function F(x,y ,z) will be zero everywhere within Ω since the scalp and skull layers do not contain generating sources. Thus there is no need to determine F(x,y,z).

It is well known that for a given function G(x,y,z), a unique solution u(x,y,z) exists which satisfies Eqs. 7, 8 and 9 . This process is often referred to as a forward problem. However, for the Deblurring process, knowledge of G(x,y,z) (one of the required boundary conditions to warrant a unique forward solution) on S_cortical is desired in addition to information about u(x,y,z) in Ω. This requirement defines the Deblurring process as an inverse problem, and the usual strategy is employed of using additional apriori information to compensate for the fact that one of the required boundary conditions is unknown. This additional information is used to help find the desired boundary information G(x,y,z) as well as the corresponding forward solution u(x,y,z).

In the EEG application, the measured potential at the scalp surface S_scalp provides additional information in the form of an extra boundary condition. That is u(x,y,z) = U(x,y ,z) on S_scalp (10) where U(x,y,z) is the potential distribution measured at the scalp. A function G^•(x,y,z) is taken as the desired soluuon to a Deblurring calculation if the corresponding forward solution u^•(x,y,z) of Eqs. 7,8 and 9 yields the given scalp potential measurement f/(x.y-z) at S_scalp. The function G^•(x,y,z) is obtained through an optimization method which produces a convergent sequence of candidate functions G_i(x,y,z) such that their corresponding forward soluuons u,(x,y,z) of Eqs.7 and 8 approach the extra boundary condition U(x,y,z) at the scalp surface. 2. Simulation Studies of the Validity of the Deblurring Method

Simulation studies demonstrate the validity of the Deblurring method. A simple conducting volume is used in the simulations for which the scalp potenual distributions satisfying Eqs. 8 and 10 and their corresponding cortical potential distributions can be exacdy determined analytically. The validity of the Deblurring method for this particular conducting volume can thus be established by comparing the deblurred data computed from these scalp data with the exact conical distribution.

A multi-layer half-shell model Ω_shall is chosen because it is a first approximation to a realistic head model and the exact soluuon u (x ,y ,z ) of Eqs. 7. 8 , 9 and 10 for the half-shell model can be obtained. For example, when Ω_shall is the upper-half two-layer shell model embedded within a three concentric sphere model, the corresponding three sphere analytical forward solution v(x ,y ,z ) due to some sources within the inner sphere but outside of Ω_shall will be the exact solution u (x ,y,z ) for Ω_shall when v(x,y ,z ) is confined to Ω_shall. Based on this observation, the following simulation procedure was designed and applied. First, the scalp, skull, and brain surfaces were modeled as top hemispheres with radii of 9. 8.5. and 8 cm. and the conductivity ratios of scalp/skull and brain/skull were taken to be 80. The scalp and skull volumes constituting the Deblurring region of interest, Ω_shall , were regularly subdivided into a total of 432 prism-shaped finite elements, each formed from two triangles located on adjacent spherical surfaces. There were 127 outer scalp points and 127 inner skull points. Second, the three-sphere analytical forward solution v(x ,y ,z ) was formulated for a dipole lying within the inner sphere. Third, the potentials for finite element vertices on the outer hemisphere were calculated by the analytical forward solution v (x ,y ,z ) to constitute the numerical boundary condition of Eq. 10 at the scalp surface. Fourth, the finite elements generated were used to discrciizc Maxwell's Equation, and the potentials for finite clement vertices on the inner hemisphere were solved by the Deblurring method. Finally, the exact potentials for these same vertices were calculated from the analytic solution v(x,y ,z ), and these potentials were compared with the deblurτcd ones.

The chosen deblurring region Ω_shall . which is embedded within a three-sphere model, allows use of the three-sphere analytic solution to generate various kinds of exact scalp data and conical data from some source located outside of the Deblurring region so that the computed data from the Deblurring method can be compared with the exact data. Table I presents the error ratio of the computed cortical potentials deblurred from the given scalp potentials and the exact conical potentials produced by six specifically chosen current dipoles. These dipoles vary in depth, orientation, and vicinity in relation to the simulated Deblurring region, and they generate a variety of scalp potential patterns. Positions and oπentauons of the dipαlcs are tabulated in spherical coordinates (p, ɸ, θ). where p is the eccentricity, o is the azimuth angle and θ is me inclination angle. The first column of Table I (dipole type) provides a descriptive summary of each dϊpole's position and orientation. In this table, if the ɸ.θ components of the dipole's position is equal to the dipole's orientation, the dipole is "radial;" otherwise it is "tangential." If the dipole is located close to the Deblurring region, a is a peripheral" dipole; otherwise it is "central." If the dipole's position is no more than 1.0 cm from the simulated conical surface it is "superficial," if greater than 1.0 cm it is "deep." The error ratio is defined as

As can be seen in the table, the error measure ranges from 0.2% lo 1.2%. Inspection of the data suggests that the most obvious differences between the deblurred and actual cortical data are due to the inexact numerical approximation of the cortical surface produced by the large finite elements. Overall, the simulation results demonstrate that the deblurred data match the theoretical data quite well.

3. Deblurring on a Realistic Head Model

The Deblurring procedure becomes more complicated when the region to be deblurred. Ω. is a realistic multi-layer head model. 3.1. Finite Element Modeling of Scalp and Skull Layers

A crucial first step is to represent the region of interest as a set of similar-shaped volumetric finite elements so that Maxwell's Equauon (Eq.7) can be discrciizcd. using these finite elements, into a system of simultaneous equations. We developed a procedure lor generating tcirahcdral finite elements as this shape is well-suited for modeling the complex geometry of the head. In practice, horizontal magnetic resonance image (MRI) slices are first inspected on a workstauon to determine the most ventral (bouom) image. Sixty-four radial lines emanating from the center point are superimposed onto the image to guide interactive placement of contour points so that each contour has consistent inter-point spacing (Figure 8 ) . Contours are generated to delineate the scalp/exterior border, the skull/scalp border, and the brain/skull border. Finite elements are automatically constructed from the contours in the following manner. For each point on a given contour for a given MRI slice, edges are formed between that point and: (a) adjacent points on the same contour for the same slice, (b) the homologous point and the next point lying on an adjacent contour of the same slice, (c) the homologous point and the next point lying on the same contour of an adjacent slice, and (d) the homologous point and the next point lying on an adjacent contour of an adjacent slice. Triangular faces are formed from the edges and tetrahedrons are formed from the faces. These tetrahedrons, which partition the entire volume, constitute the finite elements (Figure 9 ) .

22. FEM Deblurring Calculations

The Deblurring procedure is summanzed in Fig.10.Starting with an initial estimate of the cortical potentials, SG sub O (x, y, z)S. a forward solution of Maxwell's equation is computed using the FEM. This produces an estimate of the scalp potentials which is compared with the actual measured scalp data. The procedure ends if the goodness of fit reaches an acceptable criterion; otherwise, the multidimensional optimization scheme searches for a better estimate of the conical potentials and the loop contained inside the large frame in Fig.10 is repeated until convergence is reached.

Deblurring is comprised of three main computational tasks. First the FEM is applied to a conducting model composed of teirahcdral elements to find a numerical forward solution Su(x, y, z)S of Eq. 7 with boundary conditions specified in Eqs. 8 and 9 through a transformation matrix approach. Second, a sparse solver is used to find the initial numerical conical potential distribution function SG sub O (x, y, z )S. Third, an efficient multidimensional optimization scheme is employed that searches for an optimal cortical potential map SG sup * (x, y, z)S through construction and testing of a sequence of candidate functions SG sub i (x, y, z)S at SS sub conical union S sub cutoffs to find the one whose numerical forward solution Su sup * ( x, y, z)S best fits the given measured scalp potential data.

Before the FEM can be applied, the potentials measured at a relatively small number of scalp sites, say 124, must first be interpolated onto the scalp surface, which consists of some hundreds of vertices of the tetrahedrons. These interpolated potential measurements at the scalp are compared with the solution at each step of the Deblurring procedure. If the potential distribution on SS sub scalpS is denoted as SU(x, y, z)S (the extra boundary condition), the function values of SU( x, y, z )S at 124 points on SS sub scalpS are experimentally defined by the recorded data SU sub iS. A spline-like function

is then constructed and its 125 unknown coefficients, α,-, are determined using the 124 constraints U_i =P (x_i,y_i,z_i) and the constraint that

where (x_i ,y_i,z_i), i=1,2,...,124 are the channel positions, U_i, i=1,2,...,124 are the directly measured scalp potentials, and w is the averaged inter-electrode distance. Having P(x,y,z), the numerical approximation to function U(x,y,z), the potentials for every finite element vertex on the MRI-determined scalp surface S_scalp are estimated, and in so doing, the boundary condition in Eq.₁₀ is numerically established for the Deblurring procedure to proceed.

Forward solutions using the transformation matrix

Applying the finite element method to Eq. 7 with boundary conditions in 8 and 9 and using standard conductivity values for individual finite elements, the following matrix vector relation is obtained

C u = 0. (11 )

Here C is the global conductance matrix and u is a vector of dimension n_wid , which is a numerical approximation to the continuous potential distribution function u(x,y,z) defined in the Deblurring region Ω. The value of n_wid corτesponds to the total number of vcnices on all the finite elements in Ω. Next, u is partitioned into three parts, which correspond to those vertices on S_scalp, on , and the remaining vertices within the Deblurring

region, denoted by u,, u₃ and u₂ respectively. If the matrix C is partitioned correspondingly, Eq. 11 becomes:

When G(x,y-z) is given , i.e., u₃ is known, Eq.(6) is rearranged into

The resulting vectors u, and u₂ constitute the standard forward solution of Eq. 7 coupled with Eqs. 8 and 9 for a given G(x,y,z) (or for a given u₃). Eq.13 is the numerical representation of the underlying biophysical model formulated above.

During the Deblurring process, a sequence of forward solutions is produced with the same set of finite elements but different conical candidate maps (or different u₃'s). Thus, it is most efficient if the desired forward solution can be derived by simply applying a fixed matnx to a given conical potential map. Actually, only the partial forward solution at the scalp surface is needed, which is compared to the measured potential at the scalp (the extra boundary condition). This calls for a special-purpose method that reduces the time required for solving a forward solution at each step to that for a simple matrix-vector multiplication at each step. It can be shown that if the coefficient matrix of Eq. is invencd and the inverted mauix is prcmultiplicd to both sides of Eq. 13 the top half of the resulting matrix on the right side can be extracted to yield the desired partial forward solution U| for a given u₃. This transformation matrix is denoted by A. u, = Au₃. Because mauix A is derived from the matrix manipulation of the global conductivity matrix C defined in Eq.12 which is invariant for a fixed finite element model representing Ω and for a given set of conductivity values identifying various tissues within Ω, matrix A can be used to repeatedly derive forward solutions in the Deblurring process for the same subject. For a particular model of dimensionality 1886. our sparse solver uses 0.34 M flops and 0.128 Mbytes to construct matrix A. whereas a standard Cholesky Method requires 1118 Mfiops and 28.5 Mbytes to obtain the same solution.

Multidimensional optimization The Dcblurring method is essentially a multidimensional optimization scheme that searches for an optimal conical potential distribution whose forward solution best fits. the measured scalp potential distribution. In practice, the search starts with G_o(x ,y ,z). computes the corresponding forward solution u₀(x,y ,z ), and measures the goodness of fit between u_o(x,y ,z ) and the measured U (x ,y z) at the scalp surface. G ₁(x ,y ,z ) is then sought by the optimization scheme using the goodness of fit measurement and the process is repeated. The iterative process is stopped at step i₀ when G,₀(x ,y ,z ) produces a forward solution that is acceptably close to U(x,y .z).

The first step in the optimization scheme is choosing a good starting point for G₀ that is fairly close to G^a(x,y,z). The closer G₀ is to G^a (x ,y ,z ). the fewer the iterations that will be needed to find a forward solution u(x,yj) acceptably close to U(x ,y ,z ). We rearrange Eq. (6) into

and use the measured potential U (x ,y -z j in Eq.10 as

to solve for

The solution of Eq.14 is the forward solution of Eq. 7 coupled with boundary condition 10 on S_κιύp and boundary condition 8 on .

This is an auxiliary system to the forward solution system and the corresponding forward solution u^ provides a reasonable starting point for the proposed optimization scheme.

For a given transformation matrix A and a starting vector

(the numerical representation of the initial cortical potential distribution function G_o(x ,y ,z )). the goal of the iterative process is to find a vector

which is the numerical representation of funcuon G^a (x ,y ,z ). such that the resulting goodness of fit is maximized. The dimension, n , of the vector u₃ is the total number of MRI-dcrivcd finite clement vertices on the conical and cut-off surfaces. If wc define

) as the numerical representation of G,(x ,y ,z ) at step i, then each G, can be visualized geometrically as the location of a point in R*. The following optimization algorithm describes an iterative procedure for finding the numerical funcuon G, given G.., and its corresponding goodness of fit.

1) An n -dimensional hypcrcubc of radius is created with G_i-l at its center. (An initial radius value, r_D is supplied as a parameter to the algorithm. This starting radius is chosen by computing the mean of the scalp potentials.)

2) A point, chosen from the possible 2^a candidates populated on the surface of this hypercube, is found in the inner ucrauve stage such that

maximizes the goodness of fit without testing all of the 2ⁿ candidates.

is defined as

Because s₁ takes on values of either 1 or

is referred to as a bipolar vertex. The inner iterative stage consists of the following:

2.1) Given one initial bipolar vertex on the surface of this hypercube (initially chosen in the same direction as ΔG_i-2), the goodness of fit for that vertex and for the n adjacent bipolar vertices is computed.

2.2) If all the goodncss-of-fit measures for the n neighbors are worse than the one for itself, the current bipolar vertex is considered the best and the outer stage of iteration is resumed at step 3).

2.3) Otherwise, make the neighboring bipolar venex with the best goodness of fit the new initial bipolar vertex anxd go back to step 2.1).

3) If the resulting goodness of fit is better than that for G_i-1, the vector

is chosen as G_i, r_i+1 = r_i.

4) Otherwise the current hypercube is shrunk by half in radius and step 2) is repeated on a reduced hypercube. If the hypercube has not shrunk for th ree consecutive iterations, the radius of the new hypercube is increased by 75% to speed up the movement of points G, in vector space R*.

In summary, this search method proceeds by searching for a local optimal point on an intermediate domain, i.e., the surface of an n -dimensional sphere. When a better point is found on the surface, a new n -D sphere with the same radius is set up around this better point and the next better point is sought on the new surface. Otherwise the global optimal point must be contained within the current sphere, in which case the original sphere is shrunk by half and the next better point is sought on the new surface. Efficiency of this technique is achieved by exploiting the fixed structure of the transformation matrix A. and by the selection of the 2" bipolar vectors embedded in the n -dimensional sphere. The advantage of using the 2" bipolar vector structure is not only that these T vectors happen to be a set of vertices of a hypcrcubc and therefore, the maximal and minimal points within this hypercube can be estimated by the maximal and minimal points among its vertices, but also that it allows the inner search loop to be executed efficiently.

4. Application to Data from Human Subjects

4.1. Comparison of Estimated and Measured Cortical Data in an Epileptic Patient

To veπfy the results of the Deblurring procedure for a realistic head model, botii scalp and cortical averaged evoked potential (EP) data were obtained from an epileptic patient. Deblurring was performed on the scalp data, and the results were compared with the measured conical data. The MRI-dcrived finite element model for tiiis subject had 5139 ictrahcdral volumetric elements and 1300 vertices, among which 438 were scalp vertices and 393 were vertices at the conical and cut-off surfaces. The transformation matrix A. therefore, had a dimension of 438 by 393 and the vector space in which the inverse solution was sought had a dimension of 393. Figs, 11a-b show the scalp and conical brain surfaces for this subject and color isochrome maps of three different sets of EP data. The 15-Hz EP recorded from the scalp (referenced to digitally linked mastoids) is shown on the left side of Fig. 11a while the estimated cortical EP calculated via the Deblurring procedure is shown on the right side of Figure. 11a The Dcblurring results are compared with the actual cortical recording in Figure lib; the grid isochrome map is superimposed onto the (monochrome) conical surface on the right side of the figure, and the Deblurring results appear on the left. The original scalp EP distribution shows a single, near circular area of maximal potential in the mid-panetal area of the left hemisphere. The dcbluncd computed conical EP shows a more compact and irregularly shaped maximum in about the same area, with a polarity reversal towards the midline. The recorded actual cortical EP shows a single, even more compact maxima in about the same area. The computed cortical EP better matches the actual conical EP than docs the original raw EP recorded at the scalp.

Conductivity Variations

This result was based on a standard scalp/skull conductivity ratio of 80. Deblurring calculations were also performed using other values of the conductivity ratio. Figure 11c shows the Deblurring results for ratios of 40, 80 and 120. For these and other values examined. difTerent conductivity ratios produced differences in the size of active regions and in the finer detail but did not dramatically change the overall pattern of activity. 4.2. Deblurring Results for Healthy Subjects

The results of applying the Deblurring procedure to EP data obtained from a healthy subject are illustrated in Fig. 12. Fig. 12a shows the results for 15-Hz electrical stimulation of the right index, the right middle, and the left middle fingers. As would be expected, potential maxima are on the side of the coπex tiiat is coπtralateral to the stimulated hand. Stimulation of adjacent fingers of the same hand produced distinguishable patterns of cortical activation. Figure 12b compares the original EP at the scalp (recorded referenced to digitally linked mastoids) with the deblurred EP for stimulauon of this subject's left middle and right index fingers. The increased spatial detail produced by the Deblumng method is apparent.

5. Surface Integral Method

An alternative to performing dcblumng (or dipole source localization) through using numerical solutions to the differential form of Maxwell's equations is to use numerical solutions to the integral form of Maxwell's equations.

The following equation is used for the source localization problem:

where σ denotes conductivity;

is a point on the outer surface of the scalp, u is potential:

is the vector from a source to and R_j is its magnitude:

is the moment of a dipole source: / is the number of interfaces, S, , between regions differing in conductivity;

is the unit vector from the in region to the out region: and

is the vector from to the point of integration with denoting the unit vector and R denoting the magnitude. A numerical solution of this equation involves obtaining a discrete approximation to the integral. Boundary clement methods can be used to do this. An alternative is to use efficient integration methods that have been developed for the sphere and methods of computational geometry for projecting a non-spherical surface onto a sphere.

One advantage of using the integral method for source localization is that no nodes are close to the sources, a locus where potential fields vary a great deal. This is a problem for the finite element method. The following equation is used for the deblurring problem:

where the notation is that for Equauon 15 -After rcguiarization, the numerical integration techniques used for source localization are applied. An advantage of diis formulation when applied to deblurring is that it accounts for the geometry of the whole conductor whereas the finite eiement technique only accounts for the portion of the conductor that includes the electrodes on the scalp and the tissue immediately below.

6. Estimating Scalp and Skull Thickness for Models

In the finite element or surface integral methods, the local geometry of scalp and skull are each explicitly represented. If a magnetic resonance image or computerized tomography (CT) brain scan is not available, an alternative metiibd of measuring tissue thickness can be used based on low-frequency ultrasound technology.

For a low-resolution representation of skull thickness, the "end echo" in standard echoencephalography is used to estimate the thickness of the skull. This requires placing the ultrasonic probe at a position on the opposite side of the head from the place of desired measurement, normal to the scalp at the probe, and normal to the skull at the place of skull thickness measurement on the opposite side of the head from the probe. Skull thickness at each of the major cranial bones is measured in the following locations: cranial bone position probe position

left temporal T3 near T4, directly above the ear right temporal T4 nearT3. directly above the ear frontal (left) aF1-aF3 roughly 1 cm below 02

frontal (right) aF2-aF4 roughly 1 cm below 01

left occipital 01 .2-.5 cm right of FPZ

right occipital 02 .2-.5 cm left of FPZ

left parietal aP3 1 an below M2, with upward tilt right parietal aP4 1 cm below Ml, widi upward tilt

The end-echo consists primarily of the echo due to reflection from soft-tissue to inner skull interface and the echo due to outer skull to scalp interface (followed by a very small scalp-to-air echo). The delay prior to the end-echo varies, according to scalp location and inter-subject variation, roughly from 90 to 120 microseconds. The delay between the echos from either surface of the skull may vary from 1.7 to 5 microseconds. The second echo is auenuatcd with respect to the first, and requires a gain increase of 2.5 dB per MHz (typical) to be comparable in size to the first. Pulse widdis of .25 to 1 microseconds are used with burst frequencies of 2 to 10 MHz, adjusted to produce the best discrimination between the two echos.

As necessary, passive receptors arc used to scan the area of the intended skull thickness measurement to determine the accuracy of aim of the transmitted beam. When necessary, an average is taken over a small cluster of locations (3 to 6) near the listed probe site, in order to average over the variations in thickness for a given bone near the target area. In those instances when the aimed beam cannot be focused well on the parietal bones, the average thickness of these bones is estimated from normalized data on the relationship between parietal bone thickness and the average thickness of the other bones.

For finer spatial resolution in measurement of skull thickness, a near side detector can be used. It has two to four transmitters adjacent to the transmitter/receiver probe, spaced far enough away so tiiat the critical angle (27 degrees) prevents energy from the adjacent transmitters from entering the skull at an angle that could be reflected back to the probe. The phase and time delay between the adjacent transmitters and the probe is adjusted to minimize the reverberation at the probe. Thus, at the probe site the adjacent transmiuers are contributing energy soley throug rcfiections from the outer skull surface, and in a way to minimize reverberation at the probe. Only energy transmiucd by the probe itself can enter the skull and be reflected back to the probe.

The thickness of the near skull can then be measured either by subtracting independent measures of the reflection time from the outer skull surface (adjacent transmitters inactive) and reflection time of the echo from the inner skull surface (reverberation minimization active), or by inierferometry. The interferometry option sweeps the sound frequency (readjusting reverberation minimization for each frequency), and determines the change in fre quency required for one oscillation in the variation in reflected energy due to interference of the inner and outer reflections from the skull. The measured thickness of the skull is one half the speed of sound in the skull divided by the change in frequency required to move through one "interference fringe".

7. Estimating Conductivities

Up to now. estimates of conducuvity for the scalp, skull, CSF and brain have been from standard published sources. More accurate values are needed to realize the full benefit of the FEM head model. This can be done by solving for the ratio of skull-to-scalp and skull-to-brain conductivities using the scalp potential distribution produced by a single compact population of neurons, represented as an equivalent current dipole, whose location is approximately known. In one procedure, a steepest-descent, nonlinear, least-squares fit between the measured potentials and the forward solution of an equivalent current dipole in a three-sphere model is performed. The main unknown parameters are the ratios of conductivities of the tissues, as well as the dipole suength (which is not of interest here). The availability of reasonable initial guesses, as well as upper and lower bounds, from existing experimental data result in realistic solutions for the unknown parameters. A highly localized source suitable for this purpose is the somatosensory evoked potential peak at 20 msec produced by mild stimulation of the contralateral median nerve at the wrist. This potential, negative posterior to the central suicus and positive anterior to it, is generated by a single compact population of neurons onentcd taπgentially in the hand area of somatosensory area 3b and located about 2.5-3.0 cm below the surface of the scalp. Published values of the amplitude of the cortical potential and the ratio of conical to scalp potentials arc also available. The global conductivity parameter is then adjusted for each of the finite elements in the skull using local skull thickness measures determined from the low frequency ultrasound. MRI or otiicr means.

An alternative method for estimating tissue conductivity is based on measuring the reverse electromotive force produced by eddy currents induced in the tissue by an applied external field generated by a coil placed against the head of each patient. The magnitude of diis emf depends on the conductivity of the tissue and the frequency of the magnetic field, ie. higher conductivity means more eddy current is induced which gives a larger reverse emf. Very low frequencies induce eddy currents in all tissues whereas high frequencies induce currents only in the tissues closest to the coil. Scanning a spectrum of frequencies yieids data from which the conductivities of each tissue type can be computed.

Claims

CLAIMS :

1. The method of improving the spatial resolution of the data recorded from a plurality of electroencephalograph (EEG)

electrodes removably and electrically connected to the scalp of a subject, including the steps of:

(a) amplifying the brain waves detected at the electrodes, converting the amplified brain waves into digital data and entering the data into computer system means to analyze the data;

(b) measuring the thickness of the skull of the subject, the thickness of the scalp of the subject and the head shape of the subject;

(c) measuring the physical locations of the electrodes on the scalp of the subject;

(d) estimating the electrical conductivity of the skull, scalp and cerebral spinal fluid of the subject; and

(e) entering the measurements derived from (b), the

electrode locations derived from (c) and the conductivity

estimates derived from (d) into the computer system means to analyze the entered brain wave data.

2. The method of claim 1 and including the step of obtaining the measurements of (b) by magnetic resonance image (MRI).

3. The method of claim 1 and including the step of obtaining the measurements of (b) by computerized tomographic scan (CT) and mechanical measurement.

4. The method of claim 1 and including the step of obtaining the measurements of (b) by an ultrasonic scan and mechanical measurement.

5. The method of claim 1 wherein in step (e) the EEG electrode positions from (c) are registered with measurements of the external scalp surface from (b) using an iterative numerical procedure which optimally matches each EEG electrode position with the nearest scalp surface point while minimizing the distance for all other electrodes.

6. The method of claim 1 in which in step (e) the computer system means uses a numerical finite element or surface integral model of the subject's scalp, skull and brain surface.

7. The method of claim 1 in which the estimate (d) for the skull is estimated using an equivalent current dipole source model of a sensory stimulus whose major area of cortical

activation is known.

8. The method of claim 1 in which the estimates (d) are obtained by measuring the reverse electromotive force produced by eddy currents induced by an applied external field generated by a coil placed against the head of the subject.

9. The method of claim 1 and including the step of generating volumetric elements from contours drawn on medical images or scenes.

10. The method of claim 1 wherein in step (e) Maxwell's equation and a numerical method finite element method or surface integral method are used to obtain the spatial deconvolution.

11. The method of claim 1 wherein in step (e) the measured EEG potential values are interpolated into the potential values for the scalp finite element points using a local spline

interpolation function which is a linear combination of a set of basis functions which are solutions of Maxwell's equations.

12. The method of claim 1 wherein in step (e) a forward solution from the brain surface potential values to the scalp potential values is performed using a finite element based forward

transformation matrix which is solved with a sparse solver method.

13. The method of claim 1 wherein the potentials at the brain surface are estimated using a multi-dimensional optimization method to minimize the difference between the calculated and the actual scalp potentials and a bipolar-vector space implementation is employed.

14. The method of improving the spatial resolution of the data recorded from a plurality of electroencephalograph (EEG)

electrodes removably and electrically connected to the scalp of a subject, including the steps of:

(a) amplifying the brain waves detected at the electrodes, converting the amplified brain waves into digital data and entering the data into a computer system means to analyze the data;

(b) directing an ultrasonic source at the head of the subject and measuring the thickness of the skull of the subject by the reflected ultrasonic waves from the skull;

(c) measuring the physical locations of the electrodes on the subject; and

(d) entering the skull thickness derived from (b) and the electrode locations derived from (c) into the computer system means to analyze the entered brain wave data.

15. The method of improving the spatial resolution of the data recorded from a plurality of electroencephalograph (EEG)

electrodes removably and electrically connected to the scalp of a subject, including the steps of:

(a) amplifying the brain waves detected at the electrodes, converting the amplified brain waves into digital data and entering the data into a computer system means to analyze the data;

(b) measuring the physical locations of the electrodes on the subject;

(c) estimating the Laplacian Derivation (LD) proportional to the local normal current at the scalp; and

(d) entering the electrode locations derived from (b) and the estimated Laplacian Derivation derived from (c) into the computer system means to analyze the entered brain wave data.

16. The method of claim 15 wherein the estimate of step (c) is performed using a three-dimensional spline algorithm.

17. The method of claim 15 wherein the estimate of step (c) is computed using the electrical potential distribution derived from interpolation from the potentials at the electrodes.

18. The method of claim 15 wherein the estimate of step (c) is computed using a computer simulated model of the scalp of the subject.

19. The method of claim 15 wherein the estimate of step (c) is obtained using a numerical integration method of estimating

Laplacian values based upon Taylor expansions of the potential distributions and least squares optimization.

20. The method of claim 15 wherein the estimate of step (c) is computed using:

(i) the electrical potential distribution derived from interpolation from the potentials at the electrodes; and

(ii) a numerical integration method of estimating Laplacian values based upon Taylor expansions of the potential

distributions and the least squares optimization.

21. The method of measuring the positions of a plurality of electroencephalographic (EEG) electrodes electrically connected to the head of a subject by:

(a) measuring the head of the subject by measuring a plurality of standard cross-directional lines on the head to obtain a set of head measurements;

(b) entering the set of head measurements into a computer system which automatically determines the classification of the head shape according to a plurality of predetermined head shape classes;

(c) comparing the subject head size to a predetermined head model of the same class in which the position to the electrodes to be measured have been determined; and (d) determining the position of the electrodes connected to the head of the subject by scaling the measurements of the electrode positions on the head of the subject in accordance with the predetermined electrode positions of said head model.

22. The method of claim 21 wherein the measurements of step (a) are obtained by using a stretchable hat having piezoelectric ribbons embedded therein.

23. The method of claim 22 wherein the measurements are stretch measurements obtained from stretching of the ribbons and

including the steps of obtaining a vector of said stretch measurements and scaling said stretch measurements.

24. The method of claim 23 wherein said computer system is a pattern recognition neural network computer.