Using an automated approach to coregistration

This example shows how to use the coregistration functions to perform an automated MEG-MRI coregistration via scripting. Generally the results of this approach are consistent with those obtained from manual coregistration [1].

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 numpy as np

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

data_path = mne.datasets.sample.data_path()
# data_path and all paths built from it are pathlib.Path objects
subjects_dir = data_path / "subjects"
subject = "sample"

fname_raw = data_path / "MEG" / subject / f"{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.0, 0.0))

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)

    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

Initial fit with fiducials

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)

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

Refining with ICP

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.0, verbose=True)
fig = mne.viz.plot_alignment(info, trans=coreg.trans, **plot_kwargs)

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

Omitting bad points

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.0 / 1000)  # distance is in meters

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

Final coregistration fit

coreg.fit_icp(n_iterations=20, nasion_weight=10.0, 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"
)

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().

References

[1]

Jon M. Houck and Eric D. Claus. A comparison of automated and manual co-registration for magnetoencephalography. PLOS ONE, 15:e0232100, 2020. doi:10.1371/journal.pone.0232100.

Estimated memory usage: 53 MB