ENH: Add cortical ribbon workflow

smriprep/data/io_spec.json

@@ -77,6 +77,12 @@
         "hemi": ["L", "R"],
         "extension": "shape.gii",
         "suffix": ["thickness", "sulc", "curv"]
+      },
+      "anat_ribbon": {
+        "datatype": "anat",
+        "desc": "ribbon",
+        "suffix": "mask",
+        "extension": "nii.gz"
       }
     }
   },

smriprep/interfaces/workbench.py

@@ -0,0 +1,129 @@
+import os
+
+from nipype.interfaces.base import CommandLineInputSpec, File, Str, TraitedSpec, traits
+from nipype.interfaces.workbench.base import WBCommand
+
+
+class CreateSignedDistanceVolumeInputSpec(CommandLineInputSpec):
+    surf_file = File(
+        exists=True,
+        mandatory=True,
+        argstr="%s",
+        position=0,
+        desc="Input surface GIFTI file (.surf.gii)",
+    )
+    ref_file = File(
+        exists=True,
+        mandatory=True,
+        argstr="%s",
+        position=1,
+        desc="NIfTI volume in the desired output space (dims, spacing, origin)",
+    )
+    out_file = File(
+        name_source=["surf_file"],
+        name_template="%s_distvol.nii.gz",
+        argstr="%s",
+        position=2,
+        desc="Name of output volume containing signed distances",
+    )
+    out_mask = File(
+        name_source=["surf_file"],
+        name_template="%s_distmask.nii.gz",
+        argstr="-roi-out %s",
+        desc="Name of file to store a mask where the ``out_file`` has a computed value",
+    )
+    fill_value = traits.Float(
+        0.0,
+        mandatory=False,
+        usedefault=True,
+        argstr="-fill-value %f",
+        desc="value to put in all voxels that don't get assigned a distance",
+    )
+    exact_limit = traits.Float(
+        5.0,
+        usedefault=True,
+        argstr="-exact-limit %f",
+        desc="distance for exact output in mm",
+    )
+    approx_limit = traits.Float(
+        20.0,
+        usedefault=True,
+        argstr="-approx-limit %f",
+        desc="distance for approximate output in mm",
+    )
+    approx_neighborhood = traits.Int(
+        2,
+        usedefault=True,
+        argstr="-approx-neighborhood %d",
+        desc="size of neighborhood cube measured from center to face in voxels, default 2 = 5x5x5",
+    )
+    winding_method = traits.Enum(
+        "EVEN_ODD",
+        "NEGATIVE",
+        "NONZERO",
+        "NORMALS",
+        argstr="-winding %s",
+        usedefault=True,
+        desc="winding method for point inside surface test",
+    )
+
+
+class CreateSignedDistanceVolumeOutputSpec(TraitedSpec):
+    out_file = File(desc="Name of output volume containing signed distances")
+    out_mask = File(
+        desc="Name of file to store a mask where the ``out_file`` has a computed value"
+    )
+
+
+class CreateSignedDistanceVolume(WBCommand):
+    """CREATE SIGNED DISTANCE VOLUME FROM SURFACE
+    wb_command -create-signed-distance-volume
+       <surface> - the input surface
+       <refspace> - a volume in the desired output space (dims, spacing, origin)
+       <outvol> - output - the output volume
+
+       [-roi-out] - output an roi volume of where the output has a computed
+          value
+          <roi-vol> - output - the output roi volume
+
+       [-fill-value] - specify a value to put in all voxels that don't get
+          assigned a distance
+          <value> - value to fill with (default 0)
+
+       [-exact-limit] - specify distance for exact output
+          <dist> - distance in mm (default 5)
+
+       [-approx-limit] - specify distance for approximate output
+          <dist> - distance in mm (default 20)
+
+       [-approx-neighborhood] - voxel neighborhood for approximate calculation
+          <num> - size of neighborhood cube measured from center to face, in
+             voxels (default 2 = 5x5x5)
+
+       [-winding] - winding method for point inside surface test
+          <method> - name of the method (default EVEN_ODD)
+
+       Computes the signed distance function of the surface.  Exact distance is
+       calculated by finding the closest point on any surface triangle to the
+       center of the voxel.  Approximate distance is calculated starting with
+       these distances, using dijkstra's method with a neighborhood of voxels.
+       Specifying too small of an exact distance may produce unexpected results.
+       Valid specifiers for winding methods are as follows:
+
+       EVEN_ODD (default)
+       NEGATIVE
+       NONZERO
+       NORMALS
+
+       The NORMALS method uses the normals of triangles and edges, or the
+       closest triangle hit by a ray from the point.  This method may be
+       slightly faster, but is only reliable for a closed surface that does not
+       cross through itself.  All other methods count entry (positive) and exit
+       (negative) crossings of a vertical ray from the point, then counts as
+       inside if the total is odd, negative, or nonzero, respectively.
+
+    """
+
+    input_spec = CreateSignedDistanceVolumeInputSpec
+    output_spec = CreateSignedDistanceVolumeOutputSpec
+    _cmd = "wb_command -create-signed-distance-volume"

