Merge pull request #363 from effigies/enh/simplify-ribbon

RF: Replace most of anat_ribbon_wf with a Python function

.circleci/config.yml

@@ -282,6 +282,8 @@ jobs:
 
   test:
     <<: *machine_defaults
+    environment:
+      - FS_LICENSE: /tmp/fslicense/license.txt
     steps:
       - attach_workspace:
           at: /tmp
@@ -811,6 +813,7 @@ workflows:
           context:
             - nipreps-common
           requires:
+            - get_data
             - build
           filters:
             branches:

smriprep/interfaces/surf.py

@@ -157,6 +157,32 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class MakeRibbonInputSpec(TraitedSpec):
+    white_distvols = traits.List(
+        File(exists=True), minlen=2, maxlen=2, desc="White matter distance volumes"
+    )
+    pial_distvols = traits.List(
+        File(exists=True), minlen=2, maxlen=2, desc="Pial matter distance volumes"
+    )
+
+
+class MakeRibbonOutputSpec(TraitedSpec):
+    ribbon = File(desc="Binary ribbon mask")
+
+
+class MakeRibbon(SimpleInterface):
+    """Create a binary ribbon mask from white and pial distance volumes."""
+
+    input_spec = MakeRibbonInputSpec
+    output_spec = MakeRibbonOutputSpec
+
+    def _run_interface(self, runtime):
+        self._results["ribbon"] = make_ribbon(
+            self.inputs.white_distvols, self.inputs.pial_distvols, newpath=runtime.cwd
+        )
+        return runtime
+
+
 def normalize_surfs(
     in_file: str, transform_file: str | None, newpath: Optional[str] = None
 ) -> str:
