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/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 diff --git a/petprep/data/atlas/config.json b/petprep/data/atlas/config.json new file mode 100644 index 00000000..33f64651 --- /dev/null +++ b/petprep/data/atlas/config.json @@ -0,0 +1,726 @@ +{ + "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", + "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", + "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": { + "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", + "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" + } + }, + "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", + "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" + } + } +} 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)}}} diff --git a/petprep/workflows/base.py b/petprep/workflows/base.py index f71c6f7d..0a0e9b57 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 @@ -164,6 +165,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') @@ -237,21 +239,67 @@ def init_single_subject_wf(subject_id: str): spaces = config.workflow.spaces msm_sulc = config.workflow.run_msmsulc + 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 + 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( - 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) + + 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 None: + spaces.add(Reference(atlas_tpl, {})) + inputnode = pe.Node(niu.IdentityInterface(fields=['subjects_dir']), name='inputnode') bidssrc = pe.Node( @@ -364,7 +412,30 @@ 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_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( + 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 + 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( @@ -396,19 +467,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 @@ -524,29 +582,51 @@ 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: + seg_wf = init_atlas_wf( + atlas=config.workflow.atlas, + config_file=str(atlas_config), + tpl2anat_xfm=atlas_xfm, + name=f'pet_{config.workflow.atlas}_atlas_wf', + ) + seg_connect = [ + (anat_fit_wf, seg_wf, [('outputnode.t1w_preproc', 'inputnode.t1w_preproc')]), ] - ) + 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, + 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 +635,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 +699,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'), ]), @@ -638,12 +718,11 @@ def init_single_subject_wf(subject_id: str): ]), ]) # fmt:skip - if select_MNI2009c_xfm is not None: - workflow.connect([ - (select_MNI2009c_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. 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..a33d655e --- /dev/null +++ b/petprep/workflows/pet/atlas.py @@ -0,0 +1,278 @@ +# 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 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.""" + from pathlib import Path + 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", usecols=["index", "name"]) + volumes = [(data == int(idx)).sum() * voxvol for idx in labels["index"]] + 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) + + +def init_atlas_wf( + atlas: str, + config_file: str, + 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` + 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``). + + Inputs + ------ + t1w_preproc + Preprocessed T1-weighted image. + tpl2anat_xfm + Transform from template space to T1w space, typically produced by + the anatomical workflow. + + Outputs + ------- + segmentation + Atlas segmentation in T1w space. + dseg_tsv + Label table for the atlas. + """ + with open(config_file) as f: + 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) + + 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", + ) + + label_source = pe.Node(niu.IdentityInterface(fields=["dseg_tsv"]), name="label_source") + label_source.inputs.dseg_tsv = labels_tsv + + 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 + 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 + + 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, + ) + + 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")]), + (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")]), + ] + + workflow.connect(connections) # fmt:skip + + return workflow 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/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_atlas.py b/petprep/workflows/pet/tests/test_atlas.py new file mode 100644 index 00000000..47fa63af --- /dev/null +++ b/petprep/workflows/pet/tests/test_atlas.py @@ -0,0 +1,122 @@ +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 + + +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_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 = {} + + 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 + 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' + 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=None, + ) + + 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): + 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), tpl2anat_xfm=None) + + +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\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 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) 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')] diff --git a/petprep/workflows/tests/test_base.py b/petprep/workflows/tests/test_base.py index a1eeca35..4d5b0a91 100644 --- a/petprep/workflows/tests/test_base.py +++ b/petprep/workflows/tests/test_base.py @@ -143,6 +143,118 @@ 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, 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 + + 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) + + with mock_config(bids_dir=multisession_bids_root): + config.workflow.atlas = 'DKT31' + 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()) + + 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__ + + +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.""" + + 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 {'transforms': {atlas_tpl: {'reverse': 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, 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"] diff --git a/scripts/fetch_templates.py b/scripts/fetch_templates.py index a7bec82e..bee3ff8c 100755 --- a/scripts/fetch_templates.py +++ b/scripts/fetch_templates.py @@ -109,12 +109,52 @@ 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() fetch_OASIS() fetch_fsaverage() fetch_fsLR() + fetch_SUIT() if __name__ == '__main__':