1
$\begingroup$

Given input images x and h of same size, greyscale 2D and potentially complex, we seek circular cross_correlate_2d(x, h) done with fft (so no padding) in image processing context.

I offer 1000 bounty to an answer that shows it's not exactly as described here. If it is exactly as described, then answers may state as such, but no bounty.

Bounty will be awarded to a differing cross_correlate_2d that passes below validation.

Answers that seek to argue against below validation by explaining how alternative definitions are more useful for general image processing are welcome. They are also subject to the bounty, but not guaranteed as it'll depend on whether I'm convinced.

Like checking any other source, this is a verification/clarification question.


Validation description

Following Cross-correlation Wikipedia on image processing, we seek similarity of x with h:

such that out[i, j] is strength of correlation of x with h being at i, j. Hence, if h = x, then peak will be at center of out, meaning the input aligns with itself when it fully overlaps itself, just as in full (padded) cross-correlation. Here's scipy, where I added imaginary WGN to an image:

The red circle is drawn around where(out == max(out)). More generally,

that's with a flipped false positive. With unflipped but conjugated, and real WGN:

Code to reproduce is available at Github.

Validation instruction

Change cross_correlate_2d below such that there's no AssertionError. No trickery, e.g. constant offset so technically it's different. It tests that above cases are reproduced, with images replaced by complex WGN.

def cross_correlate_2d(x, h):
    h = ifftshift(ifftshift(h, axes=0), axes=1)
    return ifft2(fft2(x) * np.conj(fft2(h)))

##############################################################################
import numpy as np
from numpy.fft import fft2, ifft2, ifftshift

def crand(*s):
    return np.random.randn(*s) + 1j*np.random.randn(*s)

