Merge branch 'add_refmask_petprep' into main

Author: mnoergaard

docs/usage.rst

@@ -203,6 +203,13 @@ method with a 5 mm PSF::
     $ petprep /data/bids_root /out participant \
         --pvc-tool petpvc --pvc-method GTM --pvc-psf 5
 
+Reference region masks
+----------------------
+*PETPrep* can build masks for reference regions used in quantification.
+Use ``--ref-mask-name`` to select a predefined region and
+``--ref-mask-index`` to override the label indices.
+The presets are defined in ``petprep/data/reference_mask/config.json``.
+
 Troubleshooting
 ---------------
 Logs and crashfiles are output into the

petprep/cli/parser.py

@@ -554,6 +554,18 @@ def _bids_filter(value, parser):
         help='Segmentation method to use.',
     )
 
+    g_refmask = parser.add_argument_group('Options for reference mask generation')
+    g_refmask.add_argument(
+        '--ref-mask-name',
+        help='Predefined reference regions. Pair with --ref-mask-index to define custom labels.',
+    )
+    g_refmask.add_argument(
+        '--ref-mask-index',
+        nargs='+',
+        type=int,
+        help='List of segmentation indices to use for the reference mask.',
+    )
+
     g_pvc = parser.add_argument_group('Options for partial volume correction')
 
     try:
@@ -721,6 +733,14 @@ def parse_args(args=None, namespace=None):
             'Options --pvc-tool, --pvc-method and --pvc-psf must be used together.'
         )
 
+    if opts.ref_mask_index is not None and opts.ref_mask_name is None:
+        parser.error('Option --ref-mask-index requires --ref-mask-name.')
+
+    if opts.ref_mask_name is not None:
+        config.workflow.ref_mask_name = opts.ref_mask_name
+    if opts.ref_mask_index is not None:
+        config.workflow.ref_mask_index = tuple(opts.ref_mask_index)
+
     if opts.pvc_tool is not None:
         config.workflow.pvc_tool = opts.pvc_tool
     if opts.pvc_method is not None:

petprep/cli/tests/test_parser.py

@@ -283,3 +283,25 @@ def test_pvc_invalid_method(tmp_path, minimal_bids):
     with pytest.raises(SystemExit):
         parse_args(args=args)
     _reset_config()
+
+
+def test_reference_mask_options(tmp_path, minimal_bids, monkeypatch):
+    work_dir = tmp_path / 'work'
+    base_args = [
+        str(minimal_bids),
+        str(tmp_path / 'out'),
+        'participant',
+        '-w',
+        str(work_dir),
+        '--skip-bids-validation',
+    ]
+
+    # Missing --ref-mask-name should raise error when --ref-mask-index is used
+    with pytest.raises(SystemExit):
+        parse_args(args=base_args + ['--ref-mask-index', '1', '2'])
+    _reset_config()
+
+    parse_args(args=base_args + ['--ref-mask-name', 'cerebellum', '--ref-mask-index', '3', '4'])
+    assert config.workflow.ref_mask_name == 'cerebellum'
+    assert config.workflow.ref_mask_index == (3, 4)
+    _reset_config()

petprep/config.py

@@ -607,6 +607,12 @@ class workflow(_Config):
     pvc_psf: tuple | None = None
     """Point spread function (PSF) FWHM in mm for PVC."""
 
+    ref_mask_name: str | None = None
+    """Name of the reference region mask to generate."""
+
+    ref_mask_index: tuple | None = None
+    """Tuple of label indices used to build the reference mask."""
+
 class loggers:
     """Keep loggers easily accessible (see :py:func:`init`)."""
 

petprep/data/io_spec.json

@@ -42,9 +42,9 @@
     }
   },
   "patterns": [
-    "sub-{subject}[/ses-{session}]/{datatype<func|pet>|pet}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_res-{res}][_label-{label}][_echo-{echo}][_space-{space}][_desc-{desc}]_{suffix<bold|boldref|pet|petref|dseg|mask>}.{extension<nii|nii.gz|json>|nii.gz}",
+    "sub-{subject}[/ses-{session}]/{datatype<func|pet>|pet}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_ref-{ref}][_res-{res}][_label-{label}][_echo-{echo}][_space-{space}][_desc-{desc}]_{suffix<bold|boldref|pet|petref|dseg|mask>}.{extension<nii|nii.gz|json>|nii.gz}",
     "sub-{subject}[/ses-{session}]/{datatype<func|pet>|pet}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_run-{run}]_from-{from}_to-{to}_mode-{mode<image|points>|image}_{suffix<xfm>|xfm}.{extension<txt|h5>}",
-    "sub-{subject}[/ses-{session}]/{datatype<func|pet>|pet}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_run-{run}][_part-{part}][_desc-{desc}]_{suffix<timeseries>}.{extension<tsv|json>}",
-    "sub-{subject}[/ses-{session}]/{datatype<func|pet>|pet}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_space-{space}][_res-{res}][_den-{den}][_hemi-{hemi}[_label-{label}][_desc-{desc}]_{suffix<|boldref|petref|dseg|mask>}.{extension<dtseries.nii|dtseries.json>}"
+    "sub-{subject}[/ses-{session}]/{datatype<func|pet>|pet}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_ref-{ref}][_run-{run}][_part-{part}][_desc-{desc}]_{suffix<timeseries>}.{extension<tsv|json>}",
+    "sub-{subject}[/ses-{session}]/{datatype<func|pet>|pet}/sub-{subject}[_ses-{session}][_task-{task}][_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_ref-{ref}][_space-{space}][_res-{res}][_den-{den}][_hemi-{hemi}[_label-{label}][_desc-{desc}]_{suffix<|boldref|petref|dseg|mask>}.{extension<dtseries.nii|dtseries.json>}"
   ]
 }

petprep/data/reference_mask/config.json

@@ -0,0 +1,22 @@
+{
+  "gtm": {
+    "cerebellum": 
+    {
+    "refmask_indices": [47, 8], 
+    "erode_by_voxels": 1,
+    "exclude_indices": [7, 46, 172, 1007, 2007, 1009, 2009, 1011, 2011, 1013, 2013, 16], 
+    "dilate_by_voxels": 3
+    }
+  },
+  "wm": {
+    "semiovale":
+    {
+      "refmask_indices": [5001, 5002], 
+      "erode_by_voxels": 1,
+      "exclude_indices": [10, 49, 12, 51, 13, 52, 11, 50, 4, 43, 31, 63],
+      "dilate_by_voxels": 3, 
+      "smooth_fwhm_mm": 10,
+      "target_volume_ml": 40
+    }
+  }
+}

petprep/data/reports-spec-pet.yml

@@ -9,6 +9,10 @@ sections:
   - bids: {datatype: figures, desc: carpetplot, suffix: pet}
   - bids: {datatype: figures, desc: confoundcorr, suffix: pet}
   - bids: {datatype: figures, desc: coreg, suffix: pet}
+  - bids: {datatype: figures, desc: refmask, suffix: pet}
+    caption: Reference region mask overlaid on the PET reference and anatomical data.
+    static: false
+    subtitle: Reference mask check
 - name: About
   reportlets:
   - bids: {datatype: figures, desc: about, suffix: T1w}

petprep/data/reports-spec.yml

@@ -138,6 +138,10 @@ sections:
     caption: PET to anatomical alignment check
     static: false
     subtitle: Additional PET Visualizations
+  - bids: {datatype: figures, desc: refmask, suffix: pet}
+    caption: Reference region mask overlaid on the PET reference and anatomical.
+    static: false
+    subtitle: Reference mask check
 - name: About
   nested: true
   reportlets:

petprep/interfaces/reference_mask.py

@@ -0,0 +1,57 @@
+from nipype.interfaces.base import (
+    SimpleInterface,
+    BaseInterfaceInputSpec,
+    TraitedSpec,
+    File,
+    traits,
+    isdefined,
+)
+import os
+import nibabel as nib
+import json
+
+
+class ExtractRefRegionInputSpec(BaseInterfaceInputSpec):
+    seg_file = File(exists=True, mandatory=True, desc="Segmentation NIfTI file")
+    config_file = File(exists=True, mandatory=True, desc="Path to the config.json file")
+    segmentation_type = traits.Str(mandatory=True, desc="Type of segmentation (e.g. 'gtm', 'wm')")
+    region_name = traits.Str(
+        mandatory=True, desc="Name of the reference region (e.g. 'cerebellum')"
+    )
+    override_indices = traits.List(traits.Int, desc='Use these indices instead of configuration')
+
+
+class ExtractRefRegionOutputSpec(TraitedSpec):
+    refmask_file = File(exists=True, desc="Output reference mask NIfTI file")
+
+
+class ExtractRefRegion(SimpleInterface):
+    input_spec = ExtractRefRegionInputSpec
+    output_spec = ExtractRefRegionOutputSpec
+
+    def _run_interface(self, runtime):
+        seg_img = nib.load(self.inputs.seg_file)
+
+        if isdefined(self.inputs.override_indices):
+            cfg = {'refmask_indices': list(self.inputs.override_indices)}
+        else:
+            # Load the config
+            with open(self.inputs.config_file, 'r') as f:
+                config = json.load(f)
+
+            try:
+                cfg = config[self.inputs.segmentation_type][self.inputs.region_name]
+            except KeyError:
+                raise ValueError(
+                    f"Configuration not found for segmentation='{self.inputs.segmentation_type}' "
+                    f"and region='{self.inputs.region_name}'"
+                )
+
+        from petprep.utils.reference_mask import generate_reference_region
+
+        refmask_img = generate_reference_region(seg_img=seg_img, config=cfg)
+
+        out_file = os.path.abspath("refmask.nii.gz")
+        nib.save(refmask_img, out_file)
+        self._results["refmask_file"] = out_file
+        return runtime

petprep/interfaces/tests/test_reference_mask.py

@@ -0,0 +1,116 @@
+import json
+from pathlib import Path
+
+import nibabel as nb
+import numpy as np
+from nipype.pipeline import engine as pe
+
+from ..reference_mask import ExtractRefRegion
+
+
+def _create_seg(tmp_path: Path) -> Path:
+    data = np.zeros((5, 5, 5), dtype=np.uint8)
+    data[1, 1, 1] = 1
+    data[2, 2, 2] = 2
+    img = nb.Nifti1Image(data, np.eye(4))
+    seg_file = tmp_path / "seg.nii.gz"
+    img.to_filename(seg_file)
+    return seg_file
+
+
+def _create_config(tmp_path: Path, indices, extra=None):
+    cfg = {"testseg": {"region": {"refmask_indices": indices}}}
+    if extra:
+        cfg["testseg"]["region"].update(extra)
+    cfg_file = tmp_path / "config.json"
+    cfg_file.write_text(json.dumps(cfg))
+    return cfg_file
+
+
+def test_extract_refregion(tmp_path):
+    seg = _create_seg(tmp_path)
+    cfg = _create_config(tmp_path, [1])
+
+    node = pe.Node(
+        ExtractRefRegion(
+            seg_file=str(seg),
+            config_file=str(cfg),
+            segmentation_type="testseg",
+            region_name="region",
+        ),
+        name="er", base_dir=str(tmp_path)
+    )
+    res = node.run()
+    out = nb.load(res.outputs.refmask_file).get_fdata()
+    assert out[1, 1, 1] == 1
+    assert out.sum() == 1
+
+
+def test_extract_refregion_override(tmp_path):
+    seg = _create_seg(tmp_path)
+    cfg = _create_config(tmp_path, [1])
+
+    node = pe.Node(
+        ExtractRefRegion(
+            seg_file=str(seg),
+            config_file=str(cfg),
+            segmentation_type="testseg",
+            region_name="region",
+            override_indices=[2],
+        ),
+        name="er2", base_dir=str(tmp_path)
+    )
+    res = node.run()
+    out = nb.load(res.outputs.refmask_file).get_fdata()
+    assert out[2, 2, 2] == 1
+    assert out.sum() == 1
+
+
+def test_extract_refregion_override_ignores_config(tmp_path):
+    seg = _create_seg(tmp_path)
+    cfg = _create_config(
+        tmp_path,
+        [1],
+        {
+            "exclude_indices": [2],
+            "erode_by_voxels": 1,
+        },
+    )
+
+    node = pe.Node(
+        ExtractRefRegion(
+            seg_file=str(seg),
+            config_file=str(cfg),
+            segmentation_type="testseg",
+            region_name="region",
+            override_indices=[2],
+        ),
+        name='er3',
+        base_dir=str(tmp_path),
+    )
+    res = node.run()
+    out = nb.load(res.outputs.refmask_file).get_fdata()
+    assert out[2, 2, 2] == 1
+    assert out.sum() == 1
+
+
+def test_extract_refregion_override_missing_config(tmp_path):
+    seg = _create_seg(tmp_path)
+    cfg_file = tmp_path / 'config.json'
+    cfg_file.write_text(json.dumps({'testseg': {}}))
+
+    node = pe.Node(
+        ExtractRefRegion(
+            seg_file=str(seg),
+            config_file=str(cfg_file),
+            segmentation_type='testseg',
+            region_name='missing',
+            override_indices=[2],
+        ),
+        name='er4',
+        base_dir=str(tmp_path),
+    )
+    res = node.run()
+    out = nb.load(res.outputs.refmask_file).get_fdata()
+    assert out[2, 2, 2] == 1
+    assert out.sum() == 1