Why is my phase correlation not equal to the real correlation?

Question

If I understand the correlation theorem correctly, it states:

$ f(x,y) \unicode{x2606} \bar g(x,y) = \mathfrak{F}^{-1} \left\{ F^*(u,v) G(u,v) \right\}, $

also called a phase correlation. Above $f(x,y)$ and $\bar g(x,y)$ are 2D images of the same "size", in the sense that we define $\bar g(x,y)$ as

$\bar g(x,y) = g(x \bmod m, y \bmod n), $

yielding $g(x,y)$ as an infinite sequence of the same repeating image. In other words $\bar g(0,-1) = g(0,n-1)$, $\bar g(0,0) = g(0,n)$ and so on, where the size of $g(x,y)$ is $m \times n$.

$\mathfrak{F}^{-1}$ is the inverse discrete Fourier transform. $F(u,v)$ and $G(u,v)$ are the discrete Fourier transforms of respectively $f(x,y)$ and $g(x,y)$. $F^*(u,v)$ is the complex conjugate of $F(u,v)$.

For those unfamiliar with correlation, it is defined as:

$ f(x,y) \unicode{x2606} g(x,y) = \sum\limits_{s=-a}^a \sum\limits_{t=-b}^b f(s,t) g(x+s, y+t). $

Here $a=\frac{m-1}{2}$ and $b = \frac{n-1}{2}$, $m \times n$ being the size of $f(x,y)$ (assuming odd sizes).

If the above is correct, one should assume that the correlation and phase correlation calculations in psuedo code below yielded the same answer (give or take the imprecision of the storage unit). But its not even close (tried in matlab and python).

f = rand(3, 3)
g = rand(3, 3)
g_bar = pad(g, size=1, mode=wrap)

f_dft = fft2(f)
g_dft = fft2(g)

c = correlate(f, g_bar)
ph_c = ifft2(conjugate(f_dft) * g_dft)

diff = c - real(ph_c)

print(diff)

Why is this not the case? Is my understanding of the correlation theorem wrong?

Here are the psuedo code above in:

• python
• matlab

Answer 1

My assumption that $g(x,y)$ should be "wrap-padded" was wrong. Cross correlation is zero padded, so the image should be zero padded before taking the DFT.

In other words,

$\bar g(x,y) = \begin{cases} & g(x,y) \text{ if } x \in [0,m) \text{ and } y \in [0,n) \\ & 0 \text{ if } x \notin [0,m) \text{ and } y \notin [0,n) \\ \end{cases}$

In psuedo:

size = 3
padding = size//2    # integer division

f = rand(size, size)
g = rand(size, size)

# cross correlation is defined with zero padding!
c = correlate(f, g, mode="constant")

# zero pad to keep frequency information on edge
f_padded = pad(f, padding, mode="constant")
g_padded = pad(g, padding, mode="constant")

# phase correlation
F = fft2(f_padded)
G = fft2(g_padded)
FG = conjugate(F) * G
ph_c = ifft2(FG)

# why do we need to shift?
ph_c = fftshift(ph_c)

# clip back to original size
ph_c = ph_c[padding:-padding, padding:-padding]

# reverse axes - why?
ph_c = ph_c[::-1, ::-1]

diff = abs(c - real(ph_c))
print(diff < 1e-12)

As seen from the comments in the psuedo code, two new questions arose. So this is kind of a half answer.

Here's the implementation in:

• python
• matlab