refactor code

create-dataset/1-simulation-create-noise-free-data.py

@@ -29,22 +29,17 @@ def __init__(self, bvals, bvecs):
         self.bvals = bvals
         self.bvecs = bvecs
 
-    def add_noise(Sj, SNR, type_noise):
-        # b0vols = np.argwhere(bvals>50) # Get the b0 volumes, where there is no diffusion (pure T2 contrast)
-        # mean_b0vols = np.mean(Sj[b0vols])
-        # if mean_b0vols<0.000001: # Background voxels
-        #     mean_b0vols = 1
-        # sigma = mean_b0vols / SNR # SNR = mean(signal)/std(noise) --> mean(b0_vols)/sigma 
+    def add_noise(Sj, SNR, type_noise):        
         sigma = 1 / SNR 
 
         if type_noise == 'Gaussian':
             try:
                 random = np.random.normal(0, sigma, len(Sj))                
-                Sj_noise = Sj + random #sig * np.random.normal(0, s0/SNR, len(Sj))
+                Sj_noise = Sj + random
             except:                
                 print("Exception sigma")                
         elif type_noise == 'Rician':  # Noise in quadrature
-            noise_1 = np.random.normal(0, sigma, len(Sj))  #sigma or sigma**2?
+            noise_1 = np.random.normal(0, sigma, len(Sj)) 
             noise_2 = np.random.normal(0, sigma, len(Sj))
             Sj_noise = np.sqrt((Sj + noise_1) ** 2 + noise_2 ** 2)
 
@@ -63,8 +58,7 @@ def __call__(self, params):
             fi = params[1 + 3 * i]
             sumf += fi
             th = params[2 + 3 * i]
-            phi = params[3 + 3 * i]
-            #v = np.array([params[3 + 3 * i], params[4 + 3 * i], params[5 + 3 * i]])
+            phi = params[3 + 3 * i]            
             v = np.array([np.sin(th) * np.cos(phi), np.sin(th) * np.sin(phi), np.cos(th)])  # conversion to cartesians
             signal += s0 * (fi * np.exp(-d * self.bvals * np.power(np.dot(self.bvecs.T, v), 2)))    # sticks contribution to the signal
 
@@ -118,7 +112,7 @@ def cart2sph(x,y,z):
   return r, theta, phi
 
 #Load files
-dPath = 'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/data'
+dPath = '/Simulations/data'
 data = get_data(f'{dPath}/data.nii.gz')
 mask = get_data(f'{dPath}/nodif_brain_mask.nii.gz')
 mean_d = get_data(f'{dPath}/data.bedpostX/mean_dsamples.nii.gz')

create-dataset/2-simulation-add-noise-to-free-data.py

@@ -29,23 +29,17 @@ def __init__(self, bvals, bvecs):
         self.bvals = bvals
         self.bvecs = bvecs
 
-    def add_noise(Sj, SNR, type_noise):
-        # b0vols = np.argwhere(bvals>50) # Get the b0 volumes, where there is no diffusion (pure T2 contrast)
-        # mean_b0vols = np.mean(Sj[b0vols])
-        # if mean_b0vols<0.000001: # Background voxels
-        #     mean_b0vols = 1
-        # sigma = mean_b0vols / SNR # SNR = mean(signal)/std(noise) --> mean(b0_vols)/sigma 
+    def add_noise(Sj, SNR, type_noise):        
         sigma = 1 / SNR 
 
         if type_noise == 'Gaussian':
             try:
                 random = np.random.normal(0, sigma, len(Sj))                
-                Sj_noise = Sj + random #sig * np.random.normal(0, s0/SNR, len(Sj))
-                #Sj_noise_modificado = np.where(Sj_noise < 0, 0.0000001, Sj_noise)
+                Sj_noise = Sj + random
             except Exception as e:                
                 print("Se produjo una excepción:", e)             
         elif type_noise == 'Rician':  # Noise in quadrature
-            noise_1 = np.random.normal(0, sigma, len(Sj))  #sigma or sigma**2?
+            noise_1 = np.random.normal(0, sigma, len(Sj))
             noise_2 = np.random.normal(0, sigma, len(Sj))
             Sj_noise = np.sqrt((Sj + noise_1) ** 2 + noise_2 ** 2)
 
