Sky Coyote, 7 Dec 2006

**Download**

**New Work**

The only thing new in this update is sub-pixel estimation of cross-correlation peaks by fitting a continuous quadratic surface to the discrete cross-correlation matrix. Although in general a quadratic surface can have the form:

z = a + bx + cy + dxy + ex^2 + fy^2

I have chosen a variation of this general form which seems appropriate for fitting to image cross-correlations. The specific equation is:

z = bias + xscale * ((x - xcenter) * cos(angle) + (y - ycenter) * sin(angle))^2 + yscale * ((y - ycenter) * cos(angle) - (x - xcenter) * sin(angle))^2

where

bias = peak value xcenter, ycenter = location of peak xscale, yscale = steepness of parabolic cross-sections (negative values for a concave-down surface) angle = orientation of major/minor axes about z axis

This surface might be called an elliptical paraboloid. An example of such a surface for {bias = 3.5, xcenter = 5.1234, ycenter = 4.5123, xscale = -0.5, yscale = -0.0625, angle = 0.5} is shown here:

The portion of the cross-correlation matrix used to fit this function includes the peak element of the result, plus all elements +/-5 pixels away (i.e. an 11x11 subarray). The method used to calculate the 6 parameters of the fit is my 'standard' evolutionary Monte-Carlo search algorithm, which is able to solve general non-linear parametric problems to any number of decimal places by minimizing an error function. In this case, the error is the sum of squared differences between the cross-correlation array and the quadratic surface.

As a test of this method, an 11x11 array was initialized using the function shown above. The algorithm first makes an initial guess of the parameters, and an estimate of the range of nearby values to search for the best fit. Here are the first few lines of the output of a typical run of 1000 iterations:

Iter Indx Error Bias Xcenter Ycenter Xscale Yscale Angle ---- ---- -------------- ---------- ---------- ---------- ----------- ----------- ---------- -1: 0 5.29597879E+03 3.4773E+00 5.0000E+00 5.0000E+00 -5.7476E-02 -5.7476E-02 0.0000E+00 0: 22 2.68762358E+03 2.0832E+00 5.6952E+00 1.8668E+00 -1.1426E-01 -5.1798E-02 6.9559E-01 1: 60 2.28892940E+03 6.5706E-01 6.5163E+00 6.1689E+00 -1.6777E-01 -3.5214E-02 4.0379E-01 2: 30 1.78243257E+03 1.9412E+00 3.3874E+00 8.2996E+00 -2.1286E-01 -6.1229E-02 3.1037E-01 3: 81 8.31036881E+02 1.6040E-01 1.8572E+00 9.3490E+00 -2.5737E-01 -8.3481E-03 5.4650E-01 4: 0 8.31036881E+02 1.6040E-01 1.8572E+00 9.3490E+00 -2.5737E-01 -8.3481E-03 5.4650E-01 5: 0 8.31036881E+02 1.6040E-01 1.8572E+00 9.3490E+00 -2.5737E-01 -8.3481E-03 5.4650E-01 6: 0 8.31036881E+02 1.6040E-01 1.8572E+00 9.3490E+00 -2.5737E-01 -8.3481E-03 5.4650E-01 7: 41 6.45613084E+02 5.5094E-01 4.5331E+00 5.2408E+00 -3.0178E-01 -5.5670E-02 5.9663E-01 8: 0 6.45613084E+02 5.5094E-01 4.5331E+00 5.2408E+00 -3.0178E-01 -5.5670E-02 5.9663E-01 9: 31 5.27037754E+02 2.4709E-01 2.6705E+00 8.8894E+00 -3.1463E-01 3.1159E-02 4.7820E-01

Here are the last few lines of the output:

990: 18 3.05886830E-14 3.5000E+00 5.1234E+00 4.5123E+00 -5.0000E-01 -6.2500E-02 5.0000E-01 991: 0 3.05886830E-14 3.5000E+00 5.1234E+00 4.5123E+00 -5.0000E-01 -6.2500E-02 5.0000E-01 992: 0 3.05886830E-14 3.5000E+00 5.1234E+00 4.5123E+00 -5.0000E-01 -6.2500E-02 5.0000E-01 993: 0 3.05886830E-14 3.5000E+00 5.1234E+00 4.5123E+00 -5.0000E-01 -6.2500E-02 5.0000E-01 994: 0 3.05886830E-14 3.5000E+00 5.1234E+00 4.5123E+00 -5.0000E-01 -6.2500E-02 5.0000E-01 995: 0 3.05886830E-14 3.5000E+00 5.1234E+00 4.5123E+00 -5.0000E-01 -6.2500E-02 5.0000E-01 996: 0 3.05886830E-14 3.5000E+00 5.1234E+00 4.5123E+00 -5.0000E-01 -6.2500E-02 5.0000E-01 997: 0 3.05886830E-14 3.5000E+00 5.1234E+00 4.5123E+00 -5.0000E-01 -6.2500E-02 5.0000E-01 998: 0 3.05886830E-14 3.5000E+00 5.1234E+00 4.5123E+00 -5.0000E-01 -6.2500E-02 5.0000E-01 999: 0 3.05886830E-14 3.5000E+00 5.1234E+00 4.5123E+00 -5.0000E-01 -6.2500E-02 5.0000E-01

A run of 1000 iterations (about 2 seconds) is more than sufficient for this problem, as the values of all parameters are stable to 4 decimal places in fewer steps. Here is the entire output file. Note that due to the nature of this problem, there are several equivalent solutions, including angles +/- pi:

999: 0 1.75639276E-11 3.5000E+00 5.1234E+00 4.5123E+00 -5.0000E-01 -6.2500E-02 -2.6416E+00

and with the x and y axes switched:

999: 0 1.30632103E-13 3.5000E+00 5.1234E+00 4.5123E+00 -6.2500E-02 -5.0000E-01 -1.0708E+00

Nevertheless, all combinations represent the same surface.

In actual use, one first computes the cross-correlation using the direct or FFT methods, and then performs a quadratic fit to find the location of the peak in the surface. The two buttons in the XC window perform the following functions:

X-correlate: Perform cross-correlation using selected method and shift on indicated subimage rectangles. Plot resulting matrix autoscaling between max and min values. Quadfit: Fit a quadratic surface to the 11x11 subarray around peak of the XC matrix currently shown. Plot result using same scale as before.

Unfortunately, performing a quadratic fit does not always produce a more accurate result. For example, Here is the XC of an image with itself:

While the direct integer shift XC yields the correct result, the quadratic fit leads one to believe there is a slight fractional shift:

The Fourier transform XC yields a tighter quadratic fit in this case:

If the subimage is sufficiently large, the quadratic fit can detect fractional shifts with some accuracy (in this case, the actual shift was {10.5, -10.5}):

However, it doesn't help the Fourier XC:

Using actual data, and different images, one might believe the quadratic fit:

In this case, atan(0.928925/17.9365) = 2.96 degrees, which looks plausible. However, only time and use will tell if this method is worthwhile.

**Still To Do**

* Sub-pixel cross-correlations via pixel interpolation.

İSky Coyote and Dr. Eliot Young, 2006