Thread Subject: Cross-correlation with sub-pixel accuracy

I'm performing phase only correlation between 2 images to find the displacement using FFT's, which gives a discrete correlation plain. Currently I'm finding the peak location simply from the maximum value of the correlation plane. Can anyone help me find a method to find the peak location with sub-pixel accuracy by looking at the neighbouring pixels around the peak? I've seen method that use curve fitting but this is too slow. I'm looking for a simple method that can be implemented quickly without iteration. Perhaps some kind of averaging of neighbouring locations using the values as weights?
Peter

"Peter Bone" <peterbone@hotmail.com> wrote in message <gombji$mri$1@fred.mathworks.com>...
> I'm performing phase only correlation between 2 images to find the displacement using FFT's, which gives a discrete correlation plain. Currently I'm finding the peak location simply from the maximum value of the correlation plane. Can anyone help me find a method to find the peak location with sub-pixel accuracy by looking at the neighbouring pixels around the peak? I've seen method that use curve fitting but this is too slow. I'm looking for a simple method that can be implemented quickly without iteration. Perhaps some kind of averaging of neighbouring locations using the values as weights?
> Peter

Once you have both images in the frequency domain, do an fftshift to put DC (0Hz) in the centre of the plot. Now embed the image - without changing its scale, into a bigger array, such that the DC point is still in the centre, and set the extra frequency locations you have just created to zero.

Now do your correlation (you might want to do an ifft shift to put DC back in the 1,1 position first). You will find that your correlation output is the same size as your augmented array, with each element representing a fractional pixel shift.

Example, assume you start with a 128 x 128 pixel image. FFT it, shift, and embed it into a 512 x 512 pixel frame, then each element on your correlation should represent 0.2 of a pixel.

Hope that helps

Dave Robinson

"Dave Robinson" <dave.robinson@somewhere.biz> wrote in message <gome7u$rgb$1@fred.mathworks.com>...
> "Peter Bone" <peterbone@hotmail.com> wrote in message <gombji$mri$1@fred.mathworks.com>...
> > I'm performing phase only correlation between 2 images to find the displacement using FFT's, which gives a discrete correlation plain. Currently I'm finding the peak location simply from the maximum value of the correlation plane. Can anyone help me find a method to find the peak location with sub-pixel accuracy by looking at the neighbouring pixels around the peak? I've seen method that use curve fitting but this is too slow. I'm looking for a simple method that can be implemented quickly without iteration. Perhaps some kind of averaging of neighbouring locations using the values as weights?
> > Peter
>
> Once you have both images in the frequency domain, do an fftshift to put DC (0Hz) in the centre of the plot. Now embed the image - without changing its scale, into a bigger array, such that the DC point is still in the centre, and set the extra frequency locations you have just created to zero.
>
> Now do your correlation (you might want to do an ifft shift to put DC back in the 1,1 position first). You will find that your correlation output is the same size as your augmented array, with each element representing a fractional pixel shift.
>
> Example, assume you start with a 128 x 128 pixel image. FFT it, shift, and embed it into a 512 x 512 pixel frame, then each element on your correlation should represent 0.2 of a pixel.
>
> Hope that helps
>
> Dave Robinson

Thanks. I've developing this in Matlab currently but my aim is to get it to work on a DSP. The upscaling method will be too slow. Here's a very naive method I just wrote that gives a lot better results than simply finding the integer pixel location. Can anybody think of a reason why this is so wrong (it's fast)? It analyses a 5x5 region around the maximum peak (a). It gives (3.5,3.5) as the peak location instead of the max peak location of (3,3)

a = [1 1 2 2 2 ;
     1 5 10 10 5 ;
     2 10 100 99 10 ;
     2 10 99 99 10 ;
     2 5 10 10 5 ];
 
 x = 0;
 y = 0;
 sum = 0;
 for r = 1 : 5
     for c = 1 : 5
         sum = sum + a(r,c);
         x = x + c * a(r,c);
         y = y + r * a(r,c);
     end
 end
 
 x = x / sum;
 y = y / sum;

Peter

Apologies for posting the original 3 times. I have no idea how that happened.

Have you tried the elegant solution of Guizar et al, posted on the MATLAB Central File Share?