FITSRegister Update 12/22/06
Sky Coyote, 22 Dec 2006

Sub-Pixel Cross-Correlations

This update concerns finding the sub-pixel peak of a cross-correlation surface between two images by using image interpolation. For an detailed example of using a quadratic fit to estimate the location of this peak, see the previous update. The quadratic fit in 6 parameters does not appear to be sufficient to accurately locate the XC peak. For example, here is a pair of test images, where the right one has been shifted by (10.625, -10.375) pixels (a 1/8th pixel resolution). This is an exact shift, as both images were created using continuous functions (superimposed 2d Gaussians). Thus, each image is a continuous function sampled at a discrete regular grid of 512x512 points: A 32x32 subimage (domain) extracted from the shifted image was cross-correlated with a 96x96 subimage (range) from the first image using integer shifts and direct normalized parametric correlation. The peak of the resulting surface was found at:

```xc(65, 65) max = 0.998819 at:
11, -11
```

A 6 parameter elliptical paraboloid was fit to an 11x11 region around this peak, yielding the following values:

```11x11 quadfit: bias=0.998242, xcenter=10.2893, ycenter=-10.8163, xscale=-0.00150941, yscale=-0.00879815, angle=-1.25059
``` The innaccuracy of this result is typical for quadractic fits using smaller or larger regions around the integer peak location:

``` 3x 3 quadfit: bias=1.00009, xcenter=10.6447, ycenter=-10.5465, xscale=-0.00966501, yscale=-0.00210171, angle=0.312934
5x 5 quadfit: bias=1.00009, xcenter=10.5489, ycenter=-10.4866, xscale=-0.00172648, yscale=-0.00948908, angle=-1.26098
7x 7 quadfit: bias=1.00006, xcenter=10.474, ycenter=-10.5833, xscale=-0.00166612, yscale=-0.0093413, angle=-1.25798
9x 9 quadfit: bias=0.999622, xcenter=10.3854, ycenter=-10.6971, xscale=-0.0091127, yscale=-0.00159235, angle=0.316312
11x11 quadfit: bias=0.998242, xcenter=10.2893, ycenter=-10.8163, xscale=-0.00879815, yscale=-0.00150941, angle=-2.82138
13x13 quadfit: bias=0.995309, xcenter=10.1913, ycenter=-10.9303, xscale=-0.00142313, yscale=-0.00840713, angle=-1.2462
15x15 quadfit: bias=0.990296, xcenter=10.0962, ycenter=-11.0298, xscale=-0.00796007, yscale=-0.00133873, angle=3.47123
17x17 quadfit: bias=0.982878, xcenter=10.008, ycenter=-11.1069, xscale=-0.00748189, yscale=-0.00126005, angle=0.335491
19x19 quadfit: bias=0.972976, xcenter=9.9303, ycenter=-11.1558, xscale=-0.00118934, yscale=-0.00699593, angle=-1.22856
21x21 quadfit: bias=0.960718, xcenter=9.86598, ycenter=-11.1727, xscale=-0.00652047, yscale=-0.00112746, angle=0.349843
```

One could of course experiment with other fitting functions (e.g. asymmetric Gaussians), but there may be little to be gained from this approach. Instead, I have tried a more obvious 'brute force' approach to achieving sub-pixel estimates by enlarging the images to be cross-correlated. However, simply duplicating the pixels does not yield any increase in accuracy. Here is an example of cross-correlating two subimages which have been enlarged 8x by pixel duplication:

```xc took 0 secs
xc(65, 65) max = 0.998819 at:
11, -11
xc took 226 secs
xci(513, 513) max = 0.998819 at:
11, -11
``` Note that although the resulting XC surface appears nice and smooth, it does not contain any additional information over the non-enlarged result. Note also that while the initial XC calculation took less than one second, the enlarged calculation took almost 4 minutes. What is needed is a way to interpolate the subimages before performing the enlarged XC so as to find a new peak at a different location in the resulting surface. There are zillions of ways to interpolate an image, and a Google search will find some of them. Although some claim to be better than others for specific purposes, bear in mind that they all share a defect in that they involve creating fake data where none existed before. With that said, my attitude as an engineer is to start by trying the simplest of these (also the computationally fastest) until one is found that does the job. To this end I have created a very simple function that interpolates an image by powers of 2 (actually, by (size * 2) - 1 at each step) which works as follows: At each step, 5 new pixels are created between every 4 existing pixels. Four of the new pixels (on the X and Y gridlines) are the average of the two pixels nearest them (on the same X or Y gridline). The new center pixel is the average of the four corner pixels. This interpolation scheme has the following effect on the enlarged XC calculation:

```xc took 0 secs
xc(65, 65) max = 0.998819 at:
11, -11
xc took 213 secs
xci(513, 513) max = 0.999999 at:
10.625, -10.375
``` Note that not only does this change yield the correct shift result, but the match is now exact (0.999999... = 1). Pretty good for such a simple interpolation method! However, the 8x enlarged XC calculation still takes a very long time. To address this problem, I have created an iterated interpolation and XC calculation that performs a 2x interpolation and XC of both subimages within only the new 5x5 region around the peak of the previous XC surface. Thus, at each iteration, a new bit of information in the overall fractional shift is obtained (e.g. {1, 0.5, 0.25, 0.125, 0.0625...} pixel shifts), without the need to perform the direct XC computation for the entire interpolated subimages. This iterated calculation produces the same result in a fraction of the time:

