1""" 2======================================================= 3Denoise images using Local PCA via empirical thresholds 4======================================================= 5 6PCA-based denoising algorithms are effective denoising methods because they 7explore the redundancy of the multi-dimensional information of 8diffusion-weighted datasets. In this example, we will show how to 9perform the PCA-based denoising using the algorithm proposed by Manjon et al. 10[Manjon2013]_. 11 12This algorithm involves the following steps: 13 14* First, we estimate the local noise variance at each voxel. 15 16* Then, we apply PCA in local patches around each voxel over the gradient 17 directions. 18 19* Finally, we threshold the eigenvalues based on the local estimate of sigma 20 and then do a PCA reconstruction 21 22 23To perform PCA denoising without a prior noise standard deviation estimate 24please see the following example instead: :ref:`denoise_mppca` 25 26Let's load the necessary modules 27""" 28 29import numpy as np 30import nibabel as nib 31import matplotlib.pyplot as plt 32from time import time 33from dipy.core.gradients import gradient_table 34from dipy.denoise.localpca import localpca 35from dipy.denoise.pca_noise_estimate import pca_noise_estimate 36from dipy.data import get_fnames 37from dipy.io.image import load_nifti 38from dipy.io.gradients import read_bvals_bvecs 39 40""" 41 42Load one of the datasets. These data were acquired with 63 gradients and 1 43non-diffusion (b=0) image. 44 45""" 46 47dwi_fname, dwi_bval_fname, dwi_bvec_fname = get_fnames('isbi2013_2shell') 48data, affine = load_nifti(dwi_fname) 49bvals, bvecs = read_bvals_bvecs(dwi_bval_fname, dwi_bvec_fname) 50gtab = gradient_table(bvals, bvecs) 51 52print("Input Volume", data.shape) 53 54""" 55Estimate the noise standard deviation 56===================================== 57 58We use the ``pca_noise_estimate`` method to estimate the value of sigma to be 59used in local PCA algorithm proposed by Manjon et al. [Manjon2013]_. 60It takes both data and the gradient table object as input and returns an 61estimate of local noise standard deviation as a 3D array. We return a smoothed 62version, where a Gaussian filter with radius 3 voxels has been applied to the 63estimate of the noise before returning it. 64 65We correct for the bias due to Rician noise, based on an equation developed by 66Koay and Basser [Koay2006]_. 67 68""" 69 70t = time() 71sigma = pca_noise_estimate(data, gtab, correct_bias=True, smooth=3) 72print("Sigma estimation time", time() - t) 73 74""" 75Perform the localPCA using the function `localpca` 76================================================== 77 78The localpca algorithm takes into account the multi-dimensional information of 79the diffusion MR data. It performs PCA on local 4D patch and 80then removes the noise components by thresholding the lowest eigenvalues. 81The eigenvalue threshold will be computed from the local variance estimate 82performed by the ``pca_noise_estimate`` function, if this is inputted in the 83main ``localpca`` function. The relationship between the noise variance 84estimate and the eigenvalue threshold can be adjusted using the input parameter 85``tau_factor``. According to Manjon et al. [Manjon2013]_, this parameter is set 86to 2.3. 87""" 88 89t = time() 90 91denoised_arr = localpca(data, sigma, tau_factor=2.3, patch_radius=2) 92 93print("Time taken for local PCA (slow)", -t + time()) 94 95""" 96The ``localpca`` function returns the denoised data which is plotted below 97(middle panel) together with the original version of the data (left panel) and 98the algorithm residual (right panel) . 99""" 100 101sli = data.shape[2] // 2 102gra = data.shape[3] // 2 103orig = data[:, :, sli, gra] 104den = denoised_arr[:, :, sli, gra] 105rms_diff = np.sqrt((orig - den) ** 2) 106 107fig, ax = plt.subplots(1, 3) 108ax[0].imshow(orig, cmap='gray', origin='lower', interpolation='none') 109ax[0].set_title('Original') 110ax[0].set_axis_off() 111ax[1].imshow(den, cmap='gray', origin='lower', interpolation='none') 112ax[1].set_title('Denoised Output') 113ax[1].set_axis_off() 114ax[2].imshow(rms_diff, cmap='gray', origin='lower', interpolation='none') 115ax[2].set_title('Residual') 116ax[2].set_axis_off() 117plt.savefig('denoised_localpca.png', bbox_inches='tight') 118 119print("The result saved in denoised_localpca.png") 120 121""" 122.. figure:: denoised_localpca.png 123 :align: center 124 125Below we show how the denoised data can be saved. 126""" 127 128nib.save(nib.Nifti1Image(denoised_arr, 129 affine), 'denoised_localpca.nii.gz') 130 131print("Entire denoised data saved in denoised_localpca.nii.gz") 132 133""" 134.. [Manjon2013] Manjon JV, Coupe P, Concha L, Buades A, Collins DL "Diffusion 135 Weighted Image Denoising Using Overcomplete Local PCA" (2013). 136 PLoS ONE 8(9): e73021. doi:10.1371/journal.pone.0073021. 137 138.. [Koay2006] Koay CG, Basser PJ (2006). "Analytically exact correction scheme 139 for signal extraction from noisy magnitude MR signals". JMR 179: 140 317-322. 141 142.. include:: ../links_names.inc 143""" 144