@@ -64,8 +58,7 @@ def __call__(self, params):
             fi = params[1 + 3 * i]
             sumf += fi
             th = params[2 + 3 * i]
-            phi = params[3 + 3 * i]
-            #v = np.array([params[3 + 3 * i], params[4 + 3 * i], params[5 + 3 * i]])
+            phi = params[3 + 3 * i]            
             v = np.array([np.sin(th) * np.cos(phi), np.sin(th) * np.sin(phi), np.cos(th)])  # conversion to cartesians
             signal += s0 * (fi * np.exp(-d * self.bvals * np.power(np.dot(self.bvecs.T, v), 2)))    # sticks contribution to the signal
 
@@ -119,7 +112,7 @@ def cart2sph(x,y,z):
   return r, theta, phi
 
 #Load files
-dPath = 'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/noise-free-data'
+dPath = '/Simulations/noise-free-data'
 data = get_data(f'{dPath}/data_b1k.nii.gz')
 orig_data = nb.load(f'{dPath}/data_b1k.nii.gz', mmap=True)	# This is to capture the correct header of the nifti image when exporting
 

create-dataset/3-simulation-normalize-data.py

@@ -45,11 +45,11 @@ def export_nifti(data, orig_data, output_path, name):
 
 
 #Load files
-dPath = 'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/noise-free-data'
+dPath = '/Simulations/noise-free-data'
 data = get_data(f'{dPath}/test-noisy_data_snr10_b0b1b2k.nii.gz')
 
 
-#los b0 son las primeras 5 dimensiones de mi array
+#b0s are the first 5 elements
 ax_signal = len(data.shape) - 1
 mean_nob0vols = np.mean(data[..., :5], axis=ax_signal)
 

create-dataset/4-simulation-generate-10-datasets-from-noisefree-snr.py

@@ -29,23 +29,18 @@ def __init__(self, bvals, bvecs):
         self.bvals = bvals
         self.bvecs = bvecs
 
-    def add_noise(Sj, SNR, type_noise):
-        # b0vols = np.argwhere(bvals>50) # Get the b0 volumes, where there is no diffusion (pure T2 contrast)
-        # mean_b0vols = np.mean(Sj[b0vols])
-        # if mean_b0vols<0.000001: # Background voxels
-        #     mean_b0vols = 1
-        # sigma = mean_b0vols / SNR # SNR = mean(signal)/std(noise) --> mean(b0_vols)/sigma 
+    def add_noise(Sj, SNR, type_noise):        
         sigma = 1 / SNR 
 
         if type_noise == 'Gaussian':
             try:
                 random = np.random.normal(0, sigma, len(Sj))                
-                Sj_noise = Sj + random #sig * np.random.normal(0, s0/SNR, len(Sj))
+                Sj_noise = Sj + random
                 Sj_noise_modificado = np.where(Sj_noise < 0, 0.0000001, Sj_noise)
             except Exception as e:                
                 print("Se produjo una excepción:", e)             
         elif type_noise == 'Rician':  # Noise in quadrature
-            noise_1 = np.random.normal(0, sigma, len(Sj))  #sigma or sigma**2?
+            noise_1 = np.random.normal(0, sigma, len(Sj))
             noise_2 = np.random.normal(0, sigma, len(Sj))
             Sj_noise = np.sqrt((Sj + noise_1) ** 2 + noise_2 ** 2)
 
@@ -58,14 +53,13 @@ def __call__(self, params):
         d = params[0]
         v = np.zeros((n_fib, 3))
         sumf = 0
-        signal = torch.tensor((np.zeros((len(self.bvals)))))  # np.zeros((len(b)))
+        signal = torch.tensor((np.zeros((len(self.bvals)))))
 
         for i in range(0, n_fib):
             fi = params[1 + 3 * i]
             sumf += fi
             th = params[2 + 3 * i]
-            phi = params[3 + 3 * i]
-            #v = np.array([params[3 + 3 * i], params[4 + 3 * i], params[5 + 3 * i]])
+            phi = params[3 + 3 * i]            
             v = np.array([np.sin(th) * np.cos(phi), np.sin(th) * np.sin(phi), np.cos(th)])  # conversion to cartesians
             signal += s0 * (fi * np.exp(-d * self.bvals * np.power(np.dot(self.bvecs.T, v), 2)))    # sticks contribution to the signal
 
