added documentation and minor polishes

bin/2.2_DTIPreProcessing/preProcessing_DTI.py

@@ -10,39 +10,38 @@
 
 
 import nipype.interfaces.fsl as fsl
-import os,sys
+import os, sys
 import nibabel as nii
 import numpy as np
 import applyMICO
 import cv2
 
 # 1) Process MRI
-def applyBET(input_file,frac,radius,outputPath):
-
+def applyBET(input_file: str, frac: float, radius: int, output_path: str) -> str:
+    """
+    Performs brain extraction via the FSL Brain Extraction Tool (BET). Requires an appropriate input file (input_file), the fractional intensity threshold (frac), the head radius (radius) and the output path (output_path).
+    """
     # scale Nifti data by factor 10
     data = nii.load(input_file)
     imgTemp = data.get_data()
     scale = np.eye(4)* 10
     scale[3][3] = 1
     imgTemp = np.flip(imgTemp, 2)
 
-    # imgTemp = np.rot90(imgTemp,2)
-    # imgTemp = np.rot90(imgTemp,2)
-    # imgTemp = np.flip(imgTemp, 0)
     scaledNiiData = nii.Nifti1Image(imgTemp, data.affine * scale)
     hdrIn = scaledNiiData.header
     hdrIn.set_xyzt_units('mm')
     scaledNiiData = nii.as_closest_canonical(scaledNiiData)
     print('Orientation:' + str(nii.aff2axcodes(scaledNiiData.affine)))
 
-    fslPath = os.path.join(os.getcwd(),'fslScaleTemp.nii.gz')
-    nii.save(scaledNiiData, fslPath)
+    fsl_path = os.path.join(os.getcwd(),'fslScaleTemp.nii.gz')
+    nii.save(scaledNiiData, fsl_path)
 
     # extract brain
-    output_file = os.path.join(outputPath, os.path.basename(input_file).split('.')[0] + 'Bet.nii.gz')
-    myBet = fsl.BET(in_file=fslPath, out_file=output_file,frac=frac,radius=radius,robust=True, mask = True)
+    output_file = os.path.join(output_path, os.path.basename(input_file).split('.')[0] + 'Bet.nii.gz')
+    myBet = fsl.BET(in_file=fsl_path, out_file=output_file,frac=frac,radius=radius,robust=True, mask = True)
     myBet.run()
-    os.remove(fslPath)
+    os.remove(fsl_path)
 
 
     # unscale result data by factor 10ˆ(-1)
@@ -59,7 +58,10 @@ def applyBET(input_file,frac,radius,outputPath):
     print('Brain extraction DONE!')
     return output_file
 
-def smoothIMG(input_file,outputPath):
+def smoothIMG(input_file, output_path):
+    """
+    Smoothes image via FSL. Only input and output has do be specified. Parameters are fixed to box shape and to the kernel size of 0.1 voxel.
+    """
     data = nii.load(input_file)
 
     vol = data.get_data()
@@ -70,27 +72,31 @@ def smoothIMG(input_file,outputPath):
     hdrOut.set_xyzt_units('mm')
     output_file = os.path.join(os.path.dirname(input_file),
                                os.path.basename(input_file).split('.')[0] + 'DN.nii.gz')
-    # hdrOut['sform_code'] = 1
     nii.save(unscaledNiiData, output_file)
     input_file = output_file
-    #output_file =  os.path.join(os.path.dirname(input_file),os.path.basename(input_file).split('.')[0] + 'Smooth.nii.gz')
-    output_file = os.path.join(outputPath, os.path.basename(inputFile).split('.')[0] + 'Smooth.nii.gz')
-    myGauss =  fsl.SpatialFilter(in_file=input_file,out_file=output_file,operation='median',kernel_shape='box',kernel_size=0.1)
+    output_file = os.path.join(output_path, os.path.basename(input_file).split('.')[0] + 'Smooth.nii.gz')
+    myGauss =  fsl.SpatialFilter(
+        in_file = input_file,
+        out_file = output_file, 
+        operation = 'median',
+        kernel_shape = 'box',
+        kernel_size = 0.1
+    )
     myGauss.run()
     print('Smoothing DONE!')
     return output_file
 
