Faster way of getting 2D frequency amplitudes than DFT?
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.