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

#955)

doc/whats_new.rst

@@ -40,6 +40,8 @@ Enhancements
 
 - Add support for CNT (Neuroscan) files in :func:`mne_bids.write_raw_bids`, by `Yorguin Mantilla`_ (:gh:`924`)
 
+- Add the ability to write multiple landmarks with :func:`mne_bids.write_anat` (e.g. to have separate landmarks for different sessions) by `Alexandre Gramfort`_ (:gh:`955`)
+
 API and behavior changes
 ^^^^^^^^^^^^^^^^^^^^^^^^
 

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,17 @@ 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 anatomical landmark names in the JSON sidecar.
+        A suffix might be present e.g. to distinguish landmarks between
+        sessions. If provided, should not include a leading underscore ``_``.
+        For example, if the landmark names in the JSON sidecar file are
+        ``LPA_ses-1``, ``RPA_ses-1``, ``NAS_ses-1``, you should pass
+        ``'ses-1'`` here.
+        If ``None``, no suffix is appended, the landmarks named
+        ``Nasion`` (or ``NAS``), ``LPA``, and ``RPA`` will be used.
+
+        .. versionadded:: 0.10
     %(verbose)s
 
     Returns
@@ -863,14 +875,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,25 @@ 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
+    landmarks2.save(tmp_path / 'landmarks2.fif')
+    landmarks2 = tmp_path / 'landmarks2.fif'
+    write_anat(t1w_mgh, bids_path=t1_bids_path, overwrite=True,
+               deface=True,
+               landmarks={"coreg": landmarks, "deface": 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="coreg")
+    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="deface")
+    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, landmarks
+
+
 @verbose
 def write_anat(image, bids_path, landmarks=None, deface=False, overwrite=False,
                verbose=None):
@@ -1800,16 +1836,23 @@ 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 dictionary is passed, then the
+        values must be instances of :class:`~mne.channels.DigMontage` or
+        path-like objects pointing to a :class:`~mne.channels.DigMontage`
+        stored on disk, and the keys of the must be strings
+        (e.g. ``'session-1'``) which will be used as naming suffix for the
+        landmarks in the sidecar JSON file. If ``None``, no sidecar JSON file
+        will be created.
     deface : bool | dict
         If False, no defacing is performed.
-        If True, deface with default parameters.
-        `trans` and `raw` must not be `None` if True.
+        If ``True``, deface with default parameters using the provided
+        ``landmarks``. If multiple landmarks are provided, will
+        use the ones with the suffix ``'deface'``; if no landmarks with this
+        suffix exist, will use the first ones in the ``landmarks`` dictionary.
         If dict, accepts the following keys:
 
         - `inset`: how far back in voxels to start defacing
@@ -1876,35 +1919,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)[0]
+            )
+        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 '
@@ -1913,7 +1935,12 @@ def write_anat(image, bids_path, landmarks=None, deface=False, overwrite=False,
         _write_json(fname, img_json, overwrite)
 
         if deface:
-            image_nii = _deface(image_nii, landmarks, deface)
+            landmarks_deface = landmarks.get("deface")
+            if landmarks_deface is None:
+                # Take first one if none is specified for defacing.
+                landmarks_deface = next(iter(landmarks.items()))[1]
+            _, landmarks_deface = _get_landmarks(landmarks_deface, image_nii)
+            image_nii = _deface(image_nii, landmarks_deface, deface)
 
     # Save anatomical data
     if op.exists(bids_path):