How to perform intensity image registration with manual initialisation.¶
The goal here is not to explain how to manually define landmarks (for that look here) but to show how to use pointmatching and blockmatching algorithms to perform intensity image registration with manual initialisation.
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from timagetk.io import imread
from timagetk.io.image import _image_from_url
from timagetk.algorithms.trsf import apply_trsf
from timagetk.algorithms.trsf import compose_trsf
from timagetk.algorithms.blockmatching import blockmatching
from timagetk.algorithms.pointmatching import pointmatching
/builds/J-gEBwyb/0/mosaic/timagetk/src/timagetk/components/labelled_image.py:31: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
from tqdm.autonotebook import tqdm
Displaying the shared data¶
The landmarks have been manually defined on a top-view projection map, so we start by computing the intensity image projections (and associated altitude maps):
from timagetk.algorithms.reconstruction import image_surface_projection
ref_proj, ref_alti = image_surface_projection(ref_img, threshold=np.percentile(ref_img, 60), orientation=1)
flo_proj, flo_alti = image_surface_projection(flo_img, threshold=np.percentile(flo_img, 60), orientation=1)
import matplotlib.gridspec as gridspec
from matplotlib import colors
fig = plt.figure(figsize=(15, 10), constrained_layout=True)
# Create two panels
gs = gridspec.GridSpec(ncols=2, nrows=1, figure=fig)
ref_ax = fig.add_subplot(gs[0, 0])
flo_ax = fig.add_subplot(gs[0, 1])
# Add the projection map as images:
mpl_extent = {}
for idx, image in zip(['reference', 'floating'], [ref_proj, flo_proj]):
vxs = image.get_voxelsize()
im_extent = image.get_extent()
mpl_extent[idx] = (-vxs[0]/2, im_extent[0]+vxs[0]/2, im_extent[1]+vxs[1]/2, -vxs[1]/2)
images = []
images.append(ref_ax.imshow(ref_proj, cmap='gray', extent=mpl_extent['reference']))
images.append(flo_ax.imshow(flo_proj, cmap='gray', extent=mpl_extent['floating']))
ref_ax.xaxis.tick_top()
flo_ax.xaxis.tick_top()
# Add the landmarks points:
px, py, c = np.array([(p[0], p[1], n) for n, p in enumerate(ref_pts)]).T
ref_ax.scatter(px, py, c=c, cmap='Set3', ec='k', s=100)
px, py, c = np.array([(p[0], p[1], n) for n, p in enumerate(flo_pts)]).T
flo_ax.scatter(px, py, c=c, cmap='Set3', ec='k', s=100)
# Find the min and max of all colors for use in setting the color scale.
vmin = min(image.get_array().min() for image in images)
vmax = max(image.get_array().max() for image in images)
norm = colors.Normalize(vmin=vmin, vmax=vmax)
for im in images:
im.set_norm(norm)
# Add a colorbar at the bottom
_ = fig.colorbar(images[0], ax=[ref_ax, flo_ax], orientation='horizontal', fraction=.1)
Initial transformation from manually defined landmarks¶
rigid_init_trsf = pointmatching(flo_pts, ref_pts, template_img=ref_img, method='rigid', real=True)
registered_image = apply_trsf(flo_img, trsf=rigid_init_trsf, template_img=ref_img, interpolation='linear')
reg_proj, reg_alti = image_surface_projection(registered_image, threshold=np.percentile(registered_image, 75), orientation=1)
from timagetk.visu.mplt import grayscale_imshow
_ = grayscale_imshow([ref_proj, reg_proj])
0%| | 0/2 [00:00<?, ?image/s]
100%|██████████| 2/2 [00:00<00:00, 8371.86image/s]
Refining the registration¶
rigid_trsf = blockmatching(flo_img, ref_img, template_img=ref_img, method='rigid', left_trsf=rigid_init_trsf)
INFO [timagetk.algorithms.blockmatching] RIGID registration with a 'left' initialisation matrix
- processing level # 4
- Iteration # 1 Level # 4 Size 32x 32x 32
- Iteration # 2 Level # 4 Size 32x 32x 32
- Iteration # 3 Level # 4 Size 32x 32x 32
- Iteration # 4 Level # 4 Size 32x 32x 32
- Iteration # 5 Level # 4 Size 32x 32x 32
- Iteration # 6 Level # 4 Size 32x 32x 32
- Iteration # 7 Level # 4 Size 32x 32x 32
- Iteration # 8 Level # 4 Size 32x 32x 32
- Iteration # 9 Level # 4 Size 32x 32x 32
- Iteration # 10 Level # 4 Size 32x 32x 32
- processing level # 3
- Iteration # 1 Level # 3 Size 64x 64x 59
- Iteration # 2 Level # 3 Size 64x 64x 59
- Iteration # 3 Level # 3 Size 64x 64x 59
- Iteration # 4 Level # 3 Size 64x 64x 59
- Iteration # 5 Level # 3 Size 64x 64x 59
- Iteration # 6 Level # 3 Size 64x 64x 59
- Iteration # 7 Level # 3 Size 64x 64x 59
- Iteration # 8 Level # 3 Size 64x 64x 59
- Iteration # 9 Level # 3 Size 64x 64x 59
- Iteration # 10 Level # 3 Size 64x 64x 59
- processing level # 2
- Iteration # 1 Level # 2 Size 128x 128x 59
- Iteration # 2 Level # 2 Size 128x 128x 59
- Iteration # 3 Level # 2 Size 128x 128x 59
- Iteration # 4 Level # 2 Size 128x 128x 59
- Iteration # 5 Level # 2 Size 128x 128x 59
- Iteration # 6 Level # 2 Size 128x 128x 59
- Iteration # 7 Level # 2 Size 128x 128x 59
- Iteration # 8 Level # 2 Size 128x 128x 59
- Iteration # 9 Level # 2 Size 128x 128x 59
- Iteration # 10 Level # 2 Size 128x 128x 59
rigid_trsf = compose_trsf([rigid_init_trsf, rigid_trsf], template_img=ref_img)
registered_image = apply_trsf(flo_img, trsf=rigid_trsf, template_img=ref_img, interpolation='linear')
reg_proj, reg_alti = image_surface_projection(registered_image, threshold=np.percentile(registered_image, 75), orientation=1)
from timagetk.visu.mplt import grayscale_imshow
_ = grayscale_imshow([ref_proj, reg_proj])
0%| | 0/2 [00:00<?, ?image/s]
100%|██████████| 2/2 [00:00<00:00, 10343.54image/s]
