# Correct method for up-sampling cross-power spectrum for sub-pixel motion estimation

Let's say I have two images, one reference($g$) and another shifted($f$). I tried to get subpixel accuracy of the motion. Here are my steps:

1. Firstly I would calculate the normalized cross-power spectrum $C(u,v)$, where $F$ and $G$ are fft2-ed.

$$C(u,v)=\dfrac{F(u,v)G^*(u,v)}{\lvert F(u,v)G^*(u,v) \rvert}$$

1. Next I would like to get $1/k$ subpixel accuracy, by padarray the $C(u,v)$ by $k$ factor, resulting in $C_k(u,v)$

2. ifft2 $C_k(u,v)$, resulting in $c_k(x,y)$

3. Find the arg max of the $c_k(x,y)$

Here's my code in Matlab:

reference = fft2(reference);
shifted = fft2(shifted);
C=shifted.*conj(reference)./abs(shifted.*conj(reference));
[m, n]=size(reference);
[~,w] = max(real(ck(:)));
[dxk,dyk] = ind2sub(size(ck),w);


After run the code I can't get the accurate result. I hope I explained it clear enough(from what I understood from the attached paper). I really need help on this and much appreciate any help!

Referred paper: Here (start from section 3)

# Preface

I know the OP asked for Matlab, but if you follow our discussion here in the Mathematica chat it seems any answer will be helpful. Therefore, I'm going to present the approach, visualizing each step, so it can be repeated and debugged in Matlab.

The approach itself is simple and as far as I can see there are only two things to keep in mind. But first things first. I'm using a downscaled version of the Lena image that can be accessed in Mathematica through ExampleData. I shifted it by {-5.3, 8.7} pixels which can be done through ImageTransformation in combination with TranslationTransform.

{dx, dy} = {-5.3, 8.7};

gImg = ImageResize[
ColorConvert[ExampleData[{"TestImage", "Lena"}], "Grayscale"], 128];
{nx, ny} = ImageDimensions[gImg];
fImg = ImageTransformation[gImg, TranslationTransform[{dx, dy}/{nx, ny}]];

{gData, fData} = ImageData /@ {gImg, fImg};


Keep in mind that shifting the image by fractions of pixel units will resample it and introduce some blur that degrades the image. Therefore, I would not expect that the perfect solution can be found.

What we have in gData and fData are the pixel values of the images which look like this:

### Fourier transform

I don't know how fft2 in Matlab exactly works, but usually, taking the Fourier transform of an image does not result where the frequencies are centered. Let me give you an example, by showing how the spectrum looks if you just use FourierTransform

Image[Fourier[gData] // Abs]


To make the upscaling by padding work, we need to center the transform. You can do this, by multiplying each pixel by $(-1)^{x+y}$ which gives a chessboard-like fashion. This is explained in detail in the very easy to read Digital Imageprocessing by Gonzalez and Woods. A pdf of the 2nd edition can be found here and you should follow the discussion on page 154.

In Mathematica, a function to transform an image into this chessboard-like form can be written down as

centerFrequency[data_] := With[{d = Dimensions[data]},
Table[(-1)^(i + j)*data[[j, i]], {j, Last[d]}, {i,
First[d]}]
];


If we now look at the same spectrum, with this function applied before using the FFT, we get the expected result:

Image[Fourier[centerFrequency@gData] // Abs]


### Normalized Cross Power Spectrum

To create an upscaled version of the normalized cross-power spectrum, we follow exactly what is written in the paper

1. Calculate the centered fourier transform of fData and gData into f and g
2. Conjugate g
3. Calculate f*g/Norm[f*g]
4. Pad the result by a factor k of the original image dimension (I assume a square image in this example)

In code this looks like this

normCrossPow[ref_, tmpl_, k_] := Module[{
f = Fourier[centerFrequency@ref],
g = Conjugate[Fourier[centerFrequency@tmpl]], nx
},
nx = First[Dimensions[ref]]/2;
]


To get the resulting spectrum, we only need to apply InverseFourier, apply centerFrequencies again and look at the real part of the complex array.

spec = Re[
centerFrequency@InverseFourier[normCrossPow[fData, gData, 3]]
];


Above, I used a factor of k=3. The white spot in the upper right corner is our maximum we are interested in. Let us find the position of this maximum

max = FirstPosition[spec, Max[Flatten[spec]]]


This gives {36, 22}. Keep in mind, that the y-axis in image coordinates and matrix-coordinates is reversed and that in a matrix the y-values comes first. Additionally, if images have no shift at all, the maximum would be at {1,1} and we need to divide out our factor k again. Note that we divide by 4 (!!) because with a factor of k=1 we would divide by 1, with k=1 by 2, and so on.

({36, 22} - {1, 1})/4.0


This gives {8.75, 5.25} for {dy, dx} but remember that it is actually dy=-8.75 because of the reversed y-axis.

Now, we can use this result and back-transform our fImg

res = ImageTransformation[fImg, TranslationTransform[{5.25, -8.75}/{nx, ny}]]


and look at the combined image of gImg and res

Fits nicely and the only error we see is the black border that was introduced by initially transforming fImg.

• Thank you so much for this! It really meant a lot to me! However, when I re-implement the steps in matlab, I could get accurate dx, but not dy. Any ideas why was it? I would get something like this [dx=3.5000 dy=247.2500, which the dy is relatively large. The actual should be around dy=8.7400 This is the code, it should be understandable in both platform. Thank you once again! – Gregor Isack Mar 5 '17 at 8:29
• The same can happen to dx if the shift is in the other direction! Just think about it, when the maximum (now in the upper left corner) would move further to the left, it is reflected to the right side of the image. The same happens with dy. You need to calculate something like ny-dy if dy is larger than ny/2. Here, ny is the number of pixels in y direction. – halirutan Mar 5 '17 at 9:10
• Yes! How insensitive of me! Thank you so so much @halirutan ! You saved my day! And sorry for disturbing your Saturday :) – Gregor Isack Mar 5 '17 at 9:23