【DMRI】如何利用dipy进行dmri的patch2self去噪处理，并将去噪后的数据存储为nii的格式

patch2self是目前dMRI数据去噪比较好用的方法，去噪效果较好，并且在cpu上就可以运行，dipy提供了现成的函数接口

实现代码

1. 导入必要的库

   from dipy.denoise.patch2self import patch2self
   from dipy.io.image import load_nifti, save_nifti
   from dipy.io.gradients import read_bvals_bvecs
   import os
   import numpy as np
   import matplotlib.pyplot as plt
   

2. 定义加载数据并进行去噪的函数

   def loading_data(data_dir, saving_dir, is_rewrite=False):
       # 读取b值和b向量
       bvals, bvecs = read_bvals_bvecs('path/to/bval/file.bval', 'path/to/bvec/file.bvec')
       
       # 如果已存在去噪数据且不需要重写，则加载已存在的数据
       if os.path.exists(f"{saving_dir}/data.nii.gz") and not is_rewrite:
           print(f"loading exist denoised data")
           # denoised_data, affine = load_nifti(f"{saving_dir}")
       else:
           # 如果不存在去噪数据或者需要重写，则进行去噪处理
           if not os.path.exists(saving_dir):
               os.makedirs(saving_dir)
           print(f"loading noisy data")
           data, affine = load_nifti(f"{data_dir}")
           print(f"denosing data", data.shape)
   
           if data.shape[-1] != 61:
               return
   
           # 选择特定b值的索引
           idx_0 = np.abs(bvals - 0) <= 100
           idx_1000 = np.abs(bvals - 1000) <= 100
           idx_2000 = np.abs(bvals - 2000) <= 100
           # idx_3000 = np.abs(bvals - 3000) <= 100
   
           # 对不同b值的数据进行去噪
           denoised_data_0 = patch2self(data[..., idx_0], bvals[idx_0], model="ols", b0_threshold=10, patch_radius=(1,1,1), clip_negative_vals=False, shift_intensity=True)
           denoised_data_1000 = patch2self(data[..., idx_1000], bvals[idx_1000], model="ols", b0_threshold=100, patch_radius=(1,1,1), clip_negative_vals=False, shift_intensity=True)
           denoised_data_2000 = patch2self(data[..., idx_2000], bvals[idx_2000], model="ols", b0_threshold=100, patch_radius=(1,1,1), clip_negative_vals=False, shift_intensity=True)
           # denoised_data_3000 = patch2self(data[..., idx_3000], bvals[idx_3000], model="ols", b0_threshold=100, clip_negative_vals=False, shift_intensity=True)
   
           # 将去噪后的数据存储到一个新的4D数组中
           denoised_data = np.empty_like(data)
           denoised_data[..., idx_0] = denoised_data_0
           denoised_data[..., idx_1000] = denoised_data_1000
           denoised_data[..., idx_2000] = denoised_data_2000
           # denoised_data[..., idx_3000] = denoised_data_3000
   
           print(f"saving denoised data")
           print("saving_dir:", saving_dir)
           save_nifti(f"{saving_dir}/data.nii.gz", denoised_data, affine)
   
       return denoised_data, bvals, bvecs
   

3. 定义显示图像的函数

   def show_imgs(data, bvals, bvecs):
       print(f"############### show img #################")
       idx = np.abs(bvals - 1000) <= 100
       data_1000 = data[..., idx]
       print(f"data b=1000 shape: {data_1000.shape}")
       img = data_1000[:, :, data_1000.shape[-2] // 2, 0]
       plt.imshow(img, cmap="gray")
       plt.show()
   

4. 定义主函数

   def main():
       print(f"start running...")
       base_dir = "path/to/base/directory"
       denoised_dir = "path/to/denoised/directory"
   
       print(f"############## data ##############")
       data_type_dir = f"{base_dir}/"
       data_list = os.listdir(data_type_dir)
       for data_name in data_list:
           print(f"############## {data_name} data ##############")
           data_dir = f"{data_type_dir}/{data_name}"
           saving_dir = f"{denoised_dir}/{data_name}"
           denoised_data, bvals, bvecs = loading_data(data_dir, saving_dir)
           # show_imgs(denoised_data, bvals, bvecs)
   

5. 入口判断

   if __name__ == "__main__":
       main()