diff --git a/env.yml b/env.yml
index b62b48fb..5ff927a2 100644
--- a/env.yml
+++ b/env.yml
@@ -29,8 +29,6 @@ dependencies:
 - ants=2.6
 # Workflow dependencies: Connectome Workbench
 - connectome-workbench-cli=2.0
-# Workflow dependencies: Convert3d
-- convert3d=1.4
 # Workflow dependencies: FSL (versions pinned in 6.0.7.17.20250415.fe1c582e)
 - fsl-bet2=2111.8
 - fsl-flirt=2111.4
diff --git a/nibabies/interfaces/resampling.py b/nibabies/interfaces/resampling.py
index eb9f7998..05df737f 100644
--- a/nibabies/interfaces/resampling.py
+++ b/nibabies/interfaces/resampling.py
@@ -2,12 +2,11 @@
 
 import asyncio
 import os
-from collections.abc import Callable
 from functools import partial
-from typing import TypeVar
 
 import nibabel as nb
 import nitransforms as nt
+import nitransforms.resampling
 import numpy as np
 from nipype.interfaces.base import (
     File,
@@ -22,15 +21,8 @@
 from sdcflows.transform import grid_bspline_weights
 from sdcflows.utils.tools import ensure_positive_cosines
 
-from nibabies.utils.transforms import load_transforms
-
-R = TypeVar('R')
-
-
-async def worker(job: Callable[[], R], semaphore: asyncio.Semaphore) -> R:
-    async with semaphore:
-        loop = asyncio.get_running_loop()
-        return await loop.run_in_executor(None, job)
+from ..utils.asynctools import worker
+from ..utils.transforms import load_transforms
 
 
 class ResampleSeriesInputSpec(TraitedSpec):
@@ -38,7 +30,6 @@ class ResampleSeriesInputSpec(TraitedSpec):
     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)',
     )
     inverse = InputMultiObject(
@@ -100,7 +91,8 @@ def _run_interface(self, runtime):
 
         nvols = source.shape[3] if source.ndim > 3 else 1
 
-        transforms = load_transforms(self.inputs.transforms, self.inputs.inverse)
+        # No transforms appear Undefined, pass as empty list
+        transforms = load_transforms(self.inputs.transforms or [], self.inputs.inverse)
 
         pe_dir = self.inputs.pe_dir
         ro_time = self.inputs.ro_time
@@ -199,6 +191,13 @@ def _run_interface(self, runtime):
 class DistortionParametersInputSpec(TraitedSpec):
     in_file = File(exists=True, desc='EPI image corresponding to the metadata')
     metadata = traits.Dict(mandatory=True, desc='metadata corresponding to the inputs')
+    fallback = traits.Either(
+        None,
+        'estimated',
+        traits.Float,
+        usedefault=True,
+        desc='Fallback value for missing metadata',
+    )
 
 
 class DistortionParametersOutputSpec(TraitedSpec):
@@ -223,6 +222,8 @@ def _run_interface(self, runtime):
             self._results['readout_time'] = get_trt(
                 self.inputs.metadata,
                 self.inputs.in_file or None,
+                use_estimate=self.inputs.fallback == 'estimated',
+                fallback=self.inputs.fallback if isinstance(self.inputs.fallback, float) else None,
             )
             self._results['pe_direction'] = self.inputs.metadata['PhaseEncodingDirection']
         except (KeyError, ValueError):
@@ -686,9 +687,6 @@ def reconstruct_fieldmap(
             target.__class__(target.dataobj, projected_affine, target.header),
         )
     else:
-        # Below if statement was taken from fmriprep pull request #3439,
-        # along with the same explanation:
-        #
         # Hack. Sometimes the reference array is rotated relative to the fieldmap
         # and coefficient grids. As far as I know, coefficients are always RAS,
         # but good to check before doing this.
@@ -717,7 +715,7 @@ def reconstruct_fieldmap(
     )
 
     if not direct:
-        fmap_img = nt.apply(transforms, fmap_img, reference=target)
+        fmap_img = nt.resampling.apply(transforms, fmap_img, reference=target)
 
     fmap_img.header.set_intent('estimate', name='fieldmap Hz')
     fmap_img.header.set_data_dtype('float32')
diff --git a/nibabies/workflows/bold/base.py b/nibabies/workflows/bold/base.py
index 745cbb16..b9f20d92 100644
--- a/nibabies/workflows/bold/base.py
+++ b/nibabies/workflows/bold/base.py
@@ -734,8 +734,8 @@ def init_bold_wf(
         ]),
         (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.hmc_boldref', 'inputnode.hmc_boldref'),
+            ('outputnode.motion_xfm', 'inputnode.motion_xfm'),
             ('outputnode.boldref2anat_xfm', 'inputnode.boldref2anat_xfm'),
             ('outputnode.dummy_scans', 'inputnode.skip_vols'),
         ]),
diff --git a/nibabies/workflows/bold/confounds.py b/nibabies/workflows/bold/confounds.py
index 2b0d46eb..672ceb63 100644
--- a/nibabies/workflows/bold/confounds.py
+++ b/nibabies/workflows/bold/confounds.py
@@ -31,6 +31,11 @@
 from nipype.algorithms import confounds as nac
 from nipype.interfaces import utility as niu
 from nipype.pipeline import engine as pe
+from niworkflows.interfaces.confounds import (
+    FramewiseDisplacement,
+    FSLMotionParams,
+    FSLRMSDeviation,
+)
 
 from nibabies.config import DEFAULT_DISMISS_ENTITIES, DEFAULT_MEMORY_MIN_GB
 from nibabies.interfaces import DerivativesDataSink
