Image FFT radial component reduction

Question

I'm using MATLAB's fft2 to generate the fourier power spectrum of an image. Naturally the spectrum is two dimensional as well, showing the direction of the direction and the frequency of the spectral components. I'd like to get rid of the directional component and just look at the frequency, i.e. reduce the dimension of the spectrum by one.

The way I see it, one can either:

• sum over a ring containing all pixels with $$ r_1 \leq \sqrt{x^2 + y^2} \leq r_2 $$ (which would make it tricky to select the correct boundary radii to really count the pixels that correspond to a certain frequency)
• or one can convert the image to a polar coordinate system and sum over the axis corresponding to the radial component of the matrix.

I wonder which method would be preferrable and for what reasons and what normalization I would have to make (of course, for the first method, there would be much less pixels contained in the ring element for lower frequencies than for higher frequencies and I guess some similar scaling problem will occur for the second method)

Answer 1

The two methods are equivalent over discrete images and in both of them it is "...tricky to select the correct boundary radii to really count the pixels that correspond to a certain frequency".

Whether "selecting pixels by distance" or converting from rectangular to polar coordinates, there will always be the problem of what do you do when you have to "land" between pixels. That is, when you have to "land" at an image site, within the acquired data range but between known values.

The answer here is of course "interpolation". Whether you are trying to select all harmonics that are $\approx 14.48765\ldots$ pixels from the centre of the image or you are converting $(x,y)$ pixel $(15,5)$ to its $(\rho, \theta)$ corresponding location of $(\approx 14.41479\ldots , \approx 81.36745\ldots)$, you still have the problem of "landing" between pixels.

You could of course $round(\cdot)$ everything, which would result in a nearest neighbour interpolation but that doesn't look very attractive. You are likely to have to use more complex interpolation schemes.

By the way, by convention, in digital image processing, the two dimensional discrete fourier transform's output is shifted by $\frac{N}{2}$ in both dimensions in order to bring the low frequencies to the centre of the "image". Therefore, you are likely to have to use fftshift too.

As far as the last part of your question is concerned, there should be no need for "scaling" if float or double data types are used, both to calculate the spectrum and its sum over different frequencies. If you are asking what sort of scaling you should do in order to be able to display the image (?), then, please note that you may have to use a logarithmic transformation on the spectrum values before you convert the image to something like unit8 to prepare it for display. Happy to clarify this later part more if more information emerges.

Hope this helps.

Answer 2

A change of coordinate is always a tricky part when considering discrete signals since the transformation is not bijective on finite domains.

In your case, I would try the second method and the function interp2() in MATLAB. Of course depending on the process you have to apply on your image, you might had artifacts but it should do the job. If you need to come back in the real space, it might severely reduce the quality of your image.

Otherwise, I don't know if this could help you but usually when considering the 2 dimensional Fourier transform, especially when polar coordinates are involved, you might want to check the Radon transform and the Fourier slice theorem. The first one is basically a line integral over the radius summed over the angles and the second one is the explicit link between a Fourier transform in 1 and 2 dimensions. This two mathematical tools are linked and might be useful for polar Fourier transformations.