mne-tools · jhouck · Aug 21, 2019 · Aug 21, 2019 · Aug 21, 2019 · Aug 23, 2019
diff --git a/examples/forward/plot_auto_alignment.py b/examples/forward/plot_auto_alignment.py
@@ -0,0 +1,80 @@
+"""
+=========================================
+Use automated approach to co-registration
+=========================================
+
+This example shows how to use the get_mni_fiducials routine and the
+coregistration GUI functions to perform an automated MEG-MRI co-registration.
+
+.. warning:: The quality of the co-registration depends heavily upon the
+             quality of the head shape collected during subject prepration and
+             the quality of your T1-weighted MRI. Use with caution and check
+             the co-registration error.
+"""
+
+# License: BSD (3-clause)
+
+import os.path as op
+import numpy as np
+import mne
+from mne.coreg import get_mni_fiducials
+from mne.surface import dig_mri_distances
+from mne.io import write_fiducials
+from mne.io.constants import FIFF
+from mne.gui._file_traits import DigSource
+from mne.gui._fiducials_gui import MRIHeadWithFiducialsModel
+from mne.gui._coreg_gui import CoregModel
+
+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')
+fname_fids = op.join(subjects_dir, subject, 'bem', subject + '-fiducials.fif')
+fname_trans = op.join(data_path, 'MEG', subject,
+                      subject + 'audvis_raw_auto-trans.fif')
+src = mne.read_source_spaces(op.join(subjects_dir, subject, 'bem',
+                                     'sample-oct-6-src.fif'))
+
+fids_mri = get_mni_fiducials(subject, subjects_dir)
+
+# Save mri fiducials. This is mandatory. as fit_fiducials uses this file
+
+write_fiducials(fname_fids, fids_mri, coord_frame=FIFF.FIFFV_COORD_MRI)
+print('\nSaving estimated fiducials to %s \n' % fname_fids)
+
+# set up HSP DigSource
+hsp = DigSource()
+hsp.file = fname_raw
+
+# Set up subject MRI source space with fiducials
+mri = MRIHeadWithFiducialsModel(subjects_dir=subjects_dir, subject=subject)
+
+# Set up coreg model
+model = CoregModel(mri=mri, hsp=hsp)
+# Do initial fit to fiducials
+model.fit_fiducials()
+# Do initial coreg fit for outlier detection
+model.icp_iterations = int(6)
+model.nasion_weight = 2.  # For this fit we know the nasion is not precise
+# Overweighting at this step also seems to throw off the fit for some datasets
+model.fit_icp()
+model.omit_hsp_points(distance=5. / 1000)  # Distance is in meters
+# Do final coreg fit
+model.nasion_weight = 10.
+model.icp_iterations = int(20)
+model.fit_icp()
+model.save_trans(fname=fname_trans)
+errs_icp = model._get_point_distance()
+raw = mne.io.Raw(fname_raw)
+errs_nearest = dig_mri_distances(raw.info, fname_trans, subject, subjects_dir)
+
+fig = mne.viz.plot_alignment(raw.info, trans=fname_trans, subject=subject,
+                             subjects_dir=subjects_dir, surfaces='head-dense',
+                             dig=True, eeg=[], meg='sensors',
+                             coord_frame='meg')
+mne.viz.set_3d_view(fig, 45, 90, distance=0.6, focalpoint=(0., 0., 0.))
+
+print('Median distance from digitized points to head surface is %.3f mm'
+      % np.median(errs_icp * 1000))
+print('''Median distance from digitized points to head surface using nearest
+neighbor is %.3f mm''' % np.median(errs_nearest * 1000))
diff --git a/mne/coreg.py b/mne/coreg.py
@@ -1206,3 +1206,72 @@ def _scale_xfm(subject_to, xfm_fname, mri_name, subject_from, scale,
                 F_mri_ras, 'ras', 'ras'),
             F_ras_mni, 'ras', 'mni_tal')
     _write_fs_xfm(fname_to, T_ras_mni['trans'], kind)
+
+
+def get_mni_fiducials(subject, subjects_dir):
+    """Estimate fiducials for a subject.
+
+    Parameters
+    ----------
+    subject : str
+        Name of the mri subject
+    subjects_dir : None | str
+        Override the SUBJECTS_DIR environment variable
+        (sys.environ['SUBJECTS_DIR'])
+
+    Returns
+    -------
+    fids_mri : list
+        List of estimated fiducials (each point in a dict)
+
+    Notes
+    -----
+    For more details about the coordinate systems and transformations involved,
+    see https://surfer.nmr.mgh.harvard.edu/fswiki/CoordinateSystems and
+    :ref:`plot_source_alignment`
+    """
+    from .transforms import _read_fs_xfm
+    if not has_nibabel():
+        raise ImportError('This function requires nibabel.')
+        return
+
+    from nibabel import freesurfer as fs
+
+    subjects_dir = op.join(get_subjects_dir(subjects_dir, raise_error=True),
+                           subject)
+    fname_orig = op.join(subjects_dir, 'mri', 'T1.mgz')
+    fname_tal = op.join(subjects_dir, 'mri', 'transforms', 'talairach.xfm')
+    fname_fids_fs = os.path.join(os.path.dirname(__file__), 'data',
+                                 'fsaverage', 'fsaverage-fiducials.fif')
+
+    # Read fsaverage fiducials file and subject Talairach.
+    fids_default, coord_frame = read_fiducials(fname_fids_fs)
+    xfm_tal = np.matrix(_read_fs_xfm(fname_tal)[0])
+
+    # Get Freesurfer vox2ras and vox2ras-tkr as matrices.
+    mgh = fs.mghformat.load(fname_orig)
+    vox2ras = np.matrix(mgh.header.get_vox2ras())
+    vox2ras_tkr = np.matrix(mgh.header.get_vox2ras_tkr())
+
+    # Get transform for RAS to MNI305.
+    mni305ras = xfm_tal * vox2ras * np.linalg.inv(vox2ras_tkr)
+    mni2ras = np.linalg.inv(mni305ras)
+
+    # Split up fiducials. Convert to mm since this is Freesurfer's unit.
+    lpa = np.append(fids_default[0]['r'] * 1000, 1)
+    nas = np.append(fids_default[1]['r'] * 1000, 1)
+    rpa = np.append(fids_default[2]['r'] * 1000, 1)
+
+    # Apply transformation, to fsaverage (MNI) fiducials, convert from mm.
+    lpa_ras = np.delete(np.dot(mni2ras, lpa.T), 3) / 1000
+    nas_ras = np.delete(np.dot(mni2ras, nas.T), 3) / 1000
+    rpa_ras = np.delete(np.dot(mni2ras, rpa.T), 3) / 1000
+
+    # Build new fiducials structure from old.
+    fids_mri = fids_default
+    # Copy to fiducials.
+    fids_mri[0]['r'] = lpa_ras
+    fids_mri[1]['r'] = nas_ras
+    fids_mri[2]['r'] = rpa_ras
+
+    return fids_mri