Using MATLAB to calculate offset between successive images

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)  

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.

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.

below link will help you find transformation between 2 images and correct/recover the distorted (in your case, image with offset)

http://in.mathworks.com/help/vision/ref/estimategeometrictransform.html

index_pairs = matchFeatures(featuresOriginal,featuresDistorted, 'unique', true);
matchedPtsOriginal  = validPtsOriginal(index_pairs(:,1));
matchedPtsDistorted = validPtsDistorted(index_pairs(:,2));
[tform,inlierPtsDistorted,inlierPtsOriginal] = estimateGeometricTransform(matchedPtsDistorted,matchedPtsOriginal,'similarity');
figure; showMatchedFeatures(original,distorted,inlierPtsOriginal,inlierPtsDistorted);

The inlierPtsDistored, inlierPtsOriginal have attributes called locations. These are nothing but matching locations of one image on another. I think from that point it is very easy to calculate offset.

The function below was my attempt to compute the cross-correlation of the two images manually. Something's not quite right though. Will look at it again this weekend if I have time. You can call the function with something like:

>> oldImage = rand(64);
>> newImage = circshift(oldImage, floor(64/2)*[1 1]);
>> offset = detectOffset(oldImage, newImage, 10)

offset =

    32    -1

function offset = detectOffset(oldImage, newImage, margin)

    if size(oldImage) ~= size(newImage)
        offset = [];
        error('Test images must be the same size.');
    end

    [imageHeight, imageWidth] = size(oldImage);

    corr = zeros(2 * imageHeight - 1, 2 * imageWidth - 1);

    for yIndex = [1:2*imageHeight-1; ...
                  imageHeight:-1:1 ones(1, imageHeight-1); ...
                  imageHeight*ones(1, imageHeight) imageHeight-1:-1:1];
        oldImage = circshift(oldImage, [1 0]);
        for xIndex = [1:2*imageWidth-1; ...
                      imageWidth:-1:1 ones(1, imageWidth-1); ...
                      imageWidth*ones(1, imageWidth) imageWidth-1:-1:1];
            oldImage = circshift(oldImage, [0 1]);
            numPoint = abs(yIndex(3) - yIndex(2) + 1) * abs(xIndex(3) - xIndex(2) + 1);
            corr(yIndex(1),xIndex(1)) = sum(sum(oldImage(yIndex(2):yIndex(3),xIndex(2):xIndex(3)) .* newImage(yIndex(2):yIndex(3),xIndex(2):xIndex(3)))) * imageHeight * imageWidth / numPoint;
        end
    end

    [value, yOffset] = max(corr(margin+1:end-margin,margin+1:end-margin));
    [dummy, xOffset] = max(value);
    offset = [yOffset(xOffset)+margin-imageHeight xOffset+margin-imageWidth];