Merge pull request #90 from bilgelm/pet3d

ENH: enable running petprep on 3D PET images

Author: effigies

.maint/CONTRIBUTORS.md

@@ -19,5 +19,4 @@ Before every release, unlisted contributors will be invited again to add their n
 | Jonghyun | Bae | | 0000-0002-3987-6603 | National Institute on Aging |
 | Legarreta Gorroño | Jon Haitz | @jhlegarreta | 0000-0002-9661-1396 | Brigham and Women's Hospital, Mass General Brigham, Harvard Medical School, MA, USA |
 | Markiewicz | Christopher J. | @effigies | 0000-0002-6533-164X | Department of Psychology, Stanford University, CA, USA |
-| Murat | Bilgel | @bilgelm | 0000-0001-5042-7422 | National Institute on Aging |
 | Schiavone | Alice | @aliswh | 0009-0006-4844-2088 | Neurobiology Research Unit, Rigshospitalet, Copenhagen and Department of Computer Science, University of Copenhagen |

Dockerfile

@@ -203,7 +203,7 @@ ENV PATH="/opt/afni-latest:$PATH" \
 RUN useradd -m -s /bin/bash -G users petprep
 WORKDIR /home/petprep
 ENV HOME="/home/petprep" \
-    LD_LIBRARY_PATH="/usr/lib/x86_64-linux-gnu:$LD_LIBRARY_PATH"
+    LD_LIBRARY_PATH="/usr/lib/x86_64-linux-gnu:${LD_LIBRARY_PATH:-}"
 
 COPY --from=micromamba /bin/micromamba /bin/micromamba
 COPY --from=micromamba /opt/conda/envs/petprep /opt/conda/envs/petprep
@@ -212,7 +212,7 @@ ENV MAMBA_ROOT_PREFIX="/opt/conda"
 RUN micromamba shell init -s bash && \
     echo "micromamba activate petprep" >> $HOME/.bashrc
 ENV PATH="/opt/conda/envs/petprep/bin:$PATH" \
-    CPATH="/opt/conda/envs/petprep/include:$CPATH" \
+    CPATH="/opt/conda/envs/petprep/include:${CPATH:-}" \
     LD_LIBRARY_PATH="/opt/conda/envs/petprep/lib:$LD_LIBRARY_PATH"
 
 # Precaching atlases

petprep/utils/misc.py

@@ -78,7 +78,14 @@ def estimate_pet_mem_usage(pet_fname: str) -> tuple[int, dict]:
     nvox = int(np.prod(img.shape, dtype='u8'))
     # Assume tools will coerce to 8-byte floats to be safe
     pet_size_gb = 8 * nvox / (1024**3)
