Registration of MEG/EEG Data with 3D MRI :

Methodology and Precision Issues

D. Schwartz,

D. Lemoine,

E. Poiseau,

C. Barillot,

Corresponding Author:

Abstract

Mapping neuro-physiological functions to high resolution MRI is an effective means to evaluate localization reconstructions and to exhibit the spatio-temporal aspects of dynamic functional processes. The registration step needed between MEG/EEG and MRI is a source of error which, for the worse cases may be greater than errors related to the localization algorithms. Several registration methods can be used: those based on fiducial markers and those based on surface matching. The aim of this paper is to propose a fully automatic surface matching method and to discuss its extended theoretical and experimental evaluation. The registration procedure matches the skin surface, segmented from MRI, and a digitized description of the head performed with a 3D tracker during the MEG/EEG examination. The registration uncertainties at the edges of the MRI volume were estimated to be between 2 and 3 mm. In comparison with commonly used manual methods the improvement in accuracy is significant. Registration uncertainties are smaller than the localization uncertainties usually observed. By minimizing manual intervention, the reliability of the registration process is increased and the accuracy is stabilized. With this automatic registration method the fusion of MEG/EEG localizations with MRI anatomical data gives highly significant information. Finally the accuracy obtained allows the use of complex anatomical constraints in the localization process without introducing large modelling errors.

Key words : MEG, MRI, data fusion, localization, anatomical constraint, registration evaluation, multimodal registration

Introduction

Magneto-Encephalo-Graphy (MEG) or Electro-Encephalo-Graphy (EEG) are two means to study the human brain electro-physiological activity. There is considerable interest in using these two techniques and in studying their limits. Both of them are completely non-invasive and look at the short latency neuronal activities with high temporal resolution (comparing to Positron Emission Tomography - PET or functional Magnetic Resonance Imaging - fMRI). These methods have been applied to a large and still expanding range of applications which include :

• Epilepsy: pre-surgical exploration of patients with a drug-resistant epilepsy [Paetau et al 1993][Ebersole et al 1995][Smith et al 1995].
• Pre-surgical functional mapping : Delineation of functional areas surrounding a lesion [Roberts et al 1995].
• Study of the normal brain by means of evoked potentials [Wood et al 1985][Hari et al 1993][Liégeois et al 1995][Pantev et al 1995].

In this framework, morphological data are necessary, and they mostly come from Magnetic Resonance Imaging (MRI). Mapping neuro-physiological functions to high resolution MRI is an effective means to evaluate the localization reconstructions and to exhibit the spatio-temporal aspects of dynamic functional processes. This fusion helps for instance to validate a localization at a given latency, according to the actual physio-pathological knowledge.

Two main sources of errors can lead to unreliable results in the MEG/EEG data process.

First of all, the localization algorithm may not converge due to bad modeling assumptions (source model or/and volume conductor model). Many researchers have studied this subject, on the theoretical side [Witwer et al 1972][Weinberg et al 1986][Cuffin 1990, 1993] as well as on the clinical side [Cohen et al 1990][Ioannides et al 1993][Gallen et al 1993]. The deduced uncertainties obviously depend on the methods used and are at best within a 5 to 10 mm range for non-simulated signals.

The second source of errors is due to the registration between MEG/EEG and MRI. As above, the uncertainty depends on the registration method and may reach, for the worst cases, values greater than those related to the localization algorithms [Schwartz et al 1995]. These errors can bias the localization results in several ways: the simplest one is an MRI positioning error of the localization when using a moving-dipole fit algorithm [Scherg 1992] (the registration uncertainty is simply added to the localization one); for more complicated localization algorithms using anatomical constraints [Scherg 1992][Dale and Sereno 1993][Bouliou et al 1995], a bad registration may bias the localization algorithm due to errors in the constraint calculation (e.g. the normal of the cortical surface).

This shows that it is as important to take into account registration uncertainties as uncertainties due to the localization algorithms. A bad registration may correlate registration errors and localizations errors making even harder the analysis of results and uncertainties of a localization process.

Several registration methods can be used in the framework of MEG/EEG-MRI data fusion. As, for the PET-MRI data fusion framework [Evans et al 1989], one can identify two families of methods: those based on fiducial markers [Hill et al 1993] and those based on surface matching [Pelizzari et al 1994, Lemoine et al. 1994, West et al 1996].

The aim of this paper is to propose a fully automatic surface matching method and to discuss its extended theoretical and experimental evaluation. The registration procedure matches the skin surface, segmented from MRI, and a digitized description of the head recorded during the MEG/EEG examination (which we will call " headshape "). This algorithm is currently used routinely on the MEG/EEG platform at the University Hospital of Rennes. We will study and discuss errors only due to the registration stage within the localization process (bold arrows on Fig. 1). Our goal is to end up with a registration algorithm that is robust and that includes procedures allowing (as much as possible) an objective evaluation of the registration value and its related residual uncertainties. If the localization capabilities are known well enough, it will be then possible, in a second stage, to delineate the global uncertainty after a MEG/EEG functional signal analysis process.

The overall objective of this paper is to propose a fully automatic surface matching method to register MEG/EEG data with MRI and to discuss its extended theoretical and experimental results. We will first describe how the data acquisition is performed, then the registration method will be detailed and then a theoretical and clinical evaluation of this approach will be presented.

 

1. Acquisition Protocols

1.1. MEG/EEG Acquisition

The main issue of concern here is the sensor position relative to the patient head. The signal acquisition by itself is not relevant to the registration issue. Although Fig. 1 summarizes the different stages of the patient head localization when dealing with the 37-sensor BTi Magnes system, this scheme is quite generic for most MEG/EEG acquisition systems. The only variable parameters are the definition of an a priori fixed referential (the basis referential) and the definition of a second referential used for the localization (the patient referential). In our case, the patient referential is defined by using the three following landmarks : the Nasion, the Left PreAuricular point (LPA) and the Right PreAuricular point (RPA). The basis referential is assigned to the sensor. LPA, RPA and Nasion coordinates are recorded on the patient's head at the beginning of the acquisition procedure by drawing a 1 mm circular trace on the skin with a pen. Then, the patient referential is defined as follows:

O, the origin of the patient referential :

X axis:

Z axis:

Y axis:

