From 9778661b59f45c5e9d3402b2066228ebc97ab143 Mon Sep 17 00:00:00 2001
From: Chris Markiewicz <effigies@gmail.com>
Date: Sat, 14 Oct 2023 21:37:09 -0400
Subject: [PATCH 01/19] ENH: Check for bold length before building workflow

---
 fmriprep/workflows/base.py      |  3 +++
 fmriprep/workflows/bold/base.py | 20 ++++++++------------
 2 files changed, 11 insertions(+), 12 deletions(-)

diff --git a/fmriprep/workflows/base.py b/fmriprep/workflows/base.py
index 2c2ad6aa4..ec901e407 100644
--- a/fmriprep/workflows/base.py
+++ b/fmriprep/workflows/base.py
@@ -483,6 +483,9 @@ def init_single_subject_wf(subject_id: str):
             precomputed=functional_cache,
             fieldmap_id=fieldmap_id,
         )
+        if bold_wf is None:
+            continue
+
         bold_wf.__desc__ = func_pre_desc + (bold_wf.__desc__ or "")
 
         workflow.connect([
diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py
index 4e4ad6e4d..fcc73686d 100644
--- a/fmriprep/workflows/bold/base.py
+++ b/fmriprep/workflows/bold/base.py
@@ -140,6 +140,13 @@ def init_bold_wf(
     fmriprep_dir = config.execution.fmriprep_dir
     omp_nthreads = config.nipype.omp_nthreads
 
+    nvols, mem_gb = _create_mem_gb(bold_file)
+    if nvols <= 5 - config.execution.sloppy:
+        config.loggers.workflow.warning(
+            f"Too short BOLD series (<= 5 timepoints). Skipping processing of <{bold_file}>."
+        )
+        return
+
     functional_cache = {}
     if config.execution.derivatives:
         from fmriprep.utils.bids import collect_derivatives, extract_entities
@@ -454,17 +461,6 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
     )
     from niworkflows.interfaces.utility import KeySelect
 
-    img = nb.load(bold_file[0] if isinstance(bold_file, (list, tuple)) else bold_file)
-    nvols = 1 if img.ndim < 4 else img.shape[3]
-    if nvols <= 5 - config.execution.sloppy:
-        config.loggers.workflow.warning(
-            f"Too short BOLD series (<= 5 timepoints). Skipping processing of <{bold_file}>."
-        )
-        return
-
-    mem_gb = {"filesize": 1, "resampled": 1, "largemem": 1}
-    bold_tlen = 10
-
     # Have some options handy
     omp_nthreads = config.nipype.omp_nthreads
     freesurfer = config.workflow.run_reconall
@@ -1439,7 +1435,7 @@ def _create_mem_gb(bold_fname):
     nvox = int(np.prod(img.shape, dtype='u8'))
     # Assume tools will coerce to 8-byte floats to be safe
     bold_size_gb = 8 * nvox / (1024**3)
-    bold_tlen = img.shape[-1]
+    bold_tlen = 1 if img.ndim < 4 else img.shape[3]
     mem_gb = {
         "filesize": bold_size_gb,
         "resampled": bold_size_gb * 4,

From c67df4aab478a40ef0ac7063771f4f054592fd7f Mon Sep 17 00:00:00 2001
From: Chris Markiewicz <effigies@gmail.com>
Date: Sat, 14 Oct 2023 21:55:08 -0400
Subject: [PATCH 02/19] RF: Delete pieces already rewritten

---
 fmriprep/workflows/bold/base.py | 581 +-------------------------------
 1 file changed, 1 insertion(+), 580 deletions(-)

diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py
index fcc73686d..5ac39db00 100644
--- a/fmriprep/workflows/bold/base.py
+++ b/fmriprep/workflows/bold/base.py
@@ -469,50 +469,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
     freesurfer_spaces = spaces.get_fs_spaces()
     project_goodvoxels = config.workflow.project_goodvoxels and config.workflow.cifti_output
 
-    # Extract BIDS entities and metadata from BOLD file(s)
-    entities = extract_entities(bold_file)
-    layout = config.execution.layout
-
-    # Extract metadata
-    all_metadata = [layout.get_metadata(fname) for fname in listify(bold_file)]
-
-    # Take first file as reference
-    ref_file = pop_file(bold_file)
-    metadata = all_metadata[0]
-    # get original image orientation
-    ref_orientation = get_img_orientation(ref_file)
-
-    echo_idxs = listify(entities.get("echo", []))
-    multiecho = len(echo_idxs) > 2
-    if len(echo_idxs) == 1:
-        config.loggers.workflow.warning(
-            f"Running a single echo <{ref_file}> from a seemingly multi-echo dataset."
-        )
-        bold_file = ref_file  # Just in case - drop the list
-
-    if len(echo_idxs) == 2:
-        raise RuntimeError(
-            "Multi-echo processing requires at least three different echos (found two)."
-        )
-
-    if multiecho:
-        # Drop echo entity for future queries, have a boolean shorthand
-        entities.pop("echo", None)
-        # reorder echoes from shortest to largest
-        tes, bold_file = zip(
-            *sorted([(layout.get_metadata(bf)["EchoTime"], bf) for bf in bold_file])
-        )
-        ref_file = bold_file[0]  # Reset reference to be the shortest TE
-        shapes = [nb.load(echo).shape for echo in bold_file]
-        if len(set(shapes)) != 1:
-            diagnostic = "\n".join(
-                f"{os.path.basename(echo)}: {shape}" for echo, shape in zip(bold_file, shapes)
-            )
-            raise RuntimeError(f"Multi-echo images found with mismatching shapes\n{diagnostic}")
-
-    if os.path.isfile(ref_file):
-        bold_tlen, mem_gb = _create_mem_gb(ref_file)
-
+    ref_file = bold_file
     wf_name = _get_wf_name(ref_file, "func_preproc")
     config.loggers.workflow.debug(
         "Creating bold processing workflow for <%s> (%.2f GB / %d TRs). "
@@ -524,50 +481,6 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
         mem_gb["largemem"],
     )
 
-    # Find associated sbref, if possible
-    overrides = {
-        "suffix": "sbref",
-        "extension": [".nii", ".nii.gz"],
-    }
-    if config.execution.bids_filters:
-        overrides.update(config.execution.bids_filters.get('sbref', {}))
-    sb_ents = {**entities, **overrides}
-    sbref_files = layout.get(return_type="file", **sb_ents)
-
-    sbref_msg = f"No single-band-reference found for {os.path.basename(ref_file)}."
-    if sbref_files and "sbref" in config.workflow.ignore:
-        sbref_msg = "Single-band reference file(s) found and ignored."
-        sbref_files = []
-    elif sbref_files:
-        sbref_msg = "Using single-band reference file(s) {}.".format(
-            ",".join([os.path.basename(sbf) for sbf in sbref_files])
-        )
-    config.loggers.workflow.info(sbref_msg)
-
-    if has_fieldmap:
-        from sdcflows import fieldmaps as fm
-
-        # We may have pruned the estimator collection due to `--ignore fieldmaps`
-        estimator_key = [
-            key
-            for key in get_estimator(layout, bold_file if not multiecho else bold_file[0])
-            if key in fm._estimators
-        ]
-
-        if not estimator_key:
-            has_fieldmap = False
-            config.loggers.workflow.critical(
-                f"None of the available B0 fieldmaps are associated to <{bold_file}>"
-            )
-        else:
-            config.loggers.workflow.info(
-                f"Found usable B0-map (fieldmap) estimator(s) <{', '.join(estimator_key)}> "
-                f"to correct <{bold_file}> for susceptibility-derived distortions."
-            )
-
-    # Check whether STC must/can be run
-    run_stc = bool(metadata.get("SliceTiming")) and "slicetiming" not in config.workflow.ignore
-
     # Build workflow
     workflow = Workflow(name=wf_name)
     workflow.__postdesc__ = """\
@@ -648,127 +561,6 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
         name="outputnode",
     )
 
-    # Outline
-    # 1) Find/create reference
-    # 2) HMC
-    # 3) Apply SDC
-    # 4) BOLD-T1w coregistration
-    # 5) T2* map
-    #
-    # Notes
-    # - STC only used for bold_split (apply)
-
-    # Generate a brain-masked conversion of the t1w
-    t1w_brain = pe.Node(ApplyMask(), name="t1w_brain")
-
-    # Track echo index - this allows us to treat multi- and single-echo workflows
-    # almost identically
-    echo_index = pe.Node(niu.IdentityInterface(fields=["echoidx"]), name="echo_index")
-    if multiecho:
-        echo_index.iterables = [("echoidx", range(len(bold_file)))]
-    else:
-        echo_index.inputs.echoidx = 0
-
-    # BOLD source: track original BOLD file(s)
-    bold_source = pe.Node(niu.Select(inlist=bold_file), name="bold_source")
-
-    # BOLD buffer: an identity used as a pointer to either the original BOLD
-    # or the STC'ed one for further use.
-    boldbuffer = pe.Node(niu.IdentityInterface(fields=["bold_file"]), name="boldbuffer")
-
-    summary = pe.Node(
-        FunctionalSummary(
-            # slice_timing=run_stc,
-            registration=("FSL", "FreeSurfer")[freesurfer],
-            registration_dof=config.workflow.bold2t1w_dof,
-            registration_init=config.workflow.bold2t1w_init,
-            pe_direction=metadata.get("PhaseEncodingDirection"),
-            echo_idx=echo_idxs,
-            tr=metadata["RepetitionTime"],
-            orientation=ref_orientation,
-        ),
-        name="summary",
-        mem_gb=config.DEFAULT_MEMORY_MIN_GB,
-        run_without_submitting=True,
-    )
-    summary.inputs.dummy_scans = config.workflow.dummy_scans
-
-    func_derivatives_wf = init_func_derivatives_wf(
-        bids_root=layout.root,
-        cifti_output=config.workflow.cifti_output,
-        freesurfer=freesurfer,
-        project_goodvoxels=project_goodvoxels,
-        all_metadata=all_metadata,
-        multiecho=multiecho,
-        output_dir=fmriprep_dir,
-        spaces=spaces,
-    )
-    func_derivatives_wf.inputs.inputnode.all_source_files = bold_file
-    func_derivatives_wf.inputs.inputnode.cifti_density = config.workflow.cifti_output
-
-    # fmt:off
-    workflow.connect([
-        (outputnode, func_derivatives_wf, [
-            ("bold_t1", "inputnode.bold_t1"),
-            ("bold_t1_ref", "inputnode.bold_t1_ref"),
-            ("bold2anat_xfm", "inputnode.bold2anat_xfm"),
-            ("anat2bold_xfm", "inputnode.anat2bold_xfm"),
-            ("hmc_xforms", "inputnode.hmc_xforms"),
-            ("bold_aseg_t1", "inputnode.bold_aseg_t1"),
-            ("bold_aparc_t1", "inputnode.bold_aparc_t1"),
-            ("bold_mask_t1", "inputnode.bold_mask_t1"),
-            ("bold_native", "inputnode.bold_native"),
-            ("bold_native_ref", "inputnode.bold_native_ref"),
-            ("bold_mask_native", "inputnode.bold_mask_native"),
-            ("bold_echos_native", "inputnode.bold_echos_native"),
-            ("confounds", "inputnode.confounds"),
-            ("surfaces", "inputnode.surf_files"),
-            ("bold_cifti", "inputnode.bold_cifti"),
-            ("cifti_metadata", "inputnode.cifti_metadata"),
-            ("t2star_bold", "inputnode.t2star_bold"),
-            ("t2star_t1", "inputnode.t2star_t1"),
-            ("t2star_std", "inputnode.t2star_std"),
-            ("confounds_metadata", "inputnode.confounds_metadata"),
-            ("acompcor_masks", "inputnode.acompcor_masks"),
-            ("tcompcor_mask", "inputnode.tcompcor_mask"),
-        ]),
-    ])
-    # fmt:on
-
-    # Generate a tentative boldref
-    initial_boldref_wf = init_bold_reference_wf(
-        name="initial_boldref_wf",
-        omp_nthreads=omp_nthreads,
-        bold_file=bold_file,
-        sbref_files=sbref_files,
-        multiecho=multiecho,
-    )
-    initial_boldref_wf.inputs.inputnode.dummy_scans = config.workflow.dummy_scans
-
-    # Select validated BOLD files (orientations checked or corrected)
-    select_bold = pe.Node(niu.Select(), name="select_bold")
-
-    # Top-level BOLD splitter
-    bold_split = pe.Node(FSLSplit(dimension="t"), name="bold_split", mem_gb=mem_gb["filesize"] * 3)
-
-    # HMC on the BOLD
-    bold_hmc_wf = init_bold_hmc_wf(
-        name="bold_hmc_wf", mem_gb=mem_gb["filesize"], omp_nthreads=omp_nthreads
-    )
-
-    # calculate BOLD registration to T1w
-    bold_reg_wf = init_bold_reg_wf(
-        bold2t1w_dof=config.workflow.bold2t1w_dof,
-        bold2t1w_init=config.workflow.bold2t1w_init,
-        freesurfer=freesurfer,
-        mem_gb=mem_gb["resampled"],
-        name="bold_reg_wf",
-        omp_nthreads=omp_nthreads,
-        sloppy=config.execution.sloppy,
-        use_bbr=config.workflow.use_bbr,
-        use_compression=False,
-    )
-
     # apply BOLD registration to T1w
     bold_t1_trans_wf = init_bold_t1_trans_wf(
         name="bold_t1_trans_wf",
@@ -791,114 +583,6 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
     )
     bold_confounds_wf.get_node("inputnode").inputs.t1_transform_flags = [False]
 
-    # SLICE-TIME CORRECTION (or bypass) #############################################
-    if run_stc:
-        bold_stc_wf = init_bold_stc_wf(name="bold_stc_wf", metadata=metadata)
-        # fmt:off
-        workflow.connect([
-            (initial_boldref_wf, bold_stc_wf, [("outputnode.skip_vols", "inputnode.skip_vols")]),
-            (select_bold, bold_stc_wf, [("out", "inputnode.bold_file")]),
-            (bold_stc_wf, boldbuffer, [("outputnode.stc_file", "bold_file")]),
-        ])
-        # fmt:on
-
-    # bypass STC from original BOLD in both SE and ME cases
-    else:
-        workflow.connect([(select_bold, boldbuffer, [("out", "bold_file")])])
-
-    # MULTI-ECHO EPI DATA #############################################
-    if multiecho:  # instantiate relevant interfaces, imports
-        split_opt_comb = bold_split.clone(name="split_opt_comb")
-
-        inputnode.inputs.bold_file = ref_file  # Replace reference w first echo
-
-        join_echos = pe.JoinNode(
-            niu.IdentityInterface(fields=["bold_files"]),
-            joinsource="echo_index",
-            joinfield=["bold_files"],
-            name="join_echos",
-        )
-
-        # create optimal combination, adaptive T2* map
-        bold_t2s_wf = init_bold_t2s_wf(
-            echo_times=tes,
-            mem_gb=mem_gb["filesize"],
-            omp_nthreads=omp_nthreads,
-            name="bold_t2smap_wf",
-        )
-
-        t2s_reporting_wf = init_t2s_reporting_wf()
-
-        ds_report_t2scomp = pe.Node(
-            DerivativesDataSink(
-                desc="t2scomp",
-                datatype="figures",
-                dismiss_entities=("echo",),
-            ),
-            name="ds_report_t2scomp",
-            run_without_submitting=True,
-        )
-
-        ds_report_t2star_hist = pe.Node(
-            DerivativesDataSink(
-                desc="t2starhist",
-                datatype="figures",
-                dismiss_entities=("echo",),
-            ),
-            name="ds_report_t2star_hist",
-            run_without_submitting=True,
-        )
-
-    bold_final = pe.Node(
-        niu.IdentityInterface(fields=["bold", "boldref", "mask", "bold_echos", "t2star"]),
-        name="bold_final",
-    )
-
-    # Generate a final BOLD reference
-    # This BOLD references *does not use* single-band reference images.
-    final_boldref_wf = init_bold_reference_wf(
-        name="final_boldref_wf",
-        omp_nthreads=omp_nthreads,
-        multiecho=multiecho,
-    )
-    final_boldref_wf.__desc__ = None  # Unset description to avoid second appearance
-
-    # for standard EPI data, pass along correct file
-    if not multiecho:
-        # fmt:off
-        workflow.connect([
-            (inputnode, func_derivatives_wf, [("bold_file", "inputnode.source_file")]),
-            (bold_split, bold_t1_trans_wf, [("out_files", "inputnode.bold_split")]),
-            (bold_hmc_wf, bold_t1_trans_wf, [("outputnode.xforms", "inputnode.hmc_xforms")]),
-        ])
-        # fmt:on
-    else:  # for meepi, use optimal combination
-        # fmt:off
-        workflow.connect([
-            # update name source for optimal combination
-            (inputnode, func_derivatives_wf, [
-                (("bold_file", combine_meepi_source), "inputnode.source_file"),
-            ]),
-            (join_echos, bold_t2s_wf, [("bold_files", "inputnode.bold_file")]),
-            (join_echos, bold_final, [("bold_files", "bold_echos")]),
-            (bold_t2s_wf, split_opt_comb, [("outputnode.bold", "in_file")]),
-            (split_opt_comb, bold_t1_trans_wf, [("out_files", "inputnode.bold_split")]),
-            (bold_t2s_wf, bold_final, [("outputnode.bold", "bold"),
-                                       ("outputnode.t2star_map", "t2star")]),
-            (inputnode, t2s_reporting_wf, [("t1w_dseg", "inputnode.label_file")]),
-            (bold_reg_wf, t2s_reporting_wf, [
-                ("outputnode.itk_t1_to_bold", "inputnode.label_bold_xform")
-            ]),
-            (bold_final, t2s_reporting_wf, [("t2star", "inputnode.t2star_file"),
-                                            ("boldref", "inputnode.boldref")]),
-            (t2s_reporting_wf, ds_report_t2scomp, [('outputnode.t2s_comp_report', 'in_file')]),
-            (t2s_reporting_wf, ds_report_t2star_hist, [("outputnode.t2star_hist", "in_file")]),
-        ])
-        # fmt:on
-
-        # Already applied in bold_bold_trans_wf, which inputs to bold_t2s_wf
-        bold_t1_trans_wf.inputs.inputnode.hmc_xforms = "identity"
-
     # Map final BOLD mask into T1w space (if required)
     nonstd_spaces = set(spaces.get_nonstandard())
     if nonstd_spaces.intersection(("T1w", "anat")):
@@ -1164,269 +848,6 @@ def _last(inlist):
     ])
     # fmt:on
 
-    if not has_fieldmap:
-        # Finalize workflow without SDC connections
-        summary.inputs.distortion_correction = "None"
-
-        # Resample in native space in just one shot
-        bold_bold_trans_wf = init_bold_preproc_trans_wf(
-            mem_gb=mem_gb["resampled"],
-            omp_nthreads=omp_nthreads,
-            use_compression=not config.execution.low_mem,
-            use_fieldwarp=False,
-            name="bold_bold_trans_wf",
-        )
-        bold_bold_trans_wf.inputs.inputnode.fieldwarp = "identity"
-
-        # fmt:off
-        workflow.connect([
-            # Connect bold_bold_trans_wf
-            (bold_source, bold_bold_trans_wf, [("out", "inputnode.name_source")]),
-            (bold_split, bold_bold_trans_wf, [("out_files", "inputnode.bold_file")]),
-            (bold_hmc_wf, bold_bold_trans_wf, [
-                ("outputnode.xforms", "inputnode.hmc_xforms"),
-            ]),
-        ])
-
-        workflow.connect([
-            (bold_bold_trans_wf, bold_final, [("outputnode.bold", "bold")]),
-            (bold_bold_trans_wf, final_boldref_wf, [
-                ("outputnode.bold", "inputnode.bold_file"),
-            ]),
-        ] if not multiecho else [
-            (initial_boldref_wf, bold_t2s_wf, [
-                ("outputnode.bold_mask", "inputnode.bold_mask"),
-            ]),
-            (bold_bold_trans_wf, join_echos, [
-                ("outputnode.bold", "bold_files"),
-            ]),
-            (join_echos, final_boldref_wf, [
-                ("bold_files", "inputnode.bold_file"),
-            ]),
-        ])
-        # fmt:on
-        return workflow
-
-    from niworkflows.interfaces.utility import KeySelect
-    from sdcflows.workflows.apply.correction import init_unwarp_wf
-    from sdcflows.workflows.apply.registration import init_coeff2epi_wf
-
-    coeff2epi_wf = init_coeff2epi_wf(
-        debug="fieldmaps" in config.execution.debug,
-        omp_nthreads=config.nipype.omp_nthreads,
-        sloppy=config.execution.sloppy,
-        write_coeff=True,
-    )
-    unwarp_wf = init_unwarp_wf(
-        free_mem=config.environment.free_mem,
-        debug="fieldmaps" in config.execution.debug,
-        omp_nthreads=config.nipype.omp_nthreads,
-    )
-    unwarp_wf.inputs.inputnode.metadata = metadata
-
-    output_select = pe.Node(
-        KeySelect(fields=["fmap", "fmap_ref", "fmap_coeff", "fmap_mask", "sdc_method"]),
-        name="output_select",
-        run_without_submitting=True,
-    )
-    output_select.inputs.key = estimator_key[0]
-    if len(estimator_key) > 1:
-        config.loggers.workflow.warning(
-            f"Several fieldmaps <{', '.join(estimator_key)}> are "
-            f"'IntendedFor' <{bold_file}>, using {estimator_key[0]}"
-        )
-
-    sdc_report = pe.Node(
-        SimpleBeforeAfter(
-            before_label="Distorted",
-            after_label="Corrected",
-            dismiss_affine=True,
-        ),
-        name="sdc_report",
-        mem_gb=0.1,
-    )
-
-    ds_report_sdc = pe.Node(
-        DerivativesDataSink(
-            base_directory=fmriprep_dir,
-            desc="sdc",
-            suffix="bold",
-            datatype="figures",
-            dismiss_entities=("echo",),
-        ),
-        name="ds_report_sdc",
-        run_without_submitting=True,
-    )
-
-    # fmt:off
-    workflow.connect([
-        (inputnode, output_select, [("fmap", "fmap"),
-                                    ("fmap_ref", "fmap_ref"),
-                                    ("fmap_coeff", "fmap_coeff"),
-                                    ("fmap_mask", "fmap_mask"),
-                                    ("sdc_method", "sdc_method"),
-                                    ("fmap_id", "keys")]),
-        (output_select, coeff2epi_wf, [
-            ("fmap_ref", "inputnode.fmap_ref"),
-            ("fmap_coeff", "inputnode.fmap_coeff"),
-            ("fmap_mask", "inputnode.fmap_mask")]),
-        (output_select, summary, [("sdc_method", "distortion_correction")]),
-        (initial_boldref_wf, coeff2epi_wf, [
-            ("outputnode.ref_image", "inputnode.target_ref"),
-            ("outputnode.bold_mask", "inputnode.target_mask")]),
-        (initial_boldref_wf, unwarp_wf, [
-            ("outputnode.ref_image", "inputnode.distorted_ref"),
-        ]),
-        (coeff2epi_wf, unwarp_wf, [
-            ("outputnode.fmap_coeff", "inputnode.fmap_coeff")]),
-        (bold_hmc_wf, unwarp_wf, [
-            ("outputnode.xforms", "inputnode.hmc_xforms")]),
-        (initial_boldref_wf, sdc_report, [
-            ("outputnode.ref_image", "before")]),
-        (bold_split, unwarp_wf, [
-            ("out_files", "inputnode.distorted")]),
-        (final_boldref_wf, sdc_report, [
-            ("outputnode.ref_image", "after"),
-            ("outputnode.bold_mask", "wm_seg")]),
-        (inputnode, ds_report_sdc, [("bold_file", "source_file")]),
-        (sdc_report, ds_report_sdc, [("out_report", "in_file")]),
-
-    ])
-    # fmt:on
-
-    if "fieldmaps" in config.execution.debug:
-        # Generate additional reportlets to assess SDC
-        from sdcflows.interfaces.reportlets import FieldmapReportlet
-
-        # First, one for checking the co-registration between fieldmap and EPI
-        sdc_coreg_report = pe.Node(
-            SimpleBeforeAfter(
-                before_label="Distorted target",
-                after_label="Fieldmap ref.",
-            ),
-            name="sdc_coreg_report",
-            mem_gb=0.1,
-        )
-        ds_report_sdc_coreg = pe.Node(
-            DerivativesDataSink(
-                base_directory=fmriprep_dir,
-                datatype="figures",
-                desc="fmapCoreg",
-                dismiss_entities=("echo",),
-                suffix="bold",
-            ),
-            name="ds_report_sdc_coreg",
-            run_without_submitting=True,
-        )
-
-        # Second, showing the fieldmap reconstructed from coefficients in the EPI space
-        fmap_report = pe.Node(FieldmapReportlet(), "fmap_report")
-
-        ds_fmap_report = pe.Node(
-            DerivativesDataSink(
-                base_directory=fmriprep_dir,
-                datatype="figures",
-                desc="fieldmap",
-                dismiss_entities=("echo",),
-                suffix="bold",
-            ),
-            name="ds_fmap_report",
-            run_without_submitting=True,
-        )
-
-        # fmt:off
-        workflow.connect([
-            (initial_boldref_wf, sdc_coreg_report, [
-                ("outputnode.ref_image", "before"),
-            ]),
-            (coeff2epi_wf, sdc_coreg_report, [
-                ("coregister.inverse_warped_image", "after"),
-            ]),
-            (final_boldref_wf, sdc_coreg_report, [
-                ("outputnode.bold_mask", "wm_seg"),
-            ]),
-            (inputnode, ds_report_sdc_coreg, [("bold_file", "source_file")]),
-            (sdc_coreg_report, ds_report_sdc_coreg, [("out_report", "in_file")]),
-            (unwarp_wf, fmap_report, [(("outputnode.fieldmap", pop_file), "fieldmap")]),
-            (coeff2epi_wf, fmap_report, [
-                ("coregister.inverse_warped_image", "reference"),
-            ]),
-            (final_boldref_wf, fmap_report, [
-                ("outputnode.bold_mask", "mask"),
-            ]),
-
-            (fmap_report, ds_fmap_report, [("out_report", "in_file")]),
-            (inputnode, ds_fmap_report, [("bold_file", "source_file")]),
-        ])
-        # fmt:on
-
-    if not multiecho:
-        # fmt:off
-        workflow.connect([
-            (unwarp_wf, bold_final, [("outputnode.corrected", "bold")]),
-            # remaining workflow connections
-            (unwarp_wf, final_boldref_wf, [
-                ("outputnode.corrected", "inputnode.bold_file"),
-            ]),
-            (unwarp_wf, bold_t1_trans_wf, [
-                # TEMPORARY: For the moment we can't use frame-wise fieldmaps
-                (("outputnode.fieldwarp_ref", pop_file), "inputnode.fieldwarp"),
-            ]),
-            (unwarp_wf, bold_std_trans_wf, [
-                # TEMPORARY: For the moment we can't use frame-wise fieldmaps
-                (("outputnode.fieldwarp_ref", pop_file), "inputnode.fieldwarp"),
-            ]),
-        ])
-        # fmt:on
-        return workflow
-
-    # Finalize connections if ME-EPI
-    join_sdc_echos = pe.JoinNode(
-        niu.IdentityInterface(
-            fields=[
-                "fieldmap",
-                "fieldwarp",
-                "corrected",
-                "corrected_ref",
-                "corrected_mask",
-            ]
-        ),
-        joinsource="echo_index",
-        joinfield=[
-            "fieldmap",
-            "fieldwarp",
-            "corrected",
-            "corrected_ref",
-            "corrected_mask",
-        ],
-        name="join_sdc_echos",
-    )
-
-    def _dpop(list_of_lists):
-        return list_of_lists[0][0]
-
-    # fmt:off
-    workflow.connect([
-        (unwarp_wf, join_echos, [
-            ("outputnode.corrected", "bold_files"),
-        ]),
-        (unwarp_wf, join_sdc_echos, [
-            ("outputnode.fieldmap", "fieldmap"),
-            ("outputnode.fieldwarp", "fieldwarp"),
-            ("outputnode.corrected", "corrected"),
-            ("outputnode.corrected_ref", "corrected_ref"),
-            ("outputnode.corrected_mask", "corrected_mask"),
-        ]),
-        # remaining workflow connections
-        (join_sdc_echos, final_boldref_wf, [
-            ("corrected", "inputnode.bold_file"),
-        ]),
-        (join_sdc_echos, bold_t2s_wf, [
-            (("corrected_mask", pop_file), "inputnode.bold_mask"),
-        ]),
-    ])
-    # fmt:on
-
     return workflow
 
 

From c4799a622ba946cd0c68b9e75a45125f125c9038 Mon Sep 17 00:00:00 2001
From: Chris Markiewicz <effigies@gmail.com>
Date: Mon, 16 Oct 2023 14:38:47 -0400
Subject: [PATCH 03/19] ENH: Add generic resampling workflow, use for
 anatomical

---
 fmriprep/workflows/bold/apply.py | 103 +++++++++++++++++++++++++++++++
 fmriprep/workflows/bold/base.py  |  44 +++++++++----
 fmriprep/workflows/bold/fit.py   |  32 +++++-----
 3 files changed, 150 insertions(+), 29 deletions(-)

diff --git a/fmriprep/workflows/bold/apply.py b/fmriprep/workflows/bold/apply.py
index 449f81ae1..8b22303eb 100644
--- a/fmriprep/workflows/bold/apply.py
+++ b/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))
diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py
index 5ac39db00..c962663ba 100644
--- a/fmriprep/workflows/bold/base.py
+++ b/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()
 
diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py
index 36ec19861..3855d4256 100644
--- a/fmriprep/workflows/bold/fit.py
+++ b/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,10 +734,11 @@ 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
@@ -746,6 +746,7 @@ def init_bold_native_wf(
         ),
         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
 
 

From 9864b53c8692b6e78bc9f1ea546e7ddcf7740d2f Mon Sep 17 00:00:00 2001
From: Chris Markiewicz <effigies@gmail.com>
Date: Tue, 17 Oct 2023 10:32:01 -0400
Subject: [PATCH 04/19] ENH: Add interpolation parameters to resampling
 interface

---
 fmriprep/interfaces/resampling.py | 37 ++++++++++++++++++++++++++++++-
 1 file changed, 36 insertions(+), 1 deletion(-)

diff --git a/fmriprep/interfaces/resampling.py b/fmriprep/interfaces/resampling.py
index 8cfcd2348..50b640abf 100644
--- a/fmriprep/interfaces/resampling.py
+++ b/fmriprep/interfaces/resampling.py
@@ -48,6 +48,16 @@ class ResampleSeriesInputSpec(TraitedSpec):
         desc="the phase-encoding direction corresponding to in_data",
     )
     num_threads = traits.Int(1, usedefault=True, desc="Number of threads to use for resampling")
+    output_data_type = traits.Str("float32", usedefault=True, desc="Data type of output image")
+    order = traits.Int(3, usedefault=True, desc="Order of interpolation (0=nearest, 3=cubic)")
+    mode = traits.Str(
+        'constant',
+        usedefault=True,
+        desc="How data is extended beyond its boundaries. "
+        "See scipy.ndimage.map_coordinates for more details.",
+    )
+    cval = traits.Float(0.0, usedefault=True, desc="Value to fill past edges of data")
+    prefilter = traits.Bool(True, usedefault=True, desc="Spline-prefilter data if order > 1")
 
 
 class ResampleSeriesOutputSpec(TraitedSpec):
@@ -94,6 +104,11 @@ def _run_interface(self, runtime):
             fieldmap=fieldmap,
             pe_info=pe_info,
             nthreads=self.inputs.num_threads,
+            output_dtype=self.inputs.output_data_type,
+            order=self.inputs.order,
+            mode=self.inputs.mode,
+            cval=self.inputs.cval,
+            prefilter=self.inputs.prefilter,
         )
         resampled.to_filename(out_path)
 
@@ -458,6 +473,11 @@ def resample_bold(
     fieldmap: nb.Nifti1Image | None,
     pe_info: list[tuple[int, float]] | None,
     nthreads: int = 1,
+    output_dtype: np.dtype | str | None = 'f4',
+    order: int = 3,
+    mode: str = 'constant',
+    cval: float = 0.0,
+    prefilter: bool = True,
 ) -> nb.Nifti1Image:
     """Resample a 4D bold series into a target space, applying head-motion
     and susceptibility-distortion correction simultaneously.
@@ -480,6 +500,17 @@ def resample_bold(
         of the data array in the second dimension.
     nthreads
         Number of threads to use for parallel resampling
+    output_dtype
+        The dtype of the output array.
+    order
+        Order of interpolation (default: 3 = cubic)
+    mode
+        How ``data`` is extended beyond its boundaries. See
+        :func:`scipy.ndimage.map_coordinates` for more details.
+    cval
+        Value to fill past edges of ``data`` if ``mode`` is ``'constant'``.
+    prefilter
+        Determines if ``data`` is pre-filtered before interpolation.
 
     Returns
     -------
@@ -527,8 +558,12 @@ def resample_bold(
         pe_info=pe_info,
         hmc_xfms=hmc_xfms,
         fmap_hz=fieldmap.get_fdata(dtype='f4'),
-        output_dtype='f4',
+        output_dtype=output_dtype,
         nthreads=nthreads,
+        order=order,
+        mode=mode,
+        cval=cval,
+        prefilter=prefilter,
     )
     resampled_img = nb.Nifti1Image(resampled_data, target.affine, target.header)
     resampled_img.set_data_dtype('f4')

From 8199d4633244615332369181f9ed5dfc5c266468 Mon Sep 17 00:00:00 2001
From: Chris Markiewicz <effigies@gmail.com>
Date: Tue, 17 Oct 2023 11:00:18 -0400
Subject: [PATCH 05/19] STY: Black

---
 fmriprep/interfaces/resampling.py | 16 ++++++++++++----
 1 file changed, 12 insertions(+), 4 deletions(-)

diff --git a/fmriprep/interfaces/resampling.py b/fmriprep/interfaces/resampling.py
index 50b640abf..82e78c9ed 100644
--- a/fmriprep/interfaces/resampling.py
+++ b/fmriprep/interfaces/resampling.py
@@ -28,7 +28,9 @@ class ResampleSeriesInputSpec(TraitedSpec):
     in_file = File(exists=True, mandatory=True, desc="3D or 4D image file to resample")
     ref_file = File(exists=True, mandatory=True, desc="File to resample in_file to")
     transforms = InputMultiObject(
-        File(exists=True), mandatory=True, desc="Transform files, from in_file to ref_file (image mode)"
+        File(exists=True),
+        mandatory=True,
+        desc="Transform files, from in_file to ref_file (image mode)",
     )
     inverse = InputMultiObject(
         traits.Bool,
@@ -120,10 +122,16 @@ class ReconstructFieldmapInputSpec(TraitedSpec):
     in_coeffs = InputMultiObject(
         File(exists=True), mandatory=True, desc="SDCflows-style spline coefficient files"
     )
-    target_ref_file = File(exists=True, mandatory=True, desc="Image to reconstruct the field in alignment with")
-    fmap_ref_file = File(exists=True, mandatory=True, desc="Reference file aligned with coefficients")
+    target_ref_file = File(
+        exists=True, mandatory=True, desc="Image to reconstruct the field in alignment with"
+    )
+    fmap_ref_file = File(
+        exists=True, mandatory=True, desc="Reference file aligned with coefficients"
+    )
     transforms = InputMultiObject(
-        File(exists=True), mandatory=True, desc="Transform files, from in_file to ref_file (image mode)"
+        File(exists=True),
+        mandatory=True,
+        desc="Transform files, from in_file to ref_file (image mode)",
     )
     inverse = InputMultiObject(
         traits.Bool,

From 41bc5f2be005767975509503a213940e267eeb36 Mon Sep 17 00:00:00 2001
From: Chris Markiewicz <effigies@gmail.com>
Date: Tue, 17 Oct 2023 11:00:39 -0400
Subject: [PATCH 06/19] RF: Rename functions to match usage

---
 fmriprep/interfaces/resampling.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/fmriprep/interfaces/resampling.py b/fmriprep/interfaces/resampling.py
index 82e78c9ed..9921d6f19 100644
--- a/fmriprep/interfaces/resampling.py
+++ b/fmriprep/interfaces/resampling.py
@@ -99,7 +99,7 @@ def _run_interface(self, runtime):
 
             pe_info = [(pe_axis, -ro_time if (axis_flip ^ pe_flip) else ro_time)] * nvols
 
-        resampled = resample_bold(
+        resampled = resample_image(
             source=source,
             target=target,
             transforms=transforms,
@@ -474,7 +474,7 @@ def resample_series(
     )
 
 
-def resample_bold(
+def resample_image(
     source: nb.Nifti1Image,
     target: nb.Nifti1Image,
     transforms: nt.TransformChain,
@@ -487,13 +487,13 @@ def resample_bold(
     cval: float = 0.0,
     prefilter: bool = True,
 ) -> nb.Nifti1Image:
-    """Resample a 4D bold series into a target space, applying head-motion
+    """Resample a 3- or 4D image into a target space, applying head-motion
     and susceptibility-distortion correction simultaneously.
 
     Parameters
     ----------
     source
-        The 4D bold series to resample.
+        The 3D bold image or 4D bold series to resample.
     target
         An image sampled in the target space.
     transforms

From ded0a73f9e0181f797efcd44176970cff5b80d24 Mon Sep 17 00:00:00 2001
From: Chris Markiewicz <effigies@gmail.com>
Date: Tue, 17 Oct 2023 14:43:28 -0400
Subject: [PATCH 07/19] ENH: Add general volumetric datasink workflow

---
 fmriprep/workflows/bold/outputs.py | 165 +++++++++++++++++++++++++++++
 1 file changed, 165 insertions(+)

diff --git a/fmriprep/workflows/bold/outputs.py b/fmriprep/workflows/bold/outputs.py
index b893abf50..8c593b42f 100644
--- a/fmriprep/workflows/bold/outputs.py
+++ b/fmriprep/workflows/bold/outputs.py
@@ -683,6 +683,171 @@ def init_ds_bold_native_wf(
     return workflow
 
 
+def init_ds_volumes_wf(
+    *,
+    bids_root: str,
+    output_dir: str,
+    space: str,
+    multiecho: bool,
+    metadata: ty.List[dict],
+    name="ds_bold_native_wf",
+) -> pe.Workflow:
+    metadata = all_metadata[0]
+    timing_parameters = prepare_timing_parameters(metadata)
+
+    workflow = pe.Workflow(name=name)
+    inputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=[
+                'source_files',
+                'ref_file',
+                'bold',  # Resampled into target space
+                'bold_mask',  # boldref space
+                'bold_ref',  # boldref space
+                't2star',  # boldref space
+                # Anatomical
+                'boldref2anat_xfm',
+                # Template
+                'anat2std_xfm',
+                # Entities
+                'space',
+                'cohort',
+                'resolution',
+            ]
+        ),
+        name='inputnode',
+    )
+
+    raw_sources = pe.Node(niu.Function(function=_bids_relative), name='raw_sources')
+    raw_sources.inputs.bids_root = bids_root
+    boldref2target = pe.Node(niu.Merge(2), name='boldref2target')
+
+    # BOLD is pre-resampled
+    ds_bold = pe.Node(
+        DerivativesDataSink(
+            base_directory=output_dir,
+            desc='preproc',
+            compress=True,
+            SkullStripped=multiecho,
+            TaskName=metadata.get('TaskName'),
+            dismiss_entities=("echo",),
+            **timing_parameters,
+        ),
+        name='ds_bold',
+        run_without_submitting=True,
+        mem_gb=DEFAULT_MEMORY_MIN_GB,
+    )
+    workflow.connect([
+        (inputnode, raw_sources, [('source_files', 'in_files')]),
+        (inputnode, boldref2target, [
+            ('boldref2anat_xfm', 'in1'),
+            ('anat2std_xfm', 'in2'),
+        ]),
+        (inputnode, ds_bold, [
+            ('source_files', 'source_file'),
+            ('bold', 'in_file'),
+            ('space', 'space'),
+            ('cohort', 'cohort'),
+            ('resolution', 'resolution'),
+        ]),
+    ])  # fmt:skip
+
+    resample_ref = pe.Node(
+        ApplyTransforms(
+            dimension=3,
+            default_value=0,
+            float=True,
+            interpolation="LanczosWindowedSinc",
+        ),
+        name="resample_ref",
+    )
+    resample_mask = pe.Node(ApplyTransforms(interpolation="MultiLabel"), name="resample_mask")
+    resamplers = [resample_ref, resample_mask]
+
+    workflow.connect([
+        (inputnode, resample_ref, [('bold_ref', 'input_image')]),
+        (inputnode, resample_mask, [('bold_mask', 'input_image')]),
+    ])  # fmt:skip
+
+    ds_ref = pe.Node(
+        DerivativesDataSink(
+            base_directory=output_dir,
+            suffix='boldref',
+            compress=True,
+            dismiss_entities=("echo",),
+        ),
+        name='ds_bold',
+        run_without_submitting=True,
+        mem_gb=DEFAULT_MEMORY_MIN_GB,
+    )
+    ds_mask = pe.Node(
+        DerivativesDataSink(
+            base_directory=output_dir,
+            desc='brain',
+            suffix='mask',
+            compress=True,
+            dismiss_entities=("echo",),
+        ),
+        name='ds_bold',
+        run_without_submitting=True,
+        mem_gb=DEFAULT_MEMORY_MIN_GB,
+    )
+    datasinks = [ds_ref, ds_mask]
+
+    if multiecho:
+        resample_t2star = pe.Node(
+            ApplyTransforms(
+                dimension=3,
+                default_value=0,
+                float=True,
+                interpolation="LanczosWindowedSinc",
+            ),
+            name="resample_ref",
+        )
+        ds_t2star = pe.Node(
+            DerivativesDataSink(
+                base_directory=output_dir,
+                suffix='T2starmap',
+                compress=True,
+                dismiss_entities=("echo",),
+                **t2star_meta,
+            ),
+            name='ds_t2star_std',
+            run_without_submitting=True,
+            mem_gb=DEFAULT_MEMORY_MIN_GB,
+        )
+        resamplers.append(resample_t2star)
+        datasinks.append(ds_t2star)
+
+        workflow.connect([(inputnode, resample_t2star, [('t2star', 'input_image')])])
+
+    workflow.connect(
+        [
+            (inputnode, resampler, [('ref_file', 'reference_image')])
+            for resampler in resamplers
+        ]
+        + [
+            (boldref2target, resampler, [('out', 'transforms')])
+            for resampler in resamplers
+        ]
+        + [
+            (inputnode, datasink, [
+                ('source_files', 'source_file'),
+                ('space', 'space'),
+                ('cohort', 'cohort'),
+                ('resolution', 'resolution'),
+            ])
+            for datasink in datasinks
+        ]
+        + [
+            (resampler, datasink, [("output_image", "in_file")])
+            for resampler, datasink in zip(resamplers, datasinks)
+        ]
+    )  # fmt:skip
+
+    return workflow
+
+
 def init_func_derivatives_wf(
     bids_root: str,
     cifti_output: bool,

From c094f0094e416d4b77f877ac410d852307e374dd Mon Sep 17 00:00:00 2001
From: Chris Markiewicz <effigies@gmail.com>
Date: Wed, 18 Oct 2023 13:49:58 -0400
Subject: [PATCH 08/19] ENH: Hook up anatomical derivatives workflow

---
 fmriprep/workflows/bold/apply.py   |  3 +++
 fmriprep/workflows/bold/base.py    | 31 +++++++++++++++++++++++++++---
 fmriprep/workflows/bold/outputs.py |  1 -
 3 files changed, 31 insertions(+), 4 deletions(-)

diff --git a/fmriprep/workflows/bold/apply.py b/fmriprep/workflows/bold/apply.py
index 8b22303eb..9aff94921 100644
--- a/fmriprep/workflows/bold/apply.py
+++ b/fmriprep/workflows/bold/apply.py
@@ -54,6 +54,8 @@ def init_bold_volumetric_resample_wf(
         name='inputnode',
     )
 
+    outputnode = pe.Node(niu.IdentityInterface(fields=["bold_file"]), name='outputnode')
+
     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)
@@ -70,6 +72,7 @@ def init_bold_volumetric_resample_wf(
         ]),
         (boldref2target, bold2target, [('out', 'in2')]),
         (bold2target, resample, [('out', 'transforms')]),
+        (resample, outputnode, [('out_file', 'bold_file')]),
     ])  # fmt:skip
 
     if fieldmap_id:
diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py
index c962663ba..a849cbc22 100644
--- a/fmriprep/workflows/bold/base.py
+++ b/fmriprep/workflows/bold/base.py
@@ -49,7 +49,11 @@
 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
-from .outputs import init_ds_bold_native_wf, init_func_derivatives_wf
+from .outputs import (
+    init_ds_bold_native_wf,
+    init_ds_volumes_wf,
+    init_func_derivatives_wf,
+)
 from .registration import init_bold_reg_wf, init_bold_t1_trans_wf
 from .resampling import (
     init_bold_preproc_trans_wf,
@@ -140,6 +144,7 @@ def init_bold_wf(
 
     fmriprep_dir = config.execution.fmriprep_dir
     omp_nthreads = config.nipype.omp_nthreads
+    all_metadata = [config.execution.layout.get_metadata(file) for file in bold_series]
 
     nvols, mem_gb = _create_mem_gb(bold_file)
     if nvols <= 5 - config.execution.sloppy:
@@ -282,7 +287,7 @@ def init_bold_wf(
             bold_output=boldref_out,
             echo_output=echos_out,
             multiecho=multiecho,
-            all_metadata=[config.execution.layout.get_metadata(file) for file in bold_series],
+            all_metadata=all_metadata,
         )
         ds_bold_native_wf.inputs.inputnode.source_files = bold_series
 
@@ -297,7 +302,27 @@ def init_bold_wf(
             ]),
         ])  # fmt:skip
 
-    anat_out = bool(nonstd_spaces.intersection(('anat', 'T1w')))
+    if nonstd_spaces.intersection(('anat', 'T1w')):
+        ds_bold_t1_wf = init_ds_volumes_wf(
+            bids_root=str(config.execution.bids_dir),
+            output_dir=fmriprep_dir,
+            multiecho=multiecho,
+            metadata=metadata[0],
+        )
+        ds_bold_t1_wf.inputs.inputnode.source_files = bold_series
+
+        workflow.connect([
+            (inputnode, ds_bold_t1_wf, [
+                ('t1w_preproc', 'inputnode.ref_file'),
+            ]),
+            (bold_fit_wf, ds_bold_t1_wf, [
+                ('outputnode.bold_mask', 'inputnode.bold_mask'),
+                ('outputnode.coreg_boldref', 'inputnode.bold_ref'),
+                ('outputnode.boldref2anat_xfm', 'inputnode.boldref2anat_xfm'),
+            ]),
+            (bold_native_wf, ds_bold_t1_wf, [('outputnode.t2star_map', 'inputnode.t2star')]),
+            (bold_anat_wf, ds_bold_t1_wf, [('outputnode.bold_file', 'inputnode.bold')]),
+        ])  # fmt:skip
 
     if multiecho:
         t2s_reporting_wf = init_t2s_reporting_wf()
diff --git a/fmriprep/workflows/bold/outputs.py b/fmriprep/workflows/bold/outputs.py
index 8c593b42f..e10754127 100644
--- a/fmriprep/workflows/bold/outputs.py
+++ b/fmriprep/workflows/bold/outputs.py
@@ -687,7 +687,6 @@ def init_ds_volumes_wf(
     *,
     bids_root: str,
     output_dir: str,
-    space: str,
     multiecho: bool,
     metadata: ty.List[dict],
     name="ds_bold_native_wf",

From 7cc25d561483f5650146b34f7b66a46951fdf36a Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Wed, 18 Oct 2023 14:23:41 -0400
Subject: [PATCH 09/19] FIX: Names, connections, metadata

---
 fmriprep/workflows/bold/base.py    |  4 +++-
 fmriprep/workflows/bold/outputs.py | 14 +++++++++-----
 2 files changed, 12 insertions(+), 6 deletions(-)

diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py
index a849cbc22..b43a5c2f9 100644
--- a/fmriprep/workflows/bold/base.py
+++ b/fmriprep/workflows/bold/base.py
@@ -307,9 +307,11 @@ def init_bold_wf(
             bids_root=str(config.execution.bids_dir),
             output_dir=fmriprep_dir,
             multiecho=multiecho,
-            metadata=metadata[0],
+            metadata=all_metadata[0],
+            name='ds_bold_t1_wf',
         )
         ds_bold_t1_wf.inputs.inputnode.source_files = bold_series
+        ds_bold_t1_wf.inputs.inputnode.space = 'T1w'
 
         workflow.connect([
             (inputnode, ds_bold_t1_wf, [
diff --git a/fmriprep/workflows/bold/outputs.py b/fmriprep/workflows/bold/outputs.py
index e10754127..a46a4fbcb 100644
--- a/fmriprep/workflows/bold/outputs.py
+++ b/fmriprep/workflows/bold/outputs.py
@@ -689,9 +689,8 @@ def init_ds_volumes_wf(
     output_dir: str,
     multiecho: bool,
     metadata: ty.List[dict],
-    name="ds_bold_native_wf",
+    name="ds_volumes_wf",
 ) -> pe.Workflow:
-    metadata = all_metadata[0]
     timing_parameters = prepare_timing_parameters(metadata)
 
     workflow = pe.Workflow(name=name)
@@ -775,7 +774,7 @@ def init_ds_volumes_wf(
             compress=True,
             dismiss_entities=("echo",),
         ),
-        name='ds_bold',
+        name='ds_ref',
         run_without_submitting=True,
         mem_gb=DEFAULT_MEMORY_MIN_GB,
     )
@@ -787,13 +786,18 @@ def init_ds_volumes_wf(
             compress=True,
             dismiss_entities=("echo",),
         ),
-        name='ds_bold',
+        name='ds_mask',
         run_without_submitting=True,
         mem_gb=DEFAULT_MEMORY_MIN_GB,
     )
     datasinks = [ds_ref, ds_mask]
 
     if multiecho:
+        t2star_meta = {
+            'Units': 's',
+            'EstimationReference': 'doi:10.1002/mrm.20900',
+            'EstimationAlgorithm': 'monoexponential decay model',
+        }
         resample_t2star = pe.Node(
             ApplyTransforms(
                 dimension=3,
@@ -801,7 +805,7 @@ def init_ds_volumes_wf(
                 float=True,
                 interpolation="LanczosWindowedSinc",
             ),
-            name="resample_ref",
+            name="resample_t2star",
         )
         ds_t2star = pe.Node(
             DerivativesDataSink(

From ad3f392a7d1268b70353374a273d8bd7ca19dca2 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Mon, 23 Oct 2023 09:02:38 -0400
Subject: [PATCH 10/19] FIX volumetric

---
 fmriprep/workflows/bold/apply.py | 1 +
 fmriprep/workflows/bold/base.py  | 2 ++
 2 files changed, 3 insertions(+)

diff --git a/fmriprep/workflows/bold/apply.py b/fmriprep/workflows/bold/apply.py
index 9aff94921..b25029de8 100644
--- a/fmriprep/workflows/bold/apply.py
+++ b/fmriprep/workflows/bold/apply.py
@@ -27,6 +27,7 @@
 
 def init_bold_volumetric_resample_wf(
     *,
+    metadata: dict,
     fieldmap_id: ty.Optional[str] = None,
     omp_nthreads: int = 1,
     name: str = 'bold_volumetric_resample_wf',
diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py
index b43a5c2f9..f606f550e 100644
--- a/fmriprep/workflows/bold/base.py
+++ b/fmriprep/workflows/bold/base.py
@@ -51,6 +51,7 @@
 from .hmc import init_bold_hmc_wf
 from .outputs import (
     init_ds_bold_native_wf,
+    init_ds_registration_wf,
     init_ds_volumes_wf,
     init_func_derivatives_wf,
 )
@@ -243,6 +244,7 @@ def init_bold_wf(
         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,
         name='bold_anat_wf',

From cac0b79925ca3b7a73cadeb5883b0d9ba6bafe89 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Mon, 23 Oct 2023 10:27:25 -0400
Subject: [PATCH 11/19] FIX: fmapreg_xfm -> boldref2fmap_xfm

---
 fmriprep/workflows/bold/apply.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/fmriprep/workflows/bold/apply.py b/fmriprep/workflows/bold/apply.py
index b25029de8..c5e4f838f 100644
--- a/fmriprep/workflows/bold/apply.py
+++ b/fmriprep/workflows/bold/apply.py
@@ -99,7 +99,7 @@ def init_bold_volumetric_resample_wf(
                 ("fmap_id", "keys"),
             ]),
             (inputnode, distortion_params, [('bold_file', 'in_file')]),
-            (inputnode, fmap2target, [('fmapreg_xfm', 'in1')]),
+            (inputnode, fmap2target, [('boldref2fmap_xfm', 'in1')]),
             (boldref2target, fmap2target, [('out', 'in2')]),
             (boldref2target, inverses, [('out', 'inlist')]),
             (fmap_select, fmap_recon, [

From 143a41fc0836c98fe06c9b0db5b0a6a80a5dcd00 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Mon, 23 Oct 2023 12:05:53 -0400
Subject: [PATCH 12/19] FIX: Connect target_ref_file in volumetric resample
 workflow

---
 fmriprep/workflows/bold/apply.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/fmriprep/workflows/bold/apply.py b/fmriprep/workflows/bold/apply.py
index c5e4f838f..ca48c5004 100644
--- a/fmriprep/workflows/bold/apply.py
+++ b/fmriprep/workflows/bold/apply.py
@@ -100,6 +100,7 @@ def init_bold_volumetric_resample_wf(
             ]),
             (inputnode, distortion_params, [('bold_file', 'in_file')]),
             (inputnode, fmap2target, [('boldref2fmap_xfm', 'in1')]),
+            (inputnode, fmap_recon, [('ref_file', 'target_ref_file')]),
             (boldref2target, fmap2target, [('out', 'in2')]),
             (boldref2target, inverses, [('out', 'inlist')]),
             (fmap_select, fmap_recon, [

From 333a8fea73bd2a904d2cbe7444efdda84eb9480f Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Mon, 23 Oct 2023 12:20:46 -0400
Subject: [PATCH 13/19] RF: Move coordinate copies into threads to delay
 allocation

---
 fmriprep/interfaces/resampling.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/fmriprep/interfaces/resampling.py b/fmriprep/interfaces/resampling.py
index 9921d6f19..c052c113c 100644
--- a/fmriprep/interfaces/resampling.py
+++ b/fmriprep/interfaces/resampling.py
@@ -275,6 +275,9 @@ def resample_vol(
         coordinates = nb.affines.apply_affine(
             hmc_xfm, coordinates.reshape(coords_shape[0], -1).T
         ).T.reshape(coords_shape)
+    else:
+        # Copy coordinates to avoid interfering with other calls
+        coordinates = coordinates.copy()
 
     vsm = fmap_hz * pe_info[1]
     coordinates[pe_info[0], ...] += vsm
@@ -377,7 +380,7 @@ async def resample_series_async(
                 partial(
                     resample_vol,
                     data=volume,
-                    coordinates=coordinates.copy(),
+                    coordinates=coordinates,
                     pe_info=pe_info[volid],
                     hmc_xfm=hmc_xfms[volid] if hmc_xfms else None,
                     fmap_hz=fmap_hz,

From d2a3d958c92f14f504f2f3f3f668aa54312ed73e Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Mon, 23 Oct 2023 15:34:07 -0400
Subject: [PATCH 14/19] FIX: Add sampling reference

---
 fmriprep/workflows/bold/apply.py | 20 ++++++++++++++------
 fmriprep/workflows/bold/base.py  |  4 +++-
 2 files changed, 17 insertions(+), 7 deletions(-)

diff --git a/fmriprep/workflows/bold/apply.py b/fmriprep/workflows/bold/apply.py
index ca48c5004..d9624d97a 100644
--- a/fmriprep/workflows/bold/apply.py
+++ b/fmriprep/workflows/bold/apply.py
@@ -7,6 +7,7 @@
 import nipype.interfaces.utility as niu
 import nipype.pipeline.engine as pe
 from niworkflows.interfaces.header import ValidateImage
+from niworkflows.interfaces.nibabel import GenerateSamplingReference
 from niworkflows.interfaces.utility import KeySelect
 from niworkflows.utils.connections import listify
 
@@ -38,7 +39,9 @@ def init_bold_volumetric_resample_wf(
         niu.IdentityInterface(
             fields=[
                 "bold_file",
-                "ref_file",
+                "bold_ref_file",
+                "target_ref_file",
+                "target_mask",
                 # HMC
                 "motion_xfm",
                 # SDC
@@ -57,20 +60,25 @@ def init_bold_volumetric_resample_wf(
 
     outputnode = pe.Node(niu.IdentityInterface(fields=["bold_file"]), name='outputnode')
 
+    gen_ref = pe.Node(GenerateSamplingReference(), name='gen_ref', mem_gb=0.3)
+
     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, gen_ref, [
+            ('bold_ref_file', 'moving_image'),
+            ('target_ref_file', 'fixed_image'),
+            ('target_mask', 'fov_mask'),
+        ]),
         (inputnode, boldref2target, [
             ('boldref2anat_xfm', 'in1'),
             ('anat2std_xfm', 'in2'),
         ]),
         (inputnode, bold2target, [('motion_xfm', 'in1')]),
-        (inputnode, resample, [
-            ('bold_file', 'in_file'),
-            ('ref_file', 'ref_file'),
-        ]),
+        (inputnode, resample, [('bold_file', 'in_file')]),
+        (gen_ref, resample, [('out_file', 'ref_file')]),
         (boldref2target, bold2target, [('out', 'in2')]),
         (bold2target, resample, [('out', 'transforms')]),
         (resample, outputnode, [('out_file', 'bold_file')]),
@@ -100,7 +108,7 @@ def init_bold_volumetric_resample_wf(
             ]),
             (inputnode, distortion_params, [('bold_file', 'in_file')]),
             (inputnode, fmap2target, [('boldref2fmap_xfm', 'in1')]),
-            (inputnode, fmap_recon, [('ref_file', 'target_ref_file')]),
+            (gen_ref, fmap_recon, [('out_file', 'target_ref_file')]),
             (boldref2target, fmap2target, [('out', 'in2')]),
             (boldref2target, inverses, [('out', 'inlist')]),
             (fmap_select, fmap_recon, [
diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py
index f606f550e..f95e740ee 100644
--- a/fmriprep/workflows/bold/base.py
+++ b/fmriprep/workflows/bold/base.py
@@ -257,7 +257,8 @@ def init_bold_wf(
             ("fmap_id", "inputnode.fmap_id"),
         ]),
         (inputnode, bold_anat_wf, [
-            ("t1w_preproc", "inputnode.ref_file"),
+            ("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"),
@@ -270,6 +271,7 @@ def init_bold_wf(
             ("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"),
         ]),

From cf92a48c4a43433ecdd31612ed8b9322570a3d81 Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Mon, 23 Oct 2023 15:36:09 -0400
Subject: [PATCH 15/19] ENH: Construct resampled data array in F order

Ensures that each volume is contiguous in memory, which should keep it
more efficient for both writing resampled volumes into the array and
writing the array to disk.
---
 fmriprep/interfaces/resampling.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/fmriprep/interfaces/resampling.py b/fmriprep/interfaces/resampling.py
index c052c113c..fd521a6ce 100644
--- a/fmriprep/interfaces/resampling.py
+++ b/fmriprep/interfaces/resampling.py
@@ -372,7 +372,7 @@ async def resample_series_async(
 
     semaphore = asyncio.Semaphore(max_concurrent)
 
-    out_array = np.zeros(coordinates.shape[1:] + data.shape[-1:], dtype=output_dtype)
+    out_array = np.zeros(coordinates.shape[1:] + data.shape[-1:], dtype=output_dtype, order='F')
 
     tasks = [
         asyncio.create_task(

From 83b83064490d31c2f5a5730fd3e89e9910b88421 Mon Sep 17 00:00:00 2001
From: Chris Markiewicz <effigies@gmail.com>
Date: Wed, 18 Oct 2023 15:05:47 -0400
Subject: [PATCH 16/19] LOG: Move logging into init_bold_wf

---
 fmriprep/workflows/bold/base.py | 19 ++++++++++---------
 1 file changed, 10 insertions(+), 9 deletions(-)

diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py
index f95e740ee..4d95f1cfd 100644
--- a/fmriprep/workflows/bold/base.py
+++ b/fmriprep/workflows/bold/base.py
@@ -154,6 +154,16 @@ def init_bold_wf(
         )
         return
 
+    config.loggers.workflow.debug(
+        "Creating bold processing workflow for <%s> (%.2f GB / %d TRs). "
+        "Memory resampled/largemem=%.2f/%.2f GB.",
+        bold_file,
+        mem_gb["filesize"],
+        nvols,
+        mem_gb["resampled"],
+        mem_gb["largemem"],
+    )
+
     functional_cache = {}
     if config.execution.derivatives:
         from fmriprep.utils.bids import collect_derivatives, extract_entities
@@ -524,15 +534,6 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
 
     ref_file = bold_file
     wf_name = _get_wf_name(ref_file, "func_preproc")
-    config.loggers.workflow.debug(
-        "Creating bold processing workflow for <%s> (%.2f GB / %d TRs). "
-        "Memory resampled/largemem=%.2f/%.2f GB.",
-        ref_file,
-        mem_gb["filesize"],
-        bold_tlen,
-        mem_gb["resampled"],
-        mem_gb["largemem"],
-    )
 
     # Build workflow
     workflow = Workflow(name=wf_name)

From 9bcfdace0f988dd8f507564e6c2fee7f5b6d6136 Mon Sep 17 00:00:00 2001
From: Chris Markiewicz <effigies@gmail.com>
Date: Wed, 18 Oct 2023 15:07:11 -0400
Subject: [PATCH 17/19] RF: Remove old t1 workflow

---
 fmriprep/workflows/bold/base.py | 48 ---------------------------------
 1 file changed, 48 deletions(-)

diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py
index 4d95f1cfd..673d16e7f 100644
--- a/fmriprep/workflows/bold/base.py
+++ b/fmriprep/workflows/bold/base.py
@@ -615,16 +615,6 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
         name="outputnode",
     )
 
-    # apply BOLD registration to T1w
-    bold_t1_trans_wf = init_bold_t1_trans_wf(
-        name="bold_t1_trans_wf",
-        freesurfer=freesurfer,
-        mem_gb=mem_gb["resampled"],
-        omp_nthreads=omp_nthreads,
-        use_compression=False,
-    )
-    bold_t1_trans_wf.inputs.inputnode.fieldwarp = "identity"
-
     # get confounds
     bold_confounds_wf = init_bold_confs_wf(
         mem_gb=mem_gb["largemem"],
@@ -637,44 +627,6 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
     )
     bold_confounds_wf.get_node("inputnode").inputs.t1_transform_flags = [False]
 
-    # Map final BOLD mask into T1w space (if required)
-    nonstd_spaces = set(spaces.get_nonstandard())
-    if nonstd_spaces.intersection(("T1w", "anat")):
-        from niworkflows.interfaces.fixes import (
-            FixHeaderApplyTransforms as ApplyTransforms,
-        )
-
-        boldmask_to_t1w = pe.Node(
-            ApplyTransforms(interpolation="MultiLabel"),
-            name="boldmask_to_t1w",
-            mem_gb=0.1,
-        )
-        # fmt:off
-        workflow.connect([
-            (bold_reg_wf, boldmask_to_t1w, [("outputnode.itk_bold_to_t1", "transforms")]),
-            (bold_t1_trans_wf, boldmask_to_t1w, [("outputnode.bold_mask_t1", "reference_image")]),
-            (bold_final, boldmask_to_t1w, [("mask", "input_image")]),
-            (boldmask_to_t1w, outputnode, [("output_image", "bold_mask_t1")]),
-        ])
-        # fmt:on
-
-        if multiecho:
-            t2star_to_t1w = pe.Node(
-                ApplyTransforms(interpolation="LanczosWindowedSinc", float=True),
-                name="t2star_to_t1w",
-                mem_gb=0.1,
-            )
-            # fmt:off
-            workflow.connect([
-                (bold_reg_wf, t2star_to_t1w, [("outputnode.itk_bold_to_t1", "transforms")]),
-                (bold_t1_trans_wf, t2star_to_t1w, [
-                    ("outputnode.bold_mask_t1", "reference_image")
-                ]),
-                (bold_final, t2star_to_t1w, [("t2star", "input_image")]),
-                (t2star_to_t1w, outputnode, [("output_image", "t2star_t1")]),
-            ])
-            # fmt:on
-
     if spaces.get_spaces(nonstandard=False, dim=(3,)):
         # Apply transforms in 1 shot
         bold_std_trans_wf = init_bold_std_trans_wf(

From 26d141fe8826e20ca828abe60694c42e4b6f0f15 Mon Sep 17 00:00:00 2001
From: Chris Markiewicz <effigies@gmail.com>
Date: Thu, 19 Oct 2023 14:27:41 -0400
Subject: [PATCH 18/19] RF: Move resampling to "full" section

---
 fmriprep/workflows/bold/base.py    | 50 ++++++++++++++++--------------
 fmriprep/workflows/bold/outputs.py | 42 ++++++++++++-------------
 2 files changed, 47 insertions(+), 45 deletions(-)

diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py
index 673d16e7f..da6052cd6 100644
--- a/fmriprep/workflows/bold/base.py
+++ b/fmriprep/workflows/bold/base.py
@@ -292,6 +292,7 @@ def init_bold_wf(
     ])  # fmt:skip
 
     boldref_out = bool(nonstd_spaces.intersection(('func', 'run', 'bold', 'boldref', 'sbref')))
+    boldref_out |= config.workflow.level == 'full'
     echos_out = multiecho and config.execution.me_output_echos
 
     if boldref_out or echos_out:
@@ -316,30 +317,6 @@ def init_bold_wf(
             ]),
         ])  # fmt:skip
 
-    if nonstd_spaces.intersection(('anat', 'T1w')):
-        ds_bold_t1_wf = init_ds_volumes_wf(
-            bids_root=str(config.execution.bids_dir),
-            output_dir=fmriprep_dir,
-            multiecho=multiecho,
-            metadata=all_metadata[0],
-            name='ds_bold_t1_wf',
-        )
-        ds_bold_t1_wf.inputs.inputnode.source_files = bold_series
-        ds_bold_t1_wf.inputs.inputnode.space = 'T1w'
-
-        workflow.connect([
-            (inputnode, ds_bold_t1_wf, [
-                ('t1w_preproc', 'inputnode.ref_file'),
-            ]),
-            (bold_fit_wf, ds_bold_t1_wf, [
-                ('outputnode.bold_mask', 'inputnode.bold_mask'),
-                ('outputnode.coreg_boldref', 'inputnode.bold_ref'),
-                ('outputnode.boldref2anat_xfm', 'inputnode.boldref2anat_xfm'),
-            ]),
-            (bold_native_wf, ds_bold_t1_wf, [('outputnode.t2star_map', 'inputnode.t2star')]),
-            (bold_anat_wf, ds_bold_t1_wf, [('outputnode.bold_file', 'inputnode.bold')]),
-        ])  # fmt:skip
-
     if multiecho:
         t2s_reporting_wf = init_t2s_reporting_wf()
 
@@ -379,6 +356,31 @@ def init_bold_wf(
     if config.workflow.level == "resampling":
         return workflow
 
+    # Full derivatives, including resampled BOLD series
+    if nonstd_spaces.intersection(('anat', 'T1w')):
+        ds_bold_t1_wf = init_ds_volumes_wf(
+            bids_root=str(config.execution.bids_dir),
+            output_dir=fmriprep_dir,
+            multiecho=multiecho,
+            metadata=all_metadata[0],
+            name='ds_bold_t1_wf',
+        )
+        ds_bold_t1_wf.inputs.inputnode.source_files = bold_series
+        ds_bold_t1_wf.inputs.inputnode.space = 'T1w'
+
+        workflow.connect([
+            (inputnode, ds_bold_t1_wf, [
+                ('t1w_preproc', 'inputnode.ref_file'),
+            ]),
+            (bold_fit_wf, ds_bold_t1_wf, [
+                ('outputnode.bold_mask', 'inputnode.bold_mask'),
+                ('outputnode.coreg_boldref', 'inputnode.bold_ref'),
+                ('outputnode.boldref2anat_xfm', 'inputnode.boldref2anat_xfm'),
+            ]),
+            (bold_native_wf, ds_bold_t1_wf, [('outputnode.t2star_map', 'inputnode.t2star')]),
+            (bold_anat_wf, ds_bold_t1_wf, [('outputnode.bold_file', 'inputnode.bold')]),
+        ])  # fmt:skip
+
     # Fill-in datasinks of reportlets seen so far
     for node in workflow.list_node_names():
         if node.split(".")[-1].startswith("ds_report"):
diff --git a/fmriprep/workflows/bold/outputs.py b/fmriprep/workflows/bold/outputs.py
index a46a4fbcb..9f141f17c 100644
--- a/fmriprep/workflows/bold/outputs.py
+++ b/fmriprep/workflows/bold/outputs.py
@@ -587,6 +587,27 @@ def init_ds_bold_native_wf(
     raw_sources.inputs.bids_root = bids_root
     workflow.connect(inputnode, 'source_files', raw_sources, 'in_files')
 
+    # Masks should be output if any other derivatives are output
+    ds_bold_mask = pe.Node(
+        DerivativesDataSink(
+            base_directory=output_dir,
+            desc='brain',
+            suffix='mask',
+            compress=True,
+            dismiss_entities=("echo",),
+        ),
+        name='ds_bold_mask',
+        run_without_submitting=True,
+        mem_gb=DEFAULT_MEMORY_MIN_GB,
+    )
+    workflow.connect([
+        (inputnode, ds_bold_mask, [
+            ('source_files', 'source_file'),
+            ('bold_mask', 'in_file'),
+        ]),
+        (raw_sources, ds_bold_mask, [('out', 'RawSources')]),
+    ])  # fmt:skip
+
     if bold_output:
         ds_bold = pe.Node(
             DerivativesDataSink(
@@ -609,27 +630,6 @@ def init_ds_bold_native_wf(
             ]),
         ])  # fmt:skip
 
-    if bold_output or echo_output:
-        ds_bold_mask = pe.Node(
-            DerivativesDataSink(
-                base_directory=output_dir,
-                desc='brain',
-                suffix='mask',
-                compress=True,
-                dismiss_entities=("echo",),
-            ),
-            name='ds_bold_mask',
-            run_without_submitting=True,
-            mem_gb=DEFAULT_MEMORY_MIN_GB,
-        )
-        workflow.connect([
-            (inputnode, ds_bold_mask, [
-                ('source_files', 'source_file'),
-                ('bold_mask', 'in_file'),
-            ]),
-            (raw_sources, ds_bold_mask, [('out', 'RawSources')]),
-        ])  # fmt:skip
-
     if bold_output and multiecho:
         t2star_meta = {
             'Units': 's',

From 0d09cc051576edfecee509d411f9f2c9a90f0a58 Mon Sep 17 00:00:00 2001
From: Chris Markiewicz <markiewicz@stanford.edu>
Date: Thu, 26 Oct 2023 12:45:57 -0400
Subject: [PATCH 19/19] Apply suggestions from code review

Co-authored-by: Mathias Goncalves <goncalves.mathias@gmail.com>
---
 fmriprep/interfaces/resampling.py |  2 +
 fmriprep/workflows/bold/apply.py  | 84 ++++++++++++++++---------------
 2 files changed, 45 insertions(+), 41 deletions(-)

diff --git a/fmriprep/interfaces/resampling.py b/fmriprep/interfaces/resampling.py
index fd521a6ce..d3febb00c 100644
--- a/fmriprep/interfaces/resampling.py
+++ b/fmriprep/interfaces/resampling.py
@@ -372,6 +372,8 @@ async def resample_series_async(
 
     semaphore = asyncio.Semaphore(max_concurrent)
 
+    # Order F ensures individual volumes are contiguous in memory
+    # Also matches NIfTI, making final save more efficient
     out_array = np.zeros(coordinates.shape[1:] + data.shape[-1:], dtype=output_dtype, order='F')
 
     tasks = [
diff --git a/fmriprep/workflows/bold/apply.py b/fmriprep/workflows/bold/apply.py
index d9624d97a..6024a73d0 100644
--- a/fmriprep/workflows/bold/apply.py
+++ b/fmriprep/workflows/bold/apply.py
@@ -29,7 +29,7 @@
 def init_bold_volumetric_resample_wf(
     *,
     metadata: dict,
-    fieldmap_id: ty.Optional[str] = None,
+    fieldmap_id: str | None = None,
     omp_nthreads: int = 1,
     name: str = 'bold_volumetric_resample_wf',
 ) -> pe.Workflow:
@@ -84,46 +84,48 @@ def init_bold_volumetric_resample_wf(
         (resample, outputnode, [('out_file', 'bold_file')]),
     ])  # 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, [('boldref2fmap_xfm', 'in1')]),
-            (gen_ref, fmap_recon, [('out_file', 'target_ref_file')]),
-            (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
+    if not fieldmap_id:
+        return workflow
+
+    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, [('boldref2fmap_xfm', 'in1')]),
+        (gen_ref, fmap_recon, [('out_file', 'target_ref_file')]),
+        (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