Image Registration Visualization Example

From MATLAB Techniques for Image Processing by Steve Eddins.

Contents

This example shows how to use transparency to visualize a registered image. Here are two datasets from Massachusetts. The first is an orthorectified aerial photograph. The second is rasterized vector data of the same region, showing roads and water features.

f_photo = imread('aerial-photo-cropped-r.tif');
imshow(f_photo)
title('Orthorectified aerial photo')
xlabel('Image courtesy of Mass GIS')
Warning: Image is too big to fit on screen; displaying at 33% 
f_roads = imread('vector-gis-data.tif');
imshow(f_roads)
title('Roads and water features')
xlabel('Image courtesy of Mass GIS')
Warning: Image is too big to fit on screen; displaying at 67% 

Using an interactive tool, I selected pairs of matching control points in the two images. I saved these points in a MAT-file to use for this example.

load registration_control_points

Given pairs of point correspondences, compute an affine transform that maps vector GIS image to aerial photo.

tform = fitgeotrans(moving_points,fixed_points,'affine');

Use that affine transform with imwarp to register the f_roads image with the f_photo image. The second output, ref, is the spatial reference for f_roads_registered.

[f_roads_registered,ref] = imwarp(f_roads,tform);
subplot(1,2,1)
imshow(f_roads)
title('Roads/water image')

subplot(1,2,2)
imshow(f_roads_registered)
title('Roads/water image registered')

Superimpose transformed roads/water image over aerial photo.

clf
imshow(f_photo)
hold on
h = imshow(f_roads_registered,ref);
hold off
axis auto
Warning: Image is too big to fit on screen; displaying at 33% 

Make the superimposed image transparent.

set(h,'AlphaData',0.4)

Why does the pond appear to be cut off in the overlay? Because the original roads/water image has been rotated to match the aerial photo. We're actually seeing a lot of fill pixels. Let's use a different fill value to make it more obvious.

[f_roads_registered,ref] = imwarp(f_roads,tform,'FillValues',255);
imshow(f_photo)
hold on
h = imshow(f_roads_registered,ref);
set(h,'AlphaData',0.4)
hold off
axis auto

The "Transform an All-Ones Matrix Transparency Mask Trick"

a = 0.4*ones(size(f_roads));
a = imwarp(a,tform);
set(h,'AlphaData',a)