[ENH, MRG] Add BIDS Standard Template Transformations

doc/api.rst

@@ -31,6 +31,7 @@ mne_bids
    get_anat_landmarks
    update_anat_landmarks
    get_head_mri_trans
+   template_to_head
    get_anonymization_daysback
    search_folder_for_text
    print_dir_tree

doc/whats_new.rst

@@ -52,12 +52,14 @@ Enhancements
 
 - Add an explanation in :ref:`ieeg-example` of why it is better to have intracranial data in individual rather than template coordinates, by `Alex Rockhill`_ (:gh:`975`)
 
-- :func:`mne_bids.update_anat_landmarks` can now directly work with fiducials saved from the MNE-Python coregistration GUI or :func:`mne.io.write_fiducials`, by Richard Höchenberger`_ (:gh:`977`)
+- :func:`mne_bids.update_anat_landmarks` can now directly work with fiducials saved from the MNE-Python coregistration GUI or :func:`mne.io.write_fiducials`, by `Richard Höchenberger`_ (:gh:`977`)
 
 - All non-MNE-Python BIDS coordinate frames are now set to ``'unknown'`` on reading, by `Alex Rockhill`_ (:gh:`979`)
 
 - :func:`mne_bids.write_raw_bids` can now write to template coordinates by `Alex Rockhill`_ (:gh:`980`)
 
+- Add :func:`mne_bids.template_to_head` to transform channel locations in BIDS standard template coordinate systems to ``head`` and also provides a ``trans``, by `Alex Rockhill`_ (:gh:`983`)
+
 API and behavior changes
 ^^^^^^^^^^^^^^^^^^^^^^^^
 

mne_bids/__init__.py

@@ -14,3 +14,4 @@
                             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

mne_bids/config.py

@@ -165,7 +165,7 @@
 # Annotations to never remove during reading or writing
 ANNOTATIONS_TO_KEEP = ('BAD_ACQ_SKIP',)
 
-BIDS_STANDARD_TEMPLATE_COORDINATE_FRAMES = [
+BIDS_STANDARD_TEMPLATE_COORDINATE_SYSTEMS = [
     'ICBM452AirSpace',
     'ICBM452Warp5Space',
     'IXI549Space',
@@ -218,7 +218,7 @@
 # accepted coordinate SI units
 BIDS_COORDINATE_UNITS = ['m', 'cm', 'mm']
 coordsys_wildcard = ['Other']
-BIDS_SHARED_COORDINATE_FRAMES = (BIDS_STANDARD_TEMPLATE_COORDINATE_FRAMES +
+BIDS_SHARED_COORDINATE_FRAMES = (BIDS_STANDARD_TEMPLATE_COORDINATE_SYSTEMS +
                                  coordsys_standard_template_deprecated +
                                  coordsys_wildcard)
 

mne_bids/dig.py

@@ -1,27 +1,30 @@
 """Read/write BIDS compatible electrode/coords structures from MNE."""
 # Authors: Adam Li <[email protected]>
 #          Stefan Appelhoff <[email protected]>
+#          Alex Rockhill <[email protected]>
 #
 # License: BSD-3-Clause
 import json
 from collections import OrderedDict
 from pathlib import Path
 import re
 import warnings
+from importlib_resources import files
 
 import mne
 import numpy as np
 from mne.io.constants import FIFF
-from mne.transforms import _str_to_frame
-from mne.utils import logger, warn
+from mne.utils import logger, warn, _validate_type, _check_option
 from mne.io.pick import _picks_to_idx
 
-from mne_bids.config import (ALLOWED_SPACES,
-                             BIDS_COORDINATE_UNITS,
+from mne_bids.config import (ALLOWED_SPACES, BIDS_COORDINATE_UNITS,
                              MNE_TO_BIDS_FRAMES, BIDS_TO_MNE_FRAMES,
-                             MNE_FRAME_TO_STR, BIDS_COORD_FRAME_DESCRIPTIONS)
+                             MNE_FRAME_TO_STR, MNE_STR_TO_FRAME,
+                             BIDS_COORD_FRAME_DESCRIPTIONS,
+                             BIDS_STANDARD_TEMPLATE_COORDINATE_SYSTEMS)
 from mne_bids.tsv_handler import _from_tsv
-from mne_bids.utils import _scale_coord_to_meters, _write_json, _write_tsv
+from mne_bids.utils import (_scale_coord_to_meters, _write_json, _write_tsv,
+                            verbose)
 from mne_bids.path import BIDSPath
 
 
@@ -396,9 +399,9 @@ def _write_dig_bids(bids_path, raw, montage=None, acpc_aligned=False,
                                     message='.*nasion not found', module='mne')
             raw.set_montage(montage)
         for ch in raw.info['chs']:
-            ch['coord_frame'] = _str_to_frame[montage_coord_frame]
+            ch['coord_frame'] = MNE_STR_TO_FRAME[montage_coord_frame]
         for d in raw.info['dig']:
-            d['coord_frame'] = _str_to_frame[montage_coord_frame]
+            d['coord_frame'] = MNE_STR_TO_FRAME[montage_coord_frame]
         with raw.info._unlock():  # add back fiducials
             raw.info['dig'] = fids + raw.info['dig']
 
@@ -543,6 +546,104 @@ def _read_dig_bids(electrodes_fpath, coordsystem_fpath,
     # put back in unknown for unknown coordinate frame
     if coord_frame == 'unknown':
         for ch in raw.info['chs']:
-            ch['coord_frame'] = _str_to_frame['unknown']
+            ch['coord_frame'] = MNE_STR_TO_FRAME['unknown']
         for d in raw.info['dig']:
-            d['coord_frame'] = _str_to_frame['unknown']
+            d['coord_frame'] = MNE_STR_TO_FRAME['unknown']
+
+
+@verbose
+def template_to_head(raw, space, coord_frame='auto', unit='auto',
+                     verbose=None):
+    """Transform a BIDS standard template montage to the head coordinate frame.
+
+    Parameters
+    ----------
+    raw : mne.io.Raw
+        The MNE-Python Raw object with a montage in one of the BIDS standard
+        template coordinate frames (modified in place).
+    space : str
+        The name of the BIDS standard template. See
+        https://bids-specification.readthedocs.io/en/stable/99-appendices/08-coordinate-systems.html#standard-template-identifiers
+        for a list of acceptable spaces.
+    coord_frame : str
+        BIDS template coordinate systems do not specify a coordinate frame
+        so this must be determined by inspecting the documentation for the
+        dataset or the electrodes.tsv file. Must be ``mri``, ``mri_voxel``
+        or ``ras``. If ``auto``, the montage will be assumed to be in
+        ``mri_voxel`` if the coordinates are strictly positive and
+        ``scanner RAS`` otherwise.
+
+        .. warning::
+
+            ``scanner RAS`` and ``surface RAS`` coordinates frames are similar
+            so be very careful not to assume a BIDS dataset's coordinates are
+            in one when they are actually in the other. The only way to tell
+            for template coordinate systems, currently, is if it is specified
+            in the dataset documentation.
+
+    unit : str
+        The unit that was used in the coordinate system specification; either
+        ``m`` or ``mm``. If ``auto``, ``m`` will be inferred if the montage
+        spans less than -1 to 1 in units (meters) and ``mm`` otherwise. If the
+        ``coord_frame`` is 'mri_voxel', ``unit`` will be ignored.
+    %(verbose)s
+
+    Returns
+    -------
+    trans : mne.transforms.Transform
+        The data transformation matrix from ``head`` to ``mri`` coordinates.
+
+    """
+    _validate_type(raw, mne.io.BaseRaw)
+    _check_option('space', space, BIDS_STANDARD_TEMPLATE_COORDINATE_SYSTEMS)
+    _check_option('coord_frame', coord_frame,
+                  ('auto', 'mri', 'mri_voxel', 'ras'))
+    _check_option('unit', unit, ('auto', 'm', 'mm'))
+    montage = raw.get_montage()
+    if montage is None:
+        raise RuntimeError('No montage found in the `raw` object')
+    montage.remove_fiducials()  # we will add fiducials so remove any
+    pos = montage.get_positions()
+    if pos['coord_frame'] not in ('mni_tal', 'unknown'):
+        raise RuntimeError(
+            "Montage coordinate frame '{}' not expected for a template "
+            "montage, should be 'unknown' or 'mni_tal'".format(
+                pos['coord_frame']))
+    locs = np.array(list(pos['ch_pos'].values()))
+    locs = locs[~np.any(np.isnan(locs), axis=1)]  # only channels with loc
+    if locs.size == 0:
+        raise RuntimeError('No channel locations found in the montage')
+    if unit == 'auto':
+        unit = 'm' if abs(locs - locs.mean(axis=0)).max() < 1 else 'mm'
+    if coord_frame == 'auto':
+        coord_frame = 'mri_voxel' if locs.min() >= 0 else 'ras'
+    # transform montage to head
+    data_dir = files('mne_bids.data')
+    # set to the right coordinate frame as specified by the user
+    for d in montage.dig:  # ensure same coordinate frame
+        d['coord_frame'] = MNE_STR_TO_FRAME[coord_frame]
+    # do the transforms, first ras -> vox if needed
+    if montage.get_positions()['coord_frame'] == 'ras':
+        ras_vox_trans = mne.read_trans(
+            data_dir.joinpath(f'space-{space}_ras-vox_trans.fif'))
+        if unit == 'm':  # must be in mm here
+            for d in montage.dig:
+                d['r'] *= 1000
+        montage.apply_trans(ras_vox_trans)
+    if montage.get_positions()['coord_frame'] == 'mri_voxel':
+        vox_mri_trans = mne.read_trans(
+            data_dir.joinpath(f'space-{space}_vox-mri_trans.fif'))
+        montage.apply_trans(vox_mri_trans)
+    assert montage.get_positions()['coord_frame'] == 'mri'
+    if not (unit == 'm' and coord_frame == 'mri'):  # if so, already in m
+        for d in montage.dig:
+            d['r'] /= 1000  # mm -> m
+    # now add fiducials (in mri coordinates)
+    fids = mne.io.read_fiducials(
+        data_dir.joinpath(f'space-{space}_fiducials.fif'))[0]
+    montage.dig = fids + montage.dig  # add fiducials
+    for fid in fids:  # ensure also in mri
+        fid['coord_frame'] = MNE_STR_TO_FRAME['mri']
+    raw.set_montage(montage)  # transform to head
+    # finally return montage
+    return mne.read_trans(data_dir.joinpath(f'space-{space}_trans.fif'))

mne_bids/tests/test_dig.py

@@ -10,13 +10,16 @@
 import os.path as op
 import numpy as np
 import pytest
+from numpy.testing import assert_almost_equal
+import warnings
 
 import mne
+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
-from mne_bids.config import (BIDS_STANDARD_TEMPLATE_COORDINATE_FRAMES,
-                             BIDS_TO_MNE_FRAMES)
+from mne_bids.dig import _write_dig_bids, _read_dig_bids, template_to_head
+from mne_bids.config import (BIDS_STANDARD_TEMPLATE_COORDINATE_SYSTEMS,
+                             BIDS_TO_MNE_FRAMES, MNE_STR_TO_FRAME)
 
 base_path = op.join(op.dirname(mne.__file__), 'io')
 subject_id = '01'
@@ -101,7 +104,7 @@ def test_dig_pixels(tmp_path):
                        bids_path.datatype, raw_test)
     mnt2 = raw_test.get_montage()
     assert mnt2.get_positions()['coord_frame'] == 'unknown'
-    np.testing.assert_array_almost_equal(
+    assert_almost_equal(
         np.array(list(mnt.get_positions()['ch_pos'].values())),
         np.array(list(mnt2.get_positions()['ch_pos'].values()))
     )
@@ -118,7 +121,7 @@ def test_dig_template(tmp_path):
 
     for datatype in ('eeg', 'ieeg'):
         bids_path = _bids_path.copy().update(root=bids_root, datatype=datatype)
-        for coord_frame in BIDS_STANDARD_TEMPLATE_COORDINATE_FRAMES:
+        for coord_frame in BIDS_STANDARD_TEMPLATE_COORDINATE_SYSTEMS:
             bids_path.update(space=coord_frame)
             mnt = montage.copy()
             pos = mnt.get_positions()
@@ -156,11 +159,128 @@ def test_dig_template(tmp_path):
 
     # test MEG
     raw_test = raw.copy()
-    for coord_frame in BIDS_STANDARD_TEMPLATE_COORDINATE_FRAMES:
+    for coord_frame in BIDS_STANDARD_TEMPLATE_COORDINATE_SYSTEMS:
         bids_path = _bids_path.copy().update(root=bids_root, datatype='meg',
                                              space=coord_frame)
         write_raw_bids(raw_test, bids_path)
         raw_test2 = read_raw_bids(bids_path)
         for ch, ch2 in zip(raw.info['chs'], raw_test2.info['chs']):
             np.testing.assert_array_equal(ch['loc'], ch2['loc'])
             assert ch['coord_frame'] == ch2['coord_frame']
+
+
+def _set_montage_no_trans(raw, montage):
+    """Set the montage without transforming to 'head'."""
+    coord_frame = montage.get_positions()['coord_frame']
+    with warnings.catch_warnings():
+        warnings.filterwarnings(action='ignore', category=RuntimeWarning,
+                                message='.*nasion not found', module='mne')
+        raw.set_montage(montage, on_missing='ignore')
+    for ch in raw.info['chs']:
+        ch['coord_frame'] = MNE_STR_TO_FRAME[coord_frame]
+    for d in raw.info['dig']:
+        d['coord_frame'] = MNE_STR_TO_FRAME[coord_frame]
+
+
+def _test_montage_trans(raw, montage, pos_test, space='fsaverage',
+                        coord_frame='auto', unit='auto'):
+    """Test if a montage is transformed correctly."""
+    _set_montage_no_trans(raw, montage)
+    trans = template_to_head(
+        raw, space, coord_frame=coord_frame, unit=unit)
+    montage_test = raw.get_montage()
+    montage_test.apply_trans(trans)
+    assert_almost_equal(
+        pos_test,
+        np.array(list(montage_test.get_positions()['ch_pos'].values())))
+
+
+def test_template_to_head():
+    """Test transforming a template montage to head."""
+    # test no montage
+    raw_test = raw.copy()
+    raw_test.set_montage(None)
+    with pytest.raises(RuntimeError, match='No montage found'):
+        template_to_head(raw_test, 'fsaverage', coord_frame='auto')
+
+    # test no channels
+    montage_empty = mne.channels.make_dig_montage(hsp=[[0, 0, 0]])
+    _set_montage_no_trans(raw_test, montage_empty)
+    with pytest.raises(RuntimeError, match='No channel locations '
+                                           'found in the montage'):
+        template_to_head(raw_test, 'fsaverage', coord_frame='auto')
+
+    # test unexpected coordinate frame
+    raw_test = raw.copy()
+    with pytest.raises(RuntimeError, match='not expected for a template'):
+        template_to_head(raw_test, 'fsaverage', coord_frame='auto')
+
+    # test all coordinate frames
+    raw_test = raw.copy()
+    raw_test.set_montage(None)
+    raw_test.pick_types(eeg=True)
+    raw_test.drop_channels(raw_test.ch_names[3:])
+    montage = mne.channels.make_dig_montage(
+        ch_pos={raw_test.ch_names[0]: [0, 0, 0],
+                raw_test.ch_names[1]: [0, 0, 0.1],
+                raw_test.ch_names[2]: [0, 0, 0.2]},
+        coord_frame='unknown')
+    for space in BIDS_STANDARD_TEMPLATE_COORDINATE_SYSTEMS:
+        for cf in ('mri', 'mri_voxel', 'ras'):
+            _set_montage_no_trans(raw_test, montage)
+            trans = template_to_head(raw_test, space, cf)
+            assert trans['from'] == MNE_STR_TO_FRAME['head']
+            assert trans['to'] == MNE_STR_TO_FRAME['mri']
+            montage_test = raw_test.get_montage()
+            pos = montage_test.get_positions()
+            assert pos['coord_frame'] == 'head'
+            assert pos['nasion'] is not None
+            assert pos['lpa'] is not None
+            assert pos['rpa'] is not None
+
+    # test that we get the right transform
+    _set_montage_no_trans(raw_test, montage)
+    trans = template_to_head(raw_test, 'fsaverage', 'mri')
+    trans2 = mne.read_trans(op.join(
+        op.dirname(op.dirname(mne_bids.__file__)), 'mne_bids', 'data',
+        'space-fsaverage_trans.fif'))
+    assert_almost_equal(trans['trans'], trans2['trans'])
+
+    # test auto coordinate frame
+
+    # test auto voxels
+    montage_vox = mne.channels.make_dig_montage(
+        ch_pos={raw_test.ch_names[0]: [2, 0, 10],
+                raw_test.ch_names[1]: [0, 0, 5.5],
+                raw_test.ch_names[2]: [0, 1, 3]},
+        coord_frame='unknown')
+    pos_test = np.array([[0.126, -0.118, 0.128],
+                         [0.128, -0.1225, 0.128],
+                         [0.128, -0.125, 0.127]])
+    _test_montage_trans(raw_test, montage_vox, pos_test,
+                        coord_frame='auto', unit='mm')
+
+    # now negative values => scanner RAS
+    montage_ras = mne.channels.make_dig_montage(
+        ch_pos={raw_test.ch_names[0]: [-30.2, 20, -40],
+                raw_test.ch_names[1]: [10, 30, 53.5],
+                raw_test.ch_names[2]: [30, -21, 33]},
+        coord_frame='unknown')
+    pos_test = np.array([[-0.0302, 0.02, -0.04],
+                         [0.01, 0.03, 0.0535],
+                         [0.03, -0.021, 0.033]])
+    _set_montage_no_trans(raw_test, montage_ras)
+    _test_montage_trans(raw_test, montage_ras, pos_test,
+                        coord_frame='auto', unit='mm')
+
+    # test auto unit
+    montage_mm = montage_ras.copy()
+    _set_montage_no_trans(raw_test, montage_mm)
+    _test_montage_trans(raw_test, montage_mm, pos_test,
+                        coord_frame='ras', unit='auto')
+
+    montage_m = montage_ras.copy()
+    for d in montage_m.dig:
+        d['r'] = np.array(d['r']) / 1000
+    _test_montage_trans(raw_test, montage_m, pos_test,
+                        coord_frame='ras', unit='auto')

mne_bids/write.py

@@ -56,7 +56,7 @@
                              BIDS_VERSION, REFERENCES, _map_options, reader,
                              ALLOWED_INPUT_EXTENSIONS, CONVERT_FORMATS,
                              ANONYMIZED_JSON_KEY_WHITELIST,
-                             BIDS_STANDARD_TEMPLATE_COORDINATE_FRAMES)
+                             BIDS_STANDARD_TEMPLATE_COORDINATE_SYSTEMS)
 
 
 _FIFF_SPLIT_SIZE = '2GB'  # MNE-Python default; can be altered during debugging
@@ -1610,7 +1610,7 @@ def write_raw_bids(raw, bids_path, events_data=None, event_id=None,
             sensor_coord_system = orient
         elif orient == 'n/a':
             sensor_coord_system = bids_path.space
-        elif bids_path.space in BIDS_STANDARD_TEMPLATE_COORDINATE_FRAMES:
+        elif bids_path.space in BIDS_STANDARD_TEMPLATE_COORDINATE_SYSTEMS:
             sensor_coord_system = bids_path.space
         elif orient != bids_path.space:
             raise ValueError(f'BIDSPath.space {bids_path.space} conflicts '

setup.cfg

@@ -38,6 +38,7 @@ install_requires =
   scipy >= 1.2.0
   jinja2  # Can be removed as soon as we drop support for MNE <1.0
   setuptools
+  importlib-resources
 packages = find:
 include_package_data = True