From 808446ce487910240f4b2ee0369b0e89d5d8dd7f Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Tue, 9 Sep 2025 13:40:42 +0200 Subject: [PATCH 01/26] ENH: Add atlas field and update config --- petprep/cli/parser.py | 11 +++++++++++ petprep/cli/tests/test_parser.py | 18 ++++++++++++++++++ petprep/config.py | 3 ++- 3 files changed, 31 insertions(+), 1 deletion(-) diff --git a/petprep/cli/parser.py b/petprep/cli/parser.py index cee010d5..a291cb4b 100644 --- a/petprep/cli/parser.py +++ b/petprep/cli/parser.py @@ -554,6 +554,13 @@ def _bids_filter(value, parser): help='Segmentation method to use.', ) + g_seg.add_argument( + '--atlas', + action='store', + default=None, + help='Atlas to use for segmentation instead of FreeSurfer.', + ) + g_refmask = parser.add_argument_group('Options for reference mask generation') g_refmask.add_argument( '--ref-mask-name', @@ -727,6 +734,10 @@ def parse_args(args=None, namespace=None): config.execution.log_level = int(max(25 - 5 * opts.verbose_count, logging.DEBUG)) config.from_dict(vars(opts), init=['nipype']) + if opts.atlas is not None: + config.workflow.atlas = opts.atlas + config.workflow.seg = None + pvc_vals = (opts.pvc_tool, opts.pvc_method, opts.pvc_psf) if any(val is not None for val in pvc_vals) and not all(val is not None for val in pvc_vals): parser.error('Options --pvc-tool, --pvc-method and --pvc-psf must be used together.') diff --git a/petprep/cli/tests/test_parser.py b/petprep/cli/tests/test_parser.py index 6eeff043..590af686 100644 --- a/petprep/cli/tests/test_parser.py +++ b/petprep/cli/tests/test_parser.py @@ -158,6 +158,24 @@ def test_parse_args(tmp_path, minimal_bids): _reset_config() +def test_atlas_overrides_seg(tmp_path, minimal_bids): + """Ensure --atlas sets atlas and disables segmentation.""" + out_dir = tmp_path / 'out' + parse_args( + args=[ + str(minimal_bids), + str(out_dir), + 'participant', + '--skip-bids-validation', + '--atlas', + 'MyAtlas', + ] + ) + assert config.workflow.atlas == 'MyAtlas' + assert config.workflow.seg is None + _reset_config() + + def test_bids_filter_file(tmp_path, capsys): bids_path = tmp_path / 'data' out_path = tmp_path / 'out' diff --git a/petprep/config.py b/petprep/config.py index b7de78ba..00bda267 100644 --- a/petprep/config.py +++ b/petprep/config.py @@ -599,7 +599,8 @@ class workflow(_Config): seg = 'gtm' """Segmentation approach ('gtm', 'brainstem', 'thalamicNuclei', 'hippocampusAmygdala', 'wm', 'raphe', 'limbic').""" - + atlas: str | None = None + """Atlas defining regional segmentation. If set, FreeSurfer segmentation is skipped.""" pvc_tool: str | None = None """Tool used for partial volume correction.""" pvc_method: str | None = None From c46bc3d401e9d7cfac5904f0b8f9828a15c7afb0 Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Tue, 9 Sep 2025 13:47:22 +0200 Subject: [PATCH 02/26] ENH: Create config.json for atlas --- petprep/data/atlas/config.json | 8 ++++++++ 1 file changed, 8 insertions(+) create mode 100644 petprep/data/atlas/config.json diff --git a/petprep/data/atlas/config.json b/petprep/data/atlas/config.json new file mode 100644 index 00000000..63f0ae13 --- /dev/null +++ b/petprep/data/atlas/config.json @@ -0,0 +1,8 @@ +{ + "MIAL67ThalamicNuclei": { + "tpl": "MNI152NLin2009cAsym", + "template_file": "tpl-MNI152NLin2009cAsym_res-01_desc-brain_T1w.nii.gz", + "atlas_file": "tpl-MNI152NLin2009cAsym_res-01_atlas-MIAL67ThalamicNuclei_dseg.nii.gz", + "label_file": "tpl-MNI152NLin2009cAsym_atlas-MIAL67ThalamicNuclei_dseg.tsv" + } +} From 7cdc09e1cdbac73f46deef67fbb3effeddb9ea20 Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Tue, 9 Sep 2025 14:10:35 +0200 Subject: [PATCH 03/26] ENH: Add init_atlas_wf function in atlas.py --- petprep/workflows/pet/__init__.py | 2 + petprep/workflows/pet/atlas.py | 214 ++++++++++++++++++++++ petprep/workflows/pet/tests/test_atlas.py | 67 +++++++ 3 files changed, 283 insertions(+) create mode 100644 petprep/workflows/pet/atlas.py create mode 100644 petprep/workflows/pet/tests/test_atlas.py diff --git a/petprep/workflows/pet/__init__.py b/petprep/workflows/pet/__init__.py index e9470d18..a8aff7a4 100644 --- a/petprep/workflows/pet/__init__.py +++ b/petprep/workflows/pet/__init__.py @@ -21,6 +21,7 @@ from .registration import init_pet_reg_wf from .resampling import init_pet_surf_wf from .tacs import init_pet_tacs_wf +from .atlas import init_atlas_wf __all__ = [ 'init_pet_confs_wf', @@ -29,4 +30,5 @@ 'init_pet_surf_wf', 'init_pet_tacs_wf', 'init_pet_ref_tacs_wf', + 'init_atlas_wf', ] diff --git a/petprep/workflows/pet/atlas.py b/petprep/workflows/pet/atlas.py new file mode 100644 index 00000000..5f1bd43b --- /dev/null +++ b/petprep/workflows/pet/atlas.py @@ -0,0 +1,214 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +# +# Copyright The NiPreps Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# We support and encourage derived works from this project, please read +# about our expectations at +# +# https://www.nipreps.org/community/licensing/ +# +"""Atlas workflows.""" + +from __future__ import annotations + +import json +from pathlib import Path + +from nipype.interfaces import utility as niu +from nipype.interfaces.ants import Registration +from nipype.pipeline import engine as pe +from niworkflows.engine.workflows import LiterateWorkflow as Workflow +from niworkflows.interfaces.fixes import FixHeaderApplyTransforms as ApplyTransforms + +from ... import config +from ...interfaces import DerivativesDataSink +from ...interfaces.bids import BIDSURI +from templateflow.api import get as get_template + +DEFAULT_MEMORY_MIN_GB = config.DEFAULT_MEMORY_MIN_GB + + +def _atlas_morph_tsv(segmentation: str, labels_tsv: str) -> str: + """Generate a TSV table of region volumes from a segmentation.""" + import nibabel as nb + import numpy as np + import pandas as pd + + seg_img = nb.load(segmentation) + data = np.asanyarray(seg_img.dataobj) + voxvol = np.prod(seg_img.header.get_zooms()[:3]) + labels = pd.read_csv(labels_tsv, sep="\t") + volumes = [(data == int(idx)).sum() * voxvol for idx in labels["index"]] + out = labels.copy() + out["volume-mm3"] = volumes + out_file = Path("morph.tsv").absolute() + out.to_csv(out_file, sep="\t", index=False) + return str(out_file) + + +def init_atlas_wf(atlas: str, config_file: str, name: str = "pet_atlas_wf") -> Workflow: + """Map a template atlas into T1w space and compute regional volumes. + + Parameters + ---------- + atlas : :class:`str` + Name of atlas (used in outputs and TemplateFlow queries). + config_file : :class:`str` + JSON file with query parameters for TemplateFlow ``get`` calls. + name : :class:`str` + Workflow name (default: ``pet_atlas_wf``). + + Inputs + ------ + t1w_preproc + Preprocessed T1-weighted image. + + Outputs + ------- + segmentation + Atlas segmentation in T1w space. + dseg_tsv + Label table for the atlas. + """ + with open(config_file) as f: + cfg = json.load(f) + + template_t1w = str(get_template(**cfg["t1w"])) + atlas_img = str(get_template(atlas=atlas, **cfg["atlas"])) + labels_tsv = str(get_template(atlas=atlas, **cfg["labels"])) + + workflow = Workflow(name=name) + + inputnode = pe.Node(niu.IdentityInterface(fields=["t1w_preproc"]), name="inputnode") + outputnode = pe.Node( + niu.IdentityInterface(fields=["segmentation", "dseg_tsv"]), + name="outputnode", + ) + + label_source = pe.Node(niu.IdentityInterface(fields=["dseg_tsv"]), name="label_source") + label_source.inputs.dseg_tsv = labels_tsv + + reg = pe.Node( + Registration( + float=True, + output_inverse_warped_image=True, + output_warped_image=True, + transforms=["Rigid", "Affine"], + transform_parameters=[(0.1,), (0.1,)], + metric=["MI", "MI"], + metric_weight=[1, 1], + radius_or_number_of_bins=[32, 32], + sampling_strategy=["Regular", "Regular"], + sampling_percentage=[0.25, 0.25], + convergence_threshold=[1e-6, 1e-6], + convergence_window_size=[10, 10], + number_of_iterations=[[1000, 500, 250, 100]] * 2, + shrink_factors=[[8, 4, 2, 1]] * 2, + smoothing_sigmas=[[3, 2, 1, 0]] * 2, + use_histogram_matching=True, + ), + name="t1_to_tpl", + ) + reg.inputs.fixed_image = template_t1w + + apply_inv = pe.Node(ApplyTransforms(interpolation="NearestNeighbor"), name="apply_atlas") + apply_inv.inputs.input_image = atlas_img + + gen_morph = pe.Node( + niu.Function( + input_names=["segmentation", "labels_tsv"], + output_names=["out_file"], + function=_atlas_morph_tsv, + ), + name="gen_morph", + ) + gen_morph.inputs.labels_tsv = labels_tsv + + sources = pe.Node( + BIDSURI( + numinputs=1, + dataset_links=config.execution.dataset_links, + out_dir=str(config.execution.petprep_dir), + ), + name="sources", + ) + + ds_seg = pe.Node( + DerivativesDataSink( + base_directory=config.execution.petprep_dir, + seg=atlas, + allowed_entities=("seg",), + suffix="dseg", + extension=".nii.gz", + compress=True, + ), + name="ds_seg", + run_without_submitting=True, + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + ds_dseg_tsv = pe.Node( + DerivativesDataSink( + base_directory=config.execution.petprep_dir, + seg=atlas, + allowed_entities=("seg",), + suffix="dseg", + extension=".tsv", + datatype="anat", + check_hdr=False, + ), + name="ds_dseg_tsv", + run_without_submitting=True, + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + ds_morph_tsv = pe.Node( + DerivativesDataSink( + base_directory=config.execution.petprep_dir, + seg=atlas, + allowed_entities=("seg",), + suffix="morph", + extension=".tsv", + datatype="anat", + check_hdr=False, + ), + name="ds_morph_tsv", + run_without_submitting=True, + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + workflow.connect( + [ + (inputnode, reg, [("t1w_preproc", "moving_image")]), + (reg, apply_inv, [("inverse_composite_transform", "transforms")]), + (inputnode, apply_inv, [("t1w_preproc", "reference_image")]), + (apply_inv, outputnode, [("output_image", "segmentation")]), + (label_source, outputnode, [("dseg_tsv", "dseg_tsv")]), + (inputnode, sources, [("t1w_preproc", "in1")]), + (apply_inv, ds_seg, [("output_image", "in_file")]), + (inputnode, ds_seg, [("t1w_preproc", "source_file")]), + (sources, ds_seg, [("out", "Sources")]), + (label_source, ds_dseg_tsv, [("dseg_tsv", "in_file")]), + (inputnode, ds_dseg_tsv, [("t1w_preproc", "source_file")]), + (sources, ds_dseg_tsv, [("out", "Sources")]), + (apply_inv, gen_morph, [("output_image", "segmentation")]), + (gen_morph, ds_morph_tsv, [("out_file", "in_file")]), + (inputnode, ds_morph_tsv, [("t1w_preproc", "source_file")]), + (sources, ds_morph_tsv, [("out", "Sources")]), + ] + ) # fmt:skip + + return workflow diff --git a/petprep/workflows/pet/tests/test_atlas.py b/petprep/workflows/pet/tests/test_atlas.py new file mode 100644 index 00000000..088400f8 --- /dev/null +++ b/petprep/workflows/pet/tests/test_atlas.py @@ -0,0 +1,67 @@ +import json + +import nibabel as nb +import numpy as np +import pandas as pd + +from petprep import config +from ..atlas import init_atlas_wf, _atlas_morph_tsv + + +def test_init_atlas_wf_build(tmp_path, monkeypatch): + t1_img = nb.Nifti1Image(np.zeros((2, 2, 2)), np.eye(4)) + t1_file = tmp_path / 't1w.nii.gz' + t1_img.to_filename(t1_file) + + atlas_img = nb.Nifti1Image(np.zeros((2, 2, 2)), np.eye(4)) + atlas_file = tmp_path / 'atlas.nii.gz' + atlas_img.to_filename(atlas_file) + + labels_file = tmp_path / 'labels.tsv' + labels_file.write_text('index\tname\n0\tbackground\n') + + cfg = { + 't1w': {'suffix': 'T1w'}, + 'atlas': {'suffix': 'dseg'}, + 'labels': {'suffix': 'dseg', 'extension': '.tsv'}, + } + cfg_file = tmp_path / 'cfg.json' + cfg_file.write_text(json.dumps(cfg)) + + def fake_get(**kwargs): + if kwargs.get('suffix') == 'T1w': + return str(t1_file) + if kwargs.get('extension') == '.tsv': + return str(labels_file) + return str(atlas_file) + + monkeypatch.setattr('petprep.workflows.pet.atlas.get_template', fake_get) + + config.execution.petprep_dir = tmp_path + config.execution.dataset_links = {} + + wf = init_atlas_wf(atlas='foo', config_file=str(cfg_file)) + node_names = [n.name for n in wf._get_all_nodes()] + assert 'apply_atlas' in node_names + assert wf.get_node('label_source').inputs.dseg_tsv == str(labels_file) + assert wf.get_node('ds_seg').inputs.seg == 'foo' + assert wf.get_node('ds_dseg_tsv').inputs.seg == 'foo' + assert wf.get_node('ds_morph_tsv').inputs.seg == 'foo' + + +def test_atlas_morph_tsv(tmp_path): + data = np.array([ + [[0, 1], [1, 1]], + [[0, 0], [1, 0]], + ], dtype='int16') + seg = nb.Nifti1Image(data, np.eye(4)) + seg_file = tmp_path / 'seg.nii.gz' + seg.to_filename(seg_file) + + labels = tmp_path / 'labels.tsv' + labels.write_text('index\tname\n0\tbg\n1\tregion\n') + + out = _atlas_morph_tsv(str(seg_file), str(labels)) + df = pd.read_csv(out, sep='\t') + volume = df.loc[df['index'] == 1, 'volume-mm3'].iloc[0] + assert volume == 4 From 0df7a349932021fa93bfb68c288afa01acf489ad Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Tue, 9 Sep 2025 14:27:45 +0200 Subject: [PATCH 04/26] ENH: Update workflow logic for atlas support --- petprep/workflows/base.py | 62 +++++++++++++++++++--------- petprep/workflows/tests/test_base.py | 48 +++++++++++++++++++++ 2 files changed, 91 insertions(+), 19 deletions(-) diff --git a/petprep/workflows/base.py b/petprep/workflows/base.py index f71c6f7d..1e74b79d 100644 --- a/petprep/workflows/base.py +++ b/petprep/workflows/base.py @@ -164,6 +164,7 @@ def init_single_subject_wf(subject_id: str): ) from petprep.workflows.pet.base import init_pet_wf + from petprep.workflows.pet import init_atlas_wf from petprep.workflows.pet.segmentation import init_segmentation_wf workflow = Workflow(name=f'sub_{subject_id}_wf') @@ -524,29 +525,52 @@ def init_single_subject_wf(subject_id: str): ]), ]) # fmt:skip - segmentation_wf = init_segmentation_wf( - seg=config.workflow.seg, - name=f'pet_{config.workflow.seg}_seg_wf', - ) - workflow.connect( - [ - ( - anat_fit_wf, - segmentation_wf, - [ - ('outputnode.t1w_preproc', 'inputnode.t1w_preproc'), - ('outputnode.subjects_dir', 'inputnode.subjects_dir'), - ('outputnode.subject_id', 'inputnode.subject_id'), - ], - ), - ] - ) + if config.workflow.atlas: + try: # Py>=3.9 + from importlib.resources import files as ir_files + except Exception: # pragma: no cover - Py<3.9 fallback + from importlib_resources import files as ir_files + + seg_wf = init_atlas_wf( + atlas=config.workflow.atlas, + config_file=str(ir_files('petprep.data.atlas') / 'config.json'), + name=f'pet_{config.workflow.atlas}_atlas_wf', + ) + workflow.connect( + [ + (anat_fit_wf, seg_wf, [('outputnode.t1w_preproc', 'inputnode.t1w_preproc')]), + ] + ) + else: + seg_wf = init_segmentation_wf( + seg=config.workflow.seg, + name=f'pet_{config.workflow.seg}_seg_wf', + ) + workflow.connect( + [ + ( + anat_fit_wf, + seg_wf, + [ + ('outputnode.t1w_preproc', 'inputnode.t1w_preproc'), + ('outputnode.subjects_dir', 'inputnode.subjects_dir'), + ('outputnode.subject_id', 'inputnode.subject_id'), + ], + ), + ] + ) if config.workflow.anat_only: return clean_datasinks(workflow) # Append the PET section to the existing anatomical excerpt # That way we do not need to filter down the number of PET datasets + seg_ref = ( + f'``{config.workflow.atlas}`` atlas' + if config.workflow.atlas + else f'``{config.workflow.seg}`` segmentation workflow from FreeSurfer' + ) + pet_pre_desc = f""" PET data preprocessing @@ -555,7 +579,7 @@ def init_single_subject_wf(subject_id: str): motion estimation and correction were carried out after generating a reference image, which was subsequently coregistered to the T1-weighted anatomical image. A brain mask was computed and the structural image was -segmented with the ``{config.workflow.seg}`` segmentation workflow from FreeSurfer.""" +segmented with the {seg_ref}.""" if config.workflow.pvc_tool and config.workflow.pvc_method: pet_pre_desc += ( @@ -619,7 +643,7 @@ def init_single_subject_wf(subject_id: str): 'inputnode.sphere_reg_fsLR', ), ]), - (segmentation_wf, pet_wf, [ + (seg_wf, pet_wf, [ ('outputnode.segmentation', 'inputnode.segmentation'), ('outputnode.dseg_tsv', 'inputnode.dseg_tsv'), ]), diff --git a/petprep/workflows/tests/test_base.py b/petprep/workflows/tests/test_base.py index a1eeca35..9eb8b218 100644 --- a/petprep/workflows/tests/test_base.py +++ b/petprep/workflows/tests/test_base.py @@ -143,6 +143,54 @@ def test_segmentation_shared_across_runs(multisession_bids_root): assert all('_seg_wf' not in n for n in pet_node.list_node_names()) +def test_atlas_replaces_segmentation(monkeypatch, multisession_bids_root): + def _dummy_atlas_wf(atlas, config_file, name='pet_atlas_wf'): + from nipype.interfaces import utility as niu + from nipype.pipeline import engine as pe + from niworkflows.engine.workflows import LiterateWorkflow as Workflow + + wf = Workflow(name=name) + inputnode = pe.Node(niu.IdentityInterface(fields=['t1w_preproc']), name='inputnode') + outputnode = pe.Node( + niu.IdentityInterface(fields=['segmentation', 'dseg_tsv']), + name='outputnode', + ) + wf.add_nodes([inputnode, outputnode]) + return wf + + monkeypatch.setattr('petprep.workflows.pet.init_atlas_wf', _dummy_atlas_wf) + + with mock_config(bids_dir=multisession_bids_root): + config.workflow.atlas = 'foo' + wf = init_single_subject_wf('01') + + flatgraph = wf._create_flat_graph() + generate_expanded_graph(flatgraph) + + atlas_wf_name = f'pet_{config.workflow.atlas}_atlas_wf' + atlas_nodes = [n for n in wf.list_node_names() if n.startswith(atlas_wf_name)] + assert atlas_nodes + + pet_wf_names = [ + n + for n in {name.split('.')[0] for name in wf.list_node_names() if name.startswith('pet_')} + if n != atlas_wf_name + ] + assert len(pet_wf_names) == 2 + + atlas_node = wf.get_node(atlas_wf_name) + for name in pet_wf_names: + pet_node = wf.get_node(name) + edge = wf._graph.get_edge_data(atlas_node, pet_node) + assert ('outputnode.segmentation', 'inputnode.segmentation') in edge['connect'] + assert ('outputnode.dseg_tsv', 'inputnode.dseg_tsv') in edge['connect'] + assert all('_atlas_wf' not in n for n in pet_node.list_node_names()) + + pet_node = wf.get_node(pet_wf_names[0]) + assert config.workflow.atlas in pet_node.__desc__ + assert config.workflow.seg not in pet_node.__desc__ + + def _make_params( pet2anat_init: str = 'auto', medial_surface_nan: bool = False, From adc448fd42037799986cecd23f312f0910b25497 Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Tue, 9 Sep 2025 14:48:27 +0200 Subject: [PATCH 05/26] ENH: Update DerivativesDataSink instantiation --- petprep/workflows/pet/base.py | 19 ++++++++++++++---- petprep/workflows/pet/tests/test_base.py | 25 ++++++++++++++++++++++++ 2 files changed, 40 insertions(+), 4 deletions(-) diff --git a/petprep/workflows/pet/base.py b/petprep/workflows/pet/base.py index d40890e7..6fa30711 100644 --- a/petprep/workflows/pet/base.py +++ b/petprep/workflows/pet/base.py @@ -692,14 +692,19 @@ def init_pet_wf( Path(pet_file).with_suffix('').with_suffix('.json') ) + entity = ( + ('seg', config.workflow.atlas) + if config.workflow.atlas + else ('seg', config.workflow.seg) + ) ds_pet_tacs = pe.Node( DerivativesDataSink( base_directory=petprep_dir, suffix='tacs', - seg=config.workflow.seg, desc='preproc', - allowed_entities=('seg',), + allowed_entities=(entity[0],), TaskName=all_metadata[0].get('TaskName'), + **{entity[0]: entity[1]}, **prepare_timing_parameters(all_metadata[0]), ), name='ds_pet_tacs', @@ -726,15 +731,21 @@ def init_pet_wf( ) pet_ref_tacs_wf.inputs.inputnode.ref_mask_name = config.workflow.ref_mask_name + entity = ( + ('seg', config.workflow.atlas) + if config.workflow.atlas + else ('seg', config.workflow.seg) + ) + ds_ref_tacs = pe.Node( DerivativesDataSink( base_directory=petprep_dir, suffix='tacs', - seg=config.workflow.seg, desc='preproc', ref=config.workflow.ref_mask_name, - allowed_entities=('seg', 'ref'), + allowed_entities=(entity[0], 'ref'), TaskName=all_metadata[0].get('TaskName'), + **{entity[0]: entity[1]}, **prepare_timing_parameters(all_metadata[0]), ), name='ds_ref_tacs', diff --git a/petprep/workflows/pet/tests/test_base.py b/petprep/workflows/pet/tests/test_base.py index 97d0ec05..0042ed8d 100644 --- a/petprep/workflows/pet/tests/test_base.py +++ b/petprep/workflows/pet/tests/test_base.py @@ -254,6 +254,31 @@ def test_pet_ref_tacs_wf_connections(bids_root: Path): assert ('outputnode.timeseries', 'in_file') in edge_ds['connect'] +def test_seg_entity_freesurfer(bids_root: Path): + """TACs datasinks use FreeSurfer segmentation metadata by default.""" + pet_series = _prep_pet_series(bids_root) + + with mock_config(bids_dir=bids_root): + config.workflow.ref_mask_name = 'cerebellum' + wf = init_pet_wf(pet_series=pet_series, precomputed={}) + + assert wf.get_node('ds_pet_tacs').inputs.seg == config.workflow.seg + assert wf.get_node('ds_ref_tacs').inputs.seg == config.workflow.seg + + +def test_seg_entity_atlas(bids_root: Path): + """TACs datasinks use atlas metadata when configured.""" + pet_series = _prep_pet_series(bids_root) + + with mock_config(bids_dir=bids_root): + config.workflow.atlas = 'foo' + config.workflow.ref_mask_name = 'cerebellum' + wf = init_pet_wf(pet_series=pet_series, precomputed={}) + + assert wf.get_node('ds_pet_tacs').inputs.seg == config.workflow.atlas + assert wf.get_node('ds_ref_tacs').inputs.seg == config.workflow.atlas + + def test_psf_metadata_propagation(bids_root: Path): """PSF values should be passed to datasinks when using AGTM.""" pet_series = _prep_pet_series(bids_root) From 67076b5451f210ddfc7534c631d0136b821fba45 Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Tue, 9 Sep 2025 14:59:15 +0200 Subject: [PATCH 06/26] FIX: Update fit.py for entity handling --- petprep/workflows/pet/fit.py | 12 ++++++--- petprep/workflows/pet/tests/test_fit.py | 36 +++++++++++++++++++++++++ 2 files changed, 45 insertions(+), 3 deletions(-) diff --git a/petprep/workflows/pet/fit.py b/petprep/workflows/pet/fit.py index e5f06107..09d67e1c 100644 --- a/petprep/workflows/pet/fit.py +++ b/petprep/workflows/pet/fit.py @@ -408,8 +408,14 @@ def init_pet_fit_wf( f'PET Stage 4: Generating {config.workflow.ref_mask_name} reference mask' ) + entity = ( + ('seg', config.workflow.atlas) + if config.workflow.atlas + else ('seg', config.workflow.seg) + ) + refmask_wf = init_pet_refmask_wf( - segmentation=config.workflow.seg, + segmentation=entity[1], ref_mask_name=config.workflow.ref_mask_name, ref_mask_index=list(config.workflow.ref_mask_index) if config.workflow.ref_mask_index @@ -443,11 +449,11 @@ def init_pet_fit_wf( DerivativesDataSink( base_directory=config.execution.petprep_dir, suffix='tacs', - seg=config.workflow.seg, desc='preproc', ref=config.workflow.ref_mask_name, - allowed_entities=('seg', 'ref'), + allowed_entities=(entity[0], 'ref'), TaskName=metadata.get('TaskName'), + **{entity[0]: entity[1]}, **timing_parameters, ), name='ds_ref_tacs', diff --git a/petprep/workflows/pet/tests/test_fit.py b/petprep/workflows/pet/tests/test_fit.py index c7256170..2958a31c 100644 --- a/petprep/workflows/pet/tests/test_fit.py +++ b/petprep/workflows/pet/tests/test_fit.py @@ -255,6 +255,42 @@ def test_refmask_report_connections(bids_root: Path, tmp_path: Path, pvc_method) assert 'ds_ref_tacs' not in wf.list_node_names() +def test_seg_entity_atlas(bids_root: Path, tmp_path: Path): + """Reference outputs use atlas metadata when configured.""" + 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) + + sidecar = Path(pet_series[0]).with_suffix('').with_suffix('.json') + sidecar.write_text('{"FrameTimesStart": [0], "FrameDuration": [1]}') + + dummy_ref = str(tmp_path / 'dummy.nii') + dummy_xfm = str(tmp_path / 'dummy.txt') + img.to_filename(dummy_ref) + np.savetxt(dummy_xfm, np.eye(4)) + precomputed = { + 'petref': dummy_ref, + 'transforms': {'hmc': dummy_xfm, 'petref2anat': dummy_xfm}, + } + + with mock_config(bids_dir=bids_root): + config.workflow.atlas = 'foo' + config.workflow.ref_mask_name = 'cerebellum' + wf = init_pet_fit_wf( + pet_series=pet_series, + precomputed=precomputed, + omp_nthreads=1, + ) + + extract = wf.get_node('pet_refmask_wf').get_node('extract_refregion') + assert extract.inputs.segmentation_type == config.workflow.atlas + ds_tacs = wf.get_node('ds_ref_tacs') + assert ds_tacs.inputs.seg == config.workflow.atlas + edge = wf._graph.get_edge_data(wf.get_node('inputnode'), wf.get_node('pet_refmask_wf')) + assert ('segmentation', 'inputnode.seg_file') in edge['connect'] + + def test_pet_fit_stage1_inclusion(bids_root: Path, tmp_path: Path): """Stage 1 should run only when HMC derivatives are missing.""" pet_series = [str(bids_root / 'sub-01' / 'pet' / 'sub-01_task-rest_run-1_pet.nii.gz')] From bc9bb1d468e207986e3fe9dd1d392610e0aca22b Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Tue, 9 Sep 2025 15:03:02 +0200 Subject: [PATCH 07/26] FIX: Update atlas config and add tests --- petprep/data/atlas/config.json | 24 ++++++++++++++--- petprep/workflows/pet/atlas.py | 32 +++++++++++++++++++---- petprep/workflows/pet/tests/test_atlas.py | 31 ++++++++++++---------- 3 files changed, 64 insertions(+), 23 deletions(-) diff --git a/petprep/data/atlas/config.json b/petprep/data/atlas/config.json index 63f0ae13..8d5e15a5 100644 --- a/petprep/data/atlas/config.json +++ b/petprep/data/atlas/config.json @@ -1,8 +1,24 @@ { "MIAL67ThalamicNuclei": { - "tpl": "MNI152NLin2009cAsym", - "template_file": "tpl-MNI152NLin2009cAsym_res-01_desc-brain_T1w.nii.gz", - "atlas_file": "tpl-MNI152NLin2009cAsym_res-01_atlas-MIAL67ThalamicNuclei_dseg.nii.gz", - "label_file": "tpl-MNI152NLin2009cAsym_atlas-MIAL67ThalamicNuclei_dseg.tsv" + "t1w": { + "tpl": "MNI152NLin2009cAsym", + "desc": "brain", + "suffix": "T1w", + "res": "01", + "extension": "nii.gz" + }, + "atlas": { + "tpl": "MNI152NLin2009cAsym", + "atlas": "MIAL67ThalamicNuclei", + "suffix": "dseg", + "res": "01", + "extension": "nii.gz" + }, + "labels": { + "tpl": "MNI152NLin2009cAsym", + "atlas": "MIAL67ThalamicNuclei", + "suffix": "dseg", + "extension": "tsv" + } } } diff --git a/petprep/workflows/pet/atlas.py b/petprep/workflows/pet/atlas.py index 5f1bd43b..2494fcaa 100644 --- a/petprep/workflows/pet/atlas.py +++ b/petprep/workflows/pet/atlas.py @@ -84,11 +84,33 @@ def init_atlas_wf(atlas: str, config_file: str, name: str = "pet_atlas_wf") -> W Label table for the atlas. """ with open(config_file) as f: - cfg = json.load(f) - - template_t1w = str(get_template(**cfg["t1w"])) - atlas_img = str(get_template(atlas=atlas, **cfg["atlas"])) - labels_tsv = str(get_template(atlas=atlas, **cfg["labels"])) + data = json.load(f) + + if atlas not in data: + raise ValueError( + f"Atlas '{atlas}' not found in {config_file}. " + f"Available atlases: {', '.join(sorted(data.keys()))}" + ) + + cfg = data[atlas] + required = {"t1w", "atlas", "labels"} + missing = required - cfg.keys() + if missing: + raise ValueError( + f"Atlas '{atlas}' missing required keys: {', '.join(sorted(missing))}" + ) + + def _tf_kwargs(d: dict) -> dict: + out = dict(d) + if "tpl" in out: + out["template"] = out.pop("tpl") + if "res" in out: + out["resolution"] = out.pop("res") + return out + + template_t1w = str(get_template(**_tf_kwargs(cfg["t1w"]))) + atlas_img = str(get_template(**_tf_kwargs(cfg["atlas"]))) + labels_tsv = str(get_template(**_tf_kwargs(cfg["labels"]))) workflow = Workflow(name=name) diff --git a/petprep/workflows/pet/tests/test_atlas.py b/petprep/workflows/pet/tests/test_atlas.py index 088400f8..2506cc3b 100644 --- a/petprep/workflows/pet/tests/test_atlas.py +++ b/petprep/workflows/pet/tests/test_atlas.py @@ -1,8 +1,8 @@ -import json - import nibabel as nb import numpy as np import pandas as pd +import pytest +from importlib.resources import files as ir_files from petprep import config from ..atlas import init_atlas_wf, _atlas_morph_tsv @@ -20,18 +20,12 @@ def test_init_atlas_wf_build(tmp_path, monkeypatch): labels_file = tmp_path / 'labels.tsv' labels_file.write_text('index\tname\n0\tbackground\n') - cfg = { - 't1w': {'suffix': 'T1w'}, - 'atlas': {'suffix': 'dseg'}, - 'labels': {'suffix': 'dseg', 'extension': '.tsv'}, - } - cfg_file = tmp_path / 'cfg.json' - cfg_file.write_text(json.dumps(cfg)) + cfg_file = ir_files('petprep.data.atlas') / 'config.json' def fake_get(**kwargs): if kwargs.get('suffix') == 'T1w': return str(t1_file) - if kwargs.get('extension') == '.tsv': + if kwargs.get('extension') == 'tsv': return str(labels_file) return str(atlas_file) @@ -40,13 +34,22 @@ def fake_get(**kwargs): config.execution.petprep_dir = tmp_path config.execution.dataset_links = {} - wf = init_atlas_wf(atlas='foo', config_file=str(cfg_file)) + wf = init_atlas_wf( + atlas='MIAL67ThalamicNuclei', + config_file=str(cfg_file), + ) node_names = [n.name for n in wf._get_all_nodes()] assert 'apply_atlas' in node_names assert wf.get_node('label_source').inputs.dseg_tsv == str(labels_file) - assert wf.get_node('ds_seg').inputs.seg == 'foo' - assert wf.get_node('ds_dseg_tsv').inputs.seg == 'foo' - assert wf.get_node('ds_morph_tsv').inputs.seg == 'foo' + assert wf.get_node('ds_seg').inputs.seg == 'MIAL67ThalamicNuclei' + assert wf.get_node('ds_dseg_tsv').inputs.seg == 'MIAL67ThalamicNuclei' + assert wf.get_node('ds_morph_tsv').inputs.seg == 'MIAL67ThalamicNuclei' + + +def test_init_atlas_wf_bad_name(tmp_path): + cfg_file = ir_files('petprep.data.atlas') / 'config.json' + with pytest.raises(ValueError, match="not found"): + init_atlas_wf(atlas='notreal', config_file=str(cfg_file)) def test_atlas_morph_tsv(tmp_path): From ec1569d85ea3cb0b601d7bc22b1a76690681a949 Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Tue, 9 Sep 2025 15:51:35 +0200 Subject: [PATCH 08/26] FIX: Update atlas configuration in config.json --- petprep/data/atlas/config.json | 653 ++++++++++++++++++++++++++++++++- 1 file changed, 644 insertions(+), 9 deletions(-) diff --git a/petprep/data/atlas/config.json b/petprep/data/atlas/config.json index 8d5e15a5..22eefb3e 100644 --- a/petprep/data/atlas/config.json +++ b/petprep/data/atlas/config.json @@ -1,24 +1,659 @@ { - "MIAL67ThalamicNuclei": { - "t1w": { - "tpl": "MNI152NLin2009cAsym", + "DKT31": { + "atlas": { + "desc": "DKT31", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "desc": "DKT31", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "HOCPALth0": { + "atlas": { + "atlas": "HOCPAL", + "desc": "th0", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "HOCPAL", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin6Asym" + }, + "t1w": { "desc": "brain", + "extension": "nii.gz", + "res": "01", "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "HOCPALth25": { + "atlas": { + "atlas": "HOCPAL", + "desc": "th25", + "extension": "nii.gz", "res": "01", - "extension": "nii.gz" + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "HOCPAL", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin6Asym" }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "HOCPALth50": { "atlas": { - "tpl": "MNI152NLin2009cAsym", - "atlas": "MIAL67ThalamicNuclei", + "atlas": "HOCPAL", + "desc": "th50", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "HOCPAL", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin6Asym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "HOCPAth0": { + "atlas": { + "atlas": "HOCPA", + "desc": "th0", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "HOCPA", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin6Asym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "HOCPAth25": { + "atlas": { + "atlas": "HOCPA", + "desc": "th25", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "HOCPA", + "extension": "tsv", "suffix": "dseg", + "tpl": "MNI152NLin6Asym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "HOCPAth50": { + "atlas": { + "atlas": "HOCPA", + "desc": "th50", + "extension": "nii.gz", "res": "01", - "extension": "nii.gz" + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" }, "labels": { - "tpl": "MNI152NLin2009cAsym", + "atlas": "HOCPA", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin6Asym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "HOSPAth0": { + "atlas": { + "atlas": "HOSPA", + "desc": "th0", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "HOSPA", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin6Asym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "HOSPAth25": { + "atlas": { + "atlas": "HOSPA", + "desc": "th25", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "HOSPA", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin6Asym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "HOSPAth50": { + "atlas": { + "atlas": "HOSPA", + "desc": "th50", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "HOSPA", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin6Asym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "MIAL67ThalamicNuclei": { + "atlas": { "atlas": "MIAL67ThalamicNuclei", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "MIAL67ThalamicNuclei", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "Schaefer20181000Parcels17Networks": { + "atlas": { + "atlas": "Schaefer2018", + "desc": "1000Parcels17Networks", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "Schaefer2018", + "desc": "1000Parcels17Networks", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "Schaefer20181000Parcels7Networks": { + "atlas": { + "atlas": "Schaefer2018", + "desc": "1000Parcels7Networks", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "Schaefer2018", + "desc": "1000Parcels7Networks", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "Schaefer2018100Parcels17Networks": { + "atlas": { + "atlas": "Schaefer2018", + "desc": "100Parcels17Networks", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "Schaefer2018", + "desc": "100Parcels17Networks", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "Schaefer2018100Parcels7Networks": { + "atlas": { + "atlas": "Schaefer2018", + "desc": "100Parcels7Networks", + "extension": "nii.gz", + "res": "01", "suffix": "dseg", - "extension": "tsv" + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "Schaefer2018", + "desc": "100Parcels7Networks", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "Schaefer2018200Parcels17Networks": { + "atlas": { + "atlas": "Schaefer2018", + "desc": "200Parcels17Networks", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "Schaefer2018", + "desc": "200Parcels17Networks", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "Schaefer2018200Parcels7Networks": { + "atlas": { + "atlas": "Schaefer2018", + "desc": "200Parcels7Networks", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "Schaefer2018", + "desc": "200Parcels7Networks", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "Schaefer2018300Parcels17Networks": { + "atlas": { + "atlas": "Schaefer2018", + "desc": "300Parcels17Networks", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "Schaefer2018", + "desc": "300Parcels17Networks", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "Schaefer2018300Parcels7Networks": { + "atlas": { + "atlas": "Schaefer2018", + "desc": "300Parcels7Networks", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "Schaefer2018", + "desc": "300Parcels7Networks", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "Schaefer2018400Parcels17Networks": { + "atlas": { + "atlas": "Schaefer2018", + "desc": "400Parcels17Networks", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "Schaefer2018", + "desc": "400Parcels17Networks", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "Schaefer2018400Parcels7Networks": { + "atlas": { + "atlas": "Schaefer2018", + "desc": "400Parcels7Networks", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "Schaefer2018", + "desc": "400Parcels7Networks", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "Schaefer2018500Parcels17Networks": { + "atlas": { + "atlas": "Schaefer2018", + "desc": "500Parcels17Networks", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "Schaefer2018", + "desc": "500Parcels17Networks", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "Schaefer2018500Parcels7Networks": { + "atlas": { + "atlas": "Schaefer2018", + "desc": "500Parcels7Networks", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "Schaefer2018", + "desc": "500Parcels7Networks", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "Schaefer2018600Parcels17Networks": { + "atlas": { + "atlas": "Schaefer2018", + "desc": "600Parcels17Networks", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "Schaefer2018", + "desc": "600Parcels17Networks", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "Schaefer2018600Parcels7Networks": { + "atlas": { + "atlas": "Schaefer2018", + "desc": "600Parcels7Networks", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "Schaefer2018", + "desc": "600Parcels7Networks", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "Schaefer2018800Parcels17Networks": { + "atlas": { + "atlas": "Schaefer2018", + "desc": "800Parcels17Networks", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "Schaefer2018", + "desc": "800Parcels17Networks", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "Schaefer2018800Parcels7Networks": { + "atlas": { + "atlas": "Schaefer2018", + "desc": "800Parcels7Networks", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "atlas": "Schaefer2018", + "desc": "800Parcels7Networks", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" + } + }, + "carpet": { + "atlas": { + "desc": "carpet", + "extension": "nii.gz", + "res": "01", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "labels": { + "desc": "carpet", + "extension": "tsv", + "suffix": "dseg", + "tpl": "MNI152NLin2009cAsym" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "MNI152NLin2009cAsym" } } } From 386997a2e1bafea7895442916d8a3f50bc114c71 Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Tue, 9 Sep 2025 16:06:46 +0200 Subject: [PATCH 09/26] FIX: Add write_composite_transform to Registration --- petprep/workflows/pet/atlas.py | 1 + 1 file changed, 1 insertion(+) diff --git a/petprep/workflows/pet/atlas.py b/petprep/workflows/pet/atlas.py index 2494fcaa..6d25d647 100644 --- a/petprep/workflows/pet/atlas.py +++ b/petprep/workflows/pet/atlas.py @@ -141,6 +141,7 @@ def _tf_kwargs(d: dict) -> dict: shrink_factors=[[8, 4, 2, 1]] * 2, smoothing_sigmas=[[3, 2, 1, 0]] * 2, use_histogram_matching=True, + write_composite_transform=True, ), name="t1_to_tpl", ) From 98fb2e587f73abbf5bb6ce3f95f6fe52c6730d47 Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Tue, 9 Sep 2025 16:07:13 +0200 Subject: [PATCH 10/26] FIX: Edit config.json for SUIT atlases --- petprep/data/atlas/config.json | 67 ++++++++++++++++++++++++++++++++++ 1 file changed, 67 insertions(+) diff --git a/petprep/data/atlas/config.json b/petprep/data/atlas/config.json index 22eefb3e..33f64651 100644 --- a/petprep/data/atlas/config.json +++ b/petprep/data/atlas/config.json @@ -1,4 +1,71 @@ { + "Buckner201117n": { + "atlas": { + "atlas": "Buckner2011", + "seg": "17n", + "extension": "nii.gz", + "suffix": "dseg", + "tpl": "SUIT" + }, + "labels": { + "atlas": "Buckner2011", + "seg": "17n", + "extension": "tsv", + "suffix": "dseg", + "tpl": "SUIT" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "SUIT" + } + }, + "Buckner20117n": { + "atlas": { + "atlas": "Buckner2011", + "seg": "7n", + "extension": "nii.gz", + "suffix": "dseg", + "tpl": "SUIT" + }, + "labels": { + "atlas": "Buckner2011", + "seg": "7n", + "extension": "tsv", + "suffix": "dseg", + "tpl": "SUIT" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "SUIT" + } + }, + "Diedrichsen2009": { + "atlas": { + "atlas": "Diedrichsen2009", + "extension": "nii.gz", + "suffix": "dseg", + "tpl": "SUIT" + }, + "labels": { + "atlas": "Diedrichsen2009", + "extension": "tsv", + "suffix": "dseg", + "tpl": "SUIT" + }, + "t1w": { + "desc": "brain", + "extension": "nii.gz", + "res": "01", + "suffix": "T1w", + "tpl": "SUIT" + } + }, "DKT31": { "atlas": { "desc": "DKT31", From bf1683e593279d1bde73d36fe15ccfb7ba9b5c31 Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Tue, 9 Sep 2025 16:08:20 +0200 Subject: [PATCH 11/26] FIX: Update fetch_templates.py to implement fetch_SUIT --- scripts/fetch_templates.py | 39 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/scripts/fetch_templates.py b/scripts/fetch_templates.py index a7bec82e..97f9af7c 100755 --- a/scripts/fetch_templates.py +++ b/scripts/fetch_templates.py @@ -109,6 +109,45 @@ def fetch_fsLR(): tf.get('fsLR', density='32k') +def fetch_SUIT(): + """ + Expected templates: + + tpl-SUIT/tpl-SUIT_T1w.nii.gz + tpl-SUIT/tpl-SUIT_atlas-Diedrichsen2009_dseg.nii.gz + tpl-SUIT/tpl-SUIT_atlas-Diedrichsen2009_dseg.tsv + tpl-SUIT/tpl-SUIT_atlas-Buckner2011_seg-7n_dseg.nii.gz + tpl-SUIT/tpl-SUIT_atlas-Buckner2011_seg-7n_dseg.tsv + tpl-SUIT/tpl-SUIT_atlas-Buckner2011_seg-17n_dseg.nii.gz + tpl-SUIT/tpl-SUIT_atlas-Buckner2011_seg-17n_dseg.tsv + """ + template = 'SUIT' + + tf.get(template, suffix='T1w') + tf.get(template, atlas='Diedrichsen2009', suffix='dseg') + tf.get(template, atlas='Diedrichsen2009', suffix='dseg', extension='tsv') + + tf.get(template, atlas='Buckner2011', seg='7n', suffix='dseg', invalid_filters='allow') + tf.get( + template, + atlas='Buckner2011', + seg='7n', + suffix='dseg', + extension='tsv', + invalid_filters='allow', + ) + + tf.get(template, atlas='Buckner2011', seg='17n', suffix='dseg', invalid_filters='allow') + tf.get( + template, + atlas='Buckner2011', + seg='17n', + suffix='dseg', + extension='tsv', + invalid_filters='allow', + ) + + def fetch_all(): fetch_MNI2009() fetch_MNI6() From 815eecebc2aef5ddc39d01e0761f0be6400bee5b Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Tue, 9 Sep 2025 16:36:33 +0200 Subject: [PATCH 12/26] FIX: Edit _atlas_morph_tsv to add Path import + update registration settings --- petprep/workflows/pet/atlas.py | 40 ++++++++++++++++++++-------------- 1 file changed, 24 insertions(+), 16 deletions(-) diff --git a/petprep/workflows/pet/atlas.py b/petprep/workflows/pet/atlas.py index 6d25d647..be801af9 100644 --- a/petprep/workflows/pet/atlas.py +++ b/petprep/workflows/pet/atlas.py @@ -25,7 +25,6 @@ from __future__ import annotations import json -from pathlib import Path from nipype.interfaces import utility as niu from nipype.interfaces.ants import Registration @@ -43,6 +42,7 @@ def _atlas_morph_tsv(segmentation: str, labels_tsv: str) -> str: """Generate a TSV table of region volumes from a segmentation.""" + from pathlib import Path import nibabel as nb import numpy as np import pandas as pd @@ -125,21 +125,29 @@ def _tf_kwargs(d: dict) -> dict: reg = pe.Node( Registration( - float=True, - output_inverse_warped_image=True, - output_warped_image=True, - transforms=["Rigid", "Affine"], - transform_parameters=[(0.1,), (0.1,)], - metric=["MI", "MI"], - metric_weight=[1, 1], - radius_or_number_of_bins=[32, 32], - sampling_strategy=["Regular", "Regular"], - sampling_percentage=[0.25, 0.25], - convergence_threshold=[1e-6, 1e-6], - convergence_window_size=[10, 10], - number_of_iterations=[[1000, 500, 250, 100]] * 2, - shrink_factors=[[8, 4, 2, 1]] * 2, - smoothing_sigmas=[[3, 2, 1, 0]] * 2, + transforms=['Rigid', 'Affine', 'SyN'], + transform_parameters=[(0.1,), (0.1,), (0.1, 3, 0)], + metric=['Mattes', 'Mattes', 'CC'], + metric_weight=[1, 1, 1], + radius_or_number_of_bins=[32, 32, 4], + sampling_strategy=['Regular', 'Regular', None], + sampling_percentage=[0.25, 0.25, None], + sigma_units=['vox', 'vox', 'vox'], + number_of_iterations=[ + [1000, 500, 250, 0], + [1000, 500, 250, 0], + [100, 70, 50, 10], + ], + shrink_factors=[ + [8, 4, 2, 1], + [8, 4, 2, 1], + [8, 4, 2, 1], + ], + smoothing_sigmas=[ + [3, 2, 1, 0], + [3, 2, 1, 0], + [3, 2, 1, 0], + ], use_histogram_matching=True, write_composite_transform=True, ), From d878c9c1001bf2d0ea2f8a6ed8e06427587bae52 Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Tue, 9 Sep 2025 19:16:00 +0200 Subject: [PATCH 13/26] FIX: Implement atlas configuration loading and connections --- petprep/workflows/base.py | 34 +++++++++++++++++++++++++++- petprep/workflows/tests/test_base.py | 11 +++++++-- 2 files changed, 42 insertions(+), 3 deletions(-) diff --git a/petprep/workflows/base.py b/petprep/workflows/base.py index 1e74b79d..4ce0dd89 100644 --- a/petprep/workflows/base.py +++ b/petprep/workflows/base.py @@ -32,6 +32,7 @@ import os import sys import warnings +import json from copy import deepcopy from nipype.interfaces import utility as niu @@ -531,9 +532,10 @@ def init_single_subject_wf(subject_id: str): except Exception: # pragma: no cover - Py<3.9 fallback from importlib_resources import files as ir_files + atlas_config = ir_files('petprep.data.atlas') / 'config.json' seg_wf = init_atlas_wf( atlas=config.workflow.atlas, - config_file=str(ir_files('petprep.data.atlas') / 'config.json'), + config_file=str(atlas_config), name=f'pet_{config.workflow.atlas}_atlas_wf', ) workflow.connect( @@ -541,6 +543,36 @@ def init_single_subject_wf(subject_id: str): (anat_fit_wf, seg_wf, [('outputnode.t1w_preproc', 'inputnode.t1w_preproc')]), ] ) + + atlas_tpl = ( + json.loads(atlas_config.read_text()) + .get(config.workflow.atlas, {}) + .get('atlas', {}) + .get('tpl') + ) + if atlas_tpl and atlas_tpl in spaces.get_spaces(): + select_atlas_tpl_xfm = pe.Node( + KeySelect(fields=['std2anat_xfm'], key=atlas_tpl), + name='select_atlas_tpl_xfm', + run_without_submitting=True, + ) + workflow.connect( + [ + ( + anat_fit_wf, + select_atlas_tpl_xfm, + [ + ('outputnode.std2anat_xfm', 'std2anat_xfm'), + ('outputnode.template', 'keys'), + ], + ), + ( + select_atlas_tpl_xfm, + seg_wf, + [('std2anat_xfm', 'inputnode.tpl2anat_xfm')], + ), + ] + ) else: seg_wf = init_segmentation_wf( seg=config.workflow.seg, diff --git a/petprep/workflows/tests/test_base.py b/petprep/workflows/tests/test_base.py index 9eb8b218..3b44e304 100644 --- a/petprep/workflows/tests/test_base.py +++ b/petprep/workflows/tests/test_base.py @@ -150,7 +150,10 @@ def _dummy_atlas_wf(atlas, config_file, name='pet_atlas_wf'): from niworkflows.engine.workflows import LiterateWorkflow as Workflow wf = Workflow(name=name) - inputnode = pe.Node(niu.IdentityInterface(fields=['t1w_preproc']), name='inputnode') + inputnode = pe.Node( + niu.IdentityInterface(fields=['t1w_preproc', 'tpl2anat_xfm']), + name='inputnode', + ) outputnode = pe.Node( niu.IdentityInterface(fields=['segmentation', 'dseg_tsv']), name='outputnode', @@ -161,7 +164,7 @@ def _dummy_atlas_wf(atlas, config_file, name='pet_atlas_wf'): monkeypatch.setattr('petprep.workflows.pet.init_atlas_wf', _dummy_atlas_wf) with mock_config(bids_dir=multisession_bids_root): - config.workflow.atlas = 'foo' + config.workflow.atlas = 'DKT31' wf = init_single_subject_wf('01') flatgraph = wf._create_flat_graph() @@ -186,6 +189,10 @@ def _dummy_atlas_wf(atlas, config_file, name='pet_atlas_wf'): assert ('outputnode.dseg_tsv', 'inputnode.dseg_tsv') in edge['connect'] assert all('_atlas_wf' not in n for n in pet_node.list_node_names()) + select_node = wf.get_node('select_atlas_tpl_xfm') + edge = wf._graph.get_edge_data(select_node, atlas_node) + assert ('std2anat_xfm', 'inputnode.tpl2anat_xfm') in edge['connect'] + pet_node = wf.get_node(pet_wf_names[0]) assert config.workflow.atlas in pet_node.__desc__ assert config.workflow.seg not in pet_node.__desc__ From 549c99e9f4d94ffb61d400c8f2fa2795158c1364 Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Tue, 9 Sep 2025 19:19:09 +0200 Subject: [PATCH 14/26] FIX: Add optional input to workflow for transforms --- petprep/workflows/pet/atlas.py | 135 +++++++++++++--------- petprep/workflows/pet/tests/test_atlas.py | 40 +++++++ 2 files changed, 122 insertions(+), 53 deletions(-) diff --git a/petprep/workflows/pet/atlas.py b/petprep/workflows/pet/atlas.py index be801af9..b17fdb08 100644 --- a/petprep/workflows/pet/atlas.py +++ b/petprep/workflows/pet/atlas.py @@ -59,7 +59,12 @@ def _atlas_morph_tsv(segmentation: str, labels_tsv: str) -> str: return str(out_file) -def init_atlas_wf(atlas: str, config_file: str, name: str = "pet_atlas_wf") -> Workflow: +def init_atlas_wf( + atlas: str, + config_file: str, + tpl2anat_xfm: str | None = None, + name: str = "pet_atlas_wf", +) -> Workflow: """Map a template atlas into T1w space and compute regional volumes. Parameters @@ -68,6 +73,10 @@ def init_atlas_wf(atlas: str, config_file: str, name: str = "pet_atlas_wf") -> W Name of atlas (used in outputs and TemplateFlow queries). config_file : :class:`str` JSON file with query parameters for TemplateFlow ``get`` calls. + tpl2anat_xfm : :class:`str`, optional + Precomputed transform from template space to T1w space. When + provided, template-to-anatomical registration is skipped and this + transform is used directly to map the atlas into T1w space. name : :class:`str` Workflow name (default: ``pet_atlas_wf``). @@ -75,6 +84,9 @@ def init_atlas_wf(atlas: str, config_file: str, name: str = "pet_atlas_wf") -> W ------ t1w_preproc Preprocessed T1-weighted image. + tpl2anat_xfm + (Optional) Transform from template space to T1w space. If defined, + it will be used instead of estimating a new registration. Outputs ------- @@ -114,7 +126,12 @@ def _tf_kwargs(d: dict) -> dict: workflow = Workflow(name=name) - inputnode = pe.Node(niu.IdentityInterface(fields=["t1w_preproc"]), name="inputnode") + inputnode = pe.Node( + niu.IdentityInterface(fields=["t1w_preproc", "tpl2anat_xfm"]), + name="inputnode", + ) + if tpl2anat_xfm: + inputnode.inputs.tpl2anat_xfm = tpl2anat_xfm outputnode = pe.Node( niu.IdentityInterface(fields=["segmentation", "dseg_tsv"]), name="outputnode", @@ -123,37 +140,41 @@ def _tf_kwargs(d: dict) -> dict: label_source = pe.Node(niu.IdentityInterface(fields=["dseg_tsv"]), name="label_source") label_source.inputs.dseg_tsv = labels_tsv - reg = pe.Node( - Registration( - transforms=['Rigid', 'Affine', 'SyN'], - transform_parameters=[(0.1,), (0.1,), (0.1, 3, 0)], - metric=['Mattes', 'Mattes', 'CC'], - metric_weight=[1, 1, 1], - radius_or_number_of_bins=[32, 32, 4], - sampling_strategy=['Regular', 'Regular', None], - sampling_percentage=[0.25, 0.25, None], - sigma_units=['vox', 'vox', 'vox'], - number_of_iterations=[ - [1000, 500, 250, 0], - [1000, 500, 250, 0], - [100, 70, 50, 10], - ], - shrink_factors=[ - [8, 4, 2, 1], - [8, 4, 2, 1], - [8, 4, 2, 1], - ], - smoothing_sigmas=[ - [3, 2, 1, 0], - [3, 2, 1, 0], - [3, 2, 1, 0], - ], - use_histogram_matching=True, - write_composite_transform=True, - ), - name="t1_to_tpl", - ) - reg.inputs.fixed_image = template_t1w + # Decide whether to compute template-to-anatomical registration or reuse a provided transform + if not tpl2anat_xfm: + reg = pe.Node( + Registration( + transforms=['Rigid', 'Affine', 'SyN'], + transform_parameters=[(0.1,), (0.1,), (0.1, 3, 0)], + metric=['Mattes', 'Mattes', 'CC'], + metric_weight=[1, 1, 1], + radius_or_number_of_bins=[32, 32, 4], + sampling_strategy=['Regular', 'Regular', None], + sampling_percentage=[0.25, 0.25, None], + sigma_units=['vox', 'vox', 'vox'], + number_of_iterations=[ + [1000, 500, 250, 0], + [1000, 500, 250, 0], + [100, 70, 50, 10], + ], + shrink_factors=[ + [8, 4, 2, 1], + [8, 4, 2, 1], + [8, 4, 2, 1], + ], + smoothing_sigmas=[ + [3, 2, 1, 0], + [3, 2, 1, 0], + [3, 2, 1, 0], + ], + use_histogram_matching=True, + write_composite_transform=True, + ), + name="t1_to_tpl", + ) + reg.inputs.fixed_image = template_t1w + else: + reg = None apply_inv = pe.Node(ApplyTransforms(interpolation="NearestNeighbor"), name="apply_atlas") apply_inv.inputs.input_image = atlas_img @@ -221,25 +242,33 @@ def _tf_kwargs(d: dict) -> dict: mem_gb=DEFAULT_MEMORY_MIN_GB, ) - workflow.connect( - [ - (inputnode, reg, [("t1w_preproc", "moving_image")]), - (reg, apply_inv, [("inverse_composite_transform", "transforms")]), - (inputnode, apply_inv, [("t1w_preproc", "reference_image")]), - (apply_inv, outputnode, [("output_image", "segmentation")]), - (label_source, outputnode, [("dseg_tsv", "dseg_tsv")]), - (inputnode, sources, [("t1w_preproc", "in1")]), - (apply_inv, ds_seg, [("output_image", "in_file")]), - (inputnode, ds_seg, [("t1w_preproc", "source_file")]), - (sources, ds_seg, [("out", "Sources")]), - (label_source, ds_dseg_tsv, [("dseg_tsv", "in_file")]), - (inputnode, ds_dseg_tsv, [("t1w_preproc", "source_file")]), - (sources, ds_dseg_tsv, [("out", "Sources")]), - (apply_inv, gen_morph, [("output_image", "segmentation")]), - (gen_morph, ds_morph_tsv, [("out_file", "in_file")]), - (inputnode, ds_morph_tsv, [("t1w_preproc", "source_file")]), - (sources, ds_morph_tsv, [("out", "Sources")]), - ] - ) # fmt:skip + connections = [ + (inputnode, apply_inv, [("t1w_preproc", "reference_image")]), + (apply_inv, outputnode, [("output_image", "segmentation")]), + (label_source, outputnode, [("dseg_tsv", "dseg_tsv")]), + (inputnode, sources, [("t1w_preproc", "in1")]), + (apply_inv, ds_seg, [("output_image", "in_file")]), + (inputnode, ds_seg, [("t1w_preproc", "source_file")]), + (sources, ds_seg, [("out", "Sources")]), + (label_source, ds_dseg_tsv, [("dseg_tsv", "in_file")]), + (inputnode, ds_dseg_tsv, [("t1w_preproc", "source_file")]), + (sources, ds_dseg_tsv, [("out", "Sources")]), + (apply_inv, gen_morph, [("output_image", "segmentation")]), + (gen_morph, ds_morph_tsv, [("out_file", "in_file")]), + (inputnode, ds_morph_tsv, [("t1w_preproc", "source_file")]), + (sources, ds_morph_tsv, [("out", "Sources")]), + ] + + if reg is not None: + connections.extend( + [ + (inputnode, reg, [("t1w_preproc", "moving_image")]), + (reg, apply_inv, [("inverse_composite_transform", "transforms")]), + ] + ) + else: + connections.append((inputnode, apply_inv, [("tpl2anat_xfm", "transforms")])) + + workflow.connect(connections) # fmt:skip return workflow diff --git a/petprep/workflows/pet/tests/test_atlas.py b/petprep/workflows/pet/tests/test_atlas.py index 2506cc3b..f34e147e 100644 --- a/petprep/workflows/pet/tests/test_atlas.py +++ b/petprep/workflows/pet/tests/test_atlas.py @@ -46,6 +46,46 @@ def fake_get(**kwargs): assert wf.get_node('ds_morph_tsv').inputs.seg == 'MIAL67ThalamicNuclei' +def test_init_atlas_wf_with_xfm(tmp_path, monkeypatch): + t1_img = nb.Nifti1Image(np.zeros((2, 2, 2)), np.eye(4)) + t1_file = tmp_path / 't1w.nii.gz' + t1_img.to_filename(t1_file) + + atlas_img = nb.Nifti1Image(np.zeros((2, 2, 2)), np.eye(4)) + atlas_file = tmp_path / 'atlas.nii.gz' + atlas_img.to_filename(atlas_file) + + labels_file = tmp_path / 'labels.tsv' + labels_file.write_text('index\tname\n0\tbackground\n') + + cfg_file = ir_files('petprep.data.atlas') / 'config.json' + + def fake_get(**kwargs): + if kwargs.get('suffix') == 'T1w': + return str(t1_file) + if kwargs.get('extension') == 'tsv': + return str(labels_file) + return str(atlas_file) + + monkeypatch.setattr('petprep.workflows.pet.atlas.get_template', fake_get) + + config.execution.petprep_dir = tmp_path + config.execution.dataset_links = {} + + xfm = tmp_path / 'tpl2anat.txt' + xfm.write_text('0') + + wf = init_atlas_wf( + atlas='MIAL67ThalamicNuclei', + config_file=str(cfg_file), + tpl2anat_xfm=str(xfm), + ) + + node_names = [n.name for n in wf._get_all_nodes()] + assert 't1_to_tpl' not in node_names + assert wf.get_node('inputnode').inputs.tpl2anat_xfm == str(xfm) + + def test_init_atlas_wf_bad_name(tmp_path): cfg_file = ir_files('petprep.data.atlas') / 'config.json' with pytest.raises(ValueError, match="not found"): From c5f06cb51939ea8e59d9cc048661a6a692a7e5cf Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Tue, 9 Sep 2025 20:17:41 +0200 Subject: [PATCH 15/26] FIX: Load template from config.json --- petprep/workflows/base.py | 76 ++++++++++++++++++--------------------- 1 file changed, 34 insertions(+), 42 deletions(-) diff --git a/petprep/workflows/base.py b/petprep/workflows/base.py index 4ce0dd89..8c4897c1 100644 --- a/petprep/workflows/base.py +++ b/petprep/workflows/base.py @@ -239,11 +239,29 @@ def init_single_subject_wf(subject_id: str): spaces = config.workflow.spaces msm_sulc = config.workflow.run_msmsulc + atlas_config = None + atlas_tpl = None + if config.workflow.atlas: + try: # Py>=3.9 + from importlib.resources import files as ir_files + except Exception: # pragma: no cover - Py<3.9 fallback + from importlib_resources import files as ir_files + + atlas_config = ir_files('petprep.data.atlas') / 'config.json' + atlas_tpl = ( + json.loads(atlas_config.read_text()) + .get(config.workflow.atlas, {}) + .get('atlas', {}) + .get('tpl') + ) + anatomical_cache = {} if config.execution.derivatives: from smriprep.utils.bids import collect_derivatives as collect_anat_derivatives std_spaces = spaces.get_spaces(nonstandard=False, dim=(3,)) + if atlas_tpl and atlas_tpl not in std_spaces: + std_spaces.append(atlas_tpl) std_spaces.append('fsnative') for deriv_dir in config.execution.derivatives.values(): anatomical_cache.update( @@ -366,7 +384,19 @@ def init_single_subject_wf(subject_id: str): # Set up the template iterator once, if used template_iterator_wf = None - select_MNI2009c_xfm = None + select_atlas_tpl_xfm = None + if atlas_tpl and atlas_tpl in spaces.get_spaces(): + select_atlas_tpl_xfm = pe.Node( + KeySelect(fields=['std2anat_xfm'], key=atlas_tpl), + name='select_atlas_tpl_xfm', + run_without_submitting=True, + ) + workflow.connect([ + (anat_fit_wf, select_atlas_tpl_xfm, [ + ('outputnode.std2anat_xfm', 'std2anat_xfm'), + ('outputnode.template', 'keys'), + ]), + ]) # fmt:skip if config.workflow.level == 'full': if spaces.cached.get_spaces(nonstandard=False, dim=(3,)): template_iterator_wf = init_template_iterator_wf( @@ -398,19 +428,6 @@ def init_single_subject_wf(subject_id: str): ]), ]) # fmt:skip - if 'MNI152NLin2009cAsym' in spaces.get_spaces(): - select_MNI2009c_xfm = pe.Node( - KeySelect(fields=['std2anat_xfm'], key='MNI152NLin2009cAsym'), - name='select_MNI2009c_xfm', - run_without_submitting=True, - ) - workflow.connect([ - (anat_fit_wf, select_MNI2009c_xfm, [ - ('outputnode.std2anat_xfm', 'std2anat_xfm'), - ('outputnode.template', 'keys'), - ]), - ]) # fmt:skip - # Thread MNI152NLin6Asym standard outputs to CIFTI subworkflow, skipping # the iterator, which targets only output spaces. # This can lead to duplication in the working directory if people actually @@ -527,12 +544,6 @@ def init_single_subject_wf(subject_id: str): ]) # fmt:skip if config.workflow.atlas: - try: # Py>=3.9 - from importlib.resources import files as ir_files - except Exception: # pragma: no cover - Py<3.9 fallback - from importlib_resources import files as ir_files - - atlas_config = ir_files('petprep.data.atlas') / 'config.json' seg_wf = init_atlas_wf( atlas=config.workflow.atlas, config_file=str(atlas_config), @@ -544,28 +555,9 @@ def init_single_subject_wf(subject_id: str): ] ) - atlas_tpl = ( - json.loads(atlas_config.read_text()) - .get(config.workflow.atlas, {}) - .get('atlas', {}) - .get('tpl') - ) - if atlas_tpl and atlas_tpl in spaces.get_spaces(): - select_atlas_tpl_xfm = pe.Node( - KeySelect(fields=['std2anat_xfm'], key=atlas_tpl), - name='select_atlas_tpl_xfm', - run_without_submitting=True, - ) + if select_atlas_tpl_xfm is not None: workflow.connect( [ - ( - anat_fit_wf, - select_atlas_tpl_xfm, - [ - ('outputnode.std2anat_xfm', 'std2anat_xfm'), - ('outputnode.template', 'keys'), - ], - ), ( select_atlas_tpl_xfm, seg_wf, @@ -694,9 +686,9 @@ def init_single_subject_wf(subject_id: str): ]), ]) # fmt:skip - if select_MNI2009c_xfm is not None: + if select_atlas_tpl_xfm is not None: workflow.connect([ - (select_MNI2009c_xfm, pet_wf, [ + (select_atlas_tpl_xfm, pet_wf, [ ('std2anat_xfm', 'inputnode.mni2009c2anat_xfm'), ]), ]) # fmt:skip From dee44615e212c9a7988733a7ff8e140d441a5c42 Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Tue, 9 Sep 2025 20:26:42 +0200 Subject: [PATCH 16/26] FIX: Replace hard-coded warp lookup with dynamic name --- petprep/workflows/base.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/petprep/workflows/base.py b/petprep/workflows/base.py index 8c4897c1..b6b232db 100644 --- a/petprep/workflows/base.py +++ b/petprep/workflows/base.py @@ -272,6 +272,13 @@ def init_single_subject_wf(subject_id: str): ) ) + if atlas_tpl: + atlas_xfm = anatomical_cache.get('std2anat_xfm', {}).get(atlas_tpl) + if atlas_xfm is not None: + anatomical_cache['mni2009c2anat_xfm'] = atlas_xfm + else: + spaces.add(Reference(atlas_tpl, {})) + inputnode = pe.Node(niu.IdentityInterface(fields=['subjects_dir']), name='inputnode') bidssrc = pe.Node( From 20eb03b98e3da64515e7ab6ebca8e8b5ec3fba80 Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Tue, 9 Sep 2025 20:37:49 +0200 Subject: [PATCH 17/26] FIX: Remove conditional around select_atlas_tpl_xfm --- petprep/workflows/base.py | 46 ++++++++++++++------------------------- 1 file changed, 16 insertions(+), 30 deletions(-) diff --git a/petprep/workflows/base.py b/petprep/workflows/base.py index b6b232db..f208cf04 100644 --- a/petprep/workflows/base.py +++ b/petprep/workflows/base.py @@ -391,19 +391,17 @@ def init_single_subject_wf(subject_id: str): # Set up the template iterator once, if used template_iterator_wf = None - select_atlas_tpl_xfm = None - if atlas_tpl and atlas_tpl in spaces.get_spaces(): - select_atlas_tpl_xfm = pe.Node( - KeySelect(fields=['std2anat_xfm'], key=atlas_tpl), - name='select_atlas_tpl_xfm', - run_without_submitting=True, - ) - workflow.connect([ - (anat_fit_wf, select_atlas_tpl_xfm, [ - ('outputnode.std2anat_xfm', 'std2anat_xfm'), - ('outputnode.template', 'keys'), - ]), - ]) # fmt:skip + select_atlas_tpl_xfm = pe.Node( + KeySelect(fields=['std2anat_xfm'], key=atlas_tpl), + name='select_atlas_tpl_xfm', + run_without_submitting=True, + ) + workflow.connect([ + (anat_fit_wf, select_atlas_tpl_xfm, [ + ('outputnode.std2anat_xfm', 'std2anat_xfm'), + ('outputnode.template', 'keys'), + ]), + ]) # fmt:skip if config.workflow.level == 'full': if spaces.cached.get_spaces(nonstandard=False, dim=(3,)): template_iterator_wf = init_template_iterator_wf( @@ -556,22 +554,10 @@ def init_single_subject_wf(subject_id: str): config_file=str(atlas_config), name=f'pet_{config.workflow.atlas}_atlas_wf', ) - workflow.connect( - [ - (anat_fit_wf, seg_wf, [('outputnode.t1w_preproc', 'inputnode.t1w_preproc')]), - ] - ) - - if select_atlas_tpl_xfm is not None: - workflow.connect( - [ - ( - select_atlas_tpl_xfm, - seg_wf, - [('std2anat_xfm', 'inputnode.tpl2anat_xfm')], - ), - ] - ) + workflow.connect([ + (anat_fit_wf, seg_wf, [('outputnode.t1w_preproc', 'inputnode.t1w_preproc')]), + (select_atlas_tpl_xfm, seg_wf, [('std2anat_xfm', 'inputnode.tpl2anat_xfm')]), + ]) else: seg_wf = init_segmentation_wf( seg=config.workflow.seg, @@ -693,7 +679,7 @@ def init_single_subject_wf(subject_id: str): ]), ]) # fmt:skip - if select_atlas_tpl_xfm is not None: + if config.workflow.atlas: workflow.connect([ (select_atlas_tpl_xfm, pet_wf, [ ('std2anat_xfm', 'inputnode.mni2009c2anat_xfm'), From 8f6fbb6a3c2bd1cfe1a2ca4c8b691cc500e92e21 Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Tue, 9 Sep 2025 21:03:34 +0200 Subject: [PATCH 18/26] FIX: Update init_atlas_wf workflow inputs and connections --- petprep/workflows/base.py | 25 +++++++++++++---------- petprep/workflows/pet/atlas.py | 20 ++++++++++-------- petprep/workflows/pet/tests/test_atlas.py | 3 ++- petprep/workflows/tests/test_base.py | 2 +- 4 files changed, 29 insertions(+), 21 deletions(-) diff --git a/petprep/workflows/base.py b/petprep/workflows/base.py index f208cf04..37877852 100644 --- a/petprep/workflows/base.py +++ b/petprep/workflows/base.py @@ -391,17 +391,19 @@ def init_single_subject_wf(subject_id: str): # Set up the template iterator once, if used template_iterator_wf = None - select_atlas_tpl_xfm = pe.Node( - KeySelect(fields=['std2anat_xfm'], key=atlas_tpl), - name='select_atlas_tpl_xfm', - run_without_submitting=True, - ) - workflow.connect([ - (anat_fit_wf, select_atlas_tpl_xfm, [ - ('outputnode.std2anat_xfm', 'std2anat_xfm'), - ('outputnode.template', 'keys'), - ]), - ]) # fmt:skip + select_atlas_tpl_xfm = None + if atlas_tpl: + select_atlas_tpl_xfm = pe.Node( + KeySelect(fields=['std2anat_xfm'], key=atlas_tpl), + name='select_atlas_tpl_xfm', + run_without_submitting=True, + ) + workflow.connect([ + (anat_fit_wf, select_atlas_tpl_xfm, [ + ('outputnode.std2anat_xfm', 'std2anat_xfm'), + ('outputnode.template', 'keys'), + ]), + ]) # fmt:skip if config.workflow.level == 'full': if spaces.cached.get_spaces(nonstandard=False, dim=(3,)): template_iterator_wf = init_template_iterator_wf( @@ -552,6 +554,7 @@ def init_single_subject_wf(subject_id: str): seg_wf = init_atlas_wf( atlas=config.workflow.atlas, config_file=str(atlas_config), + tpl2anat_xfm=None, name=f'pet_{config.workflow.atlas}_atlas_wf', ) workflow.connect([ diff --git a/petprep/workflows/pet/atlas.py b/petprep/workflows/pet/atlas.py index b17fdb08..b436d0ac 100644 --- a/petprep/workflows/pet/atlas.py +++ b/petprep/workflows/pet/atlas.py @@ -62,21 +62,25 @@ def _atlas_morph_tsv(segmentation: str, labels_tsv: str) -> str: def init_atlas_wf( atlas: str, config_file: str, - tpl2anat_xfm: str | None = None, + tpl2anat_xfm: str | None, name: str = "pet_atlas_wf", ) -> Workflow: """Map a template atlas into T1w space and compute regional volumes. + This workflow expects a transform from template space to anatomical + space, typically generated by the anatomical workflow. If no transform + is provided, a template-to-anatomical registration will be estimated as + a fallback. + Parameters ---------- atlas : :class:`str` Name of atlas (used in outputs and TemplateFlow queries). config_file : :class:`str` JSON file with query parameters for TemplateFlow ``get`` calls. - tpl2anat_xfm : :class:`str`, optional - Precomputed transform from template space to T1w space. When - provided, template-to-anatomical registration is skipped and this - transform is used directly to map the atlas into T1w space. + tpl2anat_xfm : :class:`str` + Transform from template space to T1w space. If ``None``, the workflow + will estimate this registration internally. name : :class:`str` Workflow name (default: ``pet_atlas_wf``). @@ -85,8 +89,8 @@ def init_atlas_wf( t1w_preproc Preprocessed T1-weighted image. tpl2anat_xfm - (Optional) Transform from template space to T1w space. If defined, - it will be used instead of estimating a new registration. + Transform from template space to T1w space, typically produced by + the anatomical workflow. Outputs ------- @@ -140,7 +144,7 @@ def _tf_kwargs(d: dict) -> dict: label_source = pe.Node(niu.IdentityInterface(fields=["dseg_tsv"]), name="label_source") label_source.inputs.dseg_tsv = labels_tsv - # Decide whether to compute template-to-anatomical registration or reuse a provided transform + # Atlas registration is handled upstream; estimate a transform only if none is provided if not tpl2anat_xfm: reg = pe.Node( Registration( diff --git a/petprep/workflows/pet/tests/test_atlas.py b/petprep/workflows/pet/tests/test_atlas.py index f34e147e..adeefb40 100644 --- a/petprep/workflows/pet/tests/test_atlas.py +++ b/petprep/workflows/pet/tests/test_atlas.py @@ -37,6 +37,7 @@ def fake_get(**kwargs): wf = init_atlas_wf( atlas='MIAL67ThalamicNuclei', config_file=str(cfg_file), + tpl2anat_xfm=None, ) node_names = [n.name for n in wf._get_all_nodes()] assert 'apply_atlas' in node_names @@ -89,7 +90,7 @@ def fake_get(**kwargs): def test_init_atlas_wf_bad_name(tmp_path): cfg_file = ir_files('petprep.data.atlas') / 'config.json' with pytest.raises(ValueError, match="not found"): - init_atlas_wf(atlas='notreal', config_file=str(cfg_file)) + init_atlas_wf(atlas='notreal', config_file=str(cfg_file), tpl2anat_xfm=None) def test_atlas_morph_tsv(tmp_path): diff --git a/petprep/workflows/tests/test_base.py b/petprep/workflows/tests/test_base.py index 3b44e304..995c5cac 100644 --- a/petprep/workflows/tests/test_base.py +++ b/petprep/workflows/tests/test_base.py @@ -144,7 +144,7 @@ def test_segmentation_shared_across_runs(multisession_bids_root): def test_atlas_replaces_segmentation(monkeypatch, multisession_bids_root): - def _dummy_atlas_wf(atlas, config_file, name='pet_atlas_wf'): + def _dummy_atlas_wf(atlas, config_file, tpl2anat_xfm=None, name='pet_atlas_wf'): from nipype.interfaces import utility as niu from nipype.pipeline import engine as pe from niworkflows.engine.workflows import LiterateWorkflow as Workflow From f3793a787d0ec47ba420ca3c2731522195f9bdee Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Tue, 9 Sep 2025 21:24:49 +0200 Subject: [PATCH 19/26] FIX: Update init_atlas_wf argument handling --- petprep/workflows/base.py | 13 +++++--- petprep/workflows/tests/test_base.py | 48 ++++++++++++++++++++++++++++ 2 files changed, 57 insertions(+), 4 deletions(-) diff --git a/petprep/workflows/base.py b/petprep/workflows/base.py index 37877852..7c8fa51a 100644 --- a/petprep/workflows/base.py +++ b/petprep/workflows/base.py @@ -241,6 +241,7 @@ def init_single_subject_wf(subject_id: str): atlas_config = None atlas_tpl = None + atlas_xfm = None if config.workflow.atlas: try: # Py>=3.9 from importlib.resources import files as ir_files @@ -554,13 +555,17 @@ def init_single_subject_wf(subject_id: str): seg_wf = init_atlas_wf( atlas=config.workflow.atlas, config_file=str(atlas_config), - tpl2anat_xfm=None, + tpl2anat_xfm=atlas_xfm, name=f'pet_{config.workflow.atlas}_atlas_wf', ) - workflow.connect([ + seg_connect = [ (anat_fit_wf, seg_wf, [('outputnode.t1w_preproc', 'inputnode.t1w_preproc')]), - (select_atlas_tpl_xfm, seg_wf, [('std2anat_xfm', 'inputnode.tpl2anat_xfm')]), - ]) + ] + if atlas_xfm is None and select_atlas_tpl_xfm is not None: + seg_connect.append( + (select_atlas_tpl_xfm, seg_wf, [('std2anat_xfm', 'inputnode.tpl2anat_xfm')]) + ) + workflow.connect(seg_connect) else: seg_wf = init_segmentation_wf( seg=config.workflow.seg, diff --git a/petprep/workflows/tests/test_base.py b/petprep/workflows/tests/test_base.py index 995c5cac..5bcc2169 100644 --- a/petprep/workflows/tests/test_base.py +++ b/petprep/workflows/tests/test_base.py @@ -198,6 +198,54 @@ def _dummy_atlas_wf(atlas, config_file, tpl2anat_xfm=None, name='pet_atlas_wf'): assert config.workflow.seg not in pet_node.__desc__ +def test_atlas_uses_precomputed_xfm(monkeypatch, multisession_bids_root, tmp_path): + """init_atlas_wf should be initialized with a cached transform.""" + + seen = {} + + def _dummy_atlas_wf(atlas, config_file, tpl2anat_xfm=None, name='pet_atlas_wf'): + from nipype.interfaces import utility as niu + from nipype.pipeline import engine as pe + from niworkflows.engine.workflows import LiterateWorkflow as Workflow + + seen['tpl2anat_xfm'] = tpl2anat_xfm + + wf = Workflow(name=name) + inputnode = pe.Node( + niu.IdentityInterface(fields=['t1w_preproc', 'tpl2anat_xfm']), + name='inputnode', + ) + outputnode = pe.Node( + niu.IdentityInterface(fields=['segmentation', 'dseg_tsv']), + name='outputnode', + ) + wf.add_nodes([inputnode, outputnode]) + return wf + + monkeypatch.setattr('petprep.workflows.pet.init_atlas_wf', _dummy_atlas_wf) + + atlas_tpl = 'MNI152NLin2009cAsym' + xfm_file = tmp_path / 'tpl2anat_xfm.txt' + xfm_file.write_text('0') + + def fake_collect_derivatives(derivatives_dir, subject_id, std_spaces): + return {'std2anat_xfm': {atlas_tpl: str(xfm_file)}} + + monkeypatch.setattr('smriprep.utils.bids.collect_derivatives', fake_collect_derivatives) + + with mock_config(bids_dir=multisession_bids_root): + config.workflow.atlas = 'DKT31' + config.execution.derivatives = {'smriprep': tmp_path} + wf = init_single_subject_wf('01') + + assert seen['tpl2anat_xfm'] == str(xfm_file) + + atlas_wf_name = f'pet_{config.workflow.atlas}_atlas_wf' + atlas_node = wf.get_node(atlas_wf_name) + select_node = wf.get_node('select_atlas_tpl_xfm') + assert not wf._graph.has_edge(select_node, atlas_node) + + def _make_params( pet2anat_init: str = 'auto', medial_surface_nan: bool = False, From 0cc45abcb4c5fc7f316ecef9aeb1aed0a7fa5795 Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Wed, 10 Sep 2025 09:30:25 +0200 Subject: [PATCH 20/26] FIX: Make registration conditional at runtime --- petprep/workflows/pet/atlas.py | 86 +++++++++++------------ petprep/workflows/pet/tests/test_atlas.py | 16 +++-- 2 files changed, 55 insertions(+), 47 deletions(-) diff --git a/petprep/workflows/pet/atlas.py b/petprep/workflows/pet/atlas.py index b436d0ac..c3f3eec2 100644 --- a/petprep/workflows/pet/atlas.py +++ b/petprep/workflows/pet/atlas.py @@ -144,41 +144,49 @@ def _tf_kwargs(d: dict) -> dict: label_source = pe.Node(niu.IdentityInterface(fields=["dseg_tsv"]), name="label_source") label_source.inputs.dseg_tsv = labels_tsv - # Atlas registration is handled upstream; estimate a transform only if none is provided - if not tpl2anat_xfm: - reg = pe.Node( - Registration( - transforms=['Rigid', 'Affine', 'SyN'], - transform_parameters=[(0.1,), (0.1,), (0.1, 3, 0)], - metric=['Mattes', 'Mattes', 'CC'], - metric_weight=[1, 1, 1], - radius_or_number_of_bins=[32, 32, 4], - sampling_strategy=['Regular', 'Regular', None], - sampling_percentage=[0.25, 0.25, None], - sigma_units=['vox', 'vox', 'vox'], - number_of_iterations=[ - [1000, 500, 250, 0], - [1000, 500, 250, 0], - [100, 70, 50, 10], - ], - shrink_factors=[ - [8, 4, 2, 1], - [8, 4, 2, 1], - [8, 4, 2, 1], - ], - smoothing_sigmas=[ - [3, 2, 1, 0], - [3, 2, 1, 0], - [3, 2, 1, 0], - ], - use_histogram_matching=True, - write_composite_transform=True, - ), - name="t1_to_tpl", + def _get_xfm(tpl2anat_xfm: str | None, t1w_preproc: str, template_t1w: str) -> str: + if tpl2anat_xfm: + return tpl2anat_xfm + reg = Registration( + transforms=['Rigid', 'Affine', 'SyN'], + transform_parameters=[(0.1,), (0.1,), (0.1, 3, 0)], + metric=['Mattes', 'Mattes', 'CC'], + metric_weight=[1, 1, 1], + radius_or_number_of_bins=[32, 32, 4], + sampling_strategy=['Regular', 'Regular', None], + sampling_percentage=[0.25, 0.25, None], + sigma_units=['vox', 'vox', 'vox'], + number_of_iterations=[ + [1000, 500, 250, 0], + [1000, 500, 250, 0], + [100, 70, 50, 10], + ], + shrink_factors=[ + [8, 4, 2, 1], + [8, 4, 2, 1], + [8, 4, 2, 1], + ], + smoothing_sigmas=[ + [3, 2, 1, 0], + [3, 2, 1, 0], + [3, 2, 1, 0], + ], + use_histogram_matching=True, + write_composite_transform=True, ) reg.inputs.fixed_image = template_t1w - else: - reg = None + reg.inputs.moving_image = t1w_preproc + return reg.run().outputs.inverse_composite_transform + + xfm_source = pe.Node( + niu.Function( + input_names=["tpl2anat_xfm", "t1w_preproc", "template_t1w"], + output_names=["tpl2anat_xfm"], + function=_get_xfm, + ), + name="t1_to_tpl", + ) + xfm_source.inputs.template_t1w = template_t1w apply_inv = pe.Node(ApplyTransforms(interpolation="NearestNeighbor"), name="apply_atlas") apply_inv.inputs.input_image = atlas_img @@ -248,6 +256,8 @@ def _tf_kwargs(d: dict) -> dict: connections = [ (inputnode, apply_inv, [("t1w_preproc", "reference_image")]), + (inputnode, xfm_source, [("tpl2anat_xfm", "tpl2anat_xfm"), ("t1w_preproc", "t1w_preproc")]), + (xfm_source, apply_inv, [("tpl2anat_xfm", "transforms")]), (apply_inv, outputnode, [("output_image", "segmentation")]), (label_source, outputnode, [("dseg_tsv", "dseg_tsv")]), (inputnode, sources, [("t1w_preproc", "in1")]), @@ -263,16 +273,6 @@ def _tf_kwargs(d: dict) -> dict: (sources, ds_morph_tsv, [("out", "Sources")]), ] - if reg is not None: - connections.extend( - [ - (inputnode, reg, [("t1w_preproc", "moving_image")]), - (reg, apply_inv, [("inverse_composite_transform", "transforms")]), - ] - ) - else: - connections.append((inputnode, apply_inv, [("tpl2anat_xfm", "transforms")])) - workflow.connect(connections) # fmt:skip return workflow diff --git a/petprep/workflows/pet/tests/test_atlas.py b/petprep/workflows/pet/tests/test_atlas.py index adeefb40..64737ad6 100644 --- a/petprep/workflows/pet/tests/test_atlas.py +++ b/petprep/workflows/pet/tests/test_atlas.py @@ -41,6 +41,7 @@ def fake_get(**kwargs): ) node_names = [n.name for n in wf._get_all_nodes()] assert 'apply_atlas' in node_names + assert 't1_to_tpl' in node_names assert wf.get_node('label_source').inputs.dseg_tsv == str(labels_file) assert wf.get_node('ds_seg').inputs.seg == 'MIAL67ThalamicNuclei' assert wf.get_node('ds_dseg_tsv').inputs.seg == 'MIAL67ThalamicNuclei' @@ -79,12 +80,19 @@ def fake_get(**kwargs): wf = init_atlas_wf( atlas='MIAL67ThalamicNuclei', config_file=str(cfg_file), - tpl2anat_xfm=str(xfm), + tpl2anat_xfm=None, ) - node_names = [n.name for n in wf._get_all_nodes()] - assert 't1_to_tpl' not in node_names - assert wf.get_node('inputnode').inputs.tpl2anat_xfm == str(xfm) + t1_to_tpl = wf.get_node('t1_to_tpl') + t1_to_tpl.inputs.tpl2anat_xfm = str(xfm) + t1_to_tpl.inputs.t1w_preproc = str(t1_file) + + def _fail_run(self, *args, **kwargs): + raise AssertionError('Registration should not run') + + monkeypatch.setattr('petprep.workflows.pet.atlas.Registration.run', _fail_run) + result = t1_to_tpl.run() + assert result.outputs.tpl2anat_xfm == str(xfm) def test_init_atlas_wf_bad_name(tmp_path): From 185795574a02a238f6e3c16272cbd2830aa6eefb Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Wed, 10 Sep 2025 09:34:59 +0200 Subject: [PATCH 21/26] FIX: Update derivative collection and tests --- petprep/workflows/base.py | 28 ++++++++++++++++++++++------ petprep/workflows/tests/test_base.py | 2 +- 2 files changed, 23 insertions(+), 7 deletions(-) diff --git a/petprep/workflows/base.py b/petprep/workflows/base.py index 7c8fa51a..9a561d98 100644 --- a/petprep/workflows/base.py +++ b/petprep/workflows/base.py @@ -265,14 +265,30 @@ def init_single_subject_wf(subject_id: str): std_spaces.append(atlas_tpl) std_spaces.append('fsnative') for deriv_dir in config.execution.derivatives.values(): - anatomical_cache.update( - collect_anat_derivatives( - derivatives_dir=deriv_dir, - subject_id=subject_id, - std_spaces=std_spaces, - ) + derivs = collect_anat_derivatives( + derivatives_dir=deriv_dir, + subject_id=subject_id, + std_spaces=std_spaces, ) + transforms = derivs.pop('transforms', {}) + if transforms: + std2anat = anatomical_cache.setdefault('std2anat_xfm', {}) + anat2std = anatomical_cache.setdefault('anat2std_xfm', {}) + for space, xfm in transforms.items(): + fwd = xfm.get('forward') + rev = xfm.get('reverse') + if isinstance(fwd, list): + fwd = fwd[0] + if isinstance(rev, list): + rev = rev[0] + if fwd: + anat2std[space] = fwd + if rev: + std2anat[space] = rev + + anatomical_cache.update(derivs) + if atlas_tpl: atlas_xfm = anatomical_cache.get('std2anat_xfm', {}).get(atlas_tpl) if atlas_xfm is not None: diff --git a/petprep/workflows/tests/test_base.py b/petprep/workflows/tests/test_base.py index 5bcc2169..5831dd5e 100644 --- a/petprep/workflows/tests/test_base.py +++ b/petprep/workflows/tests/test_base.py @@ -229,7 +229,7 @@ def _dummy_atlas_wf(atlas, config_file, tpl2anat_xfm=None, name='pet_atlas_wf'): xfm_file.write_text('0') def fake_collect_derivatives(derivatives_dir, subject_id, std_spaces): - return {'std2anat_xfm': {atlas_tpl: str(xfm_file)}} + return {'transforms': {atlas_tpl: {'reverse': str(xfm_file)}}} monkeypatch.setattr('smriprep.utils.bids.collect_derivatives', fake_collect_derivatives) From 4d08ed46a07b7dda4023939475f23831042fce24 Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Wed, 10 Sep 2025 12:20:26 +0200 Subject: [PATCH 22/26] FIX: Refactor MNI152NLin2009cAsym transform retrieval --- petprep/workflows/base.py | 32 +++++++++++++++++++++++--------- 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/petprep/workflows/base.py b/petprep/workflows/base.py index 9a561d98..0a0e9b57 100644 --- a/petprep/workflows/base.py +++ b/petprep/workflows/base.py @@ -289,11 +289,15 @@ def init_single_subject_wf(subject_id: str): anatomical_cache.update(derivs) + mni2009c_xfm = anatomical_cache.get('std2anat_xfm', {}).get('MNI152NLin2009cAsym') + if mni2009c_xfm is not None: + anatomical_cache['mni2009c2anat_xfm'] = mni2009c_xfm + else: + spaces.add(Reference('MNI152NLin2009cAsym', {})) + if atlas_tpl: atlas_xfm = anatomical_cache.get('std2anat_xfm', {}).get(atlas_tpl) - if atlas_xfm is not None: - anatomical_cache['mni2009c2anat_xfm'] = atlas_xfm - else: + if atlas_xfm is None: spaces.add(Reference(atlas_tpl, {})) inputnode = pe.Node(niu.IdentityInterface(fields=['subjects_dir']), name='inputnode') @@ -408,6 +412,11 @@ def init_single_subject_wf(subject_id: str): # Set up the template iterator once, if used template_iterator_wf = None + select_MNI2009c_xfm = pe.Node( + KeySelect(fields=['std2anat_xfm'], key='MNI152NLin2009cAsym'), + name='select_MNI2009c_xfm', + run_without_submitting=True, + ) select_atlas_tpl_xfm = None if atlas_tpl: select_atlas_tpl_xfm = pe.Node( @@ -421,6 +430,12 @@ def init_single_subject_wf(subject_id: str): ('outputnode.template', 'keys'), ]), ]) # fmt:skip + workflow.connect([ + (anat_fit_wf, select_MNI2009c_xfm, [ + ('outputnode.std2anat_xfm', 'std2anat_xfm'), + ('outputnode.template', 'keys'), + ]), + ]) # fmt:skip if config.workflow.level == 'full': if spaces.cached.get_spaces(nonstandard=False, dim=(3,)): template_iterator_wf = init_template_iterator_wf( @@ -703,12 +718,11 @@ def init_single_subject_wf(subject_id: str): ]), ]) # fmt:skip - if config.workflow.atlas: - workflow.connect([ - (select_atlas_tpl_xfm, pet_wf, [ - ('std2anat_xfm', 'inputnode.mni2009c2anat_xfm'), - ]), - ]) # fmt:skip + workflow.connect([ + (select_MNI2009c_xfm, pet_wf, [ + ('std2anat_xfm', 'inputnode.mni2009c2anat_xfm'), + ]), + ]) # fmt:skip # Thread MNI152NLin6Asym standard outputs to CIFTI subworkflow, skipping # the iterator, which targets only output spaces. From b4f966507d60b6d64ab92b3071b3e5a3a961e6d4 Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Wed, 10 Sep 2025 12:31:40 +0200 Subject: [PATCH 23/26] FIX: Extend niworkflows with SUIT support --- petprep/__init__.py | 12 ++++++++++++ petprep/workflows/tests/test_base.py | 9 +++++++++ scripts/fetch_templates.py | 1 + 3 files changed, 22 insertions(+) diff --git a/petprep/__init__.py b/petprep/__init__.py index 2afc7fea..db8a67cf 100644 --- a/petprep/__init__.py +++ b/petprep/__init__.py @@ -26,3 +26,15 @@ from ._version import __version__ except ImportError: __version__ = '0+unknown' + +# Ensure niworkflows recognizes the SUIT template. +try: # pragma: no cover - if niworkflows is unavailable, skip patching + from niworkflows.utils import spaces as _spaces + + if 'SUIT' not in _spaces.Reference._standard_spaces: + _spaces.Reference._standard_spaces = ( + *_spaces.Reference._standard_spaces, + 'SUIT', + ) +except Exception: # pragma: no cover + pass diff --git a/petprep/workflows/tests/test_base.py b/petprep/workflows/tests/test_base.py index 5831dd5e..4d5b0a91 100644 --- a/petprep/workflows/tests/test_base.py +++ b/petprep/workflows/tests/test_base.py @@ -198,6 +198,15 @@ def _dummy_atlas_wf(atlas, config_file, tpl2anat_xfm=None, name='pet_atlas_wf'): assert config.workflow.seg not in pet_node.__desc__ +def test_accepts_suit_atlas(bids_root): + with mock_config(bids_dir=bids_root): + config.workflow.atlas = 'Buckner201117n' + try: + init_single_subject_wf('01') + except ValueError as e: # pragma: no cover - explicit failure capture + pytest.fail(f'init_single_subject_wf raised {e!r}') + + def test_atlas_uses_precomputed_xfm(monkeypatch, multisession_bids_root, tmp_path): """init_atlas_wf should be initialized with a cached transform.""" diff --git a/scripts/fetch_templates.py b/scripts/fetch_templates.py index 97f9af7c..bee3ff8c 100755 --- a/scripts/fetch_templates.py +++ b/scripts/fetch_templates.py @@ -154,6 +154,7 @@ def fetch_all(): fetch_OASIS() fetch_fsaverage() fetch_fsLR() + fetch_SUIT() if __name__ == '__main__': From a39aa5046933e75e71ca446912775b18e26eaa7a Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Wed, 10 Sep 2025 12:43:18 +0200 Subject: [PATCH 24/26] FIX: Add datalad to project dependencies --- pyproject.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index 1fafbe49..b58ca8f2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -40,6 +40,8 @@ dependencies = [ "transforms3d >= 0.4", "toml >= 0.10", "codecarbon >= 2", + "datalad", + "datalad-osf", "APScheduler >= 3.10", ] dynamic = ["version"] From e5b3c8986c3c64e2145819228a925336fe401cdd Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Wed, 10 Sep 2025 15:41:03 +0200 Subject: [PATCH 25/26] FIX: Add tpl2anat entry to collect derivatives --- petprep/data/io_spec.json | 8 +++++++ petprep/utils/bids.py | 7 ++++++ petprep/utils/tests/test_derivative_cache.py | 25 ++++++++++++++++++++ 3 files changed, 40 insertions(+) diff --git a/petprep/data/io_spec.json b/petprep/data/io_spec.json index 6ce6d4d6..4d67600d 100644 --- a/petprep/data/io_spec.json +++ b/petprep/data/io_spec.json @@ -38,6 +38,14 @@ "mode": "image", "suffix": "xfm", "extension": ".txt" + }, + "tpl2anat": { + "datatype": "anat", + "from": ["MNI152NLin2009cAsym", "MNI152NLin6Asym", "SUIT"], + "to": ["T1w", "anat"], + "mode": "image", + "suffix": "xfm", + "extension": [".h5", ".txt"] } } }, diff --git a/petprep/utils/bids.py b/petprep/utils/bids.py index 97bc8cff..41926a84 100644 --- a/petprep/utils/bids.py +++ b/petprep/utils/bids.py @@ -86,6 +86,13 @@ def collect_derivatives( item = layout.get(return_type='filename', **query) if not item: continue + if xfm == 'tpl2anat': + tpl_cache = transforms_cache.setdefault('tpl2anat', {}) + for fname in item: + tpl = layout.parse_file_entities(fname).get('from') + if tpl: + tpl_cache[tpl] = fname + continue transforms_cache[xfm] = item[0] if len(item) == 1 else item derivs_cache['transforms'] = transforms_cache return derivs_cache diff --git a/petprep/utils/tests/test_derivative_cache.py b/petprep/utils/tests/test_derivative_cache.py index aa43b10e..02185c1d 100644 --- a/petprep/utils/tests/test_derivative_cache.py +++ b/petprep/utils/tests/test_derivative_cache.py @@ -56,3 +56,28 @@ def test_transforms_found_as_str(tmp_path: Path, xfm: str): entities=entities, ) assert derivs == {'transforms': {xfm: str(to_find)}} + + +def test_tpl2anat_transform_found(tmp_path: Path): + subject = '01' + tpl = 'MNI152NLin2009cAsym' + to_find = tmp_path.joinpath( + f'sub-{subject}', + 'anat', + f'sub-{subject}_from-{tpl}_to-T1w_mode-image_xfm.h5', + ) + to_find.parent.mkdir(parents=True) + to_find.touch() + + entities = { + 'subject': subject, + 'suffix': 'pet', + 'extension': '.nii.gz', + } + + derivs = bids.collect_derivatives( + derivatives_dir=tmp_path, + entities=entities, + ) + + assert derivs == {'transforms': {'tpl2anat': {tpl: str(to_find)}}} From df591498e5c98e5bf22f8915d28e2ae5bd048d99 Mon Sep 17 00:00:00 2001 From: mnoergaard Date: Wed, 10 Sep 2025 15:58:45 +0200 Subject: [PATCH 26/26] FIX: Update _atlas_morph_tsv for column filtering --- petprep/workflows/pet/atlas.py | 6 +++--- petprep/workflows/pet/tests/test_atlas.py | 5 ++++- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/petprep/workflows/pet/atlas.py b/petprep/workflows/pet/atlas.py index c3f3eec2..a33d655e 100644 --- a/petprep/workflows/pet/atlas.py +++ b/petprep/workflows/pet/atlas.py @@ -50,10 +50,10 @@ def _atlas_morph_tsv(segmentation: str, labels_tsv: str) -> str: seg_img = nb.load(segmentation) data = np.asanyarray(seg_img.dataobj) voxvol = np.prod(seg_img.header.get_zooms()[:3]) - labels = pd.read_csv(labels_tsv, sep="\t") + labels = pd.read_csv(labels_tsv, sep="\t", usecols=["index", "name"]) volumes = [(data == int(idx)).sum() * voxvol for idx in labels["index"]] - out = labels.copy() - out["volume-mm3"] = volumes + out = labels.assign(**{"volume-mm3": volumes}) + out = out[["index", "name", "volume-mm3"]] out_file = Path("morph.tsv").absolute() out.to_csv(out_file, sep="\t", index=False) return str(out_file) diff --git a/petprep/workflows/pet/tests/test_atlas.py b/petprep/workflows/pet/tests/test_atlas.py index 64737ad6..47fa63af 100644 --- a/petprep/workflows/pet/tests/test_atlas.py +++ b/petprep/workflows/pet/tests/test_atlas.py @@ -111,9 +111,12 @@ def test_atlas_morph_tsv(tmp_path): seg.to_filename(seg_file) labels = tmp_path / 'labels.tsv' - labels.write_text('index\tname\n0\tbg\n1\tregion\n') + labels.write_text( + 'index\tname\tcolor\n0\tbg\t#000000\n1\tregion\t#ff0000\n' + ) out = _atlas_morph_tsv(str(seg_file), str(labels)) df = pd.read_csv(out, sep='\t') + assert list(df.columns) == ['index', 'name', 'volume-mm3'] volume = df.loc[df['index'] == 1, 'volume-mm3'].iloc[0] assert volume == 4