Add remove capsule (#274)

* lithresh, gcc and mask

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

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"

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)

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")

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',

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,27 +557,23 @@ 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(
         reg.RegAladin(),
         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(

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