diff --git a/docs/usage.rst b/docs/usage.rst
index 853e736f..559692d3 100644
--- a/docs/usage.rst
+++ b/docs/usage.rst
@@ -211,16 +211,37 @@ Examples: ::
     $ petprep /data/bids_root /out participant --hmc-init-frame 10 --hmc-init-frame-fix
     $ petprep /data/bids_root /out participant --hmc-off
 
+
+Reference image selection
+-------------------------
+Use :option:`--petref` to control how the reference volume is built from the
+dynamic PET series. Each strategy uses the frame timing metadata from
+``FrameTimesStart`` and ``FrameDuration`` to weight volumes; missing metadata
+will raise an error before preprocessing starts.
+
+* ``template`` (default) reuses the motion-correction template, providing a
+  consistent target for downstream registration. When :option:`--hmc-off`
+  disables motion correction, requesting ``template`` automatically falls back
+  to ``twa`` with a warning.
+* ``twa`` computes a time-weighted average, which often emphasizes later frames with
+  higher counts and longer durations.
+* ``sum`` produces a straightforward summed image.
+* ``first5min`` averages only the first 5 minutes of PET data to capture perfusion-like
+  uptake.
+
 Anatomical co-registration
 --------------------------
 *PETPrep* aligns the PET reference volume to the T1-weighted anatomy before
-deriving downstream outputs. By default, FreeSurfer's ``mri_coreg`` performs
-the alignment, with the :option:`--pet2anat-dof` flag controlling the degrees
-of freedom (rigid-body, 6 dof, is the default). When working with low
-signal-to-noise references or challenging anatomy, the
-:option:`--pet2anat-robust` flag enables ``mri_robust_register`` with an NMI
-cost function to improve robustness. This mode is restricted to rigid-body
-alignment and therefore requires ``--pet2anat-dof 6``.
+deriving downstream outputs. The anatomical image is first trimmed with
+FSL's ``robustfov`` to remove the shoulder/neck and masked to limit registration to brain voxels. Choose
+the registration backend with :option:`--pet2anat-method`: ``mri_coreg``
+(default FreeSurfer co-registration), ``robust`` (FreeSurfer
+``mri_robust_register`` with an NMI cost function), or ``ants`` (ANTs rigid
+registration that consumes the unmasked T1w and a separate mask). The
+:option:`--pet2anat-dof` flag controls the degrees of freedom; ``robust`` and
+``ants`` are limited to rigid-body alignment and therefore require
+``--pet2anat-dof 6``. All modes emit paired ITK transforms for reuse in later
+resampling steps.
 
 Segmentation
 ----------------
diff --git a/docs/workflows.rst b/docs/workflows.rst
index 12d6c83f..d6c66e48 100644
--- a/docs/workflows.rst
+++ b/docs/workflows.rst
@@ -342,29 +342,17 @@ PET reference image estimation
     wf = init_raw_petref_wf()
 
 This workflow estimates a reference image for a
-:abbr:`PET (positron emission tomography)` series as follows:
-When T1-saturation effects ("dummy scans" or non-steady state volumes) are
-detected, they are averaged and used as reference due to their
-superior tissue contrast.
-Otherwise, a median of motion corrected subset of volumes is used.
-
-This reference is used for :ref:`head-motion estimation <pet_hmc>`.
-
-For the :ref:`registration workflow <pet_reg>`, the reference image is
-either the above described reference image or a single-band reference,
-if one is found in the input dataset.
-In either case, this image is contrast-enhanced and skull-stripped
-(see :py:func:`~niworkflows.func.util.init_enhance_and_skullstrip_bold_wf`).
-If fieldmaps are present, the skull-stripped reference is corrected
-prior to registration.
-
-.. figure:: _static/sub-01_task-balloonanalogrisktask_run-1_desc-rois_bold.svg
-
-    The red contour shows the brain mask estimated for a PET reference volume.
-    The blue and magenta contours show the tCompCor and aCompCor masks,
-    respectively. (See :ref:`pet_confounds`, below.)
-
-.. _pet_hmc:
+:abbr:`PET (positron emission tomography)` series according to the strategy
+requested with :option:`--petref`:
+
+* ``template`` (default) uses the motion-correction template built during
+  head-motion estimation. When head-motion correction is disabled, ``template``
+  automatically falls back to ``twa`` while emitting a warning.
+* ``twa`` computes a time-weighted average of the motion-corrected series,
+  preserving dynamic information when later frames carry more counts.
+* ``sum`` produces a summed image of the motion-corrected series.
+* ``first5min`` weights only the first five minutes of the acquisition, which
+  can be helpful for tracers whose early dynamics resemble perfusion.
 
 Head-motion estimation
 ~~~~~~~~~~~~~~~~~~~~~~
@@ -438,17 +426,19 @@ PET to T1w registration
         omp_nthreads=1,
         pet2anat_dof=6,
         mem_gb=3,
-        use_robust_register=False,
+        pet2anat_method='mri_coreg',
     )
 
 The PET reference volume is aligned to the skull-stripped anatomical image
-using FreeSurfer's ``mri_coreg`` with the number of degrees of freedom set via
-the :option:`--pet2anat-dof` flag. The resulting affine is converted to ITK
-format for downstream application, along with its inverse.
-
-If co-registration proves challenging, the :option:`--pet2anat-robust` flag
-switches the workflow to FreeSurfer's ``mri_robust_register`` with an NMI cost function and restricted
-to rigid-body (6 dof) transforms. This method is more robust to large initial misalignments.
+using the method selected via :option:`--pet2anat-method`. The anatomical image
+is first cropped with FSL's ``robustfov`` and masked to focus the alignment on
+brain tissue (ANTs receives the unmasked cropped image together with its mask).
+By default, the workflow runs FreeSurfer's ``mri_coreg`` with the number of
+degrees of freedom set via the :option:`--pet2anat-dof` flag. Alternative modes
+include FreeSurfer's ``mri_robust_register`` (``--pet2anat-method robust``) and
+ANTs rigid registration (``--pet2anat-method ants``). Both alternatives are
+limited to rigid-body alignment (6 DoF). The resulting affine is converted to
+ITK format for downstream application, along with its inverse.
 
 Resampling PET runs onto standard spaces
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -464,11 +454,10 @@ Resampling PET runs onto standard spaces
     )
 
 This sub-workflow concatenates the transforms calculated upstream (see
-`Head-motion estimation`_, `Susceptibility Distortion Correction (SDC)`_ --if
-fieldmaps are available--, `EPI to T1w registration`_, and an anatomical-to-standard
+`Head-motion estimation`_, `PET to T1w registration`_, and an anatomical-to-standard
 transform from `Preprocessing of structural MRI`_) to map the
