[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,137 @@ 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']
+
+
+def _get_montage_from_info(info):
+    """Workaround to get montage from info."""
+    # XXX remove upon 0.11 release
+    if info['dig'] is None:
+        return None
+    # obtain coord_frame, and landmark coords
+    # (nasion, lpa, rpa, hsp, hpi) from DigPoints
+    montage_bunch = mne.io._digitization._get_data_as_dict_from_dig(
+        info['dig'])
+    coord_frame = mne.transforms._frame_to_str.get(montage_bunch.coord_frame)
+
+    # get the channel names and chs data structure
+    ch_names, chs = info['ch_names'], info['chs']
+    picks = mne.pick_types(info, meg=False, eeg=True, seeg=True,
+                           ecog=True, dbs=True, fnirs=True, exclude=[])
+
+    # channel positions from dig do not match ch_names one to one,
+    # so use loc[:3] instead
+    ch_pos = {ch_names[ii]: chs[ii]['loc'][:3] for ii in picks}
+
+    # create montage
+    montage = mne.channels.make_dig_montage(
+        ch_pos=ch_pos,
+        coord_frame=coord_frame,
+        nasion=montage_bunch.nasion,
+        lpa=montage_bunch.lpa,
+        rpa=montage_bunch.rpa,
+        hsp=montage_bunch.hsp,
+        hpi=montage_bunch.hpi,
+    )
+    return montage
+
+
+@verbose
+def template_to_head(info, space, coord_frame='auto', unit='auto',
+                     verbose=None):
+    """Transform a BIDS standard template montage to the head coordinate frame.
+
+    Parameters
+    ----------
+    %(info_not_none)s The info is 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(info, mne.io.Info)
+    _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'))
+    # XXX: change to after 0.11 release
+    # montage = info.get_montage()
+    montage = _get_montage_from_info(info)
+    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 / 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 / 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 / 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']
+    info.set_montage(montage)  # transform to head
+    # finally return montage
+    return mne.read_trans(data_dir / f'space-{space}_trans.fif')