[BUG, MRG] Fix ACPC Coordinates in surface RAS when they should be in…

… scanner RAS (#990)

* wip, helper functions don't work

* working version

* add PR number

* typo

* fix errors

* eric review

doc/api.rst

@@ -31,6 +31,8 @@ mne_bids
    get_anat_landmarks
    update_anat_landmarks
    get_head_mri_trans
+   convert_montage_to_mri
+   convert_montage_to_ras
    template_to_head
    get_anonymization_daysback
    search_folder_for_text

doc/whats_new.rst

@@ -57,7 +57,7 @@ Detailed list of changes
 🪲 Bug fixes
 ^^^^^^^^^^^^
 
-- ...
+- Fix ACPC in ``surface RAS`` instead of ``scanner RAS`` in :ref:`ieeg-example` and add convenience functions :func:`mne_bids.convert_montage_to_ras` and :func:`mne_bids.convert_montage_to_mri` to help, by `Alex Rockhill`_ :gh:`990`
 
 :doc:`Find out what was new in previous releases <whats_new_previous_releases>`
 

examples/convert_ieeg_to_bids.py

@@ -60,7 +60,8 @@
 from mne_bids import (BIDSPath, write_raw_bids, write_anat,
                       get_anat_landmarks, read_raw_bids,
                       search_folder_for_text, print_dir_tree,
-                      template_to_head)
+                      template_to_head, convert_montage_to_ras,
+                      convert_montage_to_mri)
 
 # %%
 # Step 1: Download the data
@@ -98,9 +99,10 @@
 trans = mne.coreg.estimate_head_mri_t('sample_seeg', subjects_dir)
 
 # %%
-# Now let's convert the montage to "mri"
+# Now let's convert the montage to "ras"
 montage = raw.get_montage()
 montage.apply_trans(trans)  # head->mri
+convert_montage_to_ras(montage, 'sample_seeg', subjects_dir)  # mri->ras
 
 # %%
 # BIDS vs MNE-Python Coordinate Systems
@@ -237,13 +239,16 @@
 raw2 = read_raw_bids(bids_path=bids_path)
 
 # %%
-# Now we have to go back to "head" coordinates with the head->mri transform.
+# Now we have to go back to "head" coordinates.
 #
 # .. note:: If you were downloading this from ``OpenNeuro``, you would
 #           have to run the Freesurfer ``recon-all`` to get the transforms.
 
 montage2 = raw2.get_montage()
 
+# we need to go from scanner RAS back to surface RAS (requires recon-all)
+convert_montage_to_mri(montage2, 'sample_seeg', subjects_dir=subjects_dir)
+
 # this uses Freesurfer recon-all subject directory
 montage2.add_estimated_fiducials('sample_seeg', subjects_dir=subjects_dir)
 
@@ -259,7 +264,9 @@
 # the fiducials which are estimated; as long as the head->mri trans
 # is computed with the same fiducials, the coordinates will be the same
 # in ACPC space which is what matters
-montage2 = raw.get_montage()  # get montage in 'head' coordinates
+montage = raw.get_montage()  # the original montage in 'head' coordinates
+montage.apply_trans(trans)
+montage2 = raw2.get_montage()  # the recovered montage in 'head' coordinates
 montage2.apply_trans(trans2)
 
 # compare with standard

mne_bids/__init__.py

@@ -14,4 +14,5 @@
                             get_anat_landmarks, anonymize_dataset)
 from mne_bids.sidecar_updates import update_sidecar_json, update_anat_landmarks
 from mne_bids.inspect import inspect_dataset
-from mne_bids.dig import template_to_head
+from mne_bids.dig import (template_to_head, convert_montage_to_ras,
+                          convert_montage_to_mri)

mne_bids/config.py

@@ -254,7 +254,7 @@
     'ElektaNeuromag': 'head',
     'ChietiItab': 'head',
     'CapTrak': 'head',
-    'ACPC': 'mri',  # assumes T1 is ACPC-aligned, if not the coordinates are lost  # noqa
+    'ACPC': 'ras',  # assumes T1 is ACPC-aligned, if not the coordinates are lost  # noqa
     'fsaverage': 'mni_tal',  # XXX: note fsaverage and MNI305 are the same  # noqa
     'MNI305': 'mni_tal'
 }

mne_bids/dig.py

@@ -4,6 +4,7 @@
 #          Alex Rockhill <aprockhill@mailbox.org>
 #
 # License: BSD-3-Clause
+import os.path as op
 import json
 from collections import OrderedDict
 from pathlib import Path
@@ -14,7 +15,8 @@
 import mne
 import numpy as np
 from mne.io.constants import FIFF
