Separate goodvoxels mask creation from fsLR resampling (#3170)

Closes #3169. This could just be a first step toward moving one or more of these workflows from fMRIPrep to niworkflows.


- Move `init_goodvoxels_bold_mask_wf` call out of `init_bold_fsLR_resampling_wf`. This will allow downstream tools (e.g., ASLPrep) to generate the goodvoxels mask from one file (the ASL time series) and apply it to others (e.g., CBF maps).
- Clean up `init_bold_fsLR_resampling_wf` (e.g., remove unused imports and add all parameters to docstring).
- Move BOLD-to-anatomical resampling workflow call from `resampling` level to `full` level, since it doesn't produce any outputs.

Author: effigies

fmriprep/workflows/bold/base.py

@@ -31,8 +31,6 @@
 """
 import typing as ty
 
-import nibabel as nb
-import numpy as np
 from nipype.interfaces import utility as niu
 from nipype.pipeline import engine as pe
 from niworkflows.utils.connections import listify
@@ -295,44 +293,20 @@ def init_bold_wf(
         fieldmap_id=fieldmap_id,
         omp_nthreads=omp_nthreads,
     )
-    bold_anat_wf = init_bold_volumetric_resample_wf(
-        metadata=all_metadata[0],
-        fieldmap_id=fieldmap_id if not multiecho else None,
-        omp_nthreads=omp_nthreads,
-        mem_gb=mem_gb,
-        name='bold_anat_wf',
-    )
-    bold_anat_wf.inputs.inputnode.resolution = "native"
 
     workflow.connect([
         (inputnode, bold_native_wf, [
             ("fmap_ref", "inputnode.fmap_ref"),
             ("fmap_coeff", "inputnode.fmap_coeff"),
             ("fmap_id", "inputnode.fmap_id"),
         ]),
-        (inputnode, bold_anat_wf, [
-            ("t1w_preproc", "inputnode.target_ref_file"),
-            ("t1w_mask", "inputnode.target_mask"),
-            ("fmap_ref", "inputnode.fmap_ref"),
-            ("fmap_coeff", "inputnode.fmap_coeff"),
-            ("fmap_id", "inputnode.fmap_id"),
-        ]),
         (bold_fit_wf, bold_native_wf, [
             ("outputnode.coreg_boldref", "inputnode.boldref"),
             ("outputnode.bold_mask", "inputnode.bold_mask"),
             ("outputnode.motion_xfm", "inputnode.motion_xfm"),
             ("outputnode.boldref2fmap_xfm", "inputnode.boldref2fmap_xfm"),
             ("outputnode.dummy_scans", "inputnode.dummy_scans"),
         ]),
-        (bold_fit_wf, bold_anat_wf, [
-            ("outputnode.coreg_boldref", "inputnode.bold_ref_file"),
-            ("outputnode.boldref2fmap_xfm", "inputnode.boldref2fmap_xfm"),
-            ("outputnode.boldref2anat_xfm", "inputnode.boldref2anat_xfm"),
-        ]),
-        (bold_native_wf, bold_anat_wf, [
-            ("outputnode.bold_minimal", "inputnode.bold_file"),
-            ("outputnode.motion_xfm", "inputnode.motion_xfm"),
-        ]),
     ])  # fmt:skip
 
     boldref_out = bool(nonstd_spaces.intersection(('func', 'run', 'bold', 'boldref', 'sbref')))
@@ -400,6 +374,35 @@ def init_bold_wf(
     if config.workflow.level == "resampling":
         return workflow
 
+    # Resample to anatomical space
+    bold_anat_wf = init_bold_volumetric_resample_wf(
+        metadata=all_metadata[0],
+        fieldmap_id=fieldmap_id if not multiecho else None,
+        omp_nthreads=omp_nthreads,
+        mem_gb=mem_gb,
+        name='bold_anat_wf',
+    )
+    bold_anat_wf.inputs.inputnode.resolution = "native"
+
+    workflow.connect([
+        (inputnode, bold_anat_wf, [
+            ("t1w_preproc", "inputnode.target_ref_file"),
+            ("t1w_mask", "inputnode.target_mask"),
+            ("fmap_ref", "inputnode.fmap_ref"),
+            ("fmap_coeff", "inputnode.fmap_coeff"),
+            ("fmap_id", "inputnode.fmap_id"),
+        ]),
+        (bold_fit_wf, bold_anat_wf, [
+            ("outputnode.coreg_boldref", "inputnode.bold_ref_file"),
+            ("outputnode.boldref2fmap_xfm", "inputnode.boldref2fmap_xfm"),
+            ("outputnode.boldref2anat_xfm", "inputnode.boldref2anat_xfm"),
+        ]),
+        (bold_native_wf, bold_anat_wf, [
+            ("outputnode.bold_minimal", "inputnode.bold_file"),
+            ("outputnode.motion_xfm", "inputnode.motion_xfm"),
+        ]),
+    ])  # fmt:skip
+
     # Full derivatives, including resampled BOLD series
     if nonstd_spaces.intersection(('anat', 'T1w')):
         ds_bold_t1_wf = init_ds_volumes_wf(
@@ -507,7 +510,11 @@ def init_bold_wf(
         ])  # fmt:skip
 
     if config.workflow.cifti_output:
-        from .resampling import init_bold_fsLR_resampling_wf, init_bold_grayords_wf
+        from .resampling import (
+            init_bold_fsLR_resampling_wf,
+            init_bold_grayords_wf,
+            init_goodvoxels_bold_mask_wf,
+        )
 
         bold_MNI6_wf = init_bold_volumetric_resample_wf(
             metadata=all_metadata[0],
@@ -518,12 +525,29 @@ def init_bold_wf(
         )
 
         bold_fsLR_resampling_wf = init_bold_fsLR_resampling_wf(
-            estimate_goodvoxels=config.workflow.project_goodvoxels,
             grayord_density=config.workflow.cifti_output,
             omp_nthreads=omp_nthreads,
             mem_gb=mem_gb["resampled"],
         )
 
+        if config.workflow.project_goodvoxels:
+            goodvoxels_bold_mask_wf = init_goodvoxels_bold_mask_wf(mem_gb["resampled"])
+
+            workflow.connect([
+                (inputnode, goodvoxels_bold_mask_wf, [("anat_ribbon", "inputnode.anat_ribbon")]),
+                (bold_anat_wf, goodvoxels_bold_mask_wf, [
+                    ("outputnode.bold_file", "inputnode.bold_file"),
+                ]),
+                (goodvoxels_bold_mask_wf, bold_fsLR_resampling_wf, [
+                    ("outputnode.goodvoxels_mask", "inputnode.volume_roi"),
+                ]),
+            ])  # fmt:skip
+
+            bold_fsLR_resampling_wf.__desc__ += """\
+A "goodvoxels" mask was applied during volume-to-surface sampling in fsLR space,
+excluding voxels whose time-series have a locally high coefficient of variation.
+"""
+
         bold_grayords_wf = init_bold_grayords_wf(
             grayord_density=config.workflow.cifti_output,
             mem_gb=1,
@@ -573,7 +597,6 @@ def init_bold_wf(
                 ("midthickness_fsLR", "inputnode.midthickness_fsLR"),
                 ("sphere_reg_fsLR", "inputnode.sphere_reg_fsLR"),
                 ("cortex_mask", "inputnode.cortex_mask"),
-                ("anat_ribbon", "inputnode.anat_ribbon"),
             ]),
             (bold_anat_wf, bold_fsLR_resampling_wf, [
                 ("outputnode.bold_file", "inputnode.bold_file"),

fmriprep/workflows/bold/resampling.py

@@ -503,7 +503,6 @@ def _calc_lower_thr(in_stats):
 
 def init_bold_fsLR_resampling_wf(
     grayord_density: ty.Literal['91k', '170k'],
-    estimate_goodvoxels: bool,
     omp_nthreads: int,
     mem_gb: float,
     name: str = "bold_fsLR_resampling_wf",
@@ -522,7 +521,6 @@ def init_bold_fsLR_resampling_wf(
 
             from fmriprep.workflows.bold.resampling import init_bold_fsLR_resampling_wf
             wf = init_bold_fsLR_resampling_wf(
-                estimate_goodvoxels=True,
                 grayord_density='92k',
                 omp_nthreads=1,
                 mem_gb=1,
@@ -532,9 +530,6 @@ def init_bold_fsLR_resampling_wf(
     ----------
     grayord_density : :class:`str`
         Either ``"91k"`` or ``"170k"``, representing the total *grayordinates*.
-    estimate_goodvoxels : :class:`bool`
-        Calculate mask excluding voxels with a locally high coefficient of variation to
-        exclude from surface resampling
     omp_nthreads : :class:`int`
         Maximum number of threads an individual process may use
     mem_gb : :class:`float`
@@ -546,35 +541,33 @@ def init_bold_fsLR_resampling_wf(
     ------
     bold_file : :class:`str`
         Path to BOLD file resampled into T1 space
-    surfaces : :class:`list` of :class:`str`
-        Path to left and right hemisphere white, pial and midthickness GIFTI surfaces
-    morphometrics : :class:`list` of :class:`str`
-        Path to left and right hemisphere morphometric GIFTI surfaces, which must include thickness
+    white : :class:`list` of :class:`str`
+        Path to left and right hemisphere white matter GIFTI surfaces.
+    pial : :class:`list` of :class:`str`
+        Path to left and right hemisphere pial GIFTI surfaces.
+    midthickness : :class:`list` of :class:`str`
+        Path to left and right hemisphere midthickness GIFTI surfaces.
+    midthickness_fsLR : :class:`list` of :class:`str`
+        Path to left and right hemisphere midthickness GIFTI surfaces in fsLR space.
     sphere_reg_fsLR : :class:`list` of :class:`str`
         Path to left and right hemisphere sphere.reg GIFTI surfaces, mapping from subject to fsLR
-    anat_ribbon : :class:`str`
-        Path to mask of cortical ribbon in T1w space, for calculating goodvoxels
+    cortex_mask : :class:`list` of :class:`str`
+        Path to left and right hemisphere cortical masks.
+    volume_roi : :class:`str` or Undefined
+        Pre-calculated goodvoxels mask. Not required.
 
     Outputs
     -------
     bold_fsLR : :class:`list` of :class:`str`
         Path to BOLD series resampled as functional GIFTI files in fsLR space
-    goodvoxels_mask : :class:`str`
-        Path to mask of voxels, excluding those with locally high coefficients of variation
 
     """
     import templateflow.api as tf
     from niworkflows.engine.workflows import LiterateWorkflow as Workflow
     from niworkflows.interfaces.utility import KeySelect
     from smriprep import data as smriprep_data
-    from smriprep.interfaces.workbench import SurfaceResample
 
-    from fmriprep.interfaces.gifti import CreateROI
-    from fmriprep.interfaces.workbench import (
-        MetricFillHoles,
-        MetricRemoveIslands,
-        VolumeToSurfaceMapping,
-    )
+    from fmriprep.interfaces.workbench import VolumeToSurfaceMapping
 
     fslr_density = "32k" if grayord_density == "91k" else "59k"
 
@@ -589,13 +582,13 @@ def init_bold_fsLR_resampling_wf(
         niu.IdentityInterface(
             fields=[
                 'bold_file',
-                'anat_ribbon',
                 'white',
                 'pial',
                 'midthickness',
                 'midthickness_fsLR',
                 'sphere_reg_fsLR',
                 'cortex_mask',
+                'volume_roi',
             ]
         ),
         name='inputnode',
@@ -614,7 +607,7 @@ def init_bold_fsLR_resampling_wf(
     )
 
     outputnode = pe.Node(
-        niu.IdentityInterface(fields=['bold_fsLR', 'goodvoxels_mask']),
+        niu.IdentityInterface(fields=['bold_fsLR']),
         name='outputnode',
     )
 
@@ -689,6 +682,7 @@ def init_bold_fsLR_resampling_wf(
         # Resample BOLD to native surface, dilate and mask
         (inputnode, volume_to_surface, [
             ('bold_file', 'volume_file'),
+            ('volume_roi', 'volume_roi'),
         ]),
         (select_surfaces, volume_to_surface, [
             ('midthickness', 'surface_file'),
@@ -715,27 +709,6 @@ def init_bold_fsLR_resampling_wf(
         (joinnode, outputnode, [('bold_fsLR', 'bold_fsLR')]),
     ])  # fmt:skip
 
-    if estimate_goodvoxels:
-        workflow.__desc__ += """\
-A "goodvoxels" mask was applied during volume-to-surface sampling in fsLR space,
-excluding voxels whose time-series have a locally high coefficient of variation.
-"""
-
-        goodvoxels_bold_mask_wf = init_goodvoxels_bold_mask_wf(mem_gb)
-
-        workflow.connect([
-            (inputnode, goodvoxels_bold_mask_wf, [
-                ("bold_file", "inputnode.bold_file"),
-                ("anat_ribbon", "inputnode.anat_ribbon"),
-            ]),
-            (goodvoxels_bold_mask_wf, volume_to_surface, [
-                ("outputnode.goodvoxels_mask", "volume_roi"),
-            ]),
-            (goodvoxels_bold_mask_wf, outputnode, [
-                ("outputnode.goodvoxels_mask", "goodvoxels_mask"),
-            ]),
-        ])  # fmt:skip
-
     return workflow