[ENH, MRG] Allow template coordinate frame writing

doc/whats_new.rst

@@ -56,6 +56,8 @@ Enhancements
 
 - 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`)
+
 API and behavior changes
 ^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -65,6 +67,8 @@ API and behavior changes
 
 - Passing ``fs_subject=None`` to :func:`get_head_mri_trans` has been deprecated. Please pass the FreeSurfer subject name explicitly, by Richard Höchenberger`_ (:gh:`977`)
 
+- Corrupted or missing fiducials in 'head' coordinates now raise an error instead of warning in :func:`mne_bids.write_raw_bids` by `Alex Rockhill`_ (:gh:`980`)
+
 Requirements
 ^^^^^^^^^^^^
 

examples/convert_ieeg_to_bids.py

@@ -20,6 +20,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:
 
@@ -50,6 +53,7 @@
 # %%
 
 import os.path as op
+import numpy as np
 import shutil
 
 from nilearn.plotting import plot_anat
@@ -331,3 +335,133 @@
 # ``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.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)
+
+# write to BIDS, with the montage in fsaverage coordinates
+write_raw_bids(raw, bids_path, anonymize=dict(daysback=40000),
+               montage=montage, overwrite=True)
+
+# %%
+# Now, let's see how we would work with the data in MNE. As shown below, the
+# montage has the same coordinates as when it was written, but the coordinate
+# frame is unknown.
+
+# read back in the raw data
+raw = read_raw_bids(bids_path=bids_path)
+
+# get the montage
+montage2 = raw.get_montage()
+print('Montage set to: ' + montage2.get_positions()['coord_frame'])
+
+# check that we can recover the coordinates
+print('Recovered coordinate: {recovered}\n'
+      'Saved coordinate:     {saved}'.format(
+          recovered=montage2.dig[0]['r'],
+          saved=montage.dig[0]['r']))
+
+# %%
+# To work with this data in the template coordinate frame, all the same steps
+# can be followed in :ref:`tut-working-with-seeg` and
+# :ref:`tut-working-with-ecog` after the montage is transformed to surface
+# RAS coordinates for the template MRI (if it isn't there already).
+#
+# First we'll need the ``subjects_dir`` where the recon-all for the template
+# brain is stored.
+#
+# .. code-block:: python
+#
+#    subjects_dir = op.join(misc_path, 'subjects')  # for example
+#    # add some plotting keyword arguments
+#    brain_kwargs = dict(cortex='low_contrast', alpha=0.2, background='white')
+#
+# If the montage is already in surface RAS for the template MRI, we can use:
+#
+# .. code-block:: python
+#
+#      # identity transform since 'unknown' is already 'mri' == surface RAS
+#      # for the template brain MRI
+#      trans = mne.transforms.Transform(
+#          fro='unknown', to='mri', trans=np.eye(4))
+#      brain = mne.viz.Brain('sample_seeg', subjects_dir=subjects_dir,
+#          **brain_kwargs)
+#      brain.add_sensors(raw.info, trans=trans)
+#
+# If the montage was in voxel coordinates, we'll first have to transform
+# to surface RAS:
+#
+# .. code-block:: python
+#
+#      import nibabel as nib
+#      template_T1 = nib.load(op.join(subjects_dir, 'MNI152NLin2009bAsym',
+#                                     'mri', 'T1.mgz'))
+#      trans = mne.transforms.Transform(  # use vox to surface RAS transform
+#          fro='unknown', to='mri', trans=template_T1.header.get_vox2ras_tkr())
+#      brain = mne.viz.Brain(
+#          'sample_seeg', subjects_dir=subjects_dir, **brain_kwargs)
+#      brain.add_sensors(raw.info, trans=trans)
+#
+#
+# Finally, if the montage was in scanner RAS coordinates, we'll have to
+# transform it back to voxels first before going to surface RAS:
+#
+# .. code-block:: python
+#
+#      import nibabel as nib
+#      template_T1 = nib.load(op.join(subjects_dir, 'MNI152NLin2009bAsym',
+#                                     'mri', 'T1.mgz'))
+#      ras_vox_t = template_T1.header.get_ras2vox()
+#      vox_mri_t = template_T1.header.get_vox2ras_tkr()
+#      ras_mri_t = np.dot(ras_vox_t, vox_mri_t)  # ras->vox with vox->mri
+#      trans = mne.transforms.Transform(
+#          fro='unknown', to='mri', trans=ras_mri_t)
+#      brain = mne.viz.Brain(
+#          'sample_seeg', subjects_dir=subjects_dir, **brain_kwargs)
+#      brain.add_sensors(raw.info, trans=trans)

mne_bids/config.py

@@ -163,7 +163,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',
@@ -216,7 +216,7 @@
 # accepted coordinate SI units
 BIDS_COORDINATE_UNITS = ['m', 'cm', 'mm']
 coordsys_wildcard = ['Other']
-BIDS_SHARED_COORDINATE_FRAMES = (coordsys_standard_template +
+BIDS_SHARED_COORDINATE_FRAMES = (BIDS_STANDARD_TEMPLATE_COORDINATE_FRAMES +
                                  coordsys_standard_template_deprecated +
                                  coordsys_wildcard)
 
@@ -264,11 +264,7 @@
 MNE_TO_BIDS_FRAMES = {
     'ctf_head': 'CTF',
     'head': 'CapTrak',
-    'mni_tal': 'fsaverage',
-    # 'fs_tal': 'fsaverage',  # XXX: not used
-    'unknown': 'Other',
-    'ras': 'Other',
-    'mri': 'Other'
+    'mni_tal': 'fsaverage'
 }
 
 # these coordinate frames in mne-python are related to scalp/meg
@@ -289,7 +285,7 @@
 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 '
@@ -318,7 +314,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., '