@@ -118,10 +123,10 @@ def init_bold_confs_wf(
         when available.
     bold_mask
         BOLD series mask
-    movpar_file
-        SPM-formatted motion parameters file
-    rmsd_file
-        Root mean squared deviation as measured by ``fsl_motion_outliers`` [Jenkinson2002]_.
+    motion_xfm
+        ITK-formatted head motion transforms
+    hmc_boldref
+        BOLD volume after HMC
     skip_vols
         number of non steady state volumes
     anat_mask
@@ -219,8 +224,8 @@ def init_bold_confs_wf(
             fields=[
                 'bold',
                 'bold_mask',
-                'movpar_file',
-                'rmsd_file',
+                'motion_xfm',
+                'hmc_boldref',
                 'skip_vols',
                 'anat_mask',
                 'anat_tpms',
@@ -260,8 +265,12 @@ def init_bold_confs_wf(
         mem_gb=mem_gb,
     )
 
+    # Motion parameters
+    motion_params = pe.Node(FSLMotionParams(), name='motion_params')
+
     # Frame displacement
-    fdisp = pe.Node(nac.FramewiseDisplacement(parameter_source='SPM'), name='fdisp', mem_gb=mem_gb)
+    fdisp = pe.Node(FramewiseDisplacement(), name='fdisp')
+    rmsd = pe.Node(FSLRMSDeviation(), name='rmsd')
 
     # Generate aCompCor probseg maps
     acc_masks = pe.Node(aCompCorMasks(is_aseg=freesurfer), name='acc_masks')
@@ -365,18 +374,6 @@ def init_bold_confs_wf(
         mem_gb=0.01,
         run_without_submitting=True,
     )
-    add_motion_headers = pe.Node(
-        AddTSVHeader(columns=['trans_x', 'trans_y', 'trans_z', 'rot_x', 'rot_y', 'rot_z']),
-        name='add_motion_headers',
-        mem_gb=0.01,
-        run_without_submitting=True,
-    )
-    add_rmsd_header = pe.Node(
-        AddTSVHeader(columns=['rmsd']),
-        name='add_rmsd_header',
-        mem_gb=0.01,
-        run_without_submitting=True,
-    )
     concat = pe.Node(GatherConfounds(), name='concat', mem_gb=0.01, run_without_submitting=True)
 
     # CompCor metadata
@@ -529,7 +526,11 @@ def _select_cols(table):
         # connect inputnode to each non-anatomical confound node
         (inputnode, dvars, [('bold', 'in_file'),
                             ('bold_mask', 'in_mask')]),
-        (inputnode, fdisp, [('movpar_file', 'in_file')]),
+        (inputnode, motion_params, [('motion_xfm', 'xfm_file'),
+                                    ('hmc_boldref', 'boldref_file')]),
+        (inputnode, rmsd, [('motion_xfm', 'xfm_file'),
+                           ('hmc_boldref', 'boldref_file')]),
+        (motion_params, fdisp, [('out_file', 'in_file')]),
         # Brain mask
         (inputnode, anat_mask_tfm, [('anat_mask', 'input_image'),
                                    ('bold_mask', 'reference_image'),
@@ -572,8 +573,6 @@ def _select_cols(table):
         (merge_rois, signals, [('out', 'label_files')]),
 
         # Collate computed confounds together
-        (inputnode, add_motion_headers, [('movpar_file', 'in_file')]),
-        (inputnode, add_rmsd_header, [('rmsd_file', 'in_file')]),
         (dvars, add_dvars_header, [('out_nstd', 'in_file')]),
         (dvars, add_std_dvars_header, [('out_std', 'in_file')]),
         (signals, concat, [('out_file', 'signals')]),
@@ -582,8 +581,8 @@ def _select_cols(table):
                             ('pre_filter_file', 'cos_basis')]),
         (rename_acompcor, concat, [('components_file', 'acompcor')]),
         (crowncompcor, concat, [('components_file', 'crowncompcor')]),
-        (add_motion_headers, concat, [('out_file', 'motion')]),
-        (add_rmsd_header, concat, [('out_file', 'rmsd')]),
+        (motion_params, concat, [('out_file', 'motion')]),
+        (rmsd, concat, [('out_file', 'rmsd')]),
         (add_dvars_header, concat, [('out_file', 'dvars')]),
         (add_std_dvars_header, concat, [('out_file', 'std_dvars')]),
 
diff --git a/nibabies/workflows/bold/fit.py b/nibabies/workflows/bold/fit.py
index 4b9e8da6..a79ceb00 100644
--- a/nibabies/workflows/bold/fit.py
+++ b/nibabies/workflows/bold/fit.py
@@ -21,16 +21,16 @@
 #     https://www.nipreps.org/community/licensing/
 #
 import os
+import re
 import typing as ty
 
 import nibabel as nb
 from nipype.interfaces import utility as niu
 from nipype.pipeline import engine as pe
-from niworkflows.func.util import init_enhance_and_skullstrip_bold_wf
+from niworkflows.func.util import init_enhance_and_skullstrip_bold_wf, init_skullstrip_bold_wf
 from niworkflows.interfaces.header import ValidateImage
 from niworkflows.interfaces.nitransforms import ConcatenateXFMs
 from niworkflows.interfaces.utility import KeySelect
-from sdcflows.workflows.apply.correction import init_unwarp_wf
 from sdcflows.workflows.apply.registration import init_coeff2epi_wf
 
 from nibabies import config
@@ -50,12 +50,13 @@
 # BOLD workflows
 from nibabies.workflows.bold.hmc import init_bold_hmc_wf
 from nibabies.workflows.bold.outputs import (
+    init_ds_boldmask_wf,
     init_ds_boldref_wf,
     init_ds_hmc_wf,
     init_ds_registration_wf,
     init_func_fit_reports_wf,
 )
-from nibabies.workflows.bold.reference import init_raw_boldref_wf
+from nibabies.workflows.bold.reference import init_raw_boldref_wf, init_validation_and_dummies_wf
 from nibabies.workflows.bold.registration import init_bold_reg_wf
 from nibabies.workflows.bold.stc import init_bold_stc_wf
 from nibabies.workflows.bold.t2s import init_bold_t2s_wf
@@ -100,6 +101,7 @@ def init_bold_fit_wf(
     reference_anat: Anatomical,
     precomputed: dict | None = None,
     fieldmap_id: str | None = None,
+    jacobian: bool = False,
     omp_nthreads: int = 1,
     name: str = 'bold_fit_wf',
 ) -> pe.Workflow:
@@ -181,10 +183,6 @@ 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.
 
@@ -243,8 +241,8 @@ def init_bold_fit_wf(
     # Boolean used to update workflow self-descriptions
     multiecho = len(bold_series) > 1
 
-    have_hmcref = 'hmc_boldref' in precomputed
-    have_coregref = 'coreg_boldref' in precomputed
+    hmc_boldref = precomputed.get('hmc_boldref')
+    coreg_boldref = precomputed.get('coreg_boldref')
     # Can contain
     #  1) boldref2fmap
     #  2) boldref2anat
@@ -290,8 +288,6 @@ def init_bold_fit_wf(
                 'motion_xfm',
                 'boldref2anat_xfm',
                 'boldref2fmap_xfm',
-                'movpar_file',
-                'rmsd_file',
             ],
         ),
         name='outputnode',
@@ -305,9 +301,7 @@ def init_bold_fit_wf(
         name='hmcref_buffer',
     )
     fmapref_buffer = pe.Node(niu.Function(function=_select_ref), name='fmapref_buffer')
-    hmc_buffer = pe.Node(
-        niu.IdentityInterface(fields=['hmc_xforms', 'movpar_file', 'rmsd_file']), name='hmc_buffer'
-    )
+    hmc_buffer = pe.Node(niu.IdentityInterface(fields=['hmc_xforms']), name='hmc_buffer')
     fmapreg_buffer = pe.Node(
         niu.IdentityInterface(fields=['boldref2fmap_xfm']), name='fmapreg_buffer'
     )
@@ -315,6 +309,20 @@ def init_bold_fit_wf(
         niu.IdentityInterface(fields=['boldref', 'boldmask']), name='regref_buffer'
     )
 
+    if hmc_boldref:
+        hmcref_buffer.inputs.boldref = hmc_boldref
+        config.loggers.workflow.debug(f'Reusing motion correction reference: {hmc_boldref}')
+    if hmc_xforms:
+        hmc_buffer.inputs.hmc_xforms = hmc_xforms
+        config.loggers.workflow.debug(f'Reusing motion correction transforms: {hmc_xforms}')
+    if boldref2fmap_xform:
+        fmapreg_buffer.inputs.boldref2fmap_xfm = boldref2fmap_xform
+        config.loggers.workflow.debug(f'Reusing BOLD-to-fieldmap transform: {boldref2fmap_xform}')
+    if coreg_boldref:
+        regref_buffer.inputs.boldref = coreg_boldref
+        config.loggers.workflow.debug(f'Reusing coregistration reference: {coreg_boldref}')
+    fmapref_buffer.inputs.sbref_files = sbref_files
+
     use_fs = config.workflow.surface_recon_method is not None
     summary = pe.Node(
         FunctionalSummary(
@@ -343,8 +351,7 @@ def init_bold_fit_wf(
 
     func_fit_reports_wf = init_func_fit_reports_wf(
         reference_anat=reference_anat,
-        # TODO: Enable sdc report even if we find coregref
-        sdc_correction=not (have_coregref or fieldmap_id is None),
+        sdc_correction=fieldmap_id is None,
         output_dir=config.execution.output_dir,
     )
 
@@ -353,16 +360,13 @@ def init_bold_fit_wf(
             ('boldref', 'hmc_boldref'),
             ('dummy_scans', 'dummy_scans'),
         ]),
+        (hmcref_buffer, fmapref_buffer, [('boldref', 'boldref_files')]),
         (regref_buffer, outputnode, [
             ('boldref', 'coreg_boldref'),
             ('boldmask', 'bold_mask'),
         ]),
         (fmapreg_buffer, outputnode, [('boldref2fmap_xfm', 'boldref2fmap_xfm')]),
-        (hmc_buffer, outputnode, [
-            ('hmc_xforms', 'motion_xfm'),
-            ('movpar_file', 'movpar_file'),
-            ('rmsd_file', 'rmsd_file'),
-        ]),
+        (hmc_buffer, outputnode, [('hmc_xforms', 'motion_xfm')]),
         (inputnode, func_fit_reports_wf, [
             ('bold_file', 'inputnode.source_file'),
             ('anat_preproc', 'inputnode.anat_preproc'),
@@ -385,7 +389,7 @@ def init_bold_fit_wf(
         niu.IdentityInterface(fields=['in_file']),
         name='hmc_boldref_source_buffer',
     )
-    if not have_hmcref:
+    if not hmc_boldref:
         config.loggers.workflow.info('Stage 1: Adding HMC boldref workflow')
         hmc_boldref_wf = init_raw_boldref_wf(
             name='hmc_boldref_wf',
@@ -396,19 +400,20 @@ def init_bold_fit_wf(
         hmc_boldref_wf.inputs.inputnode.dummy_scans = config.workflow.dummy_scans
 
         ds_hmc_boldref_wf = init_ds_boldref_wf(
-            output_dir=config.execution.output_dir,
+            source_file=bold_file,
+            output_dir=config.execution.nibabies_dir,
             desc='hmc',
             name='ds_hmc_boldref_wf',
         )
         ds_hmc_boldref_wf.inputs.inputnode.source_files = [bold_file]
 
         workflow.connect([
+            (hmc_boldref_wf, ds_hmc_boldref_wf, [('outputnode.boldref', 'inputnode.boldref')]),
+            (ds_hmc_boldref_wf, hmcref_buffer, [('outputnode.boldref', 'boldref')]),
             (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')]),
             (hmc_boldref_wf, func_fit_reports_wf, [
                 ('outputnode.validation_report', 'inputnode.validation_report'),
@@ -420,14 +425,16 @@ def init_bold_fit_wf(
     else:
         config.loggers.workflow.info('Found HMC boldref - skipping Stage 1')
 
-        validate_bold = pe.Node(ValidateImage(), name='validate_bold')
-        validate_bold.inputs.in_file = bold_file
-
-        hmcref_buffer.inputs.boldref = precomputed['hmc_boldref']
+        validation_and_dummies_wf = init_validation_and_dummies_wf(bold_file=bold_file)
 
         workflow.connect([
-            (validate_bold, hmcref_buffer, [('out_file', 'bold_file')]),
-            (validate_bold, func_fit_reports_wf, [('out_report', 'inputnode.validation_report')]),
+            (validation_and_dummies_wf, hmcref_buffer, [
+                ('outputnode.bold_file', 'bold_file'),
+                ('outputnode.skip_vols', 'dummy_scans'),
+            ]),
+            (validation_and_dummies_wf, func_fit_reports_wf, [
+                ('outputnode.validation_report', 'inputnode.validation_report'),
+            ]),
             (hmcref_buffer, hmc_boldref_source_buffer, [('boldref', 'in_file')]),
         ])  # fmt:skip
 
@@ -438,7 +445,10 @@ def init_bold_fit_wf(
             name='bold_hmc_wf', mem_gb=mem_gb['filesize'], omp_nthreads=omp_nthreads
         )
 
-        ds_hmc_wf = init_ds_hmc_wf(output_dir=config.execution.output_dir)
+        ds_hmc_wf = init_ds_hmc_wf(
+            source_file=bold_file,
+            output_dir=config.execution.nibabies_dir,
+        )
         ds_hmc_wf.inputs.inputnode.source_files = [bold_file]
 
         workflow.connect([
@@ -447,141 +457,207 @@ def init_bold_fit_wf(
                 ('bold_file', 'inputnode.bold_file'),
             ]),
             (bold_hmc_wf, ds_hmc_wf, [('outputnode.xforms', 'inputnode.xforms')]),
-            (bold_hmc_wf, hmc_buffer, [
-                ('outputnode.movpar_file', 'movpar_file'),
-                ('outputnode.rmsd_file', 'rmsd_file'),
-            ]),
             (ds_hmc_wf, hmc_buffer, [('outputnode.xforms', 'hmc_xforms')]),
         ])  # fmt:skip
     else:
         config.loggers.workflow.info('Found motion correction transforms - skipping Stage 2')
-        hmc_buffer.inputs.hmc_xforms = hmc_xforms
 
-    # Stage 3: Create coregistration reference
-    # Fieldmap correction only happens during fit if this stage is needed
-    if not have_coregref:
-        config.loggers.workflow.info('Stage 3: Adding coregistration boldref workflow')
+    # Stage 3: Register fieldmap to boldref and reconstruct in BOLD space
+    if fieldmap_id:
+        config.loggers.workflow.info('Stage 3: Adding fieldmap reconstruction workflow')
+        fmap_select = pe.Node(
+            KeySelect(
+                fields=['fmap_ref', 'fmap_coeff', 'fmap_mask', 'sdc_method'],
+                key=fieldmap_id,
+            ),
+            name='fmap_select',
+            run_without_submitting=True,
+        )
+
+        boldref_fmap = pe.Node(ReconstructFieldmap(inverse=[True]), name='boldref_fmap', mem_gb=1)
 
-        # Select initial boldref, enhance contrast, and generate mask
-        fmapref_buffer.inputs.sbref_files = sbref_files
+        workflow.connect([
+            (inputnode, fmap_select, [
+                ('fmap_ref', 'fmap_ref'),
+                ('fmap_coeff', 'fmap_coeff'),
+                ('fmap_mask', 'fmap_mask'),
+                ('sdc_method', 'sdc_method'),
+                ('fmap_id', 'keys'),
+            ]),
+            (fmapref_buffer, boldref_fmap, [('out', 'target_ref_file')]),
+            (fmapreg_buffer, boldref_fmap, [('boldref2fmap_xfm', 'transforms')]),
+            (fmap_select, boldref_fmap, [
+                ('fmap_coeff', 'in_coeffs'),
+                ('fmap_ref', 'fmap_ref_file'),
+            ]),
+            (fmap_select, func_fit_reports_wf, [('fmap_ref', 'inputnode.fmap_ref')]),
+            (fmap_select, summary, [('sdc_method', 'distortion_correction')]),
+            (fmapref_buffer, func_fit_reports_wf, [('out', 'inputnode.sdc_boldref')]),
+            (fmapreg_buffer, func_fit_reports_wf, [
+                ('boldref2fmap_xfm', 'inputnode.boldref2fmap_xfm'),
+            ]),
+            (boldref_fmap, func_fit_reports_wf, [('out_file', 'inputnode.fieldmap')]),
+        ])  # fmt:skip
+
+        if not boldref2fmap_xform:
+            config.loggers.workflow.info('Stage 3: Registering fieldmap to boldref')
+            fmapreg_wf = init_coeff2epi_wf(
+                debug='fieldmaps' in config.execution.debug,
+                omp_nthreads=config.nipype.omp_nthreads,
+                sloppy=config.execution.sloppy,
+                name='fmapreg_wf',
+            )
+
+            itk_mat2txt = pe.Node(ConcatenateXFMs(out_fmt='itk'), name='itk_mat2txt')
+            fmapreg_source_files = pe.Node(
+                niu.Merge(2), name='fmapreg_source_files', run_without_submitting=True
+            )
+
+            ds_fmapreg_wf = init_ds_registration_wf(
+                source_file=bold_file,
+                output_dir=config.execution.nibabies_dir,
+                source='boldref',
+                dest=re.sub(r'[^a-zA-Z0-9]', '', fieldmap_id),
+                desc='fmap',
+                name='ds_fmapreg_wf',
+            )
+
+            workflow.connect([
+                (fmap_select, fmapreg_wf, [
+                    ('fmap_ref', 'inputnode.fmap_ref'),
+                    ('fmap_mask', 'inputnode.fmap_mask'),
+                ]),
+                (fmapref_buffer, fmapreg_source_files, [('out', 'in1')]),
+                (fmap_select, fmapreg_source_files, [
+                    ('fmap_ref', 'in2'),
+                ]),
+                (fmapreg_wf, itk_mat2txt, [('outputnode.target2fmap_xfm', 'in_xfms')]),
+                (itk_mat2txt, ds_fmapreg_wf, [('out_xfm', 'inputnode.xform')]),
+                (fmapreg_source_files, ds_fmapreg_wf, [('out', 'inputnode.source_files')]),
+                (ds_fmapreg_wf, fmapreg_buffer, [('outputnode.xform', 'boldref2fmap_xfm')]),
+            ])  # fmt:skip
+        else:
+            config.loggers.workflow.info(
+                'Stage 3: Found fieldmap transform - skipping registration'
+            )
+    else:
+        config.loggers.workflow.info('No fieldmap correction - skipping Stage 3')
+
+    # Stage 4: Create coregistration reference
+    if not coreg_boldref:
+        config.loggers.workflow.info('Stage 4: Adding coregistration boldref workflow')
+
+        # If sbref files are available, add them to the list of sources
+        if sbref_files and nb.load(sbref_files[0]).ndim > 3:
+            raw_sbref_wf = init_raw_boldref_wf(
+                name='raw_sbref_wf',
+                bold_file=sbref_files[0],
+                multiecho=len(sbref_files) > 1,
+                ref_frame_start=config.workflow.hmc_bold_frame,
+            )
+            workflow.connect(raw_sbref_wf, 'outputnode.boldref', fmapref_buffer, 'sbref_files')
+
+        # TODO: Use the premask option here to avoid registering to MNI152NLin2009cAsym
         enhance_boldref_wf = init_enhance_and_skullstrip_bold_wf(omp_nthreads=omp_nthreads)
+        coreg_ref_source_files = pe.Node(
+            niu.Merge(3), name='coreg_ref_source_files', run_without_submitting=True
+        )
 
         ds_coreg_boldref_wf = init_ds_boldref_wf(
-            output_dir=config.execution.output_dir,
+            source_file=bold_file,
+            output_dir=config.execution.nibabies_dir,
             desc='coreg',
             name='ds_coreg_boldref_wf',
         )
+        ds_boldmask_wf = init_ds_boldmask_wf(
+            source_file=bold_file,
+            output_dir=config.execution.nibabies_dir,
+            desc='brain',
+            name='ds_boldmask_wf',
+        )
 
         workflow.connect([
-            (hmcref_buffer, fmapref_buffer, [('boldref', 'boldref_files')]),
             (fmapref_buffer, enhance_boldref_wf, [('out', 'inputnode.in_file')]),
-            (hmc_boldref_source_buffer, ds_coreg_boldref_wf, [
-                ('in_file', 'inputnode.source_files'),
-            ]),
+            (fmapref_buffer, coreg_ref_source_files, [('out', 'in1')]),
+            (coreg_ref_source_files, ds_coreg_boldref_wf, [('out', 'inputnode.source_files')]),
             (ds_coreg_boldref_wf, regref_buffer, [('outputnode.boldref', 'boldref')]),
-            (fmapref_buffer, func_fit_reports_wf, [('out', 'inputnode.sdc_boldref')]),
+            (ds_coreg_boldref_wf, ds_boldmask_wf, [('outputnode.boldref', 'inputnode.source_files')]),
+            (ds_boldmask_wf, regref_buffer, [('outputnode.boldmask', 'boldmask')]),
         ])  # fmt:skip
 
         if fieldmap_id:
-            fmap_select = pe.Node(
-                KeySelect(
-                    fields=['fmap_ref', 'fmap_coeff', 'fmap_mask', 'sdc_method'],
-                    key=fieldmap_id,
+            distortion_params = pe.Node(
+                DistortionParameters(
+                    metadata=metadata,
+                    in_file=bold_file,
+                    # TODO: Estimate fallback readout time
                 ),
-                name='fmap_select',
+                name='distortion_params',
                 run_without_submitting=True,
             )
-
-            if not boldref2fmap_xform:
-                fmapreg_wf = init_coeff2epi_wf(
-                    debug='fieldmaps' in config.execution.debug,
-                    omp_nthreads=config.nipype.omp_nthreads,
-                    sloppy=config.execution.sloppy,
-                    name='fmapreg_wf',
-                )
-
-                itk_mat2txt = pe.Node(ConcatenateXFMs(out_fmt='itk'), name='itk_mat2txt')
-
-                ds_fmapreg_wf = init_ds_registration_wf(
-                    output_dir=config.execution.output_dir,
-                    source='boldref',
-                    dest=fieldmap_id.replace('_', ''),
-                    name='ds_fmapreg_wf',
-                )
-                ds_fmapreg_wf.inputs.inputnode.source_files = [bold_file]
-
-                workflow.connect([
-                    (enhance_boldref_wf, fmapreg_wf, [
-                        ('outputnode.bias_corrected_file', 'inputnode.target_ref'),
-                        ('outputnode.mask_file', 'inputnode.target_mask'),
-                    ]),
-                    (fmap_select, fmapreg_wf, [
-                        ('fmap_ref', 'inputnode.fmap_ref'),
-                        ('fmap_mask', 'inputnode.fmap_mask'),
-                    ]),
-                    (fmapreg_wf, itk_mat2txt, [('outputnode.target2fmap_xfm', 'in_xfms')]),
-                    (itk_mat2txt, ds_fmapreg_wf, [('out_xfm', 'inputnode.xform')]),
-                    (ds_fmapreg_wf, fmapreg_buffer, [('outputnode.xform', 'boldref2fmap_xfm')]),
-                ])  # fmt:skip
-            else:
-                fmapreg_buffer.inputs.boldref2fmap_xfm = boldref2fmap_xform
-
-            unwarp_wf = init_unwarp_wf(
-                free_mem=config.environment.free_mem,
-                debug='fieldmaps' in config.execution.debug,
-                omp_nthreads=config.nipype.omp_nthreads,
+            unwarp_boldref = pe.Node(
+                ResampleSeries(jacobian=jacobian),
+                name='unwarp_boldref',
+                n_procs=omp_nthreads,
+                mem_gb=mem_gb['resampled'],
             )
-            unwarp_wf.inputs.inputnode.metadata = layout.get_metadata(bold_file)
+
+            skullstrip_bold_wf = init_skullstrip_bold_wf()
 
             workflow.connect([
-                (inputnode, fmap_select, [
-                    ('fmap_ref', 'fmap_ref'),
-                    ('fmap_coeff', 'fmap_coeff'),
-                    ('fmap_mask', 'fmap_mask'),
-                    ('sdc_method', 'sdc_method'),
-                    ('fmap_id', 'keys'),
-                ]),
-                (fmap_select, unwarp_wf, [
-                    ('fmap_coeff', 'inputnode.fmap_coeff'),
+                (fmapref_buffer, unwarp_boldref, [('out', 'ref_file')]),
+                (enhance_boldref_wf, unwarp_boldref, [
+                    ('outputnode.bias_corrected_file', 'in_file'),
                 ]),
-                (fmapreg_buffer, unwarp_wf, [
-                    # This looks backwards, but unwarp_wf describes transforms in
-                    # terms of points while we (and init_coeff2epi_wf) describe them
-                    # in terms of images. Mapping fieldmap coordinates into boldref
-                    # coordinates maps the boldref image onto the fieldmap image.
-                    ('boldref2fmap_xfm', 'inputnode.fmap2data_xfm'),
+                (boldref_fmap, unwarp_boldref, [('out_file', 'fieldmap')]),
+                (distortion_params, unwarp_boldref, [
+                    ('readout_time', 'ro_time'),
+                    ('pe_direction', 'pe_dir'),
                 ]),
-                (enhance_boldref_wf, unwarp_wf, [
-                    ('outputnode.bias_corrected_file', 'inputnode.distorted'),
+                (unwarp_boldref, ds_coreg_boldref_wf, [
+                    ('out_file', 'inputnode.boldref'),
                 ]),
-                (unwarp_wf, ds_coreg_boldref_wf, [
-                    ('outputnode.corrected', 'inputnode.boldref'),
+                (ds_coreg_boldref_wf, skullstrip_bold_wf, [
+                    ('outputnode.boldref', 'inputnode.in_file'),
                 ]),
-                (unwarp_wf, regref_buffer, [
-                    ('outputnode.corrected_mask', 'boldmask'),
+                (skullstrip_bold_wf, ds_boldmask_wf, [
+                    ('outputnode.mask_file', 'inputnode.boldmask'),
                 ]),
-                (fmap_select, func_fit_reports_wf, [('fmap_ref', 'inputnode.fmap_ref')]),
-                (fmap_select, summary, [('sdc_method', 'distortion_correction')]),
-                (fmapreg_buffer, func_fit_reports_wf, [
-                    ('boldref2fmap_xfm', 'inputnode.boldref2fmap_xfm'),
-                ]),
-                (unwarp_wf, func_fit_reports_wf, [('outputnode.fieldmap', 'inputnode.fieldmap')]),
+                (fmapreg_buffer, coreg_ref_source_files, [('boldref2fmap_xfm', 'in2')]),
+                (fmap_select, coreg_ref_source_files, [('fmap_coeff', 'in3')]),
             ])  # fmt:skip
+
+            if not boldref2fmap_xform:
+                workflow.connect([
+                    (enhance_boldref_wf, fmapreg_wf, [
+                        ('outputnode.bias_corrected_file', 'inputnode.target_ref'),
+                        ('outputnode.mask_file', 'inputnode.target_mask'),
+                    ]),
+                ])  # fmt:skip
+
         else:
             workflow.connect([
                 (enhance_boldref_wf, ds_coreg_boldref_wf, [
                     ('outputnode.bias_corrected_file', 'inputnode.boldref'),
                 ]),
-                (enhance_boldref_wf, regref_buffer, [
-                    ('outputnode.mask_file', 'boldmask'),
+                (enhance_boldref_wf, ds_boldmask_wf, [
+                    ('outputnode.mask_file', 'inputnode.boldmask'),
                 ]),
             ])  # fmt:skip
     else:
-        config.loggers.workflow.info('Found coregistration reference - skipping Stage 3')
-        regref_buffer.inputs.boldref = precomputed['coreg_boldref']
+        config.loggers.workflow.info('Found coregistration reference - skipping Stage 4')
+
+        # TODO: Allow precomputed bold masks to be passed
+        # Also needs consideration for how it interacts above
+        skullstrip_precomp_ref_wf = init_skullstrip_bold_wf(name='skullstrip_precomp_ref_wf')
+        skullstrip_precomp_ref_wf.inputs.inputnode.in_file = coreg_boldref
+        workflow.connect([
+            (skullstrip_precomp_ref_wf, regref_buffer, [('outputnode.mask_file', 'boldmask')])
+        ])  # fmt:skip
 
     if not boldref2anat_xform:
-        # calculate BOLD registration to anat
+        config.loggers.workflow.info('Stage 5: Adding coregistration workflow')
+        # calculate BOLD registration to T1w
         bold_reg_wf = init_bold_reg_wf(
             bold2anat_dof=config.workflow.bold2anat_dof,
             bold2anat_init=config.workflow.bold2anat_init,
@@ -593,9 +669,11 @@ def init_bold_fit_wf(
         )
 
         ds_boldreg_wf = init_ds_registration_wf(
-            output_dir=config.execution.output_dir,
+            source_file=bold_file,
+            output_dir=config.execution.nibabies_dir,
             source='boldref',
             dest=reference_anat,
+            desc='coreg',
             name='ds_boldreg_wf',
         )
 
@@ -612,11 +690,15 @@ def init_bold_fit_wf(
             (regref_buffer, bold_reg_wf, [('boldref', 'inputnode.ref_bold_brain')]),
             # Incomplete sources
             (regref_buffer, ds_boldreg_wf, [('boldref', 'inputnode.source_files')]),
-            (bold_reg_wf, ds_boldreg_wf, [('outputnode.itk_bold_to_anat', 'inputnode.xform')]),
+            (bold_reg_wf, ds_boldreg_wf, [
+                ('outputnode.itk_bold_to_anat', 'inputnode.xform'),
+                ('outputnode.metadata', 'inputnode.metadata'),
+            ]),
             (ds_boldreg_wf, outputnode, [('outputnode.xform', 'boldref2anat_xfm')]),
             (bold_reg_wf, summary, [('outputnode.fallback', 'fallback')]),
         ])  # fmt:skip
     else:
+        config.loggers.workflow.info('Found coregistration transform - skipping Stage 5')
         outputnode.inputs.boldref2anat_xfm = boldref2anat_xform
 
     return workflow
diff --git a/nibabies/workflows/bold/hmc.py b/nibabies/workflows/bold/hmc.py
index 658fad6b..171299e5 100644
--- a/nibabies/workflows/bold/hmc.py
+++ b/nibabies/workflows/bold/hmc.py
@@ -32,8 +32,6 @@
 from nipype.interfaces import utility as niu
 from nipype.pipeline import engine as pe
 
-from ...config import DEFAULT_MEMORY_MIN_GB
-
 
 def init_bold_hmc_wf(mem_gb: float, omp_nthreads: int, name: str = 'bold_hmc_wf'):
     """
@@ -73,14 +71,9 @@ def init_bold_hmc_wf(mem_gb: float, omp_nthreads: int, name: str = 'bold_hmc_wf'
     -------
     xforms
         ITKTransform file aligning each volume to ``ref_image``
-    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]_.
 
     """
     from niworkflows.engine.workflows import LiterateWorkflow as Workflow
-    from niworkflows.interfaces.confounds import NormalizeMotionParams
     from niworkflows.interfaces.itk import MCFLIRT2ITK
 
     workflow = Workflow(name=name)
@@ -94,36 +87,24 @@ def init_bold_hmc_wf(mem_gb: float, omp_nthreads: int, name: str = 'bold_hmc_wf'
     inputnode = pe.Node(
         niu.IdentityInterface(fields=['bold_file', 'raw_ref_image']), name='inputnode'
     )
-    outputnode = pe.Node(
-        niu.IdentityInterface(fields=['xforms', 'movpar_file', 'rmsd_file']), name='outputnode'
-    )
+    outputnode = pe.Node(niu.IdentityInterface(fields=['xforms']), name='outputnode')
 
     # Head motion correction (hmc)
     mcflirt = pe.Node(
-        fsl.MCFLIRT(save_mats=True, save_plots=True, save_rms=True),
+        fsl.MCFLIRT(save_mats=True),
         name='mcflirt',
         mem_gb=mem_gb * 3,
     )
 
     fsl2itk = pe.Node(MCFLIRT2ITK(), name='fsl2itk', mem_gb=0.05, n_procs=omp_nthreads)
 
-    normalize_motion = pe.Node(
-        NormalizeMotionParams(format='FSL'), name='normalize_motion', mem_gb=DEFAULT_MEMORY_MIN_GB
-    )
-
-    def _pick_rel(rms_files):
-        return rms_files[-1]
-
     workflow.connect([
         (inputnode, mcflirt, [('raw_ref_image', 'ref_file'),
                               ('bold_file', 'in_file')]),
         (inputnode, fsl2itk, [('raw_ref_image', 'in_source'),
                               ('raw_ref_image', 'in_reference')]),
         (mcflirt, fsl2itk, [('mat_file', 'in_files')]),
-        (mcflirt, normalize_motion, [('par_file', 'in_file')]),
-        (mcflirt, outputnode, [(('rms_files', _pick_rel), 'rmsd_file')]),
         (fsl2itk, outputnode, [('out_file', 'xforms')]),
-        (normalize_motion, outputnode, [('out_file', 'movpar_file')]),
     ])  # fmt:skip
 
     return workflow
diff --git a/nibabies/workflows/bold/outputs.py b/nibabies/workflows/bold/outputs.py
index 8d051f02..65e3bf8f 100644
--- a/nibabies/workflows/bold/outputs.py
+++ b/nibabies/workflows/bold/outputs.py
@@ -411,6 +411,7 @@ def init_func_fit_reports_wf(
 
 def init_ds_boldref_wf(
     *,
+    source_file: str,
     output_dir,
     desc: str,
     name='ds_boldref_wf',
@@ -427,13 +428,14 @@ def init_ds_boldref_wf(
         BIDSURI(
             numinputs=1,
             dataset_links=config.execution.dataset_links,
-            out_dir=str(config.execution.output_dir.absolute()),
+            out_dir=str(output_dir),
         ),
         name='sources',
     )
 
     ds_boldref = pe.Node(
         DerivativesDataSink(
+            source_file=source_file,
             base_directory=output_dir,
             desc=desc,
             suffix='boldref',
@@ -444,30 +446,29 @@ def init_ds_boldref_wf(
         run_without_submitting=True,
     )
 
-    # fmt:off
     workflow.connect([
         (inputnode, sources, [('source_files', 'in1')]),
-        (inputnode, ds_boldref, [('boldref', 'in_file'),
-                                 ('source_files', 'source_file')]),
+        (inputnode, ds_boldref, [('boldref', 'in_file')]),
         (sources, ds_boldref, [('out', 'Sources')]),
         (ds_boldref, outputnode, [('out_file', 'boldref')]),
-    ])
-    # fmt:on
+    ])  # fmt:skip
 
     return workflow
 
 
 def init_ds_registration_wf(
     *,
+    source_file: str,
     output_dir: str,
     source: str,
     dest: str,
     name: str,
+    desc: str | None = None,
 ) -> pe.Workflow:
     workflow = pe.Workflow(name=name)
 
     inputnode = pe.Node(
-        niu.IdentityInterface(fields=['source_files', 'xform']),
+        niu.IdentityInterface(fields=['source_files', 'xform', 'metadata']),
         name='inputnode',
     )
     outputnode = pe.Node(niu.IdentityInterface(fields=['xform']), name='outputnode')
@@ -476,15 +477,17 @@ def init_ds_registration_wf(
         BIDSURI(
             numinputs=1,
             dataset_links=config.execution.dataset_links,
-            out_dir=str(config.execution.output_dir.absolute()),
+            out_dir=str(output_dir),
         ),
         name='sources',
     )
 
     ds_xform = pe.Node(
         DerivativesDataSink(
+            source_file=source_file,
             base_directory=output_dir,
             mode='image',
+            desc=desc,
             suffix='xfm',
             extension='.txt',
             dismiss_entities=dismiss_entities(['part']),
@@ -495,21 +498,22 @@ def init_ds_registration_wf(
         mem_gb=DEFAULT_MEMORY_MIN_GB,
     )
 
-    # fmt:off
     workflow.connect([
         (inputnode, sources, [('source_files', 'in1')]),
-        (inputnode, ds_xform, [('xform', 'in_file'),
-                               ('source_files', 'source_file')]),
+        (inputnode, ds_xform, [
+            ('xform', 'in_file'),
+            ('metadata', 'meta_dict'),
+        ]),
         (sources, ds_xform, [('out', 'Sources')]),
         (ds_xform, outputnode, [('out_file', 'xform')]),
-    ])
-    # fmt:on
+    ])  # fmt:skip
 
     return workflow
 
 
 def init_ds_hmc_wf(
     *,
+    source_file: str,
     output_dir,
     name='ds_hmc_wf',
 ) -> pe.Workflow:
@@ -525,13 +529,14 @@ def init_ds_hmc_wf(
         BIDSURI(
             numinputs=1,
             dataset_links=config.execution.dataset_links,
-            out_dir=str(config.execution.output_dir.absolute()),
+            out_dir=str(output_dir),
         ),
         name='sources',
     )
 
     ds_xforms = pe.Node(
         DerivativesDataSink(
+            source_file=source_file,
             base_directory=output_dir,
             desc='hmc',
             suffix='xfm',
@@ -544,15 +549,12 @@ def init_ds_hmc_wf(
         run_without_submitting=True,
     )
 
-    # fmt:off
     workflow.connect([
         (inputnode, sources, [('source_files', 'in1')]),
-        (inputnode, ds_xforms, [('xforms', 'in_file'),
-                                ('source_files', 'source_file')]),
+        (inputnode, ds_xforms, [('xforms', 'in_file')]),
         (sources, ds_xforms, [('out', 'Sources')]),
         (ds_xforms, outputnode, [('out_file', 'xforms')]),
-    ])
-    # fmt:on
+    ])  # fmt:skip
 
     return workflow
 
@@ -964,3 +966,53 @@ def init_bold_preproc_report_wf(
     # fmt:on
 
     return workflow
+
+
+def init_ds_boldmask_wf(
+    *,
+    source_file: str,
+    output_dir,
+    desc: str,
+    name='ds_boldmask_wf',
+) -> pe.Workflow:
+    """Write out a BOLD mask."""
+    workflow = pe.Workflow(name=name)
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(fields=['source_files', 'boldmask']),
+        name='inputnode',
+    )
+    outputnode = pe.Node(niu.IdentityInterface(fields=['boldmask']), name='outputnode')
+
+    sources = pe.Node(
+        BIDSURI(
+            numinputs=1,
+            dataset_links=config.execution.dataset_links,
+            out_dir=str(output_dir),
+        ),
+        name='sources',
+    )
+
+    ds_boldmask = pe.Node(
+        DerivativesDataSink(
+            source_file=source_file,
+            base_directory=output_dir,
+            desc=desc,
+            suffix='mask',
+            compress=True,
+            dismiss_entities=DEFAULT_DISMISS_ENTITIES,
+        ),
+        name='ds_boldmask',
+        run_without_submitting=True,
+    )
+
+    workflow.connect([
+        (inputnode, sources, [('source_files', 'in1')]),
+        (inputnode, ds_boldmask, [
+            ('boldmask', 'in_file'),
+        ]),
+        (sources, ds_boldmask, [('out', 'Sources')]),
+        (ds_boldmask, outputnode, [('out_file', 'boldmask')]),
+    ])  # fmt:skip
+
+    return workflow
diff --git a/nibabies/workflows/bold/reference.py b/nibabies/workflows/bold/reference.py
index 6066e765..cbf542ef 100644
--- a/nibabies/workflows/bold/reference.py
+++ b/nibabies/workflows/bold/reference.py
@@ -83,7 +83,6 @@ def init_raw_boldref_wf(
         beginning of ``bold_file``
 
     """
-    from niworkflows.interfaces.bold import NonsteadyStatesDetector
     from niworkflows.interfaces.images import RobustAverage
 
     workflow = Workflow(name=name)
@@ -113,11 +112,7 @@ def init_raw_boldref_wf(
     if bold_file is not None:
         inputnode.inputs.bold_file = bold_file
 
-    val_bold = pe.Node(
-        ValidateImage(),
-        name='val_bold',
-        mem_gb=DEFAULT_MEMORY_MIN_GB,
-    )
+    validation_and_dummies_wf = init_validation_and_dummies_wf()
 
     # Drop frames to avoid startle when MRI begins acquiring
     select_frames = pe.Node(
@@ -126,34 +121,29 @@ def init_raw_boldref_wf(
     )
     select_frames.inputs.ref_frame_start = ref_frame_start
 
-    get_dummy = pe.Node(NonsteadyStatesDetector(), name='get_dummy')
     gen_avg = pe.Node(RobustAverage(), name='gen_avg', mem_gb=1)
 
-    calc_dummy_scans = pe.Node(
-        niu.Function(function=pass_dummy_scans, output_names=['skip_vols_num']),
-        name='calc_dummy_scans',
-        run_without_submitting=True,
-        mem_gb=DEFAULT_MEMORY_MIN_GB,
-    )
-
     workflow.connect([
-        (inputnode, val_bold, [('bold_file', 'in_file')]),
-        (inputnode, get_dummy, [('bold_file', 'in_file')]),
-        (inputnode, calc_dummy_scans, [('dummy_scans', 'dummy_scans')]),
-        (val_bold, gen_avg, [('out_file', 'in_file')]),
-        (val_bold, select_frames, [('out_file', 'in_file')]),
+        (inputnode, validation_and_dummies_wf, [
+            ('bold_file', 'inputnode.bold_file'),
+            ('dummy_scans', 'inputnode.dummy_scans'),
+        ]),
+        (validation_and_dummies_wf, outputnode, [
+            ('outputnode.bold_file', 'bold_file'),
+            ('outputnode.skip_vols', 'skip_vols'),
+            ('outputnode.algo_dummy_scans', 'algo_dummy_scans'),
+            ('outputnode.validation_report', 'validation_report'),
+        ]),
+        (validation_and_dummies_wf, gen_avg, [
+            ('outputnode.bold_file', 'in_file'),
+        ]),
+        (validation_and_dummies_wf, select_frames, [
+            ('outputnode.bold_file', 'in_file'),
+        ]),
         (inputnode, select_frames, [('dummy_scans', 'dummy_scans')]),
         (select_frames, gen_avg, [('t_mask', 't_mask')]),
-        (get_dummy, calc_dummy_scans, [('n_dummy', 'algo_dummy_scans')]),
-        (val_bold, outputnode, [
-            ('out_file', 'bold_file'),
-            ('out_report', 'validation_report'),
-        ]),
-        (calc_dummy_scans, outputnode, [('skip_vols_num', 'skip_vols')]),
         (gen_avg, outputnode, [('out_file', 'boldref')]),
-        (get_dummy, outputnode, [('n_dummy', 'algo_dummy_scans')]),
     ])  # fmt:skip
-
     return workflow
 
 
@@ -183,3 +173,102 @@ def _select_frames(
     t_mask = np.array([False] * img_len, dtype=bool)
     t_mask[start_frame:] = True
     return start_frame, list(t_mask)
+
+
+def init_validation_and_dummies_wf(
+    bold_file: str | None = None,
+    name: str = 'validation_and_dummies_wf',
+):
+    """
+    Build a workflow that validates a BOLD image and detects non-steady-state volumes.
+
+    Workflow Graph
+        .. workflow::
+            :graph2use: orig
+            :simple_form: yes
+
+            from fmriprep.workflows.bold.reference import init_validation_and_dummies_wf
+            wf = init_validation_and_dummies_wf()
+
+    Parameters
+    ----------
+    bold_file : :obj:`str`
+        BOLD series NIfTI file
+    name : :obj:`str`
+        Name of workflow (default: ``validation_and_dummies_wf``)
+
+    Inputs
+    ------
+    bold_file : str
+        BOLD series NIfTI file
+    dummy_scans : int or None
+        Number of non-steady-state volumes specified by user at beginning of ``bold_file``
+
+    Outputs
+    -------
+    bold_file : str
+        Validated BOLD series NIfTI file
+    skip_vols : int
+        Number of non-steady-state volumes selected at beginning of ``bold_file``
+    algo_dummy_scans : int
+        Number of non-steady-state volumes agorithmically detected at
+        beginning of ``bold_file``
+
+    """
+    from niworkflows.interfaces.bold import NonsteadyStatesDetector
+
+    workflow = Workflow(name=name)
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(fields=['bold_file', 'dummy_scans']),
+        name='inputnode',
+    )
+    outputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=[
+                'bold_file',
+                'skip_vols',
+                'algo_dummy_scans',
+                't_mask',
+                'validation_report',
+            ]
+        ),
+        name='outputnode',
+    )
+
+    # Simplify manually setting input image
+    if bold_file is not None:
+        inputnode.inputs.bold_file = bold_file
+
+    val_bold = pe.Node(
+        ValidateImage(),
+        name='val_bold',
+        mem_gb=DEFAULT_MEMORY_MIN_GB,
+    )
+
+    get_dummy = pe.Node(NonsteadyStatesDetector(), name='get_dummy')
+
+    calc_dummy_scans = pe.Node(
+        niu.Function(function=pass_dummy_scans, output_names=['skip_vols_num']),
+        name='calc_dummy_scans',
+        run_without_submitting=True,
+        mem_gb=DEFAULT_MEMORY_MIN_GB,
+    )
+
+    workflow.connect([
+        (inputnode, val_bold, [('bold_file', 'in_file')]),
+        (val_bold, outputnode, [
+            ('out_file', 'bold_file'),
+            ('out_report', 'validation_report'),
+        ]),
+        (inputnode, get_dummy, [('bold_file', 'in_file')]),
+        (inputnode, calc_dummy_scans, [('dummy_scans', 'dummy_scans')]),
+        (get_dummy, calc_dummy_scans, [('n_dummy', 'algo_dummy_scans')]),
+        (get_dummy, outputnode, [
+            ('n_dummy', 'algo_dummy_scans'),
+            ('t_mask', 't_mask'),
+        ]),
+        (calc_dummy_scans, outputnode, [('skip_vols_num', 'skip_vols')]),
+    ])  # fmt:skip
+
+    return workflow
diff --git a/nibabies/workflows/bold/registration.py b/nibabies/workflows/bold/registration.py
index e3e788a5..12cdc8a3 100644
--- a/nibabies/workflows/bold/registration.py
+++ b/nibabies/workflows/bold/registration.py
@@ -15,7 +15,7 @@
 import os
 import typing as ty
 
-from nipype.interfaces import c3, fsl
+from nipype.interfaces import fsl
 from nipype.interfaces import utility as niu
 from nipype.pipeline import engine as pe
 
@@ -107,6 +107,8 @@ def init_bold_reg_wf(
         Affine transform from anatomical space to BOLD space (ITK format)
     fallback
         Boolean indicating whether BBR was rejected (mri_coreg registration returned)
+    metadata
+        Output metadata from the registration workflow
 
     See Also
     --------
@@ -133,7 +135,9 @@ def init_bold_reg_wf(
     )
 
     outputnode = pe.Node(
-        niu.IdentityInterface(fields=['itk_bold_to_anat', 'itk_anat_to_bold', 'fallback']),
+        niu.IdentityInterface(
+            fields=['itk_bold_to_anat', 'itk_anat_to_bold', 'fallback', 'metadata']
+        ),
         name='outputnode',
     )
 
@@ -168,6 +172,7 @@ def init_bold_reg_wf(
             ('outputnode.itk_bold_to_anat', 'itk_bold_to_anat'),
             ('outputnode.itk_anat_to_bold', 'itk_anat_to_bold'),
             ('outputnode.fallback', 'fallback'),
+            ('outputnode.metadata', 'metadata'),
         ]),
     ])  # fmt:skip
 
@@ -249,12 +254,15 @@ def init_bbreg_wf(
         Affine transform from anatomical space to BOLD space (ITK format)
     fallback
         Boolean indicating whether BBR was rejected (mri_coreg registration returned)
+    metadata
+        Output metadata from the registration workflow
 
     """
     from nipype.interfaces.freesurfer import BBRegister
     from niworkflows.engine.workflows import LiterateWorkflow as Workflow
     from niworkflows.interfaces.nitransforms import ConcatenateXFMs
     from niworkflows.interfaces.patches import FreeSurferSource
+    from niworkflows.interfaces.utility import DictMerge
 
     from nibabies.interfaces.patches import MRICoreg
 
@@ -291,7 +299,7 @@ def init_bbreg_wf(
         name='inputnode',
     )
     outputnode = pe.Node(
-        niu.IdentityInterface(['itk_bold_to_anat', 'itk_anat_to_bold', 'fallback']),
+        niu.IdentityInterface(['itk_bold_to_anat', 'itk_anat_to_bold', 'fallback', 'metadata']),
         name='outputnode',
     )
 
@@ -339,6 +347,12 @@ def init_bbreg_wf(
     merge_ltas = pe.Node(niu.Merge(2), name='merge_ltas', run_without_submitting=True)
     concat_xfm = pe.Node(ConcatenateXFMs(inverse=True), name='concat_xfm')
 
+    # Set up GeneratedBy metadata and add a merge node for cost, if available
+    gen_by = pe.Node(niu.Merge(2), run_without_submitting=True, name='gen_by')
+    select_gen = pe.Node(niu.Select(index=0), run_without_submitting=True, name='select_gen')
+    metadata = pe.Node(niu.Merge(2), run_without_submitting=True, name='metadata')
+    merge_meta = pe.Node(DictMerge(), run_without_submitting=True, name='merge_meta')
+
     workflow.connect([
         (inputnode, merge_ltas, [('fsnative2anat_xfm', 'in2')]),
         # Wire up the co-registration alternatives
@@ -347,10 +361,20 @@ def init_bbreg_wf(
         (merge_ltas, concat_xfm, [('out', 'in_xfms')]),
         (concat_xfm, outputnode, [('out_xfm', 'itk_bold_to_anat')]),
         (concat_xfm, outputnode, [('out_inv', 'itk_anat_to_bold')]),
+        # Wire up the metadata alternatives
+        (gen_by, select_gen, [('out', 'inlist')]),
+        (select_gen, metadata, [('out', 'in1')]),
+        (metadata, merge_meta, [('out', 'in_dicts')]),
+        (merge_meta, outputnode, [('out_dict', 'metadata')]),
     ])  # fmt:skip
 
     # Do not initialize with header, use mri_coreg
     if bold2anat_init != 'header':
+        gen_by.inputs.in2 = {
+            'GeneratedBy': [
+                {'Name': 'mri_coreg', 'Version': mri_coreg.interface.version or '<unknown>'}
+            ]
+        }
         workflow.connect([
             (inputnode, mri_coreg, [('subjects_dir', 'subjects_dir'),
                                     ('subject_id', 'subject_id'),
@@ -382,6 +406,26 @@ def init_bbreg_wf(
         (bbregister, transforms, [('out_lta_file', 'in1')]),
     ])  # fmt:skip
 
+    gen_by.inputs.in1 = {
+        'GeneratedBy': [
+            {'Name': 'bbregister', 'Version': bbregister.interface.version or '<unknown>'}
+        ]
+    }
+
+    costs = pe.Node(niu.Merge(2), run_without_submitting=True, name='costs')
+    select_cost = pe.Node(niu.Select(index=0), run_without_submitting=True, name='select_cost')
+    read_cost = pe.Node(niu.Function(function=_read_cost), name='read_cost')
+
+    workflow.connect([
+        (bbregister, costs, [
+            ('min_cost_file', 'in1'),
+            ('init_cost_file', 'in2'),
+        ]),
+        (costs, select_cost, [('out', 'inlist')]),
+        (select_cost, read_cost, [('out', 'cost_file')]),
+        (read_cost, metadata, [('out', 'in2')]),
+    ])  # fmt:skip
+
     # Short-circuit workflow building, use boundary-based registration
     if use_bbr is True:
         outputnode.inputs.fallback = False
@@ -475,12 +519,15 @@ def init_fsl_bbr_wf(
         Affine transform from anatomical space to BOLD space (ITK format)
     fallback
         Boolean indicating whether BBR was rejected (rigid FLIRT registration returned)
+    metadata
+        Output metadata from the registration workflow
 
     """
     from nipype.interfaces.freesurfer import MRICoreg
     from niworkflows.engine.workflows import LiterateWorkflow as Workflow
     from niworkflows.interfaces.freesurfer import PatchedLTAConvert as LTAConvert
     from niworkflows.interfaces.nibabel import ApplyMask
+    from niworkflows.interfaces.nitransforms import ConvertAffine
     from niworkflows.utils.images import dseg_label as _dseg_label
 
     workflow = Workflow(name=name)
@@ -514,7 +561,7 @@ def init_fsl_bbr_wf(
         name='inputnode',
     )
     outputnode = pe.Node(
-        niu.IdentityInterface(['itk_bold_to_anat', 'itk_anat_to_bold', 'fallback']),
+        niu.IdentityInterface(['itk_bold_to_anat', 'itk_anat_to_bold', 'fallback', 'metadata']),
         name='outputnode',
     )
 
@@ -533,6 +580,9 @@ def init_fsl_bbr_wf(
             'T2w intermediate for FSL is not implemented, registering with T1w instead.'
         )
 
+    metadata = pe.Node(niu.Merge(2), run_without_submitting=True, name='metadata')
+    select_meta = pe.Node(niu.Select(index=0), run_without_submitting=True, name='select_meta')
+
     # Mask T1w_preproc with T1w_mask to make T1w_brain
     mask_anat_brain = pe.Node(ApplyMask(), name='mask_anat_brain')
 
@@ -543,52 +593,47 @@ def init_fsl_bbr_wf(
         mem_gb=5,
     )
 
-    lta_to_fsl = pe.Node(LTAConvert(out_fsl=True), name='lta_to_fsl', mem_gb=DEFAULT_MEMORY_MIN_GB)
-
-    invt_bbr = pe.Node(
-        fsl.ConvertXFM(invert_xfm=True), name='invt_bbr', mem_gb=DEFAULT_MEMORY_MIN_GB
-    )
-
-    # BOLD to anat transform matrix is from fsl, using c3 tools to convert to
-    # something ANTs will like.
-    fsl2itk_fwd = pe.Node(
-        c3.C3dAffineTool(fsl2ras=True, itk_transform=True),
-        name='fsl2itk_fwd',
+    xfm2itk = pe.Node(
+        ConvertAffine(in_fmt='fsl', out_fmt='itk', inverse=True),
+        name='xfm2itk',
         mem_gb=DEFAULT_MEMORY_MIN_GB,
     )
-    fsl2itk_inv = pe.Node(
-        c3.C3dAffineTool(fsl2ras=True, itk_transform=True),
-        name='fsl2itk_inv',
-        mem_gb=DEFAULT_MEMORY_MIN_GB,
-    )
-    # fmt:off
+
+    metadata.inputs.in2 = {
+        'GeneratedBy': [
+            {'Name': 'mri_coreg', 'Version': mri_coreg.interface.version or '<unknown>'}
+        ]
+    }
+
     workflow.connect([
-        (inputnode, mask_anat_brain, [('anat_preproc', 'in_file'),
-                                      ('anat_mask', 'in_mask')]),
+        (inputnode, mask_anat_brain, [
+            ('anat_preproc', 'in_file'),
+            ('anat_mask', 'in_mask'),
+        ]),
         (inputnode, mri_coreg, [('in_file', 'source_file')]),
-        (inputnode, fsl2itk_fwd, [('in_file', 'source_file')]),
-        (inputnode, fsl2itk_inv, [('in_file', 'reference_file')]),
+        (inputnode, xfm2itk, [('in_file', 'moving')]),
         (mask_anat_brain, mri_coreg, [('out_file', 'reference_file')]),
-        (mask_anat_brain, fsl2itk_fwd, [('out_file', 'reference_file')]),
-        (mask_anat_brain, fsl2itk_inv, [('out_file', 'source_file')]),
-        (mri_coreg, lta_to_fsl, [('out_lta_file', 'in_lta')]),
-        (invt_bbr, fsl2itk_inv, [('out_file', 'transform_file')]),
-        (fsl2itk_fwd, outputnode, [('itk_transform', 'itk_bold_to_anat')]),
-        (fsl2itk_inv, outputnode, [('itk_transform', 'itk_anat_to_bold')]),
+        (mask_anat_brain, xfm2itk, [('out_file', 'reference')]),
+        (xfm2itk, outputnode, [
+            ('out_xfm', 'itk_bold_to_anat'),
+            ('out_inv', 'itk_anat_to_bold'),
+        ]),
+        # Wire up the metadata alternatives
+        (metadata, select_meta, [('out', 'inlist')]),
+        (select_meta, outputnode, [('out', 'metadata')]),
     ])  # fmt:skip
 
     # Short-circuit workflow building, use rigid registration
     if use_bbr is False:
-        # fmt:off
-        workflow.connect([
-            (lta_to_fsl, invt_bbr, [('out_fsl', 'in_file')]),
-            (lta_to_fsl, fsl2itk_fwd, [('out_fsl', 'transform_file')]),
-        ])
-        # fmt:on
+        xfm2itk.inputs.in_fmt = 'fs'  # Override
+        workflow.connect(mri_coreg, 'out_lta_file', xfm2itk, 'in_xfm')
+
         outputnode.inputs.fallback = True
 
         return workflow
 
+    lta_to_fsl = pe.Node(LTAConvert(out_fsl=True), name='lta_to_fsl', mem_gb=DEFAULT_MEMORY_MIN_GB)
+
     flt_bbr = pe.Node(
         fsl.FLIRT(cost_func='bbr', dof=bold2anat_dof, args='-basescale 1'),
         name='flt_bbr',
@@ -601,13 +646,18 @@ def init_fsl_bbr_wf(
         # Should mostly be hit while building docs
         LOGGER.warning('FSLDIR unset - using packaged BBR schedule')
         flt_bbr.inputs.schedule = data.load('flirtsch/bbr.sch')
-    # fmt:off
+
+    metadata.inputs.in1 = {
+        'GeneratedBy': [{'Name': 'flirt', 'Version': flt_bbr.interface.version or '<unknown>'}]
+    }
+
     workflow.connect([
         (inputnode, wm_mask, [('anat_dseg', 'in_seg')]),
         (inputnode, flt_bbr, [('in_file', 'in_file')]),
+        (mri_coreg, lta_to_fsl, [('out_lta_file', 'in_lta')]),
         (lta_to_fsl, flt_bbr, [('out_fsl', 'in_matrix_file')]),
-    ])
-    # fmt:on
+    ])  # fmt:skip
+
     if sloppy is True:
         downsample = pe.Node(
             niu.Function(
@@ -615,35 +665,25 @@ def init_fsl_bbr_wf(
             ),
             name='downsample',
         )
-        # fmt:off
         workflow.connect([
             (mask_anat_brain, downsample, [('out_file', 'in_file')]),
             (wm_mask, downsample, [('out', 'in_mask')]),
             (downsample, flt_bbr, [('out_file', 'reference'),
                                    ('out_mask', 'wm_seg')]),
-        ])
-        # fmt:on
+        ])  # fmt:skip
     else:
-        # fmt:off
         workflow.connect([
             (mask_anat_brain, flt_bbr, [('out_file', 'reference')]),
             (wm_mask, flt_bbr, [('out', 'wm_seg')]),
-        ])
-        # fmt:on
+        ])  # fmt:skip
 
     # Short-circuit workflow building, use boundary-based registration
     if use_bbr is True:
-        # fmt:off
-        workflow.connect([
-            (flt_bbr, invt_bbr, [('out_matrix_file', 'in_file')]),
-            (flt_bbr, fsl2itk_fwd, [('out_matrix_file', 'transform_file')]),
-        ])
-        # fmt:on
+        workflow.connect(flt_bbr, 'out_matrix_file', xfm2itk, 'in_xfm')
         outputnode.inputs.fallback = False
 
         return workflow
 
-    # Reached only if use_bbr is None
     transforms = pe.Node(niu.Merge(2), run_without_submitting=True, name='transforms')
 
     compare_transforms = pe.Node(niu.Function(function=compare_xforms), name='compare_transforms')
@@ -651,7 +691,7 @@ def init_fsl_bbr_wf(
     select_transform = pe.Node(niu.Select(), run_without_submitting=True, name='select_transform')
 
     fsl_to_lta = pe.MapNode(LTAConvert(out_lta=True), iterfield=['in_fsl'], name='fsl_to_lta')
-    # fmt:off
+
     workflow.connect([
         (flt_bbr, transforms, [('out_matrix_file', 'in1')]),
         (lta_to_fsl, transforms, [('out_fsl', 'in2')]),
@@ -664,10 +704,10 @@ def init_fsl_bbr_wf(
         # Select output transform
         (transforms, select_transform, [('out', 'inlist')]),
         (compare_transforms, select_transform, [('out', 'index')]),
-        (select_transform, invt_bbr, [('out', 'in_file')]),
-        (select_transform, fsl2itk_fwd, [('out', 'transform_file')]),
-    ])
-    # fmt:on
+        (select_transform, xfm2itk, [('out', 'in_xfm')]),
+        # Select metadata
+        (compare_transforms, select_meta, [('out', 'index')]),
+    ])  # fmt:skip
 
     return workflow
 
@@ -757,3 +797,10 @@ def _conditional_downsampling(in_file, in_mask, zoom_th=4.0):
     nb.Nifti1Image(newmaskdata, newmask.affine, hdr).to_filename(out_mask)
 
     return str(out_file), str(out_mask)
+
+
+def _read_cost(cost_file) -> dict[str, float]:
+    """Read a cost from a file."""
+    # Cost file contains mincost, WM intensity, Ctx intensity, Pct Contrast
+    with open(cost_file) as fobj:
+        return {'FinalCost': float(fobj.read().split()[0])}
diff --git a/pyproject.toml b/pyproject.toml
index 738bcc7b..949a1efa 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -25,7 +25,8 @@ dependencies = [
     "nireports >= 25.3.0",
     "nitime",
     "nitransforms >= 25.0.1",
-    "niworkflows >= 1.14.1",
+    #"niworkflows >= 1.14.1",
+    "niworkflows @ git+https://github.com/nipreps/niworkflows.git@add/nitransforms-convertaffine",
     "numpy >= 2.0",
     "packaging",
     "pandas < 3",