-from mne.utils import logger, warn, _validate_type, _check_option
+from mne.utils import (logger, warn, _validate_type, _check_option,
+                       has_nibabel, get_subjects_dir)
 from mne.io.pick import _picks_to_idx
 
 from mne_bids.config import (ALLOWED_SPACES, BIDS_COORDINATE_UNITS,
@@ -416,7 +418,8 @@ def _write_dig_bids(bids_path, raw, montage=None, acpc_aligned=False,
                                "and left and right pre-auricular point "
                                "landmarks")
 
-    if bids_path.datatype == 'ieeg' and mne_coord_frame == 'mri':
+    if bids_path.datatype == 'ieeg' and bids_path.space in (None, 'ACPC') and \
+            mne_coord_frame == 'ras':
         if not acpc_aligned:
             raise RuntimeError(
                 '`acpc_aligned` is False, if your T1 is not aligned '
@@ -656,3 +659,87 @@ def template_to_head(info, space, coord_frame='auto', unit='auto',
     info.set_montage(montage)  # transform to head
     # finally return montage
     return info, mne.read_trans(data_dir / f'space-{space}_trans.fif')
+
+
+@verbose
+def convert_montage_to_ras(montage, subject, subjects_dir=None, verbose=None):
+    """Convert a montage from surface RAS (m) to scanner RAS (m).
+
+    Parameters
+    ----------
+    montage : mne.channels.DigMontage
+        The montage in the "mri" coordinate frame. Note: modified in place.
+    %(subject)s
+    %(subjects_dir)s
+    %(verbose)s
+    """
+    if not has_nibabel():  # pragma: no cover
+        raise ImportError('This function requires nibabel.')
+    import nibabel as nib
+
+    subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)
+    T1_fname = op.join(subjects_dir, subject, 'mri', 'T1.mgz')
+    if not op.isfile(T1_fname):
+        raise RuntimeError(f'Freesurfer subject ({subject}) and/or '
+                           f'subjects_dir ({subjects_dir}, incorrectly '
+                           'formatted, T1.mgz not found')
+    T1 = nib.load(T1_fname)
+
+    # transform from "mri" (Freesurfer surface RAS) to "ras" (scanner RAS)
+    mri_vox_t = np.linalg.inv(T1.header.get_vox2ras_tkr())
+    mri_vox_t[:3, :3] *= 1000  # scale from mm to m
+    mri_vox_trans = mne.transforms.Transform(
+        fro='mri', to='mri_voxel', trans=mri_vox_t)
+
+    vox_ras_t = T1.header.get_vox2ras()
+    vox_ras_t[:3] /= 1000  # scale from mm to m
+    vox_ras_trans = mne.transforms.Transform(
+        fro='mri_voxel', to='ras', trans=vox_ras_t)
+    montage.apply_trans(  # mri->vox + vox->ras = mri->ras
+        mne.transforms.combine_transforms(mri_vox_trans, vox_ras_trans,
+                                          fro='mri', to='ras'))
+
+
+@verbose
+def convert_montage_to_mri(montage, subject, subjects_dir=None, verbose=None):
+    """Convert a montage from scanner RAS (m) to surface RAS (m).
+
+    Parameters
+    ----------
+    montage : mne.channels.DigMontage
+        The montage in the "ras" coordinate frame. Note: modified in place.
+    %(subject)s
+    %(subjects_dir)s
+    %(verbose)s
+
+    Returns
+    -------
+    ras_mri_t : mne.transforms.Transform
+        The transformation matrix from ``'ras'`` (``scanner RAS``) to
+        ``'mri'`` (``surface RAS``).
+    """
+    if not has_nibabel():  # pragma: no cover
+        raise ImportError('This function requires nibabel.')
+    import nibabel as nib
+
+    subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)
+    T1_fname = op.join(subjects_dir, subject, 'mri', 'T1.mgz')
+    if not op.isfile(T1_fname):
+        raise RuntimeError(f'Freesurfer subject ({subject}) and/or '
+                           f'subjects_dir ({subjects_dir}, incorrectly '
+                           'formatted, T1.mgz not found')
+    T1 = nib.load(T1_fname)
+
+    # transform from "ras" (scanner RAS) to "mri" (Freesurfer surface RAS)
+    ras_vox_t = T1.header.get_ras2vox()
+    ras_vox_t[:3, :3] *= 1000  # scale from mm to m
+    ras_vox_trans = mne.transforms.Transform(
+        fro='ras', to='mri_voxel', trans=ras_vox_t)
+
+    vox_mri_t = T1.header.get_vox2ras_tkr()
+    vox_mri_t[:3] /= 1000  # scale from mm to m
+    vox_mri_trans = mne.transforms.Transform(
+        fro='mri_voxel', to='mri', trans=vox_mri_t)
+    montage.apply_trans(  # ras->vox + vox->mri = ras->mri
+        mne.transforms.combine_transforms(ras_vox_trans, vox_mri_trans,
+                                          fro='ras', to='mri'))