Before each acquisition, the patient referential is recorded by using a 3D digitizer device from Polhemus [Polhemus 1991]. The headshape is digitized using the same tool. The headshape contains about 2000 to 4000 points which will be used in our registration procedure. Since the geometrical relations between the Polhemus device and the basis referential (the sensor) are known, the localization process can be done in relation to the patient referential.

 

1.2. MRI Acquisition

The head skin surface is segmented from MRI and must be as accurate as possible to ensure a reliable registration result. Hence, the acquisition is performed by using a high resolution sequence (1-2 mm slice thickness) scanning from the top of the head to the mouth level and must be fast enough to reduce patient's head motions. For coherence, the geometrical distortions of the MRI system must be known and corrected as much as possible. This stage is performed by using the manufacturer's MRI phantoms, one could also correct these distortions by using post-processing procedures [Chang et al 1992]. A typical acquisition protocol on a GE Signa 1.5 T is : Sagittal reformat, TR = 33 ms, TE = 4, Thickness = 1.3 mm, Ex = 1, FOV = 25, Matrix =256*256, mean time acquisition = 18 min. It is then possible to use previously acquired MRI data when they meet these criteria. Since all of the processing is performed in 3D, the MRI acquisition is reformated by using linear interpolation in order to obtain isotropic voxels (voxel size = 0.94 mm).

2. Registration Method

The registration procedure comprises three stages. The procedure starts with the extraction of the shapes to be matched. This means that points of both acquisitions are classified as " registrable " (i.e. regions described on both acquisitions) or as " non-registrable " (i.e. regions described only on one acquisition). The second stage consists in defining the registration cost function which should lead to fast computation as well as reliable registration. Finally the third stage concerns the minimization algorithm which ends up with the registration parameters according to the (expected) global minimum of the cost function. In our case, we minimize for six parameters (3 rotations, 3 translations).

2.1.1. Headshape Description

As described in §1.1, during the MEG/EEG examination, the headshape extraction is performed according to the patient referential defined by the three LPA, RPA and Nasion points. The Polhemus 3D digitizer is used, and we have evaluated its intrinsic accuracy by studying the reproducibility of the measurements on fixed points with several orientations of the digitizer pen (see §3.1.2.2, in our context, others studies have shown that the Polhemus 3D digitizer accuracy is equivalent to the accuracy of other systems (eg. optical detection) [Rohling et al. 1995], [Galloway et al 1994] and [Young et al 1994]). The goal is to define a surface as similar as possible to the one segmented from MRI. The main dissimilarities come from the skin thickness and rigidity relative to the underneath bone. Concerning the skull, the scalp is rather thin and directly in contact with the bone, the digitizer pen does not significantly compress the skin or cause the scalp to slip on the skull. The most critical point concerns the face of the patient. The algorithm needs good data sets from the face in order to make the matching work well and to avoid mis-registration due to the rather spherical shape of the head (which leads locally to a rather flat cost function). For this purpose, the nose and eye orbits are good landmarks and can be retrieved with good confidence from the 3D headshape digitization as well as from MRI. On the other hand, the skin thickness and its elasticity with regards to the skull increases the " digitization noise "in these regions. This noise is quantified by studying the deformation depth of the epidermis under the digitizer pen pressure: this compression varies from 1 to 3 mm in regions where no important muscle masses lie between skin and skull. In detail, cheeks, cheekbones, mouth or jaw digitization should be avoided. Obviously, the patient’s head should be fixed in order to not overstrain the operator or the patient.

The last digitization constraint deals with the number of digitized points and their layout. These two issues are discussed in the evaluation of the algorithm. The accuracy of the registration is proportional to the number of points in the headshape. The points layout is also of relevance since we are using a multiresolution registration procedure which assumes that the set of points are uniformly spread out on the surface head skin.

 

2.1.2. Skin surface extraction from MRI

2.1.2.1. Segmentation

The objective of this stage is to end-up with a connected set of voxels which most accurately describes the skin surface. For this purpose, we propose a method that employs a background noise reduction scheme, able to handle holes in the MRI volume and finally able to remove non-significant regions (i.e. artefacts) which may be present in a MRI volume. This method proceeds in four steps :

Thresholding

This stage removes most of the background noise; usually the threshold is set to values between 30 and 40 for a 8-bit coded MRI sequence. This threshold relies only on the quality of the MRI acquisition and of the 16 to 8 bits transformation.

Negative Transformation and 3D Connected Components Analysis (1)

A negative image transformation is performed on pixel gray levels followed by a 3D connected components analysis using an 18 neighborhood operator in order to get rid of head interior holes (i.e. holes not connected to the skin surface).

Dilation and Erosion

Weak connections such as artefacts are removed by using a morphological opening, the number of iterations can be adapted according to the strength of these connections (typically 2). The topology of the skin surface is preserved by applying the idempotence property.

Negative Transformation and 3D Connected Components Analysis (2)

The last image processing step consists in removing the residual regions not connected to the head after the previous transformations. As before, a negative image transformation is performed on pixel gray levels followed by a 3D connected components analysis.

A wrong segmentation will influence the registration since the skin surface extraction is directly derived from this stage. For instance, the set up of the threshold level must be a trade-off between removing most background noise and keeping the significant data. Lowering the threshold, will connect background noise to the head leaving artefacts on the skin surface while increasing the threshold will remove voxels actually belonging to the surface. After idempotence, the opening procedure does not alter the surface, but if the number of iteration is not well adapted (the idempotence is not reached), surface description errors may be retained. These different problems influence the registration on several ways : i) segmentation errors may leave local deformations of the surface and ii) a one or two voxel scaling artefact may be introduced, (i.e. under- or over-estimation of the surface location), both leading to possible mismatching (See §3.1.2.2).

 

2.1.2.2. Surface detection

Once the segmentation procedure is performed, surface detection is straightforward. A simple test on each voxel neighborhood is performed : all voxels having at least one neighbor outside the head are labeled as belonging to the skin surface.

2.2. Cost Computation

Since the algorithm computes the registration matrix between the two surfaces, it is necessary to evaluate at each iteration step the value of the residual error. A sound solution should minimize the distance of each headshape point to the scalp surface. The computation of the cost function (Euclidean distance) must be done as fast as possible in order to perform the registration in a reasonable time.

 

2.2.1. Distance transform computation

We chose to use the so-called Distance Transform [Borgefors 1986]. The input of this transform is a binary volume derived from the MRI interpolated volume. The output is a volume that contains, for each voxel, the approximation of the distance from the skin surface.