@@ -119,10 +113,10 @@ def cart2sph(x,y,z):
   return r, theta, phi
 
 #Load files
-dPath = 'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/data/generated_data_noisyfree_b0_first'
+dPath = '/Simulations/data/generated_data_noisyfree_b0_first'
 data = get_data(f'{dPath}/noisyfree_data_full_b0_first.nii.gz')
 orig_data = nb.load(f'{dPath}/noisyfree_data_full_b0_first.nii.gz', mmap=True)	# This is to capture the correct header of the nifti image when exporting
-dPath_toSave = 'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/data/generated_data_noisyfree_b0_first/snrs'
+dPath_toSave = '/Simulations/data/generated_data_noisyfree_b0_first/snrs'
 
 SNRs = [40]
 
@@ -145,7 +139,7 @@ def save_noisy_dataset(noisy_dataset, snr, i):
     dataset_name = f'{i}noisy_data_snr{snr}.nii.gz'
     export_nifti(noisy_dataset, orig_data, dPath_toSave, dataset_name)
 
-for i in range(1, 6): 
+for i in range(1, 10): 
     print('Adding noise to dataset')
     noisy_data = add_noise_to_data(data, SNRs)
     print('Saving noising dataset')

metrics/2-Generate-image-simulated-distributions.py → metrics/Generate_Fig1-simulated-distributions.py

@@ -57,16 +57,16 @@ def plot_gaussian_rician_raw_noise():
     plt.legend(SNRs_leyend)
     fig.suptitle('Simulated distributions', fontsize=14)   
 
-    plt.savefig(f'{dPath_gaussian}/figures/{iter}_histogram_raw_rician_gaussian_noise_screenshots.png', dpi=600, bbox_inches='tight')
+    plt.savefig(f'{dPath_gaussian}/figures/{iter}_Fig1_histogram_raw_rician_gaussian_noise_screenshots.png', dpi=600, bbox_inches='tight')
     plt.show()
 
 
 
 SNRs = [3,5,10,20,40]
 iter = 20
 
-dPath_gaussian = f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/Exp6-data-gaussian'
-dPath_rician = f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/Exp6-data-rician'
+dPath_gaussian = f'/Simulations/Experiments/Exp6-data-gaussian'
+dPath_rician = f'/Simulations/Experiments/Exp6-data-rician'
 
 
 plot_gaussian_rician_raw_noise()

metrics/5-Generate-image-noise-distribution-ventricles.py → metrics/Generate_Fig10-11_noise-distribution-ventricles.py

@@ -18,10 +18,10 @@
 custom_palette = sns.color_palette(list_colors) # To get the same color palette than in seaborn
 
 iter = 20
-exp = 'Exp6-data-rician'
-noise_type = 'Rician'
+exp = 'Exp6-data-gaussian'
+noise_type = 'gaussian'
 
-dPath = f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}'
+dPath = f'/Simulations/Experiments/{exp}'
 
 nrows = len(methods)
 ncols = 1
@@ -32,7 +32,7 @@
 
 for i in range(0, nrows):    
     ts_csf = np.loadtxt(f'{dPath}/{methods[i]}/snr{snr}/analysis/dMRI/processed_data/ts_CSF.txt')[3:].ravel() #First 3 columns of the fslmeants output are not of interest 
-    ts_csf = ts_csf[ts_csf>0.001] 
+    ts_csf = ts_csf #[ts_csf>0.001] 
     weights_csf = np.ones_like(ts_csf) / float(len(ts_csf))
     # Histogram
     ax[i].hist(ts_csf, color=custom_palette[i], bins=100, density=True, alpha=0.5)
@@ -51,4 +51,4 @@
 plt.tight_layout()  # Ajustar el espacio entre las subtramas
 plt.show()
 
-fig.savefig(f'{dPath}/figures/{noise_type}_{iter}_fig5a_snr{snr}_noise_dist.png', dpi=600, bbox_inches='tight')
+fig.savefig(f'{dPath}/figures/{noise_type}_{iter}_fig10_snr{snr}_noisefloor_dist.png', dpi=600, bbox_inches='tight')