-def thresh(input_file,outputPath):
+def thresh(input_file, output_path):
     #output_file = os.path.join(os.path.dirname(input_file),os.path.basename(input_file).split('.')[0]+ 'Thres.nii.gz')
-    output_file = os.path.join(outputPath, os.path.basename(input_file).split('.')[0] + 'Thres.nii.gz')
+    output_file = os.path.join(output_path, os.path.basename(input_file).split('.')[0] + 'Thres.nii.gz')
     myThres = fsl.Threshold(in_file=input_file,out_file=output_file,thresh=20)#,direction='above')
     myThres.run()
     print('Thresholding DONE!')
     return output_file
 
-def cropToSmall(input_file,outputPath):
+def cropToSmall(input_file,output_path):
     #output_file = os.path.join(os.path.dirname(input_file),os.path.basename(input_file).split('.')[0]  + 'Crop.nii.gz')
-    output_file = os.path.join(outputPath, os.path.basename(input_file).split('.')[0] + 'Crop.nii.gz')
+    output_file = os.path.join(output_path, os.path.basename(input_file).split('.')[0] + 'Crop.nii.gz')
     myCrop = fsl.ExtractROI(in_file=input_file,roi_file=output_file,x_min=40,x_size=130,y_min=50,y_size=110,z_min=0,z_size=12)
     myCrop.run()
     print('Cropping DONE!')
@@ -113,23 +119,23 @@ def cropToSmall(input_file,outputPath):
     args = parser.parse_args()
 
     # set parameters
-    inputFile = None
+    input_file = None
     if args.input is not None and args.input is not None:
-        inputFile = args.input
+        input_file = args.input
 
-    if not os.path.exists(inputFile):
-        sys.exit("Error: '%s' is not an existing directory or file %s is not in directory." % (inputFile, args.file,))
+    if not os.path.exists(input_file):
+        sys.exit("Error: '%s' is not an existing directory or file %s is not in directory." % (input_file, args.file,))
 
     frac = args.frac
     radius = args.radius
     vertical_gradient = args.vertical_gradient
-    outputPath = os.path.dirname(inputFile)
+    output_path = os.path.dirname(input_file)
 
     # 1) Process DTI
     print("DTI Preprocessing  \33[5m...\33[0m (wait!)", end="\r")
 
     # generate log - file
-    sys.stdout = open(os.path.join(os.path.dirname(inputFile), 'preprocess.log'), 'w')
+    sys.stdout = open(os.path.join(os.path.dirname(input_file), 'preprocess.log'), 'w')
 
     # print parameters
     print("Frac: %s" % frac)
@@ -139,12 +145,12 @@ def cropToSmall(input_file,outputPath):
     # 1) Process MRI
     print('Start Preprocessing ...')
 
-    outputSmooth = smoothIMG(input_file=inputFile,outputPath=outputPath)
+    output_smooth = smoothIMG(input_file = input_file, output_path = output_path)
     # intensity correction using non parametric bias field correction algorithm
-    outputMICO = applyMICO.run_MICO(outputSmooth,outputPath)
+    output_mico = applyMICO.run_MICO(output_smooth, output_path)
 
     # get rid of your skull
-    outputBET = applyBET(input_file=outputMICO,frac=frac,radius=radius,outputPath=outputPath)
+    outputBET = applyBET(input_file = output_mico, frac = frac, radius = radius, output_path = output_path)
 
     sys.stdout = sys.__stdout__
     print('DTI Preprocessing  \033[0;30;42m COMPLETED \33[0m')

bin/3.2_DTIConnectivity/dsi_tools_20170214.py

@@ -79,31 +79,27 @@ def fsl_SeparateSliceMoCo(input_file, par_folder):
 
     # sparate ref and src volume in slices
     sliceFiles = findSlicesData(os.getcwd(), dataName)
-    # refFiles = findSlicesData(os.getcwd(),'ref')
     print('For all slices ... ')
 
     # start to correct motions slice by slice
     for i in range(len(sliceFiles)):
         slc = sliceFiles[i]
