ENH: Restore confound generation (#3120)

Builds on #3116. Just check the last two commits.

Author: effigies

fmriprep/workflows/base.py

@@ -503,6 +503,7 @@ def init_single_subject_wf(subject_id: str):
                 ('outputnode.t1w_preproc', 'inputnode.t1w_preproc'),
                 ('outputnode.t1w_mask', 'inputnode.t1w_mask'),
                 ('outputnode.t1w_dseg', 'inputnode.t1w_dseg'),
+                ('outputnode.t1w_tpms', 'inputnode.t1w_tpms'),
                 ('outputnode.subjects_dir', 'inputnode.subjects_dir'),
                 ('outputnode.subject_id', 'inputnode.subject_id'),
                 ('outputnode.fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'),

fmriprep/workflows/bold/base.py

@@ -188,6 +188,7 @@ def init_bold_wf(
                 "t1w_preproc",
                 "t1w_mask",
                 "t1w_dseg",
+                "t1w_tpms",
                 "subjects_dir",
                 "subject_id",
                 "fsnative2t1w_xfm",
@@ -441,6 +442,50 @@ def init_bold_wf(
             (bold_std_wf, ds_bold_std_wf, [('outputnode.bold_file', 'inputnode.bold')]),
         ])  # fmt:skip
 
+    bold_confounds_wf = init_bold_confs_wf(
+        mem_gb=mem_gb["largemem"],
+        metadata=all_metadata[0],
+        freesurfer=config.workflow.run_reconall,
+        regressors_all_comps=config.workflow.regressors_all_comps,
+        regressors_fd_th=config.workflow.regressors_fd_th,
+        regressors_dvars_th=config.workflow.regressors_dvars_th,
+        name="bold_confounds_wf",
+    )
+
+    ds_confounds = pe.Node(
+        DerivativesDataSink(
+            base_directory=fmriprep_dir,
+            desc='confounds',
+            suffix='timeseries',
+            dismiss_entities=("echo",),
+        ),
+        name="ds_confounds",
+        run_without_submitting=True,
+        mem_gb=config.DEFAULT_MEMORY_MIN_GB,
+    )
+    ds_confounds.inputs.source_file = bold_file
+
+    workflow.connect([
+        (inputnode, bold_confounds_wf, [
+            ('t1w_tpms', 'inputnode.t1w_tpms'),
+            ('t1w_mask', 'inputnode.t1w_mask'),
+        ]),
+        (bold_fit_wf, bold_confounds_wf, [
+            ('outputnode.bold_mask', 'inputnode.bold_mask'),
+            ('outputnode.movpar_file', 'inputnode.movpar_file'),
+            ('outputnode.rmsd_file', 'inputnode.rmsd_file'),
+            ('outputnode.boldref2anat_xfm', 'inputnode.boldref2anat_xfm'),
+            ('outputnode.dummy_scans', 'inputnode.skip_vols'),
+        ]),
+        (bold_native_wf, bold_confounds_wf, [
+            ('outputnode.bold_native', 'inputnode.bold'),
+        ]),
+        (bold_confounds_wf, ds_confounds, [
+            ('outputnode.confounds_file', 'in_file'),
+            ('outputnode.confounds_metadata', 'meta_dict'),
+        ]),
+    ])  # fmt:skip
+
     # Fill-in datasinks of reportlets seen so far
     for node in workflow.list_node_names():
         if node.split(".")[-1].startswith("ds_report"):
@@ -677,18 +722,6 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
         name="outputnode",
     )
 
-    # get confounds
-    bold_confounds_wf = init_bold_confs_wf(
-        mem_gb=mem_gb["largemem"],
-        metadata=metadata,
-        freesurfer=freesurfer,
-        regressors_all_comps=config.workflow.regressors_all_comps,
-        regressors_fd_th=config.workflow.regressors_fd_th,
-        regressors_dvars_th=config.workflow.regressors_dvars_th,
-        name="bold_confounds_wf",
-    )
-    bold_confounds_wf.get_node("inputnode").inputs.t1_transform_flags = [False]
-
     # SURFACES ##################################################################################
     # Freesurfer
     if freesurfer and freesurfer_spaces:

fmriprep/workflows/bold/confounds.py

@@ -121,16 +121,16 @@ def init_bold_confs_wf(
     movpar_file
         SPM-formatted motion parameters file
     rmsd_file
-        Framewise displacement as measured by ``fsl_motion_outliers``.
+        Root mean squared deviation as measured by ``fsl_motion_outliers`` [Jenkinson2002]_.
     skip_vols
         number of non steady state volumes
     t1w_mask
         Mask of the skull-stripped template image
     t1w_tpms
         List of tissue probability maps in T1w space
-    t1_bold_xform
-        Affine matrix that maps the T1w space into alignment with
-        the native BOLD space
+    boldref2anat_xfm
+        Affine matrix that maps the BOLD reference space into alignment with
+        the anatomical (T1w) space
 
     Outputs
     -------
@@ -224,7 +224,7 @@ def init_bold_confs_wf(
                 "skip_vols",
                 "t1w_mask",
                 "t1w_tpms",
-                "t1_bold_xform",
+                "boldref2anat_xfm",
             ]
         ),
         name="inputnode",
