From 4e94067abf7e0651a9ff81ad79d3bfbde1bf2282 Mon Sep 17 00:00:00 2001 From: David Meunier Date: Wed, 4 Dec 2024 14:03:28 +0100 Subject: [PATCH 01/19] lithresh --- macapype/nodes/prepare.py | 34 +++++++++++++++++++++++++++++++ macapype/pipelines/prepare.py | 38 ++++++++++++++++++++++++++++++++++- 2 files changed, 71 insertions(+), 1 deletion(-) diff --git a/macapype/nodes/prepare.py b/macapype/nodes/prepare.py index 19feeb597..f112c763d 100644 --- a/macapype/nodes/prepare.py +++ b/macapype/nodes/prepare.py @@ -445,6 +445,40 @@ def padding_cropped_img(cropped_img_file, orig_img_file, indiv_crop): return padded_img_file +def apply_li_thresh(orig_img_file): + + import os + from skimage.filters import threshold_li + import nibabel as nib + # import numpy as np + + from nipype.utils.filemanip import split_filename as split_f + + # orig image + orig_img = nib.load(orig_img_file) + + data_orig = orig_img.get_data() + header_orig = orig_img.header + affine_orig = orig_img.affine + + print("Orig img shape:", data_orig.shape) + + img_lithresh = threshold_li(data_orig) + + # saving lithresh image + fpath, fname, ext = split_f(orig_img_file) + lithr_img_file = os.path.abspath(fname + "_lithresh" + ext) + img_thresh = nib.Nifti1Image(padded_img_data, affine=affine_orig, + header=header_orig) + nib.save(img_lithresh, lithr_img_file) + + return lithr_img_file + + + + return thresh_img_file + + if __name__ == '__main__': data_path = "/hpc/meca/data/Macaques/Macaque_hiphop/results/ucdavis" diff --git a/macapype/pipelines/prepare.py b/macapype/pipelines/prepare.py index cfcad427a..776c7c2bc 100644 --- a/macapype/pipelines/prepare.py +++ b/macapype/pipelines/prepare.py @@ -13,7 +13,7 @@ from ..utils.utils_nodes import NodeParams from ..utils.misc import parse_key -from ..nodes.prepare import average_align +from ..nodes.prepare import average_align, apply_li_thresh from ..nodes.register import pad_zero_mri @@ -70,6 +70,34 @@ def _create_avg_reorient_pipeline(name="avg_reorient_pipe", params={}): return reorient_pipe + +def _create_remove_capsule_pipeline(name="remove_capsule_pipe", params={}): + + """ + By david: + li_thresholding + gcc + """ + # creating pipeline + remove_caps_pipe = pe.Workflow(name=name) + + # Creating input node + inputnode = pe.Node( + niu.IdentityInterface(fields=['av_img', 'indiv_params']), + name='inputnode' + ) + + li_thresh = pe.Node( + niu.Function(input_names=['orig_img_file'], + output_names=['lithr_img_file'], + function=apply_li_thresh), + name="li_thresholding") + + remove_caps_pipe.connect(inputnode, 'av_img', li_thresh, 'orig_img') + + return remove_caps_pipe + + ############################################################################### # main pipeline: "short_preparation_pipe" ############################################################################### @@ -392,6 +420,14 @@ def create_short_preparation_pipe(params, params_template={}, align_T2_on_T1, "out_file", outputnode, 'native_T2') + if "remove_capsule_pipe" in params: + + remove_capsule_pipe = _create_remove_capsule_pipeline( + params["remove_capsule_pipe"]) + + 0/0 + + # register T1 to stereo image crop_aladin_pipe = create_crop_aladin_pipe( "crop_aladin_pipe", From cf93e708063546ff1ec9f5f3078f487ca4fac514 Mon Sep 17 00:00:00 2001 From: David Meunier Date: Wed, 4 Dec 2024 14:42:29 +0100 Subject: [PATCH 02/19] li_thr, gcc and mask --- macapype/nodes/prepare.py | 21 +++++++++++++-------- macapype/pipelines/prepare.py | 24 +++++++++++++++++++++++- 2 files changed, 36 insertions(+), 9 deletions(-) diff --git a/macapype/nodes/prepare.py b/macapype/nodes/prepare.py index f112c763d..b6768156a 100644 --- a/macapype/nodes/prepare.py +++ b/macapype/nodes/prepare.py @@ -444,31 +444,39 @@ def padding_cropped_img(cropped_img_file, orig_img_file, indiv_crop): return padded_img_file - def apply_li_thresh(orig_img_file): import os from skimage.filters import threshold_li import nibabel as nib - # import numpy as np + import numpy as np from nipype.utils.filemanip import split_filename as split_f # orig image orig_img = nib.load(orig_img_file) - data_orig = orig_img.get_data() + data_orig = orig_img.get_fdata() header_orig = orig_img.header affine_orig = orig_img.affine print("Orig img shape:", data_orig.shape) - img_lithresh = threshold_li(data_orig) + print("Min/max values orig img :", np.min(data_orig), np.max(data_orig)) + + lithresh = threshold_li(data_orig) + + img_lithresh_data = data_orig.astype(bool) + + img_lithresh_data[data_orig <= lithresh] = 0 + img_lithresh_data[lithresh < data_orig] = 1 + + print("Min/max values after lithresh:", np.min(img_lithresh_data), np.max(img_lithresh_data)) # saving lithresh image fpath, fname, ext = split_f(orig_img_file) lithr_img_file = os.path.abspath(fname + "_lithresh" + ext) - img_thresh = nib.Nifti1Image(padded_img_data, affine=affine_orig, + img_lithresh = nib.Nifti1Image(img_lithresh_data, affine=affine_orig, header=header_orig) nib.save(img_lithresh, lithr_img_file) @@ -476,9 +484,6 @@ def apply_li_thresh(orig_img_file): - return thresh_img_file - - if __name__ == '__main__': data_path = "/hpc/meca/data/Macaques/Macaque_hiphop/results/ucdavis" diff --git a/macapype/pipelines/prepare.py b/macapype/pipelines/prepare.py index 776c7c2bc..ed559bc1b 100644 --- a/macapype/pipelines/prepare.py +++ b/macapype/pipelines/prepare.py @@ -16,6 +16,7 @@ from ..nodes.prepare import average_align, apply_li_thresh from ..nodes.register import pad_zero_mri +from ..nodes.surface import keep_gcc # should be in nipype code directly from ..nodes.prepare import Refit @@ -87,13 +88,34 @@ def _create_remove_capsule_pipeline(name="remove_capsule_pipe", params={}): name='inputnode' ) + # lithresholding li_thresh = pe.Node( niu.Function(input_names=['orig_img_file'], output_names=['lithr_img_file'], function=apply_li_thresh), name="li_thresholding") - remove_caps_pipe.connect(inputnode, 'av_img', li_thresh, 'orig_img') + remove_caps_pipe.connect(inputnode, 'av_img', li_thresh, 'orig_img_file') + + # gcc + + gcc_mask = pe.Node( + niu.Function(input_names=['nii_file'], + output_names=['gcc_nii_file'], + function=keep_gcc), + name="gcc_mask") + + + remove_caps_pipe.connect(li_thresh, 'lithr_img_file', + gcc_mask, 'nii_file') + + # masking original_image + mask_capsule = pe.Node(fsl.ApplyMask(), name='mask_capsule') + + remove_caps_pipe.connect(inputnode, 'av_img', + mask_capsule, 'in_file') + remove_caps_pipe.connect(gcc_mask, 'gcc_nii_file', + mask_capsule, 'mask_file') return remove_caps_pipe From c84ac542923f60d714aa27c9326209a205889685 Mon Sep 17 00:00:00 2001 From: David Meunier Date: Wed, 4 Dec 2024 14:51:46 +0100 Subject: [PATCH 03/19] connect --- macapype/pipelines/prepare.py | 25 +++++++++++++++++++------ 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/macapype/pipelines/prepare.py b/macapype/pipelines/prepare.py index ed559bc1b..499267eb7 100644 --- a/macapype/pipelines/prepare.py +++ b/macapype/pipelines/prepare.py @@ -447,8 +447,15 @@ def create_short_preparation_pipe(params, params_template={}, remove_capsule_pipe = _create_remove_capsule_pipeline( params["remove_capsule_pipe"]) - 0/0 + if "avg_reorient_pipe" in params.keys(): + data_preparation_pipe.connect( + av_T1, 'outputnode.std_img', + remove_capsule_pipe, 'inputnode.avg_img') + else: + data_preparation_pipe.connect( + av_T1, 'avg_img', + remove_capsule_pipe, 'inputnode.avg_img') # register T1 to stereo image crop_aladin_pipe = create_crop_aladin_pipe( @@ -461,15 +468,21 @@ def create_short_preparation_pipe(params, params_template={}, crop_aladin_pipe, 'inputnode.native_T1') else: - if "avg_reorient_pipe" in params.keys(): + if "remove_capsule_pipe" in params.keys(): data_preparation_pipe.connect( - av_T1, 'outputnode.std_img', + remove_capsule_pipe, 'mask_capsule.out_file', crop_aladin_pipe, 'inputnode.native_T1') else: - data_preparation_pipe.connect( - av_T1, 'avg_img', - crop_aladin_pipe, 'inputnode.native_T1') + if "avg_reorient_pipe" in params.keys(): + data_preparation_pipe.connect( + av_T1, 'outputnode.std_img', + crop_aladin_pipe, 'inputnode.native_T1') + + else: + data_preparation_pipe.connect( + av_T1, 'avg_img', + crop_aladin_pipe, 'inputnode.native_T1') data_preparation_pipe.connect( inputnode, 'indiv_params', From afa2437a6f05e774ee1fb964afdc8a7bea150df7 Mon Sep 17 00:00:00 2001 From: David Meunier Date: Wed, 4 Dec 2024 14:55:04 +0100 Subject: [PATCH 04/19] params --- macapype/pipelines/prepare.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/macapype/pipelines/prepare.py b/macapype/pipelines/prepare.py index 499267eb7..0dc533c1a 100644 --- a/macapype/pipelines/prepare.py +++ b/macapype/pipelines/prepare.py @@ -73,7 +73,6 @@ def _create_avg_reorient_pipeline(name="avg_reorient_pipe", params={}): def _create_remove_capsule_pipeline(name="remove_capsule_pipe", params={}): - """ By david: li_thresholding @@ -105,7 +104,6 @@ def _create_remove_capsule_pipeline(name="remove_capsule_pipe", params={}): function=keep_gcc), name="gcc_mask") - remove_caps_pipe.connect(li_thresh, 'lithr_img_file', gcc_mask, 'nii_file') @@ -445,7 +443,7 @@ def create_short_preparation_pipe(params, params_template={}, if "remove_capsule_pipe" in params: remove_capsule_pipe = _create_remove_capsule_pipeline( - params["remove_capsule_pipe"]) + params = params["remove_capsule_pipe"]) if "avg_reorient_pipe" in params.keys(): data_preparation_pipe.connect( From 1b803382b82b17f75f1d25e96e852a70cb07cee7 Mon Sep 17 00:00:00 2001 From: David Meunier Date: Wed, 4 Dec 2024 14:57:20 +0100 Subject: [PATCH 05/19] orig_img --- macapype/pipelines/prepare.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/macapype/pipelines/prepare.py b/macapype/pipelines/prepare.py index 0dc533c1a..9f01f4f33 100644 --- a/macapype/pipelines/prepare.py +++ b/macapype/pipelines/prepare.py @@ -83,7 +83,7 @@ def _create_remove_capsule_pipeline(name="remove_capsule_pipe", params={}): # Creating input node inputnode = pe.Node( - niu.IdentityInterface(fields=['av_img', 'indiv_params']), + niu.IdentityInterface(fields=['orig_img', 'indiv_params']), name='inputnode' ) @@ -94,7 +94,7 @@ def _create_remove_capsule_pipeline(name="remove_capsule_pipe", params={}): function=apply_li_thresh), name="li_thresholding") - remove_caps_pipe.connect(inputnode, 'av_img', li_thresh, 'orig_img_file') + remove_caps_pipe.connect(inputnode, 'orig_img', li_thresh, 'orig_img_file') # gcc @@ -110,7 +110,7 @@ def _create_remove_capsule_pipeline(name="remove_capsule_pipe", params={}): # masking original_image mask_capsule = pe.Node(fsl.ApplyMask(), name='mask_capsule') - remove_caps_pipe.connect(inputnode, 'av_img', + remove_caps_pipe.connect(inputnode, 'orig_img', mask_capsule, 'in_file') remove_caps_pipe.connect(gcc_mask, 'gcc_nii_file', mask_capsule, 'mask_file') @@ -441,19 +441,18 @@ def create_short_preparation_pipe(params, params_template={}, outputnode, 'native_T2') if "remove_capsule_pipe" in params: - remove_capsule_pipe = _create_remove_capsule_pipeline( - params = params["remove_capsule_pipe"]) + params=params["remove_capsule_pipe"]) if "avg_reorient_pipe" in params.keys(): data_preparation_pipe.connect( av_T1, 'outputnode.std_img', - remove_capsule_pipe, 'inputnode.avg_img') + remove_capsule_pipe, 'inputnode.orig_img') else: data_preparation_pipe.connect( av_T1, 'avg_img', - remove_capsule_pipe, 'inputnode.avg_img') + remove_capsule_pipe, 'inputnode.orig_img') # register T1 to stereo image crop_aladin_pipe = create_crop_aladin_pipe( From 86971a6d0919978a922f1e361ddcf36015fbc3d0 Mon Sep 17 00:00:00 2001 From: David Meunier Date: Thu, 5 Dec 2024 13:14:39 +0100 Subject: [PATCH 06/19] forcing bool type --- macapype/nodes/prepare.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/macapype/nodes/prepare.py b/macapype/nodes/prepare.py index b6768156a..b77c75b6f 100644 --- a/macapype/nodes/prepare.py +++ b/macapype/nodes/prepare.py @@ -444,6 +444,7 @@ def padding_cropped_img(cropped_img_file, orig_img_file, indiv_crop): return padded_img_file + def apply_li_thresh(orig_img_file): import os @@ -466,7 +467,7 @@ def apply_li_thresh(orig_img_file): lithresh = threshold_li(data_orig) - img_lithresh_data = data_orig.astype(bool) + img_lithresh_data = data_orig.astype(np.int16) img_lithresh_data[data_orig <= lithresh] = 0 img_lithresh_data[lithresh < data_orig] = 1 From 6cf98e5c1b91c9d3f0b3a554b658442e8d137dbf Mon Sep 17 00:00:00 2001 From: David Meunier Date: Thu, 5 Dec 2024 13:50:08 +0100 Subject: [PATCH 07/19] type np.int16 --- macapype/nodes/surface.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/macapype/nodes/surface.py b/macapype/nodes/surface.py index 939fddc89..77cb1b743 100644 --- a/macapype/nodes/surface.py +++ b/macapype/nodes/surface.py @@ -26,7 +26,8 @@ def getLargestCC(segmentation): # numpy data[data > 0] = 1 - binary = np.array(data, dtype="bool") + + binary = np.array(data, dtype=np.int16) # skimage GCC new_data = getLargestCC(binary) From e2f79aa2a7ca71f2ac9c640c27e18bb65d7fcf15 Mon Sep 17 00:00:00 2001 From: David Meunier Date: Thu, 5 Dec 2024 14:35:47 +0100 Subject: [PATCH 08/19] orig_native_T1 --- macapype/pipelines/prepare.py | 19 ++++++++++++ macapype/pipelines/register.py | 56 +++++----------------------------- 2 files changed, 27 insertions(+), 48 deletions(-) diff --git a/macapype/pipelines/prepare.py b/macapype/pipelines/prepare.py index 9f01f4f33..426b1d22d 100644 --- a/macapype/pipelines/prepare.py +++ b/macapype/pipelines/prepare.py @@ -464,12 +464,20 @@ def create_short_preparation_pipe(params, params_template={}, crop_T1, "roi_file", crop_aladin_pipe, 'inputnode.native_T1') + data_preparation_pipe.connect( + crop_T1, "roi_file", + crop_aladin_pipe, 'inputnode.orig_native_T1') + else: if "remove_capsule_pipe" in params.keys(): data_preparation_pipe.connect( remove_capsule_pipe, 'mask_capsule.out_file', crop_aladin_pipe, 'inputnode.native_T1') + data_preparation_pipe.connect( + av_T1, 'outputnode.std_img', + crop_aladin_pipe, 'inputnode.orig_native_T1') + else: if "avg_reorient_pipe" in params.keys(): data_preparation_pipe.connect( @@ -481,6 +489,17 @@ def create_short_preparation_pipe(params, params_template={}, av_T1, 'avg_img', crop_aladin_pipe, 'inputnode.native_T1') + if "avg_reorient_pipe" in params.keys(): + data_preparation_pipe.connect( + av_T1, 'outputnode.std_img', + crop_aladin_pipe, 'inputnode.orig_native_T1') + + else: + data_preparation_pipe.connect( + av_T1, 'avg_img', + crop_aladin_pipe, 'inputnode.orig_native_T1') + + data_preparation_pipe.connect( inputnode, 'indiv_params', crop_aladin_pipe, 'inputnode.indiv_params') diff --git a/macapype/pipelines/register.py b/macapype/pipelines/register.py index 707f40a46..96bc5b2d3 100644 --- a/macapype/pipelines/register.py +++ b/macapype/pipelines/register.py @@ -495,7 +495,7 @@ def create_crop_aladin_pipe(name="crop_aladin_pipe", params={}): # creating inputnode inputnode = pe.Node( - niu.IdentityInterface(fields=['native_T1', + niu.IdentityInterface(fields=['native_T1', 'orig_native_T1', 'stereo_template_T1', 'indiv_params']), name='inputnode') @@ -506,33 +506,15 @@ def create_crop_aladin_pipe(name="crop_aladin_pipe", params={}): 'native_to_stereo_trans']), name='outputnode') - if "pre_crop_z_T1" in params.keys(): - - print('pre_crop_z_T1') - pre_crop_z_T1 = NodeParams( - fsl.RobustFOV(), - params=parse_key(params, "pre_crop_z_T1"), - name='pre_crop_z_T1') - - reg_pipe.connect( - inputnode, 'native_T1', - pre_crop_z_T1, 'in_file') - # align T1 on template reg_T1_on_template = NodeParams( reg.RegAladin(), params=parse_key(params, "reg_T1_on_template"), name='reg_T1_on_template') - if "pre_crop_z_T1" in params.keys(): - reg_pipe.connect( - pre_crop_z_T1, "out_roi", - reg_T1_on_template, "flo_file") - else: - - reg_pipe.connect( - inputnode, 'native_T1', - reg_T1_on_template, "flo_file") + reg_pipe.connect( + inputnode, 'native_T1', + reg_T1_on_template, "flo_file") reg_pipe.connect( inputnode, 'stereo_template_T1', @@ -573,23 +555,7 @@ def create_crop_aladin_pipe(name="crop_aladin_pipe", params={}): reg_pipe.connect(reg_T1_on_template2, 'aff_file', compose_transfo, "comp_input") - # crop_z_T1 - if "crop_z_T1" in params.keys(): - crop_z_T1 = NodeParams( - fsl.RobustFOV(), - params=parse_key(params, "crop_z_T1"), - name='crop_z_T1') - - if "reg_T1_on_template2" in params.keys(): - reg_pipe.connect( - reg_T1_on_template2, 'res_file', - crop_z_T1, 'in_file') - - else: - reg_pipe.connect( - reg_T1_on_template, 'res_file', - crop_z_T1, 'in_file') - + ### padding image to have zeros and non nans pad_image_T1 = pe.Node( ants.utils.ImageMath(), name="pad_image_T1") @@ -598,15 +564,9 @@ def create_crop_aladin_pipe(name="crop_aladin_pipe", params={}): pad_image_T1.inputs.operation = "PadImage" pad_image_T1.inputs.op2 = '200' - if "pre_crop_z_T1" in params.keys(): - reg_pipe.connect( - pre_crop_z_T1, 'out_roi', - pad_image_T1, "op1") - - else: - reg_pipe.connect( - inputnode, 'native_T1', - pad_image_T1, "op1") + reg_pipe.connect( + inputnode, 'orig_native_T1', + pad_image_T1, "op1") # resampling using transfo on much bigger image reg_resample_T1 = pe.Node( From fc30c1c0e01edd8c1257917e8d88467711f163fa Mon Sep 17 00:00:00 2001 From: David Meunier Date: Thu, 5 Dec 2024 14:39:16 +0100 Subject: [PATCH 09/19] connect --- macapype/pipelines/prepare.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/macapype/pipelines/prepare.py b/macapype/pipelines/prepare.py index 426b1d22d..ac9bf21a1 100644 --- a/macapype/pipelines/prepare.py +++ b/macapype/pipelines/prepare.py @@ -469,15 +469,12 @@ def create_short_preparation_pipe(params, params_template={}, crop_aladin_pipe, 'inputnode.orig_native_T1') else: + # connect native_T1 if "remove_capsule_pipe" in params.keys(): data_preparation_pipe.connect( remove_capsule_pipe, 'mask_capsule.out_file', crop_aladin_pipe, 'inputnode.native_T1') - data_preparation_pipe.connect( - av_T1, 'outputnode.std_img', - crop_aladin_pipe, 'inputnode.orig_native_T1') - else: if "avg_reorient_pipe" in params.keys(): data_preparation_pipe.connect( @@ -489,6 +486,7 @@ def create_short_preparation_pipe(params, params_template={}, av_T1, 'avg_img', crop_aladin_pipe, 'inputnode.native_T1') + # connect orig_native_T1 if "avg_reorient_pipe" in params.keys(): data_preparation_pipe.connect( av_T1, 'outputnode.std_img', @@ -499,7 +497,6 @@ def create_short_preparation_pipe(params, params_template={}, av_T1, 'avg_img', crop_aladin_pipe, 'inputnode.orig_native_T1') - data_preparation_pipe.connect( inputnode, 'indiv_params', crop_aladin_pipe, 'inputnode.indiv_params') From cb1daedcf4660118180c403ec0f280a18d459e73 Mon Sep 17 00:00:00 2001 From: David Meunier Date: Thu, 5 Dec 2024 15:04:40 +0100 Subject: [PATCH 10/19] adding img_type to fast_T1 and fast_T2 --- macapype/pipelines/full_pipelines.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/macapype/pipelines/full_pipelines.py b/macapype/pipelines/full_pipelines.py index 387b95b55..594f10b70 100644 --- a/macapype/pipelines/full_pipelines.py +++ b/macapype/pipelines/full_pipelines.py @@ -1054,6 +1054,7 @@ def create_full_ants_subpipes( fast_T1.inputs.output_biascorrected = True fast_T1.inputs.output_biasfield = True + fast_T1.inputs.img_type = 1 seg_pipe.connect(data_preparation_pipe, 'outputnode.preproc_T1', fast_T1, "in_files") @@ -1070,6 +1071,7 @@ def create_full_ants_subpipes( fast_T2.inputs.output_biascorrected = True fast_T2.inputs.output_biasfield = True + fast_T2.inputs.img_type = 1 seg_pipe.connect(data_preparation_pipe, 'outputnode.preproc_T2', fast_T2, "in_files") From a0ffb94cb538151ded8ff2afb550e3950685e63a Mon Sep 17 00:00:00 2001 From: David Meunier Date: Thu, 5 Dec 2024 16:32:36 +0100 Subject: [PATCH 11/19] itk_debias --- macapype/nodes/correct_bias.py | 36 ++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/macapype/nodes/correct_bias.py b/macapype/nodes/correct_bias.py index a1b7e19d5..d50d4740c 100644 --- a/macapype/nodes/correct_bias.py +++ b/macapype/nodes/correct_bias.py @@ -217,3 +217,39 @@ def _list_outputs(self): outputs["t2_debiased_brain_file"] = os.path.abspath( t2_fname + self.inputs.os + "_brain.nii.gz") return outputs + +# ############## simpleITK debias + + +def itk_debias(img_file): + + import os + import SimpleITK as sitk + + from nipype.utils.filemanip import split_filename as split_f + + # Read the input image + input_image = sitk.ReadImage(img_file, sitk.sitkFloat32) + + # Initialize the N4 bias field correction filter + n4_filter = sitk.N4BiasFieldCorrectionImageFilter() + + input_image_mask = sitk.LiThreshold(input_image, 0, 1) + + # the output image + corrected_image = n4_filter.Execute(input_image, input_image_mask) + + # get the bias field + bias_field_image = n4_filter.GetLogBiasFieldAsImage(input_image) + + # Save the corrected image and bias + fpath, fname, ext = split_f(img_file) + cor_img_file = os.path.abspath(fname + "_debias" + ext) + + bias_img_file = os.path.abspath(fname + "_bias" + ext) + + sitk.WriteImage(corrected_image, cor_img_file) + + sitk.WriteImage(bias_field_image, bias_img_file) + + return cor_img_file, bias_img_file From 18e4bf3abc82e263684a3098fe49e2b2874e915c Mon Sep 17 00:00:00 2001 From: David Meunier Date: Thu, 5 Dec 2024 17:42:04 +0100 Subject: [PATCH 12/19] fast T1 --- macapype/pipelines/full_pipelines.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/macapype/pipelines/full_pipelines.py b/macapype/pipelines/full_pipelines.py index 594f10b70..8383f161a 100644 --- a/macapype/pipelines/full_pipelines.py +++ b/macapype/pipelines/full_pipelines.py @@ -1071,7 +1071,7 @@ def create_full_ants_subpipes( fast_T2.inputs.output_biascorrected = True fast_T2.inputs.output_biasfield = True - fast_T2.inputs.img_type = 1 + fast_T2.inputs.img_type = 2 seg_pipe.connect(data_preparation_pipe, 'outputnode.preproc_T2', fast_T2, "in_files") @@ -1878,6 +1878,7 @@ def create_full_T1_ants_subpipes(params_template, params_template_stereo, fast_T1.inputs.output_biascorrected = True fast_T1.inputs.output_biasfield = True + fast_T1.inputs.img_type = 1 seg_pipe.connect(data_preparation_pipe, 'outputnode.preproc_T1', fast_T1, "in_files") From 5d0b23f26c2965d1b85f06632e83e4684822f890 Mon Sep 17 00:00:00 2001 From: David Meunier Date: Fri, 6 Dec 2024 10:14:05 +0100 Subject: [PATCH 13/19] no brain-slam --- requirements.txt | 1 - 1 file changed, 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 0707769f0..9347415c0 100644 --- a/requirements.txt +++ b/requirements.txt @@ -5,4 +5,3 @@ pybids scikit-image nibabel numpy -brain-slam From 40fb19e5fc4eba8019ac37ffa9dcbfdffeeb2922 Mon Sep 17 00:00:00 2001 From: David Meunier Date: Fri, 6 Dec 2024 10:20:28 +0100 Subject: [PATCH 14/19] forcing version --- requirements.txt | 3 +++ 1 file changed, 3 insertions(+) diff --git a/requirements.txt b/requirements.txt index 9347415c0..9bc3d6149 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,7 +1,10 @@ nipype +rdflib==6.3.1 +pandas==2.2.3 nilearn networkx pybids scikit-image nibabel numpy +SimpleITK From 46d75aadc161cbdc7210939701a7284342ebddd4 Mon Sep 17 00:00:00 2001 From: David Meunier Date: Fri, 6 Dec 2024 10:53:16 +0100 Subject: [PATCH 15/19] setting best paramters --- macapype/nodes/correct_bias.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/macapype/nodes/correct_bias.py b/macapype/nodes/correct_bias.py index d50d4740c..2bbbf991d 100644 --- a/macapype/nodes/correct_bias.py +++ b/macapype/nodes/correct_bias.py @@ -234,6 +234,11 @@ def itk_debias(img_file): # Initialize the N4 bias field correction filter n4_filter = sitk.N4BiasFieldCorrectionImageFilter() + # best paramters for marmoset Randy petra + n4_filter.SetMaximumNumberOfIterations([200, 100, 50, 25]) + n4_filter.SetConvergenceThreshold(0.0001) + n4_filter.SetWienerFilterNoise(0.01) + input_image_mask = sitk.LiThreshold(input_image, 0, 1) # the output image From 16ef5c19030a790f61a8c83ab53a680035978699 Mon Sep 17 00:00:00 2001 From: David Meunier Date: Fri, 6 Dec 2024 14:00:51 +0100 Subject: [PATCH 16/19] mask_capsule in crop_aladin_pipe --- macapype/pipelines/prepare.py | 79 +--------------------------------- macapype/pipelines/register.py | 72 +++++++++++++++++++++++++++++-- 2 files changed, 69 insertions(+), 82 deletions(-) diff --git a/macapype/pipelines/prepare.py b/macapype/pipelines/prepare.py index ac9bf21a1..ac04a0f19 100644 --- a/macapype/pipelines/prepare.py +++ b/macapype/pipelines/prepare.py @@ -13,10 +13,9 @@ from ..utils.utils_nodes import NodeParams from ..utils.misc import parse_key -from ..nodes.prepare import average_align, apply_li_thresh +from ..nodes.prepare import average_align from ..nodes.register import pad_zero_mri -from ..nodes.surface import keep_gcc # should be in nipype code directly from ..nodes.prepare import Refit @@ -72,51 +71,6 @@ def _create_avg_reorient_pipeline(name="avg_reorient_pipe", params={}): return reorient_pipe -def _create_remove_capsule_pipeline(name="remove_capsule_pipe", params={}): - """ - By david: - li_thresholding - gcc - """ - # creating pipeline - remove_caps_pipe = pe.Workflow(name=name) - - # Creating input node - inputnode = pe.Node( - niu.IdentityInterface(fields=['orig_img', 'indiv_params']), - name='inputnode' - ) - - # lithresholding - li_thresh = pe.Node( - niu.Function(input_names=['orig_img_file'], - output_names=['lithr_img_file'], - function=apply_li_thresh), - name="li_thresholding") - - remove_caps_pipe.connect(inputnode, 'orig_img', li_thresh, 'orig_img_file') - - # gcc - - gcc_mask = pe.Node( - niu.Function(input_names=['nii_file'], - output_names=['gcc_nii_file'], - function=keep_gcc), - name="gcc_mask") - - remove_caps_pipe.connect(li_thresh, 'lithr_img_file', - gcc_mask, 'nii_file') - - # masking original_image - mask_capsule = pe.Node(fsl.ApplyMask(), name='mask_capsule') - - remove_caps_pipe.connect(inputnode, 'orig_img', - mask_capsule, 'in_file') - remove_caps_pipe.connect(gcc_mask, 'gcc_nii_file', - mask_capsule, 'mask_file') - - return remove_caps_pipe - ############################################################################### # main pipeline: "short_preparation_pipe" @@ -440,20 +394,6 @@ def create_short_preparation_pipe(params, params_template={}, align_T2_on_T1, "out_file", outputnode, 'native_T2') - if "remove_capsule_pipe" in params: - remove_capsule_pipe = _create_remove_capsule_pipeline( - params=params["remove_capsule_pipe"]) - - if "avg_reorient_pipe" in params.keys(): - data_preparation_pipe.connect( - av_T1, 'outputnode.std_img', - remove_capsule_pipe, 'inputnode.orig_img') - - else: - data_preparation_pipe.connect( - av_T1, 'avg_img', - remove_capsule_pipe, 'inputnode.orig_img') - # register T1 to stereo image crop_aladin_pipe = create_crop_aladin_pipe( "crop_aladin_pipe", @@ -469,23 +409,6 @@ def create_short_preparation_pipe(params, params_template={}, crop_aladin_pipe, 'inputnode.orig_native_T1') else: - # connect native_T1 - if "remove_capsule_pipe" in params.keys(): - data_preparation_pipe.connect( - remove_capsule_pipe, 'mask_capsule.out_file', - crop_aladin_pipe, 'inputnode.native_T1') - - else: - if "avg_reorient_pipe" in params.keys(): - data_preparation_pipe.connect( - av_T1, 'outputnode.std_img', - crop_aladin_pipe, 'inputnode.native_T1') - - else: - data_preparation_pipe.connect( - av_T1, 'avg_img', - crop_aladin_pipe, 'inputnode.native_T1') - # connect orig_native_T1 if "avg_reorient_pipe" in params.keys(): data_preparation_pipe.connect( diff --git a/macapype/pipelines/register.py b/macapype/pipelines/register.py index 96bc5b2d3..e7f3fa024 100644 --- a/macapype/pipelines/register.py +++ b/macapype/pipelines/register.py @@ -15,6 +15,10 @@ NMTSubjectAlign2, NwarpApplyPriors, animal_warper) +from ..nodes.surface import keep_gcc + +from ..nodes.prepare import apply_li_thresh + def create_iterative_register_pipe( template_file, template_brain_file, template_mask_file, gm_prob_file, @@ -489,13 +493,59 @@ def create_reg_seg_pipe(name="reg_seg_pipe"): return reg_pipe +# used for registration +def _create_remove_capsule_pipeline(name="remove_capsule_pipe", params={}): + """ + By david: + li_thresholding + gcc + """ + # creating pipeline + remove_caps_pipe = pe.Workflow(name=name) + + # Creating input node + inputnode = pe.Node( + niu.IdentityInterface(fields=['orig_img', 'indiv_params']), + name='inputnode' + ) + + # lithresholding + li_thresh = pe.Node( + niu.Function(input_names=['orig_img_file'], + output_names=['lithr_img_file'], + function=apply_li_thresh), + name="li_thresholding") + + remove_caps_pipe.connect(inputnode, 'orig_img', li_thresh, 'orig_img_file') + + # gcc + + gcc_mask = pe.Node( + niu.Function(input_names=['nii_file'], + output_names=['gcc_nii_file'], + function=keep_gcc), + name="gcc_mask") + + remove_caps_pipe.connect(li_thresh, 'lithr_img_file', + gcc_mask, 'nii_file') + + # masking original_image + mask_capsule = pe.Node(fsl.ApplyMask(), name='mask_capsule') + + remove_caps_pipe.connect(inputnode, 'orig_img', + mask_capsule, 'in_file') + remove_caps_pipe.connect(gcc_mask, 'gcc_nii_file', + mask_capsule, 'mask_file') + + return remove_caps_pipe + def create_crop_aladin_pipe(name="crop_aladin_pipe", params={}): reg_pipe = pe.Workflow(name=name) # creating inputnode inputnode = pe.Node( - niu.IdentityInterface(fields=['native_T1', 'orig_native_T1', + niu.IdentityInterface(fields=['native_T1', 'stereo_template_T1', 'indiv_params']), name='inputnode') @@ -506,15 +556,29 @@ def create_crop_aladin_pipe(name="crop_aladin_pipe", params={}): 'native_to_stereo_trans']), name='outputnode') + if "remove_capsule_pipe" in params: + remove_capsule_pipe = _create_remove_capsule_pipeline( + params=params["remove_capsule_pipe"]) + + reg_pipe.connect( + inputnode, 'native_T1', + remove_capsule_pipe, 'inputnode.orig_img') + # align T1 on template reg_T1_on_template = NodeParams( reg.RegAladin(), params=parse_key(params, "reg_T1_on_template"), name='reg_T1_on_template') - reg_pipe.connect( - inputnode, 'native_T1', - reg_T1_on_template, "flo_file") + if "remove_capsule_pipe" in params: + reg_pipe.connect( + remove_capsule_pipe, 'mask_capsule.out_file', + reg_T1_on_template, "flo_file") + else: + + reg_pipe.connect( + inputnode, 'native_T1', + reg_T1_on_template, "flo_file") reg_pipe.connect( inputnode, 'stereo_template_T1', From 20fd8f6a8005062b001967f124483ddaf98094a9 Mon Sep 17 00:00:00 2001 From: David Meunier Date: Fri, 6 Dec 2024 14:01:41 +0100 Subject: [PATCH 17/19] no orig --- macapype/pipelines/register.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/macapype/pipelines/register.py b/macapype/pipelines/register.py index e7f3fa024..50ca1a1b9 100644 --- a/macapype/pipelines/register.py +++ b/macapype/pipelines/register.py @@ -629,7 +629,7 @@ def create_crop_aladin_pipe(name="crop_aladin_pipe", params={}): pad_image_T1.inputs.op2 = '200' reg_pipe.connect( - inputnode, 'orig_native_T1', + inputnode, 'native_T1', pad_image_T1, "op1") # resampling using transfo on much bigger image From 61c46c64638dc552cd16e815057d2a088e617935 Mon Sep 17 00:00:00 2001 From: David Meunier Date: Fri, 6 Dec 2024 14:03:19 +0100 Subject: [PATCH 18/19] native_T1 --- macapype/pipelines/prepare.py | 8 ++------ macapype/pipelines/register.py | 1 + 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/macapype/pipelines/prepare.py b/macapype/pipelines/prepare.py index ac04a0f19..298c845f5 100644 --- a/macapype/pipelines/prepare.py +++ b/macapype/pipelines/prepare.py @@ -404,21 +404,17 @@ def create_short_preparation_pipe(params, params_template={}, crop_T1, "roi_file", crop_aladin_pipe, 'inputnode.native_T1') - data_preparation_pipe.connect( - crop_T1, "roi_file", - crop_aladin_pipe, 'inputnode.orig_native_T1') - else: # connect orig_native_T1 if "avg_reorient_pipe" in params.keys(): data_preparation_pipe.connect( av_T1, 'outputnode.std_img', - crop_aladin_pipe, 'inputnode.orig_native_T1') + crop_aladin_pipe, 'inputnode.native_T1') else: data_preparation_pipe.connect( av_T1, 'avg_img', - crop_aladin_pipe, 'inputnode.orig_native_T1') + crop_aladin_pipe, 'inputnode.native_T1') data_preparation_pipe.connect( inputnode, 'indiv_params', diff --git a/macapype/pipelines/register.py b/macapype/pipelines/register.py index 50ca1a1b9..42ca94b5b 100644 --- a/macapype/pipelines/register.py +++ b/macapype/pipelines/register.py @@ -539,6 +539,7 @@ def _create_remove_capsule_pipeline(name="remove_capsule_pipe", params={}): return remove_caps_pipe + def create_crop_aladin_pipe(name="crop_aladin_pipe", params={}): reg_pipe = pe.Workflow(name=name) From d163d282d150290ca93115c8ee94d37a56784ed9 Mon Sep 17 00:00:00 2001 From: David Meunier Date: Fri, 6 Dec 2024 15:03:35 +0100 Subject: [PATCH 19/19] flake8 --- macapype/nodes/prepare.py | 9 +++++---- macapype/pipelines/prepare.py | 1 - macapype/pipelines/register.py | 2 +- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/macapype/nodes/prepare.py b/macapype/nodes/prepare.py index b77c75b6f..1c7d7ee7b 100644 --- a/macapype/nodes/prepare.py +++ b/macapype/nodes/prepare.py @@ -472,19 +472,20 @@ def apply_li_thresh(orig_img_file): img_lithresh_data[data_orig <= lithresh] = 0 img_lithresh_data[lithresh < data_orig] = 1 - print("Min/max values after lithresh:", np.min(img_lithresh_data), np.max(img_lithresh_data)) + print("Min/max values after lithresh:", np.min(img_lithresh_data), + np.max(img_lithresh_data)) # saving lithresh image fpath, fname, ext = split_f(orig_img_file) lithr_img_file = os.path.abspath(fname + "_lithresh" + ext) - img_lithresh = nib.Nifti1Image(img_lithresh_data, affine=affine_orig, - header=header_orig) + img_lithresh = nib.Nifti1Image(img_lithresh_data, + affine=affine_orig, + header=header_orig) nib.save(img_lithresh, lithr_img_file) return lithr_img_file - if __name__ == '__main__': data_path = "/hpc/meca/data/Macaques/Macaque_hiphop/results/ucdavis" diff --git a/macapype/pipelines/prepare.py b/macapype/pipelines/prepare.py index 298c845f5..9eb0b4278 100644 --- a/macapype/pipelines/prepare.py +++ b/macapype/pipelines/prepare.py @@ -71,7 +71,6 @@ def _create_avg_reorient_pipeline(name="avg_reorient_pipe", params={}): return reorient_pipe - ############################################################################### # main pipeline: "short_preparation_pipe" ############################################################################### diff --git a/macapype/pipelines/register.py b/macapype/pipelines/register.py index 42ca94b5b..1c39efe66 100644 --- a/macapype/pipelines/register.py +++ b/macapype/pipelines/register.py @@ -620,7 +620,7 @@ def create_crop_aladin_pipe(name="crop_aladin_pipe", params={}): reg_pipe.connect(reg_T1_on_template2, 'aff_file', compose_transfo, "comp_input") - ### padding image to have zeros and non nans + # padding image to have zeros and non nans pad_image_T1 = pe.Node( ants.utils.ImageMath(), name="pad_image_T1")