There are two advantages in using the distance transform: The distance transform is computed only once and the estimation of the cost function is as simple as getting the value in the distance transform for each headshape point. We chose a 5*5*5 template (see Figure 2) that ensures a short computation time and an good approximation of the Euclidean distance with an error less than 2 %.

 

2.2.2. Cost Evaluation

A simple estimate of the position of the headshape inside the distance transform volume is given by the mean value of all headshape point distances (called mean distance). This mean value is a good approximation of the global distance between the headshape and the MRI skin surface. This method is very cost effective. Since no significant improvement has been found by using the standard deviation associated with the mean distance, we choose to use only the square mean distance as the cost function. One advantage is that it provides a greater sensitivity for small head displacements and it increases the correlation between the mean distance and the actual mean error ranging from 0.78 to 0.98. These numbers show that our cost function is a good measure of the quality of the registration. Curves 3.a & 3.b show the values of the cost function when a shift (rotation or translation) is applied to a theoretical headshape (see 3.1) and curves 3.c & 3.d show the variation of the cost function for real data.

 

2.3. Minimization Method

Since we deal with an implicit cost function, we chose a minimization method which avoids computation of the derivative of the cost function. Our method is based on the Powell algorithm [Powell 1964].

 

2.3.1. Powell Method

Let f be the cost function to minimize. In order to meet the Powell algorithm requirements, f must satisfy the following constraints around the minimum:

The constraint i) is difficult to verify since we only get an implicit function, so we just assume that the distance transform is regular enough to verify this property. The second constraint ii) is easier to verify: the distance transform function in the vicinity of its global minimum p0 may be approximated by its Taylor series:

If i) and ii) are verified, the principle of the Powell algorithm is to find minima along conjugate direction axes (for a quadratic function A, two directions p and q are conjugate if p.A.q = 0). The method is based on the two following theorems:

Conjugate directions are found by utilizing the following property of the gradient of f, at its minimum, along a given direction. The gradient may be expressed as follows:

If the minimization is performed along a given direction called , the gradient is orthogonal to at the minimum, i.e.:

and is conjugate to .

Therefore, successive estimated directions must be conjugated. If function f is fully quadratic, the minimum is obtained after a single search along each direction. Otherwise, the search has to be iterative and the minimum is found after n iterations. Then, it is considered as a "quadratic" convergence toward the minimum [Press et al 1988]. Since an MRI acquisition can be more or less contaminated by artefacts f may be locally non-quadratic. In that case f may be convex and two search parameters may be highly correlated. According to theorem 2, a new direction xm is a linear combination of the other directions {x1,...,xm-1} therefore when 2 search vectors are linearly dependent, one of the search directions is lost. In order to remove the singularity we replace the linearly dependent direction by a new direction chosen to be orthogonal to the others and assumed to be a conjugate direction.

Our implementation of Powell algorithm may be summarized as follows:

- let n be the number of parameters to be determined (here n = 6),
- let (x1,...,xn) be the n directions in the parameters space,
- let p0 be an initial approximation value expressed in the parameters space,
for i=1 to n
/* Compute a search direction */
for r=1 to n
find lr such that f(pr-1+lrxr) is minimum
then pr=pr-1+lrxr
end for r
for r=1,…,n-1 : xr=xr+1

xn=pn-p0
if there is one dependent direction from the couple then xn = Cross Product (x1,...,xn-1)
/* Minimization along the found direction */
find l such that f(pn+lxn) is minimum
/* New initial point */
p0=pn+lxn
end for i

The minimum p0 is obtained after n iterations.

2.3.2. Multiresolution Analysis

As for other minimization algorithms, we have to deal with local minima that may lead to a sub-optimal registration result. In order to reduce the influence of these local minima, we perform a multiresolution analysis (See Figure 4). At each level of the multiresolution pyramid, we sub-sample the headshape. Then, we perform a minimization and use this result as the initial value for the following step. The search interval is adapted to each resolution step as follows:

let [x-, x+] be the search interval at the resolution level n,
let Sn be the search step of the minimum,
let xmin be the estimated minimum,
then the search interval at the resolution level n-1 is [xmin - Sn, xmin + Sn]; and the step is Sn-1 = Sn/2.

Using this technique we avoid a convergence to a solution rejected in the previous minimization step. The goal of the approach is to smooth the initial volume in order to eliminate non-significant minima in local convexities of the cost function. It is far easier to refine the registration result at higher resolution steps when the initial position is close to the global minimum. It is obvious that the sub-sampling procedure must keep the spectrum of the headshape. If there is a spectral overlapping, we may observe a displacement of the minima and therefore a non-optimal registration result. This point justifies the necessity of a uniform headshape sampling as shown above. Moreover, in our context this is not really a constraint but we have to point out that in some other cases (if the headshape is not a continuous surface) a more specific sampling procedure than a simple decimation must be used.

2.4. Inconsistent points rejection

As explained in §2, some parts of the two surfaces to be registered may be unregistrable. This situation may arise from different reasons: for examples residual artefacts in the MRI volume (due to tooth filling for instance) or inconsistent headshape points (technician moves the stylus away from the head). As a result, some data in one set may have no matching part in the other data set and that must necessarily be taken into account in order to avoid the algorithm to converge toward a mathematically good but anatomically wrong solution. For that reason, we apply an inconsistent point rejection procedure during the four last steps of the multiresolution analysis. Bad points are extracted according to a threshold value of the distance transform. A point is kept if it satisfies the following condition:

a may be adapted according to the quality of the headshape. This operation eliminates non-significant points as well as significant points with no match on the MRI surface.

 

2.5. Synthesis

We have implemented this algorithm for the determination of 6 registration parameters (rotation and translation along x, y and z axes). We leave aside a scaling parameter considering the tests performed as described in §3.1.2.2. The minimization is performed along each direction with a fixed sampling step related to the resolution level (see Table I). This approach assumes that the minimum exists in the given interval. This condition is met thanks to a pre-registration task that performs a coarse registration based on the MRI acquisition orientation.

Figure 5 shows the whole registration algorithm. This version will be used in the following evaluations. This algorithm begins with the MRI pre-processing tasks described in §1.2. The initial registration matrix (MatInit) is built from the following information [Lemoine et al 1994]:

• MRI geometrical data (size and resolution) to convert voxel dimension into metric dimension,
• the knowledge of the patient referential system defined during MEG/EEG acquisition,
• the knowledge of the patient head orientation during MRI acquisition.