-:abbr:`EPI (echo-planar imaging)`
-image to the standard spaces given by the ``--output-spaces`` argument
+:abbr:`PET (positron emission tomography)` series from native space to
+the standard spaces given by the ``--output-spaces`` argument
 (see :ref:`output-spaces`).
 It also maps the T1w-based mask to each of those standard spaces.
 
@@ -478,7 +467,7 @@ step, so as little information is lost as possible.
 The output space grid can be specified using modifiers to the ``--output-spaces``
 argument.
 
-EPI sampled to FreeSurfer surfaces
+PET sampled to FreeSurfer surfaces
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 :py:func:`~petprep.workflows.pet.resampling.init_pet_surf_wf`
 
@@ -495,7 +484,7 @@ EPI sampled to FreeSurfer surfaces
         output_dir='.',
     )
 
-If FreeSurfer processing is enabled, the motion-corrected functional series
+If FreeSurfer processing is enabled, the motion-corrected PET series
 (after single shot resampling to T1w space) is sampled to the
 surface by averaging across the cortical ribbon.
 Specifically, at each vertex, the segment normal to the white-matter surface, extending to the pial
diff --git a/petprep/cli/parser.py b/petprep/cli/parser.py
index c07844aa..0610f654 100644
--- a/petprep/cli/parser.py
+++ b/petprep/cli/parser.py
@@ -358,10 +358,14 @@ def _bids_filter(value, parser):
         '6 degrees (rotation and translation) are used by default.',
     )
     g_conf.add_argument(
-        '--pet2anat-robust',
-        action='store_true',
-        help='Use FreeSurfer mri_robust_register with an NMI cost function for'
-        'PET-to-T1w co-registration. This option is limited to 6 dof.',
+        '--pet2anat-method',
+        action='store',
+        default='mri_coreg',
+        choices=['mri_coreg', 'robust', 'ants'],
+        help='Method for PET-to-anatomical registration. '
+        '"mri_coreg" (default) uses FreeSurfer mri_coreg. '
+        '"robust" uses FreeSurfer mri_robust_register (6 DoF only). '
+        '"ants" uses ANTs rigid registration (6 DoF only).',
     )
     g_conf.add_argument(
         '--force-bbr',
@@ -575,6 +579,18 @@ def _bids_filter(value, parser):
         action='store_true',
         help='Disable head-motion correction and use the uncorrected data.',
     )
