Faster way of getting 2D frequency amplitudes than DFT?

Question

I'm making a small program which gets the DFT of an image to get a general idea of the image's overall orientation. It does this by rotating a line in a radar sweep type pattern, keeping track of which angle has the highest sum of frequency amplitudes. (It actually only does $0-180^\circ$ though because the testing line is symmetrical across the origin)

Anyways, I was wondering, is there a less computationally expensive way to get frequency amplitudes of a 2D image than doing a 2D DFT?

I know that it is separable so can do better than the naive $\mathcal O(N^4)$ operation - I think it becomes $\mathcal O(N^3)$. I also know that the FFT can do a 1D DFT in $\mathcal O(N \log N)$ instead of $\mathcal O(N^2)$. Lastly, I know that it's massively parallelizable so could be (massively) multithreaded on either CPU or GPU.

Those things aside, are there any less computationally expensive ways to get frequency amplitudes of a 2D image than to do a full 2D DFT?

Answer 1

Note that the computational complexity of a 2D FFT would be $\mathcal O(N^2\log N)$ instead of $\mathcal O(N^4)$ for the naive 2D DFT.

That said, since you mention that you want to obtain spectrum information along beams of a radar-like pattern, it sounds like you are really more interested in computing the 1D FFT of each such beam. In that case, for $M$ beams the computational complexity would be $\mathcal O(MN \log N)$.

As an illustration, lets process the following figure:

from which we shall extract a set of $M$ beams with the following sample matlab code:

% Number of beams to extract
M = 200;

% Convert the input image "img" to polar coordinates
c = size(img)/2+1;
N = max(size(img));
angles = 2*pi*[0:M-1]/M;
radius = [0:floor(N/2)-1];
imgpolar = zeros(length(radius), length(angles));
for ii=1:length(radius)
  xi = min(max(1, floor(c(2)+radius(ii)*cos(angles))), size(img,2));
  yi = min(max(1, floor(c(1)-radius(ii)*sin(angles))), size(img,1));
  imgpolar(ii,:) = img(sub2ind(size(img), yi, xi));
end
% Compute the FFT for each beam angle
ImgFD = fft(imgpolar,[],1);

figure(1);
freqs = [0:size(ImgFD,1)-1]/size(ImgFD,1);
surf(angles, freqs, 10*log10(abs(ImgFD)+1), 'EdgeColor', 'None');
view(2);
colormap("gray");
xlabel('Beam angle (radians)');
ylabel('Normalized frequency');

to yield:

which can be collapsed to the sum of amplitudes as a function of the beam angles to give:

SumAmplitudes = sum(abs(ImgFD),1);
figure(2);
hold off; plot(angles, 10*log10(SumAmplitudes+1));
xlabel('beam angle (radians)');
ylabel('Sum of amplitudes (dB)');

As a side note, if you can use sum of squared amplitudes along those beam (instead of sum of amplitudes), then you can do it directly in spatial domain thanks to Parseval's theorem (which would bring the computational complexity down to $\mathcal O(MN)$, dominated by the conversion to polar coordinates). The equivalence (for the sum of squared amplitudes) can be seen using the following code:

% Compare the result of square amplitude summation in the frequency domain vs spatial domain
SumFD = sum(abs(ImgFD).^2,1)/size(ImgFD,1);
SumSD = sum(abs(imgpolar).^2,1);

figure(3);
hold off; plot(angles, 10*log10(SumFD+1), 'b');
hold on;  plot(angles, 10*log10(SumSD+1), 'r:');
xlabel('beam angle (radians)');
ylabel('Sum of squared amplitudes (dB)');
legend('Frequency domain', 'Spatial domain', "location", "southwest");

Notice the overlap of the curves computed in the spatial and frequency domains:

Update:

If you are in fact computing beams of a radar-like pattern in a 2D frequency plot as you seem to suggest in this other post, then you best bet comes back to doing a 2D FFT which would be order $\mathcal O(N^2 \log N)$. You could then perform the conversion to polar form for the sum of frequency coefficients which would add a small $\mathcal O(N^2)$, so the result is still dominated by the 2D FFT.