4
$\begingroup$

I'd like to know if there is some general Fourier transform or other signal processing algorithm, such as a discrete wavelet transform, that is rotation equivariant.

Rotational equivariance of a function means that $g:\mathbb R^{n\times n}\to \mathbb R^{n\times n}$ means that if $R$ is a rotation, then $g(RX) = Rg(X)$ - that is, rotating the input to the function is the same as rotating the output of the function applied to the unrotated input.

The reason I would like to have a transform $\mathcal F$ that is equivariant is so that I can apply an equivariant function $g$ in the frequency domain before transforming back using $\mathcal F^{-1}$ have have the whole process be equivariant - that is:

$$\mathcal F ^{-1}\circ g\circ \mathcal F \circ RX = R\circ \mathcal F ^{-1}\circ g\circ \mathcal F X$$

The closest I have gotten so far is by padding $X$, as described in this post. I demonstrate that this roughly works on a small "image" below, but the amount of padding required for larger images becomes intractable.

import numpy as np

# generate some random data
np.random.seed(1)
X = np.random.normal(size=(3, 3))
X.round(2)
Out[2]: 
array([[ 1.62, -0.61, -0.53],
       [-1.07,  0.87, -2.3 ],
       [ 1.74, -0.76,  0.32]])

# rotate X 90 degrees
X_rt = np.rot90(X)
X_rt.round(2)
Out[3]: 
array([[-0.53, -2.3 ,  0.32],
       [-0.61,  0.87, -0.76],
       [ 1.62, -1.07,  1.74]])

# calculate the transform of X and the transform of rotated X
F = np.fft.fft2(X)
F.round(2)
Out[4]: 
array([[-0.72+0.j  ,  3.81-1.73j,  3.81+1.73j],
       [ 1.09+3.3j , -1.8 +2.99j,  4.57+1.03j],
       [ 1.09-3.3j ,  4.57-1.03j, -1.8 -2.99j]])

F_rt = np.fft.fft2(X_rt)
F_rt.round(2)
Out[5]: 
array([[-0.72+0.j  ,  1.09+3.3j ,  1.09-3.3j ],
       [-3.4 +2.43j, -3.18+3.44j,  3.48-0.06j],
       [-3.4 -2.43j,  3.48+0.06j, -3.18-3.44j]])

# if transform was rotation equivariant, then this would be the zero matrix:
(np.rot90(F) - F_rt).round(2)
Out[6]: 
array([[ 4.53+1.73j,  3.49-2.27j, -2.88+0.32j],
       [ 7.21-4.16j,  1.39-0.46j,  1.09-0.97j],
       [ 2.68+2.43j, -2.4 +3.24j,  4.27+0.14j]])

# pad X and repeat
pad = 500
X = np.pad(X, pad_width=pad)

# rotate X 90 degrees
X_rt = np.rot90(X)

# calculate the transform of X and the rotated transform of X
F = np.fft.fft2(X)
F_rt = np.fft.fft2(X_rt)

# remove padding
F = F[pad:-pad, pad:-pad]
F_rt = F_rt[pad:-pad, pad:-pad]

# approximately equivariant
(np.rot90(F) - F_rt).round(2)
Out[11]: 
array([[-0.-0.04j,  0.+0.04j, -0.-0.04j],
       [-0.-0.01j,  0.+0.01j,  0.-0.01j],
       [ 0.+0.07j, -0.-0.07j,  0.+0.07j]])
$\endgroup$
4
  • $\begingroup$ You cannot rotate the frequency domain like that. The origin is the top-left element. Any transformation you apply should maintain that value in its place. $\endgroup$ Commented Nov 5, 2022 at 0:03
  • $\begingroup$ That said, F^-1 g F R X == R F^-1 g F X if g is rotation equivariant, this equality does not require applying the rotation in between the two Fourier transforms. $\endgroup$ Commented Nov 5, 2022 at 0:06
  • $\begingroup$ @CrisLuengo Are you sure about that? I have added an example where this does not seem to be true - it is possible I have made a mistake, though $\endgroup$ Commented Nov 5, 2022 at 17:21
  • $\begingroup$ Again, in the frequency domain the origin is at the top-left element, multiplying by your W is not rotation equivariant in the frequency domain. $\endgroup$ Commented Nov 5, 2022 at 18:35

1 Answer 1

3
$\begingroup$

The Fourier transform is rotation equivariant. Its discrete counterpart is only equivariant over rotations of 90 degrees, because of sampling issues and the effects of the sampling window. However, you need to realize that the origin of the DFT (which you compute with the FFT algorithm) is in the top-left element. That is, in the discrete frequency domain, the rotation needs to be around the top-left element, and this rotation is complicated by the fact that the discrete Fourier domain is periodic.

