ENH: Add generic resampling workflow, use for anatomical

fmriprep/workflows/bold/apply.py

@@ -25,6 +25,96 @@
     from niworkflows.utils.spaces import SpatialReferences
 
 
+def init_bold_volumetric_resample_wf(
+    *,
+    fieldmap_id: ty.Optional[str] = None,
+    omp_nthreads: int = 1,
+    name: str = 'bold_volumetric_resample_wf',
+) -> pe.Workflow:
+    workflow = pe.Workflow(name=name)
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=[
+                "bold_file",
+                "ref_file",
+                # HMC
+                "motion_xfm",
+                # SDC
+                "boldref2fmap_xfm",
+                "fmap_ref",
+                "fmap_coeff",
+                "fmap_id",
+                # Anatomical
+                "boldref2anat_xfm",
+                # Template
+                "anat2std_xfm",
+            ],
+        ),
+        name='inputnode',
+    )
+
+    boldref2target = pe.Node(niu.Merge(2), name='boldref2target')
+    bold2target = pe.Node(niu.Merge(2), name='bold2target')
+    resample = pe.Node(ResampleSeries(), name="resample", n_procs=omp_nthreads)
+
+    workflow.connect([
+        (inputnode, boldref2target, [
+            ('boldref2anat_xfm', 'in1'),
+            ('anat2std_xfm', 'in2'),
+        ]),
+        (inputnode, bold2target, [('motion_xfm', 'in1')]),
+        (inputnode, resample, [
+            ('bold_file', 'in_file'),
+            ('ref_file', 'ref_file'),
+        ]),
+        (boldref2target, bold2target, [('out', 'in2')]),
+        (bold2target, resample, [('out', 'transforms')]),
+    ])  # fmt:skip
+
+    if fieldmap_id:
+        fmap_select = pe.Node(
+            KeySelect(fields=["fmap_ref", "fmap_coeff"], key=fieldmap_id),
+            name="fmap_select",
+            run_without_submitting=True,
+        )
+        distortion_params = pe.Node(
+            DistortionParameters(metadata=metadata),
+            name="distortion_params",
+            run_without_submitting=True,
+        )
+        fmap2target = pe.Node(niu.Merge(2), name='fmap2target')
+        inverses = pe.Node(niu.Function(function=_gen_inverses), name='inverses')
+
+        fmap_recon = pe.Node(ReconstructFieldmap(), name="fmap_recon")
+
+        workflow.connect([
+            (inputnode, fmap_select, [
+                ("fmap_ref", "fmap_ref"),
+                ("fmap_coeff", "fmap_coeff"),
+                ("fmap_id", "keys"),
+            ]),
+            (inputnode, distortion_params, [('bold_file', 'in_file')]),
+            (inputnode, fmap2target, [('fmapreg_xfm', 'in1')]),
+            (boldref2target, fmap2target, [('out', 'in2')]),
+            (boldref2target, inverses, [('out', 'inlist')]),
+            (fmap_select, fmap_recon, [
+                ("fmap_coeff", "in_coeffs"),
+                ("fmap_ref", "fmap_ref_file"),
+            ]),
+            (fmap2target, fmap_recon, [('out', 'transforms')]),
+            (inverses, fmap_recon, [('out', 'inverse')]),
+            # Inject fieldmap correction into resample node
+            (distortion_params, resample, [
+                ("readout_time", "ro_time"),
+                ("pe_direction", "pe_dir"),
+            ]),
+            (fmap_recon, resample, [('out_file', 'fieldmap')]),
+        ])  # fmt:skip
+
+    return workflow
+
+
 def init_bold_apply_wf(
     *,
     spaces: SpatialReferences,
@@ -49,3 +139,16 @@ def init_bold_apply_wf(
         # )
 
     return workflow
+
+
+def _gen_inverses(inlist: list) -> list[bool]:
+    """Create a list indicating the first transform should be inverted.
+
+    The input list is the collection of transforms that follow the
+    inverted one.
+    """
+    from niworkflows.utils.connections import listify
+
+    if not inlist:
+        return [True]
+    return [True] + [False] * len(listify(inlist))

fmriprep/workflows/bold/base.py

@@ -45,6 +45,7 @@
 from ...utils.meepi import combine_meepi_source
 
 # BOLD workflows
+from .apply import init_bold_volumetric_resample_wf
 from .confounds import init_bold_confs_wf, init_carpetplot_wf
 from .fit import init_bold_fit_wf, init_bold_native_wf
 from .hmc import init_bold_hmc_wf
@@ -231,27 +232,46 @@ def init_bold_wf(
     #   - Save native outputs/echos only if requested
     #
 
-    bold_native_wf = init_bold_native_wf(bold_series=bold_series, fieldmap_id=fieldmap_id)
+    bold_native_wf = init_bold_native_wf(
+        bold_series=bold_series,
+        fieldmap_id=fieldmap_id,
+        omp_nthreads=omp_nthreads,
+    )
+    bold_anat_wf = init_bold_volumetric_resample_wf(
+        fieldmap_id=fieldmap_id if not multiecho else None,
+        omp_nthreads=omp_nthreads,
+        name='bold_anat_wf',
+    )
 
     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.ref_file"),
+            ("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.fmapreg_xfm"),
+            ("outputnode.boldref2fmap_xfm", "inputnode.boldref2fmap_xfm"),
             ("outputnode.dummy_scans", "inputnode.dummy_scans"),
         ]),
+        (bold_fit_wf, bold_anat_wf, [
+            ("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
 
-    if fieldmap_id:
-        workflow.connect([
-            (inputnode, bold_native_wf, [
-                ("fmap_ref", "inputnode.fmap_ref"),
-                ("fmap_coeff", "inputnode.fmap_coeff"),
-                ("fmap_id", "inputnode.fmap_id"),
-            ]),
-        ])  # fmt:skip
-
     boldref_out = bool(nonstd_spaces.intersection(('func', 'run', 'bold', 'boldref', 'sbref')))
     echos_out = multiecho and config.execution.me_output_echos
 
@@ -277,6 +297,8 @@ def init_bold_wf(
             ]),
         ])  # fmt:skip
 
+    anat_out = bool(nonstd_spaces.intersection(('anat', 'T1w')))
+
     if multiecho:
         t2s_reporting_wf = init_t2s_reporting_wf()
 

fmriprep/workflows/bold/fit.py

@@ -593,6 +593,7 @@ def init_bold_native_wf(
     *,
     bold_series: ty.List[str],
     fieldmap_id: ty.Optional[str] = None,
+    omp_nthreads: int = 1,
     name: str = "bold_native_wf",
 ) -> pe.Workflow:
     r"""
@@ -633,7 +634,7 @@ def init_bold_native_wf(
     motion_xfm
         Affine transforms from each BOLD volume to ``hmc_boldref``, written
         as concatenated ITK affine transforms.
-    fmapreg_xfm
+    boldref2fmap_xfm
         Affine transform mapping from BOLD reference space to the fieldmap
         space, if applicable.
     fmap_id
@@ -654,14 +655,12 @@ def init_bold_native_wf(
         head motion and susceptibility distortion correction (STC, HMC, SDC)
         will all be applied to each file. For multi-echo data, the echos
         are combined to form an `optimal combination`_.
+    metadata
+        Metadata dictionary of BOLD series with the shortest echo
     motion_xfm
         Motion correction transforms for further correcting bold_minimal.
         For multi-echo data, motion correction has already been applied, so
         this will be undefined.
-    fieldmap_id
-        Fieldmap ID for further correcting bold_minimal. For multi-echo data,
-        susceptibility distortion correction has already been applied, so
-        this will be undefined.
     bold_echos
         The individual, corrected echos, suitable for use in Tedana.
         (Multi-echo only.)
@@ -721,7 +720,7 @@ def init_bold_native_wf(
                 "boldref",
                 "bold_mask",
                 "motion_xfm",
-                "fmapreg_xfm",
+                "boldref2fmap_xfm",
                 "dummy_scans",
                 # Fieldmap fit
                 "fmap_ref",
@@ -735,17 +734,19 @@ def init_bold_native_wf(
     outputnode = pe.Node(
         niu.IdentityInterface(
             fields=[
-                "bold_minimal",  # Single echo: STC; Multi-echo: optimal combination
-                "bold_native",   # STC + HMC + SDC; Multi-echo: optimal combination
-                "motion_xfm",    # motion_xfms to apply to bold_minimal (none for ME)
-                "fieldmap_id",   # fieldmap to apply to bold_minimal (none for ME)
+                "bold_minimal",
+                "bold_native",
+                "metadata",
+                # Transforms
+                "motion_xfm",
                 # Multiecho outputs
                 "bold_echos",    # Individual corrected echos
                 "t2star_map",    # T2* map
             ],  # fmt:skip
         ),
         name="outputnode",
     )
+    outputnode.inputs.metadata = metadata
 
     boldbuffer = pe.Node(
         niu.IdentityInterface(fields=["bold_file", "ro_time", "pe_dir"]), name="boldbuffer"
@@ -804,9 +805,7 @@ def init_bold_native_wf(
         ])  # fmt:skip
 
     # Resample to boldref
-    boldref_bold = pe.Node(
-        ResampleSeries(), name="boldref_bold", n_procs=config.nipype.omp_nthreads
-    )
+    boldref_bold = pe.Node(ResampleSeries(), name="boldref_bold", n_procs=omp_nthreads)
 
     workflow.connect([
         (inputnode, boldref_bold, [
@@ -825,7 +824,7 @@ def init_bold_native_wf(
         workflow.connect([
             (inputnode, boldref_fmap, [
                 ("boldref", "target_ref_file"),
-                ("fmapreg_xfm", "transforms"),
+                ("boldref2fmap_xfm", "transforms"),
             ]),
             (fmap_select, boldref_fmap, [
                 ("fmap_coeff", "in_coeffs"),
@@ -850,7 +849,7 @@ def init_bold_native_wf(
             name="bold_t2smap_wf",
         )
 
-        # Do NOT set motion_xfm or fieldmap_id on outputnode
+        # Do NOT set motion_xfm on outputnode
         # This prevents downstream resamplers from double-dipping
         workflow.connect([
             (inputnode, bold_t2s_wf, [("bold_mask", "inputnode.bold_mask")]),
@@ -870,9 +869,6 @@ def init_bold_native_wf(
             (boldref_bold, outputnode, [("out_file", "bold_native")]),
         ])  # fmt:skip
 
-        if fieldmap_id:
-            outputnode.inputs.fieldmap_id = fieldmap_id
-
     return workflow