Merge pull request #324 from soazig/soazig-master

Updating master - fixed the data.type for image data?

Makefile

@@ -53,7 +53,7 @@ multi-comparison:
 	cd code/utils/scripts && python multi_beta_script.py
 	cd code/utils/scripts && python multi_comparison_script.py
 
-all-analysis:
+analysis-except-multi:
 	make eda 
 	make linear
 	make logistic
@@ -63,9 +63,10 @@ all-analysis:
 	make glm
 	make correlation
 	make noise-pca 
-	make multi-comparison
 
 report:
 	cd paper && make all
+	cd paper && make clean
+
 
 

code/utils/scripts/convolution_high_res_script.py

@@ -41,7 +41,7 @@
 data_path = project_path+'data/ds005/'
 
 #change here to get your subject !
-subject_list = ['11', '5', '1']
+subject_list = [str(i) for i in range(1,17)]
 #change here to get your run number !
 run_list = [str(i) for i in range(1,4)]
 cond_list = [str(i) for i in range(1,5)]

code/utils/scripts/convolution_normal_script.py

@@ -38,7 +38,7 @@
 data_path = project_path+'data/ds005/'
 
 #change here to get your subject !
-subject_list = ['1','5', '11']
+subject_list = [str(i) for i in range(1,17)]
 #change here to get your run number !
 run_list = [str(i) for i in range(1,4)]
 cond_list = [str(i) for i in range(1,5)]

code/utils/scripts/eda.py

@@ -36,7 +36,7 @@
 
 # Run only for subject 1 and 5 - run 1
 run_list = [str(i) for i in range(1,2)]
-subject_list = ['1', '5']
+subject_list = ['1','5']
 
 # set gray colormap and nearest neighbor interpolation by default
 plt.rcParams['image.cmap'] = 'gray'
@@ -52,6 +52,12 @@
 # Template to plot the unmasked filetered data
 template_path = project_path+'data/mni_icbm152_t1_tal_nlin_asym_09c_2mm.nii'
 
+# Progress bar
+l = len(subject_list)*len(run_list)
+sys.stdout.write("Starting EDA analysis\n")
+sys.stdout.write("EDA: ")
+sys.stdout.flush()
+
 # Loop through the data type - raw or filtered
 for dat in path_dict:
     d_path = path_dict[dat]
@@ -64,23 +70,18 @@
                      for s in subject_list]
     for image_path in images_paths:
         name = image_path[0]
+        data_int = nib.load(image_path[1]).get_data()
+        data = data_int.astype(float)
+        mean_data = np.mean(data, axis=-1)
         # Plot
         if d_path['type']=='filtered':
-	    data = nib.load(image_path[1]).get_data()
-	    data = data.astype(float)
-            mean_data = np.mean(data, axis=-1)
             Transpose=False
 	    template_data = nib.load(template_path).get_data()
 	    plt.imshow(\
 	        plot_mosaic(template_data, transpose=Transpose), \
 	        cmap='gray', alpha=1)
 	else:
-            img = nib.load(image_path[1])
-            data = img.get_data()
-	    data = data.astype(float)
-            mean_data = np.mean(data, axis=-1)
 	    in_brain_mask = mean_data > 375
-            mean_data = np.mean(data, axis=-1)
             Transpose=True
             plt.contour(\
 	        plot_mosaic(in_brain_mask, transpose=Transpose), \
@@ -93,9 +94,10 @@
                     %(d_path['type'] + str(name)))
         #plt.show()
         plt.clf()
+        plt.close()
+        sys.stdout.write("\n\b=")
+        sys.stdout.flush()
 
-
-print("======================================")
-print("\nEDA analysis done")
-print("Mosaic plots in project_epsilon/fig/BOLD/ \n\n")
-
+sys.stdout.write("======================================\n")
+sys.stdout.write("EDA analysis done\n")
+sys.stdout.write("Mosaic plots in project_epsilon/fig/BOLD/ \n")

code/utils/scripts/hist-outliers_script.py

@@ -47,7 +47,7 @@
 dir_path = project_path+'data/ds005/'
 # Uncomment to run all subjects
 #subject_list = [str(i) for i in range(1,17)]
-subject_list = ['1','5']
+subject_list = ['1', '5']
 
 # Choose run here (1,2 and/or 3)
 run_list = [str(i) for i in range(1,2)]
@@ -61,7 +61,8 @@
 for image_path in images_paths:
     name = image_path[0]
     img = nib.load(image_path[1])
-    data = img.get_data()
+    data_int = img.get_data()
+    data = data_int.astype(float)
     vol_shape = data.shape[-1]
     mean_vol = np.mean(data, axis=-1)
 #    Take the mean value over time and plot an histogram

code/utils/scripts/mask_filtered_data.py

@@ -1,4 +1,4 @@
-"""mask_functions.py
+"""mask_filtered_data.py
 
 A collection of functions to make masks on data.
 See test_* functions in this directory for nose tests
