✨Feature/update pixdim4 for the final *bold outputs

CHANGELOG.md

@@ -20,6 +20,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - `pyproject.toml` file with `[build-system]` defined.
 - [![pre-commit.ci status](https://results.pre-commit.ci/badge/github/FCP-INDI/C-PAC/main.svg)](https://results.pre-commit.ci/latest/github/FCP-INDI/C-PAC/main) badge to [`README`](./README.md).
+- validation node to match the pixdim4 of CPAC processed bold outputs with the original raw bold sources.
 
 ### Changed
 

CPAC/pipeline/engine.py

@@ -35,7 +35,12 @@
 from CPAC.pipeline import nipype_pipeline_engine as pe
 from CPAC.pipeline.check_outputs import ExpectedOutputs
 from CPAC.pipeline.nodeblock import NodeBlockFunction
-from CPAC.pipeline.utils import MOVEMENT_FILTER_KEYS, name_fork, source_set
+from CPAC.pipeline.utils import (
+    MOVEMENT_FILTER_KEYS,
+    name_fork,
+    source_set,
+    validate_outputs,
+)
 from CPAC.registration.registration import transform_derivative
 from CPAC.resources.templates.lookup_table import lookup_identifier
 from CPAC.utils.bids_utils import res_in_filename
@@ -1371,6 +1376,7 @@ def gather_pipes(self, wf, cfg, all=False, add_incl=None, add_excl=None):
                 write_json.inputs.json_data = json_info
 
                 wf.connect(id_string, "out_filename", write_json, "filename")
+
                 ds = pe.Node(DataSink(), name=f"sinker_{resource_idx}_{pipe_x}")
                 ds.inputs.parameterization = False
                 ds.inputs.base_directory = out_dct["out_dir"]
@@ -1394,7 +1400,38 @@ def gather_pipes(self, wf, cfg, all=False, add_incl=None, add_excl=None):
                         subdir=out_dct["subdir"],
                     ),
                 )
-                wf.connect(nii_name, "out_file", ds, f'{out_dct["subdir"]}.@data')
+                if resource.endswith("_bold"):
+                    # Node to validate TR (and other scan parameters)
+                    validate_bold_header = pe.Node(
+                        Function(
+                            input_names=["input_bold", "RawSource_bold"],
+                            output_names=["output_bold"],
+                            function=validate_outputs,
+                            imports=[
+                                "from CPAC.pipeline.utils import find_pixdim4, update_pixdim4"
+                            ],
+                        ),
+                        name=f"validate_bold_header_{resource_idx}_{pipe_x}",
+                    )
+                    raw_source, raw_out = self.get_data("bold")
+                    wf.connect(
+                        [
+                            (nii_name, validate_bold_header, [(out, "input_bold")]),
+                            (
+                                raw_source,
+                                validate_bold_header,
+                                [(raw_out, "RawSource_bold")],
+                            ),
+                            (
+                                validate_bold_header,
+                                ds,
+                                [("output_bold", f'{out_dct["subdir"]}.@data')],
+                            ),
+                        ]
+                    )
+                else:
+                    wf.connect(nii_name, "out_file", ds, f'{out_dct["subdir"]}.@data')
+
                 wf.connect(write_json, "json_file", ds, f'{out_dct["subdir"]}.@json')
         outputs_logger.info(expected_outputs)
 

CPAC/pipeline/utils.py

@@ -17,13 +17,143 @@
 """C-PAC pipeline engine utilities."""
 
 from itertools import chain
+import os
+import subprocess
+
+import nibabel as nib
 
 from CPAC.func_preproc.func_motion import motion_estimate_filter
 from CPAC.utils.bids_utils import insert_entity
+from CPAC.utils.monitoring import IFLOGGER
 
 MOVEMENT_FILTER_KEYS = motion_estimate_filter.outputs
 
 
