2D Cross-Correlation using 1D FFT

Question

I need to calculate cross correlation between 2 images which I read in 2 vectors, both of them uni-dimensional.

As per my understanding I need to do the following operations:

1. FFT(1D) on vector 1 and vector 2
2. Multiplication between the outputs of the FFT applied on each vector
3. IFFT(1D) of the multiplication

If I'd wanted to use 2D FFT in the first step following this method then what would the second step be?

Answer 1

Your 2nd step is wrong, it's doing circular convolution. For circular cross-correlation, it should be:

2. Multiplication between the output of the FFT applied on the first vector and the conjugate of the output of the FFT applied on the second vector.

Aside from that, the steps are the same whether in 1D or 2D:

$$\texttt{iFFT}_{2d}\bigg(\texttt{FFT}_{2d}\big(V_1\big) \cdot \overline{\texttt{FFT}_{2d}\big(V_2\big)} \bigg)$$

You can zero-pad the inputs $V_1$ and $V_2$ to calculate non-periodic cross-correlation using circular cross-correlation.

Answer 2

As Jdip correctly pointed out, your implementation is doing circular convolution, not cross-correlation. But I want to add a nuance of discrete circular cross-correlation that makes his answer slightly incorrect.

Frequency Conjugation for Discrete vs. Continuous Signals

There are two key differences between cross-correlation and convolution:

• In cross-correlation, one of the vectors is conjugated (in the time domain)
• In convolution, one of the vectors is reversed/flipped

Thus, to perform cross-correlation via FFT-implemented circular convolution, we must pre-flip and conjugate one of the vectors: cross_correlation(x,y) = convolution(flip(conjugate(x)), y).

The conjugation property of the continuous time Fourier transform says that conjugating in the frequency domain conjugates and flips in the time domain: $\mathcal{F}^{-1} (\overline{X(\omega)}) = \overline{x(-t)}$ . So you should be able to perform cross-correlation via FFT by just conjugating one of the FFTed vectors, right? Wrong!

The reason is that for discrete signals such as images, conjugation in the Fourier domain does not equate to time reversal. Rather, it performs modulo N time reversal: $ [x_1, x_2, x_3,... x_N] => [x_1, x_N, x_{N-1}, x_{N-2},... x_2] $.

Solution

Because of this, I recommend doing the conjugation and flipping in the time domain.¹ Here are the steps:

1. Flip and conjugate $V_2$: $W_2 = \overline{\texttt{flip}(V_2)}$
2. Take the FFT of $V_1$ and $W_2$
3. Multiply FFTed vectors
4. Take the IFFT of the product

Thus:
$$ \texttt{IFFT}_{2D}(\texttt{FFT}_{2D}(V_1) \cdot \texttt{FFT}_{2D}(W_2)) $$

Or:

$$ \texttt{IFFT}_{2D}(\texttt{FFT}_{2D}(V_1) \cdot \texttt{FFT}_{2D}(\overline{\texttt{flip}(V_2)})) $$

Edit

OverLordGoldDragon has argued in the comments that we shouldn't reverse/flip $V_2$ to perform cross-correlation via convolution. I offer a simple proof in MATLAB that we should:

%Create vectors and cross-correlate them
N = 8; M = 3; 
v1 = randi(9, N, 1);        %Gives a vector of random integers less than 10
v2 = v1(3:3+M-1);           %Take a small subset of v1 as our second vector
P = N + M - 1;              %Length of full convolution
y = xcorr(v1, v2);          %Cross-correlation 
y = y(N-M+1:end);           %Get rid of extra xcorr lags (zeros)
%Now, do this with convolution 
%NOTE: No time domain conjugation because v2 is real: v2 == conj(v2)
y2 = conv(flip(v2), v1);    %y2 and y are identical
%Try conjugating in the frequency domain instead of flipping
%y3 and y are not the same!
y3 = conv(ifft(conj(fft(v2))), v1); 

Footnote

¹More efficient alternative: Use a linear phase multiply in the frequency domain to circularly shift the time domain vector by 1 sample before (or after, depending on the direction of the shift) conjugating in the frequency domain.

Answer 3

Take an image:

Cross-correlating it with impulse should yield itself, and cross-correlating with itself should peak at center. Key points:

• The operating kernel must be centered about $t=0$. For a discrete sequence $h$ of length $N$, under the FFT, this means $h[0]$, and the second-half of samples are of negative time: $h[n > N/2]$
• Assuming inputs are normal images, this means time-centering the "template" image: note $x \star h \neq h \star x$.
• fft2(x) == fft(fft(x, axis=0), axis=1)

cross correlation between 2 images which I read in 2 vectors, both of them uni-dimensional.

That's a problem since the time-centering step isn't straightforward to replicate on a flattened image, I'm unsure how it'd be done. This answer will assume a 2D array, so you can just reshape it into 2D and then back. Finally, you might want to look into padding and boundary effects.

Putting it together, here's cross-correlation of COVID with image-centered unit impulse, and with itself:

additionally, we move it and its flipped copy, and check that the more intense dot spots the unflipped version; note boundary effects:

$$ \texttt{iFFT}_{2d}\bigg( \texttt{FFT}_{2d}\big(x\big) \cdot \overline{\texttt{FFT}_{2d}\big(\texttt{iFFTSHIFT}_{2d}(h)\big)} \bigg) $$

Centering / performance note

Answer assumes we're after circular cross-correlation with no padding, and $x$ and $h$ are "normal" (non-DFT centered) images. With padding/unpadding, the most compute efficient route may be different.

Complete cross-correlation, with padding, is implemented here.

Meta note

This answer was subject of controversy. If readers aren't sure which answer to use, I encourage reading responses to the followup, where I've shown further demos and offered 1000 bounty to invalidate this answer, instead of trusting the votes.

Full code

import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import fft, ifft, fft2, ifft2, ifftshift
from PIL import Image
def cross_correlate_2d(x, h):
    h = ifftshift(ifftshift(h, axes=0), axes=1)
    return ifft2(fft2(x) * np.conj(fft2(h))).real
load image as greyscale
x = np.array(Image.open("covid.png").convert("L")) / 255.
make kernels
h0 = np.zeros(x.shape, dtype=x.dtype)
h0[h0.shape[0]//2, h0.shape[1]//2] = 1
h1 = x.copy()
compute
out0 = cross_correlate_2d(x, h0)
out1 = cross_correlate_2d(x, h1)
plot
plt.imshow(out0); plt.xticks([]); plt.yticks([]); plt.show()
plt.imshow(out1); plt.xticks([]); plt.yticks([]); plt.show()

second example, same imports:

x = np.array(Image.open("covid_target.png"  ).convert("L")) / 255.
h = np.array(Image.open("covid_template.png").convert("L")) / 255.
blank regions default to 1, undo that
x[x==1] = 0
h[h==1] = 0
out = cross_correlate_2d(x, h)
plt.imshow(out, cmap='turbo'); plt.xticks([]); plt.yticks([]); plt.show()