✨ Feature/update pixdim4 for the final *bold outputs (#2169)\

.github/Dockerfiles/base-standard.Dockerfile

@@ -25,6 +25,10 @@ USER root
 # Installing FreeSurfer
 RUN apt-get update \
     && apt-get install --no-install-recommends -y bc \
+    && if [ ! -e /lib/x86_64-linux-gnu/libcrypt.so.2 ]; then \
+    # until we upgrade to Python >=3.11
+      ln -s /lib/x86_64-linux-gnu/libcrypt.so.1 /lib/x86_64-linux-gnu/libcrypt.so.2; \
+    fi \
     && yes | mamba install tcsh \
     && yes | mamba clean --all \
     && cp -l `which tcsh` /bin/tcsh \

CHANGELOG.md

@@ -22,6 +22,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Required positional parameter "wf" in input and output of `ingress_pipeconfig_paths` function, where a node to reorient templates is added to the `wf`.
 - Required positional parameter "orientation" to `resolve_resolution`.
 - Optional positional argument "cfg" to `create_lesion_preproc`.
+- Validation node to match the pixdim4 of CPAC processed bold outputs with the original raw bold sources.
 
 ### Changed
 

CPAC/pipeline/engine.py

@@ -36,11 +36,7 @@
 from CPAC.image_utils.statistical_transforms import z_score_standardize, \
     fisher_z_score_standardize
 from CPAC.pipeline.check_outputs import ExpectedOutputs
-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.utils.bids_utils import res_in_filename
 from CPAC.utils.datasource import (
@@ -63,6 +59,7 @@
 
 logger = getLogger('nipype.workflow')
 
+
 class ResourcePool:
     def __init__(self, rpool=None, name=None, cfg=None, pipe_list=None):
 
@@ -1235,6 +1232,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}_'
                                               f'{pipe_x}')
                 ds.inputs.parameterization = False
@@ -1251,8 +1249,33 @@ def gather_pipes(self, wf, cfg, all=False, add_incl=None, add_excl=None):
                     self.cfg, unique_id, resource_idx,
                     template_desc=id_string.inputs.template_desc,
                     atlas_id=atlas_id, 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)
@@ -1716,6 +1739,7 @@ def wrap_block(node_blocks, interface, wf, cfg, strat_pool, pipe_num, opt):
 
 def ingress_raw_anat_data(wf, rpool, cfg, data_paths, unique_id, part_id,
                           ses_id):
+
     if 'anat' not in data_paths:
         print('No anatomical data present.')
         return rpool
@@ -2169,8 +2193,8 @@ def ingress_pipeconfig_paths(wf, cfg, rpool, unique_id, creds_path=None):
     # ingress config file paths
     # TODO: may want to change the resource keys for each to include one level up in the YAML as well
 
-    import pandas as pd
     import pkg_resources as p
+    import pandas as pd
 
     template_csv = p.resource_filename('CPAC', 'resources/cpac_templates.csv')
     template_df = pd.read_csv(template_csv, keep_default_na=False)

CPAC/pipeline/utils.py

@@ -23,8 +23,130 @@
 from CPAC.func_preproc.func_motion import motion_estimate_filter
 from CPAC.utils.bids_utils import insert_entity
 
-IFLOGGER = logging.getLogger("nipype.interface")
 MOVEMENT_FILTER_KEYS = motion_estimate_filter.outputs
+import nibabel as nib
+import os
+import subprocess
+
+IFLOGGER = logging.getLogger('nipype.interface')
+
+
+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.
+    """
+    output_bold = input_bold
+    try:
+        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):