Note
Click here to download the full example code
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 numpy as np
import matplotlib.pyplot as plt
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('difference, SSD=%.02f' % ssd)
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')
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.e-3)
Perform the registration with equal images.
mapping = sdr.optimize(img_ref.astype(float), img_ref.astype(float))
img_warp = mapping.transform(img_ref, 'linear')
show_images(img_ref, img_warp, 'output-0')
regtools.plot_2d_diffeomorphic_map(mapping, 5, 'map-0.png')
(array([[ 0., 0., 0., ..., 0., 0., 0.],
[ 0., 127., 127., ..., 127., 127., 127.],
[ 0., 127., 127., ..., 127., 127., 127.],
...,
[ 0., 127., 127., ..., 127., 127., 127.],
[ 0., 127., 127., ..., 127., 127., 127.],
[ 0., 127., 127., ..., 127., 127., 127.]], dtype=float32), array([[ 0., 0., 0., ..., 0., 0., 0.],
[ 0., 127., 127., ..., 127., 127., 127.],
[ 0., 127., 127., ..., 127., 127., 127.],
...,
[ 0., 127., 127., ..., 127., 127., 127.],
[ 0., 127., 127., ..., 127., 127., 127.],
[ 0., 127., 127., ..., 127., 127., 127.]], dtype=float32))
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, 'linear')
show_images(img_ref, img_warp, 'output-1')
regtools.plot_2d_diffeomorphic_map(mapping, 5, 'map-1.png')
(array([[ 0. , 0. , 0. , ..., 0. , 0. ,
0. ],
[ 0. , 89.38399 , 81.57997 , ..., 122.46935 , 125.357605,
127. ],
[ 0. , 99.70613 , 127.00001 , ..., 127. , 127. ,
127. ],
...,
[ 0. , 113.511505, 127. , ..., 124.32039 , 110.232216,
127. ],
[ 0. , 127. , 127. , ..., 117.77921 , 126.14425 ,
127. ],
[ 0. , 127. , 127. , ..., 127. , 127. ,
127. ]], dtype=float32), array([[ 0. , 0. , 0. , ..., 0. , 0. ,
0. ],
[ 0. , 127. , 127. , ..., 115.19761 , 76.96389 ,
127. ],
[ 0. , 127. , 126.99999 , ..., 127. , 10.195205,
127. ],
...,
[ 0. , 124.391754, 121.773254, ..., 127. , 127. ,
127. ],
[ 0. , 127. , 127. , ..., 127. , 127. ,
127. ],
[ 0. , 127. , 127. , ..., 127. , 127. ,
127. ]], dtype=float32))
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, 'linear')
show_images(img_ref, img_warp, 'output-2')
regtools.plot_2d_diffeomorphic_map(mapping, 5, 'map-2.png')
(array([[ 0. , 0. , 0. , ..., 0. , 0. ,
0. ],
[ 0. , 107.26481 , 107.818726, ..., 94.22971 , 106.945564,
127. ],
[ 0. , 111.12815 , 127. , ..., 127. , 127. ,
127. ],
...,
[ 0. , 122.529335, 127. , ..., 127. , 127. ,
127. ],
[ 0. , 126.99735 , 127. , ..., 127. , 126.99999 ,
127. ],
[ 0. , 127. , 127. , ..., 127. , 127. ,
127. ]], dtype=float32), array([[ 0. , 0. , 0. , ..., 0. ,
0. , 0. ],
[ 0. , 127. , 127. , ..., 38.59097 ,
127. , 127. ],
[ 0. , 127. , 127. , ..., 34.39428 ,
127. , 127. ],
...,
[ 0. , 124.35901 , 122.08593 , ..., 3.9184668,
70.04836 , 127. ],
[ 0. , 127. , 127. , ..., 65.811745 ,
127. , 127. ],
[ 0. , 127. , 127. , ..., 127. ,
127. , 127. ]], dtype=float32))
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.e-3)
Perform the registration.
mapping = sdr.optimize(img_ref.astype(float), img_in.astype(float))
img_warp = mapping.transform(img_in, 'linear')
show_images(img_ref, img_warp, 'output-3')
regtools.plot_2d_diffeomorphic_map(mapping, 5, 'map-3.png')
(array([[ 0. , 0. , 0. , ..., 0. , 0. ,
0. ],
[ 0. , 113.813194, 112.31409 , ..., 69.539505, 94.62846 ,
127. ],
[ 0. , 115.748116, 127. , ..., 127. , 127. ,
127. ],
...,
[ 0. , 126.96909 , 127. , ..., 127. , 123.10664 ,
127. ],
[ 0. , 126.97683 , 127.00001 , ..., 127. , 127. ,
127. ],
[ 0. , 127. , 127. , ..., 127. , 127. ,
127. ]], dtype=float32), array([[ 0. , 0. , 0. , ..., 0. , 0. ,
0. ],
[ 0. , 126.99999 , 127. , ..., 126.99999 , 56.18689 ,
127. ],
[ 0. , 127. , 127. , ..., 86.738045, 97.76541 ,
127. ],
...,
[ 0. , 126.99135 , 126.97298 , ..., 127. , 31.17303 ,
127. ],
[ 0. , 127. , 127. , ..., 87.87969 , 36.326607,
127. ],
[ 0. , 127. , 127. , ..., 127. , 127. ,
127. ]], dtype=float32))
Total running time of the script: ( 0 minutes 2.320 seconds)