These last 3 parameters allow a coarse match of both modalities.

The algorithm proceeds in a multiresolution analysis scheme. At each level, the headshape is sub-sampled with a step function related to the resolution. The Powell minimization is performed with an initial value being the result of the previous level (M_Regis matrix is computed). The global registration (matrix) is updated from M_Regis matrix.

Let n be the 2n resolution stage, then the global registration matrix may be expressed in the form:

If necessary, inconsistent points are removed and the procedure is resumed until reaching the last resolution. The registration procedure takes only few seconds on a SUN Sparc 20. MRI pre-processing is more time consuming (about 10 minutes) but is usually performed separately.

 

3. Evaluation

The objective of the evaluation of the registration algorithm is to demonstrate the theoretical and experimental performance (in clinical context) of the method. All the particularities related to the MEG/EEG environment will be included in the evaluation. Indeed, it is necessary to precisely define all possible errors in the process and, if possible, to build simulations in order to study their influence on the global error. If some simulations are unrealisable, the only evaluation is then an evaluation under real conditions. In that case, it is necessary to define the most objective evaluation procedure as possible, knowing that of course the ideal results will remained known.

3.1. Theoretical Evaluation

3.1.1. Strategy

First, we precisely define the type of errors that can be modeled and we will construct the simulations required to independently analyze their contribution to the overall error. Once all errors will be estimated, more complex simulations closer to the reality, of the problem, will be feasible.

As shown on figure 1, one can separate the errors according to the modality (bold arrows). Errors find their origin either during the headshape digitization, or during the MRI acquisition and/or the pre-processing steps. For both modalities we identify two types of errors: local deformations and global deformations.

Concerning MRI and the associated pre-processing, the possible head surface extraction errors were described in §2.1.2.1. Those errors are the scaling effects, the local artefacts and the noise. The scaling effects are inherent to the acquisition itself and to the extraction of the surface. Local artefacts are due mostly to dental prosthesis and eye movements. If the eye movement artefacts can be correctly removed, the dental prosthesis artefact are impossible to remove. They induce a large low signal intensity area in the lower part of the head and thus, alter the outer shape of the mandible area. This artefact may create non existing edges that can mislead the algorithm.

Concerning the headshape digitization, errors can be similarly described. The scaling effect is only due to the skin mechanical properties (Cf. §2.1.1) and, in our case, it will be considered as negligible compared to other errors. Local deformations may be due to the presence of metallic object close to the head during the digitization (e.g. retention device). But most common deformations are due to manipulation errors: either the patient moves the head during the digitization or the technician moves the stylus away from the head. The digitizing device noise is small compared to the noise induced by the skin irregularities or the hair (2 to 3 min) (Cf. §3.1.2.1).

Knowing the different sources of errors, it is easily possible to simulate them in order to understand their contribution to the overall registration process. First, we construct a " theoretical " headshape that perfectly matches the MRI. This theoretical headshape is obtained by uniformly sampling the zero valued points of the distance transform. Known transformations are then applied on this headshape in order to simulate the errors described above and quantify their effects in our algorithm.

The parameters we have selected for the comparison of the theoretical registrations are the Mean Real Error (MRE), its standard deviation, its minimum and maximum values computed on the overall MRI volume (the overall head volume is a 15 cm side cube sampled every 2 mm). The mean distance (Cf. §2.2.2) and its standard deviation will be recorded in order to check that they correctly approximate the quality of the registration.

3.1.2. Results

3.1.2.1. Intrinsic Evaluation of the Algorithm

The intrinsic evaluation of the algorithm is easy. Indeed, one needs only to apply transformations composed of rotations and translations to the theoretical headshape. We have studied the performance of the algorithm for large transformations and for smaller ones (i.e. order of magnitude of the voxel size Cf. §1.1). We have also investigated the influence of the number of points in the headshape.

- Rotation / Translation

We present two types of results that show the capabilities of the algorithm in ideal conditions (the registrations are performed without adding any perturbation (e.g. noise), so the theoretical headshape should be entirely registrable):

Test 1 : (Cf. Figure 6 & Table II) : series of registration results for different initial rotations (0 to 60° with a step of 1°) and initial translations (0 to 60 mm with a step of 1 mm) of the theoretical headshape. The goal of this test is to verify that the algorithm converges within the limit given in its implementation.

Results show that the algorithm converges as long as the transformation does not exceed the minimization interval (See Table I) (MRE = 0.21±0.09 mm). As soon as the transformation exceeds the limits of the implementation defined in Table I, the algorithm does not converge, and therefore does not generate inconsistent registrations.

Test 2 : Series of 300 registrations based on uniformly distributed starting matrices. Those matrices randomly combine rotations in the range: [-19.2° +19.2°] and translations in the range: [-30 mm +30 mm]. The test informs on the general behavior of the algorithm for randomly combined translations and rotations that were not studied by the previous simulations (see table II).

The results of this test are similar to the previous tests. None of the registrations failed. The MRE remained on the order of 1/3 of the voxel size (MRE = 0.37 ± 0.154 mm with a maximum error of the order of magnitude of 1 mm). This test shows that there were no translation/rotation combinations that the algorithm could not handle.

The algorithm does not produce good registrations from a numerical point of view that are inconsistent from an anatomical point of view. Test 1 and Test 2 show that the mean distance is highly correlated to the MRE (see Figure 6.b). The correlation values lie between 0.75 and 0.99 depending on whether missed registrations of test 1 are included (when initial transformations exceeds the minimization interval). The largest correlation is observed when failures are included. This means that bad registrations can be detected, but we can hardly differentiate two registrations with a perfect minimization of the cost function. In such cases the voxel size limit, in term of resolution of the distance transform, is reached. It is also important to notice that the algorithm performance is as good for large as for small initial transformations (see figure 6.a). This property comes from the multiresolution analysis.

-Number of points

This simulation consists in a run of 100 registrations for which the number of points in the headshape varies between 100 and 3000. The initial transformation is identical for all registrations and no additional noise is added. The goal is to test the behavior of the algorithm as a function of the number of points in the headshape. All points are uniformly sampled.

Whatever the number of headshape points is, the MRE remains similar to the two previous tests. The number of points has a negligible effect on a theoretical test. Indeed, all points are perfectly registrable and they are uniformly distributed on the head surface. This test proves the robustness of the algorithm, and the good sampling of the test points on the distance transform. It will be more interesting to observe the behavior of the algorithm with added noise.

 