-        # ref = refFiles[i]
-        # take epi as ref
         output_file = os.path.join(par_folder, os.path.basename(slc))
         myMCFLIRT = fsl.preprocess.MCFLIRT(in_file=slc, out_file=output_file, save_plots=True, terminal_output='none')
         print(myMCFLIRT.cmdline)
         myMCFLIRT.run()
         os.remove(slc)
-        # os.remove(ref)
 
     # merge slices to a single volume
-
     mcf_sliceFiles = findSlicesData(par_folder, dataName)
     output_file = os.path.join(os.path.dirname(input_file),
                                os.path.basename(input_file).split('.')[0]) + '_mcf.nii.gz'
     myMerge = fsl.Merge(in_files=mcf_sliceFiles, dimension='z', merged_file=output_file)
     print(myMerge.cmdline)
     myMerge.run()
 
-    for slc in mcf_sliceFiles: os.remove(slc)
+    for slc in mcf_sliceFiles: 
+        os.remove(slc)
 
     # unscale result data by factor 10**(-1)
     output_file = scaleBy10(output_file, inv=True)
@@ -142,7 +138,7 @@ def move_files(dir_in, dir_out, pattern):
 
 def connectivity(dsi_studio, dir_in, dir_seeds, dir_out, dir_con):
     """
-    perform seed-based fiber-tracking
+    Calculates connectivity data (types: pass and end).
     """
     if not os.path.exists(dir_in):
         sys.exit("Input directory \"%s\" does not exist." % (dir_in,))
@@ -160,24 +156,24 @@ def connectivity(dsi_studio, dir_in, dir_seeds, dir_out, dir_con):
     os.chdir(os.path.dirname(dir_in))
     cmd_ana = r'%s --action=%s --source=%s --tract=%s --connectivity=%s --connectivity_value=%s --connectivity_type=%s'
 
-
     filename = glob.glob(dir_in+'/*fib.gz')[0]
     file_trk = glob.glob(dir_in+'/*trk.gz')[0]
     file_seeds = dir_seeds
 
-    #> Fixing connectivity calculation being interrupted due to lack of connectivity values // VVF 23/05/09
+    # Performs analysis on every connectivity value within the list ('qa' may not be necessary; might be removed in the future.)
     connect_vals = ['qa', 'count']
     for i in connect_vals:
         parameters = (dsi_studio, 'ana', filename, file_trk, file_seeds, i, 'pass,end')
-        #parameters = (dsi_studio, 'ana', filename, file_trk, file_seeds, 'fa,count', 'pass,end')
         os.system(cmd_ana % parameters)
 
     #move_files(dir_in, dir_con, re.escape(filename) + '\.' + re.escape(pre_seeds) + '.*(?:\.pass\.|\.end\.)')
     move_files(os.path.dirname(file_trk), dir_con, '/*.txt')
     move_files(os.path.dirname(file_trk), dir_con, '/*.mat')
 
-#? Function (mapsgen(args)) seems to be unused and remains at that state. Can it be savely removed? // VVF 23/05/09
 def mapsgen(dsi_studio, dir_in, dir_msk, b_table, pattern_in, pattern_fib):
+    """
+    FUNCTION DEPRECATED. REMOVAL PENDING.
+    """
     pre_msk = 'bet.bin.'
 
     ext_src = '.src.gz'
@@ -206,10 +202,8 @@ def mapsgen(dsi_studio, dir_in, dir_msk, b_table, pattern_in, pattern_fib):
     for index, filename in enumerate(file_list):
         # create source files
         pos = filename.rfind('_')
-        #file_in  = os.path.join(dir_in, filename)
-        #file_src = os.path.join(dir_in, filename[:pos] + ext_src)
+
         file_src = filename[:pos] + ext_src
-        #parameters = (dsi_studio, 'src', file_in, file_src, b_table)
         parameters = (dsi_studio, 'src', filename, file_src, b_table)
         print("%d of %d:" % (index + 1, len(file_list)), cmd_src % parameters)
         subprocess.call(cmd_src % parameters)
@@ -234,19 +228,15 @@ def mapsgen(dsi_studio, dir_in, dir_msk, b_table, pattern_in, pattern_fib):
         subprocess.call(cmd_exp % parameters)
 
 def srcgen(dsi_studio, dir_in, dir_msk, dir_out, b_table):
