add and modify metrics

metrics/10-Generate_image_signal_distribution_diff_execution.py

@@ -44,7 +44,7 @@ def get_data(file, mmap=True):
 plt.plot(diff_ddm2_data_voxel2, 'o-', label=f'Voxel 2 - coords {voxel2_coords}')
 
 plt.title('Individual voxel signal in the differential')
-plt.xlabel('b-values')
+plt.xlabel('Volumes')
 plt.ylabel('Signal')
 plt.legend()
 

metrics/11-Generate-image-signal-by-snr.py

@@ -0,0 +1,101 @@
+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
+
+
+
+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')
+
+# Definir la posición específica
+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
+    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
+    gaussian_signal_timeseries = gaussian_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])
+    ax[1].plot(data_rician[snr], label=f'SNR={snr}', color=colors[snr_index])
+    snr_index = snr_index + 1
+
+# Configurar los gráficos
+ax[0].set_title('Gaussian Noisy Data')
+ax[0].set_xlabel('')
+ax[0].set_ylabel('Signal Intensity')
+ax[0].legend()
+
+ax[1].set_title('Rician Noisy Data')
+ax[1].set_xlabel('')
+ax[1].legend()
+
+plt.tight_layout()
+
+plt.savefig(f'{dPath_gaussian}/figures/{iter}_signal_intensities_in_position_screenshots.png')
+plt.show()
+
+
+
+
+

metrics/12-Generate_image_signal_bymethods_invoxel.py

@@ -0,0 +1,91 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt 
+import nibabel as nb
+import seaborn as sns
+
+
+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 get_voxel_coords(voxel_type):
+    '''
+    mean_f1samples (f1), mean_f2samples (f2), mean_f3samples (f3). 
+    Para 1fib, coge cualquier voxel donde f2 y f3 <0.05
+    Para 2fib, cualquier voxel donde f3<0.05 y f2>0.05
+    Para 3fib, cualquier voxel donde f3>0.05.
+
+    Los valores que se usan cumplen los requisitos de arriba.
+    '''
+    if voxel_type == '1fib':
+        return (39, 71, 47)
+    elif voxel_type == '2fib':
+        return (36, 25, 29)
+    elif voxel_type == '3fib':
+        return (76, 35, 29)
+    else:
+        raise ValueError("Unknown voxel type")
+
+
+snr = 5 #
+methods = ['Raw','Patch2Self', 'DDM2'] #'Raw','NLMEANS','MPPCA','Patch2Self', 'DDM2', 'AVG'
+methods_plot_names = ['Noisy','Patch2Self', 'DDM2'] #'Noisy','NLMEANS','MPPCA','Patch2Self', 'DDM2', 'AVG'
+list_colors = ['gold', 'darkorange', 'lightseagreen', 'green']#'gold', 'pink', 'blue', 'darkorange', 'lightseagreen', 'green'
+
+iter = 20
+exp = 'Exp6-data-rician'
+noise_type = 'Rician'
+voxel_type = '1fib'
+
+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}'
+
+signal_timeseries_ground_truth = ground_truth[x, y, z]
+
+data_in_voxel = {}
+
+for ds, method in enumerate(methods):
+    print(f'method: {method}')
+    name_data = f'{method.lower()}-denoised_main_snr'
+    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 
+    signal_timeseries = img[x, y, z]
+    data_in_voxel[method] = signal_timeseries
+
+
+fig, ax = plt.subplots(figsize=(5, 5))
+
+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])
+    method_index = method_index +1
+
+
+plt.xlabel('Volumes')
+plt.ylabel('Signal intensity')
+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.show()
+

metrics/4-Generate-image4-bymethod-snr.py

@@ -0,0 +1,85 @@
+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
+
+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
+
+
+SNRs = [5]
+methods = ['NLMEANS','MPPCA','Patch2Self', 'DDM2', 'AVG']
+methods_names = ['NLMEANS','MPPCA','Patch2Self', 'DDM2', 'AVG']
+list_colors = ['pink', 'blue', 'darkorange', 'lightseagreen', 'green']
+
+iter = 20
+export_figures = 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')
+
+# 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}'
+
+
+for ds, dataset in enumerate(SNRs):
+    nrows = 2
+    ncols = len(methods_names)
+    fig, ax = plt.subplots(nrows=nrows, ncols=ncols, figsize=(7*ncols, 6*nrows), sharex='row', sharey='row', gridspec_kw={'height_ratios': [3, 1]})    
+
+    for i in range(0,ncols):
+        print(dataset)        
+        method = methods_names[i]
+        print(method)
+        img = get_data(f'{dPath}/{method}/snr{dataset}/differences.nii.gz')  #difference between raw and denoised data
+        img_volume = img[:,:,:,69]
+        masked_volume = img_volume[mask>0].ravel()          #535k voxels at 1.5mm
+        img_slice = img[:,:,36,69]
+
+        # Exemplar slice map masked to brain
+        masked_slice = img_slice*mask[:,:,36]
+        v_max = 3*np.std(masked_slice)
+        ax_imshow = ax[0,i]
+        ax_imshow.imshow(ndimage.rotate(masked_slice, 90), cmap='gray', vmin=-v_max, vmax=v_max)
+        ax_imshow.axis('off')
+        ax_imshow.set_title(method, fontweight='bold', size=fig.dpi*0.3, color=list_colors[i])    # Size of titles in function of the figsize
+        plt.setp(ax_imshow, xticks=[], yticks=[])
+
+        sns.kdeplot(masked_volume,color=list_colors[i], ax=ax[1,i], fill=True)#, bins=200)
+        ax[1,i].axvline(x=np.mean(masked_volume),color=list_colors[i], ls='--', lw=2)
+        ax[1,i].axvline(x=0,color='black', ls='--', lw=1)
+        #ax[1,i].set_xlim(left=x_limit, right=x_limit)
+        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')