3.1.2.2. Effects of measurement errors

In order to simulate noise we displace each point of the theoretical headshape using a random function with a uniform distribution. In the rest of this paper, " a noise of 2 mm " means that the amplitude of the displacement is less than, or equal to 2 mm per direction. Thus, we can check that the algorithm is robust to the headshape digitization noise errors and to the MRI acquisition noise. The scale effects are simulated by applying a linear scaling factor to the initial transformation. Finally, local deformations in both modalities are simulated using inconsistent points. A given percentage of the theoretical headshape points are displaced for more than a few centimeters. The random function described above is used to select the points to move.

- Effect of the noise

In order to test the effect of the noise, we chose an initial transformation equivalent to the ones regularly seen in clinic (rotations of 10° and translations of 30 mm). The noise applied to the headshape points varied from 1 to 10 mm, we performed 100 trials per noise level. There was no inconsistent points rejection in this test.

The results (see figure 7) show that the automated registration procedure can handle relatively large noise levels. For a voxel size of 1 mm and up to a noise level of 4 mm per direction, there is no noticeable degradation of the registration. Above that noise level, the algorithm start to fail, specially if it does not converge at the first step of the multiresolution analysis. The mean error is 1.15 ± 0.481 mm and the maximum is 2.388 mm for a noise level of 4 mm.

Note that if the mean distance increases regularly with the noise level (see extended tests in Table III), the distribution of the observed errors do not change in shape. In other words, the noise induces a bias in the cost function; this offset is due to the dispersion of the points in the distance transform.

-Effects of the inconsistent points

We keep the same initial transform as previously, but we insert 15 % of inconsistent points with or without the rejection procedure. When the rejection option is chosen, the rejection parameter a is set to 20.0. We will also study the effects of the rejection threshold: in that case the initial transformation will remain constant.

The addition of inconsistent points provokes repeated failures of the algorithm. This reflects the lack of homologous regions between modalities. As soon as we add a procedure to remove inconsistent points, the algorithm retrieves its efficiency and can even improve it for large noise level. The improvement is due to the rejection of very noisy points which bother the procedure and could cause the algorithm failures in the previous test (see Figure 7). We also observed a slight degradation of the results for low noise levels (1 to 2 mm). In that case the rejection threshold was too low (see below).

We have also studied the behavior of the registration as a function of the inconsistent points rejection parameter a. For a given initial transformation and noise level, we evaluate the precision of the registration for different rejection parameters a between 10 and 40. It is clearly understood that a low rejection threshold may degrade the performances of the registration by removing significant points. For a noise level between 2 and 3 mm, 1 % to 10 % of the total number of points must be rejected, at the last step of the multiresolution procedure. Figure 8 summarizes those results and shows the area where the rejection is significant (central area containing the minimum of the error). Outside this area, the cost function can lose its significance due to the loss of a large amount of points. The value of the cost function is then artificially decreased.

-Scale effects

In order to evaluate the scale effects, we kept used a constant initial transformation (rotations = 10° and translations = 30 mm) on which we added a varying scaling factor (between 0.95 and 1.05). There was no added noise. This will indicate the performance of the algorithm when the scaling factor is not taken into account during the minimization.

The registration remained correct for values of the scaling factor between 0.97 and 1.03. This approximately corresponds to a shift of 3 mm on the skin. A scaling factor becomes necessary when the scaling effect due to the MRI and to the post-processing are greater than 3 mm. In our case, based on theoretical and clinical tests, we have estimated that the scaling factor between our modalities does not require this parameter. Note that the scaling factor can deteriorate the convergence of the algorithm by increasing the number of local minima.

-Effect of the number of points

The conditions of this test are similar to the one described in §3.1.2.1, with a 3 mm noise per direction. We noticed an increase in the quality of the registration as the number of points in the headshape increases from 200 to 2000. Above 2000 points, we reach a plateau (see figure 9). This means that there is minimum number of points required to get all the information necessary for the minimization procedure. Introducing additional points removes ambiguities generated by the noise, limits the number of local minima, and thus increases the stability of the algorithm.

-Extended simulations:

In order to perform simulations as close as possible to the reality, we mix different types of errors. Uncertainty linked to each error is estimated according the knowledge we have of the process that generates them. We estimate the digitization noise, the scaling effect, the number of inconsistent points and the range of parameters (rotations and translations) of the registration matrix generally observed.

We have estimated the value of the uncertainties for every known source of error:

The noise of the Polhemus digitizer is estimated based on different trials:

1) For a given position and direction of the stylus we test the reproducibility of the measures

The mean precision is 0.16 (± 0.07) mm

2) Same operation as for 1), but now an operator is holding the stylus.

The mean precision is 0.2 (± 0.04) mm

3) The stylus is moved between each measure, the operator tries to pick the same point (BTi phantom)

The mean precision is 0.8 (± 0.07) mm

4) Same as for 3) but the point is a dot on someone's skin (see 1.1).

The mean precision is 1.5 (± 0.43) mm

Points 1 and 2 provide an estimation of the precision of the digitizing device. Note that the error due to the Polhemus is equivalent to 1/10 of a millimeter and therefore, its contribution to the digitization can be neglected since the noise linked to the skin (trial 4) is estimated to 2 to 3 mm. The errors estimated in 3) and 4) contribute to the definition of the patient coordinate system (see §1.1). In the discussion of this paper, we will discuss the possible influence of those errors which have no direct effect on the registration itself (See Conclusion).

Scaling effects are principally due to the MRI. If the MRI calibration is good (good field homogeneity at the FOV boundaries), then the size of the voxel size matches the reality with a very slight mismatch. Then, the binarisation operation used in normal conditions do not induce scaling effects. Verifications performed with real registrations showed that no scaling effects (larger than the voxel size) could be observed. Therefore, scaling was not considered in the complete simulations.

The number of inconsistent points was fixed based on the observation of the digitization performed during MEG/EEG examinations and clinical tests. We observed the number of points rejected for each successful registration: generally for headshapes consisting of 3000 to 4000 points, fewer than 1 % of the points were inaccurately digitized (stylus slipping on the head !) and between 1 to 10 % of very noisy points could possibly removed. We have therefore set the inconsistent points ratio to 5 %.

