Note
Go to the end to download the full example code
Denoise images using the Marcenko-Pastur PCA algorithm#
The PCA-based denoising algorithm exploits the redundancy across the diffusion-weighted images [Manjon2013], [Veraart2016a]. This algorithm has been shown to provide an optimal compromise between noise suppression and loss of anatomical information for different techniques such as DTI [Manjon2013], spherical deconvolution [Veraart2016a] and DKI [Henri2018].
The basic idea behind the PCA-based denoising algorithms is to remove the components of the data that are classified as noise. The Principal Components classification can be performed based on prior noise variance estimates [Manjon2013] (see denoise_localpca) # noqa E501 or automatically based on the Marchenko-Pastur distribution [Veraa2016a]. In addition to noise suppression, the PCA algorithm can be used to get the standard deviation of the noise [Veraa2016b].
In the following example, we show how to denoise diffusion MRI images and estimate the noise standard deviation using the PCA algorithm based on the Marcenko-Pastur distribution [Veraa2016a]
Let’s load the necessary modules
# load general modules
import numpy as np
import matplotlib.pyplot as plt
from time import time
# load main pca function using Marcenko-Pastur distribution
from dipy.denoise.localpca import mppca
# load functions to fetch data for this example
from dipy.data import get_fnames
# load other dipy's functions that will be used for auxiliary analysis
from dipy.core.gradients import gradient_table
from dipy.io.image import load_nifti, save_nifti
from dipy.io.gradients import read_bvals_bvecs
from dipy.segment.mask import median_otsu
import dipy.reconst.dki as dki
For this example, we use fetch to download a multi-shell dataset which was kindly provided by Hansen and Jespersen (more details about the data are provided in their paper [Hansen2016]). The total size of the downloaded data is 192 MBytes, however you only need to fetch it once.
dwi_fname, dwi_bval_fname, dwi_bvec_fname, _ = get_fnames('cfin_multib')
data, affine = load_nifti(dwi_fname)
bvals, bvecs = read_bvals_bvecs(dwi_bval_fname, dwi_bvec_fname)
gtab = gradient_table(bvals, bvecs)
For the sake of simplicity, we only select two non-zero b-values for this example.
bvals = gtab.bvals
bvecs = gtab.bvecs
sel_b = np.logical_or(np.logical_or(bvals == 0, bvals == 1000), bvals == 2000)
data = data[..., sel_b]
gtab = gradient_table(bvals[sel_b], bvecs[sel_b])
print(data.shape)
(96, 96, 19, 67)
As one can see from its shape, the selected data contains a total of 67 volumes of images corresponding to all the diffusion gradient directions of the selected b-values.
The PCA denoising using the Marchenko-Pastur distribution can be performed by
calling the function mppca:
Time taken for local MP-PCA 145.7494671344757
Internally, the mppca algorithm denoises the diffusion-weighted data
using a 3D sliding window which is defined by the variable patch_radius.
In total, this window should comprise a larger number of voxels than the
number of diffusion-weighted volumes. Since our data has a total of 67
volumes, the patch_radius is set to 2 which corresponds to a 5x5x5 sliding
window comprising a total of 125 voxels.
# To assess the performance of the Marchenko-Pastur PCA denoising algorithm,
# an axial slice of the original data, denoised data, and residuals are
# plotted below:
sli = data.shape[2] // 2
gra = data.shape[3] - 1
orig = data[:, :, sli, gra]
den = denoised_arr[:, :, sli, gra]
rms_diff = np.sqrt((orig - den) ** 2)
fig1, ax = plt.subplots(1, 3, figsize=(12, 6),
subplot_kw={'xticks': [], 'yticks': []})
fig1.subplots_adjust(hspace=0.3, wspace=0.05)
ax.flat[0].imshow(orig.T, cmap='gray', interpolation='none',
origin='lower')
ax.flat[0].set_title('Original')
ax.flat[1].imshow(den.T, cmap='gray', interpolation='none',
origin='lower')
ax.flat[1].set_title('Denoised Output')
ax.flat[2].imshow(rms_diff.T, cmap='gray', interpolation='none',
origin='lower')
ax.flat[2].set_title('Residuals')
fig1.savefig('denoised_mppca.png')
The noise suppression can be visually appreciated by comparing the original data slice (left panel) to its denoised version (middle panel). The difference between original and denoised data showing only random noise indicates that the data’s structural information is preserved by the PCA denoising algorithm (right panel).
Below we show how the denoised data can be saved.
save_nifti('denoised_mppca.nii.gz', denoised_arr, affine)
Additionally, we show how the PCA denoising algorithm affects different diffusion measurements. For this, we run the diffusion kurtosis model below on both original and denoised versions of the data:
dkimodel = dki.DiffusionKurtosisModel(gtab)
maskdata, mask = median_otsu(data, vol_idx=[0, 1], median_radius=4, numpass=2,
autocrop=False, dilate=1)
dki_orig = dkimodel.fit(data, mask=mask)
dki_den = dkimodel.fit(denoised_arr, mask=mask)
0%| | 0/57345 [00:00<?, ?it/s]
1%|▉ | 422/57345 [00:00<00:13, 4187.58it/s]
1%|█▊ | 841/57345 [00:00<00:14, 3905.14it/s]
2%|██▋ | 1233/57345 [00:00<00:14, 3773.85it/s]
3%|███▌ | 1646/57345 [00:00<00:14, 3903.95it/s]
4%|████▍ | 2062/57345 [00:00<00:13, 3980.88it/s]
4%|█████▍ | 2522/57345 [00:00<00:13, 4173.79it/s]
5%|██████▍ | 2985/57345 [00:00<00:12, 4314.73it/s]
6%|███████▍ | 3436/57345 [00:00<00:12, 4370.24it/s]
7%|████████▍ | 3915/57345 [00:00<00:11, 4492.95it/s]
8%|█████████▍ | 4389/57345 [00:01<00:11, 4562.83it/s]
8%|██████████▍ | 4851/57345 [00:01<00:11, 4574.40it/s]
9%|███████████▍ | 5309/57345 [00:01<00:11, 4516.47it/s]
10%|████████████▌ | 5795/57345 [00:01<00:11, 4615.00it/s]
11%|█████████████▌ | 6258/57345 [00:01<00:11, 4616.10it/s]
12%|██████████████▌ | 6736/57345 [00:01<00:10, 4660.23it/s]
13%|███████████████▌ | 7212/57345 [00:01<00:10, 4683.99it/s]
13%|████████████████▋ | 7692/57345 [00:01<00:10, 4712.52it/s]
14%|█████████████████▋ | 8176/57345 [00:01<00:10, 4744.25it/s]
15%|██████████████████▋ | 8664/57345 [00:01<00:10, 4778.43it/s]
16%|███████████████████▊ | 9150/57345 [00:02<00:10, 4797.93it/s]
17%|████████████████████▊ | 9640/57345 [00:02<00:09, 4815.23it/s]
18%|█████████████████████▋ | 10127/57345 [00:02<00:09, 4825.89it/s]
19%|██████████████████████▊ | 10618/57345 [00:02<00:09, 4843.96it/s]
19%|███████████████████████▊ | 11103/57345 [00:02<00:09, 4838.24it/s]
20%|████████████████████████▊ | 11587/57345 [00:02<00:09, 4803.35it/s]
21%|█████████████████████████▉ | 12081/57345 [00:02<00:09, 4837.06it/s]
22%|██████████████████████████▉ | 12567/57345 [00:02<00:09, 4836.78it/s]
23%|████████████████████████████ | 13060/57345 [00:02<00:09, 4858.33it/s]
24%|█████████████████████████████ | 13547/57345 [00:02<00:09, 4854.43it/s]
24%|██████████████████████████████ | 14033/57345 [00:03<00:08, 4849.42it/s]
25%|███████████████████████████████▏ | 14518/57345 [00:03<00:08, 4828.82it/s]
26%|████████████████████████████████▏ | 15001/57345 [00:03<00:08, 4809.23it/s]
27%|█████████████████████████████████▏ | 15497/57345 [00:03<00:08, 4847.27it/s]
28%|██████████████████████████████████▎ | 15994/57345 [00:03<00:08, 4880.08it/s]
29%|███████████████████████████████████▎ | 16483/57345 [00:03<00:08, 4874.30it/s]
30%|████████████████████████████████████▍ | 16971/57345 [00:03<00:08, 4853.52it/s]
30%|█████████████████████████████████████▍ | 17469/57345 [00:03<00:08, 4885.27it/s]
31%|██████████████████████████████████████▌ | 17962/57345 [00:03<00:08, 4892.06it/s]
32%|███████████████████████████████████████▌ | 18459/57345 [00:03<00:07, 4901.16it/s]
33%|████████████████████████████████████████▋ | 18950/57345 [00:04<00:07, 4897.16it/s]
34%|█████████████████████████████████████████▋ | 19440/57345 [00:04<00:07, 4894.84it/s]
35%|██████████████████████████████████████████▋ | 19930/57345 [00:04<00:07, 4880.20it/s]
36%|███████████████████████████████████████████▊ | 20419/57345 [00:04<00:07, 4878.60it/s]
36%|████████████████████████████████████████████▊ | 20907/57345 [00:04<00:07, 4844.19it/s]
37%|█████████████████████████████████████████████▉ | 21392/57345 [00:04<00:07, 4841.77it/s]
38%|██████████████████████████████████████████████▉ | 21882/57345 [00:04<00:07, 4853.25it/s]
39%|███████████████████████████████████████████████▉ | 22374/57345 [00:04<00:07, 4866.14it/s]
40%|█████████████████████████████████████████████████ | 22864/57345 [00:04<00:07, 4871.20it/s]
41%|██████████████████████████████████████████████████ | 23352/57345 [00:04<00:06, 4867.71it/s]
42%|███████████████████████████████████████████████████▏ | 23839/57345 [00:05<00:06, 4849.07it/s]
42%|████████████████████████████████████████████████████▏ | 24330/57345 [00:05<00:06, 4862.02it/s]
43%|█████████████████████████████████████████████████████▏ | 24820/57345 [00:05<00:06, 4867.69it/s]
44%|██████████████████████████████████████████████████████▎ | 25307/57345 [00:05<00:06, 4861.34it/s]
45%|███████████████████████████████████████████████████████▎ | 25802/57345 [00:05<00:06, 4881.41it/s]
46%|████████████████████████████████████████████████████████▍ | 26291/57345 [00:05<00:06, 4862.61it/s]
47%|█████████████████████████████████████████████████████████▍ | 26780/57345 [00:05<00:06, 4869.41it/s]
48%|██████████████████████████████████████████████████████████▍ | 27267/57345 [00:05<00:06, 4864.21it/s]
48%|███████████████████████████████████████████████████████████▌ | 27754/57345 [00:05<00:06, 4859.71it/s]
49%|████████████████████████████████████████████████████████████▌ | 28244/57345 [00:05<00:05, 4866.34it/s]
50%|█████████████████████████████████████████████████████████████▋ | 28744/57345 [00:06<00:05, 4901.09it/s]
51%|██████████████████████████████████████████████████████████████▋ | 29235/57345 [00:06<00:05, 4896.48it/s]
52%|███████████████████████████████████████████████████████████████▊ | 29725/57345 [00:06<00:05, 4861.99it/s]
53%|████████████████████████████████████████████████████████████████▊ | 30212/57345 [00:06<00:05, 4859.93it/s]
54%|█████████████████████████████████████████████████████████████████▊ | 30699/57345 [00:06<00:05, 4855.22it/s]
54%|██████████████████████████████████████████████████████████████████▉ | 31187/57345 [00:06<00:05, 4856.26it/s]
55%|███████████████████████████████████████████████████████████████████▉ | 31678/57345 [00:06<00:05, 4866.70it/s]
56%|████████████████████████████████████████████████████████████████████▉ | 32169/57345 [00:06<00:05, 4872.11it/s]
57%|██████████████████████████████████████████████████████████████████████ | 32666/57345 [00:06<00:05, 4893.49it/s]
58%|███████████████████████████████████████████████████████████████████████ | 33156/57345 [00:06<00:04, 4887.74it/s]
59%|████████████████████████████████████████████████████████████████████████▏ | 33645/57345 [00:07<00:04, 4879.50it/s]
60%|█████████████████████████████████████████████████████████████████████████▏ | 34133/57345 [00:07<00:04, 4879.59it/s]
60%|██████████████████████████████████████████████████████████████████████████▎ | 34628/57345 [00:07<00:04, 4894.62it/s]
61%|███████████████████████████████████████████████████████████████████████████▎ | 35120/57345 [00:07<00:04, 4894.95it/s]
62%|████████████████████████████████████████████████████████████████████████████▍ | 35610/57345 [00:07<00:04, 4891.09it/s]
63%|█████████████████████████████████████████████████████████████████████████████▍ | 36100/57345 [00:07<00:04, 4872.61it/s]
64%|██████████████████████████████████████████████████████████████████████████████▍ | 36590/57345 [00:07<00:04, 4875.65it/s]
65%|███████████████████████████████████████████████████████████████████████████████▌ | 37082/57345 [00:07<00:04, 4876.18it/s]
66%|████████████████████████████████████████████████████████████████████████████████▌ | 37572/57345 [00:07<00:04, 4882.34it/s]
66%|█████████████████████████████████████████████████████████████████████████████████▋ | 38066/57345 [00:07<00:03, 4894.90it/s]
67%|██████████████████████████████████████████████████████████████████████████████████▋ | 38558/57345 [00:08<00:03, 4895.19it/s]
68%|███████████████████████████████████████████████████████████████████████████████████▊ | 39048/57345 [00:08<00:03, 4891.17it/s]
69%|████████████████████████████████████████████████████████████████████████████████████▊ | 39538/57345 [00:08<00:03, 4844.89it/s]
70%|█████████████████████████████████████████████████████████████████████████████████████▊ | 40028/57345 [00:08<00:03, 4856.70it/s]
71%|██████████████████████████████████████████████████████████████████████████████████████▉ | 40514/57345 [00:08<00:03, 4851.50it/s]
71%|███████████████████████████████████████████████████████████████████████████████████████▉ | 41000/57345 [00:08<00:03, 4848.93it/s]
72%|█████████████████████████████████████████████████████████████████████████████████████████ | 41494/57345 [00:08<00:03, 4869.89it/s]
73%|██████████████████████████████████████████████████████████████████████████████████████████ | 41986/57345 [00:08<00:03, 4881.85it/s]
74%|███████████████████████████████████████████████████████████████████████████████████████████ | 42475/57345 [00:08<00:03, 4878.25it/s]
75%|████████████████████████████████████████████████████████████████████████████████████████████▏ | 42963/57345 [00:08<00:02, 4872.99it/s]
76%|█████████████████████████████████████████████████████████████████████████████████████████████▏ | 43451/57345 [00:09<00:02, 4870.52it/s]
77%|██████████████████████████████████████████████████████████████████████████████████████████████▎ | 43945/57345 [00:09<00:02, 4883.09it/s]
77%|███████████████████████████████████████████████████████████████████████████████████████████████▎ | 44434/57345 [00:09<00:02, 4864.39it/s]
78%|████████████████████████████████████████████████████████████████████████████████████████████████▎ | 44921/57345 [00:09<00:02, 4860.42it/s]
79%|█████████████████████████████████████████████████████████████████████████████████████████████████▍ | 45408/57345 [00:09<00:02, 4843.12it/s]
80%|██████████████████████████████████████████████████████████████████████████████████████████████████▍ | 45893/57345 [00:09<00:02, 4839.90it/s]
81%|███████████████████████████████████████████████████████████████████████████████████████████████████▌ | 46389/57345 [00:09<00:02, 4872.25it/s]
82%|████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 46877/57345 [00:09<00:02, 4869.48it/s]
83%|█████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 47364/57345 [00:09<00:02, 4849.09it/s]
83%|██████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 47853/57345 [00:09<00:01, 4856.48it/s]
84%|███████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 48352/57345 [00:10<00:01, 4889.07it/s]
85%|████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 48841/57345 [00:10<00:01, 4860.97it/s]
86%|█████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 49332/57345 [00:10<00:01, 4868.85it/s]
87%|██████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 49819/57345 [00:10<00:01, 4862.68it/s]
88%|███████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 50306/57345 [00:10<00:01, 4778.21it/s]
89%|████████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 50792/57345 [00:10<00:01, 4795.89it/s]
89%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 51280/57345 [00:10<00:01, 4814.92it/s]
90%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████ | 51771/57345 [00:10<00:01, 4837.30it/s]
91%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████ | 52264/57345 [00:10<00:01, 4855.59it/s]
92%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████▏ | 52754/57345 [00:10<00:00, 4863.16it/s]
93%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████▏ | 53245/57345 [00:11<00:00, 4870.45it/s]
94%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████▎ | 53733/57345 [00:11<00:00, 4866.57it/s]
95%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▎ | 54220/57345 [00:11<00:00, 4863.70it/s]
95%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▎ | 54707/57345 [00:11<00:00, 4861.32it/s]
96%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▍ | 55194/57345 [00:11<00:00, 4857.47it/s]
97%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▍ | 55680/57345 [00:11<00:00, 4853.17it/s]
98%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▍ | 56166/57345 [00:11<00:00, 4850.33it/s]
99%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 56658/57345 [00:11<00:00, 4864.33it/s]
100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▌| 57145/57345 [00:11<00:00, 4858.90it/s]
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 57345/57345 [00:11<00:00, 4798.85it/s]
0%| | 0/57345 [00:00<?, ?it/s]
1%|▉ | 452/57345 [00:00<00:12, 4496.34it/s]
2%|██ | 950/57345 [00:00<00:11, 4769.04it/s]
3%|███▏ | 1447/57345 [00:00<00:11, 4849.22it/s]
3%|████▏ | 1932/57345 [00:00<00:11, 4833.47it/s]
4%|█████▎ | 2431/57345 [00:00<00:11, 4880.68it/s]
5%|██████▎ | 2920/57345 [00:00<00:11, 4874.53it/s]
6%|███████▍ | 3422/57345 [00:00<00:10, 4915.59it/s]
7%|████████▍ | 3923/57345 [00:00<00:10, 4939.04it/s]
8%|█████████▌ | 4417/57345 [00:00<00:10, 4934.95it/s]
9%|██████████▌ | 4913/57345 [00:01<00:10, 4935.60it/s]
9%|███████████▋ | 5407/57345 [00:01<00:10, 4928.85it/s]
10%|████████████▊ | 5904/57345 [00:01<00:10, 4936.87it/s]
11%|█████████████▊ | 6398/57345 [00:01<00:10, 4934.08it/s]
12%|██████████████▉ | 6895/57345 [00:01<00:10, 4937.28it/s]
13%|███████████████▉ | 7394/57345 [00:01<00:10, 4948.82it/s]
14%|█████████████████ | 7889/57345 [00:01<00:10, 4828.37it/s]
15%|██████████████████ | 8373/57345 [00:01<00:10, 4825.97it/s]
15%|███████████████████▏ | 8865/57345 [00:01<00:10, 4841.08it/s]
16%|████████████████████▏ | 9356/57345 [00:01<00:09, 4854.04it/s]
17%|█████████████████████▎ | 9847/57345 [00:02<00:09, 4870.65it/s]
18%|██████████████████████▏ | 10335/57345 [00:02<00:09, 4792.39it/s]
19%|███████████████████████▏ | 10825/57345 [00:02<00:09, 4815.99it/s]
20%|████████████████████████▎ | 11336/57345 [00:02<00:09, 4889.78it/s]
21%|█████████████████████████▍ | 11838/57345 [00:02<00:09, 4923.04it/s]
22%|██████████████████████████▍ | 12331/57345 [00:02<00:09, 4906.26it/s]
22%|███████████████████████████▌ | 12822/57345 [00:02<00:09, 4872.59it/s]
23%|████████████████████████████▌ | 13327/57345 [00:02<00:08, 4917.36it/s]
24%|█████████████████████████████▋ | 13825/57345 [00:02<00:08, 4930.83it/s]
25%|██████████████████████████████▋ | 14319/57345 [00:03<00:12, 3460.39it/s]
26%|███████████████████████████████▊ | 14814/57345 [00:03<00:11, 3800.18it/s]
27%|████████████████████████████████▊ | 15310/57345 [00:03<00:10, 4083.48it/s]
28%|█████████████████████████████████▉ | 15811/57345 [00:03<00:09, 4320.52it/s]
28%|██████████████████████████████████▉ | 16299/57345 [00:03<00:09, 4466.82it/s]
29%|███████████████████████████████████▉ | 16771/57345 [00:03<00:09, 4470.12it/s]
30%|████████████████████████████████████▉ | 17237/57345 [00:03<00:08, 4518.67it/s]
31%|██████████████████████████████████████ | 17729/57345 [00:03<00:08, 4619.34it/s]
32%|███████████████████████████████████████ | 18206/57345 [00:03<00:08, 4657.08it/s]
33%|████████████████████████████████████████ | 18685/57345 [00:03<00:08, 4691.35it/s]
33%|█████████████████████████████████████████▏ | 19176/57345 [00:04<00:08, 4748.79it/s]
34%|██████████████████████████████████████████▏ | 19659/57345 [00:04<00:07, 4767.63it/s]
35%|███████████████████████████████████████████▏ | 20144/57345 [00:04<00:07, 4785.84it/s]
36%|████████████████████████████████████████████▏ | 20625/57345 [00:04<00:07, 4788.69it/s]
37%|█████████████████████████████████████████████▎ | 21106/57345 [00:04<00:07, 4775.55it/s]
38%|██████████████████████████████████████████████▎ | 21587/57345 [00:04<00:07, 4779.97it/s]
38%|███████████████████████████████████████████████▎ | 22066/57345 [00:04<00:07, 4733.58it/s]
39%|████████████████████████████████████████████████▎ | 22547/57345 [00:04<00:07, 4751.32it/s]
40%|█████████████████████████████████████████████████▍ | 23036/57345 [00:04<00:07, 4787.01it/s]
41%|██████████████████████████████████████████████████▍ | 23515/57345 [00:04<00:07, 4786.75it/s]
42%|███████████████████████████████████████████████████▍ | 24006/57345 [00:05<00:06, 4817.71it/s]
43%|████████████████████████████████████████████████████▌ | 24488/57345 [00:05<00:06, 4725.28it/s]
44%|█████████████████████████████████████████████████████▌ | 24961/57345 [00:05<00:06, 4721.52it/s]
44%|██████████████████████████████████████████████████████▌ | 25436/57345 [00:05<00:06, 4725.78it/s]
45%|███████████████████████████████████████████████████████▌ | 25930/57345 [00:05<00:06, 4784.01it/s]
46%|████████████████████████████████████████████████████████▋ | 26415/57345 [00:05<00:06, 4799.72it/s]
47%|█████████████████████████████████████████████████████████▋ | 26907/57345 [00:05<00:06, 4829.45it/s]
48%|██████████████████████████████████████████████████████████▊ | 27411/57345 [00:05<00:06, 4886.78it/s]
49%|███████████████████████████████████████████████████████████▊ | 27907/57345 [00:05<00:05, 4908.33it/s]
50%|████████████████████████████████████████████████████████████▉ | 28407/57345 [00:05<00:05, 4929.09it/s]
50%|█████████████████████████████████████████████████████████████▉ | 28904/57345 [00:06<00:05, 4934.29it/s]
51%|███████████████████████████████████████████████████████████████ | 29400/57345 [00:06<00:05, 4932.65it/s]
52%|████████████████████████████████████████████████████████████████▏ | 29901/57345 [00:06<00:05, 4949.12it/s]
53%|█████████████████████████████████████████████████████████████████▏ | 30396/57345 [00:06<00:05, 4943.16it/s]
54%|██████████████████████████████████████████████████████████████████▎ | 30891/57345 [00:06<00:05, 4937.31it/s]
55%|███████████████████████████████████████████████████████████████████▎ | 31393/57345 [00:06<00:05, 4959.24it/s]
56%|████████████████████████████████████████████████████████████████████▍ | 31889/57345 [00:06<00:05, 4955.04it/s]
56%|█████████████████████████████████████████████████████████████████████▍ | 32385/57345 [00:06<00:05, 4950.23it/s]
57%|██████████████████████████████████████████████████████████████████████▌ | 32883/57345 [00:06<00:04, 4954.42it/s]
58%|███████████████████████████████████████████████████████████████████████▌ | 33380/57345 [00:07<00:04, 4953.12it/s]
59%|████████████████████████████████████████████████████████████████████████▋ | 33882/57345 [00:07<00:04, 4965.99it/s]
60%|█████████████████████████████████████████████████████████████████████████▋ | 34382/57345 [00:07<00:04, 4969.42it/s]
61%|██████████████████████████████████████████████████████████████████████████▊ | 34882/57345 [00:07<00:04, 4969.26it/s]
62%|███████████████████████████████████████████████████████████████████████████▉ | 35380/57345 [00:07<00:04, 4958.59it/s]
63%|████████████████████████████████████████████████████████████████████████████▉ | 35876/57345 [00:07<00:04, 4954.62it/s]
63%|██████████████████████████████████████████████████████████████████████████████ | 36379/57345 [00:07<00:04, 4966.66it/s]
64%|███████████████████████████████████████████████████████████████████████████████ | 36877/57345 [00:07<00:04, 4963.55it/s]
65%|████████████████████████████████████████████████████████████████████████████████▏ | 37374/57345 [00:07<00:04, 4958.32it/s]
66%|█████████████████████████████████████████████████████████████████████████████████▏ | 37870/57345 [00:07<00:03, 4937.19it/s]
67%|██████████████████████████████████████████████████████████████████████████████████▎ | 38366/57345 [00:08<00:03, 4937.49it/s]
68%|███████████████████████████████████████████████████████████████████████████████████▎ | 38866/57345 [00:08<00:03, 4951.08it/s]
69%|████████████████████████████████████████████████████████████████████████████████████▍ | 39362/57345 [00:08<00:03, 4947.15it/s]
70%|█████████████████████████████████████████████████████████████████████████████████████▍ | 39858/57345 [00:08<00:03, 4944.64it/s]
70%|██████████████████████████████████████████████████████████████████████████████████████▌ | 40355/57345 [00:08<00:03, 4945.72it/s]
71%|███████████████████████████████████████████████████████████████████████████████████████▌ | 40850/57345 [00:08<00:03, 4923.35it/s]
72%|████████████████████████████████████████████████████████████████████████████████████████▋ | 41343/57345 [00:08<00:03, 4919.61it/s]
73%|█████████████████████████████████████████████████████████████████████████████████████████▋ | 41840/57345 [00:08<00:03, 4929.53it/s]
74%|██████████████████████████████████████████████████████████████████████████████████████████▊ | 42335/57345 [00:08<00:03, 4929.79it/s]
75%|███████████████████████████████████████████████████████████████████████████████████████████▉ | 42849/57345 [00:08<00:02, 4985.71it/s]
76%|████████████████████████████████████████████████████████████████████████████████████████████▉ | 43355/57345 [00:09<00:02, 5000.48it/s]
76%|██████████████████████████████████████████████████████████████████████████████████████████████ | 43861/57345 [00:09<00:02, 5013.21it/s]
77%|███████████████████████████████████████████████████████████████████████████████████████████████▏ | 44363/57345 [00:09<00:02, 4980.11it/s]
78%|████████████████████████████████████████████████████████████████████████████████████████████████▏ | 44862/57345 [00:09<00:02, 4977.18it/s]
79%|█████████████████████████████████████████████████████████████████████████████████████████████████▎ | 45360/57345 [00:09<00:02, 4958.17it/s]
80%|██████████████████████████████████████████████████████████████████████████████████████████████████▎ | 45863/57345 [00:09<00:02, 4970.68it/s]
81%|███████████████████████████████████████████████████████████████████████████████████████████████████▍ | 46361/57345 [00:09<00:02, 4953.22it/s]
82%|████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 46863/57345 [00:09<00:02, 4960.44it/s]
83%|█████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 47360/57345 [00:09<00:02, 4962.26it/s]
83%|██████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 47857/57345 [00:09<00:01, 4829.82it/s]
84%|███████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 48342/57345 [00:10<00:01, 4830.57it/s]
85%|████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 48839/57345 [00:10<00:01, 4864.91it/s]
86%|█████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 49326/57345 [00:10<00:01, 4844.95it/s]
87%|██████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 49812/57345 [00:10<00:01, 4842.93it/s]
88%|███████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 50312/57345 [00:10<00:01, 4884.74it/s]
89%|████████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 50801/57345 [00:10<00:01, 4880.70it/s]
89%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████ | 51303/57345 [00:10<00:01, 4916.58it/s]
90%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████ | 51795/57345 [00:10<00:01, 4899.18it/s]
91%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████▏ | 52295/57345 [00:10<00:01, 4923.40it/s]
92%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████▏ | 52789/57345 [00:10<00:00, 4922.49it/s]
93%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████▎ | 53286/57345 [00:11<00:00, 4932.86it/s]
94%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████▍ | 53796/57345 [00:11<00:00, 4968.85it/s]
95%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▍ | 54293/57345 [00:11<00:00, 4949.35it/s]
96%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 54788/57345 [00:11<00:00, 4944.39it/s]
96%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 55283/57345 [00:11<00:00, 4925.22it/s]
97%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 55776/57345 [00:11<00:00, 4887.40it/s]
98%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 56265/57345 [00:11<00:00, 4880.62it/s]
99%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 56763/57345 [00:11<00:00, 4903.16it/s]
100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊| 57256/57345 [00:11<00:00, 4905.16it/s]
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 57345/57345 [00:11<00:00, 4830.06it/s]
We use the following code to plot the MD, FA and MK estimates from the two data fits:
FA_orig = dki_orig.fa
FA_den = dki_den.fa
MD_orig = dki_orig.md
MD_den = dki_den.md
MK_orig = dki_orig.mk(0, 3)
MK_den = dki_den.mk(0, 3)
fig2, ax = plt.subplots(2, 3, figsize=(10, 6),
subplot_kw={'xticks': [], 'yticks': []})
fig2.subplots_adjust(hspace=0.3, wspace=0.03)
ax.flat[0].imshow(MD_orig[:, :, sli].T, cmap='gray', vmin=0, vmax=2.0e-3,
origin='lower')
ax.flat[0].set_title('MD (DKI)')
ax.flat[1].imshow(FA_orig[:, :, sli].T, cmap='gray', vmin=0, vmax=0.7,
origin='lower')
ax.flat[1].set_title('FA (DKI)')
ax.flat[2].imshow(MK_orig[:, :, sli].T, cmap='gray', vmin=0, vmax=1.5,
origin='lower')
ax.flat[2].set_title('AD (DKI)')
ax.flat[3].imshow(MD_den[:, :, sli].T, cmap='gray', vmin=0, vmax=2.0e-3,
origin='lower')
ax.flat[3].set_title('MD (DKI)')
ax.flat[4].imshow(FA_den[:, :, sli].T, cmap='gray', vmin=0, vmax=0.7,
origin='lower')
ax.flat[4].set_title('FA (DKI)')
ax.flat[5].imshow(MK_den[:, :, sli].T, cmap='gray', vmin=0, vmax=1.5,
origin='lower')
ax.flat[5].set_title('AD (DKI)')
plt.show()
fig2.savefig('denoised_dki.png')
In the above figure, the DKI maps obtained from the original data are shown in the upper panels, while the DKI maps from the denoised data are shown in the lower panels. Substantial improvements in measurement robustness can be visually appreciated, particularly for the FA and MK estimates.
Noise standard deviation estimation using the Marchenko-Pastur PCA algorithm#
As mentioned above, the Marcenko-Pastur PCA algorithm can also be used to
estimate the image’s noise standard deviation (std). The noise std
automatically computed from the mppca function can be returned by
setting the optional input parameter return_sigma to True.
denoised_arr, sigma = mppca(data, patch_radius=2, return_sigma=True)
Let’s plot the noise standard deviation estimate:
fig3 = plt.figure('PCA Noise standard deviation estimation')
plt.imshow(sigma[..., sli].T, cmap='gray', origin='lower')
plt.axis('off')
plt.show()
fig3.savefig('pca_sigma.png')
The above figure shows that the Marchenko-Pastur PCA algorithm computes a 3D spatial varying noise level map. To obtain the mean noise std across all voxels, you can use the following lines of code:
mean_sigma = np.mean(sigma[mask])
print(mean_sigma)
4.6833844
Below we use this mean noise level estimate to compute the nominal SNR of the data (i.e. SNR at b-value=0):
b0 = denoised_arr[..., 0]
mean_signal = np.mean(b0[mask])
snr = mean_signal / mean_sigma
print(snr)
34.56180377412377
References#
[Manjon2013]
(1,2,3)
Manjon JV, Coupe P, Concha L, Buades A, Collins DL “Diffusion Weighted Image Denoising Using Overcomplete Local PCA” (2013) PLoS ONE 8(9): e73021. doi:10.1371/journal.pone.0073021.
[Veraa2016a]
Veraart J, Fieremans E, Novikov DS. 2016. Diffusion MRI noise mapping using random matrix theory. Magnetic Resonance in Medicine. doi: 10.1002/mrm.26059.
[Henri2018]
Henriques, R., 2018. Advanced Methods for Diffusion MRI Data Analysis and their Application to the Healthy Ageing Brain (Doctoral thesis). https://doi.org/10.17863/CAM.29356
[Veraa2016b]
Veraart J, Novikov DS, Christiaens D, Ades-aron B, Sijbers, Fieremans E, 2016. Denoising of Diffusion MRI using random matrix theory. Neuroimage 142:394-406. doi: 10.1016/j.neuroimage.2016.08.016
Total running time of the script: (10 minutes 31.850 seconds)