1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
|
"""
==============================================================
Denoise images using Adaptive Soft Coefficient Matching (ASCM)
==============================================================
The adaptive soft coefficient matching (ASCM) as described in
:footcite:p:`Coupe2012` is an improved extension of non-local means (NLMEANS)
denoising. ASCM gives a better denoised images from two standard non-local means
denoised versions of the original data with different degrees sharpness. Here,
one denoised input is more "smooth" than the other (the easiest way to achieve
this denoising is use :ref:`non_local_means<sphx_glr_examples_built_preprocessing_denoise_nlmeans.py>`
with two different patch radii).
ASCM involves these basic steps
* Computes wavelet decomposition of the noisy as well as denoised inputs
* Combines the wavelets for the output image in a way that it takes it's
smoothness (low frequency components) from the input with larger smoothing,
and the sharp features (high frequency components) from the input with
less smoothing.
This way ASCM gives us a well denoised output while preserving the sharpness
of the image features.
Let us load the necessary modules
""" # noqa: E501
from time import time
import matplotlib.pyplot as plt
import numpy as np
from dipy.core.gradients import gradient_table
from dipy.data import get_fnames
from dipy.denoise.adaptive_soft_matching import adaptive_soft_matching
from dipy.denoise.noise_estimate import estimate_sigma
from dipy.denoise.non_local_means import non_local_means
from dipy.io.gradients import read_bvals_bvecs
from dipy.io.image import load_nifti, save_nifti
###############################################################################
# Choose one of the data from the datasets in dipy_
dwi_fname, dwi_bval_fname, dwi_bvec_fname = get_fnames(name="sherbrooke_3shell")
data, affine = load_nifti(dwi_fname)
bvals, bvecs = read_bvals_bvecs(dwi_bval_fname, dwi_bvec_fname)
gtab = gradient_table(bvals, bvecs=bvecs)
mask = data[..., 0] > 80
data = data[..., 1]
print("vol size", data.shape)
t = time()
###############################################################################
# In order to generate the two pre-denoised versions of the data we will use
# the :ref:`non_local_means denoining<sphx_glr_examples_built_preprocessing_denoise_nlmeans.py>` # noqa: E501
# For ``non_local_means`` first we need to estimate the standard deviation of
# the noise. We use N=4 since the Sherbrooke dataset was acquired on a
# 1.5T Siemens scanner with a 4 array head coil.
sigma = estimate_sigma(data, N=4)
###############################################################################
# For the denoised version of the original data which preserves sharper
# features, we perform non-local means with smaller patch size.
den_small = non_local_means(
data, sigma=sigma, mask=mask, patch_radius=1, block_radius=1, rician=True
)
###############################################################################
# For the denoised version of the original data that implies more smoothing, we
# perform non-local means with larger patch size.
den_large = non_local_means(
data, sigma=sigma, mask=mask, patch_radius=2, block_radius=1, rician=True
)
###############################################################################
# Now we perform the adaptive soft coefficient matching. Empirically we set the
# adaptive parameter in ascm to be the average of the local noise variance,
# in this case the sigma itself.
den_final = adaptive_soft_matching(data, den_small, den_large, sigma[0])
print("total time", time() - t)
###############################################################################
# To access the quality of this denoising procedure, we plot an axial slice
# of the original data, it's denoised output and residuals.
axial_middle = data.shape[2] // 2
original = data[:, :, axial_middle].T
final_output = den_final[:, :, axial_middle].T
difference = np.abs(final_output.astype(np.float64) - original.astype(np.float64))
difference[~mask[:, :, axial_middle].T] = 0
fig, ax = plt.subplots(1, 3)
ax[0].imshow(original, cmap="gray", origin="lower")
ax[0].set_title("Original")
ax[1].imshow(final_output, cmap="gray", origin="lower")
ax[1].set_title("ASCM output")
ax[2].imshow(difference, cmap="gray", origin="lower")
ax[2].set_title("Residual")
for i in range(3):
ax[i].set_axis_off()
plt.savefig("denoised_ascm.png", bbox_inches="tight")
###############################################################################
# .. rst-class:: centered small fst-italic fw-semibold
#
# Showing the axial slice without (left) and with (middle) ASCM denoising.
#
#
# From the above figure we can see that the residual is really uniform in
# nature which dictates that ASCM denoises the data while preserving the
# sharpness of the features. Now, we are Saving the entire denoised output in
# ``denoised_ascm.nii.gz`` file.
save_nifti("denoised_ascm.nii.gz", den_final, affine)
###############################################################################
# For comparison propose we also plot the outputs of the ``non_local_means``
# (both with the larger as well as with the smaller patch radius) with the ASCM
# output.
fig, ax = plt.subplots(1, 4)
ax[0].imshow(original, cmap="gray", origin="lower")
ax[0].set_title("Original")
ax[1].imshow(
den_small[..., axial_middle].T, cmap="gray", origin="lower", interpolation="none"
)
ax[1].set_title("NLMEANS small")
ax[2].imshow(
den_large[..., axial_middle].T, cmap="gray", origin="lower", interpolation="none"
)
ax[2].set_title("NLMEANS large")
ax[3].imshow(final_output, cmap="gray", origin="lower", interpolation="none")
ax[3].set_title("ASCM ")
for i in range(4):
ax[i].set_axis_off()
plt.savefig("ascm_comparison.png", bbox_inches="tight")
###############################################################################
# .. rst-class:: centered small fst-italic fw-semibold
#
# Comparing outputs of the NLMEANS and ASCM.
#
#
# From the above figure, we can observe that the information of two
# pre-denoised versions of the raw data, ASCM outperforms standard non-local
# means in suppressing noise and preserving feature sharpness.
#
# References
# ----------
#
# .. footbibliography::
#
|
