I should also clarify that the above routine is based on the
correl(...) function as used in "Numerical Recipes in C"

Hi,
 I am attempting phase correlation
 between two sets of images: both are the same size, and are powers of
 two. The first contains an image in the top left hand corner; the rest

 of the image is black.
 The second contains the same image, offset by some amount; the rest of

 the screen is either gray or black.

Using the numerical recipes approach here is my correlation method:



void correl(float data1[], float data2[], unsigned long n, float ans[])

 {
         unsigned long no2, i;
         float dum, *fft;

        fft=vector(1,n<<1);

        twofft(data1-1, data2-1, fft, ans-1, n); // fft data1 and
data2

	//data1, data2 and ans are zero indexed, so subtract one from the
pointer
       // to the array header.

         no2=n>>1;

        for(i=1; i<=n+1; i+=2)
         {
                ans[i-1] = ((fft[i-1]*(dum=ans[i-1])) +
(fft[i]*ans[i]))/no2; // Real
 part - ans[0], ans[2] etc.
                 ans[i] = ((fft[i]*dum) -
(fft[i-1]*ans[i]))/no2;  // Imaginary part -
 ans[1], ans[3] etc.

                // Divide by the magnitude - this normalises
the results:
                 // see

http://ise.stanford.edu/class/ee392j/projects/projects/liang_report.pdf


                double normalisation = sqrt(ans[i-1]*ans[i-1] +
ans[i]*ans[i]);

                ans[i-1]/= normalisation;
                 ans[i]/=normalisation;

                // So the FFT matrix to inverse FFT is "ans".
         }

        ans[1]=ans[n];
         realft(ans-1,n,-1);
         free_vector(fft,1,n<<1);

}

One of the things I am not sure about, and would like to clarify, is
the
 division by the modulus of the complex arrays (the "normalisation"
 variable above). I've read the various papers, including the
 stanford.edu one above, and the Phase Correlation page on Wikipedia,
 and this seems like the thing to do. But I'd just like a second
opinion
 if possible.

Also, the even indexed values of ans[] will be real. and the odd will
be imaginary.
Do I just use the real values? If so, this would give a solution set of
size n/2 - which doesn't
sound right!
 

Best wishes 
 


Paul