@@ -28,7 +28,8 @@ def make_bool_mask(mask_path):
         the bool array created from the binary array
     """
     mask_img = nib.load(mask_path)
-    mask_data = mask_img.get_data()
+    mask_data_int = mask_img.get_data()
+    mask_data = mask_img_int.astype(float)
     mask_bool = mask_data == 1
     return mask_bool
 
@@ -52,8 +53,10 @@ def make_mask_filtered_data(func_path, mask_path):
     """
     func_img = nib.load(func_path)
     mask_img = nib.load(mask_path)
-    mask = mask_img.get_data()
-    func_data = func_img.get_data()
+    mask_int = mask_img.get_data()
+    func_data_int = func_img.get_data()
+    mask = mask_int.astype(float)
+    func_data = func_data_int.astype(float)
     # Make data 4D to prepare for "broadcasting"
     mask = np.reshape(mask, mask.shape + (1,))
     # "Broadcasting" expands the final length 1 dimension to match the func data

code/utils/scripts/noise-pca_filtered_script.py

@@ -2,7 +2,10 @@
 This script is used to design the design matrix for our linear regression.
 We explore the influence of linear and quadratic drifts on the model 
 performance.
-
+Script for the filtered data.
+Run with:
+    python noise-pca_filtered_script.py
+    from this directory
 """
 from __future__ import print_function, division
 import sys, os, pdb
@@ -34,7 +37,7 @@
 data_path = project_path+'data/ds005/' 
 path_dict = {'data_filtered':{ 
 			      'type' : 'filtered',
-			      'feat' : '.feat/',
+			      'feat' : '.feat',
 			      'bold_img_name' : 'filtered_func_data_mni.nii.gz',
 			      'run_path' : 'model/model001/'
 			     },
@@ -51,7 +54,7 @@
 
 # Run only for subject 1 and 5 - run 1
 run_list = [str(i) for i in range(1,2)]
-subject_list = ['1','5']
+subject_list = ['1', '5']
 
 d_path = path_dict['data_filtered'] #OR original or filtered 
 
@@ -70,6 +73,7 @@
 thres = 375 #From analysis of the histograms   
 # To be used with the filtered data
 mask_path = project_path+'data/mni_icbm152_t1_tal_nlin_asym_09c_mask_2mm.nii'
+template_path = project_path+'data/mni_icbm152_t1_tal_nlin_asym_09c_2mm.nii'
 
 sm = ''
 #sm='not_smooth/'
@@ -92,24 +96,47 @@
     if not os.path.exists(d):
         os.makedirs(d)
 
+# Progress bar
+print("\nStarting noise-pca for filtered data analysis\n")
+
 for image_path in images_paths:
     name = image_path[0]
     if d_path['type']=='filtered':
-        in_brain_img = nib.load('../../../'+
-	    'data/ds005/sub001/model/model001/task001_run001.feat/'\
-	    + 'masked_filtered_func_data_mni.nii.gz')
+        #in_brain_img = nib.load('../../../'+
+	#    'data/ds005/sub001/model/model001/task001_run001.feat/'\
+	#    + 'masked_filtered_func_data_mni.nii.gz')
 	# Image shape (91, 109, 91, 240)
-	#in_brain_img = make_mask_filtered_data(image_path[1],mask_path)
-	data = in_brain_img.get_data()
+        md = data_path + 'sub%s/'%(s.zfill(3)) + d_path['run_path'] \
+        + 'task001_run%s%s/masked_%s' %(r.zfill(3),d_path['feat'],\
+        d_path['bold_img_name']) 
+        if not os.path.exists(md):
+            print("Filtering brain image for: ")
+	    print(str(name))
+            in_brain_img = make_mask_filtered_data(image_path[1],mask_path)
+	    print("brain image filtered\n")
+        else:
+	    print("Loading filtered brain image for: ")
+	    print(str(name))
+            in_brain_img = nib.load(md)
+            print("brain image loaded\n")
+        data_int = in_brain_img.get_data()
+        data = data_int.astype(float)
         mean_data = np.mean(data, axis=-1)
-        in_brain_mask = (mean_data - 0.0) < 0.01	
+        template = nib.load(template_path)
+	template_data_int = template.get_data()
+	template_data = template_data_int.astype(float)
 	Transpose = False
+        in_brain_mask = (mean_data - 0.0) < 0.01
+        plt.imshow(plot_mosaic(template_data, transpose=Transpose),\
+				        cmap='gray', alpha=1)
     else:
         img = nib.load(image_path[1])
         data = img.get_data()
         mean_data = np.mean(data, axis=-1)
         in_brain_mask = mean_data > thres
 	Transpose = True
+        plt.contour(plot_mosaic(in_brain_mask, transpose=Transpose), \
+	            cmap='gray' , alpha=1)
     # Smoothing with Gaussian filter
     smooth_data = smoothing(data,1,range(data.shape[-1])) 
 
@@ -120,9 +147,8 @@
     # Plotting the voxels in the brain
     plt.imshow(plot_mosaic(mean_data, transpose=Transpose), cmap='gray', alpha=1)
     plt.colorbar()
-    plt.contour(plot_mosaic(in_brain_mask, transpose=Transpose),colors='blue')
     plt.title('In brain voxel mean values' + '\n' +  (d_path['type'] + str(name)))
-    plt.savefig(project_path+'fig/BOLD/%s_mean_voxels_countour.png'\
+    plt.savefig(project_path+'fig/BOLD/%s_mean_voxels.png'\
                 %(d_path['type'] + str(name)))
     #plt.show()
     #plt.clf()
@@ -347,13 +373,13 @@
 #		%(name, d_path['type'] + str(name),str(reg_str[k])))
 #        pdb.set_trace() 
 #    pdb.set_trace() 
-
+    plt.close()
+    print("=")
 
 print("======================================")
-print("\n Noise and PCA analysis done")
-print("Design Matrix including drift terms \
-      stored in project_epsilon/txt_output/drifts/ \n\n")
-print("Design Matrix including PCs terms \
-      stored in project_epsilon/txt_output/pca/\n\n")
+print("\n Noise and PCA analysis for filtered data done")
+print("Design Matrix including drift terms stored in project_epsilon/txt_output/drifts/ \n\n")
+print("Design Matrix including PCs terms stored in project_epsilon/txt_output/pca/\n\n")
+print("Mean MRSS models results in  project_epsilon/txt_output/MRSS/ds005filtered_MRSS.json\n\n")
 
 

code/utils/scripts/noise-pca_script.py

@@ -2,8 +2,12 @@
 This script is used to design the design matrix for our linear regression.
 We explore the influence of linear and quadratic drifts on the model 
 performance.
-
+Script for the raw data.
+Run with:
+    python noise-pca_script.py
+    from this directory
 """
