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:

t = time()

denoised_arr = mppca(data, patch_radius=2)

print("Time taken for local MP-PCA ", -t + time())
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')
denoise mppca

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')
denoise mppca

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')
denoise mppca

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)

Gallery generated by Sphinx-Gallery