Based on those parameters, we have performed a series of 300 registrations (1 series for a 2 mm noise per direction and 1 series for a noise of 3 mm per direction). The initial registration was randomly built (rotations within [-19°, 19°] and translations within [-30 mm, 30 mm]) (see table III). The rejection parameter a was set to 20.0.

Results show that the mean real error is similar in both cases. This confirms the tests on the noise effects. We notice that the mean distance varies even though the registrations show a similar quality (the test with 3 mm noise per direction shows slightly larger errors but the difference is not significant in a statistical sense). This is the bias mentioned above. It will be necessary to take it into account for the interpretation of real registrations. In clinical conditions, it will be possible to establish a correspondence with those tests when using the verification means discussed in §3.2. In case of good registrations, a high mean distance means a noisy headshape. However, a bad registration with a small mean distance will reflect a non noisy headshape with high real errors.

-Comparison with the manual method

To complete the study, we have performed a theoretical comparison between the theoretical results of the automated registration and the results of the manual registration method that we use when the automated registration is not feasible. The manual registration consists in identifying the same reference coordinate system on both modalities. Therefore, we pick on the MRI volume the same landmarks as the ones pointed during the MEG/EEG exam. The landmarks are vitamin A capsules taped on the patient LPA, RPA and Nasion before the MRI acquisition. Knowing the position of the capsules in the MRI volume makes it possible to compute the patient coordinate system as indicated in §1.1. For manual registrations, errors essentially come from the designation of the landmarks in both modalities. It is also possible to induce large errors if the landmarks pointed in the two modalities are not exactly the same. Theoretical tests are easy to perform. We compute the theoretical coordinate system considering that we have random uncertainties on the position of the landmarks. The errors are set to 1 mm per direction for the MEG/EEG modality and to 2 to 3 mm per direction for the MRI. For each modality, we obtain a set of uncertainties. Since computation of the reference coordinate system, is independent for each modality, we can add them in order to obtain the overall uncertainties. In order to visualize them, we compute a volume in which voxels encode the value of the registration uncertainty. It is possible to superimpose it on an MRI data set in order to view the registration errors according to their localizations within the brain. Results show registration uncertainties between 5 mm and 20 mm. Those uncertainties are far larger and more complex than the automatic method ones [Schwartz et al 1995].

3.2. Evaluation in clinical context

3.2.1. Strategy

Whatever the level of complexity of the theoretical test is, it is impossible to take into account all possible errors and moreover, to correctly simulate their interactions. For example, one can not simulate complex local deformations on the MRI or on the headshape. Also, one can not study the effect of the different shapes of the head. It is also difficult to model the noise induced by small patient motions. It is therefore necessary to define how the performances of the algorithm can be evaluated in real conditions. The definition of an objective mathematical criterion is complex since we don’t know the real registration. The registration error will therefore be an estimate of the real error. In a first approximation, the evaluation can rely on a classification of the final cost function value (mean distance) and on its comparison with theoretical results (Cf. §3.1). The mean distance can indicate a success or a failure of the registration, but it can vary with the noise even for anatomically correct registrations. Thus, we have based our comparison on two additional tests. The first one consists on a visual verification of the registration. The headshape points are superimposed on the MRI and displayed using an orthosection view (axial, sagittal and coronal). This viewing mode makes it very easy to visualize residual errors that could be difficult to perceive on a single 2D or 3D view. The second test is a comparison with the results of the manual registration (see §3.1.2.2) in order to quantify the increase in the precision due to the automated method. We use the following statistical test to compare both registration methods: We compute the registration with both manual and automatic methods. Then we randomly select a set of points (called S) in the headshape (1/10 of the total number of points). For each point in S, an expert evaluates the quality of the registration of the point for both methods (0 = bad registration, 1 = good registration). We obtain two binary distributions M and A: M for the evaluation of manual method, and A for the evaluation of automatic method. We end up with a bilateral sign test between M and A to analyze the quality difference between both registration methods. The signs were obtained by computing = A minus M. The assumption H0 is true when the manual method and the automatic method give similar results (no difference between A and M). Consequently, the assumption H1 is true when, either the quality of the manual registration is inferior to the quality of the automatic registration, or, the quality of the manual registration is superior to the quality of the automatic registration. The set S contains 200 to 400 points. Then we use the central limit theorem to obtain a confidence value assuming a normal distribution:

The quality of the registration of each point is determined by an expert with the following guide-line: for a given point of S the evaluation of the registration quality is given by taking into account the distance of the point to the MRI surface and the quality of the registration of its surrounding points. This neighbor information improves the statistical significance of each point of S. With this evaluation method, a good registration with a noisy headshape is favored over a bad registration with only a part of the headshape well registered.

The clinical evaluation was achieved on subjects who underwent an MEG/EEG examination and for which an MRI was acquired as described in §1.2. The subjects were male and female ranging between 11 and 56 years old, which can be considered as representative to head shape commonly seen (See figure 10). The headshape was digitized by the technician who usually performs the MEG/EEG exams. The digitization protocol requires the following points: good sampling of the nose and the eye contours, uniform sampling of the rest of the skull, no patient head motion. The headshapes contain between 2000 and 4000 points. The MRI pre-processing parameters are set for each type of sequence and usually are not modified from one subject to another. There were a few exceptions for very noisy MRIs. The only varying parameter was the rejection parameter a, since it is linked to the quality of the two modalities. However, it is generally set to 20 and leads to good registrations. For the visual evaluation, MRI volumes are displayed in the orthosection mode (see figure 10). The expert can move the orthosection points to any point into the volume. In a first step the registration is performed with a rejection parameter a set to 20.0. After checking the registration on the display, if the expected quality is not reached, the user can restart the registration with a new rejection parameter a.

In a second step, we perform a manual registration by identifying LPA, RPA and Nasion on MRI. After the registration step we generate the two MRI volumes containing the registered headshape. Pixels belonging to the headshape have distinctive color. During the evaluation the expert didn't know which method was used to perform registration (i.e. manual or automatic).

 

3.2.2. Results