+
 from __future__ import print_function, division
 import sys, os, pdb
 from scipy import ndimage
@@ -51,7 +55,7 @@
 
 # Run only for subject 1 and 5 - run 1
 run_list = [str(i) for i in range(1,2)]
-subject_list = ['1','5']
+subject_list = ['1', '5']
 
 d_path = path_dict['data_original'] #OR original or filtered 
 
@@ -91,6 +95,7 @@
 for d in dirs:
     if not os.path.exists(d):
         os.makedirs(d)
+print("Starting noise-pca for the raw data analysis\n")
 
 for image_path in images_paths:
     name = image_path[0]
@@ -100,13 +105,15 @@
 	    + 'masked_filtered_func_data_mni.nii.gz')
 	# Image shape (91, 109, 91, 240)
 	#in_brain_img = make_mask_filtered_data(image_path[1],mask_path)
-	data = in_brain_img.get_data()
+	data_int = in_brain_img.get_data()
+        data = data_int.astype(float)
         mean_data = np.mean(data, axis=-1)
         in_brain_mask = (mean_data - 0.0) < 0.01	
 	Transpose = False
     else:
         img = nib.load(image_path[1])
-        data = img.get_data()
+        data_int = img.get_data()
+        data = data_int.astype(float)
         mean_data = np.mean(data, axis=-1)
         in_brain_mask = mean_data > thres
 	Transpose = True
@@ -350,9 +357,8 @@
 
 print("======================================")
 print("\n Noise and PCA analysis done")
-print("Design Matrix including drift terms \
-      stored in project_epsilon/txt_output/drifts/ \n\n")
-print("Design Matrix including PCs terms \
-      stored in project_epsilon/txt_output/pca/\n\n")
+print("Design Matrix including drift terms stored in project_epsilon/txt_output/drifts/ \n\n")
+print("Design Matrix including PCs terms stored in project_epsilon/txt_output/pca/\n\n")
+print("Mean MRSS models results in project_epsilon/txt_output/MRSS/ds005_MRSS.json\n\n")
 
 

code/utils/scripts/plot_mosaic.py

@@ -48,11 +48,13 @@ def plot_mosaic(img_data, transpose=False):
     #'../../../data/ds005/sub001/BOLD/task001_run001/bold.nii.gz') 
     template = nib.load(project_path+\
                'data/mni_icbm152_t1_tal_nlin_asym_09c_2mm.nii')
-    template_data = template.get_data()
+    template_data_int = template.get_data()
+    template_data = template_data_int.astype(float)
     img = nib.load(project_path+\
          'data/ds005/sub001/model/model001/task001_run001.feat/' + \
 	 'masked_filtered_func_data_mni.nii.gz')
-    img_data = img.get_data()
+    img_data_int = img.get_data()
+    img_data = img_data_int.astype(float) 
     mean_data = np.mean(img_data, axis=-1)
     plt.title('In brain voxels - mean values')
     plt.imshow(plot_mosaic(template_data, transpose=False), cmap='gray', alpha=1)