@@ -209,7 +235,7 @@ def normalize_surfs(
             for RAS in "RAS":
                 pointset.meta.pop(f"VolGeom{XYZC}_{RAS}", None)
 
-    if newpath is not None:
+    if newpath is None:
         newpath = os.getcwd()
     out_file = os.path.join(newpath, fname)
     img.to_filename(out_file)
@@ -233,8 +259,32 @@ def fix_gifti_metadata(in_file: str, newpath: Optional[str] = None) -> str:
     if pointset.meta.get("GeometricType") == "Sphere":
         pointset.meta["GeometricType"] = "Spherical"
 
-    if newpath is not None:
+    if newpath is None:
         newpath = os.getcwd()
     out_file = os.path.join(newpath, os.path.basename(in_file))
     img.to_filename(out_file)
     return out_file
+
+
+def make_ribbon(
+    white_distvols: list[str],
+    pial_distvols: list[str],
+    newpath: Optional[str] = None,
+) -> str:
+    base_img = nb.load(white_distvols[0])
+    header = base_img.header
+    header.set_data_dtype("uint8")
+
+    ribbons = [
+        (np.array(nb.load(white).dataobj) > 0) & (np.array(nb.load(pial).dataobj) < 0)
+        for white, pial in zip(white_distvols, pial_distvols)
+    ]
+
+    if newpath is None:
+        newpath = os.getcwd()
+    out_file = os.path.join(newpath, "ribbon.nii.gz")
+
+    ribbon_data = ribbons[0] | ribbons[1]
+    ribbon = base_img.__class__(ribbon_data, base_img.affine, header)
+    ribbon.to_filename(out_file)
+    return out_file

smriprep/interfaces/tests/test_surf.py

@@ -0,0 +1,36 @@
+import nibabel as nb
+from nipype.pipeline import engine as pe
+import numpy as np
+
+from ...data import load_resource
+from ..surf import MakeRibbon
+
+
+def test_MakeRibbon(tmp_path):
+    res_template = "{path}/sub-fsaverage_res-4_hemi-{hemi}_desc-cropped_{surf}dist.nii.gz"
+    white, pial = [
+        [
+            load_resource(
+                res_template.format(path="../interfaces/tests/data", hemi=hemi, surf=surf)
+            )
+            for hemi in "LR"
+        ]
+        for surf in ("wm", "pial")
+    ]
+
+    make_ribbon = pe.Node(
+        MakeRibbon(white_distvols=white, pial_distvols=pial),
+        name="make_ribbon",
+        base_dir=tmp_path,
+    )
+
+    result = make_ribbon.run()
+
+    ribbon = nb.load(result.outputs.ribbon)
+    expected = nb.load(
+        load_resource("../interfaces/tests/data/sub-fsaverage_res-4_desc-cropped_ribbon.nii.gz")
+    )
+
+    assert ribbon.shape == expected.shape
+    assert np.allclose(ribbon.affine, expected.affine)
+    assert np.array_equal(ribbon.dataobj, expected.dataobj)

smriprep/workflows/anatomical.py

@@ -310,7 +310,7 @@ def init_anat_preproc_wf(
                 ('outputnode.fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'),
             ]),
             (anat_fit_wf, anat_ribbon_wf, [
-                ('outputnode.t1w_mask', 'inputnode.t1w_mask'),
+                ('outputnode.t1w_mask', 'inputnode.ref_file'),
             ]),
             (surface_derivatives_wf, anat_ribbon_wf, [
                 ('outputnode.white', 'inputnode.white'),

smriprep/workflows/surfaces.py

@@ -31,13 +31,14 @@
 from nipype.pipeline import engine as pe
 from nipype.interfaces.base import Undefined
 from nipype.interfaces import (
-    fsl,
     io as nio,
     utility as niu,
     freesurfer as fs,
     workbench as wb,
 )
 
+from smriprep.interfaces.surf import MakeRibbon
+
 from ..data import load_resource
 from ..interfaces.freesurfer import ReconAll, MakeMidthickness
 
@@ -1045,27 +1046,39 @@ def _sel(x):
 
 
 def init_anat_ribbon_wf(name="anat_ribbon_wf"):
+    """Create anatomical ribbon mask
+
+    Workflow Graph
+
+        .. workflow::
+            :graph2use: orig
+            :simple_form: yes
+
+            from smriprep.workflows.surfaces import init_anat_ribbon_wf
+            wf = init_anat_ribbon_wf()
+
+    Inputs
+    ------
+    white
+        Left and right gray/white surfaces (as GIFTI files)
+    pial
+        Left and right pial surfaces (as GIFTI files)
+    ref_file
+        Reference image (one 3D volume) to define the target space
+
+    Outputs
+    -------
+    anat_ribbon
+        Cortical gray matter mask, sampled into ``ref_file`` space
+    """
     DEFAULT_MEMORY_MIN_GB = 0.01
     workflow = pe.Workflow(name=name)
 
     inputnode = pe.Node(
-        niu.IdentityInterface(
-            fields=[
-                "white",
-                "pial",
-                "t1w_mask",
-            ]
-        ),
+        niu.IdentityInterface(fields=["white", "pial", "ref_file"]),
         name="inputnode",
     )
-    outputnode = pe.Node(
-        niu.IdentityInterface(
-            fields=[
-                "anat_ribbon",
-            ]
-        ),
-        name="outputnode",
-    )
+    outputnode = pe.Node(niu.IdentityInterface(fields=["anat_ribbon"]), name="outputnode")
 
     create_wm_distvol = pe.MapNode(
         CreateSignedDistanceVolume(),
@@ -1079,102 +1092,22 @@ def init_anat_ribbon_wf(name="anat_ribbon_wf"):
         name="create_pial_distvol",
     )
 
-    thresh_wm_distvol = pe.MapNode(
-        fsl.maths.MathsCommand(args="-thr 0 -bin -mul 255"),
-        iterfield=["in_file"],
-        name="thresh_wm_distvol",
-        mem_gb=DEFAULT_MEMORY_MIN_GB,
-    )
-
-    uthresh_pial_distvol = pe.MapNode(
-        fsl.maths.MathsCommand(args="-uthr 0 -abs -bin -mul 255"),
-        iterfield=["in_file"],
-        name="uthresh_pial_distvol",
-        mem_gb=DEFAULT_MEMORY_MIN_GB,
-    )
-
-    bin_wm_distvol = pe.MapNode(
-        fsl.maths.UnaryMaths(operation="bin"),
-        iterfield=["in_file"],
-        name="bin_wm_distvol",
-        mem_gb=DEFAULT_MEMORY_MIN_GB,
-    )
-
-    bin_pial_distvol = pe.MapNode(
-        fsl.maths.UnaryMaths(operation="bin"),
-        iterfield=["in_file"],
-        name="bin_pial_distvol",
-        mem_gb=DEFAULT_MEMORY_MIN_GB,
-    )
-
-    split_wm_distvol = pe.Node(
-        niu.Split(splits=[1, 1]),
-        name="split_wm_distvol",
-        mem_gb=DEFAULT_MEMORY_MIN_GB,
-        run_without_submitting=True,
-    )
-
-    merge_wm_distvol_no_flatten = pe.Node(
-        niu.Merge(2),
-        no_flatten=True,
-        name="merge_wm_distvol_no_flatten",
-        mem_gb=DEFAULT_MEMORY_MIN_GB,
-        run_without_submitting=True,
-    )
-
-    make_ribbon_vol = pe.MapNode(
-        fsl.maths.MultiImageMaths(op_string="-mas %s -mul 255"),
-        iterfield=["in_file", "operand_files"],
-        name="make_ribbon_vol",
-        mem_gb=DEFAULT_MEMORY_MIN_GB,
-    )
-
-    bin_ribbon_vol = pe.MapNode(
-        fsl.maths.UnaryMaths(operation="bin"),
-        iterfield=["in_file"],
-        name="bin_ribbon_vol",
-        mem_gb=DEFAULT_MEMORY_MIN_GB,
-    )
-
-    split_squeeze_ribbon_vol = pe.Node(
-        niu.Split(splits=[1, 1], squeeze=True),
-        name="split_squeeze_ribbon",
-        mem_gb=DEFAULT_MEMORY_MIN_GB,
-        run_without_submitting=True,
-    )
-
-    combine_ribbon_vol_hemis = pe.Node(
-        fsl.maths.BinaryMaths(operation="add"),
-        name="combine_ribbon_vol_hemis",
-        mem_gb=DEFAULT_MEMORY_MIN_GB,
-    )
+    make_ribbon = pe.Node(MakeRibbon(), name="make_ribbon", mem_gb=DEFAULT_MEMORY_MIN_GB)
 
-    # make HCP-style ribbon volume in T1w space
     # fmt: off
     workflow.connect(
         [
             (inputnode, create_wm_distvol, [
                 ("white", "surf_file"),
-                ("t1w_mask", "ref_file"),
+                ("ref_file", "ref_file"),
             ]),
             (inputnode, create_pial_distvol, [
                 ("pial", "surf_file"),
-                ("t1w_mask", "ref_file"),
+                ("ref_file", "ref_file"),
             ]),
-            (create_wm_distvol, thresh_wm_distvol, [("out_file", "in_file")]),
-            (create_pial_distvol, uthresh_pial_distvol, [("out_file", "in_file")]),
-            (thresh_wm_distvol, bin_wm_distvol, [("out_file", "in_file")]),
-            (uthresh_pial_distvol, bin_pial_distvol, [("out_file", "in_file")]),
-            (bin_wm_distvol, split_wm_distvol, [("out_file", "inlist")]),
-            (split_wm_distvol, merge_wm_distvol_no_flatten, [("out1", "in1")]),
-            (split_wm_distvol, merge_wm_distvol_no_flatten, [("out2", "in2")]),
-            (bin_pial_distvol, make_ribbon_vol, [("out_file", "in_file")]),
-            (merge_wm_distvol_no_flatten, make_ribbon_vol, [("out", "operand_files")]),
-            (make_ribbon_vol, bin_ribbon_vol, [("out_file", "in_file")]),
-            (bin_ribbon_vol, split_squeeze_ribbon_vol, [("out_file", "inlist")]),
-            (split_squeeze_ribbon_vol, combine_ribbon_vol_hemis, [("out1", "in_file")]),
-            (split_squeeze_ribbon_vol, combine_ribbon_vol_hemis, [("out2", "operand_file")]),
-            (combine_ribbon_vol_hemis, outputnode, [("out_file", "anat_ribbon")]),
+            (create_wm_distvol, make_ribbon, [("out_file", "white_distvols")]),
+            (create_pial_distvol, make_ribbon, [("out_file", "pial_distvols")]),
+            (make_ribbon, outputnode, [("ribbon", "anat_ribbon")]),
         ]
     )
     # fmt: on

smriprep/workflows/tests/test_surfaces.py

@@ -0,0 +1,56 @@
+import os
+from pathlib import Path
+from shutil import which
+
+import nibabel as nb
+import numpy as np
+import pytest
+from nipype.pipeline import engine as pe
+
+from ..surfaces import init_anat_ribbon_wf, init_gifti_surfaces_wf
+from ...data import load_resource
+
+
+def test_ribbon_workflow(tmp_path: Path):
+    """Create ribbon mask for fsaverage subject"""
+
+    for command in ("mris_convert", "wb_command"):
+        if not which(command):
+            pytest.skip(f"Could not find {command} in PATH")
+
+    if not os.path.exists(os.getenv('SUBJECTS_DIR')):
+        pytest.skip("Could not find $SUBJECTS_DIR")
+
+    # Low-res file that includes the fsaverage surfaces in its bounding box
+    # We will use it both as a template and a comparison.
+    test_ribbon = load_resource(
+        "../interfaces/tests/data/sub-fsaverage_res-4_desc-cropped_ribbon.nii.gz"
+    )
+
+    gifti_surfaces_wf = init_gifti_surfaces_wf(surfaces=['white', 'pial'])
+    anat_ribbon_wf = init_anat_ribbon_wf()
+    anat_ribbon_wf.inputs.inputnode.ref_file = test_ribbon
+
+    gifti_surfaces_wf.inputs.inputnode.subjects_dir = os.getenv('SUBJECTS_DIR')
+    gifti_surfaces_wf.inputs.inputnode.subject_id = 'fsaverage'
+
+    wf = pe.Workflow(name='test_ribbon_wf', base_dir=tmp_path)
+    # fmt: off
+    wf.connect([
+        (gifti_surfaces_wf, anat_ribbon_wf, [
+            ('outputnode.white', 'inputnode.white'),
+            ('outputnode.pial', 'inputnode.pial'),
+        ]),
+    ])
+    # fmt: on
+    result = wf.run()
+
+    make_ribbon = next(node for node in result.nodes() if node.name == 'make_ribbon')
+
+    expected = nb.load(test_ribbon)
+    ribbon = nb.load(make_ribbon.result.outputs.ribbon)
+
+    assert ribbon.shape == expected.shape
+    assert np.allclose(ribbon.affine, expected.affine)
+    # Mask data is binary, so we can use np.array_equal
+    assert np.array_equal(ribbon.dataobj, expected.dataobj)