Registers two images (2-D rigid translation) within a fraction of a pixel specified by the user. Instead of computing a zero-padded FFT (fast Fourier transform), this code uses selective upsampling by a matrix-multiply DFT (discrete FT) to dramatically reduce computation time and memory without sacrificing accuracy. With this procedure all the image points are used to compute the upsampled cross-correlation in a very small neighborhood around its peak. This algorithm is referred to as the single-step DFT algorithm in .
 Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup, "Efficient subpixel image registration algorithms," Opt. Lett. 33, 156-158 (2008).
Please refer to the attached HTML for more details and a sample implementation.
Great code, I use it for all of my microparticle tracking needs and I'll definitely be citing in an upcoming paper :)
One very simple change I made that gave me a 10-20% increase in speed was changing line 115:
CCabs = abs(CC);
to the equivalent but slightly faster:
CCabs = sqrt(real(CC).^2 + imag(CC).^2);
not sure exactly why this speeds things up but any user should definitely consider this change.
imregister in the Image Processing Toolbox will do it. Or there are alternatives like Nifty or Elastix. The latter I made a really simple wrapper for: https://www.mathworks.com/matlabcentral/fileexchange/52982-matlab-elastix
Hi Manuel. Thanks for this great function; works really well.
I was wondering if you could help me understand where it might be changed to work for 3D arrays? If you can point me in the right direction I'll give it a go, or if it's easy for you could you create a 3D version?
Thanks in advance for your time.
great piece of work
Great code! Is there a way to use the output from the registration of one image to a reference image to transform a second image with the same algorithm? Thanks!
Working good for only translational error images. But is it possible to extend to cover angular rotations and scale shifts as well? (e.g. Fourier Mellin transform) as commented earlier by Philip. Any update in this code?
Does this algorithm only work for translation? How about image rotation?
is there a way to limit the rotation angle?
I want to use to to rotate square images to a shape template (a binary frame). It then randomly rotates the images n*89° instead of -1°, which would be the closest.
I'm clueless at this point how to implement that, does anyone have an idea?
This is great. I made a few minor modifications to limit translation to only the x or y axis. Do you have this under version control anywhere?
Simply fantastic. So much more efficient and accurate for 2D translation than everything else I tried. The subpixel resolution is a bonus as well.
It would be really useful if the code could do rotation as well.
Is this only for shifts? can we register rotation as well?
wow, that one is really fast and seems to be robust. thanks a lot!
Great algorithm, it has helped so much with my image stabilization project
works very well. i'm amazed with the speed and accuracy of the algorithm.
do images need to be square - as in 2^n for the fft shift to work properly?
can this be applied to colour images in some way?
ThX, I did it, hopefully some one answer me.
Nafise, post your situation/goals, your images, and what you've tried to the newsgroup (http://www.mathworks.com/matlabcentral/newsreader/). From there people can help you.
I did it, but it is not true, because my image has black regions, and it determines that part with zero pad part.
WHAT SHOULD i DO NOW?!
Nafise b, pad your smaller image with zeros so that it is the same size.
Thanks for your code.
I want to work with your code, but I have a problem with it.
In fact, my images are 2 blocks which don't have the same size, and I want to find the best place for small block in the big block.
Something that I got from your code is for your code 2 images have same size.
What should I do now, if I want to use your code for subpixel image registration?
Great work! Question about the codes " kernr=exp((-i*2*pi/(nr*usfac))*( [0:nor-1].' - roff )*( ifftshift([0:nr-1]) - floor(nr/2) )) " :
Why does the vector [0:nr-1] need be ifftshift ?
Re: my comment on 14 Jul 2010:
I owe an apology. It seems at the point when I wrote my comment I was feeding an accidentally reshaped array into the code. Because of the artifacts in my data, it caused a wrap around artifact in the output (and the registration worked in one axis only, but I didn't notice because that happened to be the axis containing almost all of the shifts).
Great code - no resrvations! Thanks much.
Great code. Any plans to extend to cover angular rotations and scale shifts as well? (e.g. Fourier Mellin transform)
Code registers images very well. Question about the supplementary output:
output = [error,diffphase,net_row_shift,net_col_shift]
On my data it seems that some times the row_shift variable seems to 'wrap around' in away that is not intuitive. For example though the data is only 512 pixels wide, I get row shifts of -512.3, 1024.3 etc which I assume are actually very small pixel shifts that have a somewhat random multiple of the image size thown on for fun.
Wonderful. Fast and accurate.
Is there a way to extend this to 3D CT images?
Thanks for your sharing.
I want to register the pair of visible and theraml images with this one, but I don't get the required resutl.
Tremendously useful! Thanks very much.
Works right out of the box. Thanks for sharing !
Great Work - fast and accurate!!!
This method is equally robust to noise as the standard FFT upsampling, but much faster. It computes the same upsampled cross correlation but only in a small neighborhood around the peak.
Thanks to everyone for the positive comments! The published M-file is intended to get the user started with the syntax and explain what the code does. For technical documentation please refer to the reference above  and references within.
forgot to rate. wow this code is brutally fast and accurate, if you define a ROI even faster. this code needs some kind of award. how much faster would it be if it was compiled in to a mex?
works but i had to re=write code for UINT16 images, very fast and accurate. why are you using double? do you really need that accuracy in the dft?
works but i had to re=write code for UINT16 images, very fast and accurate.
Is this method less robust to noise than standard upsampling?
Always nice to see a published .m, but the images tell us nothing. A zoomed in display, abs error or something would be more relevant than two identially looking images.
Elegant solution in simple understandable code + sound documentation.
Thank you for you job, I will enjoy it.
Very good and excellent documentation to get you started with the syntax. Great job and thanks for including the paper reference.
did i mention it was easy to use?
easy to use, extremely fast, and easy to use. amazing work.
Revise license text
Rewrote portions of the code for improved consistency and better readability, added BSD license.
Win cash prizes and have your live script featured on our websiteLearn more
Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.