MRG: Make update_anat_landmarks() accept fiducials file

doc/whats_new.rst

@@ -52,13 +52,17 @@ 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`)
+
 API and behavior changes
 ^^^^^^^^^^^^^^^^^^^^^^^^
 
 - :func:`mne_bids.update_anat_landmarks` will now by default raise an exception if the requested MRI landmarks do not already exist. Use the new ``on_missing`` parameter to control this behavior, by `Richard Höchenberger`_ (:gh:`957`)
 
 - :func:`mne_bids.get_head_mri_trans` now raises a warning if ``datatype`` or ``suffix`` of the provided electrophysiological :class:`mne_bids.BIDSPath` are not set. In the future, this will raise an exception, by `Richard Höchenberger`_(:gh:`969`)
 
+- 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`)
+
 Requirements
 ^^^^^^^^^^^^
 

mne_bids/read.py

@@ -804,9 +804,11 @@ def get_head_mri_trans(bids_path, extra_params=None, t1_bids_path=None,
         Use this parameter e.g. if the T1 scan was recorded during a different
         session than the MEG. It is even possible to point to a T1 image stored
         in an entirely different BIDS dataset than the MEG data.
-    fs_subject : str | None
-        The subject identifier used for FreeSurfer. If ``None``, defaults to
-        the ``subject`` entity in ``bids_path``.
+    fs_subject : str
+        The subject identifier used for FreeSurfer.
+
+        .. versionchanged:: 0.10
+           Does not default anymore to ``bids_path.subject`` if ``None``.
     fs_subjects_dir : path-like | None
         The FreeSurfer subjects directory. If ``None``, defaults to the
         ``SUBJECTS_DIR`` environment variable.
@@ -922,11 +924,15 @@ def get_head_mri_trans(bids_path, extra_params=None, t1_bids_path=None,
             f'{mri_coords_dict}'
         )
 
-    # The MRI landmarks are in "voxels". We need to convert the to the
-    # neuromag RAS coordinate system in order to compare the with MEG landmarks
-    # see also: `mne_bids.write.write_anat`
+    # The MRI landmarks are in "voxels". We need to convert them to the
+    # Neuromag RAS coordinate system in order to compare them with MEG
+    # landmarks. See also: `mne_bids.write.write_anat`
     if fs_subject is None:
+        warn('Passing "fs_subject=None" has been deprecated and will raise '
+             'an error in future versions. Please explicitly specify the '
+             'FreeSurfer subject name.', DeprecationWarning)
         fs_subject = f'sub-{meg_bids_path.subject}'
+
     fs_subjects_dir = get_subjects_dir(fs_subjects_dir, raise_error=False)
     fs_t1_path = Path(fs_subjects_dir) / fs_subject / 'mri' / 'T1.mgz'
     if not fs_t1_path.exists():

mne_bids/sidecar_updates.py

@@ -9,10 +9,15 @@
 import json
 from collections import OrderedDict
 
-from mne.channels import DigMontage
+import numpy as np
+
+from mne.channels import DigMontage, make_dig_montage
 from mne.utils import (
     logger, _validate_type, verbose, _check_on_missing, _on_missing
 )
+from mne.io import read_fiducials
+from mne.io.constants import FIFF
+
 from mne_bids import BIDSPath
 from mne_bids.utils import _write_json
 
@@ -141,7 +146,8 @@ def _update_sidecar(sidecar_fname, key, val):
 
 @verbose
 def update_anat_landmarks(
-    bids_path, landmarks, kind=None, on_missing='raise', verbose=None
+    bids_path, landmarks, *, fs_subject=None, fs_subjects_dir=None,
+    kind=None, on_missing='raise', verbose=None
 ):
     """Update the anatomical landmark coordinates of an MRI scan.
 
@@ -152,14 +158,29 @@ def update_anat_landmarks(
     ----------
     bids_path : BIDSPath
         Path of the MR image.
-    landmarks : mne.channels.DigMontage
+    landmarks : mne.channels.DigMontage | path-like
         An :class:`mne.channels.DigMontage` instance with coordinates for the
         nasion and left and right pre-auricular points in MRI voxel
-        coordinates.
+        coordinates. Alternatively, the path to a ``*-fiducials.fif`` file as
+        produced by the MNE-Python coregistration GUI or via
+        :func:`mne.io.write_fiducials`.
 
         .. note:: :func:`mne_bids.get_anat_landmarks` provides a convenient and
                   reliable way to generate the landmark coordinates in the
                   required coordinate system.
+
+        .. note:: If ``path-like``, ``fs_subject`` and ``fs_subjects_dir``
+                  must be provided as well.
+
+        .. versionchanged:: 0.10
+           Added support for ``path-like`` input.
+    fs_subject : str | None
+        The subject identifier used for FreeSurfer. Must be provided if
+        ``landmarks`` is ``path-like``; otherwise, it will be ignored.
+    fs_subjects_dir : path-like | None
+        The FreeSurfer subjects directory. If ``None``, defaults to the
+        ``SUBJECTS_DIR`` environment variable. Must be provided if
+        ``landmarks`` is ``path-like``; otherwise, it will be ignored.
     kind : str | None
         The suffix of the anatomical landmark names in the JSON sidecar.
         A suffix might be present e.g. to distinguish landmarks between
@@ -183,7 +204,9 @@ def update_anat_landmarks(
     .. versionadded:: 0.8
     """
     _validate_type(item=bids_path, types=BIDSPath, item_name='bids_path')
-    _validate_type(item=landmarks, types=DigMontage, item_name='landmarks')
+    _validate_type(
+        item=landmarks, types=(DigMontage, 'path-like'), item_name='landmarks'
+    )
     _check_on_missing(on_missing)
 
     # Do some path verifications and fill in some gaps the users might have
@@ -231,6 +254,19 @@ def update_anat_landmarks(
             f'bids_path. Tried the following filenames: '
             f'{", ".join([p.name for p in tried_paths])}')
 
+    if not isinstance(landmarks, DigMontage):  # it's pathlike
+        if fs_subject is None:
+            raise ValueError(
+                'You must provide the "fs_subject" parameter when passing the '
+                'path to fiducials'
+            )
+        landmarks = _get_landmarks_from_fiducials_file(
+            bids_path=bids_path,
+            fname=landmarks,
+            fs_subject=fs_subject,
+            fs_subjects_dir=fs_subjects_dir
+        )
+
     positions = landmarks.get_positions()
     coord_frame = positions['coord_frame']
     if coord_frame != 'mri_voxel':
@@ -252,7 +288,10 @@ def update_anat_landmarks(
         if coords is None:
             missing_points.append(name)
         else:
-            name_to_coords_map[name] = list(coords)
+            # Funnily, np.float64 is JSON-serializabe, while np.float32 is not!
+            # Thus, cast to float64 to avoid issues (which e.g. may arise when
+            # fiducials were read from disk!)
+            name_to_coords_map[name] = list(coords.astype('float64'))
 
     if missing_points:
         raise ValueError(
@@ -290,3 +329,47 @@ def update_anat_landmarks(
         mri_json['AnatomicalLandmarkCoordinates'][name] = coords
 
     update_sidecar_json(bids_path=bids_path_json, entries=mri_json)
+
+
+def _get_landmarks_from_fiducials_file(*, bids_path, fname, fs_subject,
+                                       fs_subjects_dir):
+    """Get anatomical landmarks from fiducials file, in MRI voxel space."""
+    # avoid dicrular imports
+    from mne_bids.write import (
+        _get_t1w_mgh, _mri_landmarks_to_mri_voxels, _get_fid_coords
+    )
+
+    digpoints, coord_frame = read_fiducials(fname)
+
+    # All of this should be guaranteed, but better be safe than sorry!
+    assert coord_frame == FIFF.FIFFV_COORD_MRI
+    assert digpoints[0]['ident'] == FIFF.FIFFV_POINT_LPA
+    assert digpoints[1]['ident'] == FIFF.FIFFV_POINT_NASION
+    assert digpoints[2]['ident'] == FIFF.FIFFV_POINT_RPA
+
+    montage_loaded = make_dig_montage(
+        lpa=digpoints[0]['r'],
+        nasion=digpoints[1]['r'],
+        rpa=digpoints[2]['r'],
+        coord_frame='mri'
+    )
+    landmark_coords_mri, _ = _get_fid_coords(dig_points=montage_loaded.dig)
+    landmark_coords_mri = np.asarray(
+        (landmark_coords_mri['lpa'],
+         landmark_coords_mri['nasion'],
+         landmark_coords_mri['rpa'])
+    )
+
+    t1w_mgh = _get_t1w_mgh(fs_subject, fs_subjects_dir)
+    landmark_coords_voxels = _mri_landmarks_to_mri_voxels(
+        mri_landmarks=landmark_coords_mri * 1000,  # in mm
+        t1_mgh=t1w_mgh
+    )
+    montage_voxels = make_dig_montage(
+        lpa=landmark_coords_voxels[0],
+        nasion=landmark_coords_voxels[1],
+        rpa=landmark_coords_voxels[2],
+        coord_frame='mri_voxel'
+    )
+
+    return montage_voxels

mne_bids/tests/test_update.py

@@ -251,3 +251,30 @@ def test_update_anat_landmarks(tmp_path):
                        match='did not contain all required cardinal points'):
         update_anat_landmarks(bids_path=bids_path_mri,
                               landmarks=landmarks_invalid)
+
+    # Test with path-like landmarks
+    fiducials_path = (data_path / 'subjects' / 'sample' / 'bem' /
+                      'sample-fiducials.fif')
+
+    update_anat_landmarks(
+        bids_path=bids_path_mri,
+        landmarks=fiducials_path,
+        fs_subject='sample',
+        fs_subjects_dir=data_path / 'subjects'
+    )
+    expected_coords_in_voxels = np.array(
+        [[68.38202,  45.24057,  43.439808],  # noqa: E241
+         [42.27006,  30.758774, 74.09837 ],  # noqa: E202, E241
+         [17.044853, 46.586075, 42.618504]]
+    )
+    mri_json = json.loads(
+        bids_path_mri_json.fpath.read_text(encoding='utf-8')
+    )
+    for landmark, expected_coords in zip(
+        ('LPA', 'NAS', 'RPA'),
+        expected_coords_in_voxels
+    ):
+        assert np.allclose(
+            mri_json['AnatomicalLandmarkCoordinates'][landmark],
+            expected_coords
+        )

mne_bids/write.py

@@ -186,13 +186,13 @@ def _channels_tsv(raw, fname, overwrite=False):
 }
 
 
-def _get_fid_coords(dig, raise_error=True):
+def _get_fid_coords(dig_points, raise_error=True):
     """Get the fiducial coordinates from a DigMontage.
 
     Parameters
     ----------
-    dig : mne.channels.DigMontage
-        The dig montage with the fiducial coordinates.
+    dig_points : array-like of DigPoint
+        The digitization points of the fiducial coordinates.
     raise_error : bool
         Whether to raise an error if the coordinates are missing or
         incorrectly formatted
@@ -207,7 +207,7 @@ def _get_fid_coords(dig, raise_error=True):
     fid_coords = Bunch(nasion=None, lpa=None, rpa=None)
     fid_coord_frames = dict()
 
-    for d in dig:
+    for d in dig_points:
         if d['kind'] == FIFF.FIFFV_POINT_CARDINAL:
             key = _cardinal_ident_mapping[d['ident']]
             fid_coords[key] = d['r']
@@ -614,9 +614,10 @@ def _mri_landmarks_to_mri_voxels(mri_landmarks, t1_mgh):
         The MRI voxel-space landmark data.
     """
     # Get landmarks in voxel space, using the T1 data
-    vox2ras_tkr = t1_mgh.header.get_vox2ras_tkr()
-    ras2vox_tkr = linalg.inv(vox2ras_tkr)
-    vox_landmarks = apply_trans(ras2vox_tkr, mri_landmarks)  # in vox
+    vox2ras_tkr_t = t1_mgh.header.get_vox2ras_tkr()
+    ras_tkr2vox_t = linalg.inv(vox2ras_tkr_t)
+    vox_landmarks = apply_trans(ras_tkr2vox_t, mri_landmarks)
+
     return vox_landmarks
 
 
@@ -1789,14 +1790,12 @@ def get_anat_landmarks(image, info, trans, fs_subject, fs_subjects_dir=None):
         The measurement information from an electrophysiology recording of
         the subject with the anatomical landmarks stored in its
         :class:`mne.channels.DigMontage`.
-    trans : mne.transforms.Transform | str
+    trans : mne.transforms.Transform | path-like
         The transformation matrix from head to MRI coordinates. Can
         also be a string pointing to a ``.trans`` file containing the
-        transformation matrix. If ``None`` and no ``landmarks`` parameter is
-        passed, no sidecar JSON file will be created.
+        transformation matrix.
     fs_subject : str
-        The subject identifier used for FreeSurfer. If ``None``, defaults to
-        the ``subject`` entity in ``bids_path``. Must be provided to write
+        The subject identifier used for FreeSurfer. Must be provided to write
         the anatomical landmarks if they are not provided in MRI voxel space.
         This is because the head coordinate of a
         :class:`mne.channels.DigMontage` is aligned using FreeSurfer surfaces.
@@ -1820,32 +1819,49 @@ def get_anat_landmarks(image, info, trans, fs_subject, fs_subjects_dir=None):
     landmarks = np.asarray((coords_dict['lpa'],
                             coords_dict['nasion'],
                             coords_dict['rpa']))
+
     # get trans and ensure it is from head to MRI
     trans, _ = _get_trans(trans, fro='head', to='mri')
     landmarks = _meg_landmarks_to_mri_landmarks(landmarks, trans)
-    fs_subjects_dir = get_subjects_dir(fs_subjects_dir, raise_error=True)
-    t1_fname = Path(fs_subjects_dir) / fs_subject / 'mri' / 'T1.mgz'
-    if not t1_fname.exists():
-        raise ValueError('Freesurfer recon-all subject folder '
-                         'is incorrect or improperly formatted, '
-                         f'got {Path(fs_subjects_dir) / fs_subject}')
-    t1w_img = _load_image(str(t1_fname), name='T1.mgz')
-    t1w_mgh = nib.MGHImage(t1w_img.dataobj, t1w_img.affine)
-    # go to T1 voxel space from surface RAS/TkReg RAS/freesurfer
+
+    # Get FS T1 image in MGH format
+    t1w_mgh = _get_t1w_mgh(fs_subject, fs_subjects_dir)
+
+    # FS MGH image: go to T1 voxel space from surface RAS/TkReg RAS/freesurfer
     landmarks = _mri_landmarks_to_mri_voxels(landmarks, t1w_mgh)
-    # go to T1 scanner space from T1 voxel space
+
+    # FS MGH image: go to T1 scanner space from T1 voxel space
     landmarks = _mri_voxels_to_mri_scanner_ras(landmarks, t1w_mgh)
+
+    # Input image: go to T1 voxel space from T1 scanner space
     if isinstance(image, BIDSPath):
         image = image.fpath
     img_nii = _load_image(image, name='image')
     img_mgh = nib.MGHImage(img_nii.dataobj, img_nii.affine)
     landmarks = _mri_scanner_ras_to_mri_voxels(landmarks, img_mgh)
+
     landmarks = mne.channels.make_dig_montage(
         lpa=landmarks[0], nasion=landmarks[1], rpa=landmarks[2],
         coord_frame='mri_voxel')
+
     return landmarks
 
 
+def _get_t1w_mgh(fs_subject, fs_subjects_dir):
+    """Return the T1w image in MGH format."""
+    import nibabel as nib
+
+    fs_subjects_dir = get_subjects_dir(fs_subjects_dir, raise_error=True)
+    t1_fname = Path(fs_subjects_dir) / fs_subject / 'mri' / 'T1.mgz'
+    if not t1_fname.exists():
+        raise ValueError('Freesurfer recon-all subject folder '
+                         'is incorrect or improperly formatted, '
+                         f'got {Path(fs_subjects_dir) / fs_subject}')
+    t1w_img = _load_image(str(t1_fname), name='T1.mgz')
+    t1w_mgh = nib.MGHImage(t1w_img.dataobj, t1w_img.affine)
+    return t1w_mgh
+
+
 def _get_landmarks(landmarks, image_nii, kind=''):
     import nibabel as nib
     if isinstance(landmarks, (str, Path)):