smriprep/utils/bids.py

@@ -40,7 +40,7 @@ def get_outputnode_spec():
     'anat2std_xfm', 'std2anat_xfm',
     't1w_aseg', 't1w_aparc',
     't1w2fsnative_xfm', 'fsnative2t1w_xfm',
-    'surfaces', 'morphometrics']
+    'surfaces', 'morphometrics', 'anat_ribbon']
 
     """
     spec = loads(Path(pkgrf("smriprep", "data/io_spec.json")).read_text())["queries"]

smriprep/workflows/anatomical.py

@@ -51,7 +51,7 @@
 from ..utils.misc import apply_lut as _apply_bids_lut, fs_isRunning as _fs_isRunning
 from .norm import init_anat_norm_wf
 from .outputs import init_anat_reports_wf, init_anat_derivatives_wf
-from .surfaces import init_surface_recon_wf
+from .surfaces import init_anat_ribbon_wf, init_surface_recon_wf
 
 LOGGER = logging.getLogger("nipype.workflow")
 
@@ -582,6 +582,8 @@ def _check_img(img):
         ])
         # fmt:on
 
+    # Anatomical ribbon file using HCP signed-distance volume method
+    anat_ribbon_wf = init_anat_ribbon_wf()
     # fmt:off
     workflow.connect([
         (inputnode, fs_isrunning, [
@@ -610,6 +612,11 @@ def _check_img(img):
             ('outputnode.morphometrics', 'morphometrics'),
             ('outputnode.out_aseg', 't1w_aseg'),
             ('outputnode.out_aparc', 't1w_aparc')]),
+        (surface_recon_wf, anat_ribbon_wf, [
+            ('outputnode.surfaces', 'inputnode.surfaces'),
+            ('outputnode.out_brainmask', 'inputnode.t1w_mask')]),
+        (anat_ribbon_wf, outputnode, [
+            ("outputnode.anat_ribbon", "anat_ribbon")]),
         (applyrefined, buffernode, [('out_file', 't1w_brain')]),
         (surface_recon_wf, buffernode, [
             ('outputnode.out_brainmask', 't1w_mask')]),
@@ -626,6 +633,9 @@ def _check_img(img):
             ('surfaces', 'inputnode.surfaces'),
             ('morphometrics', 'inputnode.morphometrics'),
         ]),
+        (anat_ribbon_wf, anat_derivatives_wf, [
+            ("outputnode.anat_ribbon", "inputnode.anat_ribbon"),
+        ]),
     ])
     # fmt:on
 

smriprep/workflows/outputs.py

@@ -267,6 +267,8 @@ def init_anat_derivatives_wf(
         GIFTI surfaces (gray/white boundary, midthickness, pial, inflated)
     morphometrics
         GIFTIs of cortical thickness, curvature, and sulcal depth
+    anat_ribbon
+        Cortical ribbon volume in T1w space
     t1w_fs_aseg
         FreeSurfer's aseg segmentation, in native T1w space
     t1w_fs_aparc
@@ -294,6 +296,7 @@ def init_anat_derivatives_wf(
                 "fsnative2t1w_xfm",
                 "surfaces",
                 "morphometrics",
+                "anat_ribbon",
                 "t1w_fs_aseg",
                 "t1w_fs_aparc",
             ]
@@ -637,6 +640,19 @@ def init_anat_derivatives_wf(
         name="ds_morphs",
         run_without_submitting=True,
     )
+    # Ribbon volume
+    ds_anat_ribbon = pe.Node(
+        DerivativesDataSink(
+            base_directory=output_dir,
+            desc="ribbon",
+            suffix="mask",
+            extension=".nii.gz",
+            compress=True,
+        ),
+        name="ds_anat_ribbon",
+        run_without_submitting=True,
+    )
+
     # Parcellations
     ds_t1w_fsaseg = pe.Node(
         DerivativesDataSink(
@@ -675,6 +691,9 @@ def init_anat_derivatives_wf(
                                     ('source_files', 'source_file')]),
         (inputnode, ds_t1w_fsparc, [('t1w_fs_aparc', 'in_file'),
                                     ('source_files', 'source_file')]),
+        (inputnode, ds_anat_ribbon, [('anat_ribbon', 'in_file'),
+                                     ('source_files', 'source_file')]),
+
     ])
     # fmt:on
     return workflow