ENH : add support for landmarks as dict in write_anat and kind param …

mne_bids/read.py

@@ -761,7 +761,8 @@ def read_raw_bids(bids_path, extra_params=None, verbose=None):
 
 @verbose
 def get_head_mri_trans(bids_path, extra_params=None, t1_bids_path=None,
-                       fs_subject=None, fs_subjects_dir=None, verbose=None):
+                       fs_subject=None, fs_subjects_dir=None, kind="",
+                       verbose=None):
     """Produce transformation matrix from MEG and MRI landmark points.
 
     Will attempt to read the landmarks of Nasion, LPA, and RPA from the sidecar
@@ -797,6 +798,13 @@ def get_head_mri_trans(bids_path, extra_params=None, t1_bids_path=None,
         ``SUBJECTS_DIR`` environment variable.
 
         .. versionadded:: 0.8
+    kind : str | None
+        The suffix of the fiducials names to be used. If None, no suffix
+        is appended and one uses the landmarks called Nasion/NAS, LPA, and RPA.
+        The kind should be without the ``_`` inserted between landmark name
+        and the suffix kind.
+
+        .. versionadded:: 0.10
     %(verbose)s
 
     Returns
@@ -863,14 +871,15 @@ def get_head_mri_trans(bids_path, extra_params=None, t1_bids_path=None,
     mri_coords_dict = t1w_json.get('AnatomicalLandmarkCoordinates', dict())
 
     # landmarks array: rows: [LPA, NAS, RPA]; columns: [x, y, z]
+    suffix = f"_{kind}" if kind else ""
     mri_landmarks = np.full((3, 3), np.nan)
     for landmark_name, coords in mri_coords_dict.items():
-        if landmark_name.upper() == 'LPA':
+        if landmark_name.upper() == ('LPA' + suffix).upper():
             mri_landmarks[0, :] = coords
-        elif landmark_name.upper() == 'RPA':
+        elif landmark_name.upper() == ('RPA' + suffix).upper():
             mri_landmarks[2, :] = coords
-        elif (landmark_name.upper() == 'NAS' or
-              landmark_name.lower() == 'nasion'):
+        elif (landmark_name.upper() == ('NAS' + suffix).upper() or
+              landmark_name.lower() == ('nasion' + suffix).lower()):
             mri_landmarks[1, :] = coords
         else:
             continue

mne_bids/tests/test_read.py

@@ -317,6 +317,22 @@ def test_get_head_mri_trans(tmp_path):
             bids_path=bids_path, fs_subject='bad',
             fs_subjects_dir=subjects_dir)
 
+    # Case 3: write with suffix for kind
+    landmarks2 = landmarks.copy()
+    landmarks2.dig[0]['r'] *= -1
+    write_anat(t1w_mgh, bids_path=t1_bids_path, overwrite=True,
+               landmarks={"session1": landmarks, "session2": landmarks2})
+    read_trans1 = get_head_mri_trans(
+        bids_path=meg_bids_path, t1_bids_path=t1_bids_path,
+        fs_subject='sample', fs_subjects_dir=subjects_dir,
+        kind="session1")
+    assert np.allclose(trans['trans'], read_trans1['trans'])
+    read_trans2 = get_head_mri_trans(
+        bids_path=meg_bids_path, t1_bids_path=t1_bids_path,
+        fs_subject='sample', fs_subjects_dir=subjects_dir,
+        kind="session2")
+    assert not np.allclose(trans['trans'], read_trans2['trans'])
+
 
 def test_handle_events_reading(tmp_path):
     """Test reading events from a BIDS events.tsv file."""

mne_bids/write.py

@@ -1774,6 +1774,42 @@ def get_anat_landmarks(image, info, trans, fs_subject, fs_subjects_dir=None):
     return landmarks
 
 
+def _get_landmarks(landmarks, image_nii, kind=''):
+    import nibabel as nib
+    if isinstance(landmarks, (str, Path)):
+        landmarks, coord_frame = read_fiducials(landmarks)
+        landmarks = np.array([landmark['r'] for landmark in
+                              landmarks], dtype=float)  # unpack
+    else:
+        # Prepare to write the sidecar JSON, extract MEG landmarks
+        coords_dict, coord_frame = _get_fid_coords(landmarks.dig)
+        landmarks = np.asarray((coords_dict['lpa'],
+                                coords_dict['nasion'],
+                                coords_dict['rpa']))
+
+    # check if coord frame is supported
+    if coord_frame not in (FIFF.FIFFV_MNE_COORD_MRI_VOXEL,
+                           FIFF.FIFFV_MNE_COORD_RAS):
+        raise ValueError(f'Coordinate frame not supported: {coord_frame}')
+
+    # convert to voxels from scanner RAS to voxels
+    if coord_frame == FIFF.FIFFV_MNE_COORD_RAS:
+        # Make MGH image for header properties
+        img_mgh = nib.MGHImage(image_nii.dataobj, image_nii.affine)
+        landmarks = _mri_scanner_ras_to_mri_voxels(
+            landmarks * 1e3, img_mgh)
+
+    suffix = f"_{kind}" if kind else ""
+
+    # Write sidecar.json
+    img_json = {
+        'LPA' + suffix: list(landmarks[0, :]),
+        'NAS' + suffix: list(landmarks[1, :]),
+        'RPA' + suffix: list(landmarks[2, :])
+    }
+    return img_json
+
+
 @verbose
 def write_anat(image, bids_path, landmarks=None, deface=False, overwrite=False,
                verbose=None):
@@ -1800,12 +1836,16 @@ def write_anat(image, bids_path, landmarks=None, deface=False, overwrite=False,
         **must** have the ``root`` and ``subject`` attributes set.
         The suffix is assumed to be ``'T1w'`` if not present. It can
         also be ``'FLASH'``, for example, to indicate FLASH MRI.
-    landmarks : mne.channels.DigMontage | str | None
+    landmarks : mne.channels.DigMontage | path-like | dict | None
         The montage or path to a montage with landmarks that can be
         passed to provide information for defacing. Landmarks can be determined
         from the head model using `mne coreg` GUI, or they can be determined
-        from the MRI using freeview.  If ``None`` and no ``trans`` parameter
-        is passed, no sidecar JSON file will be created.
+        from the MRI using freeview.  If a dict is passed then the values
+        must be instances of ``DigMontage`` or path-like pointing to filenames,
+        and the keys of the dict must be strings (e.g. ``'session1'``) which
+        will be used as naming suffix of the points in the sidecar JSON file.
+        If ``None`` and no ``trans`` parameter is passed, no sidecar JSON
+        file will be created.
     deface : bool | dict
         If False, no defacing is performed.
         If True, deface with default parameters.
@@ -1876,35 +1916,14 @@ def write_anat(image, bids_path, landmarks=None, deface=False, overwrite=False,
 
     # Check if we have necessary conditions for writing a sidecar JSON
     if write_sidecar:
-        if isinstance(landmarks, str):
-            landmarks, coord_frame = read_fiducials(landmarks)
-            landmarks = np.array([landmark['r'] for landmark in
-                                  landmarks], dtype=float)  # unpack
-        else:
-            # Prepare to write the sidecar JSON, extract MEG landmarks
-            coords_dict, coord_frame = _get_fid_coords(landmarks.dig)
-            landmarks = np.asarray((coords_dict['lpa'],
-                                    coords_dict['nasion'],
-                                    coords_dict['rpa']))
-
-        # check if coord frame is supported
-        if coord_frame not in (FIFF.FIFFV_MNE_COORD_MRI_VOXEL,
-                               FIFF.FIFFV_MNE_COORD_RAS):
-            raise ValueError(f'Coordinate frame not supported: {coord_frame}')
-
-        # convert to voxels from scanner RAS to voxels
-        if coord_frame == FIFF.FIFFV_MNE_COORD_RAS:
-            # Make MGH image for header properties
-            img_mgh = nib.MGHImage(image_nii.dataobj, image_nii.affine)
-            landmarks = _mri_scanner_ras_to_mri_voxels(
-                landmarks * 1e3, img_mgh)
-
-        # Write sidecar.json
-        img_json = dict()
-        img_json['AnatomicalLandmarkCoordinates'] = \
-            {'LPA': list(landmarks[0, :]),
-             'NAS': list(landmarks[1, :]),
-             'RPA': list(landmarks[2, :])}
+        if not isinstance(landmarks, dict):
+            landmarks = {"": landmarks}
+        img_json = {}
+        for kind, this_landmarks in landmarks.items():
+            img_json.update(
+                _get_landmarks(this_landmarks, image_nii, kind=kind)
+            )
+        img_json = {'AnatomicalLandmarkCoordinates': img_json}
         fname = bids_path.copy().update(extension='.json')
         if op.isfile(fname) and not overwrite:
             raise IOError('Wanted to write a file but it already exists and '