diff --git a/create-dataset/1-simulation-create-noise-free-data.py b/create-dataset/1-simulation-create-noise-free-data.py index c387bfd..735d150 100644 --- a/create-dataset/1-simulation-create-noise-free-data.py +++ b/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') diff --git a/create-dataset/2-simulation-add-noise-to-free-data.py b/create-dataset/2-simulation-add-noise-to-free-data.py index 04b0ef0..260e125 100644 --- a/create-dataset/2-simulation-add-noise-to-free-data.py +++ b/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 diff --git a/create-dataset/3-simulation-normalize-data.py b/create-dataset/3-simulation-normalize-data.py index 673cfd9..fcf1b0a 100644 --- a/create-dataset/3-simulation-normalize-data.py +++ b/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) diff --git a/create-dataset/4-simulation-generate-10-datasets-from-noisefree-snr.py b/create-dataset/4-simulation-generate-10-datasets-from-noisefree-snr.py index fb23c37..0a31259 100644 --- a/create-dataset/4-simulation-generate-10-datasets-from-noisefree-snr.py +++ b/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') diff --git a/metrics/9-Generate-image-fslcc.py b/metrics/9-Generate-image-fslcc.py deleted file mode 100644 index 5c11751..0000000 --- a/metrics/9-Generate-image-fslcc.py +++ /dev/null @@ -1,61 +0,0 @@ -import numpy as np -import pandas as pd -import matplotlib.pyplot as plt -import nibabel as nb -from scipy import ndimage, stats -from scipy.stats import ranksums -import os -import subprocess -from pathlib import Path -import warnings -import sys -import seaborn as sns - -SNRs = [3,5,10,20,40] -methods = ['Raw','NLMEANS','MPPCA','Patch2Self', 'DDM2', 'AVG'] # -methods_names_short = ['Noisy','NLM', 'MPPCA', 'P2S','DDM2', 'AVG'] -list_colors = ['gold','pink', 'blue', 'darkorange', 'lightseagreen', 'green'] -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' -metric = 'MD' - -dPath = f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}' - -data = {} - -for method in methods: - print(f'method:{method}') - data[method] = [] - for snr in SNRs: - print(f'snr:{snr}') - fsl_cc = np.loadtxt(f'{dPath}/{method}/snr{snr}/analysis/data.dti/fslcc{metric}_res.txt')[2:].ravel() #First 2 columns of the fslcc output are not of interest - fsl_val = 0 - if(fsl_cc.size > 0): - fsl_val = fsl_cc[0] - print(f'fsl_cc: {fsl_val}') - data[method].append(fsl_val) - -print(data) - -# Crear la figura y los ejes -fig, ax = plt.subplots(figsize=(5,5)) - -method_index = 0 -for method in methods: - ax.plot(SNRs, data[method], marker='o', label=methods_names_short[method_index], color=list_colors[method_index]) - method_index = method_index + 1 - -ax.set_xlabel('SNR') -ax.set_ylabel('Correlation') -plt.title(metric) -ax.legend() - -plt.tight_layout() -plt.savefig(f'{dPath}/figures/{noise_type}{iter}_fslcc{metric}_screenshots.png') -# Mostrar la gráfica -plt.show() - - diff --git a/metrics/2-Generate-image-simulated-distributions.py b/metrics/Generate_Fig1-simulated-distributions.py similarity index 83% rename from metrics/2-Generate-image-simulated-distributions.py rename to metrics/Generate_Fig1-simulated-distributions.py index ff62681..46ec9fa 100644 --- a/metrics/2-Generate-image-simulated-distributions.py +++ b/metrics/Generate_Fig1-simulated-distributions.py @@ -57,7 +57,7 @@ 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() @@ -65,8 +65,8 @@ def plot_gaussian_rician_raw_noise(): 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() diff --git a/metrics/5-Generate-image-noise-distribution-ventricles.py b/metrics/Generate_Fig10-11_noise-distribution-ventricles.py similarity index 85% rename from metrics/5-Generate-image-noise-distribution-ventricles.py rename to metrics/Generate_Fig10-11_noise-distribution-ventricles.py index b0e3e08..759fb53 100644 --- a/metrics/5-Generate-image-noise-distribution-ventricles.py +++ b/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') \ No newline at end of file +fig.savefig(f'{dPath}/figures/{noise_type}_{iter}_fig10_snr{snr}_noisefloor_dist.png', dpi=600, bbox_inches='tight') \ No newline at end of file diff --git a/metrics/Generate_Fig12-noise_floor_by_bvalues-NO_DDM2_P2S_snr3.py b/metrics/Generate_Fig12-noise_floor_by_bvalues-NO_DDM2_P2S_snr3.py new file mode 100644 index 0000000..2de587f --- /dev/null +++ b/metrics/Generate_Fig12-noise_floor_by_bvalues-NO_DDM2_P2S_snr3.py @@ -0,0 +1,149 @@ +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import nibabel as nb +from scipy import ndimage, stats +from scipy.stats import ranksums +import os +import subprocess +from pathlib import Path +import warnings +import sys +import seaborn as sns +import math +from skimage.metrics import structural_similarity as ssim +from skimage import feature, transform +from skimage.metrics import peak_signal_noise_ratio as psnr + +def get_data(file, mmap=True): + """ + Load NIfTI image data from a file. + + Parameters: + file (str): The path to the NIfTI file. + mmap (bool, optional): Whether to use memory-mapped file access. Default is True. + + Returns: + numpy.ndarray: The voxel data from the NIfTI file. + """ + import nibabel as nb + img = nb.load(file, mmap=mmap) + img_voxels = img.get_fdata() + return img_voxels + +def calculate_median_noise(denoised_data, bvals, mask): + #calcular la media por cada b-shell + indices_b1 = np.where((bvals > 100) & (bvals < 1500))[0] + indices_b2 = np.where(bvals > 1500)[0] + + masked_data_list_b1 = [denoised_data[:, :, :, i][mask == 1] for i in indices_b1] + combined_masked_data_b1 = np.concatenate(masked_data_list_b1) + median_b1 = np.median(combined_masked_data_b1) + + masked_data_list_b2 = [denoised_data[:, :, :, i][mask == 1] for i in indices_b2] + combined_masked_data_b2 = np.concatenate(masked_data_list_b2) + median_b2 = np.median(combined_masked_data_b2) + + return median_b1, median_b2 +##### + +SNRs = [3,5,10,20,40] #3,5,10,20,40 +methods = ['Raw','NLMEANS','MPPCA','Patch2Self', 'DDM2', 'AVG'] +methods_plot_names = ['Noisy','NLMEANS','MPPCA','Patch2Self', 'DDM2', 'AVG'] + +iter = 21 +exp = 'Exp6-data-gaussian' +noise_type = 'Gaussian' +list_colors = ['gold','pink', 'blue', 'darkorange', 'lightseagreen', 'green'] + + +mask = get_data(f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}/Dataset/mean_f1samples_mask_new.nii.gz') +ground_truth = get_data(f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}/Dataset/noisyfree_data_full_b0_first.nii.gz') +bvals = np.genfromtxt(f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}/Dataset/data_b0b1b2k.bval', dtype=np.float32) + + +# Identificar las posiciones donde la máscara tiene un valor de 1 (dentro de la máscara) +mask_positions = np.where(mask == 1) +dPath = f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}' + +data1k = {} +data2k = {} +for method in methods: + data1k[method] = [] + data2k[method] = [] + +for method in methods: + print(f'method:{method}') + for snr in SNRs: + print(f'snr:{snr}') + if (method == 'DDM2' or method == 'Patch2Self') and snr == 3: + continue # Omitir los métodos DDM y P2S con snr=3 + name_data = f'{method.lower()}-denoised_main_snr' + if method == 'Raw': + name_data = f'main-{noise_type}-noisy_data_snr' + + denoised_data = get_data(f'{dPath}/{method}/snr{snr}/{name_data}{snr}.nii.gz') + + median_b1, median_b2 = calculate_median_noise(denoised_data, bvals, mask) + + print(f'noisefloorb1: {median_b1}') + print(f'noisefloorb2: {median_b2}') + data1k[method].append((snr, median_b1)) + data2k[method].append((snr, median_b2)) + +print(data1k) +print(data2k) + +# Crear la figura y los ejes +fig, ax = plt.subplots(1, 2) + +# Graficar los errores de cada método como una línea +method_index = 0 +linestyle = ['solid', 'dashdot', 'solid', '--', 'solid'] +for method in methods: + print(method) + marker = 'o' + linestyle = 'solid' + #if method == 'NLMEANS': + # marker = 'X' + if method == 'Patch2Self': + linestyle = 'dashdot' + + # Filtrar los valores para no incluir snr=3 para los métodos DDM y P2S + filtered_snr_data1k = [(snr, value) for snr, value in data1k[method] if not ((method == 'DDM2' or method == 'Patch2Self') and snr == 3)] + filtered_snr_data2k = [(snr, value) for snr, value in data2k[method] if not ((method == 'DDM2' or method == 'Patch2Self') and snr == 3)] + + # Separar SNRs y valores + filtered_snr1k = [snr for snr, _ in filtered_snr_data1k] + filtered_values1k = [value for _, value in filtered_snr_data1k] + filtered_snr2k = [snr for snr, _ in filtered_snr_data2k] + filtered_values2k = [value for _, value in filtered_snr_data2k] + + # Graficar los datos filtrados + ax[0].plot(filtered_snr1k, filtered_values1k, marker=marker, label=methods_plot_names[method_index], color=list_colors[method_index], linestyle=linestyle) + ax[1].plot(filtered_snr2k, filtered_values2k, marker=marker, label=methods_plot_names[method_index], color=list_colors[method_index], linestyle=linestyle) + + method_index += 1 + + + + + +# Agregar etiquetas, título y leyenda +ax[0].set_xlabel('SNR') +ax[0].set_ylabel('noise-floor') +ax[0].set_title(r'$b=1000 \, \mathrm{mm/s^2}$') +#ax[0].legend() + +ax[1].set_xlabel('SNR') +ax[1].set_ylabel('') +ax[1].set_title(r'$b=2000 \, \mathrm{mm/s^2}$') +ax[1].legend() + +fig.suptitle('Noise-Floor', fontsize=14) + +plt.tight_layout() +plt.savefig(f'{dPath}/figures/{noise_type}_{iter}_group_noise_floor_screenshots.png') +# Mostrar la gráfica +plt.show() + diff --git a/metrics/3.1-Generate-image-noise-floor-bshell-shells-by-images.py b/metrics/Generate_Fig12-noise_floor_by_bvalues.py similarity index 85% rename from metrics/3.1-Generate-image-noise-floor-bshell-shells-by-images.py rename to metrics/Generate_Fig12-noise_floor_by_bvalues.py index 8478a64..fbacb56 100644 --- a/metrics/3.1-Generate-image-noise-floor-bshell-shells-by-images.py +++ b/metrics/Generate_Fig12-noise_floor_by_bvalues.py @@ -47,7 +47,7 @@ def calculate_median_noise(denoised_data, bvals, mask): return median_b1, median_b2 ##### -SNRs = [3,5,10,20,40] #3,5,10,20,40 +SNRs = [3,5,10,20,40] methods = ['Raw','NLMEANS','MPPCA','Patch2Self', 'DDM2', 'AVG'] methods_plot_names = ['Noisy','NLMEANS','MPPCA','Patch2Self', 'DDM2', 'AVG'] @@ -57,14 +57,14 @@ def calculate_median_noise(denoised_data, bvals, mask): list_colors = ['gold','deeppink', 'blue', 'darkorange', 'lightseagreen', 'green'] -mask = get_data(f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}/Dataset/mean_f1samples_mask_new.nii.gz') -ground_truth = get_data(f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}/Dataset/noisyfree_data_full_b0_first.nii.gz') -bvals = np.genfromtxt(f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}/Dataset/data_b0b1b2k.bval', dtype=np.float32) +mask = get_data(f'/Simulations/Experiments/{exp}/Dataset/mean_f1samples_mask_new.nii.gz') +ground_truth = get_data(f'/Simulations/Experiments/{exp}/Dataset/noisyfree_data_full_b0_first.nii.gz') +bvals = np.genfromtxt(f'/Simulations/Experiments/{exp}/Dataset/data_b0b1b2k.bval', dtype=np.float32) # Identificar las posiciones donde la máscara tiene un valor de 1 (dentro de la máscara) mask_positions = np.where(mask == 1) -dPath = f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}' +dPath = f'/Simulations/Experiments/{exp}' data1k = {} data2k = {} @@ -84,6 +84,7 @@ def calculate_median_noise(denoised_data, bvals, mask): denoised_data = get_data(f'{dPath}/{method}/snr{snr}/{name_data}{snr}.nii.gz') median_b1, median_b2 = calculate_median_noise(denoised_data, bvals, mask) + print(f'noisefloorb1: {median_b1}') print(f'noisefloorb2: {median_b2}') data1k[method].append(median_b1) @@ -136,7 +137,7 @@ def calculate_median_noise(denoised_data, bvals, mask): fig.suptitle('Noise-Floor', fontsize=14) plt.tight_layout() -plt.savefig(f'{dPath}/figures/{noise_type}_{iter}_group_noise_floor_screenshots.png') +plt.savefig(f'{dPath}/figures/Fig12-{noise_type}_{iter}_noise_floor_by_bvalues_screenshots.png') # Mostrar la gráfica plt.show() diff --git a/metrics/Generate_Fig15_correlation_FA_MD.py b/metrics/Generate_Fig15_correlation_FA_MD.py new file mode 100644 index 0000000..18f3d0c --- /dev/null +++ b/metrics/Generate_Fig15_correlation_FA_MD.py @@ -0,0 +1,73 @@ +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import nibabel as nb +from scipy import ndimage, stats +from scipy.stats import ranksums +import os +import subprocess +from pathlib import Path +import warnings +import sys +import seaborn as sns + +SNRs = [3,5,10,20,40] +methods = ['Raw','NLMEANS','MPPCA','Patch2Self', 'DDM2', 'AVG'] # +methods_names_short = ['Noisy','NLM', 'MPPCA', 'P2S','DDM2', 'AVG'] +list_colors = ['gold','pink', 'blue', 'darkorange', 'lightseagreen', 'green'] +custom_palette = sns.color_palette(list_colors) # To get the same color palette than in seaborn + +iter = 21 +exp = 'Exp6-data-gaussian' +noise_type = 'Gaussian' +metric = 'MD' + +dPath = f'/Simulations/Experiments/{exp}' + + +data = {} + +for method in methods: + print(f'method:{method}') + data[method] = [] + for snr in SNRs: + if (method == 'DDM2' or method == 'Patch2Self') and snr == 3: + continue # Omitir el método DDM2 y P2S con snr=3 + print(f'snr:{snr}') + fsl_cc = np.loadtxt(f'{dPath}/{method}/snr{snr}/analysis/data.dti/fslcc{metric}_res.txt')[2:].ravel() # First 2 columns of the fslcc output are not of interest + fsl_val = 0 + if fsl_cc.size > 0: + fsl_val = fsl_cc[0] + print(f'fsl_cc: {fsl_val}') + data[method].append((snr, fsl_val)) # Guardar pares (snr, fsl_val) + +print(data) + +# Crear la figura y los ejes +fig, ax = plt.subplots(figsize=(5, 5)) + +method_index = 0 +for method in methods: + # Filtrar los valores para no incluir snr=3 para los métodos DDM2 y P2S + if method == 'DDM2' or method == 'Patch2Self': + filtered_data = [(snr, value) for snr, value in data[method] if snr != 3] + filtered_snr = [snr for snr, _ in filtered_data] + filtered_values = [value for _, value in filtered_data] + ax.plot(filtered_snr, filtered_values, marker='o', label=methods_names_short[method_index], color=list_colors[method_index]) + else: + snr_values, values = zip(*data[method]) + ax.plot(snr_values, values, marker='o', label=methods_names_short[method_index], color=list_colors[method_index]) + method_index += 1 + + +ax.set_xlabel('SNR') +ax.set_ylabel('Correlation') +plt.title(metric) +ax.legend() + +plt.tight_layout() +plt.savefig(f'{dPath}/figures/Fig15_{noise_type}{iter}_correlation_FA_MD_fslcc{metric}_screenshots.png') +# Mostrar la gráfica +plt.show() + + diff --git a/metrics/10-Generate_image_signal_distribution_diff_execution.py b/metrics/Generate_Fig17_differences_two_exec_DDM2.py similarity index 71% rename from metrics/10-Generate_image_signal_distribution_diff_execution.py rename to metrics/Generate_Fig17_differences_two_exec_DDM2.py index 985c954..feabd34 100644 --- a/metrics/10-Generate_image_signal_distribution_diff_execution.py +++ b/metrics/Generate_Fig17_differences_two_exec_DDM2.py @@ -26,7 +26,7 @@ def get_data(file, mmap=True): iter = 20 method = 'DDM2' -dPath_gaussian = f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/Exp6-data-gaussian' +dPath_gaussian = f'/Simulations/Experiments/Exp6-data-gaussian' exec1_ddm2_data = get_data(f'{dPath_gaussian}/{method}/snr{snr}/diff_executions/exec1_hardi150_denoised.nii.gz') exec2_ddm2_data = get_data(f'{dPath_gaussian}/{method}/snr{snr}/diff_executions/exec2_hardi150_denoised.nii.gz') @@ -40,15 +40,15 @@ def get_data(file, mmap=True): plt.figure(figsize=(10, 6)) -plt.plot(diff_ddm2_data_voxel1, 'o-', label=f'Voxel 1 - coords:{voxel1_coords}') -plt.plot(diff_ddm2_data_voxel2, 'o-', label=f'Voxel 2 - coords {voxel2_coords}') +plt.plot(diff_ddm2_data_voxel1, 'o-', label=f'Voxel 1') +plt.plot(diff_ddm2_data_voxel2, 'o-', label=f'Voxel 2') -plt.title('Individual voxel signal in the differential') +plt.title('Difference between executions') plt.xlabel('Volumes') plt.ylabel('Signal') plt.legend() -plt.savefig(f'{dPath_gaussian}/figures/{iter}_{method}_voxel_signal_in_differential_screenshots.png', dpi=600, bbox_inches='tight') +plt.savefig(f'{dPath_gaussian}/figures/Fig17-{iter}_{method}_differences_two_exec_DDM2_screenshots.png', dpi=600, bbox_inches='tight') plt.show() diff --git a/metrics/11-Generate-image-signal-by-snr.py b/metrics/Generate_Fig2_signal_intensities_in_slice.py similarity index 70% rename from metrics/11-Generate-image-signal-by-snr.py rename to metrics/Generate_Fig2_signal_intensities_in_slice.py index 436d5a2..f3e2f1f 100644 --- a/metrics/11-Generate-image-signal-by-snr.py +++ b/metrics/Generate_Fig2_signal_intensities_in_slice.py @@ -36,44 +36,41 @@ def get_data(file, mmap=True): SNRs = [5, 10] 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' -ground_truth = get_data(f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/Exp6-data-gaussian/Dataset/noisyfree_data_full_b0_first.nii.gz') +dPath_gaussian = f'/Simulations/Experiments/Exp6-data-gaussian' +dPath_rician = f'/Simulations/Experiments/Exp6-data-rician' +ground_truth = get_data(f'/Simulations/Experiments/Exp6-data-gaussian/Dataset/noisyfree_data_full_b0_first.nii.gz') -# Definir la posición específica +# Position x, y, z = 51, 51, 35 ground_truth_signal_timeseries = ground_truth[x, y, z, :] -# Preparar los gráficos fig, ax = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(10, 4)) - colors = ['orange', 'green'] -# Crear diccionarios para almacenar las series temporales + data_gaussian = {} data_rician = {} -# Iterar sobre los valores de SNR + for snr in SNRs: print(f'snr: {snr}') - # Cargar los datos ruidosos + # Load noisy data gaussian_noisy_data = get_data(f'{dPath_gaussian}/RAW/snr{snr}/main-Gaussian-noisy_data_snr{snr}.nii.gz') rician_noisy_data = get_data(f'{dPath_rician}/RAW/snr{snr}/main-Rician-noisy_data_snr{snr}.nii.gz') - # Extraer la serie temporal de la señal en la posición específica + #Extract the signal in position gaussian_signal_timeseries = gaussian_noisy_data[x, y, z, :] - rician_signal_timeseries = rician_noisy_data[x, y, z, :] + rician_signal_timeseries = rician_noisy_data[x, y, z, :] - # Almacenar las series temporales en los diccionarios data_gaussian[snr] = gaussian_signal_timeseries data_rician[snr] = rician_signal_timeseries ax[0].plot(ground_truth_signal_timeseries, label='noise-free', color='black', linestyle = 'dashed') ax[1].plot(ground_truth_signal_timeseries, label='noise-free', color='black', linestyle = 'dashed') -# Graficar las series temporales + snr_index = 0 for snr in SNRs: ax[0].plot(data_gaussian[snr], label=f'SNR={snr}', color=colors[snr_index]) @@ -92,7 +89,7 @@ def get_data(file, mmap=True): plt.tight_layout() -plt.savefig(f'{dPath_gaussian}/figures/{iter}_signal_intensities_in_position_screenshots.png') +plt.savefig(f'{dPath_gaussian}/figures/{iter}_Fig2_signal_intensities_in_slice_screenshots.png') plt.show() diff --git a/metrics/1-Generate-image-central-slice-methods.py b/metrics/Generate_Fig4-5_central_view_slice.py similarity index 84% rename from metrics/1-Generate-image-central-slice-methods.py rename to metrics/Generate_Fig4-5_central_view_slice.py index 0fcdd44..4d4e2c0 100644 --- a/metrics/1-Generate-image-central-slice-methods.py +++ b/metrics/Generate_Fig4-5_central_view_slice.py @@ -91,23 +91,23 @@ def getVMinMax(method, snr): def printGroundTruth(ground_truth): plt.imshow(ndimage.rotate(ground_truth[:,:,36,65], 90), cmap='gray') - plt.xticks([]) # Oculta los valores en el eje x - plt.yticks([]) # Oculta los valores en el eje y + plt.xticks([]) + plt.yticks([]) plt.savefig(f'{dPath}/figures/{iter}_central_view_slice_ground-truth_screenshots.png') plt.show() -SNRs = [3,5,10,20,40] # -methods = ['Raw','NLMEANS','MPPCA','Patch2Self', 'DDM2', 'AVG'] # -methods_plot_names = ['Noisy','NLMEANS','MPPCA','Patch2Self', 'DDM2', 'AVG'] # +SNRs = [3,5,10,20,40] +methods = ['Raw','NLMEANS','MPPCA','Patch2Self', 'DDM2', 'AVG'] +methods_plot_names = ['Noisy','NLMEANS','MPPCA','Patch2Self', 'DDM2', 'AVG'] list_colors = ['gold', 'pink', 'blue', 'darkorange', 'lightseagreen', 'green'] iter = 20 exp = 'Exp6-data-gaussian' noise_type = 'Gaussian' -ground_truth = get_data(f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}/Dataset/noisyfree_data_full_b0_first.nii.gz') -dPath = f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}' -mask = get_data(f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}/Dataset/nodif_brain_mask.nii.gz') +ground_truth = get_data(f'/Simulations/Experiments/{exp}/Dataset/noisyfree_data_full_b0_first.nii.gz') +dPath = f'/Simulations/Experiments/{exp}' +mask = get_data(f'/Simulations/Experiments/{exp}/Dataset/nodif_brain_mask.nii.gz') #Get ground-truth #print('Printing GroundTruth') @@ -129,7 +129,7 @@ def printGroundTruth(ground_truth): if method == 'Raw': name_data = f'main-{noise_type}-noisy_data_snr' - img = get_data(f'{dPath}/{method}/snr{snr}/{name_data}{snr}.nii.gz') #esta es la imagen denoised o el raw + img = get_data(f'{dPath}/{method}/snr{snr}/{name_data}{snr}.nii.gz') #denoised dataset or raw vmin, vmax = getVMinMax(method, snr) mask_expanded = np.expand_dims(mask, axis=-1) imagex = img * mask_expanded @@ -143,5 +143,5 @@ def printGroundTruth(ground_truth): plt.tight_layout() print('Saving fig') -plt.savefig(f'{dPath}/figures/{iter}_central_view_slice_screenshots.png', dpi=600, bbox_inches='tight') +plt.savefig(f'{dPath}/figures/{iter}_Fig4-5_central_view_slice_screenshots.png', dpi=600, bbox_inches='tight') diff --git a/metrics/1-Generate-image-coronal-slice-methods.py b/metrics/Generate_Fig4-5_coronal_view_slice.py similarity index 84% rename from metrics/1-Generate-image-coronal-slice-methods.py rename to metrics/Generate_Fig4-5_coronal_view_slice.py index 4c90332..6874cf5 100644 --- a/metrics/1-Generate-image-coronal-slice-methods.py +++ b/metrics/Generate_Fig4-5_coronal_view_slice.py @@ -91,23 +91,23 @@ def getVMinMax(method, snr): def printGroundTruth(ground_truth): plt.imshow(ndimage.rotate(ground_truth[:,42,:,61], 90), cmap='gray') - plt.xticks([]) # Oculta los valores en el eje x - plt.yticks([]) # Oculta los valores en el eje y + plt.xticks([]) + plt.yticks([]) plt.savefig(f'{dPath}/figures/{iter}_coronal_view_slice_ground-truth_screenshots.png') plt.show() SNRs = [3,5,10,20,40] -methods = ['Raw','NLMEANS','MPPCA','Patch2Self', 'DDM2', 'AVG'] # -methods_plot_names = ['Noisy','NLMEANS','MPPCA','Patch2Self', 'DDM2', 'AVG'] # +methods = ['Raw','NLMEANS','MPPCA','Patch2Self', 'DDM2', 'AVG'] +methods_plot_names = ['Noisy','NLMEANS','MPPCA','Patch2Self', 'DDM2', 'AVG'] list_colors = ['gold', 'pink', 'blue', 'darkorange', 'lightseagreen', 'green'] iter = 20 exp = 'Exp6-data-rician' noise_type = 'Rician' -ground_truth = get_data(f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}/Dataset/noisyfree_data_full_b0_first.nii.gz') -dPath = f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}' -mask = get_data(f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}/Dataset/nodif_brain_mask.nii.gz') +ground_truth = get_data(f'/Simulations/Experiments/{exp}/Dataset/noisyfree_data_full_b0_first.nii.gz') +dPath = f'/Simulations/Experiments/{exp}' +mask = get_data(f'/Simulations/Experiments/{exp}/Dataset/nodif_brain_mask.nii.gz') #Get ground-truth #print('Printing GroundTruth') @@ -129,7 +129,7 @@ def printGroundTruth(ground_truth): if method == 'Raw': name_data = f'main-{noise_type}-noisy_data_snr' - img = get_data(f'{dPath}/{method}/snr{snr}/{name_data}{snr}.nii.gz') #esta es la imagen denoised o el raw + img = get_data(f'{dPath}/{method}/snr{snr}/{name_data}{snr}.nii.gz') #denoised dataset or raw dataset vmin, vmax = getVMinMax(method, snr) vmin, vmax = getVMinMax(method, snr) @@ -145,5 +145,5 @@ def printGroundTruth(ground_truth): plt.tight_layout() print('Saving fig') -plt.savefig(f'{dPath}/figures/{iter}_coronal_view_slice_screenshots.png', dpi=600, bbox_inches='tight') +plt.savefig(f'{dPath}/figures/{iter}_Fig4-5_coronal_view_slice_screenshots.png', dpi=600, bbox_inches='tight') diff --git a/metrics/7-Generate-image-group-metrics-snr-errors.py b/metrics/Generate_Fig6_reconstruction_error_metrics.py similarity index 73% rename from metrics/7-Generate-image-group-metrics-snr-errors.py rename to metrics/Generate_Fig6_reconstruction_error_metrics.py index 69cd13d..2b83dce 100644 --- a/metrics/7-Generate-image-group-metrics-snr-errors.py +++ b/metrics/Generate_Fig6_reconstruction_error_metrics.py @@ -45,9 +45,7 @@ def calculate_snr(denoised_data_mask): snr = mean_signal / noise_std # Calcula la SNR return snr -def calculate_metric(metric_name, ground_truth, noisy_data, denoised_data): - #difference = noisy_data - denoised_data - #difference_in_mask = difference[mask_positions] +def calculate_metric(metric_name, ground_truth, noisy_data, denoised_data): if metric_name == 'rmse': #se calcula con la imagen original - GroundTruth @@ -64,13 +62,7 @@ def calculate_metric(metric_name, ground_truth, noisy_data, denoised_data): ground_truth_in_mask = ground_truth[mask_positions] denoised_in_mask = denoised_data[mask_positions] difference_grount_denoise_in_mask = ground_truth_in_mask - denoised_in_mask - max_pixel_value = np.max(difference_grount_denoise_in_mask) - #difference_grount_denoise = ground_truth - denoised_data - #difference_grount_denoise_in_mask = difference_grount_denoise[mask_positions] - - #print(f'max_pixel:{max_pixel_value}') - #tipo_de_datos = difference_grount_denoise_in_mask.dtype - #print(f'tipo de datos:{tipo_de_datos}') + max_pixel_value = np.max(difference_grount_denoise_in_mask) psnr = calculate_psnr(difference_grount_denoise_in_mask, max_pixel_value) return psnr if metric_name == 'ssmi': @@ -101,58 +93,65 @@ def calculate_metric(metric_name, ground_truth, noisy_data, denoised_data): method_int = [0, 1, 2, 3] metric_names = ['rmse', 'r2', 'ssmi'] #, 'psnr', 'snr', 'ssmi', 'r2' iter = 20 -exp = 'Exp6-data-rician' -noise_type = 'Rician' +exp = 'Exp6-data-gaussian' +noise_type = 'Gaussian' list_colors = ['pink', 'blue', 'darkorange', 'lightseagreen', 'green'] -mask = get_data(f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}/Dataset/nodif_brain_mask.nii.gz') -ground_truth = get_data(f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}/Dataset/noisyfree_data_full_b0_first.nii.gz') +mask = get_data(f'/Simulations/Experiments/{exp}/Dataset/nodif_brain_mask.nii.gz') +ground_truth = get_data(f'/Simulations/Experiments/{exp}/Dataset/noisyfree_data_full_b0_first.nii.gz') -# Identificar las posiciones donde la máscara tiene un valor de 1 (dentro del cerebro) +# Identify the positions where the mask has value = 1 (inside the brain) mask_positions = np.where(mask == 1) -dPath = f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}' +dPath = f'/Simulations/Experiments/{exp}' fig, axs = plt.subplots(1, 3, figsize=(15, 5)) # Configurar tamaños de fuente específicos plt.rcParams.update({ - 'font.size': 11, # Tamaño de fuente por defecto - 'axes.titlesize': 12, # Tamaño de fuente del título del gráfico - 'axes.labelsize': 12, # Tamaño de fuente de las etiquetas de los ejes - 'xtick.labelsize': 11, # Tamaño de fuente de las etiquetas de las marcas del eje x - 'ytick.labelsize': 11 # Tamaño de fuente de las etiquetas de las marcas del eje y + 'font.size': 11, + 'axes.titlesize': 12, + 'axes.labelsize': 12, + 'xtick.labelsize': 11, + 'ytick.labelsize': 11 }) - for ds, metric in enumerate(metric_names): median_difference_methods = [] for method in methods: print(f'method:{method}') differencebymethod = [] for snr in SNRs: + if (method == 'DDM2' or method == 'Patch2Self') and snr == 3: + continue # Omitir el método DDM2 y P2S con snr=3 print(f'snr:{snr}') noisy_data = get_data(f'{dPath}/RAW/snr{snr}/main-{noise_type}-noisy_data_snr{snr}.nii.gz') denoised_data = get_data(f'{dPath}/{method}/snr{snr}/{method.lower()}-denoised_main_snr{snr}.nii.gz') difference = calculate_metric(metric, ground_truth, noisy_data, denoised_data) print(f'difference: {difference}') - differencebymethod.append(difference) + differencebymethod.append((snr, difference)) # Guardar pares (snr, diferencia) median_difference_methods.append(differencebymethod) print(differencebymethod) print(median_difference_methods) - #row = ds // 2 - #col = ds % 2 - # Graficar los errores de cada método como una línea method_index = 0 for method in methods: - axs[ds].plot(SNRs, median_difference_methods[method_index], marker='o', label=method, color=list_colors[method_index]) - method_index = method_index + 1 + if method == 'DDM2' or method == 'Patch2Self': + # Filtrar los valores para no incluir snr=3 + filtered_data = [(snr, difference) for snr, difference in median_difference_methods[method_index] if snr != 3] + filtered_snr = [snr for snr, _ in filtered_data] + filtered_differences = [difference for _, difference in filtered_data] + axs[ds].plot(filtered_snr, filtered_differences, marker='o', label=method, color=list_colors[method_index]) + else: + snr_values, differences = zip(*median_difference_methods[method_index]) + axs[ds].plot(snr_values, differences, marker='o', label=method, color=list_colors[method_index]) + method_index += 1 + + - # Agregar etiquetas, título y leyenda y_name = metric.upper() axs[ds].set_xlabel('SNR') axs[ds].set_ylabel('') @@ -161,9 +160,9 @@ def calculate_metric(metric_name, ground_truth, noisy_data, denoised_data): axs[2].legend() -# Mostrar la gráfica + plt.tight_layout() -plt.savefig(f'{dPath}/figures/{iter}_metrics_group.png', dpi=600, bbox_inches='tight') +plt.savefig(f'{dPath}/figures/{noise_type}_{iter}_Fig6_reconstruction_error_metrics.png', dpi=600, bbox_inches='tight') plt.show() diff --git a/metrics/12-Generate_image_signal_bymethods_invoxel.py b/metrics/Generate_Fig7-8_signal_intensities_voxel_fibers.py similarity index 83% rename from metrics/12-Generate_image_signal_bymethods_invoxel.py rename to metrics/Generate_Fig7-8_signal_intensities_voxel_fibers.py index 5c7c126..0c37359 100644 --- a/metrics/12-Generate_image_signal_bymethods_invoxel.py +++ b/metrics/Generate_Fig7-8_signal_intensities_voxel_fibers.py @@ -52,8 +52,8 @@ def get_voxel_coords(voxel_type): x, y, z = get_voxel_coords(voxel_type) -ground_truth = get_data(f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}/Dataset/noisyfree_data_full_b0_first.nii.gz') -dPath = f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}' +ground_truth = get_data(f'/Simulations/Experiments/{exp}/Dataset/noisyfree_data_full_b0_first.nii.gz') +dPath = f'/Simulations/Experiments/{exp}' signal_timeseries_ground_truth = ground_truth[x, y, z] @@ -65,7 +65,7 @@ def get_voxel_coords(voxel_type): if method == 'Raw': name_data = f'main-{noise_type}-noisy_data_snr' - img = get_data(f'{dPath}/{method}/snr{snr}/{name_data}{snr}.nii.gz') #esta es la imagen denoised o el raw + img = get_data(f'{dPath}/{method}/snr{snr}/{name_data}{snr}.nii.gz') #denoised dataset or raw dataset signal_timeseries = img[x, y, z] data_in_voxel[method] = signal_timeseries @@ -74,7 +74,7 @@ def get_voxel_coords(voxel_type): ax.plot(signal_timeseries_ground_truth, label='noise free', color='black', linestyle = 'dashed') -# Graficar las series temporales + method_index = 0 for method in methods: ax.plot(data_in_voxel[method], label=methods_plot_names[method_index], color=list_colors[method_index]) @@ -86,6 +86,6 @@ def get_voxel_coords(voxel_type): plt.legend() plt.title(f'SNR = {snr}') -plt.savefig(f'{dPath}/figures/{noise_type}_{voxel_type}_{iter}_signal_bymethod_invoxel_snr{snr}_screenshots.png', dpi=600, bbox_inches='tight') +plt.savefig(f'{dPath}/figures/Fig7-8_{noise_type}_{voxel_type}_{iter}_signal_intensities_voxels_fibers_bymethod_snr{snr}_screenshots.png', dpi=600, bbox_inches='tight') plt.show() diff --git a/metrics/4-Generate-image4-bymethod-snr.py b/metrics/Generate_Fig9-difference_raw_method-snr.py similarity index 85% rename from metrics/4-Generate-image4-bymethod-snr.py rename to metrics/Generate_Fig9-difference_raw_method-snr.py index cd491b6..dc1ce82 100644 --- a/metrics/4-Generate-image4-bymethod-snr.py +++ b/metrics/Generate_Fig9-difference_raw_method-snr.py @@ -38,12 +38,12 @@ def get_data(file, mmap=True): exp = 'Exp6-data-rician' noise_type = 'rician' -mask = get_data(f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}/Dataset/nodif_brain_mask.nii.gz') -ground_truth = get_data(f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}/Dataset/noisyfree_data_full_b0_first.nii.gz') +mask = get_data(f'/Simulations/Experiments/{exp}/Dataset/nodif_brain_mask.nii.gz') +ground_truth = get_data(f'/Simulations/Experiments/{exp}/Dataset/noisyfree_data_full_b0_first.nii.gz') # Identificar las posiciones donde la máscara tiene un valor de 1 (dentro del cerebro) mask_positions = np.where(mask == 1) -dPath = f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}' +dPath = f'/Simulations/Experiments/{exp}' for ds, dataset in enumerate(SNRs): @@ -76,10 +76,11 @@ def get_data(file, mmap=True): ax[1, i].set_xlim(-0.8, 0.8) plt.setp(ax[1,i], yticks=[]) + print('----') fig.suptitle(f'Dataset SNR={SNRs[ds]}', fontweight='bold', size=fig.dpi*0.4) plt.subplots_adjust(wspace=0.0, hspace=-0.25) fig.tight_layout() plt.show() - fig.savefig(f'{dPath}/figures/{noise_type}_{iter}_fig4_snr{dataset}_screenshots.png', dpi=600, bbox_inches='tight') \ No newline at end of file + fig.savefig(f'{dPath}/figures/{noise_type}_{iter}_Fig9_snr{dataset}_screenshots.png', dpi=600, bbox_inches='tight') \ No newline at end of file diff --git a/metrics/verify-correlation.py b/metrics/verify-correlation.py new file mode 100644 index 0000000..6bc76e4 --- /dev/null +++ b/metrics/verify-correlation.py @@ -0,0 +1,31 @@ +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import nibabel as nb +from scipy import ndimage, stats +from scipy.stats import ranksums +import os +import subprocess +from pathlib import Path +import warnings +import sys +import seaborn as sns + + + +iter = 21 +exp = 'Exp6-data-rician' +noise_type = 'Rician' + + +dPath = f'C:/Users/dayri/Documents/UNED/TFM/Related_projects/Simulations/Simulations/Experiments/{exp}' + +fsl_cc = np.loadtxt(f'{dPath}/AVG/snr3/analysis/dMRI/processed_data/fslcc_with_noisyfree_res.txt', usecols=2)[3:].ravel() # First 2 columns of the fslcc output are not of interest +median_value = np.median(fsl_cc) +print(f'fslcc_with_noisyfree_res: {median_value}') + + +fsl_cc_mppca = np.loadtxt(f'{dPath}/MPPCA/snr3/analysis/dMRI/processed_data/fslcc_with_noisyfree_res.txt', usecols=2)[3:].ravel() # First 2 columns of the fslcc output are not of interest +median_value_mppca = np.median(fsl_cc_mppca) +print(f'fslcc_with_noisyfree_res: {median_value_mppca}') +