diff --git a/examples/convert_ieeg_to_bids.py b/examples/convert_ieeg_to_bids.py index a52ad3623..4571c7446 100644 --- a/examples/convert_ieeg_to_bids.py +++ b/examples/convert_ieeg_to_bids.py @@ -18,10 +18,7 @@ 4. Cite MNE-BIDS. -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. +5. Repeat the process for the ``fsaverage`` template coordinate space. The iEEG data will be written by :func:`write_raw_bids` with the addition of extra metadata elements in the following files: @@ -56,6 +53,7 @@ import numpy as np import shutil +import nibabel as nib from nilearn.plotting import plot_anat import mne @@ -117,13 +115,12 @@ # - `background on FreeSurfer`_ # - `MNE-Python coordinate frames`_ # -# Currently, MNE-Python supports the ``mni_tal`` and ``mri`` coordinate frames, +# MNE-Python supports using ``mni_tal`` and ``mri`` coordinate frames, # corresponding to the ``fsaverage`` and ``ACPC`` (for an ACPC-aligned T1) BIDS # coordinate systems respectively. All other coordinate coordinate frames in -# MNE-Python if written with :func:`mne_bids.write_raw_bids` are written with -# coordinate system ``'Other'``. Note, then we suggest using -# :func:`mne_bids.update_sidecar_json` to update the sidecar -# ``*_coordsystem.json`` file to add additional information. +# MNE-Python, if written with :func:`mne_bids.write_raw_bids`, must have +# an :attr:`mne_bids.BIDSPath.space` specified, and will be read in with +# the montage channel locations set to the coordinate frame 'unknown'. # # Step 2: Formatting as BIDS # -------------------------- @@ -209,7 +206,7 @@ # %% # MNE-BIDS has created a suitable directory structure for us, and among other -# meta data files, it started an ``events.tsv``` and ``channels.tsv`` file, +# meta data files, it started an ``events.tsv`` and ``channels.tsv`` file, # and created an initial ``dataset_description.json`` file on top! # # Now it's time to manually check the BIDS directory and the meta files to add @@ -219,6 +216,7 @@ search_folder_for_text('n/a', bids_root) +# %% # Remember that there is a convenient JavaScript tool to validate all your BIDS # directories called the "BIDS-validator", available as a web version and a # command line tool: @@ -235,13 +233,7 @@ # :func:`read_raw_bids` to read in the data. # read in the BIDS dataset to plot the coordinates -raw = read_raw_bids(bids_path=bids_path) - -# compare with standard -print('Recovered coordinate: {recovered}\n' - 'Saved coordinate: {saved}'.format( - recovered=raw.info['chs'][0]['loc'][:3], - saved=montage.dig[0]['r'])) +raw2 = read_raw_bids(bids_path=bids_path) # %% # Now we have to go back to "head" coordinates with the head->mri transform. @@ -249,11 +241,31 @@ # .. note:: If you were downloading this from ``OpenNeuro``, you would # have to run the Freesurfer ``recon-all`` to get the transforms. -montage = raw.get_montage() +montage2 = raw2.get_montage() + # this uses Freesurfer recon-all subject directory -montage.add_estimated_fiducials('sample_seeg', subjects_dir=subjects_dir) +montage2.add_estimated_fiducials('sample_seeg', subjects_dir=subjects_dir) + +# get head->mri trans, invert from mri->head +trans2 = mne.transforms.invert_transform( + mne.channels.compute_native_head_t(montage2)) + # now the montage is properly in "head" and ready for analysis in MNE -raw.set_montage(montage) +raw2.set_montage(montage2) + +# get the monage, apply the trans and make sure it's the same +# note: the head coordinates may differ because they are defined by +# 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 +montage2.apply_trans(trans2) + +# compare with standard +print('Recovered coordinate: {recovered}\n' + 'Saved coordinate: {saved}'.format( + recovered=montage2.get_positions()['ch_pos']['LENT 1'], + saved=montage.get_positions()['ch_pos']['LENT 1'])) # %% # Step 4: Cite mne-bids @@ -310,13 +322,7 @@ montage=montage, overwrite=True) # read in the BIDS dataset -raw = read_raw_bids(bids_path=bids_path) - -# check that we can recover the coordinates -print('Recovered coordinate: {recovered}\n' - 'Saved coordinate: {saved}'.format( - recovered=raw.info['chs'][0]['loc'][:3], - saved=montage.dig[0]['r'])) +raw2 = read_raw_bids(bids_path=bids_path) # %% # Now we should go back to "head" coordinates. We do this with ``fsaverage`` @@ -325,152 +331,287 @@ # ``subjects_dir`` with ``fsaverage`` in it, which is accessible using # :func:`mne.datasets.fetch_fsaverage`. -montage = raw.get_montage() +montage2 = raw2.get_montage() # add fiducials for "mni_tal" (which is the coordinate frame fsaverage is in) # so that it can properly be set to "head" -montage.add_mni_fiducials(subjects_dir=subjects_dir) +montage2.add_mni_fiducials(subjects_dir=subjects_dir) -# Many other templates are included in the Freesurfer installation, -# so, for those, the fiducials can be estimated with -# ``montage.add_estimated_fiducials(template, os.environ['FREESURFER_HOME'])`` -# where ``template`` maybe be ``cvs_avg35_inMNI152`` for instance -raw.set_montage(montage) +# get the new head->mri (in this case mri == mni because fsavearge is in MNI) +mni_head_t = mne.channels.compute_native_head_t(montage2) -# %% -# Step 6: Store coordinates in another template space accepted by BIDS -# -------------------------------------------------------------------- -# As of this 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) +# set the montage transforming to the "head" coordinate frame +raw2.set_montage(montage2) +# check that we can recover the coordinates +print('Recovered coordinate head: {recovered}\n' + 'Saved coordinate head: {saved}'.format( + recovered=raw2.info['chs'][0]['loc'][:3], + saved=raw.info['chs'][0]['loc'][:3])) + +# check difference in trans +print('Recovered trans:\n{recovered}\n' + 'Original trans:\n{saved}'.format( + recovered=mni_head_t['trans'].round(3), + # combine head->mri with mri->mni to get head->mni + # and then invert to get mni->head + saved=np.linalg.inv(np.dot(trans['trans'], mri_mni_t['trans']) + ).round(3))) + +# ensure that the data in MNI coordinates is exactly the same +# (within numerical precision) +montage2 = raw2.get_montage() # get montage after transformed back to head +montage2.apply_trans(mne.transforms.invert_transform(mni_head_t)) +print('Recovered coordinate: {recovered}\n' + 'Saved coordinate: {saved}'.format( + recovered=montage2.get_positions()['ch_pos']['LENT 1'], + saved=montage.get_positions()['ch_pos']['LENT 1'])) -# 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') +# %% +# As you can see the coordinates stored in the ``raw`` object are slightly off. +# This is because the ``head`` coordinate frame is defined by the fiducials +# (nasion, left and right pre-auricular points), and, in the first case, +# the fiducials were found on the individual anatomy and then transformed +# to MNI space, whereas, in the second case, they were found directly on +# the template brain (this was done once for the template so that we could +# just load it from a file). This difference means that there are slightly +# different head->mri transforms. Once these transforms are applied, however, +# the positions are the same in MNI coordinates which is what is important. +# +# As a final step, let's go over how to assign the fiducials for a template +# brain where they are not found for you. Many template coordinate systems +# are allowed by BIDS but are not used in MNE-Python. +# +# .. note:: +# As of this 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. + +pos = montage.get_positions() +# fiducial points are included in the sample data but could be found using +# Freeview (note the coordinates are in MNE "mri" coordinates which is the +# same as Freesurfers TkRegRAS but in meters not millimeters -- +# divide the TkRegRAS values by 1000) +# or fiducial points could also be found with the MNE coregistration GUI +nas = pos['nasion'] +lpa = pos['lpa'] +rpa = pos['rpa'] + +print('Fiducial points determined from the template head anatomy:\n' + f'nasion: {nas}\nlpa: {lpa}\nrpa: {rpa}') + +# read raw in again to start over +raw2 = read_raw_bids(bids_path=bids_path) +montage2 = raw2.get_montage() + +# note: for fsaverage, the montage will be in the coordinate frame +# 'mni_tal' because it is recognized by MNE but other templates will be +# in the 'unknown' coordinate frame because they are not recognized by MNE +pos2 = montage2.get_positions() + +# here we will set the coordinate frame to be 'mri' because our channel +# positions and fiducials are in the Freesurfer surface RAS coordinate +# frame of the template T1 MRI (in this case fsaverage) +montage2 = mne.channels.make_dig_montage( # add fiducials + ch_pos=pos2['ch_pos'], + nasion=nas, + lpa=lpa, + rpa=rpa, + coord_frame='mri') + +# get head->mri trans, invert from mri->head +trans2 = mne.transforms.invert_transform( + mne.channels.compute_native_head_t(montage2)) + +# set the montage to transform back to 'head' +raw2.set_montage(montage2) + +# check that the coordinates were recovered +montage2 = raw2.get_montage() # get montage after transformed back to head +montage2.apply_trans(trans2) +print('Recovered coordinate: {recovered}\n' + 'Saved coordinate: {saved}'.format( + recovered=montage2.get_positions()['ch_pos']['LENT 1'], + saved=montage.get_positions()['ch_pos']['LENT 1'])) -# load our raw data again -raw = mne.io.read_raw_fif(op.join( +# %% +# We showed how to add the fiducials from ``montage.add_mni_fiducials`` +# yourself, which allows us to use any coordinate frame, even if it's not MNI, +# but we still assumed that the coordinates were in Freesurfer surface RAS. +# Unfortunately, this is not guranteed to be the case because the BIDS +# template descriptions only specify the anatomical space (as defined by the +# template T1 MRI) not the coordinate frame of the space; the coordinates +# could be in terms of surface RAS, voxels of the template MRI or scanner RAS +# MRI coordinates. See :ref:`tut-source-alignment` for a tutorial explaining +# the different coordinate frames. If the coordinates are in voxels or scanner +# RAS, we'll have to find the fiducials in those same coordinates. + +# get the template T1 (must be mgz to have surface RAS transforms +# the Freesurfer command `mri_convert` can be used to convert +subjects_dir = op.join(mne.datasets.sample.data_path(), 'subjects') +template_T1 = nib.load(op.join(subjects_dir, 'fsaverage', 'mri', 'T1.mgz')) + +# get vox->mri transform +vox_mri_t = template_T1.header.get_vox2ras_tkr() + +# transform the channel data to voxels just to demonstrate how to transform it +# back (this is the case where the BIDS-formatted data in the template space is +# in voxel coordinates so it would already be transformed when read in) +raw = mne.io.read_raw_fif(op.join( # load our raw data again 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) +montage = raw.get_montage() # get the original montage +montage.apply_trans(trans) # head->mri +mri_vox_trans = mne.transforms.Transform( + fro='mri', + to='mri_voxel', + trans=np.linalg.inv(vox_mri_t) +) +scale_t = np.eye(4) +scale_t[:3, :3] *= 1000 # m->mm +montage.apply_trans(mne.transforms.Transform('mri', 'mri', scale_t)) +montage.apply_trans(mri_vox_trans) # transform our original montage + +# print transformed fiducials, these could be found in the voxel locator in +# Freesurfer's `freeview` based on the template MRI since we wouldn't have +# them found for us for a template other than fsaverage +pos = montage.get_positions() +nas = pos['nasion'] +lpa = pos['lpa'] +rpa = pos['rpa'] + +print('Fiducial points determined from the template head anatomy in voxels:\n' + f'nasion: {nas}\nlpa: {lpa}\nrpa: {rpa}') + +# read raw in again to start over +raw2 = read_raw_bids(bids_path=bids_path) # %% -# 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. +# Now, it's the same as if we just got the montage from the raw in voxels +# i.e. if we would do ``montage = raw.get_montage()`` and get the positions +# directly in voxels instead of doing the transform to voxels first +montage2 = mne.channels.make_dig_montage( # add fiducials + ch_pos=pos['ch_pos'], + nasion=nas, + lpa=lpa, + rpa=rpa, + coord_frame='mri_voxel') +vox_mri_trans = mne.transforms.Transform( + fro='mri_voxel', + to='mri', + trans=vox_mri_t +) +montage2.apply_trans(vox_mri_trans) +scale_t = np.eye(4) +scale_t[:3, :3] /= 1000 # mm->m +montage2.apply_trans(mne.transforms.Transform('mri', 'mri', scale_t)) + +# get head->mri trans, invert from mri->head +trans2 = mne.transforms.invert_transform( + mne.channels.compute_native_head_t(montage2)) + +# set the montage to transform back to 'head' +raw2.set_montage(montage2) + +# check that the coordinates were recovered +montage = raw.get_montage() # get the original montage +montage.apply_trans(trans) # head->mri +montage2 = raw2.get_montage() # get montage after transformed back to head +montage2.apply_trans(trans2) +print('Recovered coordinate: {recovered}\n' + 'Saved coordinate: {saved}'.format( + recovered=montage2.get_positions()['ch_pos']['LENT 1'], + saved=montage.get_positions()['ch_pos']['LENT 1'])) -# read back in the raw data -raw = read_raw_bids(bids_path=bids_path) +# %% +# Finally, the template could also be in scanner RAS: -# get the montage -montage2 = raw.get_montage() -print('Montage set to: ' + montage2.get_positions()['coord_frame']) +# transform the channel data to scanner RAS just to demonstrate how to +# transform it back (this is the case where the BIDS-formatted data in the +# template space is in scanner RAS coordinates so it would already be +# transformed when read in) +raw = mne.io.read_raw_fif(op.join( # load our raw data again + misc_path, 'seeg', 'sample_seeg_ieeg.fif')) +montage = raw.get_montage() # get the original montage +montage.apply_trans(trans) # head->mri +scale_t = np.eye(4) +scale_t[:3, :3] *= 1000 # m->mm +montage.apply_trans(mne.transforms.Transform('mri', 'mri', scale_t)) +montage.apply_trans(mri_vox_trans) # transform our original montage +vox_ras_t = template_T1.header.get_vox2ras() # no tkr for scanner RAS +vox_ras_trans = mne.transforms.Transform( + fro='mri_voxel', + to='ras', + trans=vox_ras_t +) +montage.apply_trans(vox_ras_trans) + +# print transformed fiducials, these could be found in the RAS locator in +# Freesurfer's `freeview` based on the template MRI +# note in this case, they are in mm and should not be transformed +pos = montage.get_positions() +nas = pos['nasion'] +lpa = pos['lpa'] +rpa = pos['rpa'] + +print('Fiducial points determined from the template head anatomy in ' + f'scanner RAS:\nnasion: {nas}\nlpa: {lpa}\nrpa: {rpa}') + +# read raw in again to start over +raw2 = read_raw_bids(bids_path=bids_path) -# check that we can recover the coordinates +# %% +# Now, it's the same as if we just got the montage from the raw in scanner RAS +# i.e. we would use ``montage = raw.get_montage()`` instead of having to do +# the transforms above +montage2 = mne.channels.make_dig_montage( # add fiducials + ch_pos=pos['ch_pos'], + nasion=nas, + lpa=lpa, + rpa=rpa, + coord_frame='ras') +ras_vox_t = template_T1.header.get_ras2vox() +vox_mri_t = template_T1.header.get_vox2ras_tkr() +ras_mri_trans = mne.transforms.Transform( + fro='ras', + to='mri', + trans=np.dot(ras_vox_t, vox_mri_t) # combine ras->vox and vox->mri to get +) # ras->mri +montage2.apply_trans(ras_mri_trans) +scale_t = np.eye(4) +scale_t[:3, :3] /= 1000 # mm->m +montage2.apply_trans(mne.transforms.Transform('mri', 'mri', scale_t)) + +# get head->mri trans, invert from mri->head +trans2 = mne.transforms.invert_transform( + mne.channels.compute_native_head_t(montage2)) + +# set the montage to transform back to 'head' +raw2.set_montage(montage2) + +# check that the coordinates were recovered exactly this time +montage = raw.get_montage() # get the original montage +montage.apply_trans(trans) # head->mri +montage2 = raw2.get_montage() # get montage after transformed back to head +montage2.apply_trans(trans2) print('Recovered coordinate: {recovered}\n' 'Saved coordinate: {saved}'.format( - recovered=montage2.dig[0]['r'], - saved=montage.dig[0]['r'])) + recovered=montage2.get_positions()['ch_pos']['LENT 1'], + saved=montage.get_positions()['ch_pos']['LENT 1'])) # %% -# 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) +# In summary, as we saw, these standard template spaces that are allowable by +# BIDS are quite complicated. We therefore only cover these cases because +# datasets are allowed to be in these coordinate systems, and we want to be +# able to analyze them with MNE-Python. The template coordinate spaces +# don't specify a coordinate frame, so it is better to save the raw data in +# the individual's ACPC space, allowing the person analyzing the data to +# transform the positions to whatever template they want. Thus, we recommend +# if at all possible, saving BIDS iEEG data in ACPC coordinate space +# corresponding to the individual subject's brain, not in a template +# coordinate frame. diff --git a/mne_bids/dig.py b/mne_bids/dig.py index 35ff12f56..7c7f21fa2 100644 --- a/mne_bids/dig.py +++ b/mne_bids/dig.py @@ -21,8 +21,7 @@ MNE_TO_BIDS_FRAMES, BIDS_TO_MNE_FRAMES, MNE_FRAME_TO_STR, BIDS_COORD_FRAME_DESCRIPTIONS) from mne_bids.tsv_handler import _from_tsv -from mne_bids.utils import (_extract_landmarks, _scale_coord_to_meters, - _write_json, _write_tsv) +from mne_bids.utils import _scale_coord_to_meters, _write_json, _write_tsv from mne_bids.path import BIDSPath @@ -281,24 +280,24 @@ def _write_coordsystem_json(*, raw, unit, hpi_coord_system, Defaults to False. """ - dig = raw.info['dig'] - if dig is None: - dig = [] - - coord_frame = set([dig[ii]['coord_frame'] for ii in range(len(dig))]) - if len(coord_frame) > 1: # noqa E501 - raise ValueError('All HPI, electrodes, and fiducials must be in the ' - 'same coordinate frame. Found: "{}"' - .format(coord_frame)) + if raw.get_montage() is None: + dig = list() + coords = dict() + else: + montage = raw.get_montage() + pos = montage.get_positions() + dig = list() if montage.dig is None else montage.dig + coords = dict( + NAS=list() if pos['nasion'] is None else pos['nasion'].tolist(), + LPA=list() if pos['lpa'] is None else pos['lpa'].tolist(), + RPA=list() if pos['rpa'] is None else pos['rpa'].tolist()) # get the coordinate frame description sensor_coord_system_descr = (BIDS_COORD_FRAME_DESCRIPTIONS .get(sensor_coord_system.lower(), "n/a")) - coords = _extract_landmarks(dig) # create the coordinate json data structure based on 'datatype' if datatype == 'meg': - landmarks = dict(coords) hpi = {d['ident']: d for d in dig if d['kind'] == FIFF.FIFFV_POINT_HPI} if hpi: for ident in hpi.keys(): @@ -311,7 +310,7 @@ def _write_coordsystem_json(*, raw, unit, hpi_coord_system, 'HeadCoilCoordinates': coords, 'HeadCoilCoordinateSystem': hpi_coord_system, 'HeadCoilCoordinateUnits': unit, # XXX validate this - 'AnatomicalLandmarkCoordinates': landmarks, + 'AnatomicalLandmarkCoordinates': coords, 'AnatomicalLandmarkCoordinateSystem': sensor_coord_system, 'AnatomicalLandmarkCoordinateUnits': unit } @@ -386,18 +385,22 @@ def _write_dig_bids(bids_path, raw, montage=None, acpc_aligned=False, if montage is None: montage = raw.get_montage() - else: - # transform back to native montage coordinates if given + else: # assign montage to raw but supress any coordinate transforms + montage = montage.copy() # don't modify original + montage_coord_frame = montage.get_positions()['coord_frame'] + fids = [d for d in montage.dig # save to add back + if d['kind'] == FIFF.FIFFV_POINT_CARDINAL] + montage.remove_fiducials() # prevent coordinate transform with warnings.catch_warnings(): warnings.filterwarnings(action='ignore', category=RuntimeWarning, message='.*nasion not found', module='mne') raw.set_montage(montage) - montage_coord_frame = montage.get_positions()['coord_frame'] - if montage_coord_frame != 'head': - for ch in raw.info['chs']: - ch['coord_frame'] = _str_to_frame[montage_coord_frame] - for d in raw.info['dig']: - d['coord_frame'] = _str_to_frame[montage_coord_frame] + for ch in raw.info['chs']: + ch['coord_frame'] = _str_to_frame[montage_coord_frame] + for d in raw.info['dig']: + d['coord_frame'] = _str_to_frame[montage_coord_frame] + with raw.info._unlock(): # add back fiducials + raw.info['dig'] = fids + raw.info['dig'] # get the accepted mne-python coordinate frames coord_frame_int = int(montage.dig[0]['coord_frame']) @@ -405,9 +408,8 @@ def _write_dig_bids(bids_path, raw, montage=None, acpc_aligned=False, coord_frame = MNE_TO_BIDS_FRAMES.get(mne_coord_frame, None) if coord_frame == 'CapTrak' and bids_path.datatype in ('eeg', 'nirs'): - coords = _extract_landmarks(raw.info['dig']) - landmarks = set(['RPA', 'NAS', 'LPA']) == set(list(coords.keys())) - if not landmarks: + pos = raw.get_montage().get_positions() + if any([pos[fid_key] is None for fid_key in ('nasion', 'lpa', 'rpa')]): raise RuntimeError("'head' coordinate frame must contain nasion " "and left and right pre-auricular point " "landmarks") diff --git a/mne_bids/read.py b/mne_bids/read.py index 87834a6b4..8f72aec37 100644 --- a/mne_bids/read.py +++ b/mne_bids/read.py @@ -28,7 +28,7 @@ from mne_bids.config import (ALLOWED_DATATYPE_EXTENSIONS, ANNOTATIONS_TO_KEEP, reader, _map_options) -from mne_bids.utils import _extract_landmarks, _get_ch_type_mapping, verbose +from mne_bids.utils import _get_ch_type_mapping, verbose from mne_bids.path import (BIDSPath, _parse_ext, _find_matching_sidecar, _infer_datatype, get_bids_path_from_fname) @@ -964,10 +964,20 @@ def get_head_mri_trans(bids_path, extra_params=None, t1_bids_path=None, extra_params['allow_maxshield'] = True raw = read_raw_bids(bids_path=meg_bids_path, extra_params=extra_params) - meg_coords_dict = _extract_landmarks(raw.info['dig']) - meg_landmarks = np.asarray((meg_coords_dict['LPA'], - meg_coords_dict['NAS'], - meg_coords_dict['RPA'])) + + if (raw.get_montage() is None or + raw.get_montage().get_positions() is None or + any([raw.get_montage().get_positions()[fid_key] is None + for fid_key in ('nasion', 'lpa', 'rpa')])): + raise RuntimeError( + f'Could not extract fiducial points from ``raw`` file: ' + f'{meg_bids_path}\n\n' + f'The ``raw`` file SHOULD contain digitization points ' + 'for the nasion and left and right pre-auricular points ' + 'but none were found' + ) + pos = raw.get_montage().get_positions() + meg_landmarks = np.asarray((pos['lpa'], pos['nasion'], pos['rpa'])) # Given the two sets of points, fit the transform trans_fitted = fit_matched_points(src_pts=meg_landmarks, diff --git a/mne_bids/tests/test_dig.py b/mne_bids/tests/test_dig.py index 37d6ea6a2..2bdf5d42f 100644 --- a/mne_bids/tests/test_dig.py +++ b/mne_bids/tests/test_dig.py @@ -73,6 +73,40 @@ def test_dig_io(tmp_path): write_raw_bids(raw_test, bids_path) +def test_dig_pixels(tmp_path): + """Test dig stored correctly for the Pixels coordinate frame.""" + bids_root = tmp_path / 'bids1' + + # test coordinates in pixels + bids_path = _bids_path.copy().update( + root=bids_root, datatype='ieeg', space='Pixels') + os.makedirs(op.join(bids_root, 'sub-01', 'ses-01', bids_path.datatype), + exist_ok=True) + raw_test = raw.copy() + raw_test.pick_types(eeg=True) + raw_test.del_proj() + raw_test.set_channel_types({ch: 'ecog' for ch in raw_test.ch_names}) + + mnt = raw_test.get_montage() + # fake transform to pixel coordinates + mnt.apply_trans(mne.transforms.Transform('head', 'unknown')) + _write_dig_bids(bids_path, raw_test, mnt) + electrodes_path = bids_path.copy().update( + task=None, run=None, suffix='electrodes', extension='.tsv') + coordsystem_path = bids_path.copy().update( + task=None, run=None, suffix='coordsystem', extension='.json') + with pytest.warns(RuntimeWarning, + match='not an MNE-Python coordinate frame'): + _read_dig_bids(electrodes_path, coordsystem_path, + bids_path.datatype, raw_test) + mnt2 = raw_test.get_montage() + assert mnt2.get_positions()['coord_frame'] == 'unknown' + np.testing.assert_array_almost_equal( + np.array(list(mnt.get_positions()['ch_pos'].values())), + np.array(list(mnt2.get_positions()['ch_pos'].values())) + ) + + @pytest.mark.filterwarnings('ignore:The unit for chann*.:RuntimeWarning:mne') def test_dig_template(tmp_path): """Test that eeg and ieeg dig are stored properly.""" diff --git a/mne_bids/tests/test_read.py b/mne_bids/tests/test_read.py index c60f5f024..8de7ad43f 100644 --- a/mne_bids/tests/test_read.py +++ b/mne_bids/tests/test_read.py @@ -260,6 +260,20 @@ def test_get_head_mri_trans(tmp_path): fs_subject='sample', fs_subjects_dir=subjects_dir) + # test raw with no fiducials to provoke error + t1w_bids_path = write_anat( # put back + t1w_mgh, bids_path=t1w_bids_path, landmarks=landmarks, overwrite=True + ) + montage = raw.get_montage() + montage.remove_fiducials() + raw_test = raw.copy() + raw_test.set_montage(montage) + raw_test.save(bids_path.fpath, overwrite=True) + + with pytest.raises(RuntimeError, match='Could not extract fiducial'): + get_head_mri_trans(bids_path=bids_path, fs_subject='sample', + fs_subjects_dir=subjects_dir) + # test we are permissive for different casings of landmark names in the # sidecar, and also accept "nasion" instead of just "NAS" raw = _read_raw_fif(raw_fname) diff --git a/mne_bids/utils.py b/mne_bids/utils.py index 43320e8e1..30e116d12 100644 --- a/mne_bids/utils.py +++ b/mne_bids/utils.py @@ -15,7 +15,6 @@ import numpy as np from mne.channels import make_standard_montage -from mne.io.constants import FIFF from mne.io.kit.kit import get_kit_info from mne.io.pick import pick_types from mne.utils import warn, logger, verbose, check_version @@ -279,30 +278,6 @@ def _infer_eeg_placement_scheme(raw): return placement_scheme -def _extract_landmarks(dig): - """Extract NAS, LPA, and RPA from raw.info['dig'].""" - coords = dict() - coord_frame = dict() - - landmarks = {d['ident']: d for d in dig - if d['kind'] == FIFF.FIFFV_POINT_CARDINAL} - if landmarks: - if FIFF.FIFFV_POINT_NASION in landmarks: - coords['NAS'] = landmarks[FIFF.FIFFV_POINT_NASION]['r'].tolist() - coord_frame['NAS'] = (landmarks[FIFF.FIFFV_POINT_NASION] - ['coord_frame']) - if FIFF.FIFFV_POINT_LPA in landmarks: - coords['LPA'] = landmarks[FIFF.FIFFV_POINT_LPA]['r'].tolist() - coord_frame['LPA'] = landmarks[FIFF.FIFFV_POINT_LPA]['coord_frame'] - if FIFF.FIFFV_POINT_RPA in landmarks: - coords['RPA'] = landmarks[FIFF.FIFFV_POINT_RPA]['r'].tolist() - coord_frame['RPA'] = landmarks[FIFF.FIFFV_POINT_RPA]['coord_frame'] - - assert len(set(coord_frame.values())) < 2 # must all be the same - - return coords - - def _scale_coord_to_meters(coord, unit): """Scale units to meters (mne-python default).""" if unit == 'cm':