-    pet_tlen = img.shape[-1]
+
+    if img.ndim == 4:
+        pet_tlen = img.shape[3]
+    elif img.ndim == 3:
+        pet_tlen = 1
+    else:
+        raise ValueError('PET image must be 3D or 4D')
+
     mem_gb = {
         'filesize': pet_size_gb,
         'resampled': pet_size_gb * 4,

petprep/workflows/pet/base.py

@@ -166,11 +166,6 @@ def init_pet_wf(
     all_metadata = [config.execution.layout.get_metadata(file) for file in pet_series]
 
     nvols, mem_gb = estimate_pet_mem_usage(pet_file)
-    if nvols <= 1 - config.execution.sloppy:
-        config.loggers.workflow.warning(
-            f'Too short PET series (<= 5 timepoints). Skipping processing of <{pet_file}>.'
-        )
-        return
 
     config.loggers.workflow.debug(
         'Creating pet processing workflow for <%s> (%.2f GB / %d frames). '
@@ -755,81 +750,82 @@ def init_pet_wf(
             (pet_ref_tacs_wf, ds_ref_tacs, [('outputnode.timeseries', 'in_file')]),
         ])  # fmt:skip
 
-    pet_confounds_wf = init_pet_confs_wf(
-        mem_gb=mem_gb['largemem'],
-        metadata=all_metadata[0],
-        freesurfer=config.workflow.run_reconall,
-        regressors_all_comps=config.workflow.regressors_all_comps,
-        regressors_fd_th=config.workflow.regressors_fd_th,
-        regressors_dvars_th=config.workflow.regressors_dvars_th,
-        name='pet_confounds_wf',
-    )
-
-    ds_confounds = pe.Node(
-        DerivativesDataSink(
-            base_directory=petprep_dir,
-            desc='confounds',
-            suffix='timeseries',
-        ),
-        name='ds_confounds',
-        run_without_submitting=True,
-        mem_gb=config.DEFAULT_MEMORY_MIN_GB,
-    )
-    ds_confounds.inputs.source_file = pet_file
-
-    workflow.connect([
-        (inputnode, pet_confounds_wf, [
-            ('t1w_tpms', 'inputnode.t1w_tpms'),
-            ('t1w_mask', 'inputnode.t1w_mask'),
-        ]),
-        (pet_fit_wf, pet_confounds_wf, [
-            ('outputnode.pet_mask', 'inputnode.pet_mask'),
-            ('outputnode.petref', 'inputnode.petref'),
-            ('outputnode.motion_xfm', 'inputnode.motion_xfm'),
-            ('outputnode.petref2anat_xfm', 'inputnode.petref2anat_xfm'),
-        ]),
-        (pet_native_wf, pet_confounds_wf, [
-            ('outputnode.pet_native', 'inputnode.pet'),
-        ]),
-        (pet_confounds_wf, ds_confounds, [
-            ('outputnode.confounds_file', 'in_file'),
-            ('outputnode.confounds_metadata', 'meta_dict'),
-        ]),
-    ])  # fmt:skip
-
-    if spaces.get_spaces(nonstandard=False, dim=(3,)):
-        carpetplot_wf = init_carpetplot_wf(
-            mem_gb=mem_gb['resampled'],
+    if nvols > 1:  # run these only if 4-D PET
+        pet_confounds_wf = init_pet_confs_wf(
+            mem_gb=mem_gb['largemem'],
             metadata=all_metadata[0],
-            cifti_output=config.workflow.cifti_output,
-            name='carpetplot_wf',
+            freesurfer=config.workflow.run_reconall,
+            regressors_all_comps=config.workflow.regressors_all_comps,
+            regressors_fd_th=config.workflow.regressors_fd_th,
+            regressors_dvars_th=config.workflow.regressors_dvars_th,
+            name='pet_confounds_wf',
         )
 
-        if config.workflow.cifti_output:
-            workflow.connect(
-                pet_grayords_wf, 'outputnode.cifti_pet', carpetplot_wf, 'inputnode.cifti_pet',
-            )  # fmt:skip
-
-        def _last(inlist):
-            return inlist[-1]
+        ds_confounds = pe.Node(
+            DerivativesDataSink(
+                base_directory=petprep_dir,
+                desc='confounds',
+                suffix='timeseries',
+            ),
+            name='ds_confounds',
+            run_without_submitting=True,
+            mem_gb=config.DEFAULT_MEMORY_MIN_GB,
+        )
+        ds_confounds.inputs.source_file = pet_file
 
         workflow.connect([
-            (inputnode, carpetplot_wf, [
-                ('mni2009c2anat_xfm', 'inputnode.std2anat_xfm'),
+            (inputnode, pet_confounds_wf, [
+                ('t1w_tpms', 'inputnode.t1w_tpms'),
+                ('t1w_mask', 'inputnode.t1w_mask'),
             ]),
-            (pet_fit_wf, carpetplot_wf, [
+            (pet_fit_wf, pet_confounds_wf, [
                 ('outputnode.pet_mask', 'inputnode.pet_mask'),
+                ('outputnode.petref', 'inputnode.petref'),
+                ('outputnode.motion_xfm', 'inputnode.motion_xfm'),
                 ('outputnode.petref2anat_xfm', 'inputnode.petref2anat_xfm'),
             ]),
-            (pet_native_wf, carpetplot_wf, [
+            (pet_native_wf, pet_confounds_wf, [
                 ('outputnode.pet_native', 'inputnode.pet'),
             ]),
-            (pet_confounds_wf, carpetplot_wf, [
-                ('outputnode.confounds_file', 'inputnode.confounds_file'),
-                ('outputnode.crown_mask', 'inputnode.crown_mask'),
+            (pet_confounds_wf, ds_confounds, [
+                ('outputnode.confounds_file', 'in_file'),
+                ('outputnode.confounds_metadata', 'meta_dict'),
             ]),
         ])  # fmt:skip
 
+        if spaces.get_spaces(nonstandard=False, dim=(3,)):
+            carpetplot_wf = init_carpetplot_wf(
+                mem_gb=mem_gb['resampled'],
+                metadata=all_metadata[0],
+                cifti_output=config.workflow.cifti_output,
+                name='carpetplot_wf',
+            )
+
+            if config.workflow.cifti_output:
+                workflow.connect(
+                    pet_grayords_wf, 'outputnode.cifti_pet', carpetplot_wf, 'inputnode.cifti_pet',
+                )  # fmt:skip
+
+            def _last(inlist):
+                return inlist[-1]
+
+            workflow.connect([
+                (inputnode, carpetplot_wf, [
+                    ('mni2009c2anat_xfm', 'inputnode.std2anat_xfm'),
+                ]),
+                (pet_fit_wf, carpetplot_wf, [
+                    ('outputnode.pet_mask', 'inputnode.pet_mask'),
+                    ('outputnode.petref2anat_xfm', 'inputnode.petref2anat_xfm'),
+                ]),
+                (pet_native_wf, carpetplot_wf, [
+                    ('outputnode.pet_native', 'inputnode.pet'),
+                ]),
+                (pet_confounds_wf, carpetplot_wf, [
+                    ('outputnode.confounds_file', 'inputnode.confounds_file'),
+                    ('outputnode.crown_mask', 'inputnode.crown_mask'),
+                ]),
+            ])  # fmt:skip
+
     # Fill-in datasinks of reportlets seen so far
     for node in workflow.list_node_names():
         if node.split('.')[-1].startswith('ds_report'):

petprep/workflows/pet/fit.py

@@ -25,6 +25,7 @@
 import nibabel as nb
 from nipype.interfaces import utility as niu
 from nipype.pipeline import engine as pe
+from nitransforms.linear import Affine
 from niworkflows.interfaces.header import ValidateImage
 from niworkflows.utils.connections import listify
 
@@ -198,12 +199,18 @@ def init_pet_fit_wf(
     )
     hmc_buffer = pe.Node(niu.IdentityInterface(fields=['hmc_xforms']), name='hmc_buffer')
 
+    if pet_tlen <= 1:  # 3D PET
+        petref = pet_file
+        idmat_fname = config.execution.work_dir / 'idmat.tfm'
+        Affine().to_filename(idmat_fname, fmt='itk')
+        hmc_xforms = idmat_fname
+        config.loggers.workflow.debug('3D PET file - motion correction not needed')
     if petref:
         petref_buffer.inputs.petref = petref
-        config.loggers.workflow.debug('Reusing motion correction reference: %s', petref)
+        config.loggers.workflow.debug('(Re)using motion correction reference: %s', petref)
     if hmc_xforms:
         hmc_buffer.inputs.hmc_xforms = hmc_xforms
-        config.loggers.workflow.debug('Reusing motion correction transforms: %s', hmc_xforms)
+        config.loggers.workflow.debug('(Re)using motion correction transforms: %s', hmc_xforms)
 
     timing_parameters = prepare_timing_parameters(metadata)
     frame_durations = timing_parameters.get('FrameDuration')