Using an automated approach to coregistration

This example shows how to use the coregistration functions to perform an automated MEG-MRI coregistration via scripting.

Warning

The quality of the coregistration depends heavily upon the quality of the head shape points (HSP) collected during subject prepration and the quality of your T1-weighted MRI. Use with caution and check the coregistration error.

# Author: Jon Houck <jon.houck@gmail.com>
#         Guillaume Favelier <guillaume.favelier@gmail.com>
#
# License: BSD-3-Clause

import os.path as op
import numpy as np

import mne
from mne.coreg import Coregistration
from mne.io import read_info


data_path = mne.datasets.sample.data_path()
subjects_dir = op.join(data_path, 'subjects')
subject = 'sample'

fname_raw = op.join(data_path, 'MEG', subject, subject + '_audvis_raw.fif')
info = read_info(fname_raw)
plot_kwargs = dict(subject=subject, subjects_dir=subjects_dir,
                   surfaces="head-dense", dig=True, eeg=[],
                   meg='sensors', show_axes=True,
                   coord_frame='meg')
view_kwargs = dict(azimuth=45, elevation=90, distance=0.6,
                   focalpoint=(0., 0., 0.))

Out:

Read a total of 3 projection items:
    PCA-v1 (1 x 102)  idle
    PCA-v2 (1 x 102)  idle
    PCA-v3 (1 x 102)  idle

Set up the coregistration model

fiducials = "estimated"  # get fiducials from fsaverage
coreg = Coregistration(info, subject, subjects_dir, fiducials=fiducials)
fig = mne.viz.plot_alignment(info, trans=coreg.trans, **plot_kwargs)

Out:

    Triangle neighbors and vertex normals...
Using high resolution head model in /home/circleci/mne_data/MNE-sample-data/subjects/sample/surf/lh.seghead
    Triangle neighbors and vertex normals...
Estimating fiducials from fsaverage.
Using lh.seghead for head surface.
Channel types:: grad: 203, mag: 102

Do first a coregistration fit using only 3 fiducial points. This allows to find a good initial solution before further optimization using head shape points. This can also be useful to detect outlier head shape points which are too far from the skin surface. One can see for example that on this dataset there is one such point and we will omit it from the subsequent fit.

coreg.fit_fiducials(verbose=True)
fig = mne.viz.plot_alignment(info, trans=coreg.trans, **plot_kwargs)

Out:

Aligning using fiducials
Start median distance:  17.29 mm
End   median distance:  11.11 mm
Using lh.seghead for head surface.
Channel types:: grad: 203, mag: 102

Next we refine the transformation using a few iteration of the Iterative Closest Point (ICP) algorithm. As the initial fiducials are obtained from fsaverage and not from precise manual picking in the GUI we do a fit with reduced weight for the nasion.

coreg.fit_icp(n_iterations=6, nasion_weight=2., verbose=True)
fig = mne.viz.plot_alignment(info, trans=coreg.trans, **plot_kwargs)

Out:

Aligning using ICP
Start     median distance:  11.11 mm
  ICP  1  median distance:   6.57 mm
  ICP  2  median distance:   4.39 mm
  ICP  3  median distance:   3.48 mm
  ICP  4  median distance:   2.86 mm
  ICP  5  median distance:   2.43 mm
  ICP  6  median distance:   2.14 mm
End       median distance:   2.14 mm
Using lh.seghead for head surface.
Channel types:: grad: 203, mag: 102

It is now very clear that we have one point that is an outlier and that should be removed.

coreg.omit_head_shape_points(distance=5. / 1000)  # distance is in meters

Out:

Coregistration: Excluding 3 head shape points with distance >= 0.005 m.

Do a final coregistration fit

coreg.fit_icp(n_iterations=20, nasion_weight=10., verbose=True)
fig = mne.viz.plot_alignment(info, trans=coreg.trans, **plot_kwargs)
mne.viz.set_3d_view(fig, **view_kwargs)

dists = coreg.compute_dig_mri_distances() * 1e3  # in mm
print(
    f"Distance between HSP and MRI (mean/min/max):\n{np.mean(dists):.2f} mm "
    f"/ {np.min(dists):.2f} mm / {np.max(dists):.2f} mm"
)

Out:

Aligning using ICP
Start     median distance:   2.00 mm
  ICP  1  median distance:   1.99 mm
  ICP  2  median distance:   2.02 mm
  ICP  3  median distance:   1.79 mm
  ICP  4  median distance:   1.78 mm
  ICP  5  median distance:   1.77 mm
End       median distance:   1.77 mm
Using lh.seghead for head surface.
Channel types:: grad: 203, mag: 102
Distance between HSP and MRI (mean/min/max):
2.00 mm / 0.19 mm / 6.37 mm

Warning

Don’t forget to save the resulting trans matrix!

mne.write_trans('/path/to/filename-trans.fif', coreg.trans)

Note

The mne.coreg.Coregistration class has the ability to compute MRI scale factors using set_scale_mode() that is useful for creating surrogate MRI subjects, i.e., using a template MRI (such as one from mne.datasets.fetch_infant_template()) matched to a subject’s head digitization. When scaling is desired, a scaled surrogate MRI should be created using mne.scale_mri().

Estimated memory usage: 36 MB

previous

Source alignment and coordinate frames

next

Head model and forward computation