diff --git a/macapype/nodes/correct_bias.py b/macapype/nodes/correct_bias.py index a1b7e19d..2bbbf991 100644 --- a/macapype/nodes/correct_bias.py +++ b/macapype/nodes/correct_bias.py @@ -217,3 +217,44 @@ 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() + + # 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 + 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 diff --git a/macapype/nodes/prepare.py b/macapype/nodes/prepare.py index 19feeb59..1c7d7ee7 100644 --- a/macapype/nodes/prepare.py +++ b/macapype/nodes/prepare.py @@ -445,6 +445,47 @@ 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_fdata() + header_orig = orig_img.header + affine_orig = orig_img.affine + + print("Orig img shape:", data_orig.shape) + + 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(np.int16) + + 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_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/nodes/surface.py b/macapype/nodes/surface.py index 939fddc8..77cb1b74 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) diff --git a/macapype/pipelines/full_pipelines.py b/macapype/pipelines/full_pipelines.py index 387b95b5..8383f161 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 = 2 seg_pipe.connect(data_preparation_pipe, 'outputnode.preproc_T2', fast_T2, "in_files") @@ -1876,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") diff --git a/macapype/pipelines/prepare.py b/macapype/pipelines/prepare.py index cfcad427..9eb0b427 100644 --- a/macapype/pipelines/prepare.py +++ b/macapype/pipelines/prepare.py @@ -70,6 +70,7 @@ def _create_avg_reorient_pipeline(name="avg_reorient_pipe", params={}): return reorient_pipe + ############################################################################### # main pipeline: "short_preparation_pipe" ############################################################################### @@ -403,6 +404,7 @@ def create_short_preparation_pipe(params, params_template={}, crop_aladin_pipe, 'inputnode.native_T1') else: + # connect orig_native_T1 if "avg_reorient_pipe" in params.keys(): data_preparation_pipe.connect( av_T1, 'outputnode.std_img', diff --git a/macapype/pipelines/register.py b/macapype/pipelines/register.py index 707f40a4..1c39efe6 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,6 +493,53 @@ 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) @@ -506,17 +557,13 @@ 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') + if "remove_capsule_pipe" in params: + remove_capsule_pipe = _create_remove_capsule_pipeline( + params=params["remove_capsule_pipe"]) reg_pipe.connect( - inputnode, 'native_T1', - pre_crop_z_T1, 'in_file') + inputnode, 'native_T1', + remove_capsule_pipe, 'inputnode.orig_img') # align T1 on template reg_T1_on_template = NodeParams( @@ -524,9 +571,9 @@ def create_crop_aladin_pipe(name="crop_aladin_pipe", params={}): params=parse_key(params, "reg_T1_on_template"), name='reg_T1_on_template') - if "pre_crop_z_T1" in params.keys(): + if "remove_capsule_pipe" in params: reg_pipe.connect( - pre_crop_z_T1, "out_roi", + remove_capsule_pipe, 'mask_capsule.out_file', reg_T1_on_template, "flo_file") else: @@ -573,23 +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") - # 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 +629,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, 'native_T1', + pad_image_T1, "op1") # resampling using transfo on much bigger image reg_resample_T1 = pe.Node( diff --git a/requirements.txt b/requirements.txt index 0707769f..9bc3d614 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,8 +1,10 @@ nipype +rdflib==6.3.1 +pandas==2.2.3 nilearn networkx pybids scikit-image nibabel numpy -brain-slam +SimpleITK