diff --git a/CHANGELOG.md b/CHANGELOG.md index 0bf7516ef7..8b9af1ce55 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -29,6 +29,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - New switch `mask_sbref` under `func_input_prep` in functional registration and set to default `on`. - New resource `desc-head_bold` as non skull-stripped bold from nodeblock `bold_masking`. - `censor_file_path` from `offending_timepoints_connector` in the `build_nuisance_regressor` node. +- `mri_robust_template` for longitudinal template generation. +- `max_iter` parameter for longitudinal template generation. ### Changed @@ -45,6 +47,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Input `desc-brain_bold` to `desc-preproc_bold` for `sbref` generation nodeblock `coregistration_prep_vol`. - Turned `generate_xcpqc_files` on for all preconfigurations except `blank`. - Introduced specific switch `restore_t1w_intensity` for `correct_restore_brain_intensity_abcd` nodeblock, enabling it by default only in `abcd-options` pre-config. +- Updated GitHub Actions to run automated integration and regression tests on HPC. +- Made `mri_robust_template` default implementation for longitudinal template generation. ### Fixed @@ -54,6 +58,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - A bug in which bandpass filters always assumed 1D regressor files have exactly 5 header rows. - Removed an erroneous connection to AFNI 3dTProject in nuisance denoising that would unnecessarily send a spike regressor as a censor. This would sometimes cause TRs to unnecessarily be dropped from the timeseries as if scrubbing were being performed. - Supplied missing `subject_id` for longitudinal workflow logger and make that field optional for the logger. +- Lingering calls to `cpac_outputs.csv` (was changed to `cpac_outputs.tsv` in v1.8.1). +- A bug in the `freesurfer_abcd_preproc` nodeblock where the `Template` image was incorrectly used as `reference` during the `inverse_warp` step. Replacing it with the subject-specific `T1w` image resolved the issue of the `desc-restoreBrain_T1w` being chipped off. +- A bug in `ideal_bandpass` where the frequency mask was incorrectly applied, which caused filter to fail in certain cases. ### Removed diff --git a/CPAC/alff/alff.py b/CPAC/alff/alff.py index f8bfc1a0b8..4f8aac4e57 100644 --- a/CPAC/alff/alff.py +++ b/CPAC/alff/alff.py @@ -1,5 +1,5 @@ # -*- coding: utf-8 -*- -# Copyright (C) 2012-2024 C-PAC Developers +# Copyright (C) 2012-2025 C-PAC Developers # This file is part of C-PAC. @@ -23,9 +23,8 @@ from CPAC.alff.utils import get_opt_string from CPAC.pipeline import nipype_pipeline_engine as pe from CPAC.pipeline.nodeblock import nodeblock -from CPAC.registration.registration import apply_transform +from CPAC.registration.utils import apply_transform from CPAC.utils.interfaces import Function -from CPAC.utils.utils import check_prov_for_regtool def create_alff(wf_name="alff_workflow"): @@ -320,10 +319,7 @@ def alff_falff(wf, cfg, strat_pool, pipe_num, opt=None): def alff_falff_space_template(wf, cfg, strat_pool, pipe_num, opt=None): outputs = {} if strat_pool.check_rpool("desc-denoisedNofilt_bold"): - xfm_prov = strat_pool.get_cpac_provenance( - "from-bold_to-template_mode-image_xfm" - ) - reg_tool = check_prov_for_regtool(xfm_prov) + reg_tool = strat_pool.reg_tool("from-bold_to-template_mode-image_xfm") num_cpus = cfg.pipeline_setup["system_config"]["max_cores_per_participant"] diff --git a/CPAC/anat_preproc/anat_preproc.py b/CPAC/anat_preproc/anat_preproc.py index 751fb499d3..67764e085b 100644 --- a/CPAC/anat_preproc/anat_preproc.py +++ b/CPAC/anat_preproc/anat_preproc.py @@ -34,7 +34,7 @@ wb_command, ) from CPAC.pipeline import nipype_pipeline_engine as pe -from CPAC.pipeline.nodeblock import nodeblock +from CPAC.pipeline.nodeblock import nodeblock, NODEBLOCK_RETURN from CPAC.utils.interfaces import Function from CPAC.utils.interfaces.fsl import Merge as fslMerge @@ -1227,15 +1227,7 @@ def freesurfer_fsl_brain_connector(wf, cfg, strat_pool, pipe_num, opt): wf.connect(node, out, convert_fs_T1_to_nifti, "in_file") # 3dresample -orient RPI -inset brainmask.nii.gz -prefix brain_fs.nii.gz - reorient_fs_brainmask = pe.Node( - interface=afni.Resample(), - name=f"reorient_fs_brainmask_{node_id}", - mem_gb=0, - mem_x=(0.0115, "in_file", "t"), - ) - reorient_fs_brainmask.inputs.orientation = cfg.pipeline_setup["desired_orientation"] - reorient_fs_brainmask.inputs.outputtype = "NIFTI_GZ" - + reorient_fs_brainmask = cfg.orientation_node(f"reorient_fs_brainmask_{node_id}") wf.connect( convert_fs_brainmask_to_nifti, "out_file", reorient_fs_brainmask, "in_file" ) @@ -1249,15 +1241,7 @@ def freesurfer_fsl_brain_connector(wf, cfg, strat_pool, pipe_num, opt): wf.connect(reorient_fs_brainmask, "out_file", binarize_fs_brain, "in_file") # 3dresample -orient RPI -inset T1.nii.gz -prefix head_fs.nii.gz - reorient_fs_T1 = pe.Node( - interface=afni.Resample(), - name=f"reorient_fs_T1_{node_id}", - mem_gb=0, - mem_x=(0.0115, "in_file", "t"), - ) - reorient_fs_T1.inputs.orientation = cfg.pipeline_setup["desired_orientation"] - reorient_fs_T1.inputs.outputtype = "NIFTI_GZ" - + reorient_fs_T1 = cfg.orientation_node(f"reorient_fs_T1_{node_id}") wf.connect(convert_fs_T1_to_nifti, "out_file", reorient_fs_T1, "in_file") # flirt -in head_fs.nii.gz -ref ${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz \ @@ -1447,22 +1431,14 @@ def mask_T2(wf_name="mask_T2"): inputs=["T1w"], outputs=["desc-preproc_T1w", "desc-reorient_T1w", "desc-head_T1w"], ) -def anatomical_init(wf, cfg, strat_pool, pipe_num, opt=None): +def anatomical_init(wf, cfg, strat_pool, pipe_num, opt=None) -> NODEBLOCK_RETURN: anat_deoblique = pe.Node(interface=afni.Refit(), name=f"anat_deoblique_{pipe_num}") anat_deoblique.inputs.deoblique = True node, out = strat_pool.get_data("T1w") wf.connect(node, out, anat_deoblique, "in_file") - anat_reorient = pe.Node( - interface=afni.Resample(), - name=f"anat_reorient_{pipe_num}", - mem_gb=0, - mem_x=(0.0115, "in_file", "t"), - ) - anat_reorient.inputs.orientation = cfg.pipeline_setup["desired_orientation"] - anat_reorient.inputs.outputtype = "NIFTI_GZ" - + anat_reorient = cfg.orientation_node(f"anat_reorient_{pipe_num}") wf.connect(anat_deoblique, "out_file", anat_reorient, "in_file") outputs = { @@ -2262,15 +2238,7 @@ def anatomical_init_T2(wf, cfg, strat_pool, pipe_num, opt=None): node, out = strat_pool.get_data("T2w") wf.connect(node, out, T2_deoblique, "in_file") - T2_reorient = pe.Node( - interface=afni.Resample(), - name=f"T2_reorient_{pipe_num}", - mem_gb=0, - mem_x=(0.0115, "in_file", "t"), - ) - T2_reorient.inputs.orientation = cfg.pipeline_setup["desired_orientation"] - T2_reorient.inputs.outputtype = "NIFTI_GZ" - + T2_reorient = cfg.orientation_node(f"T2_reorient_{pipe_num}") wf.connect(T2_deoblique, "out_file", T2_reorient, "in_file") outputs = { diff --git a/CPAC/anat_preproc/lesion_preproc.py b/CPAC/anat_preproc/lesion_preproc.py index 21628c97f0..dc872483bc 100644 --- a/CPAC/anat_preproc/lesion_preproc.py +++ b/CPAC/anat_preproc/lesion_preproc.py @@ -1,5 +1,5 @@ # -*- coding: utf-8 -*- -# Copyright (C) 2019-2023 C-PAC Developers +# Copyright (C) 2019-2025 C-PAC Developers # This file is part of C-PAC. @@ -20,6 +20,7 @@ from CPAC.pipeline import nipype_pipeline_engine as pe from CPAC.utils.interfaces import Function +from CPAC.utils.nifti_utils import orientation_node def inverse_lesion(lesion_path): @@ -126,18 +127,10 @@ def create_lesion_preproc(cfg=None, wf_name="lesion_preproc"): preproc.connect(lesion_deoblique, "out_file", outputnode, "refit") # Anatomical reorientation - lesion_reorient = pe.Node( - interface=afni.Resample(), - name="lesion_reorient", - mem_gb=0, - mem_x=(0.0115, "in_file", "t"), + node_name = "lesion_reorient" + lesion_reorient = ( + cfg.orientation_node(node_name) if cfg else orientation_node(node_name, "RPI") ) - - lesion_reorient.inputs.orientation = ( - cfg.pipeline_setup["desired_orientation"] if cfg else "RPI" - ) - lesion_reorient.inputs.outputtype = "NIFTI_GZ" - preproc.connect(lesion_deoblique, "out_file", lesion_reorient, "in_file") preproc.connect(lesion_reorient, "out_file", outputnode, "reorient") diff --git a/CPAC/func_preproc/func_motion.py b/CPAC/func_preproc/func_motion.py index bea7d2e29c..5a6794c1b2 100644 --- a/CPAC/func_preproc/func_motion.py +++ b/CPAC/func_preproc/func_motion.py @@ -1,4 +1,4 @@ -# Copyright (C) 2012-2024 C-PAC Developers +# Copyright (C) 2012-2025 C-PAC Developers # This file is part of C-PAC. @@ -34,7 +34,6 @@ from CPAC.pipeline.nodeblock import nodeblock from CPAC.pipeline.schema import valid_options from CPAC.utils.interfaces.function import Function -from CPAC.utils.utils import check_prov_for_motion_tool @nodeblock( @@ -68,8 +67,7 @@ ) def calc_motion_stats(wf, cfg, strat_pool, pipe_num, opt=None): """Calculate motion statistics for motion parameters.""" - motion_prov = strat_pool.get_cpac_provenance("desc-movementParameters_motion") - motion_correct_tool = check_prov_for_motion_tool(motion_prov) + motion_correct_tool = strat_pool.motion_tool("desc-movementParameters_motion") coordinate_transformation = [ "filtered-coordinate-transformation", "coordinate-transformation", diff --git a/CPAC/func_preproc/func_preproc.py b/CPAC/func_preproc/func_preproc.py index ff626765c4..b519ac3501 100644 --- a/CPAC/func_preproc/func_preproc.py +++ b/CPAC/func_preproc/func_preproc.py @@ -1,4 +1,4 @@ -# Copyright (C) 2012-2023 C-PAC Developers +# Copyright (C) 2012-2025 C-PAC Developers # This file is part of C-PAC. @@ -718,16 +718,7 @@ def func_reorient(wf, cfg, strat_pool, pipe_num, opt=None): node, out = strat_pool.get_data("bold") wf.connect(node, out, func_deoblique, "in_file") - func_reorient = pe.Node( - interface=afni_utils.Resample(), - name=f"func_reorient_{pipe_num}", - mem_gb=0, - mem_x=(0.0115, "in_file", "t"), - ) - - func_reorient.inputs.orientation = cfg.pipeline_setup["desired_orientation"] - func_reorient.inputs.outputtype = "NIFTI_GZ" - + func_reorient = cfg.orientation_node(f"func_reorient_{pipe_num}") wf.connect(func_deoblique, "out_file", func_reorient, "in_file") outputs = { @@ -1320,16 +1311,7 @@ def bold_mask_anatomical_refined(wf, cfg, strat_pool, pipe_num, opt=None): node, out = strat_pool.get_data("bold") wf.connect(node, out, func_deoblique, "in_file") - func_reorient = pe.Node( - interface=afni_utils.Resample(), - name=f"raw_func_reorient_{pipe_num}", - mem_gb=0, - mem_x=(0.0115, "in_file", "t"), - ) - - func_reorient.inputs.orientation = cfg.pipeline_setup["desired_orientation"] - func_reorient.inputs.outputtype = "NIFTI_GZ" - + func_reorient = cfg.orientation_node(f"raw_func_reorient_{pipe_num}") wf.connect(func_deoblique, "out_file", func_reorient, "in_file") wf.connect(func_reorient, "out_file", init_bold_mask, "inputspec.func") diff --git a/CPAC/longitudinal/__init__.py b/CPAC/longitudinal/__init__.py index 4545170f29..42bc68c8c5 100644 --- a/CPAC/longitudinal/__init__.py +++ b/CPAC/longitudinal/__init__.py @@ -18,6 +18,7 @@ from CPAC.utils.docs import DOCS_URL_PREFIX +assert isinstance(__doc__, str) __doc__ += f""" See {DOCS_URL_PREFIX}/user/longitudinal diff --git a/CPAC/longitudinal/preproc.py b/CPAC/longitudinal/preproc.py index 9fbe31c6b5..a884057f41 100644 --- a/CPAC/longitudinal/preproc.py +++ b/CPAC/longitudinal/preproc.py @@ -19,7 +19,9 @@ from collections import Counter from multiprocessing.dummy import Pool as ThreadPool +from multiprocessing.pool import Pool import os +from typing import Literal, Optional import numpy as np import nibabel as nib @@ -131,27 +133,23 @@ def norm_transformation(input_mat): def template_convergence( - mat_file, mat_type="matrix", convergence_threshold=np.finfo(np.float64).eps -): + mat_file: str, + mat_type: Literal["matrix", "ITK"] = "matrix", + convergence_threshold: float | np.float64 = np.finfo(np.float64).eps, +) -> bool: """Check that the deistance between matrices is smaller than the threshold. Calculate the distance between transformation matrix with a matrix of no transformation. Parameters ---------- - mat_file : str + mat_file path to an fsl flirt matrix - mat_type : str - 'matrix'(default), 'ITK' + mat_type The type of matrix used to represent the transformations - convergence_threshold : float - (numpy.finfo(np.float64).eps (default)) threshold for the convergence + convergence_threshold The threshold is how different from no transformation is the transformation matrix. - - Returns - ------- - bool """ if mat_type == "matrix": translation, oth_transform = read_mat(mat_file) @@ -346,51 +344,68 @@ def flirt_node(in_img, output_img, output_mat): return node_list +def check_convergence(mat_list, mat_type, convergence_threshold) -> bool: + """Test if every transformation matrix has reached the convergence threshold.""" + convergence_list = [ + template_convergence(mat, mat_type, convergence_threshold) for mat in mat_list + ] + return all(convergence_list) + + +@Function.sig_imports( + [ + "from multiprocessing.pool import Pool", + "from typing import Literal, Optional", + "from nipype.pipeline import engine as pe", + "from CPAC.longitudinal.preproc import check_convergence", + ] +) def template_creation_flirt( - input_brain_list, - input_skull_list, - init_reg=None, - avg_method="median", - dof=12, - interp="trilinear", - cost="corratio", - mat_type="matrix", - convergence_threshold=-1, - thread_pool=2, - unique_id_list=None, -): + input_brain_list: list[str], + input_skull_list: list[str], + init_reg: Optional[list[pe.Node]] = None, + avg_method: Literal["median", "mean", "std"] = "median", + dof: Literal[12, 9, 7, 6] = 12, + interp: Literal["trilinear", "nearestneighbour", "sinc", "spline"] = "trilinear", + cost: Literal[ + "corratio", "mutualinfo", "normmi", "normcorr", "leastsq", "labeldiff", "bbr" + ] = "corratio", + mat_type: Literal["matrix", "ITK"] = "matrix", + convergence_threshold: float | np.float64 = -1, + max_iter: int = 5, + thread_pool: int | Pool = 2, + unique_id_list: Optional[list[str]] = None, +) -> tuple[str, str, list[str], list[str], list[str]]: """Create a temporary template from a list of images. Parameters ---------- - input_brain_list : list of str + input_brain_list list of brain images paths - input_skull_list : list of str + input_skull_list list of skull images paths - init_reg : list of Node + init_reg (default None so no initial registration performed) the output of the function register_img_list with another reference Reuter et al. 2012 (NeuroImage) section "Improved template estimation" doi:10.1016/j.neuroimage.2012.02.084 uses a ramdomly selected image from the input dataset - avg_method : str - function names from numpy library such as 'median', 'mean', 'std' ... - dof : integer (int of long) - number of transform degrees of freedom (FLIRT) (12 by default) - interp : str - ('trilinear' (default) or 'nearestneighbour' or 'sinc' or 'spline') + avg_method + function names from numpy library + dof + number of transform degrees of freedom (FLIRT) + interp final interpolation method used in reslicing - cost : str - ('mutualinfo' or 'corratio' (default) or 'normcorr' or 'normmi' or - 'leastsq' or 'labeldiff' or 'bbr') + cost cost function - mat_type : str - 'matrix'(default), 'ITK' + mat_type The type of matrix used to represent the transformations - convergence_threshold : float + convergence_threshold (numpy.finfo(np.float64).eps (default)) threshold for the convergence The threshold is how different from no transformation is the transformation matrix. + max_iter + Maximum number of iterations if transformation does not converge thread_pool : int or multiprocessing.dummy.Pool (default 2) number of threads. You can also provide a Pool so the node will be added to it to be run. @@ -463,6 +478,8 @@ def template_creation_flirt( warp_list, ) + output_brain_list = list(input_brain_list) + output_skull_list = list(input_skull_list) # Chris: I added this part because it is mentioned in the paper but I actually never used it # You could run a first register_img_list() with a selected image as starting point and # give the output to this function @@ -471,18 +488,11 @@ def template_creation_flirt( output_brain_list = [node.inputs.out_file for node in init_reg] mat_list = [node.inputs.out_matrix_file for node in init_reg] warp_list = mat_list - # test if every transformation matrix has reached the convergence - convergence_list = [ - template_convergence(mat, mat_type, convergence_threshold) - for mat in mat_list - ] - converged = all(convergence_list) + converged = check_convergence(mat_list, mat_type, convergence_threshold) else: msg = "init_reg must be a list of FLIRT nipype nodes files" raise ValueError(msg) else: - output_brain_list = input_brain_list - output_skull_list = input_skull_list converged = False temporary_brain_template = os.path.join( @@ -496,7 +506,14 @@ def template_creation_flirt( and the loop stops when this temporary template is close enough (with a transformation distance smaller than the threshold) to all the images of the precedent iteration. """ - while not converged: + iterator = 1 + iteration = 0 + if max_iter == -1: + # make iteration < max_iter always True + iterator = 0 + iteration = -2 + while not converged and iteration < max_iter: + iteration += iterator temporary_brain_template, temporary_skull_template = create_temporary_template( input_brain_list=output_brain_list, input_skull_list=output_skull_list, @@ -551,13 +568,7 @@ def template_creation_flirt( warp_list[index] = warp_list_filenames[index] output_brain_list = [node.inputs.out_file for node in reg_list_node] - - # test if every transformation matrix has reached the convergence - convergence_list = [ - template_convergence(mat, mat_type, convergence_threshold) - for mat in mat_list - ] - converged = all(convergence_list) + converged = check_convergence(mat_list, mat_type, convergence_threshold) if isinstance(thread_pool, int): pool.close() @@ -609,7 +620,7 @@ def subject_specific_template( "from collections import Counter", "from multiprocessing.dummy import Pool as ThreadPool", "from nipype.interfaces.fsl import ConvertXFM", - "from CPAC.longitudinal_pipeline.longitudinal_preproc import (" + "from CPAC.longitudinal.preproc import (" " create_temporary_template," " register_img_list," " template_convergence" @@ -628,6 +639,7 @@ def subject_specific_template( "cost", "mat_type", "convergence_threshold", + "max_iter", "thread_pool", "unique_id_list", ], diff --git a/CPAC/longitudinal/robust_template.py b/CPAC/longitudinal/robust_template.py new file mode 100644 index 0000000000..16a52e8e12 --- /dev/null +++ b/CPAC/longitudinal/robust_template.py @@ -0,0 +1,154 @@ +# -*- coding: utf-8 -*- +# Copyright (C) 2024-2025 C-PAC Developers + +# This file is part of C-PAC. + +# C-PAC is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or (at your +# option) any later version. + +# C-PAC is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +# License for more details. + +# You should have received a copy of the GNU Lesser General Public +# License along with C-PAC. If not, see . +"""Create longitudinal template using ``mri_robust_template``.""" + +import os +from typing import cast, Literal + +from nipype.interfaces.base import ( + File, + InputMultiPath, + isdefined, + OutputMultiPath, + traits, +) +from nipype.interfaces.freesurfer import longitudinal +from nipype.interfaces.freesurfer.preprocess import MRIConvert + +from CPAC.pipeline import nipype_pipeline_engine as pe +from CPAC.utils.configuration import Configuration + + +class RobustTemplateInputSpec(longitudinal.RobustTemplateInputSpec): # noqa: D101 + affine = traits.Bool(default_value=False, desc="compute 12 DOF registration") + mapmov = traits.Either( + InputMultiPath(File(exists=False)), + traits.Bool, + argstr="--mapmov %s", + desc="output images: map and resample each input to template", + ) + maxit = traits.Int( + argstr="--maxit %d", + mandatory=False, + desc="iterate max # times (if #tp>2 default 6, else 5 for 2tp reg.)", + ) + + +class RobustTemplateOutputSpec(longitudinal.RobustTemplateOutputSpec): # noqa: D101 + mapmov = OutputMultiPath( + File(), + desc="each input mapped and resampled to longitudinal template", + ) + + +class RobustTemplate(longitudinal.RobustTemplate): # noqa: D101 + # STATEMENT OF CHANGES: + # This class is derived from sources licensed under the Apache-2.0 terms, + # and this class has been changed. + + # CHANGES: + # * Added handling for `affine`, `mapmov` and `maxit`. + # * Renamed transform outputs. + + # ORIGINAL WORK'S ATTRIBUTION NOTICE: + # Copyright (c) 2009-2016, Nipype developers + + # Licensed under the Apache License, Version 2.0 (the "License"); + # you may not use this file except in compliance with the License. + # You may obtain a copy of the License at + + # http://www.apache.org/licenses/LICENSE-2.0 + + # Unless required by applicable law or agreed to in writing, software + # distributed under the License is distributed on an "AS IS" BASIS, + # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + # See the License for the specific language governing permissions and + # limitations under the License. + + # Prior to release 0.12, Nipype was licensed under a BSD license. + + # Modifications copyright (C) 2024 C-PAC Developers + input_spec = RobustTemplateInputSpec + output_spec = RobustTemplateOutputSpec + + def _format_arg(self, name, spec, value): + if name == "average_metric": + # return enumeration value + return spec.argstr % {"mean": 0, "median": 1}[value] + if name in ("mapmov", "transform_outputs", "scaled_intensity_outputs"): + value = self._list_outputs()[name] + return super()._format_arg(name, spec, value) + + def _list_outputs(self): + """:py:meth:`~nipype.interfaces.freesurfer.RobustTemplate._list_outputs` + `mapmov`.""" + outputs = self.output_spec().get() + outputs["out_file"] = os.path.abspath(self.inputs.out_file) + n_files = len(self.inputs.in_files) + fmt = "{}{:02d}.{}" if n_files > 9 else "{}{:d}.{}" # noqa: PLR2004 + for key, prefix, ext in [ + ("transform_outputs", "space-longitudinal", "lta"), + ("scaled_intensity_outputs", "is", "txt"), + ("mapmov", "space-longitudinal", "nii.gz"), + ]: + if isdefined(getattr(self.inputs, key)): + fnames = getattr(self.inputs, key) + if fnames is True: + fnames = [fmt.format(prefix, i + 1, ext) for i in range(n_files)] + outputs[key] = [os.path.abspath(x) for x in fnames] + return outputs + + +def mri_robust_template( + name: str, cfg: Configuration, num_sessions: int +) -> pe.Workflow: + """Return a subworkflow to run `mri_robust_template` with common options. + + Converts transform files to FSL format. + """ + wf = pe.Workflow(name=name) + node = pe.Node( + RobustTemplate( + affine=cfg["longitudinal_template_generation", "dof"] == 12, # noqa: PLR2004 + average_metric=cfg["longitudinal_template_generation", "average_method"], + auto_detect_sensitivity=True, + mapmov=True, + out_file=f"{name}.mgz", + transform_outputs=True, + ), + name="mri_robust_template", + ) + max_iter = cast( + int | Literal["default"], cfg["longitudinal_template_generation", "max_iter"] + ) + if isinstance(max_iter, int): + node.set_input("maxit", max_iter) + + nifti_template = pe.Node(MRIConvert(out_type="niigz"), name="NIfTI-template") + wf.connect(node, "out_file", nifti_template, "in_file") + reorient_template = cfg.orientation_node("reorient_longitudinal_template", pe.Node) + wf.connect(nifti_template, "out_file", reorient_template, "in_file") + + nifti_outputs = pe.MapNode(MRIConvert(), name="NIfTI-mapmov", iterfield=["in_file"]) + wf.connect(node, "mapmov", nifti_outputs, "in_file") + reorient_outputs = cfg.orientation_node("reorient_longitudinal_session", pe.MapNode) + wf.connect(nifti_outputs, "out_file", reorient_outputs, "in_file") + reorient_outputs.set_input( + "out_file", [f"space-longitudinal{i + 1}.nii.gz" for i in range(num_sessions)] + ) + + return wf diff --git a/CPAC/longitudinal/wf/__init__.py b/CPAC/longitudinal/wf/__init__.py new file mode 100644 index 0000000000..990d1541f5 --- /dev/null +++ b/CPAC/longitudinal/wf/__init__.py @@ -0,0 +1,25 @@ +# Copyright (C) 2024 C-PAC Developers + +# This file is part of C-PAC. + +# C-PAC is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or (at your +# option) any later version. + +# C-PAC is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +# License for more details. + +# You should have received a copy of the GNU Lesser General Public +# License along with C-PAC. If not, see . +"""Workflows for longitudinal preprocessing.""" + +from CPAC.utils.docs import DOCS_URL_PREFIX + +assert isinstance(__doc__, str) +__doc__ += f""" + +See {DOCS_URL_PREFIX}/user/longitudinal +""" # noqa: A001 diff --git a/CPAC/longitudinal/wf/anat.py b/CPAC/longitudinal/wf/anat.py index ab34b03a69..baf2b8ce3b 100644 --- a/CPAC/longitudinal/wf/anat.py +++ b/CPAC/longitudinal/wf/anat.py @@ -1,5 +1,5 @@ # -*- coding: utf-8 -*- -# Copyright (C) 2020-2024 C-PAC Developers +# Copyright (C) 2020-2025 C-PAC Developers # This file is part of C-PAC. @@ -15,27 +15,35 @@ # You should have received a copy of the GNU Lesser General Public # License along with C-PAC. If not, see . -import os +"""Longitudinal workflows for anatomical data.""" +from typing import cast, Optional + +from networkx.classes.digraph import DiGraph +from nipype import config as nipype_config from nipype.interfaces import fsl -from indi_aws import aws_utils +from nipype.interfaces.utility import Merge from CPAC.longitudinal.preproc import subject_specific_template +from CPAC.longitudinal.robust_template import mri_robust_template +from CPAC.longitudinal.wf.utils import ( + check_creds_path, + cross_graph_connections, + cross_graph_identity, + get_output_from_graph, + select_session_node, +) from CPAC.pipeline import nipype_pipeline_engine as pe from CPAC.pipeline.cpac_pipeline import ( build_anat_preproc_stack, - build_segmentation_stack, build_T1w_registration_stack, connect_pipeline, initialize_nipype_wf, ) -from CPAC.pipeline.engine import ingress_output_dir, initiate_rpool -from CPAC.pipeline.nodeblock import nodeblock -from CPAC.registration.registration import apply_transform -from CPAC.utils.configuration.configuration import Configuration -from CPAC.utils.interfaces.datasink import DataSink -from CPAC.utils.interfaces.function import Function -from CPAC.utils.utils import check_prov_for_regtool +from CPAC.pipeline.engine import ingress_output_dir, initiate_rpool, ResourcePool +from CPAC.pipeline.nodeblock import nodeblock, NODEBLOCK_RETURN +from CPAC.registration.utils import apply_transform +from CPAC.utils.configuration import Configuration @nodeblock( @@ -45,7 +53,10 @@ inputs=["desc-brain_T1w"], outputs=["space-T1w_desc-brain_mask"], ) -def mask_T1w_longitudinal_template(wf, cfg, strat_pool, pipe_num, opt=None): +def mask_T1w_longitudinal_template( + wf: pe.Workflow, cfg, strat_pool, pipe_num, opt=None +) -> NODEBLOCK_RETURN: + """Create a native-space brain mask for longitudinal template generation.""" brain_mask = pe.Node( interface=fsl.maths.MathsCommand(), name=f"longitudinal_anatomical_brain_mask_{pipe_num}", @@ -57,177 +68,51 @@ def mask_T1w_longitudinal_template(wf, cfg, strat_pool, pipe_num, opt=None): outputs = {"space-T1w_desc-brain_mask": (brain_mask, "out_file")} - return (wf, outputs) + return wf, outputs -def create_datasink( - datasink_name, - config, - subject_id, - session_id="", - strat_name="", - map_node_iterfield=None, -) -> pe.Node | pe.MapNode: - """ - Parameters - ---------- - datasink_name - config - subject_id - session_id - strat_name - map_node_iterfield - """ - encrypt_data = config.pipeline_setup["Amazon-AWS"]["s3_encryption"] - - # TODO Enforce value with schema validation - # Extract credentials path for output if it exists - try: - # Get path to creds file - creds_path = "" - if config.pipeline_setup["Amazon-AWS"]["aws_output_bucket_credentials"]: - creds_path = str( - config.pipeline_setup["Amazon-AWS"]["aws_output_bucket_credentials"] - ) - creds_path = os.path.abspath(creds_path) - - if ( - config.pipeline_setup["output_directory"]["path"] - .lower() - .startswith("s3://") - ): - # Test for s3 write access - s3_write_access = aws_utils.test_bucket_access( - creds_path, config.pipeline_setup["output_directory"]["path"] - ) - - if not s3_write_access: - msg = "Not able to write to bucket!" - raise Exception(msg) - - except Exception as e: - if ( - config.pipeline_setup["output_directory"]["path"] - .lower() - .startswith("s3://") - ): - err_msg = ( - "There was an error processing credentials or " - "accessing the S3 bucket. Check and try again.\n" - "Error: %s" % e - ) - raise Exception(err_msg) - - if map_node_iterfield is not None: - ds = pe.MapNode( - DataSink(infields=map_node_iterfield), - name=f"sinker_{datasink_name}", - iterfield=map_node_iterfield, - ) - else: - ds = pe.Node(DataSink(), name=f"sinker_{datasink_name}") - - ds.inputs.base_directory = config.pipeline_setup["output_directory"]["path"] - ds.inputs.creds_path = creds_path - ds.inputs.encrypt_bucket_keys = encrypt_data - ds.inputs.container = os.path.join( - "pipeline_%s_%s" % (config.pipeline_setup["pipeline_name"], strat_name), - subject_id, - session_id, - ) - return ds - - -def connect_anat_preproc_inputs( - strat, anat_preproc, strat_name, strat_nodes_list_list, workflow -): - """ - Parameters - ---------- - strat : Strategy - the strategy object you want to fork - anat_preproc : Workflow - the anat_preproc workflow node to be connected and added to the resource pool - strat_name : str - name of the strategy - strat_nodes_list_list : list - a list of strat_nodes_list - workflow : Workflow - main longitudinal workflow - - Returns - ------- - new_strat : Strategy - the fork of strat with the resource pool updated - strat_nodes_list_list : list - a list of strat_nodes_list - """ - new_strat = strat.fork() - - tmp_node, out_key = new_strat["anatomical"] - workflow.connect(tmp_node, out_key, anat_preproc, "inputspec.anat") - - tmp_node, out_key = new_strat["template_cmass"] - workflow.connect(tmp_node, out_key, anat_preproc, "inputspec.template_cmass") - - new_strat.append_name(anat_preproc.name) - - new_strat.update_resource_pool( - { - "anatomical_brain": (anat_preproc, "outputspec.brain"), - "anatomical_skull_leaf": (anat_preproc, "outputspec.reorient"), - "anatomical_brain_mask": (anat_preproc, "outputspec.brain_mask"), - } - ) - - try: - strat_nodes_list_list[strat_name].append(new_strat) - except KeyError: - strat_nodes_list_list[strat_name] = [new_strat] - - return new_strat, strat_nodes_list_list - - -def pick_map(file_list, index, file_type): +def pick_map( + file_list: list[list[str]] | list[str], index: str, file_type: str +) -> Optional[str]: + """Choose a file from a list of files.""" if isinstance(file_list, list): - if len(file_list) == 1: + if len(file_list) == 1 and isinstance(file_list[0], list): file_list = file_list[0] for file_name in file_list: + assert isinstance(file_name, str) if file_name.endswith(f"{file_type}_{index}.nii.gz"): return file_name return None -def select_session(session, output_brains, warps): - brain_path = None - warp_path = None - for brain_path in output_brains: - if f"{session}_" in brain_path: - break - for warp_path in warps: - if f"{session}_" in warp_path: - break - return (brain_path, warp_path) - - @nodeblock( name="mask_longitudinal_T1w_brain", config=["longitudinal_template_generation"], switch=["run"], - inputs=["space-longitudinal_desc-brain_T1w"], - outputs=["space-longitudinal_desc-brain_mask"], + inputs=["longitudinal-template_space-longitudinal_desc-brain_T1w"], + outputs=["longitudinal-template_space-longitudinal_desc-brain_mask"], ) -def mask_longitudinal_T1w_brain(wf, cfg, strat_pool, pipe_num, opt=None): +def mask_longitudinal_T1w_brain( + wf, cfg, strat_pool, pipe_num, opt=None +) -> NODEBLOCK_RETURN: + """Create brain mask for longitudinal T1w image.""" brain_mask = pe.Node( interface=fsl.maths.MathsCommand(), name=f"longitudinal_T1w_brain_mask_{pipe_num}", ) brain_mask.inputs.args = "-bin" - node, out = strat_pool.get_data("space-longitudinal_desc-brain_T1w") + node, out = strat_pool.get_data( + "longitudinal-template_space-longitudinal_desc-brain_T1w" + ) wf.connect(node, out, brain_mask, "in_file") - outputs = {"space-longitudinal_desc-brain_mask": (brain_mask, "out_file")} + outputs = { + "longitudinal-template_space-longitudinal_desc-brain_mask": ( + brain_mask, + "out_file", + ) + } return (wf, outputs) @@ -238,24 +123,24 @@ def mask_longitudinal_T1w_brain(wf, cfg, strat_pool, pipe_num, opt=None): switch=["run"], inputs=[ ( - "space-longitudinal_desc-brain_T1w", + "longitudinal-template_space-longitudinal_desc-brain_T1w", "from-longitudinal_to-template_mode-image_xfm", - ) + ), + "T1w-brain-template", ], outputs=["space-template_desc-brain_T1w"], ) -def warp_longitudinal_T1w_to_template(wf, cfg, strat_pool, pipe_num, opt=None): - xfm_prov = strat_pool.get_cpac_provenance( - "from-longitudinal_to-template_mode-image_xfm" - ) - reg_tool = check_prov_for_regtool(xfm_prov) - +def warp_longitudinal_T1w_to_template( + wf, cfg, strat_pool, pipe_num, opt=None +) -> NODEBLOCK_RETURN: + """Transform longitudinal T1w images to template space.""" + reg_tool = strat_pool.reg_tool("from-longitudinal_to-template_mode-image_xfm") num_cpus = cfg.pipeline_setup["system_config"]["max_cores_per_participant"] num_ants_cores = cfg.pipeline_setup["system_config"]["num_ants_threads"] apply_xfm = apply_transform( - f"warp_longitudinal_to_T1template_{pipe_num}", + f"warp_longitudinal_to_template_{pipe_num}", reg_tool, time_series=False, num_cpus=num_cpus, @@ -271,18 +156,25 @@ def warp_longitudinal_T1w_to_template(wf, cfg, strat_pool, pipe_num, opt=None): "anatomical_registration" ]["registration"]["FSL-FNIRT"]["interpolation"] - node, out = strat_pool.get_data("space-longitudinal_desc-brain_T1w") + node, out = strat_pool.get_data( + "longitudinal-template_space-longitudinal_desc-brain_T1w" + ) wf.connect(node, out, apply_xfm, "inputspec.input_image") - node, out = strat_pool.get_data("T1w_brain_template") + node, out = strat_pool.get_data("T1w-brain-template") wf.connect(node, out, apply_xfm, "inputspec.reference") node, out = strat_pool.get_data("from-longitudinal_to-template_mode-image_xfm") wf.connect(node, out, apply_xfm, "inputspec.transform") - outputs = {"space-template_desc-brain_T1w": (apply_xfm, "outputspec.output_image")} + outputs = { + "space-template_desc-brain_T1w": ( + apply_xfm, + "outputspec.output_image", + ) + } - return (wf, outputs) + return wf, outputs @nodeblock( @@ -291,7 +183,11 @@ def warp_longitudinal_T1w_to_template(wf, cfg, strat_pool, pipe_num, opt=None): switch=["run"], inputs=[ ( - "from-longitudinal_to-T1w_mode-image_desc-linear_xfm", + "space-longitudinal_desc-brain_T1w", + [ + "from-longitudinal_to-T1w_mode-image_desc-linear_xfm", + "from-T1w_to-longitudinal_mode-image_desc-linear_xfm", + ], "space-longitudinal_label-CSF_mask", "space-longitudinal_label-GM_mask", "space-longitudinal_label-WM_mask", @@ -301,32 +197,71 @@ def warp_longitudinal_T1w_to_template(wf, cfg, strat_pool, pipe_num, opt=None): "space-longitudinal_label-CSF_probseg", "space-longitudinal_label-GM_probseg", "space-longitudinal_label-WM_probseg", - ) - ], - outputs=[ - "label-CSF_mask", - "label-GM_mask", - "label-WM_mask", - "label-CSF_desc-preproc_mask", - "label-GM_desc-preproc_mask", - "label-WM_desc-preproc_mask", - "label-CSF_probseg", - "label-GM_probseg", - "label-WM_probseg", + ), + "T1w-brain-template", ], + outputs={ + "from-longitudinal_to-T1w_mode-image_desc-linear_xfm": {}, + "from-longitudinal_to-T1w_mode-image_desc-linear_warp": {}, + "label-CSF_mask": {}, + "label-GM_mask": {}, + "label-WM_mask": {}, + "label-CSF_desc-preproc_mask": {}, + "label-GM_desc-preproc_mask": {}, + "label-WM_desc-preproc_mask": {}, + "label-CSF_probseg": {}, + "label-GM_probseg": {}, + "label-WM_probseg": {}, + }, ) -def warp_longitudinal_seg_to_T1w(wf, cfg, strat_pool, pipe_num, opt=None): - xfm_prov = strat_pool.get_cpac_provenance( - "from-longitudinal_to-T1w_mode-image_desc-linear_xfm" - ) - reg_tool = check_prov_for_regtool(xfm_prov) +def warp_longitudinal_seg_to_T1w( + wf: pe.Workflow, + cfg: Configuration, + strat_pool: ResourcePool, + pipe_num: int, + opt: Optional[str] = None, +) -> NODEBLOCK_RETURN: + """Transform anatomical segmentation from longitudinal template to T1w space.""" + outputs = {} + if strat_pool.check_rpool("from-longitudinal_to-T1w_mode-image_desc-linear_xfm"): + reg_tool = strat_pool.reg_tool( + "from-longitudinal_to-T1w_mode-image_desc-linear_xfm" + ) + xfm: tuple[pe.Node, str] = strat_pool.get_data( + "from-longitudinal_to-T1w_mode-image_desc-linear_xfm" + ) + else: + reg_tool = strat_pool.reg_tool( + "from-T1w_to-longitudinal_mode-image_desc-linear_xfm" + ) + # create inverse xfm if we don't have it + invt = pe.Node(interface=fsl.ConvertXFM(), name=f"convert_xfm_{pipe_num}") + invt.inputs.invert_xfm = True + wf.connect( + *strat_pool.get_data("from-T1w_to-longitudinal_mode-image_desc-linear_xfm"), + invt, + "in_file", + ) + xfm = (invt, "out_file") + outputs["from-longitudinal_to-T1w_mode-image_desc-linear_xfm"] = xfm + if reg_tool == "fsl": + warp = pe.Node( + fsl.ConvertWarp(relwarp=True, out_relwarp=True), + name=f"convert_warp_{pipe_num}", + ) + wf.connect(*xfm, warp, "postmat") + wf.connect( + *strat_pool.get_data("space-longitudinal_desc-brain_T1w"), + warp, + "reference", + ) + outputs["from-longitudinal_to-T1w_mode-image_desc-linear_warp"] = ( + warp, + "out_file", + ) num_cpus = cfg.pipeline_setup["system_config"]["max_cores_per_participant"] - num_ants_cores = cfg.pipeline_setup["system_config"]["num_ants_threads"] - - outputs = {} - labels = [ "CSF_mask", "CSF_desc-preproc_mask", @@ -338,7 +273,6 @@ def warp_longitudinal_seg_to_T1w(wf, cfg, strat_pool, pipe_num, opt=None): "WM_desc-preproc_mask", "WM_probseg", ] - for label in labels: apply_xfm = apply_transform( f"warp_longitudinal_seg_to_T1w_{label}_{pipe_num}", @@ -360,73 +294,66 @@ def warp_longitudinal_seg_to_T1w(wf, cfg, strat_pool, pipe_num, opt=None): node, out = strat_pool.get_data("space-longitudinal_desc-brain_T1w") wf.connect(node, out, apply_xfm, "inputspec.input_image") - node, out = strat_pool.get_data("T1w_brain_template") + node, out = strat_pool.get_data("T1w-brain-template") wf.connect(node, out, apply_xfm, "inputspec.reference") - node, out = strat_pool.get_data("from-longitudinal_to-template_mode-image_xfm") - wf.connect(node, out, apply_xfm, "inputspec.transform") - + wf.connect(warp, "out_file", apply_xfm, "inputspec.transform") outputs[f"label-{label}"] = (apply_xfm, "outputspec.output_image") - return (wf, outputs) + return wf, outputs def anat_longitudinal_wf( - subject_id: str, sub_list: list[dict], config: Configuration + subject_id: str, sub_list: list[dict], config: Configuration, dry_run: bool = False ) -> None: """ - Create and run anatomical longitudinal workflow(s). + Create and run longitudinal workflows for anatomical data. Parameters ---------- subject_id the id of the subject sub_list - this is a list of sessions for one subject and each session if the same dictionary as the one given to - prep_workflow + a list of sessions for one subject config - a configuration object containing the information of the pipeline config. (Same as for prep_workflow) + a Configuration object containing the information for the participant pipeline + dry_run + build graph without running? """ + nipype_config.update_config( + { + "execution": { + "crashfile_format": "txt", + "stop_on_first_crash": config[ + "pipeline_setup", "system_config", "fail_fast" + ], + } + } + ) config["subject_id"] = subject_id - session_id_list: list[list] = [] + session_id_list: list[str] = [] """List of lists for every strategy""" session_wfs = {} - cpac_dirs = [] - out_dir = config.pipeline_setup["output_directory"]["path"] - - orig_pipe_name = config.pipeline_setup["pipeline_name"] + orig_pipe_name: str = config.pipeline_setup["pipeline_name"] - # Loop over the sessions to create the input for the longitudinal - # algorithm + strats_dct: dict[str, list[tuple[pe.Node, str] | str]] = { + "desc-brain_T1w": [], + "desc-head_T1w": [], + } for session in sub_list: - unique_id = session["unique_id"] + # Loop over the sessions to create the input for the longitudinal algorithm + unique_id: str = session["unique_id"] session_id_list.append(unique_id) + input_creds_path = check_creds_path(session.get("creds_path"), subject_id) - try: - creds_path = session["creds_path"] - if creds_path and "none" not in creds_path.lower(): - if os.path.exists(creds_path): - input_creds_path = os.path.abspath(creds_path) - else: - err_msg = ( - 'Credentials path: "%s" for subject "%s" ' - 'session "%s" was not found. Check this path ' - "and try again." % (creds_path, subject_id, unique_id) - ) - raise Exception(err_msg) - else: - input_creds_path = None - except KeyError: - input_creds_path = None - - workflow = initialize_nipype_wf( + workflow: pe.Workflow = initialize_nipype_wf( config, - sub_list[0], - # just grab the first one for the name + subject_id, + unique_id, name="anat_longitudinal_pre-preproc", ) - + rpool: ResourcePool workflow, rpool = initiate_rpool(workflow, config, session) pipeline_blocks = build_anat_preproc_stack(rpool, config) workflow = connect_pipeline(workflow, config, rpool, pipeline_blocks) @@ -435,231 +362,273 @@ def anat_longitudinal_wf( rpool.gather_pipes(workflow, config) - workflow.run() + for key in strats_dct.keys(): + strats_dct[key].append(cast(tuple[pe.Node, str], rpool.get_data(key))) + if not dry_run: + workflow_graph: DiGraph = workflow.run() + for key in strats_dct.keys(): # get the outputs from run-nodes + for index, data in enumerate(list(strats_dct[key])): + if isinstance(data, tuple): + strats_dct[key][index] = get_output_from_graph( + workflow, workflow_graph, *data + ) - cpac_dir = os.path.join( - out_dir, f"pipeline_{orig_pipe_name}", f"{subject_id}_{unique_id}" - ) - cpac_dirs.append(os.path.join(cpac_dir, "anat")) - - # Now we have all the anat_preproc set up for every session - # loop over the different anat preproc strategies - strats_brain_dct = {} - strats_head_dct = {} - for cpac_dir in cpac_dirs: - if os.path.isdir(cpac_dir): - for filename in os.listdir(cpac_dir): - if "T1w.nii" in filename: - for tag in filename.split("_"): - if "desc-" in tag and "brain" in tag: - if tag not in strats_brain_dct: - strats_brain_dct[tag] = [] - strats_brain_dct[tag].append( - os.path.join(cpac_dir, filename) - ) - if tag not in strats_head_dct: - strats_head_dct[tag] = [] - head_file = filename.replace(tag, "desc-reorient") - strats_head_dct[tag].append( - os.path.join(cpac_dir, head_file) - ) - - for strat in strats_brain_dct.keys(): - wf = initialize_nipype_wf( - config, - sub_list[0], - # just grab the first one for the name - name=f"template_node_{strat}", + wf = initialize_nipype_wf( + config, + subject_id, + name="template_node_brain", + ) + + num_sessions = len(strats_dct["desc-brain_T1w"]) + merge_brains = pe.Node(Merge(num_sessions), name="merge_brains") + merge_skulls = pe.Node(Merge(num_sessions), name="merge_skulls") + wf.add_nodes([merge_brains, merge_skulls]) + for i in list(range(0, num_sessions)): + wf._connect_node_or_path_for_merge( + merge_brains, strats_dct, "desc-brain_T1w", i ) + wf._connect_node_or_path_for_merge(merge_skulls, strats_dct, "desc-head_T1w", i) - config.pipeline_setup["pipeline_name"] = f"longitudinal_{orig_pipe_name}" + long_id = f"{subject_id}_desc-brain_T1w" + wf, rpool = initiate_rpool(wf, config, part_id=subject_id) - template_node_name = f"longitudinal_anat_template_{strat}" + match config["longitudinal_template_generation", "using"]: + case "C-PAC legacy": + brain_output = "brain_template" + head_output = "skull_template" - # This node will generate the longitudinal template (the functions are - # in longitudinal_preproc) - # Later other algorithms could be added to calculate it, like the - # multivariate template from ANTS - # It would just require to change it here. - template_node = subject_specific_template(workflow_name=template_node_name) + # This node will generate the longitudinal template (the functions are + # in longitudinal_preproc) + # Later other algorithms could be added to calculate it, like the + # multivariate template from ANTS + # It would just require to change it here. - template_node.inputs.set( - avg_method=config.longitudinal_template_generation["average_method"], - dof=config.longitudinal_template_generation["dof"], - interp=config.longitudinal_template_generation["interp"], - cost=config.longitudinal_template_generation["cost"], - convergence_threshold=config.longitudinal_template_generation[ - "convergence_threshold" - ], - thread_pool=config.longitudinal_template_generation["thread_pool"], - unique_id_list=list(session_wfs.keys()), - ) + # multiple variable names here for compatibility with other options later in this function + brain_template_node = wholehead_template_node = template_node = ( + subject_specific_template(workflow_name="longitudinal_anat_template") + ) - template_node.inputs.input_brain_list = strats_brain_dct[strat] - template_node.inputs.input_skull_list = strats_head_dct[strat] + template_node.inputs.set( + avg_method=config.longitudinal_template_generation["average_method"], + dof=config.longitudinal_template_generation["dof"], + interp=config.longitudinal_template_generation["legacy-specific"][ + "interp" + ], + cost=config.longitudinal_template_generation["legacy-specific"]["cost"], + convergence_threshold=config.longitudinal_template_generation[ + "legacy-specific" + ]["convergence_threshold"], + thread_pool=config.longitudinal_template_generation["legacy-specific"][ + "thread_pool" + ], + unique_id_list=list(session_wfs.keys()), + ) - long_id = f"longitudinal_{subject_id}_strat-{strat}" + wf.connect(merge_brains, "out", brain_template_node, "input_brain_list") + wf.connect(merge_skulls, "out", wholehead_template_node, "input_skull_list") - wf, rpool = initiate_rpool(wf, config, part_id=long_id) + case "mri_robust_template": + brain_output = head_output = "reorient_longitudinal_template.out_file" + brain_template_node = mri_robust_template( + f"mri_robust_template_brain_{subject_id}", config, len(sub_list) + ) + wholehead_template_node = mri_robust_template( + f"mri_robust_template_head_{subject_id}", config, len(sub_list) + ) + wf.connect( + merge_brains, "out", brain_template_node, "mri_robust_template.in_files" + ) + wf.connect( + merge_skulls, + "out", + wholehead_template_node, + "mri_robust_template.in_files", + ) - rpool.set_data( - "space-longitudinal_desc-brain_T1w", - template_node, - "brain_template", - {}, - "", - template_node_name, - ) + case _: + msg = ": ".join( + [ + "Invalid 'using' value for longitudinal template generation", + str(config["longitudinal_template_generation", "using"]), + ] + ) + raise ValueError(msg) + + rpool.set_data( + "longitudinal-template_space-longitudinal_desc-brain_T1w", + brain_template_node, + brain_output, + {}, + "", + brain_template_node.name, + ) + rpool.set_data( + "longitudinal-template_space-longitudinal_desc-head_T1w", + wholehead_template_node, + head_output, + {}, + "", + wholehead_template_node.name, + ) - rpool.set_data( - "space-longitudinal_desc-brain_T1w-template", - template_node, - "brain_template", - {}, - "", - template_node_name, - ) + pipeline_blocks = [mask_longitudinal_T1w_brain] + pipeline_blocks = build_T1w_registration_stack( + rpool, config, pipeline_blocks, space="longitudinal" + ) - rpool.set_data( - "space-longitudinal_desc-reorient_T1w", - template_node, - "skull_template", - {}, - "", - template_node_name, - ) + cross_pool_keys = [ + # "from-longitudinal_to-template_mode-image_xfm", + # "from-template_to-longitudinal_mode-image_desc-linear_xfm", + "longitudinal-template_space-longitudinal_desc-brain_mask", + "longitudinal-template_space-longitudinal_desc-brain_T1w", + "longitudinal-template_space-longitudinal_desc-head_T1w", + ] + rpool.gather_pipes(wf, config) + wf = connect_pipeline(wf, config, rpool, pipeline_blocks) + + wf_graph: DiGraph | pe.Workflow = ( + cast(DiGraph, wf.run()) if not dry_run else cast(pe.Workflow, wf) + ) + + # now, just write out a copy of the above to each session + config.pipeline_setup["pipeline_name"] = orig_pipe_name + longitudinal_rpool = rpool + for i, session in enumerate(sub_list): + unique_id = session["unique_id"] + input_creds_path = check_creds_path(session.get("creds_path"), subject_id) + + ses_wf = initialize_nipype_wf(config, subject_id, unique_id) + + ses_wf, rpool = initiate_rpool(ses_wf, config, session) + + if "derivatives_dir" in session: + ses_wf, rpool = ingress_output_dir( + ses_wf, + config, + rpool, + long_id, + data_paths=session, + part_id=subject_id, + ses_id=unique_id, + creds_path=input_creds_path, + ) + + select_sess = select_session_node(unique_id) + + match config["longitudinal_template_generation", "using"]: + case "C-PAC legacy": + cross_graph_connections( + wf, + wf_graph, + ses_wf, + merge_brains, + brain_template_node, + "out", + "input_brain_list", + ) + cross_graph_connections( + wf, + wf_graph, + ses_wf, + merge_skulls, + brain_template_node, + "out", + "input_skull_list", + ) + cross_graph_connections( + wf, + wf_graph, + ses_wf, + brain_template_node, + select_sess, + "output_brain_list", + "outputs", + ) + + case "mri_robust_template": + assert isinstance(brain_template_node, pe.Workflow) + assert isinstance(wholehead_template_node, pe.Workflow) + index = i + 1 + head_select_sess = select_session_node(unique_id, "wholehead") + select_sess.set_input("session", f"space-longitudinal{index}") + head_select_sess.set_input("session", f"space-longitudinal{index}") + input_name = "outputs" + output_name = "reorient_longitudinal_session_.out_file" + cross_graph_connections( + wf, + wf_graph, + ses_wf, + brain_template_node, + select_sess, + output_name, + input_name, + ) + cross_graph_connections( + wf, + wf_graph, + ses_wf, + wholehead_template_node, + head_select_sess, + output_name, + input_name, + ) + rpool.set_data( + "space-longitudinal_desc-head_T1w", + head_select_sess, + "path", + {}, + "", + head_select_sess.name, + ) rpool.set_data( - "space-longitudinal_desc-reorient_T1w-template", - template_node, - "skull_template", + "space-longitudinal_desc-brain_T1w", + select_sess, + "path", {}, "", - template_node_name, + select_sess.name, ) - pipeline_blocks = [mask_longitudinal_T1w_brain] - - pipeline_blocks = build_T1w_registration_stack(rpool, config, pipeline_blocks) - - pipeline_blocks = build_segmentation_stack(rpool, config, pipeline_blocks) - - wf = connect_pipeline(wf, config, rpool, pipeline_blocks) - + config.pipeline_setup["pipeline_name"] = orig_pipe_name excl = [ + # "from-T1w_to-longitudinal_mode-image_desc-linear_xfm", "space-longitudinal_desc-brain_T1w", - "space-longitudinal_desc-reorient_T1w", - "space-longitudinal_desc-brain_mask", + "space-longitudinal_desc-head_T1w", + # "space-template_desc-brain_T1w", + # "space-T1w_desc-brain_mask", ] - rpool.gather_pipes(wf, config, add_excl=excl) - - # this is going to run multiple times! - # once for every strategy! - wf.run() - - # now, just write out a copy of the above to each session - config.pipeline_setup["pipeline_name"] = orig_pipe_name - for session in sub_list: - unique_id = session["unique_id"] - + rpool.gather_pipes(ses_wf, config, add_excl=excl) + for key in cross_pool_keys: + node, out = longitudinal_rpool.get_data(key) + if isinstance(wf_graph, DiGraph): + assert isinstance(out, str) + node = cross_graph_identity(wf, wf_graph, node, out) try: - creds_path = session["creds_path"] - if creds_path and "none" not in creds_path.lower(): - if os.path.exists(creds_path): - input_creds_path = os.path.abspath(creds_path) - else: - err_msg = ( - 'Credentials path: "%s" for subject "%s" ' - 'session "%s" was not found. Check this path ' - "and try again." % (creds_path, subject_id, unique_id) - ) - raise Exception(err_msg) - else: - input_creds_path = None - except KeyError: - input_creds_path = None - - wf = initialize_nipype_wf(config, sub_list[0]) - - wf, rpool = initiate_rpool(wf, config, session) - - config.pipeline_setup["pipeline_name"] = f"longitudinal_{orig_pipe_name}" - rpool = ingress_output_dir( - config, rpool, long_id, creds_path=input_creds_path - ) - - select_node_name = f"select_{unique_id}" - select_sess = pe.Node( - Function( - input_names=["session", "output_brains", "warps"], - output_names=["brain_path", "warp_path"], - function=select_session, - ), - name=select_node_name, - ) - select_sess.inputs.session = unique_id - - wf.connect(template_node, "output_brain_list", select_sess, "output_brains") - wf.connect(template_node, "warp_list", select_sess, "warps") - - rpool.set_data( - "space-longitudinal_desc-brain_T1w", - select_sess, - "brain_path", - {}, - "", - select_node_name, - ) - + json_info: dict = longitudinal_rpool.get_json( + key, next(iter(longitudinal_rpool.rpool[key].keys())) + ) + except (AttributeError, KeyError, StopIteration): + json_info = {} rpool.set_data( - "from-T1w_to-longitudinal_mode-image_desc-linear_xfm", - select_sess, - "warp_path", - {}, + key, + node, + out, + json_info, "", - select_node_name, + f"longitudinal_{subject_id}", ) - config.pipeline_setup["pipeline_name"] = orig_pipe_name - excl = ["space-template_desc-brain_T1w", "space-T1w_desc-brain_mask"] - - rpool.gather_pipes(wf, config, add_excl=excl) - wf.run() - - # begin single-session stuff again - for session in sub_list: - unique_id = session["unique_id"] - - try: - creds_path = session["creds_path"] - if creds_path and "none" not in creds_path.lower(): - if os.path.exists(creds_path): - input_creds_path = os.path.abspath(creds_path) - else: - err_msg = ( - 'Credentials path: "%s" for subject "%s" ' - 'session "%s" was not found. Check this path ' - "and try again." % (creds_path, subject_id, unique_id) - ) - raise Exception(err_msg) - else: - input_creds_path = None - except KeyError: - input_creds_path = None - - wf = initialize_nipype_wf(config, sub_list[0]) - - wf, rpool = initiate_rpool(wf, config, session) - - pipeline_blocks = [ - warp_longitudinal_T1w_to_template, - warp_longitudinal_seg_to_T1w, - ] - - wf = connect_pipeline(wf, config, rpool, pipeline_blocks) + # pipeline_blocks = build_segmentation_stack( + # rpool, + # config, + # [warp_longitudinal_T1w_to_template] # , warp_longitudinal_seg_to_T1w], + # ) - rpool.gather_pipes(wf, config) + ses_wf = connect_pipeline(ses_wf, config, rpool, pipeline_blocks) + rpool.gather_pipes(ses_wf, config) # this is going to run multiple times! # once for every strategy! - wf.run() + if not dry_run: # check select_sess + ses_wf.run() diff --git a/CPAC/longitudinal/wf/func.py b/CPAC/longitudinal/wf/func.py index 00847073f1..3a7ca0669e 100644 --- a/CPAC/longitudinal/wf/func.py +++ b/CPAC/longitudinal/wf/func.py @@ -22,6 +22,7 @@ import nipype.interfaces.io as nio from CPAC.longitudinal.preproc import subject_specific_template +from CPAC.longitudinal.wf.utils import check_creds_path from CPAC.pipeline import nipype_pipeline_engine as pe from CPAC.registration import ( create_fsl_flirt_linear_reg, @@ -82,22 +83,7 @@ def func_preproc_longitudinal_wf(subject_id, sub_list, config): unique_id = sub_dict["unique_id"] session_id_list.append(unique_id) - try: - creds_path = sub_dict["creds_path"] - if creds_path and "none" not in creds_path.lower(): - if os.path.exists(creds_path): - input_creds_path = os.path.abspath(creds_path) - else: - err_msg = ( - 'Credentials path: "%s" for subject "%s" was not ' - "found. Check this path and try again." - % (creds_path, subject_id) - ) - raise Exception(err_msg) - else: - input_creds_path = None - except KeyError: - input_creds_path = None + input_creds_path = check_creds_path(sub_dict.get("creds_path"), subject_id) strat = Strategy() strat_list = [strat] diff --git a/CPAC/longitudinal/wf/utils.py b/CPAC/longitudinal/wf/utils.py new file mode 100644 index 0000000000..f1bdd18677 --- /dev/null +++ b/CPAC/longitudinal/wf/utils.py @@ -0,0 +1,188 @@ +# -*- coding: utf-8 -*- +# Copyright (C) 2020-2024 C-PAC Developers + +# This file is part of C-PAC. + +# C-PAC is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or (at your +# option) any later version. + +# C-PAC is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +# License for more details. + +# You should have received a copy of the GNU Lesser General Public +# License along with C-PAC. If not, see . +"""Utilities for longitudinal workflows.""" + +from pathlib import Path +from typing import Any, cast, Optional + +from networkx.classes.digraph import DiGraph +from nipype.interfaces.utility import IdentityInterface + +from CPAC.pipeline import nipype_pipeline_engine as pe +from CPAC.pipeline.utils import get_edges_with_node +from CPAC.utils.interfaces.function import Function + +_TEMPLATE_PATTERN = r"((sym)?template|longitudinal)" +LONGITUDINAL_TEMPLATE_PATTERN = f"from-{_TEMPLATE_PATTERN}_to-{_TEMPLATE_PATTERN}" + + +def check_creds_path(creds_path: Optional[str], subject_id: str) -> Optional[str]: + """Check credentials path.""" + if creds_path and "none" not in creds_path.lower(): + _creds_path = Path(creds_path) + if _creds_path.exists(): + return str(_creds_path.absolute()) + err_msg = ( + 'Credentials path: "%s" for subject "%s" was not ' + "found. Check this path and try again." % (creds_path, subject_id) + ) + raise FileNotFoundError(err_msg) + return None + + +@Function.sig_imports( + [ + "from networkx.classes.digraph import DiGraph", + "from CPAC.pipeline import nipype_pipeline_engine as pe", + "from CPAC.longitudinal.wf.utils import get_output_from_graph", + ] +) +def cross_graph_connections( + wf: pe.Workflow, + wf1: DiGraph | pe.Workflow, + wf2: pe.Workflow, + node1: pe.Node | pe.Workflow, + node2: pe.Node | pe.Workflow, + output_name: str, + input_name: str, +) -> None: + """Make cross-graph connections appropriate to dry-run status. + + Parameters + ---------- + wf + The graph that already ran + + wf1 + The results of the graph that already ran + + wf2 + The graph that runs second + + node1 + The node from ``wf1`` + + node2 + The node from ``wf2`` + + output_name + The output name from ``node1`` + + input_name + The input name from ``node2`` + """ + if isinstance(wf1, pe.Workflow): # dry run + if isinstance(node1, pe.Workflow): + sub_node_name, output_name = output_name.rsplit(".", 1) + node1 = cast(pe.Node, node1.get_node(sub_node_name)) + wf2.connect(node1, output_name, node2, input_name) + else: + setattr( + node2.inputs, input_name, get_output_from_graph(wf, wf1, node1, output_name) + ) + + +def select_session(session: str, outputs: list[str]) -> str: + """Select output brain image and warp for given session.""" + try: + return next(iter(path for path in outputs if session in path)) + except StopIteration as stop_iteration: + msg = f"{session} not found in {outputs}.\n" + raise FileExistsError(msg) from stop_iteration + + +def select_session_node(unique_id: str, suffix: str = "") -> pe.Node: + """Create a Node to select a single subject's output image and transform. + + Note + ---- + FSL is the only currenlty implemented registration tool for longitudinal template + generation, so it's hardcoded into the name of this node for + feeding :py:meth:`~CPAC.utils.utils.check_prov_for_regtool`. + """ + if suffix: + suffix = f"_{suffix.lstrip('_')}" + select_sess = pe.Node( + Function( + input_names=["session", "outputs"], + output_names=["path"], + function=select_session, + ), + name=f"longitudinalSelect{suffix.title()}_{unique_id}", + ) + select_sess.set_input("session", f"{unique_id}_") + return select_sess + + +def cross_graph_identity( + wf: pe.Workflow, graph: DiGraph, node: pe.Node | pe.Workflow, output_name: str +) -> pe.Node: + """Sever connection to ``node_1``'s workflow while maintaining value.""" + identity_interface = pe.Node( + IdentityInterface([output_name]), + name="_".join( + [str(getattr(node, "fullname", getattr(node, "name", ""))), output_name] + ).replace(".", "_"), + ) + identity_interface.set_input( + output_name, get_output_from_graph(wf, graph, node, output_name) + ) + return identity_interface + + +def get_output_from_graph( + wf: pe.Workflow, graph: DiGraph, node: pe.Node | pe.Workflow, output_name: str +) -> Any: + """Get an output from a graph that has been run.""" + nodename = str(node.fullname) + if isinstance(node, pe.Workflow): + sub_node_name, output_name = output_name.rsplit(".", 1) + nodename = f"{nodename}.{sub_node_name}" + edges = get_edges_with_node(node, output_name) + for edge in reversed(edges): + try: + return get_output_from_graph( + wf, + graph, + edge[0], + next( + iter( + connection + for connection in edge[2]["connect"] + if connection[1] == output_name + ) + )[0], + ) + except StopIteration: + continue + try: + output = getattr( + next( + iter( + _node + for _node in graph + if _node.fullname.endswith(nodename) + or _node.fullname.endswith(f"{nodename}_") + ) + ).result.outputs, + output_name, + ) + except StopIteration as stop_iteration: + msg = f"{nodename} not found in completed workflow." + raise FileNotFoundError(msg) from stop_iteration + return output diff --git a/CPAC/nuisance/bandpass.py b/CPAC/nuisance/bandpass.py index 451d4a5b9e..626fbc32ce 100644 --- a/CPAC/nuisance/bandpass.py +++ b/CPAC/nuisance/bandpass.py @@ -1,3 +1,21 @@ +# Copyright (C) 2019-2025 C-PAC Developers + +# This file is part of C-PAC. + +# C-PAC is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or (at your +# option) any later version. + +# C-PAC is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +# License for more details. + +# You should have received a copy of the GNU Lesser General Public +# License along with C-PAC. If not, see . +"""Bandpass filtering utilities.""" + import os from pathlib import Path @@ -6,43 +24,59 @@ import nibabel as nib from scipy.fftpack import fft, ifft +from CPAC.utils.monitoring import IFLOGGER + def ideal_bandpass(data, sample_period, bandpass_freqs): - # Derived from YAN Chao-Gan 120504 based on REST. + """ + Apply ideal bandpass filtering to a 1D time series data using FFT. Derived from YAN Chao-Gan 120504 based on REST. + + Parameters + ---------- + data : NDArray + 1D time series data to be filtered. + sample_period : float + Length of sampling period in seconds. + bandpass_freqs : tuple + Tuple containing the bandpass frequencies (LowCutoff, HighCutoff). + + Returns + ------- + NDArray + Filtered time series data. + + """ sample_freq = 1.0 / sample_period sample_length = data.shape[0] + nyquist_freq = sample_freq / 2.0 - data_p = np.zeros(int(2 ** np.ceil(np.log2(sample_length)))) + # Length of zero-padded data for efficient FFT + N = int(2 ** np.ceil(np.log2(len(data)))) + data_p = np.zeros(N) data_p[:sample_length] = data LowCutoff, HighCutoff = bandpass_freqs if LowCutoff is None: # No lower cutoff (low-pass filter) low_cutoff_i = 0 - elif LowCutoff > sample_freq / 2.0: + elif LowCutoff > nyquist_freq: # Cutoff beyond fs/2 (all-stop filter) - low_cutoff_i = int(data_p.shape[0] / 2) + low_cutoff_i = int(N / 2) else: - low_cutoff_i = np.ceil(LowCutoff * data_p.shape[0] * sample_period).astype( - "int" - ) + low_cutoff_i = np.ceil(LowCutoff * N * sample_period).astype("int") - if HighCutoff > sample_freq / 2.0 or HighCutoff is None: + if HighCutoff is None or HighCutoff > nyquist_freq: # Cutoff beyond fs/2 or unspecified (become a highpass filter) - high_cutoff_i = int(data_p.shape[0] / 2) + high_cutoff_i = int(N / 2) else: - high_cutoff_i = np.fix(HighCutoff * data_p.shape[0] * sample_period).astype( - "int" - ) + high_cutoff_i = np.fix(HighCutoff * N * sample_period).astype("int") freq_mask = np.zeros_like(data_p, dtype="bool") freq_mask[low_cutoff_i : high_cutoff_i + 1] = True - freq_mask[data_p.shape[0] - high_cutoff_i : data_p.shape[0] + 1 - low_cutoff_i] = ( - True - ) + freq_mask[N - high_cutoff_i : N + 1 - low_cutoff_i] = True f_data = fft(data_p) - f_data[freq_mask is not True] = 0.0 + f_data[~freq_mask] = 0.0 return np.real_if_close(ifft(f_data)[:sample_length]) @@ -63,7 +97,7 @@ def read_1D(one_D: Path | str) -> tuple[list[str], NDArray]: def bandpass_voxels(realigned_file, regressor_file, bandpass_freqs, sample_period=None): - """Performs ideal bandpass filtering on each voxel time-series. + """Perform ideal bandpass filtering on each voxel time-series. Parameters ---------- @@ -91,7 +125,9 @@ def bandpass_voxels(realigned_file, regressor_file, bandpass_freqs, sample_perio hdr = nii.header sample_period = float(hdr.get_zooms()[3]) # Sketchy check to convert TRs in millisecond units - if sample_period > 20.0: + if sample_period > 20.0: # noqa: PLR2004 + message = f"Sample period ({sample_period}) is very large. Assuming milliseconds and converting to seconds." + IFLOGGER.warning(message) sample_period /= 1000.0 Y_bp = np.zeros_like(Y) diff --git a/CPAC/nuisance/nuisance.py b/CPAC/nuisance/nuisance.py index ce4c1298da..40dc2b1585 100644 --- a/CPAC/nuisance/nuisance.py +++ b/CPAC/nuisance/nuisance.py @@ -1,4 +1,4 @@ -# Copyright (C) 2012-2024 C-PAC Developers +# Copyright (C) 2012-2025 C-PAC Developers # This file is part of C-PAC. @@ -50,7 +50,6 @@ from CPAC.utils.interfaces.function import Function from CPAC.utils.interfaces.pc import PC from CPAC.utils.monitoring import IFLOGGER -from CPAC.utils.utils import check_prov_for_regtool from .bandpass import afni_1dBandpass, bandpass_voxels @@ -2011,8 +2010,7 @@ def filtering_bold_and_regressors( outputs=["desc-preproc_bold", "desc-cleaned_bold"], ) def ICA_AROMA_FSLreg(wf, cfg, strat_pool, pipe_num, opt=None): - xfm_prov = strat_pool.get_cpac_provenance("from-T1w_to-template_mode-image_xfm") - reg_tool = check_prov_for_regtool(xfm_prov) + reg_tool = strat_pool.reg_tool("from-T1w_to-template_mode-image_xfm") if reg_tool != "fsl": return (wf, None) @@ -2058,8 +2056,7 @@ def ICA_AROMA_FSLreg(wf, cfg, strat_pool, pipe_num, opt=None): outputs=["desc-preproc_bold", "desc-cleaned_bold"], ) def ICA_AROMA_ANTsreg(wf, cfg, strat_pool, pipe_num, opt=None): - xfm_prov = strat_pool.get_cpac_provenance("from-bold_to-template_mode-image_xfm") - reg_tool = check_prov_for_regtool(xfm_prov) + reg_tool = strat_pool.reg_tool("from-bold_to-template_mode-image_xfm") if reg_tool != "ants": return (wf, None) @@ -2129,8 +2126,7 @@ def ICA_AROMA_ANTsreg(wf, cfg, strat_pool, pipe_num, opt=None): outputs=["desc-preproc_bold", "desc-cleaned_bold"], ) def ICA_AROMA_FSLEPIreg(wf, cfg, strat_pool, pipe_num, opt=None): - xfm_prov = strat_pool.get_cpac_provenance("from-bold_to-EPItemplate_mode-image_xfm") - reg_tool = check_prov_for_regtool(xfm_prov) + reg_tool = strat_pool.reg_tool("from-bold_to-EPItemplate_mode-image_xfm") if reg_tool != "fsl": return (wf, None) @@ -2182,8 +2178,7 @@ def ICA_AROMA_FSLEPIreg(wf, cfg, strat_pool, pipe_num, opt=None): outputs=["desc-preproc_bold", "desc-cleaned_bold"], ) def ICA_AROMA_ANTsEPIreg(wf, cfg, strat_pool, pipe_num, opt=None): - xfm_prov = strat_pool.get_cpac_provenance("from-bold_to-EPItemplate_mode-image_xfm") - reg_tool = check_prov_for_regtool(xfm_prov) + reg_tool = strat_pool.reg_tool("from-bold_to-EPItemplate_mode-image_xfm") if reg_tool != "ants": return (wf, None) @@ -2513,15 +2508,13 @@ def nuisance_regressors_generation( if space == "T1w": prefixes[0] = "" if strat_pool.check_rpool("from-template_to-T1w_mode-image_desc-linear_xfm"): - xfm_prov = strat_pool.get_cpac_provenance( + reg_tool = strat_pool.reg_tool( "from-template_to-T1w_mode-image_desc-linear_xfm" ) - reg_tool = check_prov_for_regtool(xfm_prov) elif space == "bold": - xfm_prov = strat_pool.get_cpac_provenance( + reg_tool = strat_pool.reg_tool( "from-EPItemplate_to-bold_mode-image_desc-linear_xfm" ) - reg_tool = check_prov_for_regtool(xfm_prov) if reg_tool is not None: use_ants = reg_tool == "ants" else: diff --git a/CPAC/nuisance/tests/test_bandpass.py b/CPAC/nuisance/tests/test_bandpass.py index 452b55d3c7..364d4605d5 100644 --- a/CPAC/nuisance/tests/test_bandpass.py +++ b/CPAC/nuisance/tests/test_bandpass.py @@ -1,4 +1,4 @@ -# Copyright (C) 2022 - 2024 C-PAC Developers +# Copyright (C) 2022-2025 C-PAC Developers # This file is part of C-PAC. @@ -19,11 +19,14 @@ from importlib.abc import Traversable from importlib.resources import files from pathlib import Path +from typing import Optional +import numpy as np from numpy.typing import NDArray import pytest +from scipy.fft import fft -from CPAC.nuisance.bandpass import read_1D +from CPAC.nuisance.bandpass import ideal_bandpass, read_1D RAW_ONE_D: Traversable = files("CPAC").joinpath("nuisance/tests/regressors.1D") @@ -46,3 +49,66 @@ def test_read_1D(start_line: int, tmp_path: Path) -> None: assert data.shape == (10, 29) # all header lines should be captured assert len(header) == 5 - start_line + + +@pytest.mark.parametrize("sample_period", [1.0, 0.1]) +@pytest.mark.parametrize( + "lowcut, highcut, in_freq, out_freq", + [ + (0.005, 0.05, 0.01, 0.2), + (0.01, 0.1, 0.02, 0.15), + (0.02, 0.08, 0.04, 0.12), + (None, 0.1, 0.02, 0.15), + (0.2, None, 0.22, 0.1), + ], +) +def test_ideal_bandpass_with_various_cutoffs( + lowcut: Optional[float], + highcut: Optional[float], + in_freq: float, + out_freq: float, + sample_period: float, +) -> None: + """Test the ideal bandpass filter with various cutoff frequencies.""" + t = np.arange(512) * sample_period + signal = np.sin(2 * np.pi * in_freq * t) + np.sin(2 * np.pi * out_freq * t) + + filtered = ideal_bandpass(signal, sample_period, (lowcut, highcut)) + + freqs = np.fft.fftfreq(len(signal), d=sample_period) + orig_fft = np.abs(fft(signal)) + filt_fft = np.abs(fft(filtered)) + + idx_in = np.argmin(np.abs(freqs - in_freq)) + idx_out = np.argmin(np.abs(freqs - out_freq)) + + assert filt_fft[idx_in] > 0.5 * orig_fft[idx_in] + assert filt_fft[idx_out] < 0.1 * orig_fft[idx_out] + + +@pytest.mark.parametrize("sample_period", [1.0, 0.1]) +def test_ideal_bandpass_cutoffs_clamped_to_nyquist(sample_period): + """Test that ideal_bandpass clamps cutoffs to Nyquist frequency.""" + N = 512 + t = np.arange(N) * sample_period + nyquist = 0.5 / sample_period + + freq_below = nyquist * 0.95 + freq_above = nyquist * 1.05 + + signal = np.sin(2 * np.pi * freq_below * t) + np.sin(2 * np.pi * freq_above * t) + + lowcut = nyquist + 0.01 + highcut = nyquist + 0.1 + + filtered = ideal_bandpass(signal, sample_period, (lowcut, highcut)) + + freqs = np.fft.fftfreq(N, d=sample_period) + filt_fft = np.abs(fft(filtered)) + + idx_below = np.argmin(np.abs(freqs - freq_below)) + idx_above = np.argmin(np.abs(freqs - freq_above)) + + acceptable_threshold = 1e-3 # threshold for numerical stability + assert filt_fft[idx_below] < acceptable_threshold + assert filt_fft[idx_above] < acceptable_threshold diff --git a/CPAC/pipeline/cpac_pipeline.py b/CPAC/pipeline/cpac_pipeline.py index 1b64b286a8..b3a1c7f1bb 100644 --- a/CPAC/pipeline/cpac_pipeline.py +++ b/CPAC/pipeline/cpac_pipeline.py @@ -1,4 +1,4 @@ -# Copyright (C) 2012-2024 C-PAC Developers +# Copyright (C) 2012-2025 C-PAC Developers # This file is part of C-PAC. @@ -25,10 +25,10 @@ import sys import time from time import strftime +from typing import Literal, Optional import yaml import nipype -from nipype import config, logging from flowdump import save_workflow_json, WorkflowJSONMeta from indi_aws import aws_utils, fetch_creds @@ -130,7 +130,7 @@ # pylint: disable=wrong-import-order from CPAC.pipeline import nipype_pipeline_engine as pe from CPAC.pipeline.check_outputs import check_outputs -from CPAC.pipeline.engine import initiate_rpool, NodeBlock +from CPAC.pipeline.engine import initiate_rpool, NodeBlock, ResourcePool from CPAC.pipeline.nipype_pipeline_engine.plugins import ( LegacyMultiProcPlugin, MultiProcPlugin, @@ -163,7 +163,6 @@ warp_deriv_mask_to_EPItemplate, warp_deriv_mask_to_T1template, warp_sbref_to_T1template, - warp_T1mask_to_template, warp_timeseries_to_EPItemplate, warp_timeseries_to_T1template, warp_timeseries_to_T1template_abcd, @@ -171,7 +170,7 @@ warp_timeseries_to_T1template_deriv, warp_tissuemask_to_EPItemplate, warp_tissuemask_to_T1template, - warp_wholeheadT1_to_template, + warp_to_template, ) from CPAC.reho.reho import reho, reho_space_template from CPAC.sca.sca import dual_regression, multiple_regression, SCA_AVG @@ -199,6 +198,7 @@ from CPAC.utils.monitoring import ( FMLOGGER, getLogger, + init_loggers, log_nodes_cb, log_nodes_initial, LOGTAIL, @@ -222,11 +222,11 @@ def run_workflow( - sub_dict, - c, - run, - pipeline_timing_info=None, - p_name=None, + sub_dict: dict, + c: Configuration, + run: bool, + pipeline_timing_info: Optional[list] = None, + p_name: Optional[str] = None, plugin="MultiProc", plugin_args=None, test_config=False, @@ -257,8 +257,6 @@ def run_workflow( 0 for success 1 for general failure """ - from CPAC.utils.datasource import bidsier_prefix - if plugin is not None and not isinstance(plugin, str): msg = ( 'CPAC.pipeline.cpac_pipeline.run_workflow requires a ' @@ -273,37 +271,7 @@ def run_workflow( subject_id, p_name, log_dir = set_subject(sub_dict, c) c["subject_id"] = subject_id - - set_up_logger( - f"{subject_id}_expectedOutputs", - filename=f'{bidsier_prefix(c["subject_id"])}_' 'expectedOutputs.yml', - level="info", - log_dir=log_dir, - mock=True, - overwrite_existing=True, - ) - if c.pipeline_setup["Debugging"]["verbose"]: - set_up_logger("CPAC.engine", level="debug", log_dir=log_dir, mock=True) - - config.update_config( - { - "logging": { - "log_directory": log_dir, - "log_to_file": bool( - getattr(c.pipeline_setup["log_directory"], "run_logging", True) - ), - }, - "execution": { - "crashfile_format": "txt", - "resource_monitor_frequency": 0.2, - "stop_on_first_crash": c[ - "pipeline_setup", "system_config", "fail_fast" - ], - }, - } - ) - config.enable_resource_monitor() - logging.update_logging(config) + init_loggers(subject_id, c, log_dir, mock=True) # Start timing here pipeline_start_time = time.time() @@ -556,6 +524,7 @@ def run_workflow( workflow_result = None exitcode = 0 + cb_log_filename = os.path.join(log_dir, "callback.log") try: subject_info["resource_pool"] = [] @@ -567,8 +536,6 @@ def run_workflow( subject_info["status"] = "Running" # Create callback logger - cb_log_filename = os.path.join(log_dir, "callback.log") - try: if not os.path.exists(os.path.dirname(cb_log_filename)): os.makedirs(os.path.dirname(cb_log_filename)) @@ -599,8 +566,9 @@ def run_workflow( plugin = MultiProcPlugin(plugin_args) try: - # Actually run the pipeline now, for the current subject - workflow_result = workflow.run(plugin=plugin, plugin_args=plugin_args) + if run: + # Actually run the pipeline now, for the current subject + workflow_result = workflow.run(plugin=plugin, plugin_args=plugin_args) except UnicodeDecodeError: msg = ( "C-PAC migrated from Python 2 to Python 3 in v1.6.2 (see " @@ -815,7 +783,7 @@ def run_workflow( run_start=pipeline_start_datetime, run_finish=strftime("%Y-%m-%d %H:%M:%S"), output_check=check_outputs( - c.pipeline_setup["output_directory"]["path"], + c["pipeline_setup"]["output_directory"]["path"], log_dir, c.pipeline_setup["pipeline_name"], c["subject_id"], @@ -857,14 +825,19 @@ def remove_workdir(wdpath: str) -> None: FMLOGGER.warning("Could not remove working directory %s", wdpath) -def initialize_nipype_wf(cfg, sub_data_dct, name=""): +def initialize_nipype_wf( + cfg: Configuration, + subject: str, + session: Optional[str] = None, + name: Optional[str] = None, +) -> pe.Workflow: """Initialize a new nipype workflow.""" - if name: - name = f"_{name}" + name = f"_{name}" if name else "" - workflow_name = ( - f'cpac{name}_{sub_data_dct["subject_id"]}_{sub_data_dct["unique_id"]}' - ) + identifier = subject + if session: + identifier = "_".join([identifier, session]) + workflow_name = f"cpac{name}_{identifier}" wf = pe.Workflow(name=workflow_name) wf.base_dir = cfg.pipeline_setup["working_directory"]["path"] wf.config["execution"] = { @@ -1062,25 +1035,33 @@ def build_anat_preproc_stack(rpool, cfg, pipeline_blocks=None): return pipeline_blocks -def build_T1w_registration_stack(rpool, cfg, pipeline_blocks=None): +def build_T1w_registration_stack( + rpool: ResourcePool, + cfg: Configuration, + pipeline_blocks: Optional[list] = None, + space: Literal["longitudinal", "T1w"] = "T1w", +): """Build the T1w registration pipeline blocks.""" if not pipeline_blocks: pipeline_blocks = [] reg_blocks = [] - if not rpool.check_rpool("from-T1w_to-template_mode-image_xfm"): + if not rpool.check_rpool(f"from-{space}_to-template_mode-image_xfm"): reg_blocks = [ [register_ANTs_anat_to_template, register_FSL_anat_to_template], overwrite_transform_anat_to_template, - warp_wholeheadT1_to_template, - warp_T1mask_to_template, ] + if space == "T1w": + reg_blocks += [ + warp_to_template("wholehead", space), + warp_to_template("mask", space), + ] if not rpool.check_rpool("desc-restore-brain_T1w"): reg_blocks.append(correct_restore_brain_intensity_abcd) if cfg.voxel_mirrored_homotopic_connectivity["run"]: - if not rpool.check_rpool("from-T1w_to-symtemplate_mode-image_xfm"): + if not rpool.check_rpool(f"from-{space}_to-symtemplate_mode-image_xfm"): reg_blocks.append( [ register_symmetric_ANTs_anat_to_template, @@ -1223,8 +1204,9 @@ def build_workflow(subject_id, sub_dict, cfg, pipeline_name=None): from CPAC.utils.datasource import gather_extraction_maps # Workflow setup - wf = initialize_nipype_wf(cfg, sub_dict, name=pipeline_name) - + wf = initialize_nipype_wf( + cfg, sub_dict["subject_id"], sub_dict.get("unique_id", None), name=pipeline_name + ) # Extract credentials path if it exists try: creds_path = sub_dict["creds_path"] diff --git a/CPAC/pipeline/cpac_runner.py b/CPAC/pipeline/cpac_runner.py index e5eef08138..166ff36faf 100644 --- a/CPAC/pipeline/cpac_runner.py +++ b/CPAC/pipeline/cpac_runner.py @@ -18,22 +18,23 @@ from multiprocessing import Process import os +from pathlib import Path from time import strftime import warnings from voluptuous.error import Invalid import yaml -from CPAC.longitudinal_pipeline.longitudinal_workflow import anat_longitudinal_wf +from CPAC.longitudinal.wf.anat import anat_longitudinal_wf from CPAC.pipeline.utils import get_shell from CPAC.utils.configuration import check_pname, Configuration, set_subject from CPAC.utils.configuration.yaml_template import upgrade_pipeline_to_1_8 from CPAC.utils.ga import track_run -from CPAC.utils.monitoring import failed_to_start, log_nodes_cb, WFLOGGER +from CPAC.utils.monitoring import failed_to_start, init_loggers, log_nodes_cb, WFLOGGER -# Run condor jobs def run_condor_jobs(c, config_file, subject_list_file, p_name): + """Run condor jobs.""" # Import packages import subprocess from time import strftime @@ -236,7 +237,7 @@ def run_cpac_on_cluster(config_file, subject_list_file, cluster_files_dir): f.write(pid) -def run_T1w_longitudinal(sublist, cfg): +def run_T1w_longitudinal(sublist, cfg: Configuration, dry_run: bool = False): subject_id_dict = {} for sub in sublist: @@ -249,7 +250,11 @@ def run_T1w_longitudinal(sublist, cfg): # sessions for each participant as value for subject_id, sub_list in subject_id_dict.items(): if len(sub_list) > 1: - anat_longitudinal_wf(subject_id, sub_list, cfg) + log_dir: str + _, _, log_dir = set_subject(sub_list[0], cfg) + log_dir = str(Path(log_dir).parent / f"{subject_id}_longitudinal") + init_loggers(subject_id, cfg, log_dir, mock=True) + anat_longitudinal_wf(subject_id, sub_list, cfg, dry_run=dry_run) elif len(sub_list) == 1: warnings.warn( "\n\nThere is only one anatomical session " @@ -491,161 +496,10 @@ def run( """ # BEGIN LONGITUDINAL TEMPLATE PIPELINE - if ( - hasattr(c, "longitudinal_template_generation") - and c.longitudinal_template_generation["run"] - ): - run_T1w_longitudinal(sublist, c) + if c["longitudinal_template_generation", "run"]: + run_T1w_longitudinal(sublist, c, dry_run=test_config) # TODO functional longitudinal pipeline - """ - if valid_longitudinal_data: - rsc_file_list = [] - for dirpath, dirnames, filenames in os.walk(c.pipeline_setup[ - 'output_directory']['path']): - for f in filenames: - # TODO is there a better way to check output folder name? - if f != '.DS_Store' and 'T1w_longitudinal_pipeline' in dirpath: - rsc_file_list.append(os.path.join(dirpath, f)) - - subject_specific_dict = {subj: [] for subj in subject_id_dict.keys()} - session_specific_dict = {os.path.join(session['subject_id'], session['unique_id']): [] for session in sublist} - for rsc_path in rsc_file_list: - key = [s for s in session_specific_dict.keys() if s in rsc_path] - if key: - session_specific_dict[key[0]].append(rsc_path) - else: - subj = [s for s in subject_specific_dict.keys() if s in rsc_path] - if subj: - subject_specific_dict[subj[0]].append(rsc_path) - - # update individual-specific outputs: - # anatomical_brain, anatomical_brain_mask and anatomical_reorient - for key in session_specific_dict.keys(): - for f in session_specific_dict[key]: - sub, ses = key.split('/') - ses_list = [subj for subj in sublist if sub in subj['subject_id'] and ses in subj['unique_id']] - if len(ses_list) > 1: - raise Exception("There are several files containing " + f) - if len(ses_list) == 1: - ses = ses_list[0] - subj_id = ses['subject_id'] - tmp = f.split(c.pipeline_setup['output_directory']['path'])[-1] - keys = tmp.split(os.sep) - if keys[0] == '': - keys = keys[1:] - if len(keys) > 1: - if ses.get('resource_pool') is None: - ses['resource_pool'] = { - keys[0].split(c.pipeline_setup['pipeline_name'] + '_')[-1]: { - keys[-2]: f - } - } - else: - strat_key = keys[0].split(c.pipeline_setup['pipeline_name'] + '_')[-1] - if ses['resource_pool'].get(strat_key) is None: - ses['resource_pool'].update({ - strat_key: { - keys[-2]: f - } - }) - else: - ses['resource_pool'][strat_key].update({ - keys[-2]: f - }) - - for key in subject_specific_dict: - for f in subject_specific_dict[key]: - ses_list = [subj for subj in sublist if key in subj['anat']] - for ses in ses_list: - tmp = f.split(c.pipeline_setup['output_directory']['path'])[-1] - keys = tmp.split(os.sep) - if keys[0] == '': - keys = keys[1:] - if len(keys) > 1: - if ses.get('resource_pool') is None: - ses['resource_pool'] = { - keys[0].split(c.pipeline_setup['pipeline_name'] + '_')[-1]: { - keys[-2]: f - } - } - else: - strat_key = keys[0].split(c.pipeline_setup['pipeline_name'] + '_')[-1] - if ses['resource_pool'].get(strat_key) is None: - ses['resource_pool'].update({ - strat_key: { - keys[-2]: f - } - }) - else: - if keys[-2] == 'anatomical_brain' or keys[-2] == 'anatomical_brain_mask' or keys[-2] == 'anatomical_skull_leaf': - pass - elif 'apply_warp_anat_longitudinal_to_standard' in keys[-2] or 'fsl_apply_xfm_longitudinal' in keys[-2]: - # TODO update!!! - # it assumes session id == last key (ordered by session count instead of session id) + 1 - # might cause problem if session id is not continuous - def replace_index(target1, target2, file_path): - index1 = file_path.index(target1)+len(target1) - index2 = file_path.index(target2)+len(target2) - file_str_list = list(file_path) - file_str_list[index1] = "*" - file_str_list[index2] = "*" - file_path_updated = "".join(file_str_list) - file_list = glob.glob(file_path_updated) - file_list.sort() - return file_list - if ses['unique_id'] == str(int(keys[-2][-1])+1): - if keys[-3] == 'seg_probability_maps': - f_list = replace_index('seg_probability_maps_', 'segment_prob_', f) - ses['resource_pool'][strat_key].update({ - keys[-3]: f_list - }) - elif keys[-3] == 'seg_partial_volume_files': - f_list = replace_index('seg_partial_volume_files_', 'segment_pve_', f) - ses['resource_pool'][strat_key].update({ - keys[-3]: f_list - }) - else: - ses['resource_pool'][strat_key].update({ - keys[-3]: f # keys[-3]: 'anatomical_to_standard' - }) - elif keys[-2] != 'warp_list': - ses['resource_pool'][strat_key].update({ - keys[-2]: f - }) - elif keys[-2] == 'warp_list': - if 'ses-'+ses['unique_id'] in tmp: - ses['resource_pool'][strat_key].update({ - keys[-2]: f - }) - for key in subject_specific_dict: - ses_list = [subj for subj in sublist if key in subj['anat']] - for ses in ses_list: - for reg_strat in strat_list: - try: - ss_strat_list = list(ses['resource_pool']) - for strat_key in ss_strat_list: - try: - ses['resource_pool'][strat_key].update({ - 'registration_method': reg_strat['registration_method'] - }) - except KeyError: - pass - except KeyError: - pass - - yaml.dump(sublist, open(os.path.join(c.pipeline_setup['working_directory']['path'],'data_config_longitudinal.yml'), 'w'), default_flow_style=False) - WFLOGGER.info("\n\nLongitudinal pipeline completed.\n\n") - - # skip main preprocessing - if ( - not c.anatomical_preproc['run'] and - not c.functional_preproc['run'] - ): - sys.exit() - """ - # END LONGITUDINAL TEMPLATE PIPELINE - # If it only allows one, run it linearly if c.pipeline_setup["system_config"]["num_participants_at_once"] == 1: for sub in sublist: diff --git a/CPAC/pipeline/engine.py b/CPAC/pipeline/engine.py index 8b3c726b14..c8631d4d68 100644 --- a/CPAC/pipeline/engine.py +++ b/CPAC/pipeline/engine.py @@ -22,7 +22,8 @@ import json import os import re -from typing import Optional +from types import NotImplementedType +from typing import cast, Generator, Literal, Optional import warnings import pandas as pd @@ -35,6 +36,7 @@ fisher_z_score_standardize, z_score_standardize, ) +from CPAC.longitudinal.wf.utils import LONGITUDINAL_TEMPLATE_PATTERN from CPAC.pipeline import nipype_pipeline_engine as pe from CPAC.pipeline.check_outputs import ExpectedOutputs from CPAC.pipeline.nodeblock import NodeBlockFunction @@ -43,11 +45,12 @@ name_fork, source_set, ) -from CPAC.registration.registration import transform_derivative +from CPAC.registration.utils import transform_derivative from CPAC.resources.templates.lookup_table import lookup_identifier from CPAC.utils.bids_utils import res_in_filename from CPAC.utils.configuration import Configuration from CPAC.utils.datasource import ( + bidsier_prefix, create_anat_datasource, create_func_datasource, create_general_datasource, @@ -62,6 +65,7 @@ WARNING_FREESURFER_OFF_WITH_DATA, WFLOGGER, ) +from CPAC.utils.monitoring.custom_logging import MockLogger from CPAC.utils.outputs import Outputs from CPAC.utils.utils import ( check_prov_for_regtool, @@ -73,7 +77,9 @@ class ResourcePool: - def __init__(self, rpool=None, name=None, cfg=None, pipe_list=None): + def __init__( + self, rpool=None, name: Optional[str] = None, cfg=None, pipe_list=None + ): if not rpool: self.rpool = {} else: @@ -84,7 +90,7 @@ def __init__(self, rpool=None, name=None, cfg=None, pipe_list=None): else: self.pipe_list = pipe_list - self.name = name + self.name = name or "" self.info = {} if cfg: @@ -148,6 +154,32 @@ def __str__(self) -> str: return f"ResourcePool({self.name}): {list(self.rpool)}" return f"ResourcePool: {list(self.rpool)}" + def _set_id_parts(self) -> None: + """Set part_id and ses_id.""" + unique_id = self.name + setattr(self, "_part_id", unique_id.split("_")[0]) + if "_" not in unique_id: + setattr(self, "_ses_id", None) + return + ses_id = unique_id.split("_")[1] + if "ses-" not in ses_id: + ses_id = f"ses-{ses_id}" + setattr(self, "_ses_id", ses_id) + + @property + def part_id(self) -> str: + """Access participant ID.""" + if not hasattr(self, "_part_id"): + self._set_id_parts() + return getattr(self, "_part_id") + + @property + def ses_id(self) -> str: + """Access session ID.""" + if not hasattr(self, "_part_id"): + self._set_id_parts() + return getattr(self, "_ses_id") + def append_name(self, name): self.name.append(name) @@ -196,7 +228,9 @@ def back_propogate_template_name( pass return - def get_name(self): + def get_name(self) -> str: + if not hasattr(self, "_part_id"): + self._set_id_parts() return self.name def check_rpool(self, resource): @@ -506,6 +540,28 @@ def get_cpac_provenance(self, resource, strat=None): json_data = self.get_json(resource, strat) return json_data["CpacProvenance"] + def motion_tool( + self, resource, strat=None + ) -> Optional[Literal["3dvolreg", "mcflirt"]]: + """Check provenance for motion correction tool.""" + prov = self.get_cpac_provenance(resource, strat) + last_entry = get_last_prov_entry(prov) + last_node = last_entry.split(":")[1] + if "3dvolreg" in last_node.lower(): + return "3dvolreg" + if "mcflirt" in last_node.lower(): + return "mcflirt" + # check entire prov + if "3dvolreg" in str(prov): + return "3dvolreg" + if "mcflirt" in str(prov): + return "mcflirt" + return None + + def reg_tool(self, resource, strat=None) -> Optional[Literal["ants", "fsl"]]: + """Check provenance for registration tool.""" + return check_prov_for_regtool(self.get_cpac_provenance(resource, strat)) + @staticmethod def generate_prov_string(prov): # this will generate a string from a SINGLE RESOURCE'S dictionary of @@ -855,7 +911,6 @@ def derivative_xfm(self, wf, label, connection, json_info, pipe_idx, pipe_x): self.num_ants_cores, ants_interp=self.ants_interp, fsl_interp=self.fsl_interp, - opt=None, ) wf.connect(connection[0], connection[1], xfm, "inputspec.in_file") @@ -1080,10 +1135,20 @@ def post_process(self, wf, label, connection, json_info, pipe_idx, pipe_x, outs) def gather_pipes(self, wf, cfg, all=False, add_incl=None, add_excl=None): excl = [] substring_excl = [] - outputs_logger = getLogger( - f'{cfg.get("subject_id", getattr(wf, "name", ""))}_expectedOutputs' - ) - expected_outputs = ExpectedOutputs() + try: + unique_id = re.match(r"(.*_)(sub-.*)", wf.name).group(2) # pyright: ignore[reportOptionalMemberAccess] + except (AttributeError, IndexError): + unique_id = cfg.get("subject_id", getattr(wf, "name", "")) + unique_id = bidsier_prefix(unique_id) + outputs_logger = getLogger(f"{unique_id}_expectedOutputs") + expected = {} + if isinstance(outputs_logger, MockLogger): + try: + # load already-expected outputs + expected = outputs_logger.yaml_contents() + except (FileNotFoundError, TypeError): + pass + expected_outputs = ExpectedOutputs(expected) if add_excl: excl += add_excl @@ -1101,13 +1166,20 @@ def gather_pipes(self, wf, cfg, all=False, add_incl=None, add_excl=None): excl += Outputs.debugging for resource in self.rpool.keys(): - if resource not in Outputs.any: + output_resource: str = ( + resource[22:] + if resource.startswith("longitudinal-template_") + else resource + ) + + if output_resource not in Outputs.any: continue if resource in excl: continue drop = False + for substring_list in substring_excl: bool_list = [] for substring in substring_list: @@ -1134,22 +1206,26 @@ def gather_pipes(self, wf, cfg, all=False, add_incl=None, add_excl=None): # TODO: other stuff like acq- etc. for pipe_idx in self.rpool[resource]: - unique_id = self.get_name() - part_id = unique_id.split("_")[0] - ses_id = unique_id.split("_")[1] - - if "ses-" not in ses_id: - ses_id = f"ses-{ses_id}" - out_dir = cfg.pipeline_setup["output_directory"]["path"] pipe_name = cfg.pipeline_setup["pipeline_name"] - container = os.path.join(f"pipeline_{pipe_name}", part_id, ses_id) - filename = f"{unique_id}_{res_in_filename(self.cfg, resource)}" + longitudinal_xfm = bool( + re.search(LONGITUDINAL_TEMPLATE_PATTERN, resource) + ) + if self.ses_id and not longitudinal_xfm: + container = os.path.join( + f"pipeline_{pipe_name}", self.part_id, self.ses_id + ) + else: + container = os.path.join(f"pipeline_{pipe_name}", self.part_id) + resource_name = self.get_name() + if resource_name.startswith("longitudinal-template_"): + resource_name = resource_name[22:] + filename = f"{resource_name}_{res_in_filename(self.cfg, resource)}" out_path = os.path.join(out_dir, container, subdir, filename) out_dct = { - "unique_id": unique_id, + "unique_id": self.part_id if longitudinal_xfm else self.get_name(), "out_dir": out_dir, "container": container, "subdir": subdir, @@ -1160,7 +1236,6 @@ def gather_pipes(self, wf, cfg, all=False, add_incl=None, add_excl=None): # TODO: have to link the pipe_idx's here. and call up 'desc-preproc_T1w' from a Sources in a json and replace. here. # TODO: can do the pipeline_description.json variants here too! - for resource in self.rpool.keys(): if resource not in Outputs.any: continue @@ -1236,7 +1311,11 @@ def gather_pipes(self, wf, cfg, all=False, add_incl=None, add_excl=None): unlabelled.remove(key) # del all_forks for pipe_idx in self.rpool[resource]: - pipe_x = self.get_pipe_number(pipe_idx) + try: + pipe_x = self.get_pipe_number(pipe_idx) + except ValueError: + # already gone + continue json_info = self.rpool[resource][pipe_idx]["json"] out_dct = self.rpool[resource][pipe_idx]["out"] @@ -1293,7 +1372,7 @@ def gather_pipes(self, wf, cfg, all=False, add_incl=None, add_excl=None): output_names=["out_filename"], function=create_id_string, ), - name=f"id_string_{resource_idx}_{pipe_x}", + name=f"id_string_{unique_id}_{resource_idx}_{pipe_x}", ) id_string.inputs.cfg = self.cfg id_string.inputs.unique_id = unique_id @@ -1347,7 +1426,9 @@ def gather_pipes(self, wf, cfg, all=False, add_incl=None, add_excl=None): ) ) ) - nii_name = pe.Node(Rename(), name=f"nii_{resource_idx}_{pipe_x}") + nii_name = pe.Node( + Rename(), name=f"nii_{unique_id}_{resource_idx}_{pipe_x}" + ) nii_name.inputs.keep_ext = True if resource in Outputs.ciftis: @@ -1366,14 +1447,7 @@ def gather_pipes(self, wf, cfg, all=False, add_incl=None, add_excl=None): wf.connect(id_string, "out_filename", nii_name, "format_string") node, out = self.rpool[resource][pipe_idx]["data"] - if not node: - msg = f"Resource {resource} not found in resource pool." - raise FileNotFoundError(msg) - try: - wf.connect(node, out, nii_name, "in_file") - except OSError as os_error: - WFLOGGER.warning(os_error) - continue + wf.connect(node, out, nii_name, "in_file") write_json_imports = ["import os", "import json"] write_json = pe.Node( @@ -2635,7 +2709,14 @@ def _set_nested(attr, keys): return wf, rpool -def initiate_rpool(wf, cfg, data_paths=None, part_id=None): +def initiate_rpool( + wf: pe.Workflow, + cfg: Configuration, + data_paths=None, + part_id=None, + *, + rpool: Optional[ResourcePool] = None, +) -> tuple[pe.Workflow, ResourcePool]: """ Initialize a new ResourcePool. @@ -2674,7 +2755,7 @@ def initiate_rpool(wf, cfg, data_paths=None, part_id=None): unique_id = part_id creds_path = None - rpool = ResourcePool(name=unique_id, cfg=cfg) + rpool = ResourcePool(rpool=rpool.rpool if rpool else None, name=unique_id, cfg=cfg) if data_paths: # ingress outdir @@ -2707,7 +2788,7 @@ def initiate_rpool(wf, cfg, data_paths=None, part_id=None): # output files with 4 different scans - return (wf, rpool) + return wf, rpool def run_node_blocks(blocks, data_paths, cfg=None): @@ -2768,13 +2849,13 @@ class NodeData: -------- >>> rp = ResourcePool() >>> rp.node_data(None) - NotImplemented (NotImplemented) + NodeData(NotImplemented, NotImplemented) >>> rp.set_data('test', ... pe.Node(Function(input_names=[]), 'test'), ... 'b', [], 0, 'test') >>> rp.node_data('test') - test (b) + NodeData(test, b) >>> rp.node_data('test').out 'b' @@ -2787,10 +2868,25 @@ class NodeData: # pylint: disable=too-few-public-methods def __init__(self, strat_pool=None, resource=None, **kwargs): + """Initialize NodeData.""" self.node = NotImplemented self.out = NotImplemented if strat_pool is not None and resource is not None: - self.node, self.out = strat_pool.get_data(resource, **kwargs) + self.node, self.out = cast( + tuple[pe.Node, str], strat_pool.get_data(resource, **kwargs) + ) - def __repr__(self): # noqa: D105 + def __iter__( + self, + ) -> Generator[pe.Node | NotImplementedType, str | NotImplementedType, None]: + """Expand NodeData into node, data.""" + yield self.node + yield self.out + + def __repr__(self) -> str: + """Return reproducible string representation of NodeData.""" + return f"NodeData({getattr(self.node, 'name', str(self.node))}, {self.out})" + + def __str__(self) -> str: + """Return string representation of NodeData.""" return f'{getattr(self.node, "name", str(self.node))} ({self.out})' diff --git a/CPAC/pipeline/nipype_pipeline_engine/engine.py b/CPAC/pipeline/nipype_pipeline_engine/engine.py index 743285ae9d..455e38280c 100644 --- a/CPAC/pipeline/nipype_pipeline_engine/engine.py +++ b/CPAC/pipeline/nipype_pipeline_engine/engine.py @@ -8,6 +8,7 @@ # * Applies a random seed # * Supports overriding memory estimates via a log file and a buffer # * Adds quotation marks around strings in dotfiles +# * Adds methods for cross-graph connections # ORIGINAL WORK'S ATTRIBUTION NOTICE: # Copyright (c) 2009-2016, Nipype developers @@ -50,16 +51,18 @@ for Nipype's documentation. """ # pylint: disable=line-too-long +from collections.abc import Mapping, Sequence from copy import deepcopy from inspect import Parameter, Signature, signature import os import re -from typing import Any, ClassVar, Optional +from typing import Any, cast, ClassVar, Optional, TYPE_CHECKING from numpy import prod from traits.trait_base import Undefined from traits.trait_handlers import TraitListObject from nibabel import load +from nipype.interfaces.base.support import InterfaceResult from nipype.interfaces.utility import Function from nipype.pipeline import engine as pe from nipype.pipeline.engine.utils import ( @@ -76,6 +79,9 @@ from CPAC.utils.monitoring import getLogger, WFLOGGER +if TYPE_CHECKING: + pass + # set global default mem_gb DEFAULT_MEM_GB = 2.0 UNDEFINED_SIZE = (42, 42, 42, 1200) @@ -527,6 +533,25 @@ def __init__(self, name, base_dir=None, debug=False): self._nodes_cache = set() self._nested_workflows_cache = set() + def copy_input_connections(self, node1: pe.Node, node2: pe.Node) -> None: + """Copy input connections from ``node1`` to ``node2``.""" + new_connections: list[tuple[pe.Node, str, pe.Node, str]] = [] + for connection in self._graph.edges: + _out: pe.Node + _in: pe.Node + _out, _in = connection + if _in == node1: + details = self._graph.get_edge_data(*connection) + if "connect" in details: + for connect in details["connect"]: + new_connections.append((_out, connect[0], node2, connect[1])) + for connection in new_connections: + try: + self.connect(*connection) + except Exception: + # connection already exists + continue + def _configure_exec_nodes(self, graph): """Ensure that each node knows where to get inputs from.""" for node in graph.nodes(): @@ -565,6 +590,20 @@ def _configure_exec_nodes(self, graph): except (FileNotFoundError, KeyError, TypeError): self._handle_just_in_time_exception(node) + def _connect_node_or_path_for_merge( + self, + node: pe.Node, + strats_dct: Mapping[str, Sequence[tuple[pe.Node, str] | str]], + key: str, + index: int, + ) -> None: + """Set input to either a Node or a path string for cross-graph Merge Nodes.""" + _input: str = f"in{index + 1}" + if isinstance(strats_dct[key][index], str): + node.set_input(_input, strats_dct[key][index]) + else: + self.connect(*strats_dct[key][index], node, _input) + def _get_dot( self, prefix=None, hierarchy=None, colored=False, simple_form=True, level=0 ): @@ -678,6 +717,41 @@ def _get_dot( WFLOGGER.debug("cross connection: %s", dotlist[-1]) return ("\n" + prefix).join(dotlist) + def get_output(self, node: pe.Node, out: str) -> Any: + """Get an output path from an already-run Node.""" + result_nodes = cast(list[pe.Node], self.run(updatehash=True).nodes) + orig_wd = os.getcwd() + output = Undefined + try: + # look for exact match + _run_node: pe.Node = next( + iter(_ for _ in result_nodes if _.fullname == node.fullname) + ) + except StopIteration as stop_interation: + # look for match in subgraph + try: + _run_node: pe.Node = next( + iter( + _ + for _ in result_nodes + if node.fullname + and _.fullname + and _.fullname.endswith(node.fullname) + ) + ) + except StopIteration: + msg = f"Could not find {node.fullname} in {self}'s run Nodes." + raise LookupError(msg) from stop_interation + try: + os.chdir(_run_node.output_dir()) + _res: InterfaceResult = _run_node.run() + output = getattr(_res.outputs, out) + if output is Undefined: + output = _run_node.interface._list_outputs().get(out, Undefined) + finally: + os.chdir(orig_wd) + return output + def _handle_just_in_time_exception(self, node): # pylint: disable=protected-access if hasattr(self, "_local_func_scans"): diff --git a/CPAC/pipeline/nodeblock.py b/CPAC/pipeline/nodeblock.py index 53b9db1330..86d1da9d54 100644 --- a/CPAC/pipeline/nodeblock.py +++ b/CPAC/pipeline/nodeblock.py @@ -1,4 +1,4 @@ -# Copyright (C) 2023-2024 C-PAC Developers +# Copyright (C) 2023-2025 C-PAC Developers # This file is part of C-PAC. @@ -16,7 +16,13 @@ # License along with C-PAC. If not, see . """Class and decorator for NodeBlock functions.""" -from typing import Any, Callable, Optional +from typing import Any, Callable, Optional, TypeAlias + +from nipype.pipeline import engine as pe + +NODEBLOCK_RETURN: TypeAlias = tuple[ + pe.Workflow, dict[str, tuple[pe.Node | pe.Workflow, str]] +] class NodeBlockFunction: @@ -55,9 +61,9 @@ def __init__( """ self.option_val: Optional[str | list[str]] = option_val """Indicates values for which this NodeBlock should be active.""" - self.inputs: Optional[list[str | list | tuple]] = inputs + self.inputs: list[str | list | tuple] = inputs or [] """ResourcePool keys indicating resources needed for the NodeBlock's functionality.""" - self.outputs: Optional[list[str] | dict[str, Any]] = outputs + self.outputs: list[str] | dict[str, Any] = outputs or [] """ ResourcePool keys indicating resources generated or updated by the NodeBlock, optionally including metadata for the outputs' respective sidecars. @@ -78,7 +84,7 @@ def __init__( ).rstrip() # all node block functions have this signature - def __call__(self, wf, cfg, strat_pool, pipe_num, opt=None): + def __call__(self, wf, cfg, strat_pool, pipe_num, opt=None) -> NODEBLOCK_RETURN: """ Parameters diff --git a/CPAC/pipeline/schema.py b/CPAC/pipeline/schema.py index 828c0b1aec..112492ae4d 100644 --- a/CPAC/pipeline/schema.py +++ b/CPAC/pipeline/schema.py @@ -21,6 +21,7 @@ from itertools import chain, permutations import re from subprocess import CalledProcessError +from typing import Any as AnyType import numpy as np from pathvalidate import sanitize_filename @@ -852,22 +853,38 @@ def sanitize(filename): }, "longitudinal_template_generation": { "run": bool1_1, + "using": In({"mri_robust_template", "C-PAC legacy"}), "average_method": In({"median", "mean", "std"}), "dof": In({12, 9, 7, 6}), - "interp": In({"trilinear", "nearestneighbour", "sinc", "spline"}), - "cost": In( - { - "corratio", - "mutualinfo", - "normmi", - "normcorr", - "leastsq", - "labeldiff", - "bbr", - } + "max_iter": Any( + All(Number, Range(min=0, min_included=False)), In([-1, "default"]) + ), + "legacy-specific": Maybe( + Schema( + { + "convergence_threshold": Any( + All(Number, Range(min=0, max=1, min_included=False)), -1 + ), + "interp": Maybe( + In({"trilinear", "nearestneighbour", "sinc", "spline"}) + ), + "cost": Maybe( + In( + { + "corratio", + "mutualinfo", + "normmi", + "normcorr", + "leastsq", + "labeldiff", + "bbr", + } + ) + ), + "thread_pool": Maybe(int), + } + ) ), - "thread_pool": int, - "convergence_threshold": Number, }, "functional_preproc": { "run": bool1_1, @@ -1267,6 +1284,20 @@ def sanitize(filename): ) +def check_unimplemented( + to_check: dict[str, AnyType], k_v_pairs: list[tuple[str, AnyType]], category: str +) -> None: + """Check for unimplemented combinations in subschema. + + Raise NotImplementedError if any found. + """ + error_msg = "`{value}` is not implemented for {category} `{key}`." + for key, value in k_v_pairs: + if to_check[key] == value: + msg = error_msg.format(category=category, key=key, value=value) + raise NotImplementedError(msg) + + def schema(config_dict): """Validate a participant-analysis pipeline configuration. @@ -1330,7 +1361,9 @@ def schema(config_dict): "``nuisance_corrections: 2-nuisance_regression: space`` " f"to ``template`` {or_else}" ) - raise ExclusiveInvalid(msg) + raise ExclusiveInvalid( + msg, path=["nuisance_corrections", "2-nuisance_regression", "space"] + ) if any( registration != "ANTS" for registration in partially_validated["registration_workflows"][ @@ -1343,7 +1376,15 @@ def schema(config_dict): "``registration_workflows: anatomical_registration: " f"registration: using`` to ``ANTS`` {or_else}" ) - raise ExclusiveInvalid(msg) + raise ExclusiveInvalid( + msg, + path=[ + "registration_workflows", + "anatomical_registration", + "registration", + "using", + ], + ) except KeyError: pass try: @@ -1398,7 +1439,9 @@ def schema(config_dict): "[!] Ingress_regressors and create_regressors can't both run! " " Try turning one option off.\n " ) - raise ExclusiveInvalid(msg) + raise ExclusiveInvalid( + msg, path=["nuisance_corrections", "2-nuisance_regression"] + ) overwrite = partially_validated["registration_workflows"][ "anatomical_registration" @@ -1411,9 +1454,20 @@ def schema(config_dict): "anatomical_registration" ]["registration"]["using"] ): + msg = ( + "[!] Overwrite transform method is the same as the anatomical" + " registration method!\nNo need to overwrite transform with the same" + " registration method. Please turn it off or use a different" + " registration method." + ) raise ExclusiveInvalid( - "[!] Overwrite transform method is the same as the anatomical registration method! " - "No need to overwrite transform with the same registration method. Please turn it off or use a different registration method." + msg, + path=[ + "registration_workflows", + "anatomical_registration", + "overwrite_transform", + "run", + ], ) except KeyError: pass @@ -1445,6 +1499,21 @@ def schema(config_dict): raise OSError(msg) from error except KeyError: pass + try: + # check for incompatible longitudinal options + lgt = partially_validated["longitudinal_template_generation"] + if lgt["using"] == "mri_robust_template": + check_unimplemented( + lgt, + [("average_method", "std"), ("dof", 9), ("dof", 7), ("max_iter", -1)], + "longitudinal `mri_robust_template`", + ) + if lgt["using"] == "C-PAC legacy": + check_unimplemented( + lgt, [("max_iter", "default")], "C-PAC legacy longitudinal" + ) + except KeyError: + pass return partially_validated diff --git a/CPAC/pipeline/test/test_engine.py b/CPAC/pipeline/test/test_engine.py index 25b16d9e44..677c3489a0 100644 --- a/CPAC/pipeline/test/test_engine.py +++ b/CPAC/pipeline/test/test_engine.py @@ -49,11 +49,11 @@ def test_ingress_func_raw_data(pipe_config, bids_dir, test_dir): cfg.pipeline_setup["output_directory"]["path"] = os.path.join(test_dir, "out") cfg.pipeline_setup["working_directory"]["path"] = os.path.join(test_dir, "work") - wf = initialize_nipype_wf(cfg, sub_data_dct) - part_id = sub_data_dct["subject_id"] ses_id = sub_data_dct["unique_id"] + wf = initialize_nipype_wf(cfg, part_id, ses_id) + unique_id = f"{part_id}_{ses_id}" rpool = ResourcePool(name=unique_id, cfg=cfg) @@ -76,11 +76,11 @@ def test_ingress_anat_raw_data(pipe_config, bids_dir, test_dir): cfg.pipeline_setup["output_directory"]["path"] = os.path.join(test_dir, "out") cfg.pipeline_setup["working_directory"]["path"] = os.path.join(test_dir, "work") - wf = initialize_nipype_wf(cfg, sub_data_dct) - part_id = sub_data_dct["subject_id"] ses_id = sub_data_dct["unique_id"] + wf = initialize_nipype_wf(cfg, part_id, ses_id) + unique_id = f"{part_id}_{ses_id}" rpool = ResourcePool(name=unique_id, cfg=cfg) @@ -103,11 +103,11 @@ def test_ingress_pipeconfig_data(pipe_config, bids_dir, test_dir): cfg.pipeline_setup["working_directory"]["path"] = os.path.join(test_dir, "work") cfg.pipeline_setup["log_directory"]["path"] = os.path.join(test_dir, "logs") - wf = initialize_nipype_wf(cfg, sub_data_dct) - part_id = sub_data_dct["subject_id"] ses_id = sub_data_dct["unique_id"] + wf = initialize_nipype_wf(cfg, part_id, ses_id) + unique_id = f"{part_id}_{ses_id}" rpool = ResourcePool(name=unique_id, cfg=cfg) @@ -128,7 +128,9 @@ def test_build_anat_preproc_stack(pipe_config, bids_dir, test_dir): cfg.pipeline_setup["working_directory"]["path"] = os.path.join(test_dir, "work") cfg.pipeline_setup["log_directory"]["path"] = os.path.join(test_dir, "logs") - wf = initialize_nipype_wf(cfg, sub_data_dct) + wf = initialize_nipype_wf( + cfg, sub_data_dct["subject_id"], sub_data_dct["unique_id"] + ) wf, rpool = initiate_rpool(wf, cfg, sub_data_dct) @@ -149,7 +151,9 @@ def test_build_workflow(pipe_config, bids_dir, test_dir): cfg.pipeline_setup["working_directory"]["path"] = os.path.join(test_dir, "work") cfg.pipeline_setup["log_directory"]["path"] = os.path.join(test_dir, "logs") - wf = initialize_nipype_wf(cfg, sub_data_dct) + wf = initialize_nipype_wf( + cfg, sub_data_dct["subject_id"], sub_data_dct["unique_id"] + ) wf, rpool = initiate_rpool(wf, cfg, sub_data_dct) diff --git a/CPAC/pipeline/utils.py b/CPAC/pipeline/utils.py index d135addc41..ffa5700a9d 100644 --- a/CPAC/pipeline/utils.py +++ b/CPAC/pipeline/utils.py @@ -16,14 +16,18 @@ # License along with C-PAC. If not, see . """C-PAC pipeline engine utilities.""" +from collections.abc import Sequence from itertools import chain import os import subprocess -from typing import Optional +from typing import Optional, TYPE_CHECKING from CPAC.func_preproc.func_motion import motion_estimate_filter from CPAC.utils.bids_utils import insert_entity +if TYPE_CHECKING: + from CPAC.pipeline.nipype_pipeline_engine import Node, Workflow + MOVEMENT_FILTER_KEYS = motion_estimate_filter.outputs @@ -241,3 +245,14 @@ def _update_resource_idx(resource_idx, out_dct, key, value): resource_idx = insert_entity(resource_idx, key, value) out_dct["filename"] = insert_entity(out_dct["filename"], key, value) return resource_idx, out_dct + + +def get_edges_with_node( + wf: "Workflow", parameter: str +) -> Sequence[tuple["Node", "Node", dict[str, list[tuple[str, str]]]]]: + """Get all edges containing a given parameter.""" + return [ + edge + for edge in wf._graph.edges(data=True) + if any(parameter in _ for _ in edge[2].get("connect")) + ] diff --git a/CPAC/registration/longitudinal.py b/CPAC/registration/longitudinal.py new file mode 100644 index 0000000000..74dacd6504 --- /dev/null +++ b/CPAC/registration/longitudinal.py @@ -0,0 +1,208 @@ +# Copyright (C) 2025 C-PAC Developers + +# This file is part of C-PAC. + +# C-PAC is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or (at your +# option) any later version. + +# C-PAC is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +# License for more details. + +# You should have received a copy of the GNU Lesser General Public +# License along with C-PAC. If not, see . +# pylint: disable=too-many-lines,ungrouped-imports,wrong-import-order +"""Longitudial registration workflows and utilities.""" + +from nipype.interfaces.utility import IdentityInterface +from nipype.pipeline.engine import Node, Workflow + +from CPAC.pipeline import nipype_pipeline_engine as pe +from CPAC.pipeline.engine import ResourcePool +from CPAC.registration.utils import ( + apply_transform, + CommonRegistrationInputs, + compose_ants_warp, + RegistrationTemplates, +) + + +def get_common_reg_inputs( + strat_pool: ResourcePool, + registration_templates: RegistrationTemplates = RegistrationTemplates(), +) -> CommonRegistrationInputs: + """Get common longitudinal registration inputs.""" + orig = "longitudinal" + has_longitudinal = strat_pool.check_rpool( + "longitudinal-template_space-longitudinal_desc-brain_T1w" + ) + is_longitudinal = False + if has_longitudinal and not strat_pool.check_rpool("desc-preproc_T1w"): + orig = "longitudinal" + has_longitudinal = False + is_longitudinal = True + input_brain = strat_pool.node_data( + "longitudinal-template_space-longitudinal_desc-brain_T1w" + ) + input_head = strat_pool.node_data( + "longitudinal-template_space-longitudinal_desc-head_T1w" + ) + reference_mask = strat_pool.node_data( + "longitudinal-template_space-longitudinal_desc-brain_mask" + ) + lesion_mask = None + else: + orig = "T1w" + input_brain = strat_pool.node_data("desc-preproc_T1w") + input_head = strat_pool.node_data( + [ # TODO: check the order of T1w + "desc-restore_T1w", + "desc-head_T1w", + "desc-preproc_T1w", + ] + ) + reference_mask = ( + strat_pool.node_data(registration_templates.reference_mask) + if strat_pool.check_rpool(registration_templates.reference_mask) + else None + ) + lesion_mask = ( + strat_pool.node_data("label-lesion_mask") + if strat_pool.check_rpool("label-lesion_mask") + else None + ) + if has_longitudinal: + t1w_brain_template = strat_pool.node_data( + "longitudinal-template_space-longitudinal_desc-brain_T1w" + ) + t1w_template = strat_pool.node_data( + "longitudinal-template_space-longitudinal_desc-head_T1w" + ) + else: + t1w_brain_template = strat_pool.node_data( + registration_templates.reference_brain + ) + t1w_template = strat_pool.node_data(registration_templates.reference_head) + if is_longitudinal: + brain_mask = strat_pool.node_data( + "longitudinal-template_space-longitudinal_desc-brain_mask" + ) + else: + brain_mask = strat_pool.node_data( + [ + "space-T1w_desc-brain_mask", + "space-T1w_desc-acpcbrain_mask", + ] + ) + return CommonRegistrationInputs( + orig, + has_longitudinal, + input_brain, + input_head, + reference_mask, + lesion_mask, + t1w_brain_template, + t1w_template, + brain_mask, + ) + + +def t1w_to_longitudinal_to_template_ants( + wf: Workflow, strat_pool: ResourcePool, template: str = "template" +) -> dict[str, tuple[Node | Workflow, str]]: + """Combine and apply T1w transforms from native to longitudinal to template.""" + t1_to_long_to_template = pe.Node( + IdentityInterface( + fields=[ + "input_brain", + "reference_brain", + "input_head", + "reference_head", + "input_mask", + "reference_mask", + "transform", + "interpolation", + ], + interpolation=strat_pool.ants_interp, + ), + name=f"t1_to_long_to_{template}_inputspec", + ) + list_of_xfms = pe.Node( + IdentityInterface(fields=["T1w_to_longitudinal", "longitudinal_to_template"]), + "list_of_xfms", + ) + wf.connect( + *strat_pool.node_data("from-T1w_to-longitudinal_mode-image_xfm"), + list_of_xfms, + "T1w_to_longitudinal", + ) + wf.connect( + *strat_pool.node_data(f"from-longitudinal_to-{template}_mode-image_xfm"), + list_of_xfms, + "longitudinal_to_template", + ) + composite_xfm = compose_ants_warp( + wf=wf, + name="T1wLinearTemplate_xfm", + input_node=t1_to_long_to_template, + warp_from="input_brain", + warp_to="reference_brain", + inputs=[(list_of_xfms, ["longitudinal_to_template", "T1w_to_longitudinal"])], + ) + list_of_invxfms = pe.Node( + IdentityInterface(fields=["longitudinal_to_T1w", "template_to_longitudinal"]), + "list_of_inv_xfms", + ) + wf.connect( + *strat_pool.node_data("from-longitudinal_to-T1w_mode-image_xfm"), + list_of_invxfms, + "longitudinal_to_T1w", + ) + wf.connect( + *strat_pool.node_data(f"from-{template}_to-longitudinal_mode-image_xfm"), + list_of_invxfms, + "template_to_longitudinal", + ) + inverse_composite_xfm = compose_ants_warp( + wf=wf, + name="T1wLinearTemplateInv_xfm", + input_node=t1_to_long_to_template, + warp_from="reference_brain", + warp_to="input_brain", + inputs=[(list_of_invxfms, ["longitudinal_to_T1w", "template_to_longitudinal"])], + ) + preproc = apply_transform(f"warp_t1_to_longitudinal_to_{template}", "ants") + preproc.inputs.inputspec.interpolation = strat_pool.ants_interp + wf.connect( + [ + ( + t1_to_long_to_template, + preproc, + [ + ("input_brain", "inputspec.input_image"), + ("reference_brain", "reference"), + ], + ), + (composite_xfm, preproc, [("output_image", "transform")]), + ] + ) + outputs: dict[str, tuple[Node | Workflow, str]] = { + f"space-{template}_desc-preproc_T1w": (preproc, "output_image"), + # "from-T1w_to-template_mode-image_desc-linear_xfm": None, + # "from-template_to-T1w_mode-image_desc-linear_xfm": None, + # "from-T1w_to-template_mode-image_desc-nonlinear_xfm": None, + # "from-template_to-T1w_mode-image_desc-nonlinear_xfm": None, + f"from-T1w_to-{template}_mode-image_xfm": (composite_xfm, "output_image"), + f"from-{template}_to-T1w_mode-image_xfm": ( + inverse_composite_xfm, + "output_image", + ), + } + + return outputs + + +t1w_to_longitudinal_to_template = {"ants": t1w_to_longitudinal_to_template_ants} diff --git a/CPAC/registration/registration.py b/CPAC/registration/registration.py index 258cb9712d..9103c3690a 100644 --- a/CPAC/registration/registration.py +++ b/CPAC/registration/registration.py @@ -1,4 +1,4 @@ -# Copyright (C) 2012-2024 C-PAC Developers +# Copyright (C) 2012-2025 C-PAC Developers # This file is part of C-PAC. @@ -17,7 +17,7 @@ # pylint: disable=too-many-lines,ungrouped-imports,wrong-import-order """Workflows for registration.""" -from typing import Optional +from typing import Any, cast, Literal, Optional, TYPE_CHECKING from voluptuous import RequiredFieldInvalid from nipype.interfaces import afni, ants, c3, fsl, utility as util @@ -25,344 +25,34 @@ from CPAC.anat_preproc.lesion_preproc import create_lesion_preproc from CPAC.func_preproc.func_preproc import fsl_afni_subworkflow -from CPAC.func_preproc.utils import chunk_ts, split_ts_chunks from CPAC.pipeline import nipype_pipeline_engine as pe -from CPAC.pipeline.nodeblock import nodeblock +from CPAC.pipeline.nodeblock import nodeblock, NODEBLOCK_RETURN, NodeBlockFunction +from CPAC.registration.longitudinal import ( + get_common_reg_inputs, + t1w_to_longitudinal_to_template, +) from CPAC.registration.utils import ( + apply_transform, change_itk_transform_type, - check_transforms, - generate_inverse_transform_flags, + compose_ants_warp, + convert_pedir, hardcoded_reg, - interpolation_string, one_d_to_mat, + prep_reg_connector, + prepend_space, + REGISTRATION_SPACE, + RegistrationTemplates, run_c3d, run_c4d, seperate_warps_list, - single_ants_xfm_to_list, + xfm_outputs, ) from CPAC.utils.interfaces import Function from CPAC.utils.interfaces.fsl import Merge as fslMerge -from CPAC.utils.utils import check_prov_for_motion_tool, check_prov_for_regtool - -def apply_transform( - wf_name, - reg_tool, - time_series=False, - multi_input=False, - num_cpus=1, - num_ants_cores=1, -): - """Apply transform.""" - if not reg_tool: - msg = ( - "\n[!] Developer info: the 'reg_tool' parameter sent to the" - f" 'apply_transform' node for '{wf_name}' is empty.\n" - ) - raise RequiredFieldInvalid(msg) - - wf = pe.Workflow(name=wf_name) - - inputNode = pe.Node( - util.IdentityInterface( - fields=["input_image", "reference", "transform", "interpolation"] - ), - name="inputspec", - ) - - outputNode = pe.Node( - util.IdentityInterface(fields=["output_image"]), name="outputspec" - ) - - if int(num_cpus) > 1 and time_series: - # parallelize time series warp application - # we need the node to be a MapNode to feed in the list of functional - # time series chunks - multi_input = True - - if reg_tool == "ants": - if multi_input: - apply_warp = pe.MapNode( - interface=ants.ApplyTransforms(), - name=f"apply_warp_{wf_name}", - iterfield=["input_image"], - mem_gb=0.7, - mem_x=(1708448960473801 / 151115727451828646838272, "input_image"), - ) - else: - apply_warp = pe.Node( - interface=ants.ApplyTransforms(), - name=f"apply_warp_{wf_name}", - mem_gb=0.7, - mem_x=(1708448960473801 / 151115727451828646838272, "input_image"), - ) - - apply_warp.inputs.dimension = 3 - apply_warp.interface.num_threads = int(num_ants_cores) - - if time_series: - apply_warp.inputs.input_image_type = 3 - - wf.connect(inputNode, "reference", apply_warp, "reference_image") - - interp_string = pe.Node( - Function( - input_names=["interpolation", "reg_tool"], - output_names=["interpolation"], - function=interpolation_string, - ), - name="interp_string", - mem_gb=2.5, - ) - interp_string.inputs.reg_tool = reg_tool - - wf.connect(inputNode, "interpolation", interp_string, "interpolation") - wf.connect(interp_string, "interpolation", apply_warp, "interpolation") - - ants_xfm_list = pe.Node( - Function( - input_names=["transform"], - output_names=["transform_list"], - function=single_ants_xfm_to_list, - ), - name="single_ants_xfm_to_list", - mem_gb=2.5, - ) - - wf.connect(inputNode, "transform", ants_xfm_list, "transform") - wf.connect(ants_xfm_list, "transform_list", apply_warp, "transforms") - - # parallelize the apply warp, if multiple CPUs, and it's a time - # series! - if int(num_cpus) > 1 and time_series: - chunk_imports = ["import nibabel as nib"] - chunk = pe.Node( - Function( - input_names=["func_file", "n_chunks", "chunk_size"], - output_names=["TR_ranges"], - function=chunk_ts, - imports=chunk_imports, - ), - name=f"chunk_{wf_name}", - mem_gb=2.5, - ) - - # chunk.inputs.n_chunks = int(num_cpus) - - # 10-TR sized chunks - chunk.inputs.chunk_size = 10 - - wf.connect(inputNode, "input_image", chunk, "func_file") - - split_imports = ["import os", "import subprocess"] - split = pe.Node( - Function( - input_names=["func_file", "tr_ranges"], - output_names=["split_funcs"], - function=split_ts_chunks, - imports=split_imports, - ), - name=f"split_{wf_name}", - mem_gb=2.5, - ) - - wf.connect(inputNode, "input_image", split, "func_file") - wf.connect(chunk, "TR_ranges", split, "tr_ranges") - - wf.connect(split, "split_funcs", apply_warp, "input_image") - - func_concat = pe.Node( - interface=afni_utils.TCat(), name=f"func_concat_{wf_name}", mem_gb=2.5 - ) - func_concat.inputs.outputtype = "NIFTI_GZ" - - wf.connect(apply_warp, "output_image", func_concat, "in_files") - - wf.connect(func_concat, "out_file", outputNode, "output_image") - - else: - wf.connect(inputNode, "input_image", apply_warp, "input_image") - wf.connect(apply_warp, "output_image", outputNode, "output_image") - - elif reg_tool == "fsl": - if multi_input: - apply_warp = pe.MapNode( - interface=fsl.ApplyWarp(), - name="fsl_apply_warp", - iterfield=["in_file"], - mem_gb=2.5, - ) - else: - apply_warp = pe.Node( - interface=fsl.ApplyWarp(), name="fsl_apply_warp", mem_gb=2.5 - ) - - interp_string = pe.Node( - Function( - input_names=["interpolation", "reg_tool"], - output_names=["interpolation"], - function=interpolation_string, - ), - name="interp_string", - mem_gb=2.5, - ) - interp_string.inputs.reg_tool = reg_tool - - wf.connect(inputNode, "interpolation", interp_string, "interpolation") - wf.connect(interp_string, "interpolation", apply_warp, "interp") - - # mni to t1 - wf.connect(inputNode, "reference", apply_warp, "ref_file") - - # NOTE: C-PAC now converts all FSL xfm's to .nii, so even if the - # inputNode 'transform' is a linear xfm, it's a .nii and must - # go in as a warpfield file - wf.connect(inputNode, "transform", apply_warp, "field_file") - - # parallelize the apply warp, if multiple CPUs, and it's a time - # series! - if int(num_cpus) > 1 and time_series: - chunk_imports = ["import nibabel as nib"] - chunk = pe.Node( - Function( - input_names=["func_file", "n_chunks", "chunk_size"], - output_names=["TR_ranges"], - function=chunk_ts, - imports=chunk_imports, - ), - name=f"chunk_{wf_name}", - mem_gb=2.5, - ) - - # chunk.inputs.n_chunks = int(num_cpus) - - # 10-TR sized chunks - chunk.inputs.chunk_size = 10 - - wf.connect(inputNode, "input_image", chunk, "func_file") - - split_imports = ["import os", "import subprocess"] - split = pe.Node( - Function( - input_names=["func_file", "tr_ranges"], - output_names=["split_funcs"], - function=split_ts_chunks, - imports=split_imports, - ), - name=f"split_{wf_name}", - mem_gb=2.5, - ) - - wf.connect(inputNode, "input_image", split, "func_file") - wf.connect(chunk, "TR_ranges", split, "tr_ranges") - - wf.connect(split, "split_funcs", apply_warp, "in_file") - - func_concat = pe.Node( - interface=afni_utils.TCat(), name=f"func_concat{wf_name}" - ) - func_concat.inputs.outputtype = "NIFTI_GZ" - - wf.connect(apply_warp, "out_file", func_concat, "in_files") - - wf.connect(func_concat, "out_file", outputNode, "output_image") - - else: - wf.connect(inputNode, "input_image", apply_warp, "in_file") - wf.connect(apply_warp, "out_file", outputNode, "output_image") - - return wf - - -def transform_derivative( - wf_name, - label, - reg_tool, - num_cpus, - num_ants_cores, - ants_interp=None, - fsl_interp=None, - opt=None, -): - """Transform output derivatives to template space. - - This function is designed for use with the NodeBlock connection engine. - """ - wf = pe.Workflow(name=wf_name) - - inputnode = pe.Node( - util.IdentityInterface(fields=["in_file", "reference", "transform"]), - name="inputspec", - ) - - multi_input = False - if "statmap" in label: - multi_input = True - - stack = False - if "correlations" in label: - stack = True - - apply_xfm = apply_transform( - f"warp_{label}_to_template", - reg_tool, - time_series=stack, - multi_input=multi_input, - num_cpus=num_cpus, - num_ants_cores=num_ants_cores, - ) - - if reg_tool == "ants": - apply_xfm.inputs.inputspec.interpolation = ants_interp - elif reg_tool == "fsl": - apply_xfm.inputs.inputspec.interpolation = fsl_interp - - wf.connect(inputnode, "in_file", apply_xfm, "inputspec.input_image") - wf.connect(inputnode, "reference", apply_xfm, "inputspec.reference") - wf.connect(inputnode, "transform", apply_xfm, "inputspec.transform") - - outputnode = pe.Node(util.IdentityInterface(fields=["out_file"]), name="outputspec") - - wf.connect(apply_xfm, "outputspec.output_image", outputnode, "out_file") - - return wf - - -def convert_pedir(pedir, convert="xyz_to_int"): - """FSL Flirt requires pedir input encoded as an int.""" - if convert == "xyz_to_int": - conv_dct = { - "x": 1, - "y": 2, - "z": 3, - "x-": -1, - "y-": -2, - "z-": -3, - "i": 1, - "j": 2, - "k": 3, - "i-": -1, - "j-": -2, - "k-": -3, - "-x": -1, - "-i": -1, - "-y": -2, - "-j": -2, - "-z": -3, - "-k": -3, - } - elif convert == "ijk_to_xyz": - conv_dct = {"i": "x", "j": "y", "k": "z", "i-": "x-", "j-": "y-", "k-": "z-"} - - if isinstance(pedir, bytes): - pedir = pedir.decode() - if not isinstance(pedir, str): - msg = f"\n\nPhase-encoding direction must be a string value.\n\nValue: {pedir}\n\n" - raise ValueError(msg) - if pedir not in conv_dct.keys(): - msg = f"\n\nInvalid phase-encoding direction entered: {pedir}\n\n" - raise ValueError(msg) - return conv_dct[pedir] +if TYPE_CHECKING: + from CPAC.pipeline.engine import ResourcePool + from CPAC.utils.configuration import Configuration def create_fsl_flirt_linear_reg(name="fsl_flirt_linear_reg"): @@ -1106,23 +796,24 @@ def bbreg_args(bbreg_target): def create_wf_calculate_ants_warp( - name="create_wf_calculate_ants_warp", num_threads=1, reg_ants_skull=1 -): + name: Optional[str] = "create_wf_calculate_ants_warp", + num_threads: int = 1, + reg_ants_skull: int = 1, +) -> pe.Workflow: """Calculate the nonlinear ANTS registration transform. This workflow employs the antsRegistration tool: http://stnava.github.io/ANTs/ - Parameters ---------- - name : string, optional + name Name of the workflow. Returns ------- - calc_ants_warp_wf : nipype.pipeline.engine.Workflow + calc_ants_warp_wf Notes ----- @@ -1303,149 +994,87 @@ def create_wf_calculate_ants_warp( calculate_ants_warp.interface.num_threads = num_threads - select_forward_initial = pe.Node( - Function( - input_names=["warp_list", "selection"], - output_names=["selected_warp"], - function=seperate_warps_list, - ), - name="select_forward_initial", - ) - - select_forward_initial.inputs.selection = "Initial" - - select_forward_rigid = pe.Node( - Function( - input_names=["warp_list", "selection"], - output_names=["selected_warp"], - function=seperate_warps_list, - ), - name="select_forward_rigid", - ) - - select_forward_rigid.inputs.selection = "Rigid" - - select_forward_affine = pe.Node( - Function( - input_names=["warp_list", "selection"], - output_names=["selected_warp"], - function=seperate_warps_list, - ), - name="select_forward_affine", - ) - - select_forward_affine.inputs.selection = "Affine" - - select_forward_warp = pe.Node( - Function( - input_names=["warp_list", "selection"], - output_names=["selected_warp"], - function=seperate_warps_list, - ), - name="select_forward_warp", - ) - - select_forward_warp.inputs.selection = "Warp" - - select_inverse_warp = pe.Node( - Function( - input_names=["warp_list", "selection"], - output_names=["selected_warp"], - function=seperate_warps_list, - ), - name="select_inverse_warp", - ) - - select_inverse_warp.inputs.selection = "Inverse" - - calc_ants_warp_wf.connect( - inputspec, "moving_brain", calculate_ants_warp, "moving_brain" - ) + select_forward = {} + for selection in ["Initial", "Rigid", "Affine", "Warp", "Inverse"]: + select_forward[selection] = pe.Node( + Function( + input_names=["warp_list", "selection"], + output_names=["selected_warp"], + function=seperate_warps_list, + ), + name=f"select_forward_{selection.lower()}", + ) + select_forward[selection].inputs.selection = selection calc_ants_warp_wf.connect( - inputspec, "reference_brain", calculate_ants_warp, "reference_brain" + [ + ( + inputspec, + calculate_ants_warp, + [ + ("moving_brain", "moving_brain"), + ("reference_brain", "reference_brain"), + ("fixed_image_mask", "fixed_image_mask"), + ("reference_mask", "reference_mask"), + ("moving_mask", "moving_mask"), + ("ants_para", "ants_para"), + ("interp", "interp"), + ], + ) + ] ) if reg_ants_skull == 1: calculate_ants_warp.inputs.reg_with_skull = 1 - calc_ants_warp_wf.connect( - inputspec, "moving_skull", calculate_ants_warp, "moving_skull" - ) - - calc_ants_warp_wf.connect( - inputspec, "reference_skull", calculate_ants_warp, "reference_skull" + [ + ( + inputspec, + calculate_ants_warp, + [ + ("moving_skull", "moving_skull"), + ( + "reference_skull", + "reference_skull", + ), + ], + ) + ] ) - else: calc_ants_warp_wf.connect( - inputspec, "moving_brain", calculate_ants_warp, "moving_skull" - ) - - calc_ants_warp_wf.connect( - inputspec, "reference_brain", calculate_ants_warp, "reference_skull" + [ + ( + inputspec, + calculate_ants_warp, + [ + ("moving_brain", "moving_skull"), + ("reference_brain", "reference_skull"), + ], + ) + ] ) - calc_ants_warp_wf.connect( - inputspec, "fixed_image_mask", calculate_ants_warp, "fixed_image_mask" - ) - - calc_ants_warp_wf.connect( - inputspec, "reference_mask", calculate_ants_warp, "reference_mask" - ) - - calc_ants_warp_wf.connect( - inputspec, "moving_mask", calculate_ants_warp, "moving_mask" - ) - - calc_ants_warp_wf.connect(inputspec, "ants_para", calculate_ants_warp, "ants_para") - - calc_ants_warp_wf.connect(inputspec, "interp", calculate_ants_warp, "interp") - # inter-workflow connections - - calc_ants_warp_wf.connect( - calculate_ants_warp, "warp_list", select_forward_initial, "warp_list" - ) - - calc_ants_warp_wf.connect( - calculate_ants_warp, "warp_list", select_forward_rigid, "warp_list" - ) - - calc_ants_warp_wf.connect( - calculate_ants_warp, "warp_list", select_forward_affine, "warp_list" - ) - - calc_ants_warp_wf.connect( - calculate_ants_warp, "warp_list", select_forward_warp, "warp_list" - ) - - calc_ants_warp_wf.connect( - calculate_ants_warp, "warp_list", select_inverse_warp, "warp_list" - ) + for select in select_forward.values(): + calc_ants_warp_wf.connect(calculate_ants_warp, "warp_list", select, "warp_list") # connections to outputspec - calc_ants_warp_wf.connect( - select_forward_initial, "selected_warp", outputspec, "ants_initial_xfm" + select_forward["Initial"], "selected_warp", outputspec, "ants_initial_xfm" ) - calc_ants_warp_wf.connect( - select_forward_rigid, "selected_warp", outputspec, "ants_rigid_xfm" + select_forward["Rigid"], "selected_warp", outputspec, "ants_rigid_xfm" ) - calc_ants_warp_wf.connect( - select_forward_affine, "selected_warp", outputspec, "ants_affine_xfm" + select_forward["Affine"], "selected_warp", outputspec, "ants_affine_xfm" ) - calc_ants_warp_wf.connect( - select_forward_warp, "selected_warp", outputspec, "warp_field" + select_forward["Warp"], "selected_warp", outputspec, "warp_field" ) - calc_ants_warp_wf.connect( - select_inverse_warp, "selected_warp", outputspec, "inverse_warp_field" + select_forward["Inverse"], "selected_warp", outputspec, "inverse_warp_field" ) - calc_ants_warp_wf.connect( calculate_ants_warp, "warped_image", outputspec, "normalized_output_brain" ) @@ -1454,8 +1083,13 @@ def create_wf_calculate_ants_warp( def FSL_registration_connector( - wf_name, cfg, orig="T1w", opt=None, symmetric=False, template="T1w" -): + wf_name: str, + cfg: "Configuration", + orig: REGISTRATION_SPACE = "T1w", + opt: Optional[Literal["FSL", "FSL-linear"]] = None, + symmetric: bool = False, + template: REGISTRATION_SPACE = "T1w", +) -> NODEBLOCK_RETURN: """Transform raw data to template with FSL.""" wf = pe.Workflow(name=wf_name) @@ -1474,102 +1108,110 @@ def FSL_registration_connector( ] ), name="inputspec", + fnirt_config=cfg[ + "registration_workflows", + "functional_registration", + "EPI_registration", + "FSL-FNIRT", + "fnirt_config", + ], + interpolation=cfg[ + "registration_workflows", + "anatomical_registration", + "registration", + "FSL-FNIRT", + "interpolation", + ], ) - sym = "" - symm = "" - if symmetric: - sym = "sym" - symm = "_symmetric" + sym, symm, tmpl, template = prep_reg_connector(symmetric, template) - tmpl = "" - if template == "EPI": - tmpl = "EPI" + flirt_reg_anat_mni = create_fsl_flirt_linear_reg(f"anat_mni_flirt_register{symm}") - if opt in ("FSL", "FSL-linear"): - flirt_reg_anat_mni = create_fsl_flirt_linear_reg( - f"anat_mni_flirt_register{symm}" - ) + # Input registration parameters + wf.connect( + [ + ( + inputNode, + flirt_reg_anat_mni, + [ + ("interpolation", "inputspec.interp"), + ("input_brain", "inputspec.input_brain"), + ("reference_brain", "inputspec.reference_brain"), + ], + ) + ] + ) - # Input registration parameters - wf.connect(inputNode, "interpolation", flirt_reg_anat_mni, "inputspec.interp") + write_lin_composite_xfm = pe.Node( + interface=fsl.ConvertWarp(), name=f"fsl_lin-warp_to_nii{symm}" + ) - wf.connect( - inputNode, "input_brain", flirt_reg_anat_mni, "inputspec.input_brain" - ) + wf.connect(inputNode, "reference_brain", write_lin_composite_xfm, "reference") + wf.connect( + flirt_reg_anat_mni, + "outputspec.linear_xfm", + write_lin_composite_xfm, + "premat", + ) - wf.connect( - inputNode, - "reference_brain", - flirt_reg_anat_mni, - "inputspec.reference_brain", - ) + write_invlin_composite_xfm = pe.Node( + interface=fsl.ConvertWarp(), name=f"fsl_invlin-warp_to_nii{symm}" + ) - write_lin_composite_xfm = pe.Node( - interface=fsl.ConvertWarp(), name=f"fsl_lin-warp_to_nii{symm}" - ) + wf.connect(inputNode, "reference_brain", write_invlin_composite_xfm, "reference") - wf.connect(inputNode, "reference_brain", write_lin_composite_xfm, "reference") + wf.connect( + flirt_reg_anat_mni, + "outputspec.invlinear_xfm", + write_invlin_composite_xfm, + "premat", + ) - wf.connect( + outputs = { + f"space-{sym}{template}_desc-preproc_{orig}": ( flirt_reg_anat_mni, - "outputspec.linear_xfm", + "outputspec.output_brain", + ), + f"from-{orig}_to-{sym}{tmpl}{template}_mode-image_desc-linear_xfm": ( write_lin_composite_xfm, - "premat", - ) - - write_invlin_composite_xfm = pe.Node( - interface=fsl.ConvertWarp(), name=f"fsl_invlin-warp_to_nii{symm}" - ) - - wf.connect( - inputNode, "reference_brain", write_invlin_composite_xfm, "reference" - ) - - wf.connect( - flirt_reg_anat_mni, - "outputspec.invlinear_xfm", + "out_file", + ), + f"from-{sym}{tmpl}{template}_to-{orig}_mode-image_desc-linear_xfm": ( write_invlin_composite_xfm, - "premat", - ) - - outputs = { - f"space-{sym}template_desc-preproc_{orig}": ( - flirt_reg_anat_mni, - "outputspec.output_brain", - ), - f"from-{orig}_to-{sym}{tmpl}template_mode-image_desc-linear_xfm": ( - write_lin_composite_xfm, - "out_file", - ), - f"from-{sym}{tmpl}template_to-{orig}_mode-image_desc-linear_xfm": ( - write_invlin_composite_xfm, - "out_file", - ), - f"from-{orig}_to-{sym}{tmpl}template_mode-image_xfm": ( - write_lin_composite_xfm, - "out_file", - ), - } + "out_file", + ), + f"from-{orig}_to-{sym}{tmpl}{template}_mode-image_xfm": ( + write_lin_composite_xfm, + "out_file", + ), + } - if opt == "FSL": + if opt == "FSL": # as opposed to "FSL-linear" fnirt_reg_anat_mni = create_fsl_fnirt_nonlinear_reg_nhp( f"anat_mni_fnirt_register{symm}" ) wf.connect( - inputNode, "input_brain", fnirt_reg_anat_mni, "inputspec.input_brain" - ) - - wf.connect( - inputNode, - "reference_brain", - fnirt_reg_anat_mni, - "inputspec.reference_brain", + [ + ( + inputNode, + fnirt_reg_anat_mni, + [ + ("input_brain", "inputspec.input_brain"), + ("reference_brain", "inputspec.reference_brain"), + ("input_head", "inputspec.input_skull"), + ("reference_head", "inputspec.reference_skull"), + ("reference_mask", "inputspec.ref_mask"), + ( + "fnirt_config", + "inputspec.fnirt_config", + ), # assign the FSL FNIRT config file specified in pipeline config.yml + ], + ) + ] ) - wf.connect(inputNode, "input_head", fnirt_reg_anat_mni, "inputspec.input_skull") - # NOTE: crossover from above opt block wf.connect( flirt_reg_anat_mni, @@ -1578,54 +1220,46 @@ def FSL_registration_connector( "inputspec.linear_aff", ) - wf.connect( - inputNode, "reference_head", fnirt_reg_anat_mni, "inputspec.reference_skull" - ) - - wf.connect( - inputNode, "reference_mask", fnirt_reg_anat_mni, "inputspec.ref_mask" - ) - - # assign the FSL FNIRT config file specified in pipeline config.yml - wf.connect( - inputNode, "fnirt_config", fnirt_reg_anat_mni, "inputspec.fnirt_config" - ) - # NOTE: this is an UPDATE because of the opt block above added_outputs = { - f"space-{sym}template_desc-preproc_{orig}": ( + f"space-{sym}{template}_desc-preproc_{orig}": ( fnirt_reg_anat_mni, "outputspec.output_brain", ), - f"space-{sym}template_desc-head_{orig}": ( + f"space-{sym}{template}_desc-head_{orig}": ( fnirt_reg_anat_mni, "outputspec.output_head", ), - f"space-{sym}template_desc-{'brain' if orig == 'T1w' else orig}_mask": ( + f"space-{sym}{template}_desc-{'brain' if orig == 'T1w' else orig}_mask": ( fnirt_reg_anat_mni, "outputspec.output_mask", ), - f"space-{sym}template_desc-T1wT2w_biasfield": ( + f"space-{sym}{template}_desc-T1wT2w_biasfield": ( fnirt_reg_anat_mni, "outputspec.output_biasfield", ), - f"from-{orig}_to-{sym}{tmpl}template_mode-image_xfm": ( + f"from-{orig}_to-{sym}{tmpl}{template}_mode-image_xfm": ( fnirt_reg_anat_mni, "outputspec.nonlinear_xfm", ), - f"from-{orig}_to-{sym}{tmpl}template_mode-image_warp": ( + f"from-{orig}_to-{sym}{tmpl}{template}_mode-image_warp": ( fnirt_reg_anat_mni, "outputspec.nonlinear_warp", ), } outputs.update(added_outputs) - return (wf, outputs) + return wf, outputs def ANTs_registration_connector( - wf_name, cfg, params, orig="T1w", symmetric=False, template="T1w" -): + wf_name: str, + cfg: "Configuration", + params: list[dict[str, Any]], + orig: REGISTRATION_SPACE = "T1w", + symmetric: bool = False, + template: REGISTRATION_SPACE = "T1w", +) -> NODEBLOCK_RETURN: """Transform raw data to template with ANTs.""" wf = pe.Workflow(name=wf_name) @@ -1643,17 +1277,15 @@ def ANTs_registration_connector( ] ), name="inputspec", + interpolation=cfg[ + "registration_workflows", + "anatomical_registration", + "registration", + "ANTs", + "interpolation", + ], ) - - sym = "" - symm = "" - if symmetric: - sym = "sym" - symm = "_symmetric" - - tmpl = "" - if template == "EPI": - tmpl = "EPI" + sym, symm, tmpl, template = prep_reg_connector(symmetric, template) if params is None: err_msg = ( @@ -1672,31 +1304,27 @@ def ANTs_registration_connector( ) ants_reg_anat_mni.inputs.inputspec.ants_para = params - wf.connect(inputNode, "interpolation", ants_reg_anat_mni, "inputspec.interp") - # calculating the transform with the skullstripped is # reported to be better, but it requires very high # quality skullstripping. If skullstripping is imprecise # registration with skull is preferred - - wf.connect(inputNode, "input_brain", ants_reg_anat_mni, "inputspec.moving_brain") - - wf.connect( - inputNode, "reference_brain", ants_reg_anat_mni, "inputspec.reference_brain" - ) - - wf.connect(inputNode, "input_head", ants_reg_anat_mni, "inputspec.moving_skull") - - wf.connect( - inputNode, "reference_head", ants_reg_anat_mni, "inputspec.reference_skull" - ) - - wf.connect(inputNode, "input_mask", ants_reg_anat_mni, "inputspec.moving_mask") - wf.connect( - inputNode, "reference_mask", ants_reg_anat_mni, "inputspec.reference_mask" + [ + ( + inputNode, + ants_reg_anat_mni, + [ + ("interpolation", "inputspec.interp"), + ("input_brain", "inputspec.moving_brain"), + ("reference_brain", "inputspec.reference_brain"), + ("input_head", "inputspec.moving_skull"), + ("reference_head", "inputspec.reference_skull"), + ("input_mask", "inputspec.moving_mask"), + ("reference_mask", "inputspec.reference_mask"), + ], + ) + ] ) - ants_reg_anat_mni.inputs.inputspec.fixed_image_mask = None if orig == "T1w": @@ -1714,341 +1342,121 @@ def ANTs_registration_connector( ) # combine the linear xfm's into one - makes it easier downstream - write_composite_linear_xfm = pe.Node( - interface=ants.ApplyTransforms(), - name=f"write_composite_linear{symm}_xfm", - mem_gb=1.155, - mem_x=(1708448960473801 / 1208925819614629174706176, "input_image"), - ) - write_composite_linear_xfm.inputs.print_out_composite_warp_file = True - write_composite_linear_xfm.inputs.output_image = ( - f"from-{orig}_to-{sym}{tmpl}template_mode-image_desc-linear_xfm.nii.gz" - ) - - wf.connect(inputNode, "input_brain", write_composite_linear_xfm, "input_image") - - wf.connect( - inputNode, "reference_brain", write_composite_linear_xfm, "reference_image" - ) - - wf.connect(inputNode, "interpolation", write_composite_linear_xfm, "interpolation") - - write_composite_linear_xfm.inputs.input_image_type = 0 - write_composite_linear_xfm.inputs.dimension = 3 - - collect_transforms = pe.Node( - util.Merge(3), - name=f"collect_transforms{symm}", - mem_gb=0.8, - mem_x=(263474863123069 / 37778931862957161709568, "in1"), - ) - - wf.connect( - ants_reg_anat_mni, "outputspec.ants_affine_xfm", collect_transforms, "in1" - ) - - wf.connect( - ants_reg_anat_mni, "outputspec.ants_rigid_xfm", collect_transforms, "in2" - ) - - wf.connect( - ants_reg_anat_mni, "outputspec.ants_initial_xfm", collect_transforms, "in3" - ) - - # check transform list to exclude Nonetype (missing) init/rig/affine - check_transform = pe.Node( - Function( - input_names=["transform_list"], - output_names=["checked_transform_list", "list_length"], - function=check_transforms, - ), - name="check_transforms", - mem_gb=6, - ) - - wf.connect(collect_transforms, "out", check_transform, "transform_list") - - wf.connect( - check_transform, - "checked_transform_list", - write_composite_linear_xfm, - "transforms", + write_composite_linear_xfm = compose_ants_warp( + wf=wf, + name=f"linear{sym}_xfm", + input_node=inputNode, + warp_from="input_brain", + warp_to="reference_brain", + inputs=[ + ( + ants_reg_anat_mni, + [ + "outputspec.ants_affine_xfm", + "outputspec.ants_rigid_xfm", + "outputspec.ants_initial_xfm", + ], + ) + ], ) - # combine the linear xfm's into one - makes it easier downstream - write_composite_invlinear_xfm = pe.Node( - interface=ants.ApplyTransforms(), - name=f"write_composite_invlinear{symm}_xfm", + # combine the inverse linear xfm's into one - makes it easier downstream + write_composite_invlinear_xfm = compose_ants_warp( + wf=wf, + name=f"invlinear{sym}_xfm", + input_node=inputNode, + warp_from="reference_brain", + warp_to="input_brain", + inputs=[ + ( + ants_reg_anat_mni, + [ + "outputspec.ants_initial_xfm", + "outputspec.ants_rigid_xfm", + "outputspec.ants_affine_xfm", + ], + ) + ], + inv=True, mem_gb=1.05, mem_x=(1367826948979337 / 151115727451828646838272, "input_image"), ) - write_composite_invlinear_xfm.inputs.print_out_composite_warp_file = True - write_composite_invlinear_xfm.inputs.output_image = ( - f"from-{sym}{tmpl}template_to-{orig}_mode-image_desc-linear_xfm.nii.gz" - ) - - wf.connect( - inputNode, "reference_brain", write_composite_invlinear_xfm, "input_image" - ) - - wf.connect( - inputNode, "input_brain", write_composite_invlinear_xfm, "reference_image" - ) - - wf.connect( - inputNode, "interpolation", write_composite_invlinear_xfm, "interpolation" - ) - - write_composite_invlinear_xfm.inputs.input_image_type = 0 - write_composite_invlinear_xfm.inputs.dimension = 3 - - collect_inv_transforms = pe.Node( - util.Merge(3), name=f"collect_inv_transforms{symm}" - ) - - wf.connect( - ants_reg_anat_mni, "outputspec.ants_initial_xfm", collect_inv_transforms, "in1" - ) - - wf.connect( - ants_reg_anat_mni, "outputspec.ants_rigid_xfm", collect_inv_transforms, "in2" - ) - - wf.connect( - ants_reg_anat_mni, "outputspec.ants_affine_xfm", collect_inv_transforms, "in3" - ) - - # check transform list to exclude Nonetype (missing) init/rig/affine - check_invlinear_transform = pe.Node( - Function( - input_names=["transform_list"], - output_names=["checked_transform_list", "list_length"], - function=check_transforms, - ), - name="check_inv_transforms", - ) - - wf.connect( - collect_inv_transforms, "out", check_invlinear_transform, "transform_list" - ) - - wf.connect( - check_invlinear_transform, - "checked_transform_list", - write_composite_invlinear_xfm, - "transforms", - ) - - # generate inverse transform flags, which depends on the - # number of transforms - inverse_transform_flags = pe.Node( - Function( - input_names=["transform_list"], - output_names=["inverse_transform_flags"], - function=generate_inverse_transform_flags, - ), - name="inverse_transform_flags", - ) - - wf.connect( - check_invlinear_transform, - "checked_transform_list", - inverse_transform_flags, - "transform_list", - ) - - wf.connect( - inverse_transform_flags, - "inverse_transform_flags", - write_composite_invlinear_xfm, - "invert_transform_flags", - ) # combine ALL xfm's into one - makes it easier downstream - write_composite_xfm = pe.Node( - interface=ants.ApplyTransforms(), name=f"write_composite_{symm}xfm", mem_gb=1.5 - ) - write_composite_xfm.inputs.print_out_composite_warp_file = True - write_composite_xfm.inputs.output_image = ( - f"from-{orig}_to-{sym}{tmpl}template_mode-image_xfm.nii.gz" - ) - - wf.connect(inputNode, "input_brain", write_composite_xfm, "input_image") - - wf.connect(inputNode, "reference_brain", write_composite_xfm, "reference_image") - - wf.connect(inputNode, "interpolation", write_composite_xfm, "interpolation") - - write_composite_xfm.inputs.input_image_type = 0 - write_composite_xfm.inputs.dimension = 3 - - collect_all_transforms = pe.Node( - util.Merge(4), name=f"collect_all_transforms{symm}" - ) - - wf.connect( - ants_reg_anat_mni, "outputspec.warp_field", collect_all_transforms, "in1" - ) - - wf.connect( - ants_reg_anat_mni, "outputspec.ants_affine_xfm", collect_all_transforms, "in2" - ) - - wf.connect( - ants_reg_anat_mni, "outputspec.ants_rigid_xfm", collect_all_transforms, "in3" - ) - - wf.connect( - ants_reg_anat_mni, "outputspec.ants_initial_xfm", collect_all_transforms, "in4" - ) - - # check transform list to exclude Nonetype (missing) init/rig/affine - check_all_transform = pe.Node( - Function( - input_names=["transform_list"], - output_names=["checked_transform_list", "list_length"], - function=check_transforms, - ), - name="check_all_transforms", - ) - - wf.connect(collect_all_transforms, "out", check_all_transform, "transform_list") - - wf.connect( - check_all_transform, "checked_transform_list", write_composite_xfm, "transforms" + write_composite_xfm = compose_ants_warp( + wf=wf, + name=f"write_composite_{sym}_xfm", + input_node=inputNode, + warp_from="input_brain", + warp_to="reference_brain", + inputs=[ + ( + ants_reg_anat_mni, + [ + "outputspec.warp_field", + "outputspec.ants_affine_xfm", + "outputspec.ants_rigid_xfm", + "outputspec.ants_initial_xfm", + ], + ) + ], ) # combine ALL xfm's into one - makes it easier downstream - write_composite_inv_xfm = pe.Node( - interface=ants.ApplyTransforms(), - name=f"write_composite_inv_{symm}xfm", + write_composite_inv_xfm = compose_ants_warp( + wf=wf, + name=f"write_composite_inv_{sym}_xfm", + input_node=inputNode, + warp_from="reference_brain", + warp_to="input_brain", + inputs=[ + ( + ants_reg_anat_mni, + [ + "outputspec.ants_initial_xfm", + "outputspec.ants_rigid_xfm", + "outputspec.ants_affine_xfm", + "outputspec.inverse_warp_field", + ], + ) + ], + inv=True, mem_gb=0.3, mem_x=(6278549929741219 / 604462909807314587353088, "input_image"), ) - write_composite_inv_xfm.inputs.print_out_composite_warp_file = True - write_composite_inv_xfm.inputs.output_image = ( - f"from-{sym}{tmpl}template_to-{orig}_mode-image_xfm.nii.gz" - ) - - wf.connect(inputNode, "reference_brain", write_composite_inv_xfm, "input_image") - - wf.connect(inputNode, "input_brain", write_composite_inv_xfm, "reference_image") - - wf.connect(inputNode, "interpolation", write_composite_inv_xfm, "interpolation") - - write_composite_inv_xfm.inputs.input_image_type = 0 - write_composite_inv_xfm.inputs.dimension = 3 - - collect_all_inv_transforms = pe.Node( - util.Merge(4), name=f"collect_all_inv_transforms{symm}" - ) - - wf.connect( - ants_reg_anat_mni, - "outputspec.ants_initial_xfm", - collect_all_inv_transforms, - "in1", - ) - - wf.connect( - ants_reg_anat_mni, - "outputspec.ants_rigid_xfm", - collect_all_inv_transforms, - "in2", - ) - - wf.connect( - ants_reg_anat_mni, - "outputspec.ants_affine_xfm", - collect_all_inv_transforms, - "in3", - ) - - wf.connect( - ants_reg_anat_mni, - "outputspec.inverse_warp_field", - collect_all_inv_transforms, - "in4", - ) - - # check transform list to exclude Nonetype (missing) init/rig/affine - check_all_inv_transform = pe.Node( - Function( - input_names=["transform_list"], - output_names=["checked_transform_list", "list_length"], - function=check_transforms, - ), - name="check_all_inv_transforms", - ) - - wf.connect( - collect_all_inv_transforms, "out", check_all_inv_transform, "transform_list" - ) - - wf.connect( - check_all_inv_transform, - "checked_transform_list", - write_composite_inv_xfm, - "transforms", - ) - - # generate inverse transform flags, which depends on the - # number of transforms - inverse_all_transform_flags = pe.Node( - Function( - input_names=["transform_list"], - output_names=["inverse_transform_flags"], - function=generate_inverse_transform_flags, - ), - name="inverse_all_transform_flags", - ) - - wf.connect( - check_all_inv_transform, - "checked_transform_list", - inverse_all_transform_flags, - "transform_list", - ) - - wf.connect( - inverse_all_transform_flags, - "inverse_transform_flags", - write_composite_inv_xfm, - "invert_transform_flags", - ) outputs = { - f"space-{sym}template_desc-preproc_{orig}": ( + f"space-{sym}{template}_desc-preproc_T1w": ( ants_reg_anat_mni, "outputspec.normalized_output_brain", ), - f"from-{orig}_to-{sym}{tmpl}template_mode-image_xfm": ( + f"from-{orig}_to-{sym}{tmpl}{template}_mode-image_xfm": ( write_composite_xfm, "output_image", ), - f"from-{sym}{tmpl}template_to-{orig}_mode-image_xfm": ( + f"from-{sym}{tmpl}{template}_to-{orig}_mode-image_xfm": ( write_composite_inv_xfm, "output_image", ), - f"from-{orig}_to-{sym}{tmpl}template_mode-image_desc-linear_xfm": ( + f"from-{orig}_to-{sym}{tmpl}{template}_mode-image_desc-linear_xfm": ( write_composite_linear_xfm, "output_image", ), - f"from-{sym}{tmpl}template_to-{orig}_mode-image_desc-linear_xfm": ( + f"from-{sym}{tmpl}{template}_to-{orig}_mode-image_desc-linear_xfm": ( write_composite_invlinear_xfm, "output_image", ), - f"from-{orig}_to-{sym}{tmpl}template_mode-image_desc-nonlinear_xfm": ( + f"from-{orig}_to-{sym}{tmpl}{template}_mode-image_desc-nonlinear_xfm": ( ants_reg_anat_mni, "outputspec.warp_field", ), - f"from-{sym}{tmpl}template_to-{orig}_mode-image_desc-nonlinear_xfm": ( + f"from-{sym}{tmpl}{template}_to-{orig}_mode-image_desc-nonlinear_xfm": ( ants_reg_anat_mni, "outputspec.inverse_warp_field", ), } - return (wf, outputs) + return wf, outputs def bold_to_T1template_xfm_connector( @@ -2227,8 +1635,14 @@ def bold_to_T1template_xfm_connector( option_val=["FSL", "FSL-linear"], inputs=[ ( - ["desc-preproc_T1w", "space-longitudinal_desc-reorient_T1w"], - ["desc-brain_T1w", "space-longitudinal_desc-brain_T1w"], + [ + "desc-preproc_T1w", + "longitudinal-template_space-longitudinal_desc-head_T1w", + ], + [ + "desc-brain_T1w", + "longitudinal-template_space-longitudinal_desc-brain_T1w", + ], ), "T1w-template", "T1w-brain-template", @@ -2239,19 +1653,16 @@ def bold_to_T1template_xfm_connector( outputs={ "space-template_desc-preproc_T1w": {"Template": "T1w-brain-template"}, "space-template_desc-head_T1w": {"Template": "T1w-template"}, - "space-template_desc-brain_mask": {"Template": "T1w-template"}, "space-template_desc-T1wT2w_biasfield": {"Template": "T1w-template"}, - "from-T1w_to-template_mode-image_desc-linear_xfm": {"Template": "T1w-template"}, - "from-template_to-T1w_mode-image_desc-linear_xfm": {"Template": "T1w-template"}, - "from-T1w_to-template_mode-image_xfm": {"Template": "T1w-template"}, "from-T1w_to-template_mode-image_warp": {"Template": "T1w-template"}, - "from-longitudinal_to-template_mode-image_desc-linear_xfm": { - "Template": "T1w-template" - }, - "from-template_to-longitudinal_mode-image_desc-linear_xfm": { - "Template": "T1w-template" - }, - "from-longitudinal_to-template_mode-image_xfm": {"Template": "T1w-template"}, + **xfm_outputs( + spaces={ + "longitudinal": "longitudinal-T1w-template", + "template": "T1w-template", + "T1w": "native T1w", + }, + template="template", + ), }, ) def register_FSL_anat_to_template(wf, cfg, strat_pool, pipe_num, opt=None): @@ -2260,16 +1671,9 @@ def register_FSL_anat_to_template(wf, cfg, strat_pool, pipe_num, opt=None): f"register_{opt}_anat_to_template_{pipe_num}", cfg, orig="T1w", opt=opt ) - fsl.inputs.inputspec.interpolation = cfg.registration_workflows[ - "anatomical_registration" - ]["registration"]["FSL-FNIRT"]["interpolation"] - - fsl.inputs.inputspec.fnirt_config = cfg.registration_workflows[ - "anatomical_registration" - ]["registration"]["FSL-FNIRT"]["fnirt_config"] - connect, brain = strat_pool.get_data( - ["desc-brain_T1w", "space-longitudinal_desc-brain_T1w"], report_fetched=True + ["desc-brain_T1w", "longitudinal-template_space-longitudinal_desc-brain_T1w"], + report_fetched=True, ) node, out = connect wf.connect(node, out, fsl, "inputspec.input_brain") @@ -2281,7 +1685,10 @@ def register_FSL_anat_to_template(wf, cfg, strat_pool, pipe_num, opt=None): wf.connect(node, out, fsl, "inputspec.reference_head") node, out = strat_pool.get_data( - ["desc-preproc_T1w", "space-longitudinal_desc-reorient_T1w"] + [ + "desc-preproc_T1w", + "longitudinal-template_space-longitudinal_desc-head_T1w", + ] ) wf.connect(node, out, fsl, "inputspec.input_head") @@ -2289,7 +1696,9 @@ def register_FSL_anat_to_template(wf, cfg, strat_pool, pipe_num, opt=None): wf.connect(node, out, fsl, "inputspec.reference_mask") if "space-longitudinal" in brain: - for key in outputs.keys(): + for key in list( # `list` for a copy, to make outputs mutable in this loop + outputs.keys() + ): if "from-T1w" in key: new_key = key.replace("from-T1w", "from-longitudinal") outputs[new_key] = outputs[key] @@ -2299,7 +1708,7 @@ def register_FSL_anat_to_template(wf, cfg, strat_pool, pipe_num, opt=None): outputs[new_key] = outputs[key] del outputs[key] - return (wf, outputs) + return wf, outputs @nodeblock( @@ -2310,8 +1719,14 @@ def register_FSL_anat_to_template(wf, cfg, strat_pool, pipe_num, opt=None): option_val=["FSL", "FSL-linear"], inputs=[ ( - ["desc-preproc_T1w", "space-longitudinal_desc-reorient_T1w"], - ["desc-brain_T1w", "space-longitudinal_desc-brain_T1w"], + [ + "desc-preproc_T1w", + "longitudinal-template_space-longitudinal_desc-head_T1w", + ], + [ + "desc-brain_T1w", + "longitudinal-template_space-longitudinal_desc-brain_T1w", + ], ), "T1w-template-symmetric", "T1w-brain-template-symmetric", @@ -2327,20 +1742,14 @@ def register_FSL_anat_to_template(wf, cfg, strat_pool, pipe_num, opt=None): "brain_mask", ] }, - **{ - output: {"Template": "T1w-template-symmetric"} - for output in [ - "space-symtemplate_desc-head_T1w", - "from-T1w_to-symtemplate_mode-image_desc-linear_xfm", - "from-symtemplate_to-T1w_mode-image_desc-linear_xfm", - "from-T1w_to-symtemplate_mode-image_warp", - "from-T1w_to-symtemplate_mode-image_xfm", - "from-longitudinal_to-symtemplate_mode-image_desc-linear_xfm", - "from-symtemplate_to-longitudinal_mode-image_desc-linear_xfm", - "from-longitudinal_to-symtemplate_mode-image_xfm", - "space-symtemplate_desc-T1wT2w_biasfield", - ] - }, + **xfm_outputs( + spaces={ + "longitudinal": "longitudinal-T1w-template", + "symtemplate": "T1w-template-symmetric", + "T1w": "native T1w", + }, + template="template", + ), }, ) def register_symmetric_FSL_anat_to_template(wf, cfg, strat_pool, pipe_num, opt=None): @@ -2353,16 +1762,9 @@ def register_symmetric_FSL_anat_to_template(wf, cfg, strat_pool, pipe_num, opt=N symmetric=True, ) - fsl.inputs.inputspec.interpolation = cfg.registration_workflows[ - "anatomical_registration" - ]["registration"]["FSL-FNIRT"]["interpolation"] - - fsl.inputs.inputspec.fnirt_config = cfg.registration_workflows[ - "anatomical_registration" - ]["registration"]["FSL-FNIRT"]["fnirt_config"] - connect, brain = strat_pool.get_data( - ["desc-brain_T1w", "space-longitudinal_desc-brain_T1w"], report_fetched=True + ["desc-brain_T1w", "longitudinal-template_space-longitudinal_desc-brain_T1w"], + report_fetched=True, ) node, out = connect wf.connect(node, out, fsl, "inputspec.input_brain") @@ -2371,7 +1773,10 @@ def register_symmetric_FSL_anat_to_template(wf, cfg, strat_pool, pipe_num, opt=N wf.connect(node, out, fsl, "inputspec.reference_brain") node, out = strat_pool.get_data( - ["desc-preproc_T1w", "space-longitudinal_desc-reorient_T1w"] + [ + "desc-preproc_T1w", + "longitudinal-template_space-longitudinal_desc-head_T1w", + ] ) wf.connect(node, out, fsl, "inputspec.input_head") @@ -2382,7 +1787,9 @@ def register_symmetric_FSL_anat_to_template(wf, cfg, strat_pool, pipe_num, opt=N wf.connect(node, out, fsl, "inputspec.reference_mask") if "space-longitudinal" in brain: - for key in outputs.keys(): + for key in list( # `list` for a copy, to make outputs mutable in this loop + outputs.keys() + ): if "from-T1w" in key: new_key = key.replace("from-T1w", "from-longitudinal") outputs[new_key] = outputs[key] @@ -2427,14 +1834,6 @@ def register_FSL_EPI_to_template(wf, cfg, strat_pool, pipe_num, opt=None): template="EPI", ) - fsl.inputs.inputspec.interpolation = cfg["registration_workflows"][ - "functional_registration" - ]["EPI_registration"]["FSL-FNIRT"]["interpolation"] - - fsl.inputs.inputspec.fnirt_config = cfg["registration_workflows"][ - "functional_registration" - ]["EPI_registration"]["FSL-FNIRT"]["fnirt_config"] - node, out = strat_pool.get_data("sbref") wf.connect(node, out, fsl, "inputspec.input_brain") @@ -2461,18 +1860,22 @@ def register_FSL_EPI_to_template(wf, cfg, strat_pool, pipe_num, opt=None): option_val="ANTS", inputs=[ ( - ["desc-preproc_T1w", "space-longitudinal_desc-brain_T1w"], + "desc-preproc_T1w", + "longitudinal-template_space-longitudinal_desc-brain_T1w", [ "space-T1w_desc-brain_mask", - "space-longitudinal_desc-brain_mask", "space-T1w_desc-acpcbrain_mask", ], + "longitudinal-template_space-longitudinal_desc-brain_mask", [ "desc-restore_T1w", "desc-head_T1w", "desc-preproc_T1w", - "space-longitudinal_desc-reorient_T1w", ], + "longitudinal-template_space-longitudinal_desc-head_T1w", + "space-longitudinal_desc-head_T1w", + "space-longitudinal_desc-brain_T1w", + "space-longitudinal_desc-preproc_T1w", "space-template_desc-head_T1w", "space-template_desc-preproc_T1w", ), @@ -2487,143 +1890,75 @@ def register_FSL_EPI_to_template(wf, cfg, strat_pool, pipe_num, opt=None): "template space.", "Template": "T1w-template", }, - "from-T1w_to-template_mode-image_desc-linear_xfm": { - "Description": "Linear (affine) transform from T1w native space " - "to T1w-template space.", - "Template": "T1w-template", - }, - "from-template_to-T1w_mode-image_desc-linear_xfm": { - "Description": "Linear (affine) transform from T1w-template space " - "to T1w native space.", - "Template": "T1w-template", - }, - "from-T1w_to-template_mode-image_desc-nonlinear_xfm": { - "Description": "Nonlinear (warp field) transform from T1w native " - "space to T1w-template space.", - "Template": "T1w-template", - }, - "from-template_to-T1w_mode-image_desc-nonlinear_xfm": { - "Description": "Nonlinear (warp field) transform from " - "T1w-template space to T1w native space.", - "Template": "T1w-template", - }, - "from-T1w_to-template_mode-image_xfm": { - "Description": "Composite (affine + warp field) transform from " - "T1w native space to T1w-template space.", - "Template": "T1w-template", - }, - "from-template_to-T1w_mode-image_xfm": { - "Description": "Composite (affine + warp field) transform from " - "T1w-template space to T1w native space.", - "Template": "T1w-template", - }, - "from-longitudinal_to-template_mode-image_desc-linear_xfm": { - "Description": "Linear (affine) transform from " - "longitudinal-template space to T1w-template " - "space.", - "Template": "T1w-template", - }, - "from-template_to-longitudinal_mode-image_desc-linear_xfm": { - "Description": "Linear (affine) transform from T1w-template " - "space to longitudinal-template space.", - "Template": "T1w-template", - }, - "from-longitudinal_to-template_mode-image_desc-nonlinear_xfm": { - "Description": "Nonlinear (warp field) transform from " - "longitudinal-template space to T1w-template " - "space.", - "Template": "T1w-template", - }, - "from-template_to-longitudinal_mode-image_desc-nonlinear_xfm": { - "Description": "Nonlinear (warp field) transform from " - "T1w-template space to longitudinal-template " - "space.", - "Template": "T1w-template", - }, - "from-longitudinal_to-template_mode-image_xfm": { - "Description": "Composite (affine + warp field) transform from " - "longitudinal-template space to T1w-template " - "space.", - "Template": "T1w-template", - }, - "from-template_to-longitudinal_mode-image_xfm": { - "Description": "Composite (affine + warp field) transform from " - "T1w-template space to longitudinal-template " - "space.", - "Template": "T1w-template", - }, + **xfm_outputs( + spaces={ + "longitudinal": "longitudinal-T1w-template", + "template": "T1w-template", + "T1w": "native T1w", + }, + template="template", + ), }, ) -def register_ANTs_anat_to_template(wf, cfg, strat_pool, pipe_num, opt=None): +def register_ANTs_anat_to_template( + wf: pe.Workflow, + cfg: "Configuration", + strat_pool: "ResourcePool", + pipe_num: int, + opt: Optional[str] = None, +) -> NODEBLOCK_RETURN: """Register T1w to template with ANTs.""" - params = cfg.registration_workflows["anatomical_registration"]["registration"][ - "ANTs" - ]["T1_registration"] - - ants_rc, outputs = ANTs_registration_connector( - f"ANTS_T1_to_template_{pipe_num}", cfg, params, orig="T1w" - ) - - ants_rc.inputs.inputspec.interpolation = cfg.registration_workflows[ - "anatomical_registration" - ]["registration"]["ANTs"]["interpolation"] - - connect, brain = strat_pool.get_data( - ["desc-preproc_T1w", "space-longitudinal_desc-brain_T1w"], report_fetched=True + params = cast( + list[dict[str, Any]], + cfg[ + "registration_workflows", + "anatomical_registration", + "registration", + "ANTs", + "T1_registration", + ], ) - node, out = connect - wf.connect(node, out, ants_rc, "inputspec.input_brain") + ( + orig, + has_longitudinal, + input_brain, + input_head, + reference_mask, + lesion_mask, + t1w_brain_template, + t1w_template, + brain_mask, + ) = get_common_reg_inputs(strat_pool) + _rc_params = { + "wf_name": f"ANTS_T1_to_template_{pipe_num}", + "cfg": cfg, + "params": params, + "orig": orig, + } + if has_longitudinal: + _rc_params["wf_name"] = f"ANTS_T1_to_longitudinal_{pipe_num}" + _rc_params["template"] = "longitudinal" + ants_rc, outputs = ANTs_registration_connector(**_rc_params) - t1w_brain_template = strat_pool.node_data("T1w-brain-template") + wf.connect(input_brain.node, input_brain.out, ants_rc, "inputspec.input_brain") wf.connect( t1w_brain_template.node, t1w_brain_template.out, ants_rc, "inputspec.reference_brain", ) - - # TODO check the order of T1w - node, out = strat_pool.get_data( - [ - "desc-restore_T1w", - "desc-head_T1w", - "desc-preproc_T1w", - "space-longitudinal_desc-reorient_T1w", - ] - ) - wf.connect(node, out, ants_rc, "inputspec.input_head") - - t1w_template = strat_pool.node_data("T1w-template") + wf.connect(input_head.node, input_head.out, ants_rc, "inputspec.input_head") wf.connect(t1w_template.node, t1w_template.out, ants_rc, "inputspec.reference_head") - - brain_mask = strat_pool.node_data( - [ - "space-T1w_desc-brain_mask", - "space-longitudinal_desc-brain_mask", - "space-T1w_desc-acpcbrain_mask", - ] - ) wf.connect(brain_mask.node, brain_mask.out, ants_rc, "inputspec.input_mask") - - if strat_pool.check_rpool("T1w-brain-template-mask"): - node, out = strat_pool.get_data("T1w-brain-template-mask") - wf.connect(node, out, ants_rc, "inputspec.reference_mask") - - if strat_pool.check_rpool("label-lesion_mask"): - node, out = strat_pool.get_data("label-lesion_mask") - wf.connect(node, out, ants_rc, "inputspec.lesion_mask") - - if "space-longitudinal" in brain: - for key in outputs: - for direction in ["from", "to"]: - if f"{direction}-T1w" in key: - new_key = key.replace( - f"{direction}-T1w", f"{direction}-longitudinal" - ) - outputs[new_key] = outputs[key] - del outputs[key] - - return (wf, outputs) + if reference_mask: + wf.connect( + reference_mask.node, reference_mask.out, ants_rc, "inputspec.reference_mask" + ) + if lesion_mask: + wf.connect(lesion_mask.node, lesion_mask.out, ants_rc, "inputspec.lesion_mask") + if has_longitudinal: + outputs.update(t1w_to_longitudinal_to_template["ants"](wf, strat_pool)) + return wf, outputs @nodeblock( @@ -2634,13 +1969,19 @@ def register_ANTs_anat_to_template(wf, cfg, strat_pool, pipe_num, opt=None): option_val="ANTS", inputs=[ ( - ["desc-preproc_T1w", "space-longitudinal_desc-brain_T1w"], - ["space-T1w_desc-brain_mask", "space-longitudinal_desc-brain_mask"], + "desc-preproc_T1w", + "longitudinal-template_space-longitudinal_desc-brain_T1w", [ + "space-T1w_desc-brain_mask", + "space-T1w_desc-acpcbrain_mask", + ], + "longitudinal-template_space-longitudinal_desc-brain_mask", + [ + "desc-restore_T1w", "desc-head_T1w", "desc-preproc_T1w", - "space-longitudinal_desc-reorient_T1w", ], + "longitudinal-template_space-longitudinal_desc-head_T1w", ), "T1w-template-symmetric", "T1w-brain-template-symmetric", @@ -2651,103 +1992,79 @@ def register_ANTs_anat_to_template(wf, cfg, strat_pool, pipe_num, opt=None): "space-symtemplate_desc-preproc_T1w": { "Template": "T1w-brain-template-symmetric" }, - "from-T1w_to-symtemplate_mode-image_desc-linear_xfm": { - "Template": "T1w-template-symmetric" - }, - "from-symtemplate_to-T1w_mode-image_desc-linear_xfm": { - "Template": "T1w-template-symmetric" - }, - "from-T1w_to-symtemplate_mode-image_desc-nonlinear_xfm": { - "Template": "T1w-template-symmetric" - }, - "from-symtemplate_to-T1w_mode-image_desc-nonlinear_xfm": { - "Template": "T1w-template-symmetric" - }, - "from-T1w_to-symtemplate_mode-image_xfm": { - "Template": "T1w-template-symmetric" - }, - "from-symtemplate_to-T1w_mode-image_xfm": { - "Template": "T1w-template-symmetric" - }, - "from-longitudinal_to-symtemplate_mode-image_desc-linear_xfm": { - "Template": "T1w-template-symmetric" - }, - "from-symtemplate_to-longitudinal_mode-image_desc-linear_xfm": { - "Template": "T1w-template-symmetric" - }, - "from-longitudinal_to-symtemplate_mode-image_desc-nonlinear_xfm": { - "Template": "T1w-template-symmetric" - }, - "from-symtemplate_to-longitudinal_mode-image_desc-nonlinear_xfm": { - "Template": "T1w-template-symmetric" - }, - "from-longitudinal_to-symtemplate_mode-image_xfm": { - "Template": "T1w-template-symmetric" - }, - "from-symtemplate_to-longitudinal_mode-image_xfm": { - "Template": "T1w-template-symmetric" - }, + **xfm_outputs( + spaces={ + "longitudinal": "longitudinal-T1w-template", + "symtemplate": "T1w-template-symmetric", + "T1w": "native T1w", + }, + template="symtemplate", + ), }, ) -def register_symmetric_ANTs_anat_to_template(wf, cfg, strat_pool, pipe_num, opt=None): +def register_symmetric_ANTs_anat_to_template( + wf: pe.Workflow, + cfg: "Configuration", + strat_pool: "ResourcePool", + pipe_num: int, + opt: Optional[str] = None, +) -> NODEBLOCK_RETURN: """Register T1 to symmetric template with ANTs.""" - params = cfg.registration_workflows["anatomical_registration"]["registration"][ - "ANTs" - ]["T1_registration"] - - ants, outputs = ANTs_registration_connector( - f"ANTS_T1_to_template_symmetric_{pipe_num}", - cfg, - params, - orig="T1w", - symmetric=True, - ) - - ants.inputs.inputspec.interpolation = cfg.registration_workflows[ - "anatomical_registration" - ]["registration"]["ANTs"]["interpolation"] - - connect, brain = strat_pool.get_data( - ["desc-preproc_T1w", "space-longitudinal_desc-brain_T1w"], report_fetched=True - ) - node, out = connect - wf.connect(node, out, ants, "inputspec.input_brain") - - node, out = strat_pool.get_data("T1w-brain-template-symmetric") - wf.connect(node, out, ants, "inputspec.reference_brain") - - node, out = strat_pool.get_data( - ["desc-head_T1w", "desc-preproc_T1w", "space-longitudinal_desc-reorient_T1w"] + params = cast( + list[dict[str, Any]], + cfg[ + "registration_workflows", + "anatomical_registration", + "registration", + "ANTs", + "T1_registration", + ], ) - wf.connect(node, out, ants, "inputspec.input_head") - - node, out = strat_pool.get_data("T1w-template-symmetric") - wf.connect(node, out, ants, "inputspec.reference_head") - - node, out = strat_pool.get_data( - ["space-T1w_desc-brain_mask", "space-longitudinal_desc-brain_mask"] + ( + orig, + has_longitudinal, + input_brain, + input_head, + reference_mask, + lesion_mask, + t1w_brain_template, + t1w_template, + brain_mask, + ) = get_common_reg_inputs( + strat_pool, + registration_templates=RegistrationTemplates( + reference_brain="T1w-brain-template-symmetric", + reference_head="T1w-template-symmetric", + reference_mask="dilated-symmetric-brain-mask", + ), ) - wf.connect(node, out, ants, "inputspec.input_mask") - - node, out = strat_pool.get_data("dilated-symmetric-brain-mask") - wf.connect(node, out, ants, "inputspec.reference_mask") - - if strat_pool.check_rpool("label-lesion_mask"): - node, out = strat_pool.get_data("label-lesion_mask") - wf.connect(node, out, ants, "inputspec.lesion_mask") - - if "space-longitudinal" in brain: - for key in outputs.keys(): - if "from-T1w" in key: - new_key = key.replace("from-T1w", "from-longitudinal") - outputs[new_key] = outputs[key] - del outputs[key] - if "to-T1w" in key: - new_key = key.replace("to-T1w", "to-longitudinal") - outputs[new_key] = outputs[key] - del outputs[key] - - return (wf, outputs) + _rc_params = { + "wf_name": f"ANTS_T1_to_template_symmetric_{pipe_num}", + "cfg": cfg, + "params": params, + "orig": orig, + "symmetric": True, + } + if has_longitudinal: + _rc_params["wf_name"] = f"ANTS_longitudinal_to_template_symmetric_{pipe_num}" + _rc_params["template"] = "longitudinal" + ants, outputs = ANTs_registration_connector(**_rc_params) + wf.connect(*input_brain, ants, "inputspec.input_brain") + wf.connect(*t1w_brain_template, ants, "inputspec.reference_brain") + wf.connect(*input_head, ants, "inputspec.input_head") + wf.connect(*t1w_template, ants, "inputspec.reference_head") + wf.connect(*brain_mask, ants, "inputspec.input_mask") + if reference_mask: + wf.connect(*reference_mask, ants, "inputspec.reference_mask") + if lesion_mask: + wf.connect(*lesion_mask, ants, "inputspec.lesion_mask") + if has_longitudinal: + outputs.update( + t1w_to_longitudinal_to_template["ants"]( + wf, strat_pool, template="symtemplate" + ) + ) + return wf, outputs @nodeblock( @@ -2840,7 +2157,10 @@ def register_ANTs_EPI_to_template(wf, cfg, strat_pool, pipe_num, opt=None): inputs=[ ( "desc-restore-brain_T1w", - ["desc-preproc_T1w", "space-longitudinal_desc-brain_T1w"], + [ + "desc-preproc_T1w", + "space-longitudinal_desc-brain_T1w", + ], ["desc-restore_T1w", "desc-preproc_T1w", "desc-reorient_T1w", "T1w"], ["desc-preproc_T1w", "desc-reorient_T1w", "T1w"], "space-T1w_desc-brain_mask", @@ -2860,9 +2180,7 @@ def register_ANTs_EPI_to_template(wf, cfg, strat_pool, pipe_num, opt=None): ) def overwrite_transform_anat_to_template(wf, cfg, strat_pool, pipe_num, opt=None): """Overwrite ANTs transforms with FSL transforms.""" - xfm_prov = strat_pool.get_cpac_provenance("from-T1w_to-template_mode-image_xfm") - - reg_tool = check_prov_for_regtool(xfm_prov) + reg_tool = strat_pool.reg_tool("from-T1w_to-template_mode-image_xfm") if opt.lower() == "fsl" and reg_tool.lower() == "ants": # Apply head-to-head transforms on brain using ABCD-style registration @@ -3080,7 +2398,7 @@ def overwrite_transform_anat_to_template(wf, cfg, strat_pool, pipe_num, opt=None else: outputs = {} - return (wf, outputs) + return wf, outputs @nodeblock( @@ -3233,7 +2551,7 @@ def coregistration_prep_fmriprep(wf, cfg, strat_pool, pipe_num, opt=None): else: outputs["sbref"] = strat_pool.get_data("desc-unifized_bold") - return (wf, outputs) + return wf, outputs @nodeblock( @@ -3489,8 +2807,7 @@ def create_func_to_T1template_xfm(wf, cfg, strat_pool, pipe_num, opt=None): Condense the BOLD-to-T1 coregistration transform and the T1-to-template transform into one transform matrix. """ - xfm_prov = strat_pool.get_cpac_provenance("from-T1w_to-template_mode-image_xfm") - reg_tool = check_prov_for_regtool(xfm_prov) + reg_tool = strat_pool.reg_tool("from-T1w_to-template_mode-image_xfm") xfm, outputs = bold_to_T1template_xfm_connector( f"create_func_to_T1wtemplate_xfm_{pipe_num}", cfg, reg_tool, symmetric=False @@ -3568,8 +2885,7 @@ def create_func_to_T1template_symmetric_xfm(wf, cfg, strat_pool, pipe_num, opt=N Condense the BOLD-to-T1 coregistration transform and the T1-to-symmetric-template transform into one transform matrix. """ - xfm_prov = strat_pool.get_cpac_provenance("from-T1w_to-symtemplate_mode-image_xfm") - reg_tool = check_prov_for_regtool(xfm_prov) + reg_tool = strat_pool.reg_tool("from-T1w_to-symtemplate_mode-image_xfm") xfm, outputs = bold_to_T1template_xfm_connector( f"create_func_to_T1wsymtemplate_xfm_{pipe_num}", @@ -3774,8 +3090,7 @@ def apply_phasediff_to_timeseries_separately(wf, cfg, strat_pool, pipe_num, opt= ) def apply_blip_to_timeseries_separately(wf, cfg, strat_pool, pipe_num, opt=None): """Apply blip to timeseries.""" - xfm_prov = strat_pool.get_cpac_provenance("from-bold_to-template_mode-image_xfm") - reg_tool = check_prov_for_regtool(xfm_prov) + reg_tool = strat_pool.reg_tool("from-bold_to-template_mode-image_xfm") outputs = {"desc-preproc_bold": strat_pool.get_data("desc-preproc_bold")} if strat_pool.check_rpool("ants-blip-warp"): @@ -3836,115 +3151,112 @@ def apply_blip_to_timeseries_separately(wf, cfg, strat_pool, pipe_num, opt=None) return (wf, outputs) -@nodeblock( - name="transform_whole_head_T1w_to_T1template", - config=["registration_workflows", "anatomical_registration"], - switch=["run"], - inputs=[ - ( - "desc-head_T1w", - "from-T1w_to-template_mode-image_xfm", - "space-template_desc-head_T1w", - ), - "T1w-template", - ], - outputs={"space-template_desc-head_T1w": {"Template": "T1w-template"}}, -) -def warp_wholeheadT1_to_template(wf, cfg, strat_pool, pipe_num, opt=None): - """Warp T1 head to template.""" - xfm_prov = strat_pool.get_cpac_provenance("from-T1w_to-template_mode-image_xfm") - reg_tool = check_prov_for_regtool(xfm_prov) - - num_cpus = cfg.pipeline_setup["system_config"]["max_cores_per_participant"] - - num_ants_cores = cfg.pipeline_setup["system_config"]["num_ants_threads"] - - apply_xfm = apply_transform( - f"warp_wholehead_T1w_to_T1template_{pipe_num}", - reg_tool, - time_series=False, - num_cpus=num_cpus, - num_ants_cores=num_ants_cores, - ) - - if reg_tool == "ants": - apply_xfm.inputs.inputspec.interpolation = cfg.registration_workflows[ - "functional_registration" - ]["func_registration_to_template"]["ANTs_pipelines"]["interpolation"] - elif reg_tool == "fsl": - apply_xfm.inputs.inputspec.interpolation = cfg.registration_workflows[ - "functional_registration" - ]["func_registration_to_template"]["FNIRT_pipelines"]["interpolation"] - - connect = strat_pool.get_data("desc-head_T1w") - node, out = connect - wf.connect(node, out, apply_xfm, "inputspec.input_image") - - node, out = strat_pool.get_data("T1w-template") - wf.connect(node, out, apply_xfm, "inputspec.reference") - - node, out = strat_pool.get_data("from-T1w_to-template_mode-image_xfm") - wf.connect(node, out, apply_xfm, "inputspec.transform") +def warp_to_template( + warp_what: Literal["mask", "wholehead"], space_from: Literal["longitudinal", "T1w"] +) -> NodeBlockFunction: + """Get a NodeBlockFunction to transform a resource from ``space`` to template. - outputs = {"space-template_desc-head_T1w": (apply_xfm, "outputspec.output_image")} + The resource being warped needs to be the first list or string in the tuple + in the first position of the decorator's "inputs". + """ + _decorators = { + "mask": { + "name": f"transform_{space_from}-mask_to_T1-template", + "switch": [ + ["registration_workflows", "anatomical_registration", "run"], + ["anatomical_preproc", "run"], + ["anatomical_preproc", "brain_extraction", "run"], + ], + "inputs": [ + ( + f"space-{space_from}_desc-brain_mask", + f"from-{space_from}_to-template_mode-image_xfm", + ), + "T1w-template", + ], + "outputs": {"space-template_desc-brain_mask": {"Template": "T1w-template"}}, + }, + "wholehead": { + "name": f"transform_wholehead_{space_from}_to_T1template", + "config": ["registration_workflows", "anatomical_registration"], + "switch": ["run"], + "inputs": [ + ( + ["desc-head_T1w", "desc-reorient_T1w"], + [ + f"from-{space_from}_to-template_mode-image_xfm", + f"from-{space_from}_to-template_mode-image_xfm", + ], + "space-template_desc-head_T1w", + ), + "T1w-template", + ], + "outputs": {"space-template_desc-head_T1w": {"Template": "T1w-template"}}, + }, + } + if space_from != "T1w": + _decorators[warp_what]["inputs"][0] = ( + prepend_space(_decorators[warp_what]["inputs"][0][0], space_from), + *_decorators[warp_what]["inputs"][0][1:], + ) - return (wf, outputs) + @nodeblock(**_decorators[warp_what]) + def warp_to_template_fxn( + wf: pe.Workflow, + cfg: "Configuration", + strat_pool: "ResourcePool", + pipe_num: int, + opt: Optional[str] = None, + ) -> NODEBLOCK_RETURN: + """Transform a resource to template space.""" + reg_tool = strat_pool.reg_tool(f"from-{space_from}_to-template_mode-image_xfm") + num_cpus = cfg.pipeline_setup["system_config"]["max_cores_per_participant"] -@nodeblock( - name="transform_T1mask_to_T1template", - switch=[ - ["registration_workflows", "anatomical_registration", "run"], - ["anatomical_preproc", "run"], - ["anatomical_preproc", "brain_extraction", "run"], - ], - inputs=[ - ("space-T1w_desc-brain_mask", "from-T1w_to-template_mode-image_xfm"), - "T1w-template", - ], - outputs={"space-template_desc-brain_mask": {"Template": "T1w-template"}}, -) -def warp_T1mask_to_template(wf, cfg, strat_pool, pipe_num, opt=None): - """Warp T1 mask to template.""" - xfm_prov = strat_pool.get_cpac_provenance("from-T1w_to-template_mode-image_xfm") - reg_tool = check_prov_for_regtool(xfm_prov) + num_ants_cores = cfg.pipeline_setup["system_config"]["num_ants_threads"] - num_cpus = cfg.pipeline_setup["system_config"]["max_cores_per_participant"] + apply_xfm = apply_transform( + f"warp_{space_from}{warp_what}_to_T1template_{pipe_num}", + reg_tool, + time_series=False, + num_cpus=num_cpus, + num_ants_cores=num_ants_cores, + ) - num_ants_cores = cfg.pipeline_setup["system_config"]["num_ants_threads"] + if warp_what == "mask": + apply_xfm.inputs.inputspec.interpolation = "NearestNeighbor" + else: + tool = ( + "ANTs" if reg_tool == "ants" else "FNIRT" if reg_tool == "fsl" else None + ) + if not tool: + msg = f"Warp {warp_what} to template not implemented for {reg_tool}." + raise NotImplementedError(msg) + apply_xfm.inputs.inputspec.interpolation = cfg.registration_workflows[ + "functional_registration" + ]["func_registration_to_template"][f"{tool}_pipelines"]["interpolation"] - apply_xfm = apply_transform( - f"warp_T1mask_to_T1template_{pipe_num}", - reg_tool, - time_series=False, - num_cpus=num_cpus, - num_ants_cores=num_ants_cores, - ) + # the resource being warped needs to be inputs[0][0] for this + node, out = strat_pool.get_data(_decorators[warp_what]["inputs"][0][0]) + wf.connect(node, out, apply_xfm, "inputspec.input_image") - apply_xfm.inputs.inputspec.interpolation = "NearestNeighbor" - """ - if reg_tool == 'ants': - apply_xfm.inputs.inputspec.interpolation = cfg.registration_workflows[ - 'functional_registration']['func_registration_to_template'][ - 'ANTs_pipelines']['interpolation'] - elif reg_tool == 'fsl': - apply_xfm.inputs.inputspec.interpolation = cfg.registration_workflows[ - 'functional_registration']['func_registration_to_template'][ - 'FNIRT_pipelines']['interpolation'] - """ - connect = strat_pool.get_data("space-T1w_desc-brain_mask") - node, out = connect - wf.connect(node, out, apply_xfm, "inputspec.input_image") + node, out = strat_pool.get_data("T1w-template") + wf.connect(node, out, apply_xfm, "inputspec.reference") - node, out = strat_pool.get_data("T1w-template") - wf.connect(node, out, apply_xfm, "inputspec.reference") + node, out = strat_pool.get_data(f"from-{space_from}_to-template_mode-image_xfm") + wf.connect(node, out, apply_xfm, "inputspec.transform") - node, out = strat_pool.get_data("from-T1w_to-template_mode-image_xfm") - wf.connect(node, out, apply_xfm, "inputspec.transform") + outputs = { + # there's only one output, so that's what we give here + next(iter(_decorators[warp_what]["outputs"].keys())): ( + apply_xfm, + "outputspec.output_image", + ) + } - outputs = {"space-template_desc-brain_mask": (apply_xfm, "outputspec.output_image")} + return wf, outputs - return (wf, outputs) + return warp_to_template_fxn @nodeblock( @@ -3967,8 +3279,7 @@ def warp_T1mask_to_template(wf, cfg, strat_pool, pipe_num, opt=None): ) def warp_timeseries_to_T1template(wf, cfg, strat_pool, pipe_num, opt=None): """Warp timeseries to T1 template.""" - xfm_prov = strat_pool.get_cpac_provenance("from-bold_to-template_mode-image_xfm") - reg_tool = check_prov_for_regtool(xfm_prov) + reg_tool = strat_pool.reg_tool("from-bold_to-template_mode-image_xfm") num_cpus = cfg.pipeline_setup["system_config"]["max_cores_per_participant"] @@ -4030,8 +3341,7 @@ def warp_timeseries_to_T1template(wf, cfg, strat_pool, pipe_num, opt=None): ) def warp_timeseries_to_T1template_deriv(wf, cfg, strat_pool, pipe_num, opt=None): """Warp timeseries to T1 template at derivative resolution.""" - xfm_prov = strat_pool.get_cpac_provenance("from-bold_to-template_mode-image_xfm") - reg_tool = check_prov_for_regtool(xfm_prov) + reg_tool = strat_pool.reg_tool("from-bold_to-template_mode-image_xfm") num_cpus = cfg.pipeline_setup["system_config"]["max_cores_per_participant"] @@ -4833,7 +4143,7 @@ def warp_timeseries_to_T1template_dcan_nhp(wf, cfg, strat_pool, pipe_num, opt=No }, ) def single_step_resample_timeseries_to_T1template( - wf, cfg, strat_pool, pipe_num, opt=None + wf, cfg, strat_pool: "ResourcePool", pipe_num, opt=None ): """Apply motion correction, coreg, anat-to-template transforms... @@ -4872,9 +4182,8 @@ def single_step_resample_timeseries_to_T1template( # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED # OF THE POSSIBILITY OF SUCH DAMAGE. - # Modifications copyright (C) 2021 - 2024 C-PAC Developers - xfm_prov = strat_pool.get_cpac_provenance("from-T1w_to-template_mode-image_xfm") - reg_tool = check_prov_for_regtool(xfm_prov) + # Modifications copyright (C) 2021 - 2025 C-PAC Developers + reg_tool = strat_pool.reg_tool("from-T1w_to-template_mode-image_xfm") bbr2itk = pe.Node( Function( @@ -4932,9 +4241,7 @@ def single_step_resample_timeseries_to_T1template( wf.connect(node, out, motionxfm2itk, "source_file") node, out = strat_pool.get_data("coordinate-transformation") - motion_correct_tool = check_prov_for_motion_tool( - strat_pool.get_cpac_provenance("coordinate-transformation") - ) + motion_correct_tool = strat_pool.motion_tool("coordinate-transformation") if motion_correct_tool == "mcflirt": wf.connect(node, out, motionxfm2itk, "transform_file") elif motion_correct_tool == "3dvolreg": @@ -5371,6 +4678,7 @@ def warp_deriv_mask_to_EPItemplate(wf, cfg, strat_pool, pipe_num, opt=None): "label-WM_mask", "label-GM_mask", "from-T1w_to-template_mode-image_xfm", + "from-longitudinal_to-template_mode-image_xfm", ), "T1w-template", ], @@ -5387,7 +4695,9 @@ def warp_tissuemask_to_T1template(wf, cfg, strat_pool, pipe_num, opt=None): cfg, strat_pool, pipe_num, - xfm="from-T1w_to-template_mode-image_xfm", + xfm="from-longitudinal_to-template_mode-image_xfm" + if strat_pool.check_rpool("from-longitudinal_to-template_mode-image_xfm") + else "from-T1w_to-template_mode-image_xfm", template_space="T1", ) @@ -5467,14 +4777,15 @@ def warp_tissuemask_to_template(wf, cfg, strat_pool, pipe_num, xfm, template_spa "outputspec.output_image", ) for tissue in tissue_types + if apply_xfm[tissue] } return _warp_return(wf, apply_xfm, outputs) def warp_resource_to_template( wf: pe.Workflow, - cfg, - strat_pool, + cfg: "Configuration", + strat_pool: "ResourcePool", pipe_num: int, input_resource: list[str] | str, xfm: str, @@ -5522,8 +4833,7 @@ def warp_resource_to_template( if template_space == "": template_space = "T1w" # determine tool used for registration - xfm_prov = strat_pool.get_cpac_provenance(xfm) - reg_tool = check_prov_for_regtool(xfm_prov) + reg_tool = strat_pool.reg_tool(xfm) # set 'resource' if strat_pool.check_rpool(input_resource): resource, input_resource = strat_pool.get_data( diff --git a/CPAC/registration/utils.py b/CPAC/registration/utils.py index 4e0dc4421e..90f78011b6 100644 --- a/CPAC/registration/utils.py +++ b/CPAC/registration/utils.py @@ -1,4 +1,4 @@ -# Copyright (C) 2014-2024 C-PAC Developers +# Copyright (C) 2014-2025 C-PAC Developers # This file is part of C-PAC. @@ -18,9 +18,367 @@ import os import subprocess +from typing import Literal, NamedTuple, Optional, overload, TYPE_CHECKING, TypeAlias import numpy as np from voluptuous import RequiredFieldInvalid +from nipype.interfaces.afni.utils import TCat +from nipype.interfaces.ants import ApplyTransforms as AntsApplyTransforms +from nipype.interfaces.fsl import ApplyWarp +from nipype.interfaces.utility import IdentityInterface, Merge +from nipype.pipeline.engine import Node as NipypeNode, Workflow + +from CPAC.func_preproc.utils import chunk_ts, split_ts_chunks +from CPAC.pipeline.nipype_pipeline_engine import MapNode, Node +from CPAC.utils.interfaces import Function + +if TYPE_CHECKING: + from CPAC.pipeline.engine import NodeData + +REGISTRATION_SPACE: TypeAlias = Literal[ + "bold", "EPI", "T1", "T1w", "longitudinal", "template" +] + + +class CommonRegistrationInputs(NamedTuple): + """Input nodes registration methods take in common.""" + + orig: Literal["longitudinal", "T1w"] + has_longitudinal: bool + input_brain: "NodeData" + input_head: "NodeData" + reference_mask: Optional["NodeData"] + lesion_mask: Optional["NodeData"] + t1w_brain_template: "NodeData" + t1w_template: "NodeData" + brain_mask: "NodeData" + + +class RegistrationTemplates(NamedTuple): + """Keys for registration templates in strat pool.""" + + reference_brain: list[str] | str = "T1w-brain-template" + reference_head: list[str] | str = "T1w-template" + reference_mask: list[str] | str = "T1w-brain-template-mask" + + +@overload +def convert_pedir(pedir: bytes | str, convert: Literal["xyz_to_int"]) -> int: ... +@overload +def convert_pedir(pedir: bytes | str, convert: Literal["ijk_to_xyz"]) -> str: ... +def convert_pedir( + pedir: bytes | str, convert: Literal["xyz_to_int", "ijk_to_xyz"] = "xyz_to_int" +) -> int | str: + """FSL Flirt requires pedir input encoded as an int.""" + if convert == "xyz_to_int": + conv_dct = { + "x": 1, + "y": 2, + "z": 3, + "x-": -1, + "y-": -2, + "z-": -3, + "i": 1, + "j": 2, + "k": 3, + "i-": -1, + "j-": -2, + "k-": -3, + "-x": -1, + "-i": -1, + "-y": -2, + "-j": -2, + "-z": -3, + "-k": -3, + } + elif convert == "ijk_to_xyz": + conv_dct = {"i": "x", "j": "y", "k": "z", "i-": "x-", "j-": "y-", "k-": "z-"} + + if isinstance(pedir, bytes): + pedir = pedir.decode() + if not isinstance(pedir, str): + msg = f"\n\nPhase-encoding direction must be a string value.\n\nValue: {pedir}\n\n" + raise ValueError(msg) + if pedir not in conv_dct.keys(): + msg = f"\n\nInvalid phase-encoding direction entered: {pedir}\n\n" + raise ValueError(msg) + return conv_dct[pedir] + + +def apply_transform( + wf_name: str, + reg_tool: str, + time_series: bool = False, + multi_input: bool = False, + num_cpus: int = 1, + num_ants_cores: int = 1, +): + """Apply transform.""" + if not reg_tool: + msg = ( + "\n[!] Developer info: the 'reg_tool' parameter sent to the" + f" 'apply_transform' node for '{wf_name}' is empty.\n" + ) + raise RequiredFieldInvalid(msg) + + wf = Workflow(name=wf_name) + + inputNode = Node( + IdentityInterface( + fields=["input_image", "reference", "transform", "interpolation"] + ), + name="inputspec", + ) + + outputNode = Node(IdentityInterface(fields=["output_image"]), name="outputspec") + + if int(num_cpus) > 1 and time_series: + # parallelize time series warp application + # we need the node to be a MapNode to feed in the list of functional + # time series chunks + multi_input = True + + if reg_tool == "ants": + if multi_input: + apply_warp = MapNode( + interface=AntsApplyTransforms(), + name=f"apply_warp_{wf_name}", + iterfield=["input_image"], + mem_gb=0.7, + mem_x=(1708448960473801 / 151115727451828646838272, "input_image"), + ) + else: + apply_warp = Node( + interface=AntsApplyTransforms(), + name=f"apply_warp_{wf_name}", + mem_gb=0.7, + mem_x=(1708448960473801 / 151115727451828646838272, "input_image"), + ) + + apply_warp.inputs.dimension = 3 + apply_warp.interface.num_threads = int(num_ants_cores) + + if time_series: + apply_warp.inputs.input_image_type = 3 + + wf.connect(inputNode, "reference", apply_warp, "reference_image") + + interp_string = Node( + Function( + input_names=["interpolation", "reg_tool"], + output_names=["interpolation"], + function=interpolation_string, + ), + name="interp_string", + mem_gb=2.5, + ) + interp_string.inputs.reg_tool = reg_tool + + wf.connect(inputNode, "interpolation", interp_string, "interpolation") + wf.connect(interp_string, "interpolation", apply_warp, "interpolation") + + ants_xfm_list = Node( + Function( + input_names=["transform"], + output_names=["transform_list"], + function=single_ants_xfm_to_list, + ), + name="single_ants_xfm_to_list", + mem_gb=2.5, + ) + + wf.connect(inputNode, "transform", ants_xfm_list, "transform") + wf.connect(ants_xfm_list, "transform_list", apply_warp, "transforms") + + # parallelize the apply warp, if multiple CPUs, and it's a time + # series! + if int(num_cpus) > 1 and time_series: + chunk_imports = ["import nibabel as nib"] + chunk = Node( + Function( + input_names=["func_file", "n_chunks", "chunk_size"], + output_names=["TR_ranges"], + function=chunk_ts, + imports=chunk_imports, + ), + name=f"chunk_{wf_name}", + mem_gb=2.5, + ) + + # chunk.inputs.n_chunks = int(num_cpus) + + # 10-TR sized chunks + chunk.inputs.chunk_size = 10 + + wf.connect(inputNode, "input_image", chunk, "func_file") + + split_imports = ["import os", "import subprocess"] + split = Node( + Function( + input_names=["func_file", "tr_ranges"], + output_names=["split_funcs"], + function=split_ts_chunks, + imports=split_imports, + ), + name=f"split_{wf_name}", + mem_gb=2.5, + ) + + wf.connect(inputNode, "input_image", split, "func_file") + wf.connect(chunk, "TR_ranges", split, "tr_ranges") + + wf.connect(split, "split_funcs", apply_warp, "input_image") + + func_concat = Node( + interface=TCat(), name=f"func_concat_{wf_name}", mem_gb=2.5 + ) + func_concat.inputs.outputtype = "NIFTI_GZ" + + wf.connect(apply_warp, "output_image", func_concat, "in_files") + + wf.connect(func_concat, "out_file", outputNode, "output_image") + + else: + wf.connect(inputNode, "input_image", apply_warp, "input_image") + wf.connect(apply_warp, "output_image", outputNode, "output_image") + + elif reg_tool == "fsl": + if multi_input: + apply_warp = MapNode( + interface=ApplyWarp(), + name="fsl_apply_warp", + iterfield=["in_file"], + mem_gb=2.5, + ) + else: + apply_warp = Node(interface=ApplyWarp(), name="fsl_apply_warp", mem_gb=2.5) + + interp_string = Node( + Function( + input_names=["interpolation", "reg_tool"], + output_names=["interpolation"], + function=interpolation_string, + ), + name="interp_string", + mem_gb=2.5, + ) + interp_string.inputs.reg_tool = reg_tool + + wf.connect(inputNode, "interpolation", interp_string, "interpolation") + wf.connect(interp_string, "interpolation", apply_warp, "interp") + + # mni to t1 + wf.connect(inputNode, "reference", apply_warp, "ref_file") + + # NOTE: C-PAC now converts all FSL xfm's to .nii, so even if the + # inputNode 'transform' is a linear xfm, it's a .nii and must + # go in as a warpfield file + wf.connect(inputNode, "transform", apply_warp, "field_file") + + # parallelize the apply warp, if multiple CPUs, and it's a time + # series! + if int(num_cpus) > 1 and time_series: + chunk_imports = ["import nibabel as nib"] + chunk = Node( + Function( + input_names=["func_file", "n_chunks", "chunk_size"], + output_names=["TR_ranges"], + function=chunk_ts, + imports=chunk_imports, + ), + name=f"chunk_{wf_name}", + mem_gb=2.5, + ) + + # chunk.inputs.n_chunks = int(num_cpus) + + # 10-TR sized chunks + chunk.inputs.chunk_size = 10 + + wf.connect(inputNode, "input_image", chunk, "func_file") + + split_imports = ["import os", "import subprocess"] + split = Node( + Function( + input_names=["func_file", "tr_ranges"], + output_names=["split_funcs"], + function=split_ts_chunks, + imports=split_imports, + ), + name=f"split_{wf_name}", + mem_gb=2.5, + ) + + wf.connect(inputNode, "input_image", split, "func_file") + wf.connect(chunk, "TR_ranges", split, "tr_ranges") + + wf.connect(split, "split_funcs", apply_warp, "in_file") + + func_concat = Node(interface=TCat(), name=f"func_concat{wf_name}") + func_concat.inputs.outputtype = "NIFTI_GZ" + + wf.connect(apply_warp, "out_file", func_concat, "in_files") + + wf.connect(func_concat, "out_file", outputNode, "output_image") + + else: + wf.connect(inputNode, "input_image", apply_warp, "in_file") + wf.connect(apply_warp, "out_file", outputNode, "output_image") + + return wf + + +def transform_derivative( + wf_name: str, + label: str, + reg_tool: Optional[str], + num_cpus: int, + num_ants_cores: int, + ants_interp: Optional[str] = None, + fsl_interp: Optional[str] = None, +): + """Transform output derivatives to template space. + + This function is designed for use with the NodeBlock connection engine. + """ + wf = Workflow(name=wf_name) + + inputnode = Node( + IdentityInterface(fields=["in_file", "reference", "transform"]), + name="inputspec", + ) + + multi_input = False + if "statmap" in label: + multi_input = True + + stack = False + if "correlations" in label: + stack = True + + apply_xfm = apply_transform( + f"warp_{label}_to_template", + reg_tool, + time_series=stack, + multi_input=multi_input, + num_cpus=num_cpus, + num_ants_cores=num_ants_cores, + ) + + if reg_tool == "ants": + apply_xfm.inputs.inputspec.interpolation = ants_interp + elif reg_tool == "fsl": + apply_xfm.inputs.inputspec.interpolation = fsl_interp + + wf.connect(inputnode, "in_file", apply_xfm, "inputspec.input_image") + wf.connect(inputnode, "reference", apply_xfm, "inputspec.reference") + wf.connect(inputnode, "transform", apply_xfm, "inputspec.transform") + + outputnode = Node(IdentityInterface(fields=["out_file"]), name="outputspec") + + wf.connect(apply_xfm, "outputspec.output_image", outputnode, "out_file") + + return wf def single_ants_xfm_to_list(transform): @@ -808,3 +1166,191 @@ def run_c4d(input_name, output_name): os.system(cmd) return output1, output2, output3 + + +@overload +def prepend_space(resource: list[str], space: str) -> list[str]: ... +@overload +def prepend_space(resource: str, space: str) -> str: ... +def prepend_space(resource: str | list[str], space: str) -> str | list[str]: + """Given a resource or list of resources, return same but with updated space.""" + if isinstance(resource, list): + return [prepend_space(_, space) for _ in resource] + prefix = "longitudinal-template_" if space == "longitudinal" else "" + if "space" not in resource: + return f"{prefix}space-{space}_{resource}" + pre, post = resource.split("space-") + _old_space, post = post.split("_", 1) + return f"{prefix}space-{space}_".join([pre, post]) + + +def collect_xfms( + wf: Workflow, + name: str, + node_inputs: list[tuple[Node | Workflow, list[str]]], + *, + mem_gb: float = 0.8, + mem_x: tuple[float, str] = (263474863123069 / 37778931862957161709568, "in1"), +) -> Node: + """Create a node to collect transforms to compose into one.""" + if len(node_inputs) == 1: + input_node, inputs = node_inputs[0] + else: + msg = "Combining transforms from multiple nodes not yet implemented." + raise NotImplementedError(msg) + numinputs = len(inputs) + node = Node(interface=Merge(numinputs), name=name, mem_gb=mem_gb, mem_x=mem_x) + # check transform list to exclude Nonetype (missing) init/rig/affine + check_transform = Node( + Function( + input_names=["transform_list"], + output_names=["checked_transform_list", "list_length"], + function=check_transforms, + ), + name=f"check_transforms_{name}", + mem_gb=6, + ) + wf.connect( + [ + ( + input_node, + node, + [(node_input, f"in{i + 1}") for i, node_input in enumerate(inputs)], + ), + (node, check_transform, [("out", "transform_list")]), + ] + ) + return check_transform + + +def compose_ants_warp( # noqa: PLR0913 + wf: Workflow, + name: str, + input_node: NipypeNode | Workflow, + warp_from: str, + warp_to: str, + inputs: list[tuple[Node | Workflow, list[str]]], + inv: bool = False, + *, + dimension: int = 3, + input_image_type: int = 0, + mem_gb: float = 1.155, + mem_x: tuple[float, str] = ( + 1708448960473801 / 1208925819614629174706176, + "input_image", + ), +) -> Node: + """Create a Node to combine xfms.""" + node = Node( + interface=AntsApplyTransforms( + dimension=dimension, + input_image_type=input_image_type, + print_out_composite_warp_file=True, + output_image=f"{name}.nii.gz", + ), + name=f"write_composite_{name}", + mem_gb=mem_gb, + mem_x=mem_x, + ) + wf.connect( + [ + ( + input_node, + node, + [ + (warp_from, "input_image"), + (warp_to, "reference_image"), + ("interpolation", "interpolation"), + ], + ) + ] + ) + collect_transforms = collect_xfms(wf, f"collect_{name}", inputs) + wf.connect(collect_transforms, "checked_transform_list", node, "transforms") + if inv: + # generate inverse transform flags, which depends on the + # number of transforms + inverse_transform_flags = Node( + Function( + input_names=["transform_list"], + output_names=["inverse_transform_flags"], + function=generate_inverse_transform_flags, + ), + name=f"inverse_transform_flags_{name}", + ) + wf.connect( + collect_transforms, + "checked_transform_list", + inverse_transform_flags, + "transform_list", + ) + wf.connect( + inverse_transform_flags, + "inverse_transform_flags", + node, + "invert_transform_flags", + ) + return node + + +def prep_reg_connector( + symmetric: bool, template: REGISTRATION_SPACE +) -> tuple[ + Literal["", "sym"], + Literal["", "_symmetric"], + Literal["", "EPI"], + REGISTRATION_SPACE, +]: + """Return some formatted strings. + + Returns + ------- + sym + String to indicate if symmetric in resource names + + symm + String to indicate if symmetric in node names + + tmpl + String to indicate EPI template space in resource names + + template + String to indicate template space in resource names + """ + sym = "" + symm = "" + if symmetric: + sym = "sym" + symm = "_symmetric" + + tmpl = "" + match template: + case "EPI": + tmpl = "EPI" + template = "template" + case "longitudinal": + template = template # noqa: PLW0127 + case _: + template = "template" + return sym, symm, tmpl, template + + +def xfm_outputs(spaces: dict[str, str], template: str) -> dict[str, dict[str, str]]: + """Build dictionary for XFM output specs.""" + transform_types = { + "": "Composite (affine + warp field)", + "_desc-linear": "Linear (affine)", + "_desc-nonlinear": "Nonlinear (warp field)", + } + return { + f"from-{origin}_to-{destination}_mode-image{transform_type}_xfm": { + "Description": f"{transform_type_desc} transform from {origin_desc} space to {destination_desc} space.", + "Template": spaces.get( + template, destination_desc if destination != "T1w" else origin_desc + ), + } + for origin, origin_desc in spaces.items() + for destination, destination_desc in spaces.items() + for transform_type, transform_type_desc in transform_types.items() + if origin != destination + } diff --git a/CPAC/resources/configs/1.7-1.8-nesting-mappings.yml b/CPAC/resources/configs/1.7-1.8-nesting-mappings.yml index dd83685bc1..8ad79125b0 100644 --- a/CPAC/resources/configs/1.7-1.8-nesting-mappings.yml +++ b/CPAC/resources/configs/1.7-1.8-nesting-mappings.yml @@ -322,12 +322,15 @@ longitudinal_template_dof: - dof longitudinal_template_interp: - longitudinal_template_generation + - legacy-specific - interp longitudinal_template_cost: - longitudinal_template_generation + - legacy-specific - cost longitudinal_template_thread_pool: - longitudinal_template_generation + - legacy-specific - thread_pool longitudinal_template_convergence_threshold: - longitudinal_template_generation diff --git a/CPAC/resources/configs/pipeline_config_blank.yml b/CPAC/resources/configs/pipeline_config_blank.yml index 59442c47dc..d629534e5b 100644 --- a/CPAC/resources/configs/pipeline_config_blank.yml +++ b/CPAC/resources/configs/pipeline_config_blank.yml @@ -1510,29 +1510,43 @@ longitudinal_template_generation: # at once. run: Off - # Freesurfer longitudinal template algorithm using FSL FLIRT + # Implementation to use + # Options: mri_robust_template, C-PAC legacy + using: mri_robust_template + # Method to average the dataset at each iteration of the template creation - # Options: median, mean or std + # Options: median, mean + # Additional option if using "C-PAC legacy": std average_method: median # Degree of freedom for FLIRT in the template creation - # Options: 12 (affine), 9 (traditional), 7 (global rescale) or 6 (rigid body) + # Options: 12 (affine) or 6 (rigid body) + # Additional options if using "C-PAC legacy": 9 (traditional), 7 (global rescale) dof: 12 - # Interpolation parameter for FLIRT in the template creation - # Options: trilinear, nearestneighbour, sinc or spline - interp: trilinear + # Maximum iterations + # Stop after this many iterations, even if still above convergence_threshold + # Additional option if using "mri_robust_template": "default" means 5 for 2 sessions, 6 for more than 2 sessions + # Additional option if using "C-PAC legacy": -1 means loop forever until reaching convergence threshold + max_iter: default + + # Options for C-PAC legacy implementation that are not configurable in mri_robust_template + legacy-specific: + + # Threshold of transformation distance to consider that the loop converged + # (-1 means numpy.finfo(np.float64).eps and is the default) + convergence_threshold: -1 - # Cost function for FLIRT in the template creation - # Options: corratio, mutualinfo, normmi, normcorr, leastsq, labeldiff or bbr - cost: corratio + # Interpolation parameter for FLIRT in the template creation + # Options: trilinear, nearestneighbour, sinc or spline + interp: - # Number of threads used for one run of the template generation algorithm - thread_pool: 2 + # Cost function for FLIRT in the template creation + # Options: corratio, mutualinfo, normmi, normcorr, leastsq, labeldiff or bbr + cost: - # Threshold of transformation distance to consider that the loop converged - # (-1 means numpy.finfo(np.float64).eps and is the default) - convergence_threshold: -1 + # Number of threads used for one run of the template generation algorithm + thread_pool: # OUTPUTS AND DERIVATIVES # ----------------------- diff --git a/CPAC/resources/configs/pipeline_config_default.yml b/CPAC/resources/configs/pipeline_config_default.yml index 070bce1196..16d3f7af83 100644 --- a/CPAC/resources/configs/pipeline_config_default.yml +++ b/CPAC/resources/configs/pipeline_config_default.yml @@ -257,30 +257,43 @@ longitudinal_template_generation: # at once. run: Off - # Freesurfer longitudinal template algorithm using FSL FLIRT + # Implementation to use + # Options: mri_robust_template, C-PAC legacy + using: mri_robust_template + # Method to average the dataset at each iteration of the template creation - # Options: median, mean or std + # Options: median, mean + # Additional option if using "C-PAC legacy": std average_method: median # Degree of freedom for FLIRT in the template creation - # Options: 12 (affine), 9 (traditional), 7 (global rescale) or 6 (rigid body) + # Options: 12 (affine) or 6 (rigid body) + # Additional options if using "C-PAC legacy": 9 (traditional), 7 (global rescale) dof: 12 - # Interpolation parameter for FLIRT in the template creation - # Options: trilinear, nearestneighbour, sinc or spline - interp: trilinear + # Maximum iterations + # Stop after this many iterations, even if still above convergence_threshold + # Additional option if using "mri_robust_template": "default" means 5 for 2 sessions, 6 for more than 2 sessions + # Additional option if using "C-PAC legacy": -1 means loop forever until reaching convergence threshold + max_iter: 6 + + # Options for C-PAC legacy implementation that are not configurable in mri_robust_template + legacy-specific: - # Cost function for FLIRT in the template creation - # Options: corratio, mutualinfo, normmi, normcorr, leastsq, labeldiff or bbr - cost: corratio + # Threshold of transformation distance to consider that the loop converged + # (-1 means numpy.finfo(np.float64).eps and is the default) + convergence_threshold: -1 - # Number of threads used for one run of the template generation algorithm - thread_pool: 2 + # Interpolation parameter for FLIRT in the template creation + # Options: trilinear, nearestneighbour, sinc or spline + interp: trilinear - # Threshold of transformation distance to consider that the loop converged - # (-1 means numpy.finfo(np.float64).eps and is the default) - convergence_threshold: -1 + # Cost function for FLIRT in the template creation + # Options: corratio, mutualinfo, normmi, normcorr, leastsq, labeldiff or bbr + cost: corratio + # Number of threads used for one run of the template generation algorithm + thread_pool: 2 anatomical_preproc: diff --git a/CPAC/resources/cpac_outputs.tsv b/CPAC/resources/cpac_outputs.tsv index 8fe4cd284f..5b9e0a208b 100644 --- a/CPAC/resources/cpac_outputs.tsv +++ b/CPAC/resources/cpac_outputs.tsv @@ -76,6 +76,11 @@ space-bold_label-GM_desc-eroded_mask mask functional func NIfTI space-bold_label-GM_mask mask functional func NIfTI space-bold_label-WM_desc-eroded_mask mask functional func NIfTI space-bold_label-WM_mask mask functional func NIfTI +longitudinal-template_space-longitudinal_desc-brain_T1w T1w longitudinal T1w anat NIfTI +longitudinal-template_space-longitudinal_desc-head_T1w T1w longitudinal T1w anat NIfTI +longitudinal-template_space-longitudinal_desc-brain_mask mask longitudinal T1w anat NIfTI +space-longitudinal_desc-brain_T1w T1w longitudinal T1w anat NIfTI +space-longitudinal_desc-head_T1w T1w longitudinal T1w anat NIfTI space-longitudinal_desc-brain_mask mask longitudinal T1w anat NIfTI space-longitudinal_label-CSF_desc-preproc_mask mask longitudinal T1w anat NIfTI space-longitudinal_label-CSF_mask mask longitudinal T1w anat NIfTI @@ -186,6 +191,7 @@ from-EPItemplate_to-bold_mode-image_desc-nonlinear_xfm xfm func NIfTI from-longitudinal_to-symtemplate_mode-image_desc-linear_xfm xfm anat NIfTI from-longitudinal_to-symtemplate_mode-image_desc-nonlinear_xfm xfm anat NIfTI from-longitudinal_to-symtemplate_mode-image_xfm xfm anat NIfTI +from-longitudinal_to-T1w_mode-image_desc-linear_xfm xfm anat NIfTI from-longitudinal_to-template_mode-image_desc-linear_xfm xfm anat NIfTI from-longitudinal_to-template_mode-image_desc-nonlinear_xfm xfm anat NIfTI from-longitudinal_to-template_mode-image_xfm xfm anat NIfTI @@ -196,6 +202,7 @@ from-symtemplate_to-longitudinal_mode-image_xfm xfm anat NIfTI from-symtemplate_to-T1w_mode-image_desc-linear_xfm xfm anat NIfTI from-symtemplate_to-T1w_mode-image_desc-nonlinear_xfm xfm anat NIfTI from-symtemplate_to-T1w_mode-image_xfm xfm anat NIfTI +from-T1w_to-longitudinal_mode-image_desc-linear_xfm xfm anat NIfTI from-T1w_to-symtemplate_mode-image_desc-linear_xfm xfm anat NIfTI from-T1w_to-symtemplate_mode-image_desc-nonlinear_xfm xfm anat NIfTI from-T1w_to-symtemplate_mode-image_xfm xfm anat NIfTI @@ -258,7 +265,6 @@ hemi-R_space-native_veryinflated surface_derived func GIFTI surf hemi-L_space-fsLR_den-164k_midthickness surface_derived func GIFTI surf hemi-R_space-fsLR_den-164k_midthickness surface_derived func GIFTI surf hemi-L_space-fsLR_den-32k_midthickness surface_derived func GIFTI surf -hemi-L_space-fsLR_den-32k_midthickness surface_derived func GIFTI surf hemi-L_space-native_midthickness surface_derived func GIFTI surf hemi-R_space-native_midthickness surface_derived func GIFTI surf hemi-L_space-fsLR_den-32k_pial surface_derived func GIFTI surf diff --git a/CPAC/seg_preproc/seg_preproc.py b/CPAC/seg_preproc/seg_preproc.py index f769cf14b3..472f6d308e 100644 --- a/CPAC/seg_preproc/seg_preproc.py +++ b/CPAC/seg_preproc/seg_preproc.py @@ -1,4 +1,4 @@ -# Copyright (C) 2012-2023 C-PAC Developers +# Copyright (C) 2012-2025 C-PAC Developers # This file is part of C-PAC. @@ -20,8 +20,11 @@ from CPAC.anat_preproc.utils import mri_convert from CPAC.pipeline import nipype_pipeline_engine as pe from CPAC.pipeline.nodeblock import nodeblock -from CPAC.registration.registration import apply_transform -from CPAC.registration.utils import check_transforms, generate_inverse_transform_flags +from CPAC.registration.utils import ( + apply_transform, + check_transforms, + generate_inverse_transform_flags, +) from CPAC.seg_preproc.utils import ( check_if_file_is_empty, hardcoded_antsJointLabelFusion, @@ -35,7 +38,6 @@ from CPAC.utils.interfaces.function.seg_preproc import ( pick_tissue_from_labels_file_interface, ) -from CPAC.utils.utils import check_prov_for_regtool def process_segment_map(wf_name, use_priors, use_custom_threshold, reg_tool): @@ -495,8 +497,14 @@ def create_seg_preproc_antsJointLabel_method(wf_name="seg_preproc_templated_base option_val="FSL-FAST", inputs=[ ( - ["desc-brain_T1w", "space-longitudinal_desc-brain_T1w"], - ["space-T1w_desc-brain_mask", "space-longitudinal_desc-brain_mask"], + [ + "desc-brain_T1w", + "space-longitudinal_desc-brain_T1w", + ], + [ + "space-T1w_desc-brain_mask", + "longitudinal-template_space-longitudinal_desc-brain_mask", + ], [ "from-template_to-T1w_mode-image_desc-linear_xfm", "from-template_to-longitudinal_mode-image_desc-linear_xfm", @@ -507,27 +515,10 @@ def create_seg_preproc_antsJointLabel_method(wf_name="seg_preproc_templated_base "WM-path", ], outputs=[ - "label-CSF_mask", - "label-GM_mask", - "label-WM_mask", - "label-CSF_desc-preproc_mask", - "label-GM_desc-preproc_mask", - "label-WM_desc-preproc_mask", - "label-CSF_probseg", - "label-GM_probseg", - "label-WM_probseg", - "label-CSF_pveseg", - "label-GM_pveseg", - "label-WM_pveseg", - "space-longitudinal_label-CSF_mask", - "space-longitudinal_label-GM_mask", - "space-longitudinal_label-WM_mask", - "space-longitudinal_label-CSF_desc-preproc_mask", - "space-longitudinal_label-GM_desc-preproc_mask", - "space-longitudinal_label-WM_desc-preproc_mask", - "space-longitudinal_label-CSF_probseg", - "space-longitudinal_label-GM_probseg", - "space-longitudinal_label-WM_probseg", + f"{long}label-{tissue}_{entity}" + for long in ["", "space-longitudinal_"] + for tissue in ["CSF", "GM", "WM"] + for entity in ["mask", "desc-preproc_mask", "probseg", "pveseg"] ], ) def tissue_seg_fsl_fast(wf, cfg, strat_pool, pipe_num, opt=None): @@ -536,7 +527,6 @@ def tissue_seg_fsl_fast(wf, cfg, strat_pool, pipe_num, opt=None): # triggered by 'segments' boolean input (-g or --segments) # 'probability_maps' output is a list of individual probability maps # triggered by 'probability_maps' boolean input (-p) - segment = pe.Node( interface=fsl.FAST(), name=f"segment_{pipe_num}", @@ -574,7 +564,8 @@ def tissue_seg_fsl_fast(wf, cfg, strat_pool, pipe_num, opt=None): ) connect, resource = strat_pool.get_data( - ["desc-brain_T1w", "space-longitudinal_desc-brain_T1w"], report_fetched=True + ["desc-brain_T1w", "space-longitudinal_desc-brain_T1w"], + report_fetched=True, ) node, out = connect wf.connect(node, out, segment, "in_files") @@ -596,10 +587,9 @@ def tissue_seg_fsl_fast(wf, cfg, strat_pool, pipe_num, opt=None): xfm = "from-template_to-T1w_mode-image_desc-linear_xfm" if "space-longitudinal" in resource: xfm = "from-template_to-longitudinal_mode-image_desc-linear_xfm" - xfm_prov = strat_pool.get_cpac_provenance(xfm) - reg_tool = check_prov_for_regtool(xfm_prov) + + reg_tool = strat_pool.reg_tool(xfm) else: - xfm_prov = None reg_tool = None xfm = None @@ -662,7 +652,10 @@ def tissue_seg_fsl_fast(wf, cfg, strat_pool, pipe_num, opt=None): wf.connect(node, out, process_wm, "inputspec.brain") node, out = strat_pool.get_data( - ["space-T1w_desc-brain_mask", "space-longitudinal_desc-brain_mask"] + [ + "space-T1w_desc-brain_mask", + "longitudinal-template_space-longitudinal_desc-brain_mask", + ] ) wf.connect(node, out, process_csf, "inputspec.brain_mask") wf.connect(node, out, process_gm, "inputspec.brain_mask") @@ -752,10 +745,7 @@ def tissue_seg_fsl_fast(wf, cfg, strat_pool, pipe_num, opt=None): outputs=["label-CSF_mask", "label-GM_mask", "label-WM_mask"], ) def tissue_seg_T1_template_based(wf, cfg, strat_pool, pipe_num, opt=None): - xfm_prov = strat_pool.get_cpac_provenance( - "from-template_to-T1w_mode-image_desc-linear_xfm" - ) - reg_tool = check_prov_for_regtool(xfm_prov) + reg_tool = strat_pool.reg_tool("from-template_to-T1w_mode-image_desc-linear_xfm") use_ants = reg_tool == "ants" csf_template2t1 = tissue_mask_template_to_t1(f"CSF_{pipe_num}", use_ants) @@ -806,10 +796,9 @@ def tissue_seg_T1_template_based(wf, cfg, strat_pool, pipe_num, opt=None): ], ) def tissue_seg_EPI_template_based(wf, cfg, strat_pool, pipe_num, opt=None): - xfm_prov = strat_pool.get_cpac_provenance( + reg_tool = strat_pool.reg_tool( "from-EPItemplate_to-bold_mode-image_desc-linear_xfm" ) - reg_tool = check_prov_for_regtool(xfm_prov) use_ants = reg_tool == "ants" csf_template2t1 = tissue_mask_template_to_t1("CSF", use_ants) diff --git a/CPAC/unet/tests/test_torch.py b/CPAC/unet/tests/test_torch.py index 30195e5e36..b2a99c1090 100644 --- a/CPAC/unet/tests/test_torch.py +++ b/CPAC/unet/tests/test_torch.py @@ -45,13 +45,11 @@ def test_import_torch(monkeypatch, readonly, tmp_path, workdir): @pytest.mark.parametrize("error", [ImportError, ModuleNotFoundError, None]) def test_validate_unet(error): - """Test that pipeline validation throws error if torch is not - installable. - """ + """Test that pipeline validation throws error if torch is not installable.""" if error: import_module = MagicMock(side_effect=error()) - with patch("importlib.import_module", import_module): - with pytest.raises(OSError) as os_error: + with pytest.raises((OSError, error)) as os_error: + with patch("importlib.import_module", import_module): from CPAC.utils.configuration import Preconfiguration monkey = Preconfiguration("monkey") diff --git a/CPAC/utils/configuration/__init__.py b/CPAC/utils/configuration/__init__.py index b4a69671ce..8a07d2b2aa 100644 --- a/CPAC/utils/configuration/__init__.py +++ b/CPAC/utils/configuration/__init__.py @@ -20,6 +20,7 @@ from .configuration import ( check_pname, Configuration, + NestedKeyMixin, preconfig_yaml, Preconfiguration, set_subject, @@ -30,6 +31,7 @@ "Configuration", "configuration", "diff", + "NestedKeyMixin", "Preconfiguration", "preconfig_yaml", "set_subject", diff --git a/CPAC/utils/configuration/configuration.py b/CPAC/utils/configuration/configuration.py index 23904d8e40..5fdcacce47 100644 --- a/CPAC/utils/configuration/configuration.py +++ b/CPAC/utils/configuration/configuration.py @@ -1,4 +1,4 @@ -# Copyright (C) 2022-2024 C-PAC Developers +# Copyright (C) 2022-2025 C-PAC Developers # This file is part of C-PAC. @@ -16,18 +16,21 @@ # License along with C-PAC. If not, see . """C-PAC Configuration class and related functions.""" +from collections.abc import Iterable, KeysView +from importlib.resources import files import os import re -from typing import Any, Optional +from typing import Any, cast, Literal, Optional, overload from warnings import warn from click import BadParameter -import pkg_resources as p import yaml -from .diff import dct_diff +from CPAC.pipeline.nipype_pipeline_engine import MapNode, Node +from .diff import dct_diff, DiffDict CONFIG_KEY_TYPE = str | list[str] +_DICT = dict SPECIAL_REPLACEMENT_STRINGS = {r"${resolution_for_anat}", r"${func_resolution}"} @@ -42,7 +45,343 @@ def __init__(self): super().__init__() -class Configuration: +class NestedKeyMixin: + """Provide methods for getting and setting nested keys.""" + + def dict(self) -> dict[Any, Any]: + """Show contents as a dict.""" + return {k: v for k, v in self.__dict__.items() if not callable(v)} + + def __contains__(self, item: str | list[Any]) -> bool: + """Check if an item is in the Configuration.""" + if isinstance(item, str): + return item in self.keys() + try: + self.get_nested(self, item) + return True + except KeyError: + return False + + def __getitem__(self, key: Iterable) -> Any: + """Get an item from a nested dictionary.""" + self._check_keys(key) + if isinstance(key, str): + return getattr(self, key) + if isinstance(key, (list, tuple)): + return self.get_nested(self, key) + self.key_type_error(key) + return None + + def __setitem__(self, key: Iterable, value: Any) -> None: + """Set an item in a nested dictionary.""" + self._check_keys(key) + if isinstance(key, str): + setattr(self, key, value) + elif isinstance(key, (list, tuple)): + self.set_nested(self, key, value) + else: + self.key_type_error(key) + + def __sub__(self: "NestedKeyMixin", other: "NestedKeyMixin") -> DiffDict: + """Return the set difference between two nested dictionaries. + + Examples + -------- + >>> diff = (Preconfiguration('fmriprep-options') + ... - Preconfiguration('default')) + >>> diff['pipeline_setup']['pipeline_name'] + ('cpac_fmriprep-options', 'cpac-default-pipeline') + >>> diff['pipeline_setup']['pipeline_name'].s_value + 'cpac_fmriprep-options' + >>> diff['pipeline_setup']['pipeline_name'].t_value + 'cpac-default-pipeline' + >>> diff.s_value['pipeline_setup']['pipeline_name'] + 'cpac_fmriprep-options' + >>> diff.t_value['pipeline_setup']['pipeline_name'] + 'cpac-default-pipeline' + >>> diff['pipeline_setup']['pipeline_name'].left + 'cpac_fmriprep-options' + >>> diff.left['pipeline_setup']['pipeline_name'] + 'cpac_fmriprep-options' + >>> diff['pipeline_setup']['pipeline_name'].minuend + 'cpac_fmriprep-options' + >>> diff.minuend['pipeline_setup']['pipeline_name'] + 'cpac_fmriprep-options' + >>> diff['pipeline_setup']['pipeline_name'].right + 'cpac-default-pipeline' + >>> diff.right['pipeline_setup']['pipeline_name'] + 'cpac-default-pipeline' + >>> diff['pipeline_setup']['pipeline_name'].subtrahend + 'cpac-default-pipeline' + >>> diff.subtrahend['pipeline_setup']['pipeline_name'] + 'cpac-default-pipeline' + """ + return dct_diff(self.dict(), other.dict()) + + def _nonestr_to_None(self, d): + """Recursive method to type convert 'None' to None in nested config. + + Parameters + ---------- + d : any + config item to check + + Returns + ------- + d : any + same item, same type, but with 'none' strings converted to + Nonetypes + """ + if isinstance(d, str) and d.lower() == "none": + return None + if isinstance(d, list): + return [self._nonestr_to_None(i) for i in d] + if isinstance(d, set): + return {self._nonestr_to_None(i) for i in list(d)} + if isinstance(d, dict): + return {i: self._nonestr_to_None(d[i]) for i in d} + return d + + @staticmethod + def _check_keys(keys: Iterable) -> None: + """Check that keys are iterable and at least 1 key is provided.""" + if not keys: + if isinstance(keys, Iterable): + error = KeyError + msg = "No keys provided to `set_nested`." + else: + error = TypeError + msg = f"`set_nested` keys must be iterable, got {type(keys)}." + raise error(msg) + + def get_nested(self, _d: "NestedKeyMixin | _DICT", keys: Iterable) -> Any: + """Get a value from a Configuration dictionary given a nested key.""" + self._check_keys(keys) + if _d is None: + _d = {} + if isinstance(keys, str): + return _d[keys] + if isinstance(keys, (list, tuple)): + if len(keys) > 1: + return self.get_nested(_d[keys[0]], keys[1:]) + assert len(keys) == 1 + return _d[keys[0]] + return _d + + def keys(self) -> KeysView[Any]: + """Show toplevel keys of a nested dict.""" + return self.dict().keys() + + @overload + def set_nested( + self, d: "NestedKeyMixin", keys: Iterable, value: Any + ) -> "NestedKeyMixin": ... + @overload + def set_nested(self, d: _DICT, keys: Iterable, value: Any) -> _DICT: ... + def set_nested( + self, d: "NestedKeyMixin | _DICT", keys: Iterable, value: Any + ) -> "NestedKeyMixin | _DICT": + """Set a nested key in a nested dictionary.""" + self._check_keys(keys) + if isinstance(keys, str): + d[keys] = value + elif isinstance(keys, (list, tuple)): + if len(keys) > 1: + d[keys[0]] = self.set_nested(d[keys[0]], keys[1:], value) + else: + assert len(keys) == 1 + d[keys[0]] = value + return d + + def _check_if_switch(self, key: CONFIG_KEY_TYPE, error: bool = False) -> bool: + """Check if a given entity is a switch. + + Parameters + ---------- + key : str or list of str + key to check + + error : bool + raise a TypeError if not a switch + + Returns + ------- + bool + True if the given key is a switch, False otherwise + + Examples + -------- + >>> c = Configuration() + >>> c._check_if_switch('anatomical_preproc') + False + >>> c._check_if_switch(['anatomical_preproc']) + False + >>> c._check_if_switch(['anatomical_preproc', 'run']) + True + """ + _maybe_switch = self[key] + if isinstance(_maybe_switch, bool): + return True + if isinstance(_maybe_switch, list): + _answer = all(isinstance(_, bool) for _ in _maybe_switch) + if _answer: + return _answer + if error: + msg = f"`{key}` is not a switch in {self!s}." + raise TypeError(msg) + return False + + def _switch_bool(self, key: CONFIG_KEY_TYPE, value: bool, exclusive: bool) -> bool: + """Return True if the key is set to the given value or False otherwise. + + Parameters + ---------- + key : str or list of str + key to check + + value : bool + value to check for + + exclusive : bool + return False if forking (both True and False) + + Returns + ------- + bool + True if the given key is set to the given value or False + otherwise. If exclusive is True, return False if the key + is set to both True and False. + """ + if not (exclusive and self.switch_is_on_off(key)): + if isinstance(self[key], bool): + return self[key] is value + if isinstance(self[key], list): + return value in self[key] + return False + + def switch_is_off(self, key: CONFIG_KEY_TYPE, exclusive: bool = False) -> bool: + """Return True if the key is set to 'off' OR 'on' and 'off' or False otherwise. + + Used for tracking forking. + + Parameters + ---------- + key : str or list of str + key to check + + exclusive : bool, optional, default: False + return False if the key is set to 'on' and 'off' + + Returns + ------- + bool + True if key is set to 'off', False if not set to 'off'. + If exclusive is set to True, return False if the key is + set to 'on' and 'off'. + + Examples + -------- + >>> c = Configuration() + >>> c.switch_is_off(['nuisance_corrections', '2-nuisance_regression', + ... 'run']) + True + >>> c = Configuration({'nuisance_corrections': { + ... '2-nuisance_regression': {'run': [True, False]}}}) + >>> c.switch_is_off(['nuisance_corrections', '2-nuisance_regression', + ... 'run']) + True + >>> c.switch_is_off(['nuisance_corrections', '2-nuisance_regression', + ... 'run'], exclusive=True) + False + """ + self._check_if_switch(key, True) + return self._switch_bool(key, False, exclusive) + + def switch_is_on(self, key: CONFIG_KEY_TYPE, exclusive: bool = False) -> bool: + """Return True if the key is set to 'on' OR 'on' and 'off' or False otherwise. + + Used for tracking forking. + + Parameters + ---------- + key : str or list of str + key to check + + exclusive : bool, optional, default: False + return False if the key is set to 'on' and 'off' + + Returns + ------- + bool + True if key is set to 'on', False if not set to 'on'. + If exclusive is set to True, return False if the key is + set to 'on' and 'off'. + + Examples + -------- + >>> c = Configuration() + >>> c.switch_is_on(['nuisance_corrections', '2-nuisance_regression', + ... 'run']) + False + >>> c = Configuration({'nuisance_corrections': { + ... '2-nuisance_regression': {'run': [True, False]}}}) + >>> c.switch_is_on(['nuisance_corrections', '2-nuisance_regression', + ... 'run']) + True + >>> c.switch_is_on(['nuisance_corrections', '2-nuisance_regression', + ... 'run'], exclusive=True) + False + """ + self._check_if_switch(key, True) + return self._switch_bool(key, True, exclusive) + + def switch_is_on_off(self, key: CONFIG_KEY_TYPE) -> bool: + """Return True if the key is set to both 'on' and 'off' or False otherwise. + + Used for tracking forking. + + Parameters + ---------- + key : str or list of str + key to check + + Returns + ------- + bool + True if key is set to 'on' and 'off', False otherwise + + Examples + -------- + >>> c = Configuration() + >>> c.switch_is_on_off(['nuisance_corrections', + ... '2-nuisance_regression', 'run']) + False + >>> c = Configuration({'nuisance_corrections': { + ... '2-nuisance_regression': {'run': [True, False]}}}) + >>> c.switch_is_on_off(['nuisance_corrections', + ... '2-nuisance_regression', 'run']) + True + """ + self._check_if_switch(key, True) + if isinstance(self[key], list): + return True in self[key] and False in self[key] + return False + + def key_type_error(self, key): + """Raise a KeyError if an inappropriate type of key is attempted.""" + raise KeyError( + " ".join( + [ + "Configuration key must be a string, list, or tuple;", + type(key).__name__, + f"`{key!s}`", + "was given.", + ] + ) + ) + + +class Configuration(NestedKeyMixin): """ Class to set dictionary keys as map attributes. @@ -104,6 +443,25 @@ class Configuration: 'slack_420349_preconfig' """ + amplitude_low_frequency_fluctuation: dict + anatomical_preproc: dict + FROM: str + functional_preproc: dict + longitudinal_template_generation: dict + network_centrality: dict + nuisance_corrections: dict + pipeline_setup: dict + post_processing: dict + PyPEER: dict + regional_homogeneity: dict + registration_workflows: dict + seed_based_correlation_analysis: dict + segmentation: dict + skip_env_check: bool + surface_analysis: dict + timeseries_extraction: dict + voxel_mirrored_homotopic_connectivity: dict + def __init__( self, config_map: Optional[dict] = None, skip_env_check: bool = False ) -> None: @@ -159,6 +517,7 @@ def __init__( regressor["Name"] = nipype_friendly_name(regressor["Name"]) config_map = schema(config_map) + assert isinstance(config_map, dict) # remove 'skip env check' now that the config is validated if "skip env check" in config_map: @@ -171,88 +530,34 @@ def __init__( for key in config_map: # set attribute setattr(self, key, self.set_without_ENV(config_map[key])) - else: - # set FSLDIR to the environment $FSLDIR if the user sets it to - # 'FSLDIR' in the pipeline config file - _FSLDIR = config_map.get("FSLDIR") - if _FSLDIR and bool(re.match(r"^[\$\{]{0,2}?FSLDIR[\}]?$", _FSLDIR)): - config_map["FSLDIR"] = os.environ["FSLDIR"] - for key in config_map: - # set attribute - setattr(self, key, self.set_from_ENV(config_map[key])) - self._update_attr() - - # set working directory as an environment variable - os.environ["CPAC_WORKDIR"] = self["pipeline_setup", "working_directory", "path"] - - def __str__(self): - return f"C-PAC Configuration ('{self['pipeline_setup', 'pipeline_name']}')" - - def __repr__(self): - """Show Configuration as a dict when accessed directly.""" - return str(self.dict()) - - def __copy__(self): - newone = type(self)({}) - newone.__dict__.update(self.__dict__) - newone._update_attr() - return newone - - def __getitem__(self, key): - if isinstance(key, str): - return getattr(self, key) - if isinstance(key, (list, tuple)): - return self.get_nested(self, key) - self.key_type_error(key) - return None - - def __setitem__(self, key, value): - if isinstance(key, str): - setattr(self, key, value) - elif isinstance(key, (list, tuple)): - self.set_nested(self, key, value) - else: - self.key_type_error(key) - - def __sub__(self: "Configuration", other: "Configuration"): - """Return the set difference between two Configurations. - - Examples - -------- - >>> diff = (Preconfiguration('fmriprep-options') - ... - Preconfiguration('default')) - >>> diff['pipeline_setup']['pipeline_name'] - ('cpac_fmriprep-options', 'cpac-default-pipeline') - >>> diff['pipeline_setup']['pipeline_name'].s_value - 'cpac_fmriprep-options' - >>> diff['pipeline_setup']['pipeline_name'].t_value - 'cpac-default-pipeline' - >>> diff.s_value['pipeline_setup']['pipeline_name'] - 'cpac_fmriprep-options' - >>> diff.t_value['pipeline_setup']['pipeline_name'] - 'cpac-default-pipeline' - >>> diff['pipeline_setup']['pipeline_name'].left - 'cpac_fmriprep-options' - >>> diff.left['pipeline_setup']['pipeline_name'] - 'cpac_fmriprep-options' - >>> diff['pipeline_setup']['pipeline_name'].minuend - 'cpac_fmriprep-options' - >>> diff.minuend['pipeline_setup']['pipeline_name'] - 'cpac_fmriprep-options' - >>> diff['pipeline_setup']['pipeline_name'].right - 'cpac-default-pipeline' - >>> diff.right['pipeline_setup']['pipeline_name'] - 'cpac-default-pipeline' - >>> diff['pipeline_setup']['pipeline_name'].subtrahend - 'cpac-default-pipeline' - >>> diff.subtrahend['pipeline_setup']['pipeline_name'] - 'cpac-default-pipeline' - """ - return dct_diff(self.dict(), other.dict()) + else: + # set FSLDIR to the environment $FSLDIR if the user sets it to + # 'FSLDIR' in the pipeline config file + _FSLDIR = config_map.get("FSLDIR") + if _FSLDIR and bool(re.match(r"^[\$\{]{0,2}?FSLDIR[\}]?$", _FSLDIR)): + config_map["FSLDIR"] = os.environ["FSLDIR"] + for key in config_map: + # set attribute + setattr(self, key, self.set_from_ENV(config_map[key])) + self._update_attr() - def dict(self) -> dict[Any, Any]: - """Show contents of a C-PAC configuration as a dict.""" - return {k: v for k, v in self.__dict__.items() if not callable(v)} + # set working directory as an environment variable + os.environ["CPAC_WORKDIR"] = self["pipeline_setup", "working_directory", "path"] + + def __str__(self): + """Return string representation of a Configuration instance.""" + return f"C-PAC Configuration ('{self['pipeline_setup', 'pipeline_name']}')" + + def __repr__(self): + """Show Configuration as a dict when accessed directly.""" + return str(self.dict()) + + def __copy__(self): + """Copy a pipeline Configuration.""" + newone = type(self)({}) + newone.__dict__.update(self.__dict__) + newone._update_attr() + return newone def get(self, key: Any, default: Any = None, /) -> Any: """Provide convenience access from `Configuration` to :meth:`dict.get` . @@ -269,34 +574,6 @@ def get(self, key: Any, default: Any = None, /) -> Any: """ return self.dict().get(key, default) - def keys(self): - """Show toplevel keys of a C-PAC configuration dict.""" - return self.dict().keys() - - def _nonestr_to_None(self, d): - """Recursive method to type convert 'None' to None in nested config. - - Parameters - ---------- - d : any - config item to check - - Returns - ------- - d : any - same item, same type, but with 'none' strings converted to - Nonetypes - """ - if isinstance(d, str) and d.lower() == "none": - return None - if isinstance(d, list): - return [self._nonestr_to_None(i) for i in d] - if isinstance(d, set): - return {self._nonestr_to_None(i) for i in d} - if isinstance(d, dict): - return {i: self._nonestr_to_None(d[i]) for i in d} - return d - def set_from_ENV(self, conf): # pylint: disable=invalid-name """Replace strings like $VAR and ${VAR} with environment variable values. @@ -369,9 +646,11 @@ def set_without_ENV(self, conf): # pylint: disable=invalid-name return conf def sub_pattern(self, pattern, orig_key): + """Make a defined pattern substitution.""" return orig_key.replace(pattern, self[pattern[2:-1].split(".")]) def check_pattern(self, orig_key, tags=None): + """Make defined pattern substitutions.""" if tags is None: tags = [] if isinstance(orig_key, dict): @@ -425,222 +704,31 @@ def check_path(key): setattr(self, attr_key, new_key) def update(self, key, val=ConfigurationDictUpdateConflation()): + """Update a C-PAC pipeline Configuration.""" if isinstance(key, dict): raise ConfigurationDictUpdateConflation if isinstance(val, Exception): raise val setattr(self, key, val) - def get_nested(self, _d, keys): - if _d is None: - _d = {} - if isinstance(keys, str): - return _d[keys] - if isinstance(keys, (list, tuple)): - if len(keys) > 1: - return self.get_nested(_d[keys[0]], keys[1:]) - return _d[keys[0]] - return _d - - def set_nested(self, d, keys, value): # pylint: disable=invalid-name - if isinstance(keys, str): - d[keys] = value - elif isinstance(keys, (list, tuple)): - if len(keys) > 1: - d[keys[0]] = self.set_nested(d[keys[0]], keys[1:], value) - else: - d[keys[0]] = value - return d - - def _check_if_switch(self, key: CONFIG_KEY_TYPE, error: bool = False) -> bool: - """Check if a given entity is a switch. - - Parameters - ---------- - key : str or list of str - key to check - - error : bool - raise a TypeError if not a switch - - Returns - ------- - bool - True if the given key is a switch, False otherwise - - Examples - -------- - >>> c = Configuration() - >>> c._check_if_switch('anatomical_preproc') - False - >>> c._check_if_switch(['anatomical_preproc']) - False - >>> c._check_if_switch(['anatomical_preproc', 'run']) - True - """ - _maybe_switch = self[key] - if isinstance(_maybe_switch, bool): - return True - if isinstance(_maybe_switch, list): - _answer = all(isinstance(_, bool) for _ in _maybe_switch) - if _answer: - return _answer - if error: - msg = f"`{key}` is not a switch in {self!s}." - raise TypeError(msg) - return False - - def _switch_bool(self, key: CONFIG_KEY_TYPE, value: bool, exclusive: bool) -> bool: - """Return True if the key is set to the given value or False otherwise. - - Parameters - ---------- - key : str or list of str - key to check - - value : bool - value to check for - - exclusive : bool - return False if forking (both True and False) - - Returns - ------- - bool - True if the given key is set to the given value or False - otherwise. If exclusive is True, return False if the key - is set to both True and False. - """ - if not (exclusive and self.switch_is_on_off(key)): - if isinstance(self[key], bool): - return self[key] is value - if isinstance(self[key], list): - return value in self[key] - return False - - def switch_is_off(self, key: CONFIG_KEY_TYPE, exclusive: bool = False) -> bool: - """Return True if the key is set to 'off' OR 'on' and 'off' or False otherwise. - - Used for tracking forking. - - Parameters - ---------- - key : str or list of str - key to check - - exclusive : bool, optional, default: False - return False if the key is set to 'on' and 'off' - - Returns - ------- - bool - True if key is set to 'off', False if not set to 'off'. - If exclusive is set to True, return False if the key is - set to 'on' and 'off'. - - Examples - -------- - >>> c = Configuration() - >>> c.switch_is_off(['nuisance_corrections', '2-nuisance_regression', - ... 'run']) - True - >>> c = Configuration({'nuisance_corrections': { - ... '2-nuisance_regression': {'run': [True, False]}}}) - >>> c.switch_is_off(['nuisance_corrections', '2-nuisance_regression', - ... 'run']) - True - >>> c.switch_is_off(['nuisance_corrections', '2-nuisance_regression', - ... 'run'], exclusive=True) - False - """ - self._check_if_switch(key, True) - return self._switch_bool(key, False, exclusive) - - def switch_is_on(self, key: CONFIG_KEY_TYPE, exclusive: bool = False) -> bool: - """Return True if the key is set to 'on' OR 'on' and 'off' or False otherwise. - - Used for tracking forking. - - Parameters - ---------- - key : str or list of str - key to check - - exclusive : bool, optional, default: False - return False if the key is set to 'on' and 'off' - - Returns - ------- - bool - True if key is set to 'on', False if not set to 'on'. - If exclusive is set to True, return False if the key is - set to 'on' and 'off'. - - Examples - -------- - >>> c = Configuration() - >>> c.switch_is_on(['nuisance_corrections', '2-nuisance_regression', - ... 'run']) - False - >>> c = Configuration({'nuisance_corrections': { - ... '2-nuisance_regression': {'run': [True, False]}}}) - >>> c.switch_is_on(['nuisance_corrections', '2-nuisance_regression', - ... 'run']) - True - >>> c.switch_is_on(['nuisance_corrections', '2-nuisance_regression', - ... 'run'], exclusive=True) - False - """ - self._check_if_switch(key, True) - return self._switch_bool(key, True, exclusive) - - def switch_is_on_off(self, key: CONFIG_KEY_TYPE) -> bool: - """Return True if the key is set to both 'on' and 'off' or False otherwise. - - Used for tracking forking. - - Parameters - ---------- - key : str or list of str - key to check - - Returns - ------- - bool - True if key is set to 'on' and 'off', False otherwise - - Examples - -------- - >>> c = Configuration() - >>> c.switch_is_on_off(['nuisance_corrections', - ... '2-nuisance_regression', 'run']) - False - >>> c = Configuration({'nuisance_corrections': { - ... '2-nuisance_regression': {'run': [True, False]}}}) - >>> c.switch_is_on_off(['nuisance_corrections', - ... '2-nuisance_regression', 'run']) - True - """ - self._check_if_switch(key, True) - if isinstance(self[key], list): - return True in self[key] and False in self[key] - return False - - def key_type_error(self, key): - """Raise a KeyError if an inappropriate type of key is attempted.""" - raise KeyError( - " ".join( - [ - "Configuration key must be a string, list, or tuple;", - type(key).__name__, - f"`{key!s}`", - "was given.", - ] - ) + @overload + def orientation_node(self, name: str, node_type: type[MapNode]) -> MapNode: ... + @overload + def orientation_node(self, name: str, node_type: type[Node]) -> Node: ... + def orientation_node( + self, name: str, node_type: type[Node | MapNode] = Node + ) -> Node | MapNode: + """Return a node configured to resample an input with AFNI 3dresample.""" + from CPAC.utils.nifti_utils import orientation_node + + orientation = cast( + Literal["RPI", "LPI", "RAI", "LAI", "RAS", "LAS", "RPS", "LPS"], + self["pipeline_setup", "desired_orientation"], ) + return orientation_node(name=name, orientation=orientation, node_type=node_type) -def check_pname(p_name: str, pipe_config: Configuration) -> str: +def check_pname(p_name: Optional[str], pipe_config: Configuration) -> str: """Check / set `p_name`, the str representation of a pipeline for use in filetrees. Parameters @@ -752,9 +840,8 @@ def preconfig_yaml(preconfig_name="default", load=False): if load: with open(preconfig_yaml(preconfig_name), "r", encoding="utf-8") as _f: return yaml.safe_load(_f) - return p.resource_filename( - "CPAC", - os.path.join("resources", "configs", f"pipeline_config_{preconfig_name}.yml"), + return files("CPAC").joinpath( + f"resources/configs/pipeline_config_{preconfig_name}.yml" ) diff --git a/CPAC/utils/configuration/yaml_template.py b/CPAC/utils/configuration/yaml_template.py index f460f90e5d..a62efc70ff 100755 --- a/CPAC/utils/configuration/yaml_template.py +++ b/CPAC/utils/configuration/yaml_template.py @@ -27,14 +27,19 @@ from click import BadParameter import yaml -from CPAC.utils.configuration import Configuration, preconfig_yaml, Preconfiguration +from CPAC.utils.configuration import ( + Configuration, + NestedKeyMixin, + preconfig_yaml, + Preconfiguration, +) from CPAC.utils.monitoring import UTLOGGER from CPAC.utils.utils import update_config_dict, update_pipeline_values_1_8, YAML_BOOLS YAML_LOOKUP = {yaml_str: key for key, value in YAML_BOOLS.items() for yaml_str in value} -class YamlTemplate: # pylint: disable=too-few-public-methods +class YamlTemplate(NestedKeyMixin): """A class to link YAML comments to the contents of a YAML file. Attributes @@ -81,8 +86,6 @@ def __init__(self, original_yaml, base_config=None): self._dict = base_config.dict() self._parse_comments() - get_nested = Configuration.get_nested - def dump(self, new_dict, parents=None): """Dump YAML from a new dictionary with comments from template dictionary. diff --git a/CPAC/utils/interfaces/function/function.py b/CPAC/utils/interfaces/function/function.py index 34d01373d5..db09d30839 100644 --- a/CPAC/utils/interfaces/function/function.py +++ b/CPAC/utils/interfaces/function/function.py @@ -110,7 +110,9 @@ def get_function_name_from_source(function_source: str) -> str: def create_function_from_source( - function_source: str, imports: Optional[list[str]] = None, ns: Optional[dict] = None + function_source: str | bytes, + imports: Optional[list[str]] = None, + ns: Optional[dict] = None, ): """Return a function object from a function source. @@ -156,28 +158,28 @@ class Function(NipypeFunction): def __init__( self, - input_names=None, - output_names="out", - function=None, - imports=None, - as_module=False, + input_names: Optional[list[str] | str] = None, + output_names: list[str] | str = "out", + function: Optional[Callable] = None, + imports: Optional[list[str]] = None, + as_module: bool = False, **inputs, ): """Initialize a :py:func`~CPAC.utils.interfaces.function.Function` interface. Parameters ---------- - input_names : single str or list or None + input_names names corresponding to function inputs if ``None``, derive input names from function argument names - output_names : single str or list + output_names names corresponding to function outputs (default: 'out'). if list of length > 1, has to match the number of outputs - function : callable + function callable python object. must be able to execute in an isolated namespace (possibly in concert with the ``imports`` parameter) - imports : list of strings + imports list of import statements that allow the function to execute in an otherwise empty namespace. If these collide with imports defined via the :py:meth:`Function.sig_imports` diff --git a/CPAC/utils/monitoring/__init__.py b/CPAC/utils/monitoring/__init__.py index 552d5fa488..f681bf7ded 100644 --- a/CPAC/utils/monitoring/__init__.py +++ b/CPAC/utils/monitoring/__init__.py @@ -26,6 +26,7 @@ FMLOGGER, getLogger, IFLOGGER, + init_loggers, set_up_logger, UTLOGGER, WFLOGGER, @@ -44,6 +45,7 @@ "FMLOGGER", "getLogger", "IFLOGGER", + "init_loggers", "LoggingHTTPServer", "LoggingRequestHandler", "log_nodes_cb", diff --git a/CPAC/utils/monitoring/custom_logging.py b/CPAC/utils/monitoring/custom_logging.py index 3d8d1b842a..0b6a6f5ea9 100644 --- a/CPAC/utils/monitoring/custom_logging.py +++ b/CPAC/utils/monitoring/custom_logging.py @@ -1,4 +1,4 @@ -# Copyright (C) 2022-2024 C-PAC Developers +# Copyright (C) 2022-2025 C-PAC Developers # This file is part of C-PAC. @@ -18,16 +18,46 @@ import logging import os +from pathlib import Path import subprocess from sys import exc_info as sys_exc_info from traceback import print_exception -from typing import Optional, Sequence +from typing import Literal, Optional, Sequence, TYPE_CHECKING, TypeAlias -from nipype import logging as nipype_logging +import yaml +from nipype import config as nipype_config, logging as nipype_logging from CPAC.utils.docs import docstring_parameter from CPAC.utils.monitoring.config import MOCK_LOGGERS +if TYPE_CHECKING: + from CPAC.utils.configuration import Configuration +LogLevel: TypeAlias = ( + Literal[ + "CRITICAL", + "critical", + "Critical", + "DEBUG", + "debug", + "Debug", + "ERROR", + "error", + "Error", + "INFO", + "info", + "Info", + "NOTSET", + "notset", + "Notset", + "NotSet", + "notSet", + "WARNING", + "warning", + "Warning", + ] + | int +) + def failed_to_start(log_dir, exception): """Launch a failed-to-start logger for a run that failed to start. @@ -46,17 +76,8 @@ def failed_to_start(log_dir, exception): logger.exception(exception) -def getLogger(name): # pylint: disable=invalid-name - """Get a mock logger if one exists, falling back on real loggers. - - Parameters - ---------- - name : str - - Returns - ------- - logger : CPAC.utils.monitoring.custom_logging.MockLogger or logging.Logger - """ +def getLogger(name: str) -> "logging.Logger | MockLogger": # pylint: disable=invalid-name + """Get a mock logger if one exists, falling back on real loggers.""" if name in MOCK_LOGGERS: return MOCK_LOGGERS[name] logger = nipype_logging.getLogger(name) @@ -143,6 +164,22 @@ def log_subprocess(cmd, *args, raise_error=True, **kwargs): return output, 0 +class ListToSetYamlLoader(yaml.Loader): + """Custom YAML loader to convert lists to sets.""" + + def construct_sequence( # pyright: ignore[reportIncompatibleMethodOverride] + self, node, deep=False + ) -> set[str]: + """Convert YAML sequence to a set.""" + return set(super().construct_sequence(node, deep)) + + +ListToSetYamlLoader.add_constructor( + yaml.resolver.BaseResolver.DEFAULT_SEQUENCE_TAG, + ListToSetYamlLoader.construct_sequence, +) + + # pylint: disable=too-few-public-methods class MockHandler: """Handler for MockLogger.""" @@ -155,7 +192,10 @@ def __init__(self, filename): class MockLogger: """Mock logging.Logger to provide API without keeping the logger in memory.""" - def __init__(self, name, filename, level, log_dir): + def __init__( + self, name: str, filename: str, level: LogLevel, log_dir: Path | str + ) -> None: + """Initialize a mock logger.""" self.name = name self.level = level self.handlers = [MockHandler(os.path.join(log_dir, filename))] @@ -210,6 +250,24 @@ def _get_first_file_handler( return handler return None + def yaml_contents(self) -> dict: + """If the logger's first handler is a YAML file, return the contents and delete them from the logger.""" + file = self._get_first_file_handler(self.handlers) + if hasattr(file, "baseFilename"): + file = Path(getattr(file, "baseFilename")) + if file.suffix == ".yml": + with file.open("r", encoding="utf-8") as f: + contents = yaml.load(f.read(), Loader=ListToSetYamlLoader) + with file.open("w", encoding="utf-8") as f: + f.write("") + return contents + error = TypeError + msg = f"Could not load YAML contents from {file}" + else: + error = FileNotFoundError + msg = f"Could not find file handler for {self.name}" + raise error(msg) + def _lazy_sub(message, *items): """Given lazy-logging syntax, return string with substitutions. @@ -240,34 +298,37 @@ def _lazy_sub(message, *items): def set_up_logger( - name, filename=None, level=None, log_dir=None, mock=False, overwrite_existing=False -): + name: str, + filename: Optional[str] = None, + level: Optional[LogLevel] = None, + log_dir: Optional[Path | str] = None, + mock: bool = False, +) -> logging.Logger | MockLogger: r"""Initialize a logger. Parameters ---------- - name : str + name logger name (for subsequent calls to ``logging.getLogger``) to write to the same log file) - filename : str, optional + filename filename to write log to. If not specified, filename will be the same as ``name`` with the extension ``log`` - level : str, optional - one of ``{critical, error, warning, info, debug, notset}``, - case-insensitive + level + https://docs.python.org/3/library/logging.html#levels - log_dir : str, optional + log_dir - mock : bool, optional + mock if ``True``, return a ``CPAC.utils.monitoring.MockLogger`` instead of a ``logging.Logger`` Returns ------- - logger : logging.Handler - initialized logging Handler + logger + initialized logger Examples -------- @@ -292,20 +353,73 @@ def set_up_logger( """ if filename is None: filename = f"{name}.log" - try: - level = getattr(logging, level.upper()) - except AttributeError: + if isinstance(level, str): + try: + level = getattr(logging, level.upper()) + except AttributeError: + pass + if not level: level = logging.NOTSET - if log_dir is None: - log_dir = os.getcwd() - filepath = os.path.join(log_dir, filename) - if overwrite_existing and os.path.exists(filepath): - with open(filepath, "w") as log_file: - log_file.write("") + log_dir = Path(log_dir) if log_dir else Path.cwd() + filepath = log_dir / filename + if not filepath.exists(): + filepath.parent.mkdir(parents=True, exist_ok=True) if mock: return MockLogger(name, filename, level, log_dir) logger = getLogger(name) + if isinstance(logger, MockLogger): + return logger logger.setLevel(level) handler = logging.FileHandler(filepath) logger.addHandler(handler) return logger + + +def init_loggers( + subject_id: str, + cpac_config: "Configuration", + log_dir: str, + mock: bool = True, + longitudinal: bool = False, +) -> None: + """Set up and configure loggers.""" + from CPAC.utils.datasource import bidsier_prefix + + if "subject_id" not in cpac_config: + cpac_config["subject_id"] = subject_id + set_up_logger( + f"{cpac_config['subject_id']}_expectedOutputs", + filename=f"{bidsier_prefix(cpac_config['subject_id'])}_expectedOutputs.yml", + level="info", + log_dir=log_dir, + mock=mock, + ) + + if cpac_config["pipeline_setup", "Debugging", "verbose"]: + set_up_logger("CPAC.engine", level="debug", log_dir=log_dir, mock=True) + + nipype_config.update_config( + { + "logging": { + "log_directory": log_dir, + "log_to_file": bool( + getattr( + cpac_config["pipeline_setup", "log_directory"], + "run_logging", + True, + ) + ), + }, + "execution": { + "crashfile_format": "txt", + "resource_monitor_frequency": 0.2, + "stop_on_first_crash": cpac_config[ + "pipeline_setup", "system_config", "fail_fast" + ], + }, + } + ) + + nipype_config.enable_resource_monitor() + + nipype_logging.update_logging(nipype_config) diff --git a/CPAC/utils/nifti_utils.py b/CPAC/utils/nifti_utils.py index 04db0d25ac..dd06fa9cd0 100644 --- a/CPAC/utils/nifti_utils.py +++ b/CPAC/utils/nifti_utils.py @@ -1,4 +1,4 @@ -# Copyright (C) 2019-2024 C-PAC Developers +# Copyright (C) 2019-2025 C-PAC Developers # This file is part of C-PAC. @@ -17,9 +17,13 @@ """Utlities for NIfTI images.""" import os +from typing import Literal, overload import numpy as np import nibabel as nib +from nipype.interfaces.afni import utils as afni_utils + +from CPAC.pipeline import nipype_pipeline_engine as pe def nifti_image_input( @@ -93,3 +97,35 @@ def inverse_nifti_values(image): out_data[zeros] = 0 return nib.nifti1.Nifti1Image(out_data, img.affine) + + +@overload +def orientation_node( + name: str, + orientation: Literal["RPI", "LPI", "RAI", "LAI", "RAS", "LAS", "RPS", "LPS"], + node_type: type[pe.MapNode], +) -> pe.MapNode: ... +@overload +def orientation_node( + name: str, + orientation: Literal["RPI", "LPI", "RAI", "LAI", "RAS", "LAS", "RPS", "LPS"], + node_type: type[pe.Node], +) -> pe.Node: ... +def orientation_node( + name: str, + orientation: Literal["RPI", "LPI", "RAI", "LAI", "RAS", "LAS", "RPS", "LPS"], + node_type: type[pe.Node | pe.MapNode] = pe.Node, +) -> pe.Node | pe.MapNode: + """Return a node configured to resample an input with AFNI 3dresample.""" + kwargs = { + "interface": afni_utils.Resample( + orientation=orientation, + outputtype="NIFTI_GZ", + ), + "name": name, + "mem_gb": 0, + "mem_x": (0.0115, "in_file", "t"), + } + if node_type == pe.MapNode: + kwargs["iterfield"] = ["in_file", "out_file"] + return node_type(**kwargs) diff --git a/CPAC/utils/utils.py b/CPAC/utils/utils.py index b459262993..0505c28c9c 100644 --- a/CPAC/utils/utils.py +++ b/CPAC/utils/utils.py @@ -1,4 +1,4 @@ -# Copyright (C) 2012-2024 C-PAC Developers +# Copyright (C) 2012-2025 C-PAC Developers # This file is part of C-PAC. @@ -73,7 +73,7 @@ def get_last_prov_entry(prov): return prov[-1] -def check_prov_for_regtool(prov): +def check_prov_for_regtool(prov) -> Optional[Literal["ants", "fsl"]]: """Check provenance for registration tool.""" last_entry = get_last_prov_entry(prov) last_node = last_entry.split(":")[1] @@ -101,22 +101,6 @@ def check_prov_for_regtool(prov): return None -def check_prov_for_motion_tool(prov): - """Check provenance for motion correction tool.""" - last_entry = get_last_prov_entry(prov) - last_node = last_entry.split(":")[1] - if "3dvolreg" in last_node.lower(): - return "3dvolreg" - if "mcflirt" in last_node.lower(): - return "mcflirt" - # check entire prov - if "3dvolreg" in str(prov): - return "3dvolreg" - if "mcflirt" in str(prov): - return "mcflirt" - return None - - def _get_flag(in_flag): return in_flag @@ -159,6 +143,7 @@ def create_id_string( fwhm=None, subdir=None, extension=None, + subject_level: bool = False, ): """Create the unique key-value identifier string for BIDS-Derivatives file names. @@ -177,6 +162,9 @@ def create_id_string( from CPAC.utils.bids_utils import combine_multiple_entity_instances, res_in_filename + if "longitudinal-template" in resource: + resource = resource.replace("longitudinal-template", "").replace("__", "") + if atlas_id: if "_desc-" in atlas_id: atlas, desc = atlas_id.split("_desc-") @@ -186,16 +174,23 @@ def create_id_string( atlas_id = atlas_id.replace("_desc-", "") resource = f"atlas-{atlas_id}_{resource}" - part_id = unique_id.split("_")[0] - ses_id = unique_id.split("_")[1] + id_parts = [] + if "_" in unique_id: + part_id, ses_id = unique_id.split("_", 1) + if "ses-" not in ses_id: + ses_id = f"ses-{ses_id}" + id_parts.append(ses_id) + else: + part_id = unique_id if "sub-" not in part_id: part_id = f"sub-{part_id}" - if "ses-" not in ses_id: - ses_id = f"ses-{ses_id}" + id_parts.insert(0, part_id) if scan_id: - out_filename = f"{part_id}_{ses_id}_task-{scan_id}_{resource}" - else: - out_filename = f"{part_id}_{ses_id}_{resource}" + if "task-" not in scan_id: + scan_id = f"task-{scan_id}" + id_parts.append(scan_id) + id_parts.append(resource) + out_filename = "_".join(id_parts) template_tag = template_desc.split(" -")[0] if template_desc else "*" for prefix in ["space-", "from-", "to-"]: diff --git a/CPAC/vmhc/vmhc.py b/CPAC/vmhc/vmhc.py index 3c547a8e2f..38499f511c 100644 --- a/CPAC/vmhc/vmhc.py +++ b/CPAC/vmhc/vmhc.py @@ -1,11 +1,28 @@ +# Copyright (C) 2012-2025 C-PAC Developers + +# This file is part of C-PAC. + +# C-PAC is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or (at your +# option) any later version. + +# C-PAC is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +# License for more details. + +# You should have received a copy of the GNU Lesser General Public +# License along with C-PAC. If not, see . +"""Voxel-Mirrored Homotopic Connectivity.""" + from nipype.interfaces import fsl from nipype.interfaces.afni import preprocess from CPAC.image_utils import spatial_smoothing from CPAC.pipeline import nipype_pipeline_engine as pe from CPAC.pipeline.nodeblock import nodeblock -from CPAC.registration.registration import apply_transform -from CPAC.utils.utils import check_prov_for_regtool +from CPAC.registration.utils import apply_transform from CPAC.vmhc import * from .utils import * @@ -60,8 +77,7 @@ def smooth_func_vmhc(wf, cfg, strat_pool, pipe_num, opt=None): outputs=["space-symtemplate_desc-sm_bold"], ) def warp_timeseries_to_sym_template(wf, cfg, strat_pool, pipe_num, opt=None): - xfm_prov = strat_pool.get_cpac_provenance("from-bold_to-symtemplate_mode-image_xfm") - reg_tool = check_prov_for_regtool(xfm_prov) + reg_tool = strat_pool.reg_tool("from-bold_to-symtemplate_mode-image_xfm") num_cpus = cfg.pipeline_setup["system_config"]["max_cores_per_participant"]