+def find_pixdim4(file_path):
+    """Find the pixdim4 value of a NIfTI file.
+
+    Parameters
+    ----------
+    file_path : str
+        Path to the NIfTI file.
+
+    Returns
+    -------
+    float
+        The pixdim4 value of the NIfTI file.
+
+    Raises
+    ------
+    FileNotFoundError
+        If the file does not exist.
+    nibabel.filebasedimages.ImageFileError
+        If there is an error loading the NIfTI file.
+    IndexError
+        If pixdim4 is not found in the header.
+    """
+    if not os.path.isfile(file_path):
+        error_message = f"File not found: {file_path}"
+        raise FileNotFoundError(file_path)
+
+    try:
+        nii = nib.load(file_path)
+        header = nii.header
+        pixdim = header.get_zooms()
+        return pixdim[3]
+    except nib.filebasedimages.ImageFileError as e:
+        error_message = f"Error loading the NIfTI file: {e}"
+        raise nib.filebasedimages.ImageFileError(error_message)
+    except IndexError as e:
+        error_message = f"pixdim4 not found in the header: {e}"
+        raise IndexError(error_message)
+
+
+def update_pixdim4(file_path, new_pixdim4):
+    """Update the pixdim4 value of a NIfTI file using 3drefit.
+
+    Parameters
+    ----------
+    file_path : str
+        Path to the NIfTI file.
+    new_pixdim4 : float
+        New pixdim4 value to update the NIfTI file with.
+
+    Raises
+    ------
+    FileNotFoundError
+        If the file does not exist.
+    subprocess.CalledProcessError
+        If there is an error running the subprocess.
+
+    Notes
+    -----
+    The pixdim4 value is the Repetition Time (TR) of the NIfTI file.
+
+    """
+    if not os.path.isfile(file_path):
+        error_message = f"File not found: {file_path}"
+        raise FileNotFoundError(error_message)
+
+    # Print the current pixdim4 value for verification
+    IFLOGGER.info(f"Updating {file_path} with new pixdim[4] value: {new_pixdim4}")
+
+    # Construct the command to update the pixdim4 value using 3drefit
+    command = ["3drefit", "-TR", str(new_pixdim4), file_path]
+
+    try:
+        subprocess.run(command, check=True)
+        IFLOGGER.info(f"Successfully updated TR to {new_pixdim4} seconds.")
+    except subprocess.CalledProcessError as e:
+        error_message = f"Error occurred while updating the file: {e}"
+        raise subprocess.CalledProcessError(error_message)
+
+
+def validate_outputs(input_bold, RawSource_bold):
+    """Match pixdim4/TR of the input_bold with RawSource_bold.
+
+    Parameters
+    ----------
+    input_bold : str
+        Path to the input BOLD file.
+    RawSource_bold : str
+        Path to the RawSource BOLD file.
+
+    Returns
+    -------
+    output_bold : str
+        Path to the output BOLD file.
+
+    Raises
+    ------
+    Exception
+        If there is an error in finding or updating pixdim4.
+    """
+    try:
+        output_bold = input_bold
+        output_pixdim4 = find_pixdim4(output_bold)
+        source_pixdim4 = find_pixdim4(RawSource_bold)
+
+        if output_pixdim4 != source_pixdim4:
+            IFLOGGER.info(
+                "TR mismatch detected between output_bold and RawSource_bold."
+            )
+            IFLOGGER.info(f"output_bold TR: {output_pixdim4} seconds")
+            IFLOGGER.info(f"RawSource_bold TR: {source_pixdim4} seconds")
+            IFLOGGER.info(
+                "Attempting to update the TR of output_bold to match RawSource_bold."
+            )
+            update_pixdim4(output_bold, source_pixdim4)
+        else:
+            IFLOGGER.debug("TR match detected between output_bold and RawSource_bold.")
+            IFLOGGER.debug(f"output_bold TR: {output_pixdim4} seconds")
+            IFLOGGER.debug(f"RawSource_bold TR: {source_pixdim4} seconds")
+        return output_bold
+    except Exception as e:
+        error_message = f"Error in validating outputs: {e}"
+        IFLOGGER.error(error_message)
+        return output_bold
+
+
 def name_fork(resource_idx, cfg, json_info, out_dct):
     """Create and insert entities for forkpoints.