1
$\begingroup$

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:

$\endgroup$
  • $\begingroup$ Welcome to Math.SE! Could you explain some of the symbols you are using, such as the star and the frak $F$? For Math.SE, it might also be useful to not post your algorithm as code, but use formulas or pseudo-code instead. In this way, mathematicians who can't read code can still help you and allows more people to benefit from your question. $\endgroup$ – Hrodelbert May 26 '15 at 12:50
  • $\begingroup$ Its updated to reflect your notes, hoping it is clearer now! $\endgroup$ – arve0 May 26 '15 at 16:09
1
$\begingroup$

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:

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.