Check for MCRIBS sphere file

xcp_d/tests/test_utils_bids.py

@@ -1,6 +1,7 @@
 """Tests for the xcp_d.utils.bids module."""
 import json
 import os
+from pathlib import Path
 
 import pytest
 from bids.layout import BIDSLayout
@@ -215,7 +216,14 @@ def test_get_freesurfer_dir(datasets):
     assert os.path.isfile(sphere_file)
 
     with pytest.raises(FileNotFoundError, match="Sphere file not found at"):
-        sphere_file = xbids.get_freesurfer_sphere(fs_dir, "fail", "L")
+        xbids.get_freesurfer_sphere(fs_dir, "fail", "L")
+
+    mcribs_sphere_file = sphere_file + "2"
+    Path(mcribs_sphere_file).touch()
+    sphere_file = xbids.get_freesurfer_sphere(fs_dir, "sub-1648798153", "L")
+    assert os.path.isfile(sphere_file)
+    assert sphere_file == mcribs_sphere_file
+    os.remove(mcribs_sphere_file)
 
 
 def test_get_entity(datasets):

xcp_d/utils/bids.py

@@ -765,7 +765,7 @@ def _get_tr(img):
 
 
 def get_freesurfer_dir(fmri_dir):
-    """Find FreeSurfer derivatives associated with preprocessing pipeline.
+    """Find FreeSurfer/MCRIBS derivatives associated with preprocessing pipeline.
 
     NOTE: This is a Node function.
 
@@ -777,42 +777,47 @@ def get_freesurfer_dir(fmri_dir):
     Returns
     -------
     freesurfer_path : :obj:`str`
-        Path to FreeSurfer derivatives.
+        Path to FreeSurfer or MCRIBS derivatives.
 
     Raises
     ------
     ValueError
-        If more than one potential FreeSurfer derivative folder is found.
+        If more than one potential FreeSurfer/MCRIBS derivative folder is found.
     NotADirectoryError
