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

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 convence 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 <[email protected]>
 #
 # 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,7 @@ 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 mne_coord_frame == 'ras':
         if not acpc_aligned:
             raise RuntimeError(
                 '`acpc_aligned` is False, if your T1 is not aligned '
@@ -656,3 +658,107 @@ 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)
+
+    # scale from mm to m
+    scale_t = np.eye(4)
+    scale_t[:3, :3] *= 1000
+    scale_trans = mne.transforms.Transform(fro='mri', to='mri', trans=scale_t)
+    montage.apply_trans(scale_trans)
+
+    # transform from "mri" (Freesurfer surface RAS) to "ras" (scanner RAS)
+    mri_vox_t = np.linalg.inv(T1.header.get_vox2ras_tkr())
+    mri_vox_trans = mne.transforms.Transform(
+        fro='mri', to='mri_voxel', trans=mri_vox_t)
+    montage.apply_trans(mri_vox_trans)  # mri->vox
+
+    vox_ras_t = T1.header.get_vox2ras()
+    vox_ras_trans = mne.transforms.Transform(
+        fro='mri_voxel', to='ras', trans=vox_ras_t)
+    montage.apply_trans(vox_ras_trans)  # vox->ras
+
+    # finally, need to put back in m
+    scale_inv_t = np.eye(4)
+    scale_inv_t[:3, :3] /= 1000
+    scale_inv_trans = mne.transforms.Transform(
+        fro='ras', to='ras', trans=scale_inv_t)
+    montage.apply_trans(scale_inv_trans)
+
+
+@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)
+
+    # scale from mm to m
+    scale_t = np.eye(4)
+    scale_t[:3, :3] *= 1000
+    scale_trans = mne.transforms.Transform(fro='ras', to='ras', trans=scale_t)
+    montage.apply_trans(scale_trans)
+
+    # transform from "ras" (scanner RAS) to "mri" (Freesurfer surface RAS)
+    ras_vox_t = T1.header.get_ras2vox()
+    ras_vox_trans = mne.transforms.Transform(
+        fro='ras', to='mri_voxel', trans=ras_vox_t)
+    montage.apply_trans(ras_vox_trans)  # ras->vox
+
+    vox_mri_t = T1.header.get_vox2ras_tkr()
+    vox_mri_trans = mne.transforms.Transform(
+        fro='mri_voxel', to='mri', trans=vox_mri_t)
+    montage.apply_trans(vox_mri_trans)  # vox->mri
+
+    # finally, need to put back in m
+    scale_inv_t = np.eye(4)
+    scale_inv_t[:3, :3] /= 1000
+    scale_inv_trans = mne.transforms.Transform(
+        fro='mri', to='mri', trans=scale_inv_t)
+    montage.apply_trans(scale_inv_trans)

mne_bids/tests/test_dig.py

@@ -17,7 +17,8 @@
 import mne_bids
 from mne.datasets import testing
 from mne_bids import BIDSPath, write_raw_bids, read_raw_bids
-from mne_bids.dig import _write_dig_bids, _read_dig_bids, template_to_head
+from mne_bids.dig import (_write_dig_bids, _read_dig_bids, template_to_head,
+                          convert_montage_to_mri, convert_montage_to_ras)
 from mne_bids.config import (BIDS_STANDARD_TEMPLATE_COORDINATE_SYSTEMS,
                              BIDS_TO_MNE_FRAMES, MNE_STR_TO_FRAME)
 
@@ -36,6 +37,8 @@
 data_path = testing.data_path()
 raw_fname = op.join(data_path, 'MEG', 'sample',
                     'sample_audvis_trunc_raw.fif')
+trans = mne.read_trans(op.join(data_path, 'MEG', 'sample',
+                               'sample_audvis_trunc-trans.fif'))
 raw = mne.io.read_raw(raw_fname)
 raw.drop_channels(raw.info['bads'])
 raw.info['line_freq'] = 60
@@ -284,3 +287,33 @@ def test_template_to_head():
         d['r'] = np.array(d['r']) / 1000
     _test_montage_trans(raw_test, montage_m, pos_test,
                         coord_frame='ras', unit='auto')
+
+
+def test_convert_montage():
+    """Test the montage RAS conversion."""
+    raw_test = raw.copy()
+    montage = raw_test.get_montage()
+    montage.apply_trans(trans)
+
+    subjects_dir = op.join(data_path, 'subjects')
+    # test read
+    with pytest.raises(RuntimeError, match='incorrectly formatted'):
+        convert_montage_to_mri(montage, 'foo', subjects_dir)
+
+    # test write
+    with pytest.raises(RuntimeError, match='incorrectly formatted'):
+        convert_montage_to_ras(montage, 'foo', subjects_dir)
+
+    # test mri to ras
+    convert_montage_to_ras(montage, 'sample', subjects_dir)
+    pos = montage.get_positions()
+    assert pos['coord_frame'] == 'ras'
+    assert_almost_equal(pos['ch_pos']['EEG 001'],
+                        [-0.0366405, 0.063066, 0.0676311])
+
+    # test ras to mri
+    convert_montage_to_mri(montage, 'sample', subjects_dir)
+    pos = montage.get_positions()
+    assert pos['coord_frame'] == 'mri'
+    assert_almost_equal(pos['ch_pos']['EEG 001'],
+                        [-0.0313669, 0.0540269, 0.0949191])

mne_bids/tests/test_read.py

@@ -22,8 +22,7 @@
 from mne_bids import BIDSPath
 from mne_bids.config import (MNE_STR_TO_FRAME, BIDS_SHARED_COORDINATE_FRAMES,
                              BIDS_TO_MNE_FRAMES)
-from mne_bids.read import (read_raw_bids,
-                           _read_raw, get_head_mri_trans,
+from mne_bids.read import (read_raw_bids, _read_raw, get_head_mri_trans,
                            _handle_events_reading)
 from mne_bids.tsv_handler import _to_tsv, _from_tsv
 from mne_bids.utils import (_write_json)
@@ -881,7 +880,7 @@ def test_handle_ieeg_coords_reading(bids_path, tmp_path):
     # ACPC should be read in as RAS for iEEG
     _update_sidecar(coordsystem_fname, 'iEEGCoordinateSystem', 'ACPC')
     raw_test = read_raw_bids(bids_path=bids_fname, verbose=False)
-    coord_frame_int = MNE_STR_TO_FRAME['mri']
+    coord_frame_int = MNE_STR_TO_FRAME['ras']
     for digpoint in raw_test.info['dig']:
         assert digpoint['coord_frame'] == coord_frame_int
 

mne_bids/tests/test_write.py

@@ -1507,7 +1507,7 @@ def test_eegieeg(dir_name, fname, reader, _bids_validate, tmp_path):
 
     # test writing to ACPC
     ecog_montage = mne.channels.make_dig_montage(ch_pos=ch_pos,
-                                                 coord_frame='mri')
+                                                 coord_frame='ras')
     bids_root = tmp_path / 'bids4'
     bids_path.update(root=bids_root, datatype='ieeg')
     # test works if ACPC-aligned is specified