Using MATLAB to calculate offset between successive images

Question

I'm taking images using a tunneling microscope. However, the scope is drifting between successive images. I'm trying to use MatLab to calculate the offset between images. The code below calculates in seconds for small images (e.g. 64x64 pixels), but takes >2 hrs to handle the 512x512 pixel images I'm dealing with. Do you have any suggestions for speeding up this code? Or do you know of better ways to track images in MatLab? Thanks for your help!

%Test templates
template = .5*ones(32);
template(25:32,:) = 0;
template(:,25:64) = 0;
data_A = template;
close all
imshow(data_A);
template(9:32,41:64) = .5;
template(:,1:24) = 0;
data_B = template;
figure, imshow(data_B);

tic

[m n] = size(data_B);
z = [];

% Loop over all possible displacements
for x = -n:n

for y = -m:m

    paddata_B = data_B;
    ax = abs(x);
    zerocols = zeros(m,ax);

    if x > 0
        paddata_B(:,1:ax) = [];
        paddata_B = [paddata_B zerocols];

    else
        paddata_B(:,(n-ax+1):n) = [];
        paddata_B = [zerocols paddata_B];

    end

    ay = abs(y);
    zerorows = zeros(ay,n);


    if y < 0
        paddata_B(1:ay,:) = [];
        paddata_B = vertcat(paddata_B, zerorows);

    else
        paddata_B((m-ay+1):m,:) = [];
        paddata_B = vertcat(zerorows, paddata_B);

    end

% Full matrix sum after array multiplication
C = paddata_B.*data_A;        
matsum = sum(sum(C));

% Populate array of matrix sums for each displacement    
z(x+n+1, y+m+1) = matsum;

end
end

toc

% Plot matrix sums
figure, surf(z), shading flat

% Find maximum value of z matrix
[max_z, imax] = max(abs(z(:)));
[xpeak, ypeak] = ind2sub(size(z),imax(1))

% Calculate displacement in pixels
corr_offset = [(xpeak-n-1) (ypeak-m-1)];
xoffset = corr_offset(1)
yoffset = corr_offset(2)  

Answer 1

What you're calculating is known as the cross-correlation of the two images. You can calculate the cross-correlation of all offsets at once using Discrete Fourier Transforms (DFT or FFT). So try something like

z = ifft2( fft2(dataA) .* fft2(dataB).' );

If you pad with zeros in the Fourier domain, you can even use this sort of math to get offsets in fractions of a pixel, and apply offsets of fractions of a pixel to an image.

Answer 2

A typical approach to this kind of problem is to use the fact that it works quickly for small images to your advantage. When you have large images, decimate them to make small images. Register the small images quickly and use the computed offset as your initial value for the next iteration. In the next iteration, you don't decimate the images as much, but you're starting with a good initial estimate of the offset so you can constrain your search for solutions to a small neighborhood near your initial estimate.

Although not written with tunneling microscopes in mind, a review paper that may be of some assistance is: "Mutual Information-Based Registration of Medical Images: A Survey" by Pluim, Maintz, and Viergever published in IEEE Transactions on Medical Imaging, Vol. 22, No. 8, p. 986.