From 60a9da70b952965bf686c5e0da2e19b3c4b819d1 Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Wed, 17 Jun 2026 16:03:24 +0200 Subject: [PATCH 1/3] ENH: Fix HMC-off identity transform collisions across PET series --- petprep/interfaces/resampling.py | 9 +++ petprep/interfaces/tests/test_resampling.py | 14 ++++ petprep/workflows/pet/fit.py | 18 ++++- petprep/workflows/pet/tests/test_fit.py | 85 ++++++++++++++++++++- 4 files changed, 120 insertions(+), 6 deletions(-) create mode 100644 petprep/interfaces/tests/test_resampling.py diff --git a/petprep/interfaces/resampling.py b/petprep/interfaces/resampling.py index 174d9f18..fa1b6c85 100644 --- a/petprep/interfaces/resampling.py +++ b/petprep/interfaces/resampling.py @@ -213,6 +213,15 @@ async def resample_series_async( The resampled array, with shape ``coordinates.shape[1:] + (N,)``, where N is the number of volumes in ``data``. """ + if hmc_xfms: + n_volumes = data.shape[-1] if data.ndim > 3 else 1 + if len(hmc_xfms) != n_volumes: + raise ValueError( + 'Head-motion transform count ' + f'({len(hmc_xfms)}) does not match the number of PET volumes ' + f'({n_volumes}).' + ) + if data.ndim == 3: return resample_vol( data, diff --git a/petprep/interfaces/tests/test_resampling.py b/petprep/interfaces/tests/test_resampling.py new file mode 100644 index 00000000..9cc0b985 --- /dev/null +++ b/petprep/interfaces/tests/test_resampling.py @@ -0,0 +1,14 @@ +import numpy as np +import pytest + +from ..resampling import resample_series + + +def test_resample_series_rejects_mismatched_hmc_transforms(): + """HMC transform mappings must contain one affine per PET volume.""" + + data = np.zeros((2, 2, 2, 2), dtype=np.float32) + coordinates = np.indices(data.shape[:3], dtype=np.float32) + + with pytest.raises(ValueError, match='Head-motion transform count .* PET volumes'): + resample_series(data, coordinates, [np.eye(4)], output_dtype=np.float32) diff --git a/petprep/workflows/pet/fit.py b/petprep/workflows/pet/fit.py index 09fc2d95..834a6040 100644 --- a/petprep/workflows/pet/fit.py +++ b/petprep/workflows/pet/fit.py @@ -276,6 +276,17 @@ def _write_identity_xforms(num_frames: int, filename: Path) -> Path: return filename +def _identity_xforms_path(pet_file: str, work_dir: Path) -> Path: + """Return a PET-specific identity transform path.""" + + pet_path = Path(pet_file) + pet_stem = pet_path + while pet_stem.suffix: + pet_stem = pet_stem.with_suffix('') + + return Path(work_dir) / f'{pet_stem.name}_idmat.tfm' + + def _construct_nu_path(subjects_dir: str, subject_id: str) -> str: """Return the expected path to FreeSurfer's ``nu.mgz`` for ``subject_id``.""" @@ -716,16 +727,15 @@ def init_pet_fit_wf( if hmc_disabled: config.execution.work_dir.mkdir(parents=True, exist_ok=True) 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) + idmat_fname = _identity_xforms_path(pet_file, config.execution.work_dir) + hmc_xforms = _write_identity_xforms(pet_tlen, idmat_fname) config.loggers.workflow.info('Head motion correction disabled; using identity transforms.') if petref_strategy == 'auto' and petref_candidates is not None: petref_candidates.inputs.template = petref if pet_tlen <= 1: # 3D PET petref = pet_file - idmat_fname = config.execution.work_dir / 'idmat.tfm' + idmat_fname = _identity_xforms_path(pet_file, config.execution.work_dir) hmc_xforms = _write_identity_xforms(pet_tlen, idmat_fname) config.loggers.workflow.debug('3D PET file - motion correction not needed') if petref: diff --git a/petprep/workflows/pet/tests/test_fit.py b/petprep/workflows/pet/tests/test_fit.py index 8dc58d60..e247451d 100644 --- a/petprep/workflows/pet/tests/test_fit.py +++ b/petprep/workflows/pet/tests/test_fit.py @@ -22,6 +22,7 @@ _extract_first5min_image, _extract_sum_image, _extract_twa_image, + _identity_xforms_path, _select_anatomical_reference, _select_best_petref, _write_identity_xforms, @@ -1495,7 +1496,7 @@ def _record_identity_xforms(num_frames, filename): assert not any(name.startswith('pet_hmc_wf') for name in wf.list_node_names()) hmc_buffer = wf.get_node('hmc_buffer') - assert str(hmc_buffer.inputs.hmc_xforms).endswith('idmat.tfm') + assert str(hmc_buffer.inputs.hmc_xforms).endswith('_idmat.tfm') assert Path(hmc_buffer.inputs.hmc_xforms).exists() assert identity_xform_frames == [data.shape[-1]] petref_buffer = wf.get_node('petref_buffer') @@ -1506,6 +1507,76 @@ def _record_identity_xforms(num_frames, filename): assert np.allclose(petref_img.get_fdata(), 14.0 / 6.0) +def test_pet_fit_hmc_off_identity_xforms_use_image_frame_count( + bids_root: Path, tmp_path: Path, monkeypatch +): + """Identity HMC transforms should follow image shape, not timing metadata length.""" + + pet_series = [str(bids_root / 'sub-01' / 'pet' / 'sub-01_task-rest_run-1_pet.nii.gz')] + data = np.zeros((2, 2, 2, 2), dtype=np.float32) + 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], "FrameDuration": [1]}') + + identity_xform_frames = [] + write_identity_xforms = pet_fit._write_identity_xforms + + def _record_identity_xforms(num_frames, filename): + identity_xform_frames.append(num_frames) + return write_identity_xforms(num_frames, filename) + + monkeypatch.setattr(pet_fit, '_write_identity_xforms', _record_identity_xforms) + + with mock_config(bids_dir=bids_root): + config.workflow.hmc_off = True + config.workflow.petref = 'sum' + init_pet_fit_wf(pet_series=pet_series, precomputed={}, omp_nthreads=1) + + assert identity_xform_frames == [data.shape[-1]] + + +def test_pet_fit_hmc_off_identity_xforms_are_pet_specific( + bids_root: Path, tmp_path: Path, monkeypatch +): + """Identity HMC transforms for one PET series should not overwrite another.""" + + pet_dir = bids_root / 'sub-01' / 'pet' + dynamic_pet = pet_dir / 'sub-01_task-rest_run-1_pet.nii.gz' + static_pet = pet_dir / 'sub-01_task-rest_run-2_pet.nii.gz' + nb.Nifti1Image(np.zeros((2, 2, 2, 32), dtype=np.float32), np.eye(4)).to_filename( + dynamic_pet + ) + nb.Nifti1Image(np.zeros((2, 2, 2), dtype=np.float32), np.eye(4)).to_filename(static_pet) + dynamic_pet.with_suffix('').with_suffix('.json').write_text( + json.dumps( + { + 'FrameTimesStart': list(range(32)), + 'FrameDuration': [1] * 32, + } + ) + ) + static_pet.with_suffix('').with_suffix('.json').write_text( + '{"FrameTimesStart": [0], "FrameDuration": [1]}' + ) + + with mock_config(bids_dir=bids_root): + config.workflow.hmc_off = True + dynamic_wf = init_pet_fit_wf(pet_series=[str(dynamic_pet)], precomputed={}, omp_nthreads=1) + static_wf = init_pet_fit_wf(pet_series=[str(static_pet)], precomputed={}, omp_nthreads=1) + + dynamic_xforms = Path(dynamic_wf.get_node('hmc_buffer').inputs.hmc_xforms) + static_xforms = Path(static_wf.get_node('hmc_buffer').inputs.hmc_xforms) + + assert dynamic_xforms != static_xforms + assert dynamic_xforms.name == 'sub-01_task-rest_run-1_pet_idmat.tfm' + assert static_xforms.name == 'sub-01_task-rest_run-2_pet_idmat.tfm' + assert np.asarray(nt.linear.load(dynamic_xforms).matrix).shape == (32, 4, 4) + assert np.asarray(nt.linear.load(static_xforms).matrix).shape == (4, 4) + + @pytest.mark.parametrize( ('frame_start_times', 'frame_durations', 'message'), [ @@ -1652,7 +1723,7 @@ def test_pet_fit_hmc_off_ignores_precomputed(bids_root: Path, tmp_path: Path): assert petref_buffer.inputs.petref != str(precomputed_petref) assert Path(petref_buffer.inputs.petref).name.endswith('_timeavgref.nii.gz') assert hmc_buffer.inputs.hmc_xforms != str(precomputed_hmc) - assert Path(hmc_buffer.inputs.hmc_xforms).name == 'idmat.tfm' + assert Path(hmc_buffer.inputs.hmc_xforms).name.endswith('_idmat.tfm') def test_pet_fit_picks_single_precomputed_derivative(bids_root: Path, tmp_path: Path): @@ -1711,6 +1782,16 @@ def test_write_identity_xforms_minimum(tmp_path: Path): assert np.allclose(matrices[0], np.eye(4)) +def test_identity_xforms_path_is_pet_specific(tmp_path: Path): + """Identity transform filenames should distinguish PET series.""" + + pet_file = '/data/sub-01/ses-closed/pet/sub-01_ses-closed_task-rest_pet.nii.gz' + + assert _identity_xforms_path(pet_file, tmp_path) == ( + tmp_path / 'sub-01_ses-closed_task-rest_pet_idmat.tfm' + ) + + def test_select_anatomical_reference_prefers_nu(tmp_path: Path): """Selecting ``anatref='nu'`` should return the FreeSurfer nu image when present.""" From d31993a94517fa249bd3a648450194cddd6301e5 Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Wed, 17 Jun 2026 16:06:24 +0200 Subject: [PATCH 2/3] FIX: style --- petprep/workflows/pet/tests/test_fit.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/petprep/workflows/pet/tests/test_fit.py b/petprep/workflows/pet/tests/test_fit.py index e247451d..c30bcb06 100644 --- a/petprep/workflows/pet/tests/test_fit.py +++ b/petprep/workflows/pet/tests/test_fit.py @@ -1546,9 +1546,7 @@ def test_pet_fit_hmc_off_identity_xforms_are_pet_specific( pet_dir = bids_root / 'sub-01' / 'pet' dynamic_pet = pet_dir / 'sub-01_task-rest_run-1_pet.nii.gz' static_pet = pet_dir / 'sub-01_task-rest_run-2_pet.nii.gz' - nb.Nifti1Image(np.zeros((2, 2, 2, 32), dtype=np.float32), np.eye(4)).to_filename( - dynamic_pet - ) + nb.Nifti1Image(np.zeros((2, 2, 2, 32), dtype=np.float32), np.eye(4)).to_filename(dynamic_pet) nb.Nifti1Image(np.zeros((2, 2, 2), dtype=np.float32), np.eye(4)).to_filename(static_pet) dynamic_pet.with_suffix('').with_suffix('.json').write_text( json.dumps( From 27ac758b0f7d0c1304e38de9f33cfde18f55ddb0 Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Wed, 17 Jun 2026 23:41:00 +0200 Subject: [PATCH 3/3] TST: Cover 3D resample_series path --- petprep/interfaces/tests/test_resampling.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/petprep/interfaces/tests/test_resampling.py b/petprep/interfaces/tests/test_resampling.py index 9cc0b985..23cd8d1d 100644 --- a/petprep/interfaces/tests/test_resampling.py +++ b/petprep/interfaces/tests/test_resampling.py @@ -4,6 +4,24 @@ from ..resampling import resample_series +def test_resample_series_handles_3d_input(): + """3D inputs should be resampled as a single volume.""" + + data = np.arange(8, dtype=np.float32).reshape((2, 2, 2)) + coordinates = np.indices(data.shape, dtype=np.float32) + + resampled = resample_series( + data, + coordinates, + [np.eye(4)], + output_dtype=np.float32, + order=0, + ) + + assert resampled.shape == data.shape + assert np.allclose(resampled, data) + + def test_resample_series_rejects_mismatched_hmc_transforms(): """HMC transform mappings must contain one affine per PET volume."""