Core function to extract a timepoint from a 4D file (Inria-Empenn#239)

* Core function to extract a timepoint from a 4D file

* Codespell

* New versions for actions

* New versions for actions

* Extract timepoints range

* Use core function instead of FSL's ExtractROI for SPM pipelines

Author: bclenet

.github/workflows/code_quality.yml

@@ -32,7 +32,7 @@ jobs:
       with:
         python-version: 3.9
 
-    - uses: actions/cache@v3
+    - uses: actions/cache@v4
       with:
         path: ~/.cache/pip
         key: ${{ runner.os }}-pip-${{ hashFiles('setup.py') }}

.github/workflows/pipeline_status.yml

@@ -37,7 +37,7 @@ jobs:
       with:
         python-version: 3.9
 
-    - uses: actions/cache@v3
+    - uses: actions/cache@v4
       with:
         path: ~/.cache/pip
         key: ${{ runner.os }}-pip-${{ hashFiles('setup.py') }}

docs/core.md

@@ -124,6 +124,17 @@ This module contains a set of functions dedicated to computations on images.
 # Get dimensions of voxels along x, y, and z in mm (returns e.g.: [1.0, 1.0, 1.0]).
 get_voxel_dimensions('/path/to/the/image.nii.gz')
 ```
+
+ * `get_image_timepoint` : extracts the 3D volume in a given time point range from a 4D Nifti image
+
+```python
+# Create a nifti 3D file containing time point 5 of the input image (returns a path to the generated 3D image)
+get_image_timepoints('/path/to/the/image.nii.gz', 5, 5):
+
+# Create a nifti 3D file containing time points 5 to 10 of the input image (returns a path to the generated 3D image)
+get_image_timepoints('/path/to/the/image.nii.gz', 5, 10):
+```
+
 ## narps_open.core.interfaces
 
 This module contains a set of interface creators inheriting form the `narps_open.core.interfaces.InterfaceCreator` abstract class.

narps_open/core/image.py

@@ -23,3 +23,32 @@ def get_voxel_dimensions(image: str) -> list:
         float(voxel_dimensions[1]),
         float(voxel_dimensions[2])
         ]
+
+def get_image_timepoints(in_file: str, start_time_point: int, end_time_point: int) -> str:
+    """
+    Extract the 3D volume at a given time point from a 4D Nifti image.
+    Create a Nifti file containing the 3D volume.
+    Return the filename of the created 3D output file.
+
+    Arguments:
+        in_file: str, string that represent an absolute path to a Nifti 4D image.
+        start_time_point: int, zero-based index of the first volume to extract from the 4D image.
+        end_time_point: int, zero-based index of the last volume to extract from the 4D image.
+
+    Returns:
+        str, path to the created file
+    """
+    # These imports must stay inside the function, as required by Nipype
+    from os import extsep
+    from os.path import abspath, basename
+    import nibabel as nib
+
+    # The output filename is based on the base filename of the input file.
+    # Note that we use abspath to write the file in the base_directory of a nipype Node.
+    out_file = abspath(
+        basename(in_file).split(extsep)[0]+f'_timepoint-{start_time_point}-{end_time_point}.nii')
+
+    # Perform timepoints extraction
+    nib.save(nib.load(in_file).slicer[..., start_time_point:end_time_point+1], out_file)
+
+    return out_file

narps_open/pipelines/team_0H5E.py

@@ -17,16 +17,14 @@
    EstimateModel, EstimateContrast, Threshold
    )
 from nipype.interfaces.spm.base import Info as SPMInfo
-from nipype.interfaces.fsl import (
-    ExtractROI
-    )
 
 from narps_open.pipelines import Pipeline
 from narps_open.data.task import TaskInformation
 from narps_open.data.participants import get_group
 from narps_open.core.common import (
     remove_parent_directory, list_intersection, elements_in_string, clean_list
     )
+from narps_open.core.image import get_image_timepoints
 from narps_open.utils.configuration import Configuration
 
 class PipelineTeam0H5E(Pipeline):
@@ -93,10 +91,13 @@ def get_preprocessing(self):
 
         # EXTRACTROI - remove first image of func
         # > Removal of "dummy" scans (deleting the first four volumes from each run)
-        remove_first_image = Node(ExtractROI(), name = 'remove_first_image')
-        remove_first_image.inputs.t_min = 4
-        remove_first_image.inputs.t_size = 449 # number of time points - 4
-        remove_first_image.inputs.output_type='NIFTI'
+        remove_first_image = Node(Function(
+            function = get_image_timepoints,
+            input_names = ['in_file', 'start_time_point', 'end_time_point'],
+            output_names = ['roi_file']
+            ), name = 'remove_first_image')
+        remove_first_image.inputs.start_time_point = 4
+        remove_first_image.inputs.end_time_point = 453 # last image
         preprocessing.connect(gunzip_func, 'out_file', remove_first_image, 'in_file')
 
         # REALIGN - motion correction

narps_open/pipelines/team_98BT.py

@@ -16,7 +16,6 @@
     DARTELNorm2MNI, FieldMap, Threshold
     )
 from nipype.interfaces.spm.base import Info as SPMInfo
-from nipype.interfaces.fsl import ExtractROI
 from nipype.algorithms.modelgen import SpecifySPMModel
 from niflow.nipype1.workflows.fmri.spm import create_DARTEL_template
 
@@ -26,6 +25,7 @@
 from narps_open.core.common import (
     remove_parent_directory, list_intersection, elements_in_string, clean_list
     )
+from narps_open.core.image import get_image_timepoints
 from narps_open.utils.configuration import Configuration
 
 class PipelineTeam98BT(Pipeline):
@@ -264,10 +264,13 @@ def get_preprocessing_sub_workflow(self):
         preprocessing.connect(slice_timing, 'timecorrected_files', motion_correction, 'in_files')
 
         # Intrasubject coregistration
-        extract_first = Node(ExtractROI(), name = 'extract_first')
-        extract_first.inputs.t_min = 1
-        extract_first.inputs.t_size = 1
-        extract_first.inputs.output_type = 'NIFTI'
+        extract_first = Node(Function(
+            function = get_image_timepoints,
+            input_names = ['in_file', 'start_time_point', 'end_time_point'],
+            output_names = ['roi_file']
+            ), name = 'extract_first')
+        extract_first.inputs.start_time_point = 1
+        extract_first.inputs.end_time_point = 1
         preprocessing.connect(
             motion_correction, 'realigned_unwarped_files', extract_first, 'in_file')
 

narps_open/pipelines/team_V55J.py

@@ -14,7 +14,6 @@
     Coregister, Smooth, OneSampleTTestDesign, EstimateModel, EstimateContrast,
     Level1Design, TwoSampleTTestDesign, RealignUnwarp,
     Normalize12, NewSegment, FieldMap, Threshold)
-from nipype.interfaces.fsl import ExtractROI
 from nipype.algorithms.modelgen import SpecifySPMModel
 from nipype.interfaces.spm.base import Info as SPMInfo
 
@@ -24,6 +23,7 @@
 from narps_open.core.common import (
     remove_parent_directory, list_intersection, elements_in_string, clean_list
     )
+from narps_open.core.image import get_image_timepoints
 from narps_open.utils.configuration import Configuration
 
 class PipelineTeamV55J(Pipeline):
@@ -100,10 +100,13 @@ def get_preprocessing(self):
 
         # EXTRACTROI - get the image 10 in func file
         # "For each run, we selected image 10 to distortion correct and asked to match the VDM file"
-        extract_tenth_image = Node(ExtractROI(), name = 'extract_tenth_image')
-        extract_tenth_image.inputs.t_min = 10
-        extract_tenth_image.inputs.t_size = 1
-        extract_tenth_image.inputs.output_type='NIFTI'
+        extract_tenth_image = Node(Function(
+            function = get_image_timepoints,
+            input_names = ['in_file', 'start_time_point', 'end_time_point'],
+            output_names = ['roi_file']
+            ), name = 'extract_tenth_image')
+        extract_tenth_image.inputs.start_time_point = 9 # 0-based 10th image
+        extract_tenth_image.inputs.end_time_point = 9
         preprocessing.connect(gunzip_func, 'out_file', extract_tenth_image, 'in_file')
 
         # FIELDMAP - Calculate VDM routine of the FieldMap tool in SPM12

tests/core/test_image.py

@@ -11,7 +11,10 @@
     pytest -q test_image.py -k <selected_test>
 """
 
-from os.path import abspath, join
+from os.path import abspath, join, basename, exists
+from filecmp import cmp
+from shutil import copyfile
+
 from numpy import isclose
 
 from pytest import mark
@@ -46,3 +49,67 @@ def test_get_voxel_dimensions():
 
         # Check voxel sizes
         assert isclose(outputs.voxel_dimensions, [8.0, 8.0, 9.6]).all()
+
+    @staticmethod
+    @mark.unit_test
+    def test_get_image_timepoints(temporary_data_dir):
+        """ Test the get_image_timepoints function """
+
+        # Path to the test image
+        test_file_path = abspath(join(
+            Configuration()['directories']['test_data'],
+            'core',
+            'image',
+            'test_timepoint_input.nii.gz'))
+
+        # Create a Nipype Node using get_image_timepoint
+        test_get_image_timepoint = Node(Function(
+            function = im.get_image_timepoints,
+            input_names = ['in_file', 'start_time_point', 'end_time_point'],
+            output_names = ['out_file']
+            ), name = 'test_get_image_timepoint')
+        test_get_image_timepoint.inputs.in_file = test_file_path
+        test_get_image_timepoint.inputs.start_time_point = 10
+        test_get_image_timepoint.inputs.end_time_point = 10
+        test_get_image_timepoint.base_dir = temporary_data_dir
+        outputs = test_get_image_timepoint.run().outputs
+
+        # Check output file name
+        assert basename(outputs.out_file) == 'test_timepoint_input_timepoint-10-10.nii'
+
+        # Check if file exists
+        assert exists(outputs.out_file)
+
+        # Compare with expected
+        out_test_file_path = abspath(join(
+            Configuration()['directories']['test_data'],
+            'core',
+            'image',
+            'test_timepoint_output.nii'))
+        assert cmp(outputs.out_file, out_test_file_path)
+
+        # Try a new time range
+        test_get_image_timepoint = Node(Function(
+            function = im.get_image_timepoints,
+            input_names = ['in_file', 'start_time_point', 'end_time_point'],
+            output_names = ['out_file']
+            ), name = 'test_get_image_timepoint')
+        test_get_image_timepoint.inputs.in_file = test_file_path
+        test_get_image_timepoint.inputs.start_time_point = 2
+        test_get_image_timepoint.inputs.end_time_point = 4
+        test_get_image_timepoint.base_dir = temporary_data_dir
+        outputs = test_get_image_timepoint.run().outputs
+
+        # Check output file name
+        assert basename(outputs.out_file) == 'test_timepoint_input_timepoint-2-4.nii'
+
+        # Check if file exists
+        assert exists(outputs.out_file)
+
+        # Compare with expected
+        out_test_file_path = abspath(join(
+            Configuration()['directories']['test_data'],
+            'core',
+            'image',
+            'test_timepoint_output_2.nii'))
+        assert cmp(outputs.out_file, out_test_file_path)