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
|
"""
=======================================================
Diffeomorphic Registration with binary and fuzzy images
=======================================================
This example demonstrates registration of a binary and a fuzzy image.
This could be seen as aligning a fuzzy (sensed) image to a binary
(e.g., template) image.
"""
import matplotlib.pyplot as plt
import numpy as np
from skimage import draw, filters
from dipy.align.imwarp import SymmetricDiffeomorphicRegistration
from dipy.align.metrics import SSDMetric
from dipy.viz import regtools
###############################################################################
# Let's generate a sample template image as the combination of three ellipses.
# We will generate the fuzzy (sensed) version of the image by smoothing
# the reference image.
def draw_ellipse(img, center, axis):
rr, cc = draw.ellipse(center[0], center[1], axis[0], axis[1], shape=img.shape)
img[rr, cc] = 1
return img
img_ref = np.zeros((64, 64))
img_ref = draw_ellipse(img_ref, (25, 15), (10, 5))
img_ref = draw_ellipse(img_ref, (20, 45), (15, 10))
img_ref = draw_ellipse(img_ref, (50, 40), (7, 15))
img_in = filters.gaussian(img_ref, sigma=3)
###############################################################################
# Let's define a small visualization function.
def show_images(img_ref, img_warp, fig_name):
fig, axarr = plt.subplots(ncols=2, figsize=(12, 5))
axarr[0].set_title("warped image & reference contour")
axarr[0].imshow(img_warp)
axarr[0].contour(img_ref, colors="r")
ssd = np.sum((img_warp - img_ref) ** 2)
axarr[1].set_title(f"difference, SSD={ssd:.02f}")
im = axarr[1].imshow(img_warp - img_ref)
plt.colorbar(im)
fig.tight_layout()
fig.savefig(fig_name + ".png")
show_images(img_ref, img_in, "input")
###############################################################################
# .. rst-class:: centered small fst-italic fw-semibold
#
# Input images before alignment.
#
#
#
# Let's use the general Registration function with some naive parameters,
# such as set `step_length` as 1 assuming maximal step 1 pixel and a
# reasonably small number of iterations since the deformation with already
# aligned images should be minimal.
sdr = SymmetricDiffeomorphicRegistration(
metric=SSDMetric(img_ref.ndim),
step_length=1.0,
level_iters=[50, 100],
inv_iter=50,
ss_sigma_factor=0.1,
opt_tol=1.0e-3,
)
###############################################################################
# Perform the registration with equal images.
mapping = sdr.optimize(img_ref.astype(float), img_ref.astype(float))
img_warp = mapping.transform(img_ref, interpolation="linear")
show_images(img_ref, img_warp, "output-0")
regtools.plot_2d_diffeomorphic_map(mapping, delta=5, fname="map-0.png")
###############################################################################
# .. rst-class:: centered small fst-italic fw-semibold
#
# Registration results for default parameters and equal images.
#
#
#
# Perform the registration with binary and fuzzy images.
mapping = sdr.optimize(img_ref.astype(float), img_in.astype(float))
img_warp = mapping.transform(img_in, interpolation="linear")
show_images(img_ref, img_warp, "output-1")
regtools.plot_2d_diffeomorphic_map(mapping, delta=5, fname="map-1.png")
###############################################################################
# .. rst-class:: centered small fst-italic fw-semibold
#
# Registration results for a naive parameter configuration.
#
#
#
# Note, we are still using a multi-scale approach which makes `step_length`
# in the upper level multiplicatively larger.
# What happens if we set `step_length` to a rather small value?
sdr.step_length = 0.1
###############################################################################
# Perform the registration and examine the output.
mapping = sdr.optimize(img_ref.astype(float), img_in.astype(float))
img_warp = mapping.transform(img_in, interpolation="linear")
show_images(img_ref, img_warp, "output-2")
regtools.plot_2d_diffeomorphic_map(mapping, delta=5, fname="map-2.png")
###############################################################################
# .. rst-class:: centered small fst-italic fw-semibold
#
# Registration results for decreased step size.
#
#
#
# An alternative scenario is to use just a single-scale level.
# Even though the warped image may look fine, the estimated deformations show
# that it is off the mark.
sdr = SymmetricDiffeomorphicRegistration(
metric=SSDMetric(img_ref.ndim),
step_length=1.0,
level_iters=[100],
inv_iter=50,
ss_sigma_factor=0.1,
opt_tol=1.0e-3,
)
###############################################################################
# Perform the registration.
mapping = sdr.optimize(img_ref.astype(float), img_in.astype(float))
img_warp = mapping.transform(img_in, interpolation="linear")
show_images(img_ref, img_warp, "output-3")
regtools.plot_2d_diffeomorphic_map(mapping, delta=5, fname="map-3.png")
###############################################################################
# .. rst-class:: centered small fst-italic fw-semibold
#
# Registration results for single level.
|