```xc took 0 secs
xc(65, 65) max = 0.998819 at:
11, -11
n = 0, sub2i(63, 63), sub1i(67, 67)
xci(5, 5) max = 0.999752 at:
-1, 1
n = 1, sub2i(125, 125), sub1i(129, 129)
xci(5, 5) max = 0.999886 at:
1, 0
n = 2, sub2i(249, 249), sub1i(253, 253)
xci(5, 5) max = 0.999999 at:
-1, 1
n = 3, sub2i(497, 497), sub1i(501, 501)
xci(5, 5) max = 0.999999 at:
0, 0
xc took 1 secs
total shift = (10.625, -10.375)
``` In this case, 4 iterations (1/16th pixel resolution) took about 1 second (which also included the time to plot the images). At each iteration, the subimages are interpolated by a factor of 2x, and a new XC is performed in the 5x5 region around the previous peak. There is no limit to the number of iterations that could be performed, each yielding another factor of 2x in the resulting resolution. However, the time and memory required at each step increases at least as O(n^2) (and possibly O(n^4)), so that in practice the effective resolution is currently limited to about 6-7 steps.

To test the robustness of this method, I created a set of multiple tests that extracted random subimages of size 16x16 up to 64x64 from the shifted image and performed iterated XCs wrt subimages of 3x that size in the unshifted image. Here are two excerpts from a 4 iteration test:

```  n    x0    y0   dx   dy       shiftx       shifty
---   ---   ---   --   --   ----------   ----------
1   150    25   20   20    10.625000   -10.375000
2   308   110   62   62    10.625000   -10.375000
3   269   156   56   56    10.625000   -10.375000
4   193   110   54   54    10.625000   -10.375000
5   213    88   39   39    10.625000   -10.375000
6   108   288   51   51    10.625000   -10.375000
7   286   367   62   62    10.625000   -10.375000
8    91   413   22   22    10.625000   -10.375000
9   295   237   22   22    10.625000   -10.375000
10   173   168   50   50    10.625000   -10.375000
...
91   387   453   20   20    10.625000   -10.375000
92   186   216   38   38    10.625000   -10.375000
93   325    71   36   36    10.625000   -10.375000
94   152   218   64   64    10.625000   -10.375000
95   267   453   28   28    10.625000   -10.375000
96   181   329   31   31    10.625000   -10.375000
97   467    78   16   16    10.625000   -10.375000
98   180   347   45   45    10.625000   -10.375000
99   374   294   62   62    10.625000   -10.375000
100   299   273   48   48    10.625000   -10.375000
run took 53 sec, 0.53 secs / xc
```

The entire results can be seen here. In all cases the correct shift was found. To further test the method, I created an additional set of slightly more complicated (multiple Gaussians of varying size) shifted and unshifted images as shown below. The right image is shifted by (10 + 63/64, -(10 + 1/64)) pixels: A 6 iteration set of XC tests on these two images produced the following output:

```  n    x0    y0   dx   dy       shiftx       shifty
---   ---   ---   --   --   ----------   ----------
1    88   162   60   60    10.984375   -10.015625
2   329   254   37   37    10.984375   -10.015625
3   422   159   24   24    10.984375   -10.015625
4   213   400   53   53    10.984375   -10.015625
5    38   206   21   21    10.984375   -10.015625
6   274    30   30   30    10.984375   -10.015625
7   157   332   60   60    10.984375   -10.015625
8   312   127   55   55    10.984375   -10.015625
9   264   196   42   42    10.984375   -10.015625
10   352   215   60   60    10.984375   -10.015625
...
91   138   340   52   52    10.984375   -10.015625
92    72   398   42   42    10.984375   -10.015625
93   365    48   45   45    10.984375   -10.015625
94   128    92   45   45    10.984375   -10.015625
95   196   326   56   56    10.984375   -10.015625
96   220   355   50   50    10.984375   -10.015625
97   355   170   33   33    10.984375   -10.015625
98   314    80   59   59    10.984375   -10.015625
99   123    38   25   25    10.984375   -10.015625
100    56    74   33   33    10.984375   -10.015625
run took 472 sec, 4.72 secs / xc
```

The entire results can be seen here. Again, in all cases the correct shift was found. As a final test, I extracted a subset of one of the processed composite Venus images and shifted it by an irrational amount of (3.14159265358979, -2.71828182846) using FITSRegister and the FFT sinc shift: I then performed multiple tests on the shifted/unshifted subsets using only 32x32 domains and 96x96 ranges, with 6 levels of iteration:

```  n    x0    y0   dx   dy       shiftx       shifty
---   ---   ---   --   --   ----------   ----------
1    53    42   32   32     3.140625    -2.718750
2    42    98   32   32     3.140625    -2.718750
3    62    48   32   32     3.140625    -2.718750
4    83    93   32   32     3.140625    -2.718750
5    45    62   32   32     3.140625    -2.718750
6    86    89   32   32     3.140625    -2.718750
7    64    61   32   32     3.140625    -2.718750
8    73    52   32   32     3.140625    -2.718750
9    48    95   32   32     3.140625    -2.718750
10    63    46   32   32     3.140625    -2.718750
run took 24 sec, 2.4 secs / xc
```

These results are well within 1/64th pixel of accuracy. Of course the real test will be to try this method on different Venus images which are separated in time. However, in that case I don't know how the result will be judged, since I doubt it is possible to align image subsets to 1/64th pixel by eye.