-        If no FreeSurfer derivatives are found.
+        If no FreeSurfer or MCRIBS derivatives are found.
     """
     import glob
     import os
 
-    # for fMRIPrep/Nibabies versions >=20.2.1
-    freesurfer_paths = sorted(glob.glob(os.path.join(fmri_dir, "sourcedata/*freesurfer*")))
-    if len(freesurfer_paths) == 0:
-        # for fMRIPrep/Nibabies versions <20.2.1
-        freesurfer_paths = sorted(
-            glob.glob(os.path.join(os.path.dirname(fmri_dir), "*freesurfer*"))
-        )
+    from nipype import logging
 
-    if len(freesurfer_paths) == 1:
-        freesurfer_path = freesurfer_paths[0]
+    LOGGER = logging.getLogger("nipype.utils")
 
-    elif len(freesurfer_paths) > 1:
-        freesurfer_paths_str = "\n\t".join(freesurfer_paths)
-        raise ValueError(
-            "More than one candidate for FreeSurfer derivatives found. "
-            "We recommend mounting only one FreeSurfer directory in your Docker/Singularity "
-            "image. "
-            f"Detected candidates:\n\t{freesurfer_paths_str}"
-        )
+    targets = {
+        "MCRIBS BIDS-format": os.path.join(fmri_dir, "sourcedata/mcribs"),
+        "MCRIBS legacy-format": os.path.join(os.path.dirname(fmri_dir), "mcribs"),
+        "FreeSurfer BIDS-format": os.path.join(fmri_dir, "sourcedata/*freesurfer*"),
+        "FreeSurfer legacy-format": os.path.join(os.path.dirname(fmri_dir), "*freesurfer*"),
+    }
 
-    else:
-        raise NotADirectoryError("No FreeSurfer derivatives found.")
+    for key, value in targets.items():
+        path_candidates = sorted(glob.glob(value))
+        if len(path_candidates) == 1:
+            path = path_candidates[0]
+            LOGGER.debug(f"{key} derivatives found at {path}")
+            return path
+        elif len(path_candidates) > 1:
+            path_candidates_str = "\n\t".join(path_candidates)
+            raise ValueError(
+                f"More than one candidate for {key} derivatives found. "
+                "We recommend mounting only one directory in your Docker/Singularity image. "
+                f"Detected candidates:\n\t{path_candidates_str}"
+            )
 
-    return freesurfer_path
+    searches_str = "\n\t".join(targets.values())
+    raise NotADirectoryError(
+        f"No FreeSurfer or MCRIBS derivatives found with the following searches:\n\t{searches_str}"
+    )
 
 
 def get_freesurfer_sphere(freesurfer_path, subject_id, hemisphere):
@@ -846,17 +851,26 @@ def get_freesurfer_sphere(freesurfer_path, subject_id, hemisphere):
     if not subject_id.startswith("sub-"):
         subject_id = f"sub-{subject_id}"
 
-    sphere_raw = os.path.join(
+    # Nibabies + MCRIBS will have a sphere.reg2 file instead of sphere.reg.
+    sphere_file = os.path.join(
         freesurfer_path,
         subject_id,
         "surf",
-        f"{hemisphere.lower()}h.sphere.reg",
+        f"{hemisphere.lower()}h.sphere.reg2",
     )
+    if not os.path.isfile(sphere_file):
+        # Go with FreeSurfer file if MCRIBS file doesn't exist.
+        sphere_file = os.path.join(
+            freesurfer_path,
+            subject_id,
+            "surf",
+            f"{hemisphere.lower()}h.sphere.reg",
+        )
 
-    if not os.path.isfile(sphere_raw):
-        raise FileNotFoundError(f"Sphere file not found at '{sphere_raw}'")
+    if not os.path.isfile(sphere_file):
+        raise FileNotFoundError(f"Sphere file not found at '{sphere_file}'")
 
-    return sphere_raw
+    return sphere_file
 
 
 def get_entity(filename, entity):

xcp_d/workflows/anatomical.py

@@ -1266,28 +1266,40 @@ def init_warp_one_hemisphere_wf(
     )
     inputnode.inputs.participant_id = participant_id
 
-    # Load the fsaverage-164k sphere
-    # NOTE: Why do we need the fsaverage mesh?
-    fsaverage_mesh = str(
-        get_template(
-            template="fsaverage",
-            space=None,
-            hemi=hemisphere,
-            density="164k",
-            desc=None,
-            suffix="sphere",
+    source_space = "fsaverage"
+    if source_space == "fsaverage":
+        # Load the fsaverage-164k sphere
+        # This links triangles to locations in fsaverage space
+        source_sphere = str(
+            get_template(
+                template="fsaverage",
+                space=None,
+                hemi=hemisphere,
+                density="164k",
+                desc=None,
+                suffix="sphere",
+            )
         )
-    )
 
-    # NOTE: Can we upload these to templateflow?
-    fs_hemisphere_to_fsLR = pkgrf(
-        "xcp_d",
-        (
-            f"data/standard_mesh_atlases/fs_{hemisphere}/"
-            f"fs_{hemisphere}-to-fs_LR_fsaverage.{hemisphere}_LR.spherical_std."
-            f"164k_fs_{hemisphere}.surf.gii"
-        ),
-    )
+        # Load the fsaverage sphere deformed to fsLR space.
+        # The triangles are the same as source_sphere, but the vertices point to locations in fsLR.
+        # NOTE: Can we upload these to templateflow?
+        source_sphere_deformed_to_target = pkgrf(
+            "xcp_d",
+            (
+                f"data/standard_mesh_atlases/fs_{hemisphere}/"
+                f"fs_{hemisphere}-to-fs_LR_fsaverage.{hemisphere}_LR.spherical_std."
+                f"164k_fs_{hemisphere}.surf.gii"
+            ),
+        )
+
+    elif source_space == "dhcp":
+        ...
+
+    elif source_space == "mcribs":
+        ...
+
+    # Load the subject's fsaverage-space sphere file
     get_freesurfer_sphere_node = pe.Node(
         Function(
             function=get_freesurfer_sphere,
@@ -1307,34 +1319,26 @@ def init_warp_one_hemisphere_wf(
     ])
     # fmt:on
 
-    # NOTE: What does this step do?
-    sphere_to_surf_gii = pe.Node(
+    # Convert the Freesurfer sphere.reg file to gifti format
+    sphere_to_gii = pe.Node(
         MRIsConvert(out_datatype="gii"),
-        name="sphere_to_surf_gii",
+        name="sphere_to_gii",
         mem_gb=mem_gb,
         n_procs=omp_nthreads,
     )
+    workflow.connect([(get_freesurfer_sphere_node, sphere_to_gii, [("sphere_raw", "in_file")])])
 
-    # fmt:off
-    workflow.connect([
-        (get_freesurfer_sphere_node, sphere_to_surf_gii, [("sphere_raw", "in_file")]),
-    ])
-    # fmt:on
-
-    # NOTE: What does this step do?
-    surface_sphere_project_unproject = pe.Node(
+    # Project the fsaverage-space Freesurfer file to fsLR space.
+    # It will retain the same triangles (topology) as before,
+    # but the vertices (coordinates) will be updated.
+    project_sphere_to_fslr = pe.Node(
         SurfaceSphereProjectUnproject(
-            sphere_project_to=fsaverage_mesh,
-            sphere_unproject_from=fs_hemisphere_to_fsLR,
+            sphere_project_to=source_sphere,
+            sphere_unproject_from=source_sphere_deformed_to_target,
         ),
-        name="surface_sphere_project_unproject",
+        name="project_sphere_to_fslr",
     )
-
-    # fmt:off
-    workflow.connect([
-        (sphere_to_surf_gii, surface_sphere_project_unproject, [("converted", "in_file")]),
-    ])
-    # fmt:on
+    workflow.connect([(sphere_to_gii, project_sphere_to_fslr, [("converted", "in_file")])])
 
     fsLR_sphere = str(
         get_template(
@@ -1347,8 +1351,7 @@ def init_warp_one_hemisphere_wf(
         )
     )
 
-    # resample the surfaces to fsLR-32k
-    # NOTE: Does that mean the data are in fsLR-164k before this?
+    # Resample the surfaces to fsLR-32k, using the subject's new fsLR-32k-deformed sphere.
     resample_to_fsLR32k = pe.MapNode(
         CiftiSurfaceResample(
             new_sphere=fsLR_sphere,
@@ -1363,7 +1366,7 @@ def init_warp_one_hemisphere_wf(
     # fmt:off
     workflow.connect([
         (inputnode, resample_to_fsLR32k, [("hemi_files", "in_file")]),
-        (surface_sphere_project_unproject, resample_to_fsLR32k, [("out_file", "current_sphere")]),
+        (project_sphere_to_fslr, resample_to_fsLR32k, [("out_file", "current_sphere")]),
     ])
     # fmt:on