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