+    g_hmc.add_argument(
+        '--petref',
+        default='template',
+        choices=['template', 'twa', 'sum', 'first5min'],
+        help=(
+            "Strategy for generating the PET reference. 'template' uses the "
+            "motion correction template, while 'twa' computes a time-weighted "
+            "average, 'sum' produces a summed image of the motion-corrected "
+            "series, and 'first5min' averages the early (0-5 minute) portion "
+            'of the acquisition.'
+        ),
+    )
 
     g_seg = parser.add_argument_group('Segmentation options')
     g_seg.add_argument(
@@ -758,8 +774,9 @@ def parse_args(args=None, namespace=None):
     parser = _build_parser()
     opts = parser.parse_args(args, namespace)
 
-    if getattr(opts, 'pet2anat_robust', False) and opts.pet2anat_dof != 6:
-        parser.error('--pet2anat-robust requires --pet2anat-dof=6.')
+    # Validate DoF constraints for registration methods
+    if opts.pet2anat_method in ('robust', 'ants') and opts.pet2anat_dof != 6:
+        parser.error(f'--pet2anat-method {opts.pet2anat_method} requires --pet2anat-dof=6.')
 
     if opts.config_file:
         skip = {} if opts.reports_only else {'execution': ('run_uuid',)}
diff --git a/petprep/config.py b/petprep/config.py
index 8624b384..e4b85e7b 100644
--- a/petprep/config.py
+++ b/petprep/config.py
@@ -556,8 +556,10 @@ class workflow(_Config):
     """Degrees of freedom of the PET-to-anatomical registration steps."""
     pet2anat_init = 'auto'
     """Initial transform for PET-to-anatomical registration."""
-    pet2anat_robust = False
-    """Use ``mri_robust_register`` for PET-to-anatomical alignment."""
+    pet2anat_method: str = 'mri_coreg'
+    """PET-to-anatomical registration method (mri_coreg, robust, or ants)."""
+    petref: str = 'template'
+    """Strategy for building the PET reference (``'template'``, ``'twa'``, ``'sum'`` or ``'first5min'``)."""
     cifti_output = None
     """Generate HCP Grayordinates, accepts either ``'91k'`` (default) or ``'170k'``."""
     hires = None
diff --git a/petprep/interfaces/reports.py b/petprep/interfaces/reports.py
index 4f465d4b..41122315 100644
--- a/petprep/interfaces/reports.py
+++ b/petprep/interfaces/reports.py
@@ -98,6 +98,7 @@ def get_world_pedir(orientation: str, pe_dir: str) -> str:
 \t\t<ul class="elem-desc">
 \t\t\t<li>Original orientation: {ornt}</li>
 \t\t\t<li>Registration: {registration}</li>
+\t\t\t<li>Reference image: {reference}</li>
 \t\t\t<li>Time zero: {time_zero}</li>
 \t\t\t<li>Radiotracer: {radiotracer}</li>
 \t\t\t<li>Injected dose: {dose} {dose_units}</li>
@@ -228,6 +229,7 @@ class FunctionalSummaryInputSpec(TraitedSpec):
     registration = traits.Enum(
         'mri_coreg',
         'mri_robust_register',
+        'ants_registration',
         'Precomputed',
         mandatory=True,
         desc='PET/anatomical registration method',
@@ -237,6 +239,18 @@ class FunctionalSummaryInputSpec(TraitedSpec):
     )
     orientation = traits.Str(mandatory=True, desc='Orientation of the voxel axes')
     metadata = traits.Dict(desc='PET metadata dictionary')
+    petref_strategy = traits.Enum(
+        'template',
+        'twa',
+        'sum',
+        'first5min',
+        mandatory=True,
+        desc='PET reference generation strategy',
+    )
+    requested_petref_strategy = traits.Enum(
+        'template', 'twa', 'sum', 'first5min', desc='User-requested PET reference strategy'
+    )
+    hmc_disabled = traits.Bool(False, desc='Head motion correction disabled')
 
 
 class FunctionalSummary(SummaryInterface):
@@ -249,8 +263,25 @@ def _generate_segment(self):
             reg = 'Precomputed affine transformation'
         elif self.inputs.registration == 'mri_coreg':
             reg = f'FreeSurfer <code>mri_coreg</code> - {dof} dof'
+        elif self.inputs.registration == 'ants_registration':
+            reg = f'ANTs <code>ants_registration</code> ({dof} DoF)'
+        elif self.inputs.registration == 'mri_robust_register':
+            reg = 'FreeSurfer <code>mri_robust_register</code> (NMI cost)'
         else:
-            reg = 'FreeSurfer <code>mri_robust_register</code> (ROBENT cost)'
+            reg = f'Unknown registration method: {self.inputs.registration}'
+
+        reference_map = {
+            'template': 'Motion correction template',
+            'twa': 'Time-weighted average of motion-corrected series',
+            'sum': 'Summed motion-corrected series',
+            'first5min': 'Early (0-5 minute) average of motion-corrected series',
+        }
+        petref_strategy = reference_map.get(self.inputs.petref_strategy, 'Unknown')
+        requested = getattr(self.inputs, 'requested_petref_strategy', None)
+        if requested and requested != self.inputs.petref_strategy:
+            petref_strategy += f" (requested '{requested}')"
+        if self.inputs.hmc_disabled:
+            petref_strategy += ' (head motion correction disabled)'
 
         meta = self.inputs.metadata or {}
         time_zero = meta.get('TimeZero', None)
@@ -288,6 +319,7 @@ def _generate_segment(self):
 
         return FUNCTIONAL_TEMPLATE.format(
             registration=reg,
+            reference=petref_strategy,
             ornt=self.inputs.orientation,
             # Use the metadata dictionary to fill in the details
             time_zero=time_zero,
diff --git a/petprep/interfaces/tests/test_reports.py b/petprep/interfaces/tests/test_reports.py
index 97f7522b..90b16389 100644
--- a/petprep/interfaces/tests/test_reports.py
+++ b/petprep/interfaces/tests/test_reports.py
@@ -76,7 +76,7 @@ def test_subject_summary_handles_missing_task(tmp_path):
 
 @pytest.mark.parametrize(
     'registration',
-    ['mri_coreg', 'mri_robust_register'],
+    ['mri_coreg', 'mri_robust_register', 'ants_registration'],
 )
 def test_functional_summary_with_metadata(registration):
     from ..reports import FunctionalSummary
@@ -85,6 +85,7 @@ def test_functional_summary_with_metadata(registration):
         registration=registration,
         registration_dof=6,
         orientation='RAS',
+        petref_strategy='template',
         metadata={
             'TracerName': 'DASB',
             'TracerRadionuclide': '[11C]',
@@ -97,6 +98,7 @@ def test_functional_summary_with_metadata(registration):
 
     segment = summary._generate_segment()
     assert registration in segment
+    assert 'Reference image: Motion correction template' in segment
     assert 'Radiotracer: [11C]DASB' in segment
     assert 'Injected dose: 100 MBq' in segment
     assert 'Number of frames: 2' in segment
diff --git a/petprep/workflows/pet/fit.py b/petprep/workflows/pet/fit.py
index 06118a41..bcb03e66 100644
--- a/petprep/workflows/pet/fit.py
+++ b/petprep/workflows/pet/fit.py
@@ -24,7 +24,6 @@
 from pathlib import Path
 
 import nibabel as nb
-import numpy as np
 from nipype.interfaces import utility as niu
 from nipype.pipeline import engine as pe
 from nitransforms.linear import Affine, LinearTransformsMapping
@@ -57,12 +56,18 @@
 
 def _extract_twa_image(
     pet_file: str,
-    output_dir: Path,
-    frame_start_times: Sequence[float] | None,
-    frame_durations: Sequence[float] | None,
+    output_dir: 'Path',
+    frame_start_times: 'Sequence[float] | None',
+    frame_durations: 'Sequence[float] | None',
 ) -> str:
     """Return a time-weighted average (twa) reference image from a 4D PET series."""
 
+    from pathlib import Path
+
+    import nibabel as nb
+    import numpy as np
+
+    output_dir = Path(output_dir)
     output_dir.mkdir(parents=True, exist_ok=True)
     img = nb.load(pet_file)
     if img.ndim < 4 or img.shape[-1] == 1:
@@ -109,6 +114,100 @@ def _extract_twa_image(
     return str(out_file)
 
 
+def _extract_first5min_image(
+    pet_file: str,
+    output_dir: 'Path',
+    frame_start_times: 'Sequence[float] | None',
+    frame_durations: 'Sequence[float] | None',
+    window_sec: float = 300.0,
+) -> str:
+    """Average the early (0-``window_sec``) portion of a PET series."""
+
+    from pathlib import Path
+
+    import nibabel as nb
+    import numpy as np
+
+    output_dir = Path(output_dir)
+    output_dir.mkdir(parents=True, exist_ok=True)
+    img = nb.load(pet_file)
+    if img.ndim < 4 or img.shape[-1] == 1:
+        return pet_file
+
+    if frame_start_times is None or frame_durations is None:
+        raise ValueError('Frame timing metadata are required to compute an early reference image.')
+
+    frame_start_times = np.asarray(frame_start_times, dtype=float)
+    frame_durations = np.asarray(frame_durations, dtype=float)
+
+    if frame_start_times.ndim != 1 or frame_durations.ndim != 1:
+        raise ValueError('Frame timing metadata must be one-dimensional sequences.')
+
+    if len(frame_start_times) != len(frame_durations):
+        raise ValueError('FrameTimesStart and FrameDuration must have the same length.')
+
+    if len(frame_durations) != img.shape[-1]:
+        raise ValueError(
+            'Frame timing metadata must match the number of frames in the PET series.'
+        )
+
+    if np.any(frame_durations <= 0):
+        raise ValueError('FrameDuration values must all be positive.')
+
+    if np.any(np.diff(frame_start_times) < 0):
+        raise ValueError('FrameTimesStart values must be non-decreasing.')
+
+    frame_ends = frame_start_times + frame_durations
+    included_durations = np.clip(frame_ends, 0.0, window_sec) - np.clip(
+        frame_start_times, 0.0, window_sec
+    )
+
+    if not np.any(included_durations > 0):
+        raise ValueError('No frames overlap with the first 5 minutes of the acquisition.')
+
+    hdr = img.header.copy()
+    data = np.asanyarray(img.dataobj)
+    weighted_average = np.average(data, axis=-1, weights=included_durations).astype(np.float32)
+    hdr.set_data_shape(weighted_average.shape)
+
+    pet_path = Path(pet_file)
+    while pet_path.suffix:
+        pet_path = pet_path.with_suffix('')
+
+    out_file = output_dir / f'{pet_path.name}_first5minref.nii.gz'
+    img.__class__(weighted_average, img.affine, hdr).to_filename(out_file)
+    return str(out_file)
+
+
+def _extract_sum_image(pet_file: str, output_dir: 'Path') -> str:
+    """Return a summed reference image from a 4D PET series."""
+
+    from pathlib import Path
+
+    import nibabel as nb
+    import numpy as np
+
+    output_dir = Path(output_dir)
+    output_dir.mkdir(parents=True, exist_ok=True)
+    img = nb.load(pet_file)
+    if img.ndim < 4 or img.shape[-1] == 1:
+        return pet_file
+
+    hdr = img.header.copy()
+    data = np.asanyarray(img.dataobj)
+    summed = data.sum(axis=-1, dtype=np.float64).astype(np.float32)
+    hdr.set_data_shape(summed.shape)
+
+    pet_path = Path(pet_file)
+    pet_stem = pet_path
+    while pet_stem.suffix:
+        pet_stem = pet_stem.with_suffix('')
+
+    out_file = output_dir / f'{pet_stem.name}_sumref.nii.gz'
+    img.__class__(summed, img.affine, hdr).to_filename(out_file)
+    return str(out_file)
+
+
 def _write_identity_xforms(num_frames: int, filename: Path) -> Path:
     """Write ``num_frames`` identity transforms to ``filename``."""
 
@@ -286,20 +385,93 @@ def init_pet_fit_wf(
             'Please check your BIDS JSON sidecar.'
         )
 
+    requested_petref_strategy = getattr(config.workflow, 'petref', 'template')
+    hmc_disabled = bool(config.workflow.hmc_off)
+    petref_strategy = requested_petref_strategy
+    if hmc_disabled and petref_strategy == 'template':
+        config.loggers.workflow.warning(
+            'Head motion correction disabled (--hmc-off); using a time-weighted average '
+            'reference instead of the motion correction template.'
+        )
+        petref_strategy = 'twa'
+
+    use_corrected_reference = petref_strategy in {'twa', 'sum', 'first5min'}
+    reference_function = _extract_twa_image
+    reference_kwargs: dict[str, object] = {
+        'output_dir': config.execution.work_dir,
+        'frame_start_times': frame_start_times,
+        'frame_durations': frame_durations,
+    }
+    reference_node_name = 'twa_reference'
+    reference_input_names = ['pet_file', 'output_dir', 'frame_start_times', 'frame_durations']
+    report_reference_input_names = reference_input_names
+
+    if petref_strategy == 'sum':
+        reference_function = _extract_sum_image
+        reference_kwargs = {'output_dir': config.execution.work_dir}
+        reference_node_name = 'sum_reference'
+        reference_input_names = ['pet_file', 'output_dir']
+
+    if petref_strategy == 'first5min':
+        reference_function = _extract_first5min_image
+        reference_node_name = 'first5min_reference'
+
+    corrected_pet_for_report = None
+    corrected_reference = None
+
+    requires_report_reference = pet_tlen > 1
+    report_pet_for_coreg = None
+    report_pet_reference = None
+
+    if requires_report_reference:
+        report_pet_for_coreg = pe.Node(
+            ResampleSeries(),
+            name='report_pet_for_coreg',
+            n_procs=omp_nthreads,
+            mem_gb=mem_gb['resampled'],
+        )
+        report_pet_reference = pe.Node(
+            niu.Function(
+                function=_extract_twa_image,
+                input_names=report_reference_input_names,
+                output_names=['out_file'],
+            ),
+            name='report_petref',
+        )
+        report_pet_reference.inputs.output_dir = config.execution.work_dir
+        report_pet_reference.inputs.frame_start_times = frame_start_times
+        report_pet_reference.inputs.frame_durations = frame_durations
+
+    if use_corrected_reference:
+        corrected_pet_for_report = pe.Node(
+            ResampleSeries(),
+            name='corrected_pet_for_report',
+            n_procs=omp_nthreads,
+            mem_gb=mem_gb['resampled'],
+        )
+        corrected_reference = pe.Node(
+            niu.Function(
+                function=reference_function,
+                input_names=reference_input_names,
+                output_names=['out_file'],
+            ),
+            name=reference_node_name,
+        )
+        corrected_reference.inputs.output_dir = config.execution.work_dir
+        if petref_strategy in {'twa', 'first5min'}:
+            corrected_reference.inputs.frame_start_times = frame_start_times
+            corrected_reference.inputs.frame_durations = frame_durations
+
     registration_method = 'Precomputed'
     if not petref2anat_xform:
-        registration_method = (
-            'mri_robust_register' if config.workflow.pet2anat_robust else 'mri_coreg'
-        )
-    hmc_disabled = bool(config.workflow.hmc_off)
+        registration_method = {
+            'mri_coreg': 'mri_coreg',
+            'robust': 'mri_robust_register',
+            'ants': 'ants_registration',
+        }[config.workflow.pet2anat_method]
     if hmc_disabled:
         config.execution.work_dir.mkdir(parents=True, exist_ok=True)
-        petref = petref or _extract_twa_image(
-            pet_file,
-            config.execution.work_dir,
-            frame_start_times,
-            frame_durations,
-        )
+        petref = petref or reference_function(pet_file, **reference_kwargs)
         idmat_fname = config.execution.work_dir / 'idmat.tfm'
         n_frames = len(frame_durations)
         hmc_xforms = _write_identity_xforms(n_frames, idmat_fname)
@@ -323,6 +495,9 @@ def init_pet_fit_wf(
             registration_dof=config.workflow.pet2anat_dof,
             orientation=orientation,
             metadata=metadata,
+            petref_strategy=petref_strategy,
+            requested_petref_strategy=requested_petref_strategy,
+            hmc_disabled=hmc_disabled,
         ),
         name='summary',
         mem_gb=config.DEFAULT_MEMORY_MIN_GB,
@@ -402,9 +577,27 @@ def init_pet_fit_wf(
             ]),
             (pet_hmc_wf, ds_hmc_wf, [('outputnode.xforms', 'inputnode.xforms')]),
             (ds_hmc_wf, hmc_buffer, [('outputnode.xforms', 'hmc_xforms')]),
-            (pet_hmc_wf, petref_buffer, [('outputnode.petref', 'petref')]),
             (pet_hmc_wf, ds_petref_wf, [('outputnode.petref', 'inputnode.petref')]),
         ])  # fmt:skip
+
+        if use_corrected_reference:
+            workflow.connect([
+                (pet_hmc_wf, corrected_pet_for_report, [('outputnode.petref', 'ref_file')]),
+                (val_pet, corrected_pet_for_report, [('out_file', 'in_file')]),
+                (hmc_buffer, corrected_pet_for_report, [('hmc_xforms', 'transforms')]),
+                (corrected_pet_for_report, corrected_reference, [('out_file', 'pet_file')]),
+                (corrected_reference, petref_buffer, [('out_file', 'petref')]),
+            ])  # fmt:skip
+        else:
+            workflow.connect([(pet_hmc_wf, petref_buffer, [('outputnode.petref', 'petref')])])
+
+        if report_pet_reference:
+            workflow.connect([
+                (pet_hmc_wf, report_pet_for_coreg, [('outputnode.petref', 'ref_file')]),
+                (val_pet, report_pet_for_coreg, [('out_file', 'in_file')]),
+                (hmc_buffer, report_pet_for_coreg, [('hmc_xforms', 'transforms')]),
+                (report_pet_for_coreg, report_pet_reference, [('out_file', 'pet_file')]),
+            ])  # fmt:skip
     else:
         config.loggers.workflow.info(
             'PET Stage 1: Found head motion correction transforms and petref - skipping Stage 1'
@@ -418,9 +611,45 @@ def init_pet_fit_wf(
 
         ])  # fmt:skip
         val_pet.inputs.in_file = pet_file
-        petref_buffer.inputs.petref = petref
+        if use_corrected_reference:
+            corrected_pet_for_report.inputs.ref_file = petref
+            workflow.connect(
+                [
+                    (val_pet, corrected_pet_for_report, [('out_file', 'in_file')]),
+                    (hmc_buffer, corrected_pet_for_report, [('hmc_xforms', 'transforms')]),
+                    (corrected_pet_for_report, corrected_reference, [('out_file', 'pet_file')]),
+                    (corrected_reference, petref_buffer, [('out_file', 'petref')]),
+                ]
+            )  # fmt:skip
+        else:
+            petref_buffer.inputs.petref = petref
+
+        if report_pet_reference:
+            workflow.connect([
+                (petref_buffer, report_pet_for_coreg, [('petref', 'ref_file')]),
+                (val_pet, report_pet_for_coreg, [('out_file', 'in_file')]),
+                (hmc_buffer, report_pet_for_coreg, [('hmc_xforms', 'transforms')]),
+                (report_pet_for_coreg, report_pet_reference, [('out_file', 'pet_file')]),
+            ])  # fmt:skip
+
+    if report_pet_reference:
+        workflow.connect([
+            (report_pet_reference, func_fit_reports_wf, [('out_file', 'inputnode.report_pet')])
+        ])  # fmt:skip
+    else:
+        func_fit_reports_wf.inputs.inputnode.report_pet = pet_file
 
     # Stage 2: Coregistration
+    rerun_coreg = petref2anat_xform and (
+        config.workflow.petref != 'template' or config.workflow.pet2anat_method != 'mri_coreg'
+    )
+    if rerun_coreg:
+        config.loggers.workflow.info(
+            'PET Stage 2: Re-running co-registration because --petref or --pet2anat-method '
+            'were set to non-default values.'
+        )
+        petref2anat_xform = None
+
     if not petref2anat_xform:
         config.loggers.workflow.info('PET Stage 2: Adding co-registration workflow of PET to T1w')
         # calculate PET registration to T1w
@@ -428,7 +657,7 @@ def init_pet_fit_wf(
             pet2anat_dof=config.workflow.pet2anat_dof,
             omp_nthreads=omp_nthreads,
             mem_gb=mem_gb['resampled'],
-            use_robust_register=config.workflow.pet2anat_robust,
+            pet2anat_method=config.workflow.pet2anat_method,
             sloppy=config.execution.sloppy,
         )
 
diff --git a/petprep/workflows/pet/outputs.py b/petprep/workflows/pet/outputs.py
index d4c726c5..389660e3 100644
--- a/petprep/workflows/pet/outputs.py
+++ b/petprep/workflows/pet/outputs.py
@@ -148,6 +148,7 @@ def init_func_fit_reports_wf(
     inputfields = [
         'source_file',
         'petref',
+        'report_pet',
         'pet_mask',
         'petref2anat_xfm',
         't1w_preproc',
@@ -261,16 +262,16 @@ def init_func_fit_reports_wf(
         ]),
         (inputnode, t1w_petref, [
             ('t1w_preproc', 'input_image'),
-            ('petref', 'reference_image'),
+            ('report_pet', 'reference_image'),
             ('petref2anat_xfm', 'transforms'),
         ]),
         (inputnode, t1w_wm, [('t1w_dseg', 'in_seg')]),
         (inputnode, petref_wm, [
-            ('petref', 'reference_image'),
+            ('report_pet', 'reference_image'),
             ('petref2anat_xfm', 'transforms'),
         ]),
         (t1w_wm, petref_wm, [('out', 'input_image')]),
-        (inputnode, crop_petref, [('petref', 'in_file'), ('pet_mask', 'mask_file')]),
+        (inputnode, crop_petref, [('report_pet', 'in_file'), ('pet_mask', 'mask_file')]),
         (t1w_petref, crop_t1w_petref, [('output_image', 'in_file')]),
         (inputnode, crop_t1w_petref, [('pet_mask', 'mask_file')]),
         (petref_wm, crop_petref_wm, [('output_image', 'in_file')]),
@@ -284,7 +285,7 @@ def init_func_fit_reports_wf(
             ]),
             (inputnode, petref_refmask, [
                 ('refmask', 'input_image'),
-                ('petref', 'reference_image'),
+                ('report_pet', 'reference_image'),
                 ('petref2anat_xfm', 'transforms'),
             ]),
             (petref_refmask, crop_petref_refmask, [('output_image', 'in_file')]),
diff --git a/petprep/workflows/pet/registration.py b/petprep/workflows/pet/registration.py
index 370da756..2d4a2d63 100644
--- a/petprep/workflows/pet/registration.py
+++ b/petprep/workflows/pet/registration.py
@@ -36,12 +36,17 @@
 AffineDOF = ty.Literal[6, 9, 12]
 
 
+def _get_first(in_list):
+    """Extract first element from a list (for ANTs transform output)."""
+    return in_list[0]
+
+
 def init_pet_reg_wf(
     *,
     pet2anat_dof: AffineDOF,
     mem_gb: float,
     omp_nthreads: int,
-    use_robust_register: bool = False,
+    pet2anat_method: str = 'mri_coreg',
     name: str = 'pet_reg_wf',
     sloppy: bool = False,
 ):
@@ -70,10 +75,11 @@ def init_pet_reg_wf(
         Size of PET file in GB
     omp_nthreads : :obj:`int`
         Maximum number of threads an individual process may use
-    use_robust_register : :obj:`bool`
-        Run FreeSurfer ``mri_robust_register`` with an NMI cost function for
-        PET-to-anatomical alignment. Only rigid-body (6 dof) alignment is
-        supported in this mode.
+    pet2anat_method : :obj:`str`
+        Method for PET-to-anatomical registration. Options are 'mri_coreg'
+        (default FreeSurfer co-registration), 'robust' (uses FreeSurfer
+        mri_robust_register with NMI, 6 DoF only), or 'ants' (uses ANTs rigid
+        registration, 6 DoF only).
     name : :obj:`str`
         Name of workflow (default: ``pet_reg_wf``)
 
@@ -95,7 +101,9 @@ def init_pet_reg_wf(
         Affine transform from anatomical space to PET space (ITK format)
 
     """
-    from nipype.interfaces.freesurfer import MRICoreg, RobustRegister
+    from nipype.interfaces.ants import Registration
+    from nipype.interfaces.freesurfer import MRIConvert, MRICoreg, RobustRegister
+    from nipype.interfaces.fsl import RobustFOV
     from niworkflows.engine.workflows import LiterateWorkflow as Workflow
     from niworkflows.interfaces.nibabel import ApplyMask
     from niworkflows.interfaces.nitransforms import ConcatenateXFMs
@@ -112,7 +120,45 @@ def init_pet_reg_wf(
     )
 
     mask_brain = pe.Node(ApplyMask(), name='mask_brain')
-    if use_robust_register:
+    crop_anat_mask = pe.Node(MRIConvert(out_type='niigz'), name='crop_anat_mask')
+    robust_fov = pe.Node(RobustFOV(output_type='NIFTI_GZ'), name='robust_fov')
+
+    if pet2anat_method == 'ants':
+        coreg = pe.Node(
+            Registration(
+                dimension=3,
+                float=True,
+                output_transform_prefix='pet2anat_',
+                output_warped_image='pet2anat_Warped.nii.gz',
+                transforms=['Rigid'],
+                transform_parameters=[(0.1,)],
+                metric=['MI'],
+                metric_weight=[1],
+                radius_or_number_of_bins=[32],
+                sampling_strategy=['Regular'],
+                sampling_percentage=[0.25],
+                number_of_iterations=[[1000, 500, 250]],
+                convergence_threshold=[1e-6],
+                convergence_window_size=[10],
+                shrink_factors=[[4, 2, 1]],
+                smoothing_sigmas=[[2, 1, 0]],
+                sigma_units=['vox'],
+                use_histogram_matching=False,
+                initial_moving_transform_com=1,
+                winsorize_lower_quantile=0.001,
+                winsorize_upper_quantile=0.999,
+            ),
+            name='ants_registration',
+            n_procs=omp_nthreads,
+            mem_gb=mem_gb * 2,
+        )
+        coreg_target = 'fixed_image'
+        coreg_mask = 'fixed_image_masks'
+        coreg_moving = 'moving_image'
+        coreg_output = 'forward_transforms'
+        coreg_output_is_list = True
+
+    elif pet2anat_method == 'robust':
         coreg = pe.Node(
             RobustRegister(
                 auto_sens=False,
@@ -128,8 +174,11 @@ def init_pet_reg_wf(
             mem_gb=5,
         )
         coreg_target = 'target_file'
+        coreg_moving = 'source_file'
         coreg_output = 'out_reg_file'
-    else:
+        coreg_output_is_list = False
+
+    else:  # mri_coreg (default)
         coreg = pe.Node(
             MRICoreg(dof=pet2anat_dof, sep=[4], ftol=0.0001, linmintol=0.01),
             name='mri_coreg',
@@ -137,31 +186,63 @@ def init_pet_reg_wf(
             mem_gb=5,
         )
         coreg_target = 'reference_file'
+        coreg_moving = 'source_file'
         coreg_output = 'out_lta_file'
+        coreg_output_is_list = False
+
     convert_xfm = pe.Node(ConcatenateXFMs(inverse=True), name='convert_xfm')
 
-    connections = [
-        (
-            inputnode,
-            mask_brain,
-            [
-                ('anat_preproc', 'in_file'),
-                ('anat_mask', 'in_mask'),
-            ],
-        ),
-        (inputnode, coreg, [('ref_pet_brain', 'source_file')]),
-        (mask_brain, coreg, [('out_file', coreg_target)]),
-        (coreg, convert_xfm, [(coreg_output, 'in_xfms')]),
-        (
-            convert_xfm,
-            outputnode,
-            [
-                ('out_xfm', 'itk_pet_to_t1'),
-                ('out_inv', 'itk_t1_to_pet'),
-            ],
-        ),
-    ]
-
-    workflow.connect(connections)  # fmt:skip
+    # Build connections dynamically based on output type
+    if coreg_output_is_list:
+        # ANTs outputs a list of transforms; take the first (and only) one
+        # ANTs gets unmasked T1W + separate mask (not pre-masked image)
+        connections = [
+            (robust_fov, mask_brain, [('out_roi', 'in_file')]),
+            (crop_anat_mask, mask_brain, [('out_file', 'in_mask')]),
+            (inputnode, coreg, [('ref_pet_brain', coreg_moving)]),
+            (
+                robust_fov,
+                coreg,
+                [
+                    ('out_roi', coreg_target),
+                ],
+            ),
+            (crop_anat_mask, coreg, [('out_file', coreg_mask)]),
+            (coreg, convert_xfm, [((coreg_output, _get_first), 'in_xfms')]),
+            (
+                convert_xfm,
+                outputnode,
+                [
+                    ('out_xfm', 'itk_pet_to_t1'),
+                    ('out_inv', 'itk_t1_to_pet'),
+                ],
+            ),
+        ]
+    else:
+        # mri_coreg and mri_robust_register output single transform file
+        connections = [
+            (robust_fov, mask_brain, [('out_roi', 'in_file')]),
+            (crop_anat_mask, mask_brain, [('out_file', 'in_mask')]),
+            (inputnode, coreg, [('ref_pet_brain', coreg_moving)]),
+            (mask_brain, coreg, [('out_file', coreg_target)]),
+            (coreg, convert_xfm, [(coreg_output, 'in_xfms')]),
+            (
+                convert_xfm,
+                outputnode,
+                [
+                    ('out_xfm', 'itk_pet_to_t1'),
+                    ('out_inv', 'itk_t1_to_pet'),
+                ],
+            ),
+        ]
+
+    workflow.connect(
+        [
+            (inputnode, robust_fov, [('anat_preproc', 'in_file')]),
+            (inputnode, crop_anat_mask, [('anat_mask', 'in_file')]),
+            (robust_fov, crop_anat_mask, [('out_roi', 'reslice_like')]),
+        ]
+        + connections
+    )  # fmt:skip
 
     return workflow
diff --git a/petprep/workflows/pet/tests/test_fit.py b/petprep/workflows/pet/tests/test_fit.py
index ee95cac9..05c60ad9 100644
--- a/petprep/workflows/pet/tests/test_fit.py
+++ b/petprep/workflows/pet/tests/test_fit.py
@@ -5,6 +5,7 @@
 import numpy as np
 import pytest
 import yaml
+from nipype.interfaces.base import Undefined
 from nipype.pipeline.engine.utils import generate_expanded_graph
 from niworkflows.utils.testing import generate_bids_skeleton
 
@@ -12,7 +13,14 @@
 from ....utils import bids
 from ...tests import mock_config
 from ...tests.test_base import BASE_LAYOUT
-from ..fit import _extract_twa_image, init_pet_fit_wf, init_pet_native_wf
+from ..fit import (
+    _extract_first5min_image,
+    _extract_sum_image,
+    _extract_twa_image,
+    _write_identity_xforms,
+    init_pet_fit_wf,
+    init_pet_native_wf,
+)
 from ..outputs import init_refmask_report_wf
 
 
@@ -176,20 +184,51 @@ def test_pet_fit_mask_connections(bids_root: Path, tmp_path: Path):
     assert ('out', 'inputnode.petmask') in ds_edge['connect']
 
 
-def test_petref_report_connections(bids_root: Path, tmp_path: Path):
-    """Ensure the PET reference is passed to the reports workflow."""
+def test_reports_use_motion_corrected_average(bids_root: Path, tmp_path: Path):
+    """Co-registration report should show the motion corrected time-weighted average."""
+
     pet_series = [str(bids_root / 'sub-01' / 'pet' / 'sub-01_task-rest_run-1_pet.nii.gz')]
-    img = nb.Nifti1Image(np.zeros((2, 2, 2, 1)), np.eye(4))
+    data = np.stack((np.ones((2, 2, 2)), np.full((2, 2, 2), 2.0)), axis=-1)
+    img = nb.Nifti1Image(data, np.eye(4))
+    for path in pet_series:
+        img.to_filename(path)
+
+    sidecar = Path(pet_series[0]).with_suffix('').with_suffix('.json')
+    sidecar.write_text('{"FrameTimesStart": [0, 1], "FrameDuration": [1, 1]}')
+
+    with mock_config(bids_dir=bids_root):
+        wf = init_pet_fit_wf(pet_series=pet_series, precomputed={}, omp_nthreads=1)
+
+    assert 'report_petref' in wf.list_node_names()
+    reports_node = wf.get_node('func_fit_reports_wf')
+    report_petref = wf.get_node('report_petref')
+    edge = wf._graph.get_edge_data(report_petref, reports_node)
+    assert ('out_file', 'inputnode.report_pet') in edge['connect']
+
 
+def test_petref_default_twa_when_hmc_disabled(bids_root: Path, tmp_path: Path):
+    """Disabling HMC should fall back to TWA references and note it in reports."""
+
+    pet_series = [str(bids_root / 'sub-01' / 'pet' / 'sub-01_task-rest_run-1_pet.nii.gz')]
+    data = np.stack((np.ones((2, 2, 2)), np.full((2, 2, 2), 2.0)), axis=-1)
+    img = nb.Nifti1Image(data, np.eye(4))
     for path in pet_series:
         img.to_filename(path)
 
+    sidecar = Path(pet_series[0]).with_suffix('').with_suffix('.json')
+    sidecar.write_text('{"FrameTimesStart": [0, 1], "FrameDuration": [1, 1]}')
+
     with mock_config(bids_dir=bids_root):
+        config.workflow.hmc_off = True
+        config.workflow.petref = 'template'
         wf = init_pet_fit_wf(pet_series=pet_series, precomputed={}, omp_nthreads=1)
 
-    petref_buffer = wf.get_node('petref_buffer')
-    edge = wf._graph.get_edge_data(petref_buffer, wf.get_node('func_fit_reports_wf'))
-    assert ('petref', 'inputnode.petref') in edge['connect']
+    assert 'twa_reference' in wf.list_node_names()
+
+    summary = wf.get_node('summary')
+    assert summary.inputs.petref_strategy == 'twa'
+    assert summary.inputs.requested_petref_strategy == 'template'
+    assert summary.inputs.hmc_disabled is True
 
 
 @pytest.mark.parametrize('pvc_method', [None, 'gtm'])
@@ -299,13 +338,35 @@ def test_pet_fit_robust_registration(bids_root: Path, tmp_path: Path):
         )
 
     with mock_config(bids_dir=bids_root):
-        config.workflow.pet2anat_robust = True
+        config.workflow.pet2anat_method = 'robust'
         config.workflow.pet2anat_dof = 6
         wf = init_pet_fit_wf(pet_series=pet_series, precomputed={}, omp_nthreads=1)
 
     node_names = wf.list_node_names()
     assert 'pet_reg_wf.mri_robust_register' in node_names
     assert 'pet_reg_wf.mri_coreg' not in node_names
+    assert 'pet_reg_wf.ants_registration' not in node_names
+
+
+def test_init_pet_fit_wf_ants_registration(bids_root: Path, tmp_path: Path):
+    """Test PET fit workflow with ANTs registration."""
+    pet_series = [str(bids_root / 'sub-01' / 'pet' / 'sub-01_task-rest_run-1_pet.nii.gz')]
+    img = nb.Nifti1Image(np.zeros((2, 2, 2, 1)), np.eye(4))
+    for path in pet_series:
+        img.to_filename(path)
+        Path(path).with_suffix('').with_suffix('.json').write_text(
+            '{"FrameTimesStart": [0], "FrameDuration": [1]}'
+        )
+
+    with mock_config(bids_dir=bids_root):
+        config.workflow.pet2anat_method = 'ants'
+        config.workflow.pet2anat_dof = 6
+        wf = init_pet_fit_wf(pet_series=pet_series, precomputed={}, omp_nthreads=1)
+
+    node_names = wf.list_node_names()
+    assert 'pet_reg_wf.ants_registration' in node_names
+    assert 'pet_reg_wf.mri_coreg' not in node_names
+    assert 'pet_reg_wf.mri_robust_register' not in node_names
 
 
 def test_pet_fit_requires_both_derivatives(bids_root: Path, tmp_path: Path):
@@ -371,6 +432,50 @@ def test_pet_fit_stage1_with_cached_baseline(bids_root: Path, tmp_path: Path):
     assert not any(name.startswith('pet_hmc_wf') for name in wf.list_node_names())
 
 
+def test_pet_fit_reruns_coreg_with_custom_options(bids_root: Path, tmp_path: Path):
+    """Custom co-registration options should ignore cached transforms."""
+
+    pet_series = [str(bids_root / 'sub-01' / 'pet' / 'sub-01_task-rest_run-1_pet.nii.gz')]
+    img = nb.Nifti1Image(np.zeros((2, 2, 2, 1)), np.eye(4))
+    for path in pet_series:
+        img.to_filename(path)
+
+    deriv_root = tmp_path / 'derivs'
+    petref = deriv_root / 'sub-01' / 'pet' / 'sub-01_task-rest_run-1_desc-hmc_petref.nii.gz'
+    hmc_xfm = (
+        deriv_root
+        / 'sub-01'
+        / 'pet'
+        / 'sub-01_task-rest_run-1_from-orig_to-petref_mode-image_xfm.txt'
+    )
+    petref2anat_xfm = (
+        deriv_root
+        / 'sub-01'
+        / 'pet'
+        / 'sub-01_task-rest_run-1_from-petref_to-anat_mode-image_xfm.txt'
+    )
+
+    petref.parent.mkdir(parents=True)
+    img.to_filename(petref)
+    np.savetxt(hmc_xfm, np.eye(4))
+    np.savetxt(petref2anat_xfm, np.eye(4))
+
+    sidecar = Path(pet_series[0]).with_suffix('').with_suffix('.json')
+    sidecar.write_text('{"FrameTimesStart": [0], "FrameDuration": [1]}')
+
+    entities = bids.extract_entities(pet_series)
+    precomputed = bids.collect_derivatives(derivatives_dir=deriv_root, entities=entities)
+
+    with mock_config(bids_dir=bids_root):
+        config.workflow.petref = 'twa'
+        config.workflow.pet2anat_method = 'robust'
+        wf = init_pet_fit_wf(pet_series=pet_series, precomputed=precomputed, omp_nthreads=1)
+
+    node_names = wf.list_node_names()
+    assert 'pet_reg_wf.mri_robust_register' in node_names
+    assert wf.get_node('outputnode').inputs.petref2anat_xfm is Undefined
+
+
 def test_pet_fit_hmc_off_disables_stage1(bids_root: Path, tmp_path: Path):
     """Disabling HMC should skip Stage 1 and use identity transforms."""
     pet_series = [str(bids_root / 'sub-01' / 'pet' / 'sub-01_task-rest_run-1_pet.nii.gz')]
@@ -436,6 +541,67 @@ def test_extract_twa_image_validation(
         )
 
 
+def test_extract_sum_image(tmp_path: Path):
+    """Summed references are written out with the expected contents."""
+
+    data = np.stack((np.ones((2, 2, 2)), np.full((2, 2, 2), 3.0)), axis=-1)
+    pet_img = nb.Nifti1Image(data.astype(np.float32), np.eye(4))
+    pet_file = tmp_path / 'pet.nii.gz'
+    pet_img.to_filename(pet_file)
+
+    out_file = _extract_sum_image(str(pet_file), tmp_path / 'out')
+
+    summed = nb.load(out_file).get_fdata()
+    assert np.allclose(summed, 4.0)
+    assert Path(out_file).name == 'pet_sumref.nii.gz'
+
+    # 3D inputs should round-trip without creating a new file
+    pet_3d = tmp_path / 'pet3d.nii.gz'
+    nb.Nifti1Image(np.zeros((2, 2, 2), dtype=np.float32), np.eye(4)).to_filename(pet_3d)
+    assert _extract_sum_image(str(pet_3d), tmp_path / 'out') == str(pet_3d)
+
+
+def test_extract_first5min_image(tmp_path: Path):
+    """Early reference averages only the first 5 minutes of data."""
+
+    data = np.stack((np.ones((2, 2, 2)), np.full((2, 2, 2), 3.0)), axis=-1)
+    pet_img = nb.Nifti1Image(data.astype(np.float32), np.eye(4))
+    pet_file = tmp_path / 'pet.nii.gz'
+    pet_img.to_filename(pet_file)
+
+    out_file = _extract_first5min_image(
+        str(pet_file),
+        tmp_path / 'out',
+        frame_start_times=[0, 400],
+        frame_durations=[400, 200],
+    )
+
+    averaged = nb.load(out_file).get_fdata()
+    expected = (1.0 * 300 + 3.0 * 0) / 300
+    assert np.allclose(averaged, expected)
+    assert Path(out_file).name == 'pet_first5minref.nii.gz'
+
+
+def test_report_petref_receives_frame_metadata(bids_root: Path, tmp_path: Path):
+    """Report reference node always receives timing metadata."""
+
+    pet_series = [str(bids_root / 'sub-01' / 'pet' / 'sub-01_task-rest_run-1_pet.nii.gz')]
+    img = nb.Nifti1Image(np.zeros((2, 2, 2, 2)), np.eye(4))
+    for path in pet_series:
+        img.to_filename(path)
+        Path(path).with_suffix('').with_suffix('.json').write_text(
+            '{"FrameTimesStart": [0, 1], "FrameDuration": [1, 1]}'
+        )
+
+    with mock_config(bids_dir=bids_root):
+        config.workflow.petref = 'sum'
+        wf = init_pet_fit_wf(pet_series=pet_series, precomputed={}, omp_nthreads=1)
+
+    report_petref = wf.get_node('report_petref')
+    assert report_petref.inputs.frame_start_times == [0, 1]
+    assert report_petref.inputs.frame_durations == [1, 1]
+
+
 def test_pet_fit_hmc_off_ignores_precomputed(bids_root: Path, tmp_path: Path):
     """Precomputed derivatives are ignored when ``--hmc-off`` is set."""
 
@@ -473,6 +639,20 @@ def test_pet_fit_hmc_off_ignores_precomputed(bids_root: Path, tmp_path: Path):
     assert Path(hmc_buffer.inputs.hmc_xforms).name == 'idmat.tfm'
 
 
+def test_write_identity_xforms_minimum(tmp_path: Path):
+    """At least one identity transform should always be written."""
+
+    xfm_file = _write_identity_xforms(0, tmp_path / 'idmat.tfm')
+
+    xforms = nt.linear.load(xfm_file)
+    matrices = np.asarray(xforms.matrix)
+    if matrices.ndim == 2:
+        matrices = matrices[np.newaxis, ...]
+
+    assert matrices.shape[0] == 1
+    assert np.allclose(matrices[0], np.eye(4))
+
+
 def test_init_refmask_report_wf(tmp_path: Path):
     """Ensure the refmask report workflow initializes without errors."""
     wf = init_refmask_report_wf(output_dir=str(tmp_path), ref_name='test')