-    #dir_src = r'..\src'
-    #dir_fib = r'..\fib_map'
+    """
+    Sources and creates fib files. Diffusivity and aniosotropy metrics are exported from data.
+    """
     dir_src = r'src'
     dir_fib = r'fib_map'
-    #dir_map = r'maps'
-    #dir_gfa = r'gfa'
     dir_qa  = r'DSI_studio'
     dir_con = r'connectivity'
-
-
     ext_src = '.src.gz'
 
-
     if not os.path.exists(dir_in):
         sys.exit("Input directory \"%s\" does not exist." % (dir_in,))
 
@@ -261,12 +251,8 @@ def srcgen(dsi_studio, dir_in, dir_msk, dir_out, b_table):
     if not os.path.isfile(b_table):
         sys.exit("File \"%s\" does not exist." % (b_table,))
 
-
-
     dir_src = make_dir(os.path.dirname(dir_out), dir_src)
     dir_fib = make_dir(os.path.dirname(dir_out), dir_fib)
-    #dir_map = make_dir(os.path.dirname(dir_fib), dir_map)
-    #dir_gfa = make_dir(os.path.dirname(dir_map), dir_gfa)
     dir_qa  = make_dir(os.path.dirname(dir_out), dir_qa)
 
     # change to input directory
@@ -278,9 +264,7 @@ def srcgen(dsi_studio, dir_in, dir_msk, dir_out, b_table):
     # create source files
     filename = os.path.basename(dir_in)
     pos = filename.rfind('.')
-    #file_in  = os.path.join(dir_in, filename)
     file_src = os.path.join(dir_src, filename[:pos] + ext_src)
-    #parameters = (dsi_studio, 'src', file_in, file_src, b_table)
     parameters = (dsi_studio, 'src', filename, file_src, b_table)
     print("Generate src-File %s:" % cmd_src % parameters)
     os.system(cmd_src % parameters)
@@ -291,32 +275,31 @@ def srcgen(dsi_studio, dir_in, dir_msk, dir_out, b_table):
     print("Generate fib-File %s:" % cmd_rec % parameters)
     os.system(cmd_rec % parameters)
 
-
     # move fib to corresponding folders
     move_files(dir_src, dir_fib, '/*fib.gz')
 
-    # # extracts maps: 2 ways:
+    # extracts maps: 2 ways:
     cmd_exp = r'%s --action=%s --source=%s --export=%s'
     file_fib = glob.glob(dir_fib+'/*fib.gz')[0]
     parameters = (dsi_studio, 'exp', file_fib, 'fa')
     print("Generate two maps %s:" % cmd_exp % parameters)
     os.system(cmd_exp % parameters)
 
-    # # extracts maps: 2 ways:
+    # extracts maps: 2 ways:
     cmd_exp = r'%s --action=%s --source=%s --export=%s'
     file_fib = glob.glob(dir_fib + '/*fib.gz')[0]
     parameters = (dsi_studio, 'exp', file_fib, 'md')
     print("Generate two maps %s:" % cmd_exp % parameters)
     os.system(cmd_exp % parameters)
 
-    # # extracts maps: 2 ways:
+    # extracts maps: 2 ways:
     cmd_exp = r'%s --action=%s --source=%s --export=%s'
     file_fib = glob.glob(dir_fib + '/*fib.gz')[0]
     parameters = (dsi_studio, 'exp', file_fib, 'ad')
     print("Generate two maps %s:" % cmd_exp % parameters)
     os.system(cmd_exp % parameters)
 
-    # # extracts maps: 2 ways:
+    # extracts maps: 2 ways:
     cmd_exp = r'%s --action=%s --source=%s --export=%s'
     file_fib = glob.glob(dir_fib + '/*fib.gz')[0]
     parameters = (dsi_studio, 'exp', file_fib, 'rd')
@@ -331,7 +314,7 @@ def srcgen(dsi_studio, dir_in, dir_msk, dir_out, b_table):
 
 def tracking(dsi_studio, dir_in):
     """
-    perform seed-based fiber-tracking
+    Performs seed-based fiber-tracking.
     """
     if not os.path.exists(dir_in):
         sys.exit("Input directory \"%s\" does not exist." % (dir_in,))