Is there an analogue to the 2D DFT that is rotation equivariant?

Question

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]])

Answer 1

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.