@@ -244,7 +244,7 @@ def init_bold_confs_wf(
 
     # Project T1w mask into BOLD space and merge with BOLD brainmask
     t1w_mask_tfm = pe.Node(
-        ApplyTransforms(interpolation="MultiLabel"),
+        ApplyTransforms(interpolation="MultiLabel", invert_transform_flags=[True]),
         name="t1w_mask_tfm",
     )
     union_mask = pe.Node(niu.Function(function=_binary_union), name="union_mask")
@@ -266,9 +266,9 @@ def init_bold_confs_wf(
     # Generate aCompCor probseg maps
     acc_masks = pe.Node(aCompCorMasks(is_aseg=freesurfer), name="acc_masks")
 
-    # Resample probseg maps in BOLD space via T1w-to-BOLD transform
+    # Resample probseg maps in BOLD space via BOLD-to-T1w transform
     acc_msk_tfm = pe.MapNode(
-        ApplyTransforms(interpolation="Gaussian"),
+        ApplyTransforms(interpolation="Gaussian", invert_transform_flags=[True]),
         iterfield=["input_image"],
         name="acc_msk_tfm",
         mem_gb=0.1,
@@ -500,7 +500,7 @@ def _select_cols(table):
         # Brain mask
         (inputnode, t1w_mask_tfm, [("t1w_mask", "input_image"),
                                    ("bold_mask", "reference_image"),
-                                   ("t1_bold_xform", "transforms")]),
+                                   ("boldref2anat_xfm", "transforms")]),
         (inputnode, union_mask, [("bold_mask", "mask1")]),
         (t1w_mask_tfm, union_mask, [("output_image", "mask2")]),
         (union_mask, dilated_mask, [("out", "in_mask")]),
@@ -512,7 +512,7 @@ def _select_cols(table):
                                ("skip_vols", "ignore_initial_volumes")]),
         (inputnode, acc_masks, [("t1w_tpms", "in_vfs"),
                                 (("bold", _get_zooms), "bold_zooms")]),
