[ENH, WIP] Add support for BIDS template standard coordinates

examples/convert_ieeg_to_bids.py

@@ -18,6 +18,9 @@
 
 5. Repeat the process for the ``fsaverage`` template coordinate frame.
 
+6. Repeat the process for one of the other standard template coordinate frames
+   allowed by BIDS.
+
 The iEEG data will be written by :func:`write_raw_bids` with
 the addition of extra metadata elements in the following files:
 
@@ -48,6 +51,7 @@
 # %%
 
 import os.path as op
+import numpy as np
 import shutil
 
 from nilearn.plotting import plot_anat
@@ -70,6 +74,7 @@
 # which is easily read in using numpy
 raw = mne.io.read_raw_fif(op.join(
     misc_path, 'seeg', 'sample_seeg_ieeg.fif'))
+raw.set_channel_types({ch: 'ecog' for ch in raw.ch_names})  # fake intracranial
 raw.info['line_freq'] = 60  # specify power line frequency as required by BIDS
 subjects_dir = op.join(misc_path, 'seeg')  # Freesurfer recon-all directory
 
@@ -195,7 +200,7 @@
 #
 # `acpc_aligned=True` affirms that our MRI is aligned to ACPC
 # if this is not true, convert to `fsaverage` (see below)!
-write_raw_bids(raw, bids_path, anonymize=dict(daysback=30000),
+write_raw_bids(raw, bids_path, anonymize=dict(daysback=40000),
                montage=montage, acpc_aligned=True, overwrite=True)
 
 # check our output
@@ -271,7 +276,14 @@
 # from "mri" aka surface RAS space.
 # ``fsaverage`` is very useful for group analysis as shown in
 # `Working with SEEG
-# <https://mne.tools/stable/auto_tutorials/misc/plot_seeg.html>`_.
+# <https://mne.tools/stable/auto_tutorials/misc/plot_seeg.html>`_. Note, this
+# is only a linear transform and so loses quite a bit of accuracy relative
+# to the needs of intracranial researchers so it is quite suboptimal. A better
+# option is to use a symmetric diffeomorphic transform to create a one-to-one
+# mapping of brain voxels from the individual's brain to the template. Even so,
+# it's better to provide the coordinates in the individual's brain space, as
+# was done above, so that the researcher who uses the coordinates has the
+# ability to tranform them to a template of their choice (see below).
 
 # ensure the output path doesn't contain any leftover files from previous
 # tests and example runs
@@ -282,6 +294,7 @@
 raw = mne.io.read_raw_fif(op.join(
     misc_path, 'seeg', 'sample_seeg_ieeg.fif'))
 raw.info['line_freq'] = 60  # specify power line frequency as required by BIDS
+raw.set_channel_types({ch: 'ecog' for ch in raw.ch_names})
 
 # get Talairach transform
 mri_mni_t = mne.read_talxfm('sample_seeg', subjects_dir)
@@ -293,7 +306,7 @@
 montage.apply_trans(mri_mni_t)
 
 # write to BIDS, this time with a template coordinate system
-write_raw_bids(raw, bids_path, anonymize=dict(daysback=30000),
+write_raw_bids(raw, bids_path, anonymize=dict(daysback=40000),
                montage=montage, overwrite=True)
 
 # read in the BIDS dataset
@@ -322,3 +335,87 @@
 # ``montage.add_estimated_fiducials(template, os.environ['FREESURFER_HOME'])``
 # where ``template`` maybe be ``cvs_avg35_inMNI152`` for instance
 raw.set_montage(montage)
+
+# %%
+# Step 6: Store coordinates in another template space accepted by BIDS
+# --------------------------------------------------------------------
+# As of thid writing, BIDS accepts channel coordinates in reference to the
+# the following template spaces: ``ICBM452AirSpace``, ``ICBM452Warp5Space``,
+# ``IXI549Space``, ``fsaverage``, ``fsaverageSym``, ``fsLR``, ``MNIColin27``,
+# ``MNI152Lin``, ``MNI152NLin2009[a-c][Sym|Asym]``, ``MNI152NLin6Sym``,
+# ``MNI152NLin6ASym``, ``MNI305``, ``NIHPD``, ``OASIS30AntsOASISAnts``,
+# ``OASIS30Atropos``, ``Talairach`` and ``UNCInfant``. As discussed above,
+# it is recommended to share the coordinates in the individual subject's
+# anatomical reference frame so that researchers who use the data can
+# transform the coordinates to any of these templates that they choose. If
+# BIDS-formatted data is shared in one of these template coordinate frames,
+# the data can still be analyzed in MNE-Python. However, MNE-Python only
+# recognizes a few coordinate frames (so that coordinate frames that are
+# not regularly used by the MNE community don't misleadingly appear to be
+# being fully maintained when they are not) so you'll need a bit more
+# know-how, which we will go over below.
+
+# ensure the output path doesn't contain any leftover files from previous
+# tests and example runs
+if op.exists(bids_root):
+    shutil.rmtree(bids_root)
+
+
+# first we'll write our data as if it were MNI152NLin2009bAsym coordinates
+# (you would need to transform your coordinates to this template first)
+bids_path.update(datatype='ieeg', space='MNI152NLin2009bAsym')
+
+# load our raw data again
+raw = mne.io.read_raw_fif(op.join(
+    misc_path, 'seeg', 'sample_seeg_ieeg.fif'))
+raw.set_channel_types({ch: 'seeg' for ch in raw.ch_names})
+raw.info['line_freq'] = 60  # specify power line frequency as required by BIDS
+
+# get the montage as stored in raw
+# MNE stores coordinates in raw objects in "head" coordinates for consistency
+montage = raw.get_montage()
+
+# define a transform to MNI152NLin2009bAsym (fake it)
+# MNE-Python doesn't recognize MNI152NLin2009bAsym, so we have to use
+# the unknown coordinate frame
+head_template_t = np.array([[1.0, 0.1, 0.2, 10.1],
+                            [-0.1, 1.0, 0.1, -20.3],
+                            [0.2, -0.1, 1.0, -30.7],
+                            [0.0, 0.0, 0.0, 1.0]])
+head_template_trans = mne.transforms.Transform(
+    fro='head', to='unknown', trans=head_template_t)
+montage.apply_trans(head_template_trans)
+
+# get fiducials in the template coordinate frame for later
+pos = montage.get_positions()
+nasion, lpa, rpa = (pos[fid] for fid in ('nasion', 'lpa', 'rpa'))
+
+# write to BIDS, with the montage in fsaverage coordinates
+write_raw_bids(raw, bids_path, anonymize=dict(daysback=40000),
+               montage=montage, overwrite=True)
+
+# read back in the raw data
+raw = read_raw_bids(bids_path=bids_path)
+
+# get the montage
+montage = raw.get_montage()
+print('Montage set to: ' + montage.get_positions()['coord_frame'])
+
+# now we can use the fiducials that we stored before to store the
+# coordinates in the "head" coordinate frame that is prefered
+# you would have to find the fiducials for the template its
+# documentation or using Freeview on the template MRI or with
+# the MNE coregistration GUI if you didn't have them already
+pos = montage.get_positions()
+montage = mne.channels.make_dig_montage(
+    ch_pos=pos['ch_pos'], nasion=nasion, lpa=lpa, rpa=rpa,
+    coord_frame=pos['coord_frame'])
+
+# now check that we can recover our transform
+head_template_trans2 = mne.channels.compute_native_head_t(montage)
+print(f'Original: {head_template_trans}\nRecovered: {head_template_trans2}')
+
+# finally we have a montage in head coordinates that can be
+# transformed back the template coordinate frame stored as
+# unknown in MNE
+raw.set_montage(montage)  # transforms to head

mne_bids/config.py

@@ -153,7 +153,7 @@
 # Annotations to never remove during reading or writing
 ANNOTATIONS_TO_KEEP = ('BAD_ACQ_SKIP',)
 
-coordsys_standard_template = [
+BIDS_STANDARD_TEMPLATE_COORDINATE_FRAMES = [
     'ICBM452AirSpace',
     'ICBM452Warp5Space',
     'IXI549Space',
@@ -198,18 +198,23 @@
 coordsys_meg = ['CTF', 'ElektaNeuromag', '4DBti', 'KitYokogawa', 'ChietiItab']
 coordsys_eeg = ['CapTrak']
 coordsys_ieeg = ['Pixels', 'ACPC']
-coordsys_wildcard = ['Other']
-coordsys_shared = (coordsys_standard_template +
-                   coordsys_standard_template_deprecated +
-                   coordsys_wildcard)
+coordsys_shared = (BIDS_STANDARD_TEMPLATE_COORDINATE_FRAMES +
+                   coordsys_standard_template_deprecated)
 
 ALLOWED_SPACES = dict()
-ALLOWED_SPACES['meg'] = coordsys_shared + coordsys_meg + coordsys_eeg
-ALLOWED_SPACES['eeg'] = coordsys_shared + coordsys_meg + coordsys_eeg
+ALLOWED_SPACES['meg'] = ALLOWED_SPACES['eeg'] = \
+    coordsys_shared + coordsys_meg + coordsys_eeg
 ALLOWED_SPACES['ieeg'] = coordsys_shared + coordsys_ieeg
 ALLOWED_SPACES['anat'] = None
 ALLOWED_SPACES['beh'] = None
 
+# write only
+ALLOWED_SPACES_WRITE = dict()
+ALLOWED_SPACES_WRITE['meg'] = ALLOWED_SPACES_WRITE['eeg'] = \
+    BIDS_STANDARD_TEMPLATE_COORDINATE_FRAMES + coordsys_meg + coordsys_eeg
+ALLOWED_SPACES_WRITE['ieeg'] = \
+    BIDS_STANDARD_TEMPLATE_COORDINATE_FRAMES + coordsys_ieeg
+
 # See: https://bids-specification.readthedocs.io/en/latest/99-appendices/04-entity-table.html#encephalography-eeg-ieeg-and-meg  # noqa
 ENTITY_VALUE_TYPE = {
     'subject': 'label',
@@ -259,9 +264,6 @@
     'head': 'CapTrak',
     'mni_tal': 'fsaverage',
     # 'fs_tal': 'fsaverage',  # XXX: not used
-    'unknown': 'Other',
-    'ras': 'Other',
-    'mri': 'Other'
 }
 
 # these coordinate frames in mne-python are related to scalp/meg
@@ -282,11 +284,18 @@
 MNE_FRAME_TO_STR = {val: key for key, val in MNE_STR_TO_FRAME.items()}
 
 # see BIDS specification for description we copied over from each
-BIDS_COORD_FRAME_DESCRIPTIONS = {
+BIDS_COORD_FRAME_DESCRIPTIONS = dict(**{
     'acpc': 'The origin of the coordinate system is at the Anterior '
             'Commissure and the negative y-axis is passing through the '
             'Posterior Commissure. The positive z-axis is passing through '
             'a mid-hemispheric point in the superior direction.',
+    'pixels': 'If electrodes are localized in 2D space (only x and y are '
+              'specified and z is n/a), then the positions in this file '
+              'must correspond to the locations expressed in pixels on '
+              'the photo/drawing/rendering of the electrodes on the brain. '
+              'In this case, coordinates must be (row,column) pairs, with '
+              '(0,0) corresponding to the upper left pixel and (N,0) '
+              'corresponding to the lower left pixel.',
     'ctf': 'ALS orientation and the origin between the ears',
     'elektaneuromag': 'RAS orientation and the origin between the ears',
     '4dbti': 'ALS orientation and the origin between the ears',
@@ -304,7 +313,59 @@
     'fsaverage': 'Defined by FreeSurfer, the MRI (surface RAS) origin is '
                  'at the center of a 256×256×256 mm^3 anisotropic volume '
                  '(may not be in the center of the head).',
-}
+    'icbm452airspace': 'Reference space defined by the "average of 452 '
+                       'T1-weighted MRIs of normal young adult brains" '
+                       'with "linear transforms of the subjects into the '
+                       'atlas space using a 12-parameter affine '
+                       'transformation"',
+    'icbm452warp5space': 'Reference space defined by the "average of 452 '
+                         'T1-weighted MRIs of normal young adult brains" '
+                         '"based on a 5th order polynomial transformation '
+                         'into the atlas space"',
+    'ixi549space': 'Reference space defined by the average of the "549 (...) '
+                   'subjects from the IXI dataset" linearly transformed to '
+                   'ICBM MNI 452.',
+    'fsaveragesym': 'The fsaverage is a dual template providing both '
+                    'volumetric and surface coordinates references. The '
+                    'volumetric template corresponds to a FreeSurfer variant '
+                    'of MNI305 space. The fsaverageSym atlas also defines a '
+                    'symmetric surface reference system (formerly described '
+                    'as fsaveragesym).',
+    'fslr': 'The fsLR is a dual template providing both volumetric and '
+            'surface coordinates references. The volumetric template '
+            'corresponds to MNI152NLin6Asym. Surface templates are given '
+            'at several sampling densities: 164k (used by HCP pipelines '
+            'for 3T and 7T anatomical analysis), 59k (used by HCP pipelines '
+            'for 7T MRI bold and DWI analysis), 32k (used by HCP pipelines '
+            'for 3T MRI bold and DWI analysis), or 4k (used by HCP '
+            'pipelines for MEG analysis) fsaverage_LR surface '
+            'reconstructed from the T1w image.',
+    'mnicolin27': 'Average of 27 T1 scans of a single subject.',
+    'mni152lin': 'Also known as ICBM (version with linear coregistration).',
+    'mni152nlin6sym': 'Also known as symmetric ICBM 6th generation '
+                      '(non-linear coregistration).',
+    'mni152nlin6asym': 'A variation of MNI152NLin6Sym built by A. Janke that '
+                       'is released as the MNI template of FSL. Volumetric '
+                       'templates included with HCP-Pipelines correspond to '
+                       'this template too.',
+    'mni305': 'Also known as avg305.',
+    'nihpd': 'Pediatric templates generated from the NIHPD sample. Available '
+             'for different age groups (4.5–18.5 y.o., 4.5–8.5 y.o., '
+             '7–11 y.o., 7.5–13.5 y.o., 10–14 y.o., 13–18.5 y.o. This '
+             'template also comes in either -symmetric or -asymmetric flavor.',
+    'oasis30antsoasisants':
+        'See https://figshare.com/articles/ANTs_ANTsR_Brain_Templates/915436',
+    'oasis30atropos':
+        'See https://mindboggle.info/data.html',
+    'talairach': 'Piecewise linear scaling of the brain is implemented as '
+                 'described in TT88.',
+    'uncinfant': 'Infant Brain Atlases from Neonates to 1- and 2-year-olds.'
+},
+    **{f'mni152nlin2009{letter}{sym}': 'Also known as ICBM '
+       '(non-linear coregistration with 40 iterations, '
+       'released in 2009). It comes in either three different flavours '
+       'each in symmetric or asymmetric version.'
+       for letter in ('a', 'b', 'c') for sym in ('Sym', 'Asym')})
 
 REFERENCES = {'mne-bids':
               'Appelhoff, S., Sanderson, M., Brooks, T., Vliet, M., '