But note that your expression $\mathcal F ^{-1}\circ g\circ \mathcal F \circ RX = R\circ \mathcal F ^{-1}\circ g\circ \mathcal F X$ doesn't require the Fourier transform to be equivariant, it requires the operation $F ^{-1}\circ g\circ \mathcal F$ to be equivariant. If $g$ is a multiplication, as you had in the now deleted section of your question, then the operation $F^{-1}\circ g\circ \mathcal F$ is the convolution with $F(g)$. A convolution is equivariant to rotation if the convolution kernel is invariant to rotation (i.e. circularly symmetric). For your case of rotations over 90 degrees, the kernel only needs to be invariant to 90 degree rotations, so a square structuring element will suffice.

Let's do the experiment with a Gaussian convolution kernel. In the frequency domain, this is a Gaussian as well. I'll apply it to an actual image, which is more interesting than random data:

import numpy as np
import imageio.v2 as imageio

img = imageio.imread("cameraman.tif")

f = np.fft.fftfreq(img.shape[0])
g = np.exp(-0.5 * f**2 / 0.1**2)
g = g[:, None] * g[None, :]

out_1 = np.rot90(np.fft.ifft2(g * np.fft.fft2(img)))
out_2 = np.fft.ifft2(g * np.fft.fft2(np.rot90(img)))
np.allclose(out_1, out_2)

In fact, since the Gaussian is rotation invariant, we can do the same experiment with an arbitrary rotation. The two results will be identical away from the image edges, where things get muddy.

import diplib as dip

def rotation(img):
    phi = 0.3
    img = dip.Image(img)
    out = dip.Rotation2D(img, phi)
    out.Crop(img.Sizes())
    return np.asarray(out)

out_1 = rotation(np.fft.ifft2(g * np.fft.fft2(img)))
out_2 = np.fft.ifft2(g * np.fft.fft2(rotation(img)))

Here is out_1, out_2 and 10 * np.abs(out_1 - out_2). I've multiplied the difference by 10 to demonstrate that the results are identical away from the image edges.

two identical images and an image that is black in the central part

$\endgroup$
12
  • $\begingroup$ This is very helpful. Just a couple of follow ups: 1. "then the operation $\mathcal F^{−1}\circ g \circ\mathcal F$ is the convolution with $\mathcal F\circ g$". I interpret this to mean that $\mathcal F^{−1}\circ g \circ\mathcal F(x) = \mathcal F(x) * g(x)$, where $*$ is convolution - could you correct me or elaborate as to why this is true? 2. "convolution is equivariant to rotation if the convolution kernel is invariant to rotation" - Could you direct me to where I might look to find a proof of this? So far, I only see ML results $\endgroup$ Commented Nov 6, 2022 at 6:27
  • $\begingroup$ 3. I think that the origin being the top left pixel is because of the ordering of the sample frequencies (i.e., np.fft.fftfreq(n)==[0, 1, ..., (n-1)/2, -(n-1)/2,..., -1] / n ). I am trying to connect this to the DFT definition: $$DFT(x)_k = \sum_{j=0}^{N-1}x_j\left[\cos 2\pi j\frac kN - i\sin2\pi j\frac kN\right]$$ From this, it seems like the sample frequencies should be $\frac kN,k\in0,...,N-1$ - why is this not the case? $\endgroup$ Commented Nov 6, 2022 at 6:42
  • $\begingroup$ 4. A comment - I am surprised you were able to achieve rotational equivariance, since the discretized version of your kernel was not exactly circularly symmetric ((np.abs(np.rot90(g) - g) / np.abs(g)).mean() gives 0.098). This difference goes to 0 as the standard deviation of the Gaussian is increased from 0.1. Also sort of surprising that the kernel had such little effect on the resulting image, and basically just turned it green and purple - decreasing the standard deviation by a factor of 10 turns the cameraman into a purple blob $\endgroup$ Commented Nov 6, 2022 at 7:26
  • $\begingroup$ @JacobHelwig (1) I’m going to express the property this way, it’s easier notation for me: $\mathcal{F}^{-1}(g \cdot \mathcal{F}(f)) = \mathcal{F}(g) * f$. This is the convolution property of the Fourier transform, and is the reason that the Fourier transform is useful for analyzing linear systems. (2) The proof for equivariance should be in any signals and systems book. I learned from Oppenheim, Willsky and Young at the time. It is easily provable for the continuous domain using the convolution property of the Fourier transform. For the discrete case, rotation is always awkward. $\endgroup$ Commented Nov 6, 2022 at 13:48
  • $\begingroup$ @JacobHelwig (3) Add to that definition of the DFT that $k in [0, N-1]$. But the DFT is periodic, such that $DFT(x)_k = DFT(x)_{k+i N}$ for any integer i. People are just used to looking at the range from -N/2 to +N/2. For a real signal, the Fourier transform is conjugate symmetric, putting the zero in the middle makes it prettier. But the DFT really only reports positive frequencies, such that k matches the index into the array. $\endgroup$ Commented Nov 6, 2022 at 13:57

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.