ENH: Add option to exclude projecting high variance voxels to surface

.circleci/config.yml

@@ -64,7 +64,7 @@ _check_skip_job: &check_skip_job
 
 version: 2.1
 orbs:
-  docker: circleci/docker@1.6.0
+  docker: circleci/docker@2.1.4
 
 jobs:
 

nibabies/cli/parser.py

@@ -333,6 +333,15 @@ def _slice_time_ref(value, parser):
         help="Replace medial wall values with NaNs on functional GIFTI files. Only "
         "performed for GIFTI files mapped to a freesurfer subject (fsaverage or fsnative).",
     )
+    g_conf.add_argument(
+        "--project-goodvoxels",
+        required=False,
+        action="store_true",
+        default=False,
+        help="Exclude voxels whose timeseries have locally high coefficient of variation "
+        "from surface resampling. Only performed for GIFTI files mapped to a freesurfer subject "
+        "(fsaverage or fsnative).",
+    )
     g_conf.add_argument(
         "--slice-time-ref",
         required=False,

nibabies/config.py

@@ -553,6 +553,8 @@ class workflow(_Config):
     """Run FreeSurfer ``recon-all`` with the ``-logitudinal`` flag."""
     medial_surface_nan = None
     """Fill medial surface with :abbr:`NaNs (not-a-number)` when sampling."""
+    project_goodvoxels = False
+    """Exclude voxels with locally high coefficient of variation from sampling."""
     regressors_all_comps = None
     """Return all CompCor components."""
     regressors_dvars_th = None

nibabies/data/tests/config.toml

@@ -52,6 +52,7 @@ hires = true
 ignore = [ "slicetiming",]
 longitudinal = false
 medial_surface_nan = false
+project_goodvoxels = false
 regressors_all_comps = false
 regressors_dvars_th = 1.5
 regressors_fd_th = 0.5

nibabies/interfaces/workbench.py

@@ -1445,3 +1445,128 @@ class VolumeLabelImport(WBCommand):
     input_spec = VolumeLabelImportInputSpec
     output_spec = VolumeLabelImportOutputSpec
     _cmd = "wb_command -volume-label-import"
+
+
+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"

nibabies/workflows/anatomical/base.py

@@ -92,6 +92,7 @@ def init_infant_anat_wf(
     from .preproc import init_anat_average_wf
     from .registration import init_coregistration_wf
     from .segmentation import init_anat_segmentations_wf
+    from .surfaces import init_anat_ribbon_wf
 
     # for now, T1w only
     num_t1w = len(t1w) if t1w else 0
@@ -143,6 +144,7 @@ def init_infant_anat_wf(
                 "surfaces",
                 "anat_aseg",
                 "anat_aparc",
+                "anat_ribbon",
                 "template",
             ]
         ),
@@ -250,6 +252,10 @@ def init_infant_anat_wf(
         templates=spaces.get_spaces(nonstandard=False, dim=(3,)),
     )
 
+    # Anatomical ribbon file using HCP signed-distance volume method
+    # if config.workflow.project_goodvoxels:
+    anat_ribbon_wf = init_anat_ribbon_wf()
+
     # fmt:off
     wf.connect([
         (inputnode, t1w_template_wf, [("t1w", "inputnode.in_files")]),
@@ -354,6 +360,9 @@ def init_infant_anat_wf(
             ("outputnode.anat2std_xfm", "inputnode.anat2std_xfm"),
             ("outputnode.std2anat_xfm", "inputnode.std2anat_xfm"),
         ]),
+        (anat_ribbon_wf, anat_derivatives_wf, [
+            ("outputnode.anat_ribbon", "inputnode.anat_ribbon"),
+        ]),
         (anat_seg_wf, anat_derivatives_wf, [
             ("outputnode.anat_dseg", "inputnode.t1w_dseg"),
             ("outputnode.anat_tpms", "inputnode.t1w_tpms"),
@@ -404,6 +413,15 @@ def init_infant_anat_wf(
             ("outputnode.out_aparc", "anat_aparc"),
             ("outputnode.out_aseg", "anat_aseg"),
         ]),
+        (t1w_preproc_wf, anat_ribbon_wf, [
+            ("outputnode.t1w_mask", "inputnode.t1w_mask"),
+        ]),
+        (surface_recon_wf, anat_ribbon_wf, [
+            ("outputnode.surfaces", "inputnode.surfaces"),
+        ]),
+        (anat_ribbon_wf, outputnode, [
+            ("outputnode.anat_ribbon", "anat_ribbon")
+        ]),
         (surface_recon_wf, anat_reports_wf, [
             ("outputnode.subject_id", "inputnode.subject_id"),
             ("outputnode.subjects_dir", "inputnode.subjects_dir"),

nibabies/workflows/anatomical/outputs.py

@@ -338,6 +338,7 @@ def init_anat_derivatives_wf(
                 "surfaces",
                 "t1w_fs_aseg",
                 "t1w_fs_aparc",
+                "anat_ribbon",
             ]
         ),
         name="inputnode",
@@ -353,6 +354,12 @@ def init_anat_derivatives_wf(
     )
     ds_t1w_preproc.inputs.SkullStripped = False
 
+    ds_anat_ribbon = pe.Node(
+        DerivativesDataSink(base_directory=output_dir, suffix="ribbon", compress=True),
+        name="ds_anat_ribbon",
+        run_without_submitting=True,
+    )
+
     ds_t1w_mask = pe.Node(
         DerivativesDataSink(base_directory=output_dir, desc="brain", suffix="mask", compress=True),
         name="ds_t1w_mask",
@@ -384,6 +391,8 @@ def init_anat_derivatives_wf(
                                   ('source_files', 'source_file')]),
         (inputnode, ds_t1w_dseg, [('t1w_dseg', 'in_file'),
                                   ('source_files', 'source_file')]),
+        (inputnode, ds_anat_ribbon, [('anat_ribbon', 'in_file'),
+                                     ('source_files', 'source_file')]),
         (raw_sources, ds_t1w_mask, [('out', 'RawSources')]),
     ])
     # fmt: on