From cb4ce9eb6dc54b8493d33cce807e4f395198b741 Mon Sep 17 00:00:00 2001 From: NicolasGensollen Date: Fri, 31 Mar 2023 14:42:31 +0200 Subject: [PATCH] Split EPI pipeline in two, improve names, add docs.... --- .../dwi_preprocessing_using_t1_pipeline.py | 6 +- .../dwi_preprocessing_using_t1_workflows.py | 349 +++++++++++++----- 2 files changed, 259 insertions(+), 96 deletions(-) diff --git a/clinica/pipelines/dwi_preprocessing_using_t1/dwi_preprocessing_using_t1_pipeline.py b/clinica/pipelines/dwi_preprocessing_using_t1/dwi_preprocessing_using_t1_pipeline.py index bbc47f8151..91a4982511 100644 --- a/clinica/pipelines/dwi_preprocessing_using_t1/dwi_preprocessing_using_t1_pipeline.py +++ b/clinica/pipelines/dwi_preprocessing_using_t1/dwi_preprocessing_using_t1_pipeline.py @@ -366,9 +366,9 @@ def build_core_nodes(self): ("out_reference_b0", "inputnode.ref_b0")]), # TODO: check if really needed... (mask_b0_pre, eddy_fsl, [("mask_file", "inputnode.in_mask")]), # Magnetic susceptibility correction - (init_node, sdc, [("t1w", "inputnode.T1")]), - (eddy_fsl, sdc, [("outputnode.out_corrected", "inputnode.DWI")]), - (eddy_fsl, sdc, [("outputnode.out_rotated_bvecs", "inputnode.bvec")]), + (init_node, sdc, [("t1w", "inputnode.t1_filename")]), + (eddy_fsl, sdc, [("outputnode.out_corrected", "inputnode.dwi_filename")]), + (eddy_fsl, sdc, [("outputnode.out_rotated_bvecs", "inputnode.b_vectors_filename")]), # Bias correction (prepare_b0, bias, [("out_updated_bval", "in_bval")]), (sdc, bias, [("outputnode.DWIs_epicorrected", "in_file"), diff --git a/clinica/pipelines/dwi_preprocessing_using_t1/dwi_preprocessing_using_t1_workflows.py b/clinica/pipelines/dwi_preprocessing_using_t1/dwi_preprocessing_using_t1_workflows.py index eec0801ec3..da918a2ddf 100644 --- a/clinica/pipelines/dwi_preprocessing_using_t1/dwi_preprocessing_using_t1_workflows.py +++ b/clinica/pipelines/dwi_preprocessing_using_t1/dwi_preprocessing_using_t1_workflows.py @@ -153,7 +153,7 @@ def epi_pipeline( base_dir: str, delete_cache: bool = False, name="susceptibility_distortion_correction_using_t1", -): +) -> Workflow: """Perform EPI correction. This workflow allows to correct for echo-planar induced susceptibility artifacts without fieldmap @@ -162,6 +162,11 @@ def epi_pipeline( function (SyN algorithm). This workflow allows also a coregistration of DWIs with their respective baseline T1-weighted structural scans in order to latter combine tracks and cortex parcellation. + This pipeline needs: + - FSL + - ANTS + - c3d + Parameters ---------- base_dir: str @@ -173,16 +178,108 @@ def epi_pipeline( name: str, optional Name of the pipeline. + Returns + ------- + Workflow : + The Nipype workflow. + This workflow has the following inputs: + - "t1_filename": The path to the T1w image + - "dwi_filename": The path to the corrected DWI image + This corresponds to the 'out_corrected' output from the eddy_fsl pipeline. + - "b_vectors_filename": The path to the file holding the rotated B-vectors + This corresponds to the 'out_rotated_bvecs' output from the eddy_fsl pipeline. + + And the following outputs: + - "DWI_2_T1_Coregistration_matrix": ??? + - "epi_correction_deformation_field": ??? + - "epi_correction_affine_transform": ??? + - "epi_correction_image_warped": ??? + - "DWIs_epicorrected": ??? + - "warp_epi": ??? + - "out_bvec": ??? + Warnings -------- - This workflow rotates the b-vectors. + This workflow rotates the b-vectors. Notes ----- - Nir et al. (2015): Connectivity network measures predict volumetric atrophy in mild cognitive impairment - Leow et al. (2007): Statistical Properties of Jacobian Maps and the Realization of - Unbiased Large Deformation Nonlinear Image Registration + Nir et al. (2015): Connectivity network measures predict volumetric atrophy in mild cognitive impairment + Leow et al. (2007): Statistical Properties of Jacobian Maps and the Realization of + Unbiased Large Deformation Nonlinear Image Registration """ + import nipype.interfaces.utility as niu + import nipype.pipeline.engine as pe + + workflow_inputs = ["t1_filename", "dwi_filename", "b_vectors_filename"] + + inputnode = pe.Node( + niu.IdentityInterface(fields=workflow_inputs), + name="inputnode", + ) + ants_registration = perform_ants_registration() + ants_registration_outputs = [ + "merged_transforms", + "dwi_to_t1_co_registration_matrix", + "epi_correction_deformation_field", + "epi_correction_affine_transform", + "epi_correction_image_warped", + "rotated_b_vectors", + ] + + epi_correction = perform_dwi_epi_correction( + base_dir=base_dir, + delete_cache=delete_cache, + ) + epi_correction_outputs = ["epi_corrected_dwi_image"] + + outputnode = pe.Node( + niu.IdentityInterface( + fields=ants_registration_outputs + epi_correction_outputs, + ), + name="outputnode", + ) + + wf = pe.Workflow(name=name) + + connections = [ + ( + inputnode, + ants_registration, + [(inpt, f"inputnode.{inpt}") for inpt in workflow_inputs], + ), + ( + inputnode, + epi_correction, + [ + ("dwi_filename", "inputnode.dwi_filename"), + ("t1_filename", "inputnode.t1_filename"), + ], + ), + ( + ants_registration, + epi_correction, + [("outputnode.merged_transforms", "inputnode.merged_transforms")], + ), + ( + ants_registration, + outputnode, + [(f"outputnode.{outpt}", outpt) for outpt in ants_registration_outputs], + ), + ( + epi_correction, + outputnode, + [(f"outputnode.{outpt}", outpt) for outpt in epi_correction_outputs], + ), + ] + + wf.connect(connections) + + return wf + + +def perform_ants_registration(name="perform_ants_registration") -> Workflow: + """Step 1 of EPI pipeline.""" import nipype.interfaces.ants as ants import nipype.interfaces.c3 as c3 import nipype.interfaces.fsl as fsl @@ -190,47 +287,51 @@ def epi_pipeline( import nipype.pipeline.engine as pe from .dwi_preprocessing_using_t1_utils import ( - ants_apply_transforms, + broadcast_matrix_filename_to_match_b_vector_length, change_itk_transform_type, - delete_temp_dirs, - expend_matrix_list, - rotate_bvecs, + rotate_b_vectors, ) inputnode = pe.Node( - niu.IdentityInterface(fields=["T1", "DWI", "bvec"]), name="inputnode" + niu.IdentityInterface( + fields=["t1_filename", "dwi_filename", "b_vectors_filename"] + ), + name="inputnode", ) + split_dwi_volumes = pe.Node(fsl.Split(dimension="t"), name="split_dwi_volumes") + pick_reference_b0 = pe.Node(niu.Select(), name="pick_reference_b0") + pick_reference_b0.inputs.index = [0] - split = pe.Node(fsl.Split(dimension="t"), name="SplitDWIs") - pick_ref = pe.Node(niu.Select(), name="Pick_b0") - pick_ref.inputs.index = [0] - - flirt_b0_2_t1 = pe.Node(interface=fsl.FLIRT(dof=6), name="flirt_B0_2_T1") - flirt_b0_2_t1.inputs.interp = "spline" - flirt_b0_2_t1.inputs.cost = "normmi" - flirt_b0_2_t1.inputs.cost_func = "normmi" + co_register_reference_b0_to_t1 = pe.Node( + interface=fsl.FLIRT(dof=6), + name="co_register_reference_b0_to_t1", + ) + co_register_reference_b0_to_t1.inputs.interp = "spline" + co_register_reference_b0_to_t1.inputs.cost = "normmi" + co_register_reference_b0_to_t1.inputs.cost_func = "normmi" expend_matrix = pe.Node( interface=niu.Function( - input_names=["in_matrix", "in_bvec"], + input_names=["matrix_filename", "b_vectors_filename"], output_names=["out_matrix_list"], - function=expend_matrix_list, + function=broadcast_matrix_filename_to_match_b_vector_length, ), name="expend_matrix", ) - rot_bvec = pe.Node( + b_vectors_rotation = pe.Node( niu.Function( - input_names=["in_matrix", "in_bvec"], - output_names=["out_file"], - function=rotate_bvecs, + input_names=["matrix_filenames", "b_vectors_filename"], + output_names=["rotated_b_vectors_filename"], + function=rotate_b_vectors, ), - name="Rotate_Bvec", + name="b_vectors_rotation", ) ants_registration = pe.Node( interface=ants.registration.RegistrationSynQuick( - transform_type="br", dimension=3 + transform_type="br", + dimension=3, ), name="antsRegistrationSyNQuick", ) @@ -239,7 +340,7 @@ def epi_pipeline( c3d_flirt2ants.inputs.itk_transform = True c3d_flirt2ants.inputs.fsl2ras = True - change_transform = pe.Node( + change_transform_type = pe.Node( niu.Function( input_names=["input_affine_file"], output_names=["updated_affine_file"], @@ -248,8 +349,108 @@ def epi_pipeline( name="change_transform_type", ) - merge_transform = pe.Node(niu.Merge(3), name="MergeTransforms") + merge_transforms = pe.Node(niu.Merge(3), name="merge_transforms") + + outputnode = pe.Node( + niu.IdentityInterface( + fields=[ + "merged_transforms", + "dwi_to_t1_co_registration_matrix", + "epi_correction_deformation_field", + "epi_correction_affine_transform", + "epi_correction_image_warped", + "rotated_b_vectors", + ] + ), + name="outputnode", + ) + + wf = pe.Workflow(name=name) + + connections = [ + (inputnode, split_dwi_volumes, [("dwi_filename", "in_file")]), + (split_dwi_volumes, pick_reference_b0, [("out_files", "inlist")]), + (pick_reference_b0, co_register_reference_b0_to_t1, [("out", "in_file")]), + (inputnode, co_register_reference_b0_to_t1, [("t1_filename", "reference")]), + (inputnode, b_vectors_rotation, [("b_vectors_filename", "b_vectors_filename")]), + ( + co_register_reference_b0_to_t1, + expend_matrix, + [("out_matrix_file", "matrix_filename")], + ), + (inputnode, expend_matrix, [("b_vectors_filename", "b_vectors_filename")]), + (expend_matrix, b_vectors_rotation, [("out_matrix_list", "matrix_filenames")]), + (inputnode, ants_registration, [("t1_filename", "fixed_image")]), + ( + co_register_reference_b0_to_t1, + ants_registration, + [("out_file", "moving_image")], + ), + (inputnode, c3d_flirt2ants, [("t1_filename", "reference_file")]), + (pick_reference_b0, c3d_flirt2ants, [("out", "source_file")]), + ( + co_register_reference_b0_to_t1, + c3d_flirt2ants, + [("out_matrix_file", "transform_file")], + ), + ( + c3d_flirt2ants, + change_transform_type, + [("itk_transform", "input_affine_file")], + ), + (change_transform_type, merge_transforms, [("updated_affine_file", "in1")]), + (ants_registration, merge_transforms, [("out_matrix", "in2")]), + (ants_registration, merge_transforms, [("forward_warp_field", "in3")]), + (merge_transforms, outputnode, [("out", "merged_transforms")]), + ( + co_register_reference_b0_to_t1, + outputnode, + [("out_matrix_file", "dwi_to_t1_co_registration_matrix")], + ), + ( + ants_registration, + outputnode, + [ + ("forward_warp_field", "epi_correction_deformation_field"), + ("out_matrix", "epi_correction_affine_transform"), + ("warped_image", "epi_correction_image_warped"), + ], + ), + ( + b_vectors_rotation, + outputnode, + [("rotated_b_vectors_filename", "rotated_b_vectors")], + ), + ] + + wf.connect(connections) + + return wf + + +def perform_dwi_epi_correction( + base_dir: str, + delete_cache: bool = False, + name="perform_dwi_epi_correction", +) -> Workflow: + """Step 2 of EPI pipeline.""" + import nipype.interfaces.ants as ants + import nipype.interfaces.fsl as fsl + import nipype.interfaces.utility as niu + import nipype.pipeline.engine as pe + + from .dwi_preprocessing_using_t1_utils import ( + ants_apply_transforms, + delete_temp_dirs, + ) + inputnode = pe.Node( + niu.IdentityInterface( + fields=["t1_filename", "dwi_filename", "merged_transforms"] + ), + name="inputnode", + ) + split_dwi_volumes = pe.Node(fsl.Split(dimension="t"), name="split_dwi_volumes") # the nipype function is not used since the results it gives are not # as good as the ones we get by using the command directly. apply_transform_image = pe.MapNode( @@ -293,7 +494,6 @@ def epi_pipeline( iterfield=["deformationField"], name="jacobian", ) - jacobian.inputs.imageDimension = 3 jacobian.inputs.outputImage = "Jacobian_image.nii.gz" @@ -303,11 +503,13 @@ def epi_pipeline( name="ModulateDWIs", ) - thres = pe.MapNode( - fsl.Threshold(thresh=0.0), iterfield=["in_file"], name="RemoveNegative" + threshold_negative = pe.MapNode( + fsl.Threshold(thresh=0.0), + iterfield=["in_file"], + name="threshold_negative", ) - merge = pe.Node(fsl.Merge(dimension="t"), name="MergeDWIs") + merge_dwi_volumes = pe.Node(fsl.Merge(dimension="t"), name="merge_dwi_volumes") # Delete the temporary directory that takes too much place delete_warp_field_tmp = pe.Node( @@ -327,74 +529,35 @@ def epi_pipeline( ] outputnode = pe.Node( - niu.IdentityInterface( - fields=[ - "DWI_2_T1_Coregistration_matrix", - "epi_correction_deformation_field", - "epi_correction_affine_transform", - "epi_correction_image_warped", - "DWIs_epicorrected", - "warp_epi", - "out_bvec", - ] - ), + niu.IdentityInterface(fields=["epi_corrected_dwi_image"]), name="outputnode", ) - wf = pe.Workflow(name="epi_pipeline") - # fmt: off - wf.connect( - [ - (inputnode, split, [("DWI", "in_file")]), - (split, pick_ref, [("out_files", "inlist")]), - (pick_ref, flirt_b0_2_t1, [("out", "in_file")]), - (inputnode, flirt_b0_2_t1, [("T1", "reference")]), - (inputnode, rot_bvec, [("bvec", "in_bvec")]), - (flirt_b0_2_t1, expend_matrix, [("out_matrix_file", "in_matrix")]), - (inputnode, expend_matrix, [("bvec", "in_bvec")]), - (expend_matrix, rot_bvec, [("out_matrix_list", "in_matrix")]), - (inputnode, ants_registration, [("T1", "fixed_image")]), - (flirt_b0_2_t1, ants_registration, [("out_file", "moving_image")]), - (inputnode, c3d_flirt2ants, [("T1", "reference_file")]), - (pick_ref, c3d_flirt2ants, [("out", "source_file")]), - (flirt_b0_2_t1, c3d_flirt2ants, [("out_matrix_file", "transform_file")]), - (c3d_flirt2ants, change_transform, [("itk_transform", "input_affine_file")]), - (change_transform, merge_transform, [("updated_affine_file", "in1")]), - (ants_registration, merge_transform, [("out_matrix", "in2")]), - (ants_registration, merge_transform, [("forward_warp_field", "in3")]), - (inputnode, apply_transform_image, [("T1", "fixed_image")]), - (split, apply_transform_image, [("out_files", "moving_image")]), - (merge_transform, apply_transform_image, [("out", "transforms")]), - - (inputnode, apply_transform_field, [("T1", "fixed_image")]), - (split, apply_transform_field, [("out_files", "moving_image")]), - (merge_transform, apply_transform_field, [("out", "transforms")]), - - (apply_transform_field, jacobian, [("warped_image", "deformationField")]), - (apply_transform_image, jacmult, [("warped_image", "operand_files")]), - (jacobian, jacmult, [("jacobian_image", "in_file")]), - (jacmult, thres, [("out_file", "in_file")]), - (thres, merge, [("out_file", "in_files")]), + wf = pe.Workflow(name=name) - (merge, outputnode, [("merged_file", "DWIs_epicorrected")]), - (flirt_b0_2_t1, outputnode, [("out_matrix_file", "DWI_2_T1_Coregistration_matrix")]), - (ants_registration, outputnode, [("forward_warp_field", "epi_correction_deformation_field"), - ("out_matrix", "epi_correction_affine_transform"), - ("warped_image", "epi_correction_image_warped")]), - (merge_transform, outputnode, [("out", "warp_epi")]), - (rot_bvec, outputnode, [("out_file", "out_bvec")]), + connections = [ + (inputnode, split_dwi_volumes, [("dwi_filename", "in_file")]), + (inputnode, apply_transform_image, [("t1_filename", "fixed_image")]), + (split_dwi_volumes, apply_transform_image, [("out_files", "moving_image")]), + (inputnode, apply_transform_image, [("merged_transforms", "transforms")]), + (inputnode, apply_transform_field, [("t1_filename", "fixed_image")]), + (split_dwi_volumes, apply_transform_field, [("out_files", "moving_image")]), + (inputnode, apply_transform_field, [("merged_transforms", "transforms")]), + (apply_transform_field, jacobian, [("warped_image", "deformationField")]), + (apply_transform_image, jacmult, [("warped_image", "operand_files")]), + (jacobian, jacmult, [("jacobian_image", "in_file")]), + (jacmult, threshold_negative, [("out_file", "in_file")]), + (threshold_negative, merge_dwi_volumes, [("out_file", "in_files")]), + (merge_dwi_volumes, outputnode, [("merged_file", "epi_corrected_dwi_image")]), + ] + if delete_cache: + connections += [ + (merge_dwi_volumes, delete_warp_field_tmp, [("merged_file", "checkpoint")]) + ] + wf.connect(connections) - ] - ) - if delete_cache: - wf.connect( - [ - (merge, delete_warp_field_tmp, [("merged_file", "checkpoint")]) - ] - ) - # fmt: on return wf