for M in (64, 65, 99):
    for N in (64, 65, 99):
        # case 1 -------------------------------------------------------------
        x = crand(M, N)

        o = cross_correlate_2d(x, x)
        hmax, wmax = np.where(abs(o) == abs(o).max())
        assert hmax == M//2, (hmax, M//2)
        assert wmax == N//2, (wmax, N//2)

        # case 2 -------------------------------------------------------------
        x = np.zeros((M, N), dtype='complex128')
        h = np.zeros((M, N), dtype='complex128')

        hctr, wctr = M//8, N//8
        hsize, wsize = hctr*2, wctr*2
        target_loc = slice(0, hsize), slice(0, wsize)
        false_positive_loc = slice(M//2, M//2 + hsize), slice(-wsize, None)
        target = crand(hsize, wsize)

        x[target_loc] = target
        x[false_positive_loc] = target[::-1, ::-1]
        h[M//2 - hctr:M//2 + hctr, N//2 - wctr:N//2 + wctr] = target

        o = cross_correlate_2d(x, h)
        hmax, wmax = np.where(abs(o) == abs(o).max())
        assert hmax == hctr, (hmax, hctr)
        assert wmax == wctr, (wmax, wctr)

        # case 3 -------------------------------------------------------------
        x[false_positive_loc] = np.conj(target)

        o = cross_correlate_2d(x, h)
        hmax, wmax = np.where(abs(o) == abs(o).max())
        assert hmax == hctr, (hmax, hctr)
        assert wmax == wctr, (wmax, wctr)
$\endgroup$

2 Answers 2

6
$\begingroup$

For the question's approach of shifting the h input image spatially using ifftshift before frequency-domain complex-conjugation, a possible motivation is to consider that the circular cross-correlation output is also an image, that the input images have anchors at their (non-periodic) centers, and that for identical input images, to indicate no displacement between the images, a peak of the circular cross-correlation output should appear at the anchor point of the reference input image h. Note that this answer assumes periodic input images, no padding, throughout, except for the notion of a non-periodic center of an image which is indirectly implied in the question.

Let's first have a look at what the question's shifting approach does in the even and odd period cases when the input images are identical impulses or impulse trains if considering the periodicity, compared to no ifftshift which would result in circular cross-correlation peaking at the first output sample. It suffices to look at the 1-dimensional case:

import numpy as np
from numpy.fft import fft, ifft, fftshift, ifftshift
x = [1, 0, 0, 0]
h = [1, 0, 0, 0]
print(np.real(ifft(fft(x) * np.conj(fft(ifftshift(h))))).astype(int))
[0 0 1 0]

That was a shift of N/2 for the even period N case, compared to [1 0 0 0].

x = [1, 0, 0, 0, 0]
h = [1, 0, 0, 0, 0]
print(np.real(ifft(fft(x) * np.conj(fft(ifftshift(h))))).astype(int))
[0 0 1 0 0]

That was a shift of (N - 1)/2 for the odd period N case, compared to [1 0 0 0 0].

Note that identical behavior would be obtained by removing the ifftshift and by fftshift of the circular cross-correlation output:

fftshift(ifft(fft(x) * np.conj(fft(h))))

It may be useful to look at the even and odd period cases in combination with 0-based and ½-based indexing schemes, both commonly used in image processing. Spoiler: It will turn out that the indexing scheme doesn't matter, so for simplicity you may wish to look at just the 0-based indexing scheme in the following:

enter image description here

The pixel indexes are shown as red dots, and the shared anchor index of the input images is shown as a cross. Now we test the question's proposed circular cross-correlation shifting scheme, on whether the peak of the circular cross-correlation output for identical inputs will be at the anchor index or not:

  • With 0-based indexing and an even period 4 the output of the question's circular cross-correlation method peaks at index 2, which is unfortunately not at the anchor index 1.5. ❌
  • With ½-based indexing and an even period 4 the output of the question's circular cross-correlation method peaks at index 2.5, which is unfortunately not at the anchor index 2. ❌
  • With 0-based indexing and an odd period 5 the output of the question's circular cross-correlation method peaks at index 2, which is at the anchor location 2 as desired. ✅
  • With ½-based indexing and an odd period 5 the output of the question's circular cross-correlation method peaks at index 2.5, which is at the anchor location 2.5 as desired. ✅

In summary, the question's approach of using ifftshift results in a small indexing error of half a sampling period, for the even period case. This problem cannot be solved by any non-interpolating shifting scheme. Cross-correlation of band-limited signals is identically bandlimited so complex exponential modulation in frequency domain could be used to do an interpolative shift assuming that the input images are band-limited horizontally and vertically below half the horizontal and vertical spatial sampling frequencies.

An alternative, and I would think the standard approach, is to do no shifting at all in the circular cross-correlation calculation, as in this answer. Then a peak at the first sample in the circular cross-correlation output is quite intuitively, it being at index zero, understood as indicating zero displacement between the input images. This is in agreement with the mathematical definition of circular cross-correlation for finite discrete functions:

$$(f \star g)[n]\ \triangleq \sum_{m=0}^{N-1} \overline{f[m]} g[(m+n)_{\text{mod}~N}]$$

As an extension, an anchor coordinate based correction to the calculated argmax of circular cross-correlation can be added if displaced anchor coordinates are needed. While also allowing non-centered anchors, this however doesn't enable as intuitive (albeit slightly erroneous) visualization of the full circular cross-correlation output as the question's approach.

Another motivation behind the question's approach may be to detect displacements using local methods that do not handle wraparound because they were designed for non-periodic data. The shift by the question's approach approximately symmetrically maximizes the amount of displacement that can be detected by such methods.

After reading the original poster's (OP's) answer, I understand that a motivation to do it as OP suggests may be to maintain compatibility when extending a library cross-correlation function to support circular padding in addition to any existing padding modes. There, if the existing implementation assumes a certain origin and support of the inputs (machine learning libraries typically provide convolutional neural network functionality with "centered" convolution kernels), it would be natural and of least surprise to a user to maintain those for the circular padding mode. Documenting this does not necessitate redefining what is meant, mathematically, by circular cross-correlation.

I can't think of any actual image processing application for circular cross-correlation with equal size x and h, other than acceleration of non-circular cross-correlation, which to me is ruled out by the question's statement "we seek circular cross_correlate_2d(x, h)".

$\endgroup$
16
  • 2
    $\begingroup$ The whole thing is about anchor points (Basically the origin). While the intuition of the OP is correct, we like functions centered, the math of DFT is having the origin at the first sample. Since all operators are shift invariant, it doesn't really matter where you shift. +1. $\endgroup$
    – Royi
    May 9, 2023 at 9:38
  • 1
    $\begingroup$ @OverLordGoldDragon, I really have no idea what you're talking about. I only know that if we talk on the valid part of the correlation / convolution, there is no problem to achieve it without any need of the shifts. It is all a game of indices. I thought Olli answers is good from your side. If it is not, I hope it will be for others. $\endgroup$
    – Royi
    May 9, 2023 at 19:25
  • 1
    $\begingroup$ @OverLordGoldDragon, I don't want to get into things. We agreed not to agree. I think there is no problem to apply any kind of cross correlation on frequency domain without any need for shifts. You think differently. It is OK. I have never ran into issues with my approach, yours works well for you. All good. Concentrate on creating more usable information. It seems you have plenty to share. $\endgroup$
    – Royi
    May 10, 2023 at 7:12
  • 1
    $\begingroup$ @Jdip Indeed I should've said "can read and don't run from drama". $\endgroup$ May 10, 2023 at 17:13
  • 3
    $\begingroup$ @OverLordGoldDragon, I suggest you stop evaluate the skills of other to have conversations. Specifically, It might be the issue with my English, but sometimes I am not sure I get the little details of what you mean. All I said, and I stand behind it, no need for fftshift() in order to implement cross correlation in frequency domain. Even for 2D. You said I'm wrong, I offered you open a question. That's it. Let's focus on creating knowledge. $\endgroup$
    – Royi
    May 10, 2023 at 17:45
0
$\begingroup$

OP's approach follows three basic premises:

  1. Output at i, j corresponds to similarity of kernel with input, with kernel centered at i, j, so that the input and output axes are perfectly aligned, as would be the case in 1D continuous time with inputs centered at $t=0$.
  2. Centering indexing conventions are followed. The "center" can either be at index 0, or its fftshift which is N//2. Usual image inputs are never centered at index 0, so N//2.
  3. The images, as-is, are ground truth. Pixels are discrete, so there's no such thing as "fractional index".

As will be shown, all standard filtering follows these premises.

Summary

  1. Only OP's approach yields valid images as outputs, for valid image inputs
  2. Only OP's approach yields outputs that can be used as features for vast majority of cases: distance algorithms, edge detection, contour tracing, etc
  3. Only OP's approach yields outputs subject to subsequent filtering, e.g. after nonlinearity
  4. Only OP's approach has the property that shifting the target in the image will shift the peak in the output by the same amount, for all shifts
  5. Only OP's approach can utilize DFT properties to yield commonly expected behavior, like time reversal
  6. Only OP's approach is consistent with filtering in vast majority of other contexts, including basic spectral manipulation (e.g. lowpass), STFT and CWT, and all of machine learning.

The following disclaimer is specific to this network's regular users: if any of these points are disputed, I'll only reply in comments if the entire answer has been read first.

Minimal version

When doing usual cross-correlation with same output mode, do you input the left or right image?

Right is what non-OP's approaches do, with the conj(fft2(h)) variant being within 1-sample shift:

import numpy as np; from numpy.fft import fft2, ifft2, fftshift
import scipy.signal
x = np.random.randn(30, 30) + 1j*np.random.randn(30, 30)
h = np.random.randn(30, 30) + 1j*np.random.randn(30, 30)

xs = fftshift(fftshift(x, axes=0), axes=1)
assert np.allclose(ifft2(fft2(x) * fft2(np.conj(h[::-1, ::-1]))),
                   scipy.signal.correlate2d(xs, h, mode='same', boundary='circular'))

Standard filtering

  1. Low-pass, high-pass, band-pass, and much of rudimentary spectral manipulations are mostly done with zero-phase filters, so that: (a) filtering doesn't shift/delay the signal, or misalign its spectral contents; (b) the filter isn't asymmetric so that left or right side of windowed input isn't favored. Zero-phase means real-valued in Fourier; this guarantees that the filter is centered at either 0 or N//2.

  2. CWT, STFT, and all time-frequency analysis is done per premise 1.

  3. Convolution and cross-correlation with 'same' output mode, which ensures that output size equals input size, unpad in a way that ensures premise 1, in every major implementation (NumPy, scipy, MATLAB, all of machine learning). The unpadding is done with surgical precision to ensure this: we unpad such that at each boundary edge, the kernel overlaps (is centered over) input as much as possible.

All of this follows an extremely basic objective of a time-to-time transform: alignment. Output at $t_0$ corresponds to input being modified around $t_0$. For an LTI system with impulse response $h$, the output at $t_0$ is a modification of input with $h$ centered over $x(t_0)$, such that $y(t)$ and $x(t)$ are directly comparable, and that if $h$ is the unit impulse, then $y(t) = x(t)$.

Standard filter formulation

In continuous time, virtually all filters will look like what's on the left rather than on the right:

It's done for same reason all of mathematics is done oriented around the origin, whenever it's the case: it's the reference coordinate. Even if we were to center elsewhere, where would it be for IIR filters? It'd be just as arbitrary as centering at $t=0$, and in the end we'd still have to align it with input to get what we want.

Yes, this is all very basic. And partly for same basic reasons, when we sample a time-domain filter in code, we get something that looks like this:

which is also what all major implementations will return windows as, and indeed this imports from scipy.signal.windows. Yet, they won't do filtering this way; all filtering is done zero-phase, 0-centered. For images, it means this MATLAB example would be:

Onto the premises:

P1: Input-output alignment

This agrees with the origin definition of cross-correlation:

$$ (x \star h)(\tau) = \int_{-\infty}^{\infty} x(t) \overline{h(t - \tau)} dt \tag{1} $$

Output at shift $\tau$ corresponds to the kernel, $h$, being shifted by $\tau$. Combined with the standard of centering $h$ over $t=0$, this becomes "$h$ being centered over $\tau$".

P2: Center index is either 0 or N//2

All major library-generated filters, if visually centered, will be centered with respect to N//2. Scipy example:

from scipy.signal.windows import hann

for N in range(64, 129):
    for sym in (True, False):
        w = hann(N, sym=sym)
        assert w[N//2] == w.max()

For sym=True, it also passes for N//2 - 1, since by definition it means left half of array equals right, and putting the second peak at N//2 + 1 would be farther off center. sym=False, the intended option for windowed FFT, only passes for N//2. ifftshift moves the peak to index 0.

More pertinently, we're working with FFT, which is centered about 0. fftshift moves this center to N//2, which is why ifftshift moves to 0. These conventions enforce DFT properties; if centers were anywhere else, we'd lose, for example, $x[n] = x[-n] \Leftrightarrow \text{zero phase}$, where negative indexing is circular per DFT's definition. Most properties would lose duality without intermediate and domain-specific shift manipulations. Conversely, if a time-domain manipulation doesn't follow said centering with standard FFT, it won't affect FFT as intended: flip(x) is not time-reversal.

P3: Discrete is ground-truth

Suppose this is false. Then every major library and centering conventions are invalidated.

Of course, we can insist on being centered over fractions, which would require bandlimited half-index-shift interpolation for the even case. This is a valid, yet completely separate sentiment, and completely separate operation that is regardless forced to build (operate) upon a non-fractional indexed output.

Other useful properties

OP's proposal is the only one that satisfies

$$ \tau + 1 > \tau $$

For ifft(fft(x) * conj(fft(h))), index N//2 and N//2 + 1 are apart by N//2: there's a jump-discontinuity, and indexing is no longer ordinal. This invalidates discrete-continuous domain correspondence, as continuous domains (in any signal processing) are by definition ordinal.

This discontinuity invalidates all spatially local operations upon the output. Any further filtering (e.g. after an intermediate operation), distance algorithms, edge detection, contour tracing - gone. OP's approach preserves spatial contiguity.

Re: Olli's answer

In summary, the question's approach of using ifftshift results in a small indexing error of half a sampling period, for the even period case. This problem cannot be solved by any non-interpolating shifting scheme.

i.e., if it's considered a problem, OP's approach suffers the same as the other alternatives in question.

Then a peak at the first sample in the circular cross-correlation output is quite intuitively, it being at index zero, understood as indicating zero displacement between the input images

Quite unintuitively, index N//2 - 1 is displacement by N//2 - 1, and index N//2 is displacement by -N//2. Shifting target in image by T left or right won't always shift output's peak by T left or right - likewise for up and down.

Another motivation behind the question's approach may be to detect displacements using local methods that do not handle wraparound because they were designed for non-periodic data.

In other words: if your inputs are valid images and you want a valid image output, use OP's approach. Indeed this refers to spatial contiguity.

Re: other approaches

We contrast OP's approach with others' by exploring latter explicitly. This section primarily concerns advocacy for fft(flip(conj(h)))-based approach with the goal of being exact, which was described in this answer by @Gillispie.

Following the answer's setup, Gillispie seeks to show that, in order to obtain equivalence between convolution and cross-correlation, we must time reverse correctly. Regardless of whether or not that point is achieved correctly, there's no mention of the critical practical fact that input images are visually centered, and that not accounting for it will run into all the problems described in this answer. Turns out, it's also not achieved correctly.

To be safe, we examine the relevant portions that show the objective is to be exact when it comes to circular time reversal:

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]$.

This completely correct observation is, unfortunately, spinned into the opposite conclusion. We re-iterate the motivation of equivalence, also of correspondence between time and freq domains:

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!

So the goal is $\mathcal{F}^{-1} (\overline{X(\omega)}) = \overline{x(-t)}$, in discrete. If readers were careful with this answer, they already see the problem. If not, let's be explicit: "$x(-t)$" assumes the existence of $x(0)$. Yet, $x(0)$ is not at index 0 - it is at index N//2, as Gillispie is aware by stating that flip is the solution. So, we seek $x[n] \rightarrow x[-n]$, with $x[0]$ at N//2. This is not achieved:

It's not achieved for any other choice of $x[0]$ centering either - not N//2 + 1, not N//2 - 1, not 0, or anything else. It is not achieved with Jdip's approach either. It is, ironically, achieved only with OP's approach:

A separate but related problem is described by Olli:

Gillispie is right to emphasize the concern, as it has numerous implications. For a zero-phase complex filter, for example, conjugation is identical to time reversal, but only if it's done right.

The motivation for OP's answer to the linked Q&A was two-fold: to correct Gillispie on his own goal, and to pose the most useful variant of cross-correlation for general image processing. After reading this answer of mine, let the reader judge whether OP succeeded.

$\endgroup$
16
  • $\begingroup$ All the properties you mention are habits and make things easier for us to work. They have no need in the Math form and actually are not enforced. For instance same in MATLAB just work by where to start / stop the iteration (By the way, filter() of MATLAB is just the shifted indices). There is no need to do things in the way you present above. No need to do fftshift() in order to apply correctly cross correlation in frequency domain. $\endgroup$
    – Royi
    May 10, 2023 at 8:09
  • 1
    $\begingroup$ @OverLordGoldDragon. You write Output at i, j corresponds to similarity of kernel with input, with kernel centered at i, j. This has no advantage over any other convention. You build system of constraints which has not advantage on any other set of coherent convention. Moreover, you say, accept my constraints and then try prove me wrong. The only issue is there is no necessity to follow those conventions. As I wrote, when you use the DFT, it is much easier to handle indices if you accept the convention of the origins being on the top left (For 2D data). $\endgroup$
    – Royi
    May 10, 2023 at 17:52
  • 1
    $\begingroup$ @OverLordGoldDragon, There is a simpler criterion which much more sense. Use the top left as the anchor. Apply the the property of the DFT: Everything is periodic. If you want, for some reason, to have the top left at a specific place, shift it. When you shift keep in mind you're on a cyclic space, so when you shift, add the periodic samples. same and full has nothing to do with this. full is the only valid convolution operator in the sense it keeps the commutative property. The same is just cropping the beginning or the tail. It is actually just a convention MATLAB created. $\endgroup$
    – Royi
    May 12, 2023 at 7:22
  • 1
    $\begingroup$ As I wrote, MATLAB itself, with filter() crops to the other side of conv() with same. Which is correct? None. Cross correlation is the normalized inner product in a vector space. It has nothing to do with shifted indices. You're mixing conventions of 2D image filtering with Math conventions. There are many ways to make things work as our intuition. There is nothing special in your code, people have been doing this or the other way (As I described above) for years. $\endgroup$
    – Royi
    May 12, 2023 at 7:25
  • 2
    $\begingroup$ Try this in MATLAB: vA = [1, 2, 3]; vB = [4, 5]; conv(vA, vB, 'full'), conv(vA, vB, 'same'), filter(vB, 1, vA). conv() with same just crops the beginning (1 sample as vB is length 2). The filter() crops the end (1 sample, same reasoning). Neither have any knowledge whether the filter is symmetric or not. $\endgroup$
    – Royi
    May 12, 2023 at 7:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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