Figures 10, 11 and table IV summarize the results of the comparison of the automatic registration algorithm with the manual registration. Most of the automated registrations studied for this comparison (85 %) were visually considered as accurate (except for LE and PL). Accurate means here that the automatic method provides a better matching of the headshape points with the MRI head surface. Even if the registration is not perfect, there is a better homogeneity of the points position. The failure of case (LE) is explained by failure to respect the headshape acquisition protocol. The face was not correctly sampled and the subject moved during the digitization step. For case (PL), the failure is explained by the bad quality of the MRI: movement of the patient during the acquisition and low contrast. This result in a poor head surface detection from the MRI. The distance transform was therefore less significant. The difference between both methods is less pronounced for three other cases (DF, CD, and SD). For the first on (DF), the automatic registration is not perfect (residual rotations), while the manual registration is good. For the second case (CD) we observe a good registration with the two methods. Finally, for the third case (SD), both methods give excellent results. Two other cases (DMT and DL) required an additional pre-processing of the MRI in order to remove strong dental artefacts. Case (QA) is interesting to some extent. (QA) was 10 years old when undergoing the MEG/EEG examination while MRI was performed 1 year earlier. There was therefore a scaling effect between the two modalities ( 1 mm). However, the algorithm could converge and the registration was accepted as correct. There was no possible manual registration for this case since the MRI procedure did not comply with the MEG/EEG protocol and no capsules were applied.

Those results show that the automatic method improves the accuracy and the robustness of the registration. The method is fast and easy to use and the requirements imposed on the modalities are not drastic. The algorithm is efficient even if there are small scaling effects or if the headshape is noisy (case DMT and DF). Up to now, only one subject could not be registered (PL). If one considers two trials with different values for the rejection parameter a (20 and 30), the registration procedure based on the chamfer volume can be performed in 5 minutes.

 

Conclusion

We have developed a registration method that minimizes manual interventions. This method has been evaluated from both theoretical and clinical points of view. The quality of the registration in terms of accuracy was evaluated for each registration. Using theoretical tests, we found mean registration errors smaller than the voxel size (0.7±0.3 mm). The uncertainties at the boundaries of the MRI volume were estimated between 2 and 3 mm. Those uncertainties are smaller than the commonly observed localization uncertainties. In comparison with the manual method, the improvement in terms of accuracy was significant. By minimizing manual intervention we increased the reliability of the registration process by stabilizing the accuracy. The last manual interaction is the recording of the patient referential during the MEG/EEG acquisition. This has been addressed by using coils to detect landmarks [Fuchs et al 1995] which allows a continuous checking of the patient referential during the acquisition and therefore removes a part of the error due to head motion. Two points need further development. First we have to find a modified cost function which can differentiate between two good registrations. Secondly, we will improve the MRI pre-processing in order to remove small artefacts (e.g. dental artefact, acquisition artefact) which decrease the reliability of the registration. With this automatic registration method, fusion of MEG/EEG localizations on MRI data gives highly significant information for the analyzis of spatio-temporal aspects of dynamic functional processes. Moreover, the accuracy obtained allows the use of complex anatomical constraints (such as accurate description of the cortical surface) in the localization process, without introducing large modelling errors. This algorithm was implemented on the Rennes’ MEG platform. Once the MRI preprocessing parameters are fixed with respect to the acquisition protocol, the successful registration rate is superior to 80% after the first trial. The method is used daily by the medical staff.

Acknowledgments

We thank Pr. P. Toulouse, responsible of the Rennes MEG platform, R.M. Beurel, technician, Pr M. Carsin, head of the radiology department, Pr J.M. Scarabin, head of the neurosurgical department and B. Black from Bti for their valuable contributions in this study.

References

Borgefors, G. Distance Transformations in Digital Images. Computer Vision, Graphics and Image Processing, 1986, 34:344-371.

Bouliou, A., Bihoué, P., Toulouse, P., Poiseau, E., Scarabin, J.M. Modelisation of the cerebral activity by a probility distribution upon the cortical surface from spatio-temporal MEG data. Human Brain Mapping, 1995, (Supp. 1):105.

Chang, H. and Fitzpatick, J.M. A technique for accurate magnetic resonance imaging in the presence of field inhomogeneities. IEEE Trans. Med. Imaging, 1992, 11:319-329.

Cohen, D., Cuffin, B. N., Yunokuchi, K., Maniewski, R., Purcell, C., Cosgrove, G. R., Ives, J., Kennedy, J. and Schoner, D. MEG versus EEG localization test using implanted sources in the human brain. Ann. Neurol., 1990,28:811-817.

Cuffin,N. B. Effects of Head Shape on EEG's and MEG's. IEEE Trans. Biomed. Eng., 1990, 37(1):44-52.

Cuffin,N. B. Effects of Local Variations in Skull and Scalp Thickness on EEG's and MEG's. IEEE Trans. Biomed. Eng., 1993, 40(1):42-48.

Dale, A. M., Sereno, M. I. Improved Localization of Cortical Activity by Combining EEG and MEG with MRI Cortical Surface Reconstruction : A Linear Approach. J. of Cognitive Neuroscience, 1993, 5(2):162-176.

Ebersole, J.S., Squires, K.C., Eliashiv, S.D., Smith, J.R. Applications of Magnetic Source Imaging in Evaluation of Candidates for Epilepsy Surgery. Neuroimaging Clin. of North America, 1995,5(2):267-287

Evans, A.C., Marret, S., Collins, L., Peters, T.M. Anatomical-functional correlative analysis of the human brain using three dimensional imaging systems. In: Schneider R.H., Dwyer III S. J, Gilbert Jost R. (Eds): Proc SPIE Vol 1092, Med Im III: Image Processing, SPIE Press, Bellingham, WA, 1989: 264-274.

Fuchs, M., Wischmann, H.-A., Wagner, M., Krüger, J. Coordinate System for Neuromagnetic and Morphological Reconstruction Overlay. IEEE Trans. Biomed. Eng., 1995, 42(4):416-420.

Gallen, C. C, Schwartz, B., Rieke, K., Pantev, C., Sobel, D., Hirschkoff, E., Bloom, F. E. Intrasubject Reliability and Validity of Somatosensory Source Localization using a large Array Biomagnetometer. Electroenceph. clin. Neurophysiol., 1994, 90:145-156.

Galloway, R.L., Maciunas, R.J., Bass, W.A., Carpini, W.J. Optical localization for interactive image-guided neurosurgery. SPIE Physics of Medical Imaging Newport Beach CA, 1994, 2164:137-146.

Hari, R., Karhu, J., Hämäläinen, M., Knuutila, J., Salonen, O., Sams, M. and Vilkman, V. Functional organization of the human first and second somatosensory cortices: a neuromagnetic study. European Journal of Neuroscience, 1993, 5:724-734.

