Using MATLAB to calculate the offset between consecutive images

Question

I take images using a tunnel microscope. However, the scope drifts between successive images. I am 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 hours to process the 512x512 pixel images I'm dealing with. Do you have any suggestions for speeding up this code? Or do you know how to better track images in MatLab? Thank you 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 [mn] = 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 are calculating is called cross-correlation of 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 fill with zeros in the Fourier domain, you can use this math to get offsets in fractions of a pixel and apply offsets of fractions of a pixel to the image.

Answer 2

A typical approach to this problem is to take advantage of the fact that it works quickly for small images to your advantage. When you have large images, reduce them to make small images. Register small images quickly and use the calculated offset as your initial value for the next iteration. At the next iteration, you do not reduce the image as much, but start with a good initial estimate of the bias so that you can limit the search for solutions to a small one next to your initial estimate.

Although not written using tunneling microscopes, a review document that may be of some help is: “Reciprocal Information Recording of Medical Images: A Review” of Plume, Mainz, and Viergever published in IEEE Transactions on Medical Imaging, Vol. 22, No. 8, p. 986.

Answer 3

the link below will help you find the conversion between the two images and fix / restore the distorted ones (in your case, the image is 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); 

InlierPtsDistored, inlierPtsOriginal have attributes called locations. This is nothing more than matching the locations of one image with another. I think it’s very easy to calculate the offset from now on.

Answer 4

The following is my attempt to calculate the cross-correlation of two images manually. Something is not quite right. 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];