-
Notifications
You must be signed in to change notification settings - Fork 429
/
contextual_enhancement.py
285 lines (222 loc) · 10.3 KB
/
contextual_enhancement.py
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
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
# -*- coding: utf-8 -*-
"""
==============================================
Crossing-preserving contextual enhancement
==============================================
This demo presents an example of crossing-preserving contextual enhancement of
FOD/ODF fields [Meesters2016]_, implementing the contextual PDE framework
of [Portegies2015a]_ for processing HARDI data. The aim is to enhance the
alignment of elongated structures in the data such that crossing/junctions are
maintained while reducing noise and small incoherent structures. This is
achieved via a hypo-elliptic 2nd order PDE in the domain of coupled positions
and orientations :math:`\mathbb{R}^3 \rtimes S^2`. This domain carries a
non-flat geometrical differential structure that allows including a notion of
alignment between neighboring points.
Let :math:`({\bf y},{\bf n}) \in \mathbb{R}^3\rtimes S^2` where
:math:`{\bf y} \in \mathbb{R}^{3}` denotes the spatial part, and
:math:`{\bf n} \in S^2` the angular part.
Let :math:`W:\mathbb{R}^3\rtimes S^2\times \mathbb{R}^{+} \to \mathbb{R}` be
the function representing the evolution of FOD/ODF field. Then, the contextual
PDE with evolution time :math:`t\geq 0` is given by:
.. math::
\begin{cases}
\frac{\partial}{\partial t} W({\bf y},{\bf n},t) &= ((D^{33}({\bf n} \cdot
\nabla)^2 + D^{44} \Delta_{S^2})W)({\bf y},{\bf n},t)
\\ W({\bf y},{\bf n},0) &= U({\bf y},{\bf n})
\end{cases},
where:
* :math:`D^{33}>0` is the coefficient for the spatial smoothing (which goes only in the direction of :math:`n`);
* :math:`D^{44}>0` is the coefficient for the angular smoothing (here :math:`\Delta_{S^2}` denotes the Laplace-Beltrami operator on the sphere :math:`S^2`);
* :math:`U:\mathbb{R}^3\rtimes S^2 \to \mathbb{R}` is the initial condition given by the noisy FOD/ODF’s field.
This equation is solved via a shift-twist convolution (denoted by :math:`\ast_{\mathbb{R}^3\rtimes S^2}`) with its corresponding kernel :math:`P_t:\mathbb{R}^3\rtimes S^2 \to \mathbb{R}^+`:
.. math::
W({\bf y},{\bf n},t) = (P_t \ast_{\mathbb{R}^3 \rtimes S^2} U)({\bf y},{\bf n})
= \int_{\mathbb{R}^3} \int_{S^2} P_t (R^T_{{\bf n}^\prime}({\bf y}-{\bf y}^\prime),
R^T_{{\bf n}^\prime} {\bf n} ) U({\bf y}^\prime, {\bf n}^\prime)
Here, :math:`R_{\bf n}` is any 3D rotation that maps the vector :math:`(0,0,1)`
onto :math:`{\bf n}`.
Note that the shift-twist convolution differs from a Euclidean convolution and
takes into account the non-flat structure of the space :math:`\mathbb{R}^3\rtimes S^2`.
The kernel :math:`P_t` has a stochastic interpretation [DuitsAndFranken2011]_.
It can be seen as the limiting distribution obtained by accumulating random
walks of particles in the position/orientation domain, where in each step the
particles can (randomly) move forward/backward along their current orientation,
and (randomly) change their orientation. This is an extension to the 3D case of
the process for contour enhancement of 2D images.
.. figure:: _static/stochastic_process.jpg
:scale: 60 %
:align: center
The random motion of particles (a) and its corresponding probability map
(b) in 2D. The 3D kernel is shown on the right. Adapted from
[Portegies2015a]_.
In practice, as the exact analytical formulas for the kernel :math:`P_t`
are unknown, we use the approximation given in [Portegies2015b]_.
"""
"""
The enhancement is evaluated on the Stanford HARDI dataset
(150 orientations, b=2000 $s/mm^2$) where Rician noise is added. Constrained
spherical deconvolution is used to model the fiber orientations.
"""
import numpy as np
from dipy.data import fetch_stanford_hardi, read_stanford_hardi
from dipy.sims.voxel import add_noise
# Read data
fetch_stanford_hardi()
img, gtab = read_stanford_hardi()
data = img.get_data()
# Add Rician noise
from dipy.segment.mask import median_otsu
b0_slice = data[:, :, :, 1]
b0_mask, mask = median_otsu(b0_slice)
np.random.seed(1)
data_noisy = add_noise(data, 10.0, np.mean(b0_slice[mask]), noise_type='rician')
# Select a small part of it.
padding = 3 # Include a larger region to avoid boundary effects
data_small = data[25-padding:40+padding, 65-padding:80+padding, 35:42]
data_noisy_small = data_noisy[25-padding:40+padding, 65-padding:80+padding, 35:42]
"""
Enables/disables interactive visualization
"""
interactive = False
"""
Fit an initial model to the data, in this case Constrained Spherical
Deconvolution is used.
"""
# Perform CSD on the original data
from dipy.reconst.csdeconv import auto_response
from dipy.reconst.csdeconv import ConstrainedSphericalDeconvModel
response, ratio = auto_response(gtab, data, roi_radius=10, fa_thr=0.7)
csd_model_orig = ConstrainedSphericalDeconvModel(gtab, response)
csd_fit_orig = csd_model_orig.fit(data_small)
csd_shm_orig = csd_fit_orig.shm_coeff
# Perform CSD on the original data + noise
response, ratio = auto_response(gtab, data_noisy, roi_radius=10, fa_thr=0.7)
csd_model_noisy = ConstrainedSphericalDeconvModel(gtab, response)
csd_fit_noisy = csd_model_noisy.fit(data_noisy_small)
csd_shm_noisy = csd_fit_noisy.shm_coeff
"""
Inspired by [Rodrigues2010]_, a lookup-table is created, containing
rotated versions of the kernel :math:`P_t` sampled over a discrete set of
orientations. In order to ensure rotationally invariant processing, the discrete
orientations are required to be equally distributed over a sphere. By default,
a sphere with 100 directions is used.
"""
from dipy.denoise.enhancement_kernel import EnhancementKernel
from dipy.denoise.shift_twist_convolution import convolve
# Create lookup table
D33 = 1.0
D44 = 0.02
t = 1
k = EnhancementKernel(D33, D44, t)
"""
Visualize the kernel
"""
from dipy.viz import window, actor
from dipy.data import get_sphere
from dipy.reconst.shm import sf_to_sh, sh_to_sf
ren = window.Renderer()
# convolve kernel with delta spike
spike = np.zeros((7, 7, 7, k.get_orientations().shape[0]), dtype=np.float64)
spike[3, 3, 3, 0] = 1
spike_shm_conv = convolve(sf_to_sh(spike, k.get_sphere(), sh_order=8), k,
sh_order=8, test_mode=True)
sphere = get_sphere('symmetric724')
spike_sf_conv = sh_to_sf(spike_shm_conv, sphere, sh_order=8)
model_kernel = actor.odf_slicer(spike_sf_conv * 6,
sphere=sphere,
norm=False,
scale=0.4)
model_kernel.display(x=3)
ren.add(model_kernel)
ren.set_camera(position=(30, 0, 0), focal_point=(0, 0, 0), view_up=(0, 0, 1))
window.record(ren, out_path='kernel.png', size=(900, 900))
if interactive:
window.show(ren)
"""
.. figure:: kernel.png
:align: center
Visualization of the contour enhancement kernel.
"""
"""
Shift-twist convolution is applied on the noisy data
"""
# Perform convolution
csd_shm_enh = convolve(csd_shm_noisy, k, sh_order=8)
"""
The Sharpening Deconvolution Transform is applied to sharpen the ODF field.
"""
# Sharpen via the Sharpening Deconvolution Transform
from dipy.reconst.csdeconv import odf_sh_to_sharp
csd_shm_enh_sharp = odf_sh_to_sharp(csd_shm_enh, sphere, sh_order=8, lambda_=0.1)
# Convert raw and enhanced data to discrete form
csd_sf_orig = sh_to_sf(csd_shm_orig, sphere, sh_order=8)
csd_sf_noisy = sh_to_sf(csd_shm_noisy, sphere, sh_order=8)
csd_sf_enh = sh_to_sf(csd_shm_enh, sphere, sh_order=8)
csd_sf_enh_sharp = sh_to_sf(csd_shm_enh_sharp, sphere, sh_order=8)
# Normalize the sharpened ODFs
csd_sf_enh_sharp = csd_sf_enh_sharp * np.amax(csd_sf_orig)/np.amax(csd_sf_enh_sharp) * 1.25
"""
The end results are visualized. It can be observed that the end result after
diffusion and sharpening is closer to the original noiseless dataset.
"""
ren = window.Renderer()
# original ODF field
fodf_spheres_org = actor.odf_slicer(csd_sf_orig,
sphere=sphere,
scale=0.4,
norm=False)
fodf_spheres_org.display(z=3)
fodf_spheres_org.SetPosition(0, 25, 0)
ren.add(fodf_spheres_org)
# ODF field with added noise
fodf_spheres = actor.odf_slicer(csd_sf_noisy,
sphere=sphere,
scale=0.4,
norm=False,)
fodf_spheres.SetPosition(0, 0, 0)
ren.add(fodf_spheres)
# Enhancement of noisy ODF field
fodf_spheres_enh = actor.odf_slicer(csd_sf_enh,
sphere=sphere,
scale=0.4,
norm=False)
fodf_spheres_enh.SetPosition(25, 0, 0)
ren.add(fodf_spheres_enh)
# Additional sharpening
fodf_spheres_enh_sharp = actor.odf_slicer(csd_sf_enh_sharp,
sphere=sphere,
scale=0.4,
norm=False)
fodf_spheres_enh_sharp.SetPosition(25, 25, 0)
ren.add(fodf_spheres_enh_sharp)
window.record(ren, out_path='enhancements.png', size=(900, 900))
if interactive:
window.show(ren)
"""
.. figure:: enhancements.png
:align: center
The results after enhancements. Top-left: original noiseless data.
Bottom-left: original data with added Rician noise (SNR=10). Bottom-right:
After enhancement of noisy data. Top-right: After enhancement and sharpening
of noisy data.
References
----------
.. [Meesters2016] S. Meesters, G. Sanguinetti, E. Garyfallidis, J. Portegies,
R. Duits. (2016) Fast implementations of contextual PDE’s for HARDI data
processing in DIPY. ISMRM 2016 conference.
.. [Portegies2015a] J. Portegies, R. Fick, G. Sanguinetti, S. Meesters,
G.Girard, and R. Duits. (2015) Improving Fiber Alignment in HARDI by
Combining Contextual PDE flow with Constrained Spherical Deconvolution.
PLoS One.
.. [Portegies2015b] J. Portegies, G. Sanguinetti, S. Meesters, and R. Duits.
(2015) New Approximation of a Scale Space Kernel on SE(3) and Applications
in Neuroimaging. Fifth International Conference on Scale Space and
Variational Methods in Computer Vision.
.. [DuitsAndFranken2011] R. Duits and E. Franken (2011) Left-invariant
diffusions on the space of positions and orientations and their application
to crossing-preserving smoothing of HARDI images. International Journal of
Computer Vision, 92:231-264.
.. [Rodrigues2010] P. Rodrigues, R. Duits, B. Romeny, A. Vilanova (2010).
Accelerated Diffusion Operators for Enhancing DW-MRI. Eurographics Workshop
on Visual Computing for Biology and Medicine. The Eurographics Association.
"""