Hill D.L.G., Hawkes, D.J., Crossman, J.E., Gleeson, M.J., Cox, T.C.S., Bracey, E.C.M.L., Strong, A.J., Graves P. Registration of MRI and CT images for skull base surgery using point-like anatomical features. Brit. J. Radiol., 1991, 35(5):1030-1035.

Ioannides, A.A, Muratore, R., Balish, M., Sato, S. In Vivo Validation of Distributed Source Solutions for the Biomegnetic Inverse Problem. Brain Topography, 1993, 5(3):263-273.

Lemoine, D., Lussot, E., Legeard, D., Barillot, C. Multimodal Registration System for the Fusion of MRI, CT, MEG and 3D or Stereotactic Angiographie Data. SPIE Medical Imaging 1994, 2164:46-56.

Liégois-Chauvel, C., Laguitton, V., Badier, J.M., Schwartz, D. and Chauvel, P. Rev. Neurol., in press, 1995,151.

Paetau, R., Hämäläinen, M., Hari, R., Kajola, M., Karhu, J., Larsen, T.A., Lindahl, E., Salonen, O. Magnetoencephalographic Evaluation of Children And Adolescents with Intractable Epilepsy. Epilepsia, 1994, 35(2):275-284.

Pellizari,C.A., Levin, D.N., Chen, G.T.Y., and Chen, C.T. Image registration based on anatomical surface matching. In: Maciunas, R. (Ed.), Interactive image-guided neurosurgery, American Association of Neurological Surgeons, 1993: 47-62.

Powell, M. J. D. An Efficient Method for Finding the Minimum of a Function of Several Variables Without Calculating Derivatives. The Computer Journal, 1964:155-162.

Press, H. P., Flannery, B. P., Teukolsky, S. A., Vetterling, W. T. Numerical Recipes in C. Cambridge University Press, Cambridge, 1988:309-315.

Polhemus. 3 SPACE® User's Manual. Kaiser aerospace & electronics company, 1991.

Roberts, T., Rowley, H., Kucharczyk, J. Applications of Magnetic Source Imaging to Presurgical Brain Mapping. Neuroimaging Clin. of North America. 1995, 5(2):251-265.

Rohling, R., Munger, P., Hollerbach, J.M., Peters, T. Comparison of relative accuracy between a mechanical and an optical position tracker for image-guided neurosurgery. Journal of Image Guided Surgery, 1995, pre-print.

Scherg, M. Functional imaging and localization of electromagnetic brain activity. Brain Topography, 1992, 5(2):361-364.

Schwartz, D., Poiseau, E., Lemoine, D. MEG/EEG and MRI registration: Errors and uncertainties evaluation. Human Brain Mapping. 1995, (Supp. 1):83

Smith, J. R., Schwartz, B. J., Gallen, C., Orrison, W., Lewine, J., Murro, A. M., King, D. W., Park, Y. D. Utilization of Multichannel Magnetoencephalography in the Guidance of Ablative Seizure Surgery. J. Epilepsy, 1995, 8:119-130.

Weinberg, H., Brickett, P., Coolsma, F., Baff, M. Magnetic Localisation of Intracranial Dipoles : Simulation with a Physical Model. Electroenceph. clin. Neurophysiol., 1986, 64:159-170.

Witwer, J. G., Trezek, G. J., Jewet, Don L. The Effect of Media Inhomogeneities Upon Intracranial Electrical Fields. IEEE Trans. Biomed. Eng., 1972, 19(5):352-362.

West, J., Fittzpatrick, M., Wang, M.Y., Dawant, M.B., Maurer, C.R., Kessler, R.M., Maciunas, R.J., Barillot, C., Lemoine, D., Collignon, A., Maes, F., Suetens, P., Vandermeulen, D., van den Elsen, P., Hemler, P.F., Napel, T., Sumanaweera, T.S., Harkness, B., Hill, D.L.G., Studholme, C., Malandain, G., Pennec, X., Noz, M.E., Maguire, G.Q., Pollack, M., Pelizzari, C.A., Robb, R.A., Hanson, D., Woods, R.P. Comparison and evaluation of retrospective intermodality image registration techniques. Will appear in : Medical Imaging 1996: Image Processing, SPIE Proc. 1996, 2710.

Wood, C. C., Cohen, D., Cuffin, B. N., & Allison, T. Electrical Sources in Human Somataosensory Cortex : Identification by Combined Magnetic and Potential Recordings. Science, 1985, 227:1051- 1053.

Young, S. T., Yip, S.W., Cheng, H.C., Shieh, D. Three-Dimensional Surface Digitizer for facial contour. IEEE Engineering in Medecine and Biology, 1994, 125:141.

Figures captions

Figure 1: Framework of data fusion between MEG/EEG data and MRI data. The bold arrows represent errors propagation pathways investigated in this paper.

Figures :

k = 2	k =1	k = 0

Tables Captions

Table I : Rotation and translation ranges and search steps during minimization for each resolution level.

Table II : Results of intrinsic tests 1 and 2 for initial transformation of the theoretical headshape in the minimization interval range (See Table I). The first column shows the mean of the Mean Real Error (MRE) for each test. The second column represents the maximum of MRE for each test. And the third column represents the mean of the final Mean Distance (cost function value) for each test. The last column shows the number of trials per test.

 

Table III : Results of extended simulations for two noise levels (2 mm and 3 mm per direction). The first column shows the mean of the Mean Real Error (MRE) for each noise level. The second column shows the maximum of MRE for each noise level.. The third column shows the mean of the final Mean Distance (cost function value) for each noise level. The last column shows the number of trials per test

 

Table IV: Results of clinical evaluation. The second column shows the final Mean Distance (cost function value) for each registration. The column "Orthosection Check" shows the expert evaluation of the registration based on orthosection display (+++: perfect, ++: good with few mis-registered points of the headshape, +: a set of points of the headshape was slightly mis-registered, -: a set of points of the headshape was strongly mis-registered, --: mis-registration, registration useless in a clinical context, ---: the automatic algorithm does not converge). Column "Sign Test": if the sign test is greater than 2.53 the automatic method is considered more efficient, if the sign test is smaller than -2.53 the automatic registration is considered less efficient and if the sign test belonged to -2.53 to 2.53 range the both methods give similar results. *: Epileptic patients.

Tables: