Merge pull request #6 from mnoergaard/main

refactor fMRIPrep to PETPrep

Still need modifications to niworkflows such as the BIDSDataGrabber and also the nipreps.json file.

Author: mnoergaard

fmriprep/cli/workflow.py

@@ -86,7 +86,7 @@ def build_workflow(config_file, retval):
     if config.execution.reports_only:
         build_log.log(25, 'Running --reports-only on participants %s', ', '.join(subject_list))
         session_list = (
-            config.execution.bids_filters.get('bold', {}).get('session')
+            config.execution.bids_filters.get('pet', config.execution.bids_filters.get('bold', {})).get('session')
             if config.execution.bids_filters
             else None
         )

fmriprep/data/reports-spec-pet.yml

@@ -0,0 +1,14 @@
+package: fmriprep
+title: PET report for participant '{subject}', session '{session}' - fMRIPrep
+sections:
+- name: PET
+  ordering: session
+  reportlets:
+  - bids: {datatype: figures, desc: summary, suffix: pet}
+  - bids: {datatype: figures, desc: validation, suffix: pet}
+  - bids: {datatype: figures, desc: carpetplot, suffix: pet}
+  - bids: {datatype: figures, desc: confoundcorr, suffix: pet}
+  - bids: {datatype: figures, desc: coreg, suffix: pet}
+- name: About
+  reportlets:
+  - bids: {datatype: figures, desc: about, suffix: T1w}

fmriprep/data/reports-spec.yml

@@ -1,5 +1,5 @@
 package: fmriprep
-title: Visual report for participant '{subject}' - fMRIPrep
+title: Visual report for participant '{subject}' - PETPrep
 sections:
 - name: Summary
   reportlets:
@@ -106,6 +106,38 @@ sections:
       effects and can inform decisions about feature orthogonalization prior to
       confound regression.
     subtitle: Correlations among nuisance regressors
+- name: PET
+  ordering: session,task,acquisition,ceagent,reconstruction,direction,run
+  reportlets:
+  - bids: {datatype: figures, desc: summary, suffix: pet}
+    caption: Summary of PET data acquisition parameters and processing workflow overview, including details such as injected dose, radiotracer used, and scan duration.
+    static: true
+    subtitle: PET Acquisition and Workflow Summary
+
+  - bids: {datatype: figures, desc: validation, suffix: pet}
+    caption: Validation of PET images against BIDS specifications and initial quality assessment including checks for missing slices, artifacts, and alignment issues.
+    static: true
+    subtitle: PET Data Validation
+
+  - bids: {datatype: figures, desc: carpetplot, suffix: pet}
+    caption: |
+      Summary statistics and global PET signal measures are presented.
+      A carpet plot displays voxel-level PET tracer uptake over time within the brain mask. Global signals calculated across the whole-brain (GS), white matter (WM), and cerebrospinal fluid (CSF) regions are plotted, along with DVARS and framewise displacement (FD) to visualize potential motion or acquisition artifacts.
+      "Ctx" = cortex, "Cb" = cerebellum, "WM" = white matter, "CSF" = cerebrospinal fluid.
+    static: false
+    subtitle: PET Summary and Carpet Plot
+
+  - bids: {datatype: figures, desc: confoundcorr, suffix: pet}
+    caption: |
+      Left: Correlation heatmap illustrating relationships among PET-derived confound variables (e.g., motion parameters, global signal).
+      Right: Magnitude of correlation between each PET confound time series and the global PET signal. High correlations suggest potential partial volume effects or motion-induced artifacts, informing subsequent confound regression strategies.
+    static: false
+    subtitle: PET Confound Correlation
+
+  - bids: {datatype: figures, desc: coreg, suffix: pet}
+    caption: PET to anatomical alignment check
+    static: false
+    subtitle: Additional PET Visualizations
 - name: About
   nested: true
   reportlets:

fmriprep/reports/core.py

@@ -119,7 +119,7 @@ def generate_reports(
             # we separate the functional reports per session
             if session_list is None:
                 all_filters = config.execution.bids_filters or {}
-                filters = all_filters.get('bold', {})
+                filters = all_filters.get("pet", all_filters.get("bold", {}))
                 session_list = config.execution.layout.get_sessions(
                     subject=subject_label, **filters
                 )
@@ -145,4 +145,21 @@ def generate_reports(
                 if report_error is not None:
                     errors.append(report_error)
 
+                bootstrap_file = data.load('reports-spec-pet.yml')
+                html_report = f'sub-{subject_label}_ses-{session_label}_pet.html'
+
+                report_error = run_reports(
+                    output_dir,
+                    subject_label,
+                    run_uuid,
+                    bootstrap_file=bootstrap_file,
+                    out_filename=html_report,
+                    reportlets_dir=reportlets_dir,
+                    errorname=f'report-{run_uuid}-{subject_label}-pet.err',
+                    subject=subject_label,
+                    session=session_label,
+                )
+                if report_error is not None:
+                    errors.append(report_error)
+
     return errors

fmriprep/reports/tests/test_reports.py

@@ -23,6 +23,10 @@
                 'sub-001_ses-003_func.html',
                 'sub-001_ses-004_func.html',
                 'sub-001_ses-005_func.html',
+                'sub-001_ses-001_pet.html',
+                'sub-001_ses-003_pet.html',
+                'sub-001_ses-004_pet.html',
+                'sub-001_ses-005_pet.html',
             ],
         ),
         (4, ['sub-001.html']),
@@ -109,3 +113,30 @@ def mock_session_list(*args, **kwargs):
         assert 'One or more execution steps failed' in html_content, (
             f'The file {expected_files[0]} did not contain the reported error.'
         )
+
+
+def test_pet_report(tmp_path, monkeypatch):
+    fake_uuid = 'fake_uuid'
+
+    pet_source = data_dir / 'work/reportlets/fmriprep'
+    sub_dir = tmp_path / 'sub-01' / 'figures'
+    sub_dir.mkdir(parents=True)
+
+    shutil.copy2(pet_source / 'sub-001/figures/sub-001_desc-about_T1w.html', sub_dir / 'sub-01_desc-about_T1w.html')
+    shutil.copy2(pet_source / 'sub-001/figures/sub-001_ses-001_task-qct_dir-LR_part-mag_desc-summary_bold.html', sub_dir / 'sub-01_ses-baseline_desc-summary_pet.html')
+    shutil.copy2(pet_source / 'sub-001/figures/sub-001_ses-001_task-qct_dir-LR_part-mag_desc-validation_bold.html', sub_dir / 'sub-01_ses-baseline_desc-validation_pet.html')
+    shutil.copy2(pet_source / 'sub-001/figures/sub-001_ses-001_task-qct_dir-LR_part-mag_desc-carpetplot_bold.svg', sub_dir / 'sub-01_ses-baseline_desc-carpetplot_pet.svg')
+    shutil.copy2(pet_source / 'sub-001/figures/sub-001_ses-001_task-qct_dir-LR_part-mag_desc-confoundcorr_bold.svg', sub_dir / 'sub-01_ses-baseline_desc-confoundcorr_pet.svg')
+    shutil.copy2(pet_source / 'sub-01/func/sub-01_task-mixedgamblestask_run-01_bold_bbr.svg', sub_dir / 'sub-01_ses-baseline_pet.svg')
+
+    config.execution.aggr_ses_reports = 4
+    config.execution.layout = BIDSLayout(data_dir / 'pet')
+    monkeypatch.setattr(config.execution, 'bids_filters', {'pet': {'session': ['baseline']}})
+
+    failed_reports = generate_reports(['01'], tmp_path, fake_uuid)
+
+    assert not failed_reports
+    html_file = tmp_path / 'sub-01.html'
+    assert html_file.is_file()
+    html_content = html_file.read_text()
+    assert '<div id="PET"' in html_content

fmriprep/workflows/pet/base.py

@@ -152,7 +152,7 @@ def init_pet_wf(
         precomputed = {}
     pet_file = pet_series
 
-    fmriprep_dir = config.execution.petprep_dir
+    petprep_dir = config.execution.petprep_dir
     omp_nthreads = config.nipype.omp_nthreads
     all_metadata = [config.execution.layout.get_metadata(file) for file in pet_series]
 
@@ -272,7 +272,7 @@ def init_pet_wf(
     if petref_out:
         ds_pet_native_wf = init_ds_pet_native_wf(
             bids_root=str(config.execution.bids_dir),
-            output_dir=fmriprep_dir,
+            output_dir=petprep_dir,
             pet_output=petref_out,
             all_metadata=all_metadata,
         )
@@ -288,7 +288,7 @@ def init_pet_wf(
         # Fill-in datasinks of reportlets seen so far
         for node in workflow.list_node_names():
             if node.split('.')[-1].startswith('ds_report'):
-                workflow.get_node(node).inputs.base_directory = fmriprep_dir
+                workflow.get_node(node).inputs.base_directory = petprep_dir
                 workflow.get_node(node).inputs.source_file = pet_file
         return workflow
 
@@ -319,7 +319,7 @@ def init_pet_wf(
     if nonstd_spaces.intersection(('anat', 'T1w')):
         ds_pet_t1_wf = init_ds_volumes_wf(
             bids_root=str(config.execution.bids_dir),
-            output_dir=fmriprep_dir,
+            output_dir=petprep_dir,
             metadata=all_metadata[0],
             name='ds_pet_t1_wf',
         )
@@ -350,7 +350,7 @@ def init_pet_wf(
         )
         ds_pet_std_wf = init_ds_volumes_wf(
             bids_root=str(config.execution.bids_dir),
-            output_dir=fmriprep_dir,
+            output_dir=petprep_dir,
             metadata=all_metadata[0],
             name='ds_pet_std_wf',
         )
@@ -401,7 +401,7 @@ def init_pet_wf(
             surface_spaces=freesurfer_spaces,
             medial_surface_nan=config.workflow.medial_surface_nan,
             metadata=all_metadata[0],
-            output_dir=fmriprep_dir,
+            output_dir=petprep_dir,
             name='pet_surf_wf',
         )
         pet_surf_wf.inputs.inputnode.source_file = pet_file
@@ -459,7 +459,7 @@ def init_pet_wf(
 
         ds_pet_cifti = pe.Node(
             DerivativesDataSink(
-                base_directory=fmriprep_dir,
+                base_directory=petprep_dir,
                 space='fsLR',
                 density=config.workflow.cifti_output,
                 suffix='pet',
@@ -523,7 +523,7 @@ def init_pet_wf(
 
     ds_confounds = pe.Node(
         DerivativesDataSink(
-            base_directory=fmriprep_dir,
+            base_directory=petprep_dir,
             desc='confounds',
             suffix='timeseries',
         ),
@@ -589,7 +589,7 @@ def _last(inlist):
     # Fill-in datasinks of reportlets seen so far
     for node in workflow.list_node_names():
         if node.split('.')[-1].startswith('ds_report'):
-            workflow.get_node(node).inputs.base_directory = fmriprep_dir
+            workflow.get_node(node).inputs.base_directory = petprep_dir
             workflow.get_node(node).inputs.source_file = pet_file
 
     return workflow

fmriprep/workflows/pet/confounds.py

@@ -274,7 +274,7 @@ def init_pet_confs_wf(
         iterfield=['in_file'],
     )
     merge_rois = pe.Node(
-        niu.Merge(3, ravel_inputs=True), name='merge_rois', run_without_submitting=True
+        niu.Merge(4, ravel_inputs=True), name='merge_rois', run_without_submitting=True
     )
     signals = pe.Node(
         SignalExtraction(class_labels=signals_class_labels), name='signals', mem_gb=mem_gb
@@ -350,18 +350,6 @@ def init_pet_confs_wf(
         mem_gb=DEFAULT_MEMORY_MIN_GB,
     )
 
-    def _last(inlist):
-        return inlist[-1]
-
-    def _select_cols(table):
-        import pandas as pd
-
-        return [
-            col
-            for col in pd.read_table(table, nrows=2).columns
-            if not col.startswith(('a_comp_cor_', 't_comp_cor_', 'std_dvars'))
-        ]
-
     workflow.connect([
         # connect inputnode to each non-anatomical confound node
         (inputnode, dvars, [('pet', 'in_file'),
@@ -392,10 +380,11 @@ def _select_cols(table):
             ('petref2anat_xfm', 'transforms'),
         ]),
         (acompcor_tfm, acompcor_bin, [('output_image', 'in_file')]),
+        (union_mask, merge_rois, [('out', 'in1')]),
         (acompcor_bin, merge_rois, [
-            (('out_mask', _last), 'in3'),
-            (('out_mask', lambda masks: masks[0]), 'in1'),
-            (('out_mask', lambda masks: masks[1]), 'in2'),
+            (('out_mask', _first), 'in2'),
+            (('out_mask', _second), 'in3'),
+            (('out_mask', _last), 'in4'),
         ]),
         (merge_rois, signals, [('out', 'label_files')]),
 
@@ -590,6 +579,34 @@ def _binary_union(mask1, mask2):
     return str(out_name)
 
 
+def _smooth_binarize(in_file, fwhm=10.0, thresh=0.2):
+    """Smooth ``in_file`` with a Gaussian kernel, binarize and keep largest cluster."""
+    from pathlib import Path
+
+    import nibabel as nb
+    import numpy as np
+    from scipy.ndimage import gaussian_filter, label
+
+    img = nb.load(in_file)
+    data = img.get_fdata(dtype=np.float32)
+    zooms = np.array(img.header.get_zooms()[:3], dtype=float)
+    sigma = (fwhm / 2.3548) / zooms
+    smoothed = gaussian_filter(data, sigma=sigma)
+    mask = smoothed > (thresh * smoothed.max())
+
+    labeled, n_labels = label(mask)
+    if n_labels > 1:
+        sizes = np.bincount(labeled.ravel())
+        sizes[0] = 0  # ignore background
+        mask = labeled == sizes.argmax()
+
+    out_img = img.__class__(mask.astype('uint8'), img.affine, img.header)
+    out_img.set_data_dtype('uint8')
+    out_name = Path('smoothed_bin_mask.nii.gz').absolute()
+    out_img.to_filename(out_name)
+    return str(out_name)
+
+
 def _carpet_parcellation(segmentation, crown_mask, nifti=False):
     """Generate a segmentation for carpet plot visualization."""
     from pathlib import Path
@@ -619,3 +636,31 @@ def _get_zooms(in_file):
     import nibabel as nb
 
     return tuple(nb.load(in_file).header.get_zooms()[:3])
+
+
+def _last(inlist):
+    """Return the last element of a list."""
+
+    return inlist[-1]
+
+
+def _first(inlist):
+    """Return the first element of a list."""
+
+    return inlist[0]
+
+
+def _second(inlist):
+    """Return the second element of a list."""
+
+    return inlist[1]
+
+def _select_cols(table):
+    """Return confound columns excluding a/tCompCor and std_dvars."""
+    import pandas as pd
+
+    return [
+        col
+        for col in pd.read_table(table, nrows=2).columns
+        if not col.startswith(('a_comp_cor_', 't_comp_cor_', 'std_dvars'))
+    ]

fmriprep/workflows/pet/fit.py

@@ -228,6 +228,7 @@ def init_pet_fit_wf(
             ('subjects_dir', 'inputnode.subjects_dir'),
             ('subject_id', 'inputnode.subject_id'),
         ]),
+        (petref_buffer, func_fit_reports_wf, [('petref', 'inputnode.petref')]),
         (outputnode, func_fit_reports_wf, [
             ('pet_mask', 'inputnode.pet_mask'),
             ('petref2anat_xfm', 'inputnode.petref2anat_xfm'),
@@ -356,15 +357,18 @@ def init_pet_fit_wf(
 
     # Stage 4: Estimate PET brain mask
     from niworkflows.interfaces.fixes import FixHeaderApplyTransforms as ApplyTransforms
-    from niworkflows.interfaces.nibabel import Binarize
 
-    from .confounds import _binary_union
+    from .confounds import _binary_union, _smooth_binarize
 
     t1w_mask_tfm = pe.Node(
         ApplyTransforms(interpolation='MultiLabel', invert_transform_flags=[True]),
         name='t1w_mask_tfm',
     )
-    petref_mask = pe.Node(Binarize(thresh_low=0.2), name='petref_mask')
+    petref_mask = pe.Node(
+        niu.Function(function=_smooth_binarize), name='petref_mask'
+    )
+    petref_mask.inputs.fwhm = 10.0
+    petref_mask.inputs.thresh = 0.2
     merge_mask = pe.Node(niu.Function(function=_binary_union), name='merge_mask')
 
     if not petref2anat_xform:
@@ -379,7 +383,7 @@ def init_pet_fit_wf(
             (inputnode, t1w_mask_tfm, [('t1w_mask', 'input_image')]),
             (petref_buffer, t1w_mask_tfm, [('petref', 'reference_image')]),
             (petref_buffer, petref_mask, [('petref', 'in_file')]),
-            (petref_mask, merge_mask, [('out_mask', 'mask1')]),
+            (petref_mask, merge_mask, [('out', 'mask1')]),
             (t1w_mask_tfm, merge_mask, [('output_image', 'mask2')]),
             (merge_mask, outputnode, [('out', 'pet_mask')]),
         ]

fmriprep/workflows/pet/outputs.py

@@ -307,7 +307,7 @@ def init_func_fit_reports_wf(
     pet_t1_report = pe.Node(
         SimpleBeforeAfter(
             before_label='T1w',
-            after_label='EPI',
+            after_label='PET',
             dismiss_affine=True,
         ),
         name='pet_t1_report',
@@ -317,6 +317,7 @@ def init_func_fit_reports_wf(
     ds_pet_t1_report = pe.Node(
         DerivativesDataSink(
             base_directory=output_dir,
+            desc='coreg',
             suffix='pet',
             datatype='figures',
         ),