-        (inputnode, acc_msk_tfm, [("t1_bold_xform", "transforms"),
+        (inputnode, acc_msk_tfm, [("boldref2anat_xfm", "transforms"),
                                   ("bold_mask", "reference_image")]),
         (inputnode, acc_msk_brain, [("bold_mask", "in_mask")]),
         (acc_masks, acc_msk_tfm, [("out_masks", "input_image")]),
@@ -624,9 +624,9 @@ def init_carpetplot_wf(
         BOLD series mask
     confounds_file
         TSV of all aggregated confounds
-    t1_bold_xform
-        Affine matrix that maps the T1w space into alignment with
-        the native BOLD space
+    boldref2anat_xfm
+        Affine matrix that maps the BOLD reference space into alignment with
+        the anatomical (T1w) space
     std2anat_xfm
         ANTs-compatible affine-and-warp transform file
     cifti_bold
@@ -653,7 +653,7 @@ def init_carpetplot_wf(
                 "bold",
                 "bold_mask",
                 "confounds_file",
-                "t1_bold_xform",
+                "boldref2anat_xfm",
                 "std2anat_xfm",
                 "cifti_bold",
                 "crown_mask",
@@ -708,6 +708,7 @@ def init_carpetplot_wf(
                     extension=[".nii", ".nii.gz"],
                 )
             ),
+            invert_transform_flags=[True, False],
             interpolation="MultiLabel",
             args="-u int",
         ),
@@ -720,7 +721,7 @@ def init_carpetplot_wf(
 
     # fmt:off
     workflow.connect([
-        (inputnode, mrg_xfms, [("t1_bold_xform", "in1"),
+        (inputnode, mrg_xfms, [("boldref2anat_xfm", "in1"),
                                ("std2anat_xfm", "in2")]),
         (inputnode, resample_parc, [("bold_mask", "reference_image")]),
         (inputnode, parcels, [("crown_mask", "crown_mask")]),

fmriprep/workflows/bold/fit.py

@@ -177,6 +177,12 @@ def init_bold_fit_wf(
     boldref2fmap_xfm
         Affine transform mapping from BOLD reference space to the fieldmap
         space, if applicable.
+    movpar_file
+        MCFLIRT motion parameters, normalized to SPM format (X, Y, Z, Rx, Ry, Rz)
+    rmsd_file
+        Root mean squared deviation as measured by ``fsl_motion_outliers`` [Jenkinson2002]_.
+    dummy_scans
+        The number of dummy scans declared or detected at the beginning of the series.
 
     See Also
     --------
@@ -267,6 +273,8 @@ def init_bold_fit_wf(
                 "motion_xfm",
                 "boldref2anat_xfm",
                 "boldref2fmap_xfm",
+                "movpar_file",
+                "rmsd_file",
             ],
         ),
         name="outputnode",
@@ -276,7 +284,8 @@ def init_bold_fit_wf(
     workflow.add_nodes([inputnode])
 
     hmcref_buffer = pe.Node(
-        niu.IdentityInterface(fields=["boldref", "bold_file"]), name="hmcref_buffer"
+        niu.IdentityInterface(fields=["boldref", "bold_file", "dummy_scans"]),
+        name="hmcref_buffer",
     )
     fmapref_buffer = pe.Node(niu.Function(function=_select_ref), name="fmapref_buffer")
     hmc_buffer = pe.Node(niu.IdentityInterface(fields=["hmc_xforms"]), name="hmc_buffer")
@@ -313,13 +322,20 @@ def init_bold_fit_wf(
 
     # fmt:off
     workflow.connect([
-        (hmcref_buffer, outputnode, [("boldref", "hmc_boldref")]),
+        (hmcref_buffer, outputnode, [
+            ("boldref", "hmc_boldref"),
+            ("dummy_scans", "dummy_scans"),
+        ]),
         (regref_buffer, outputnode, [
             ("boldref", "coreg_boldref"),
             ("boldmask", "bold_mask"),
         ]),
         (fmapreg_buffer, outputnode, [("boldref2fmap_xfm", "boldref2fmap_xfm")]),
-        (hmc_buffer, outputnode, [("hmc_xforms", "motion_xfm")]),
+        (hmc_buffer, outputnode, [
+            ("hmc_xforms", "motion_xfm"),
+            ("movpar_file", "movpar_file"),
+            ("rmsd_file", "rmsd_file"),
+        ]),
         (inputnode, func_fit_reports_wf, [
             ("bold_file", "inputnode.source_file"),
             ("t1w_preproc", "inputnode.t1w_preproc"),
@@ -360,6 +376,7 @@ def init_bold_fit_wf(
             (hmc_boldref_wf, hmcref_buffer, [
                 ("outputnode.bold_file", "bold_file"),
                 ("outputnode.boldref", "boldref"),
+                ("outputnode.skip_vols", "dummy_scans"),
             ]),
             (hmcref_buffer, ds_hmc_boldref_wf, [("boldref", "inputnode.boldref")]),
             (hmc_boldref_wf, summary, [("outputnode.algo_dummy_scans", "algo_dummy_scans")]),
@@ -403,7 +420,11 @@ def init_bold_fit_wf(
                 ("bold_file", "inputnode.bold_file"),
             ]),
             (bold_hmc_wf, ds_hmc_wf, [("outputnode.xforms", "inputnode.xforms")]),
-            (ds_hmc_wf, hmc_buffer, [("outputnode.xforms", "hmc_xforms")]),
+            (bold_hmc_wf, hmc_buffer, [
+                ("outputnode.xforms", "hmc_xforms"),
+                ("outputnode.movpar_file", "movpar_file"),
+                ("outputnode.rmsd_file", "rmsd_file"),
+            ]),
         ])
         # fmt:on
     else:

fmriprep/workflows/bold/hmc.py

@@ -75,8 +75,8 @@ def init_bold_hmc_wf(mem_gb: float, omp_nthreads: int, name: str = 'bold_hmc_wf'
         ITKTransform file aligning each volume to ``ref_image``
     movpar_file
         MCFLIRT motion parameters, normalized to SPM format (X, Y, Z, Rx, Ry, Rz)
-    rms_file
-        Framewise displacement as measured by ``fsl_motion_outliers`` [Jenkinson2002]_.
+    rmsd_file
+        Root mean squared deviation as measured by ``fsl_motion_outliers`` [Jenkinson2002]_.
 
     """
     from niworkflows.engine.workflows import LiterateWorkflow as Workflow

fmriprep/workflows/bold/outputs.py

@@ -961,17 +961,6 @@ def init_func_derivatives_wf(
     raw_sources = pe.Node(niu.Function(function=_bids_relative), name='raw_sources')
     raw_sources.inputs.bids_root = bids_root
 
-    ds_confounds = pe.Node(
-        DerivativesDataSink(
-            base_directory=output_dir,
-            desc='confounds',
-            suffix='timeseries',
-            dismiss_entities=("echo",),
-        ),
-        name="ds_confounds",
-        run_without_submitting=True,
-        mem_gb=DEFAULT_MEMORY_MIN_GB,
-    )
     ds_ref_t1w_xfm = pe.Node(
         DerivativesDataSink(
             base_directory=output_dir,
@@ -1001,9 +990,6 @@ def init_func_derivatives_wf(
     # fmt:off
     workflow.connect([
         (inputnode, raw_sources, [('all_source_files', 'in_files')]),
-        (inputnode, ds_confounds, [('source_file', 'source_file'),
-                                   ('confounds', 'in_file'),
-                                   ('confounds_metadata', 'meta_dict')]),
         (inputnode, ds_ref_t1w_xfm, [('source_file', 'source_file'),
                                      ('bold2anat_xfm', 'in_file')]),
         (inputnode, ds_ref_t1w_inv_xfm, [('source_file', 'source_file'),