Spatial frequency of pixels in an FFT transformed image

[lams]

This question is about converting from spatial to frequency domains when performing an FFT on 2D image data. I suspect my problem has a painfully obvious solution that I'm just not seeing.

If I have some image data g(x,y) where I know the pixel resolution, say 256x256 pels at a resolution of 0.1 μm/pel, then I know my image is 25.6μm across in both the X and Y axes.

When I take the FFT of this image and rotate quadrants to centre the low-frequency components (using, for example, MATLAB's fft() and fftshift() functions), how can I obtain the frequency associated with each pixel in the image? I'm quite confused about this.

The temptation is to simply compute the radial distance of each pixel relative to the centre of the spatial image g(x,y) in units of μm, and then take the reciprocal of these values. Obviously this would give me an infinity at the origin at the centre of the spatial image before rotating quadrants to centre low-frequency components. This seems wrong, but I might be wrong about being wrong.

I need to apply a low-pass filter to this image data at a specific cutoff frequency, and I'm not an optics whiz so some help would be great.

 

[Blibbler]

An idea...

For the low pass filter cutoff frequency, think about Nyquist and Shannon.

If this were an audio problem and I was looking to get rid of clicks and aliasing I would set the lpf to be just under half the sample rate and would make it a shelf filter - 8th order n-pole cascade with a little positive feedback to keep the 3db point (-93db if you're in 16 bit and counting down from 0db as your normalised signal input level) top end of the response curve below half the sample rate. Keep the order of the filter as a multiple of 4 to prevent phase shift from skewing your results.

This is an agricultural approach but it would probably work without you having to scramble your brains in the z plane.

 

[lams]

Thanks for your response. Unfortunately, my cutoff frequency relates to physical parameters I have no control over. I'm not just cleaning up signals. I really do need to be able to whack some spatial frequency axes onto my phase and magnitude plots.

For the low pass filter cutoff frequency, think about Nyquist and Shannon.

If this were an audio problem and I was looking to get rid of clicks and aliasing I would set the lpf to be just under half the sample rate and would make it a shelf filter - 8th order n-pole cascade with a little positive feedback to keep the 3db point (-93db if you're in 16 bit and counting down from 0db as your normalised signal input level) top end of the response curve below half the sample rate. Keep the order of the filter as a multiple of 4 to prevent phase shift from skewing your results.

This is an agricultural approach but it would probably work without you having to scramble your brains in the z plane.

 

[olivermsun]

If the FFT length is N=2^n and the spatial length of the record is L, then the FFT frequencies after you fftshift() are just (in Matlab code): [-N/2:N/2-1] / L .

 

[Dale]

If I have some image data g(x,y) where I know the pixel resolution, say 256x256 pels at a resolution of 0.1 μm/pel, then I know my image is 25.6μm across in both the X and Y axes.

If your image is 25.6 μm across then each "pixel" in the Fourier domain represents a certain number of cycles per 25.6 μm. The "pixel" at the origin is 0 cycles/25.6 μm = 0 m^-1 and the next "pixel" is 1 cycle/25.6 μm = 39062.5 m^-1 and so forth. Your highest Fourier "pixel" will be 128 cycles/25.6 μm = 5E6 m^-1.

 

[AlephZero]

For image data, it may be easier to think in terms of wavelength than frequency,.

From the Nyquist criterion, the shortest wavelength (i.e. highest "frequency") in the FFT has a length of 2 pixels. From the image size you know the number of pixels per inch, meter, or whatever.

 

[lams]

@olivermsun/DaleSpam - thank you for clarifying this point — I am currently computing the radial position of each point on my image in units of pixels, and then simply dividing each radial position by the width of the image (it's a square image).

@DaleSpam - you say the "highest Fourier "pixel" will be 128 cycles/25.6 μm = 5E6 m^-1". Since 128 cycles is half of the total number of pixels, i.e. the span form the origin, then shouldn't the displacement used be 25.6μm/2 = 12.8μm? Also, when I'm off-axis I assume I should still be dividing the radial position by the same displacement. Is this correct?

Thank you both for the clarification.

 

[Dale]

@DaleSpam - you say the "highest Fourier "pixel" will be 128 cycles/25.6 μm = 5E6 m^-1". Since 128 cycles is half of the total number of pixels, i.e. the span form the origin, then shouldn't the displacement used be 25.6μm/2 = 12.8μm? Also, when I'm off-axis I assume I should still be dividing the radial position by the same displacement. Is this correct?

No, you need two pixels minimum to fit one cycle. So 128 cycles is the maximum number of cycles that can fit in a 256 pixel image. So the maximum wave number in the image is 128 cycles/25.6 μm = 5E6 m^-1. You would never divide by 12.8 μm nor any other distance, you always divide by the FOV = 25.6 μm.

You mention displacement a couple of times. Displacement is irrelevant, don't think in terms of displacement. The frequency domain is in cycles per m or cycles per FOV, depending on if you are using physical wavelengths or image-based wavelength measures. The "position" in the Fourier domain is unrelated to any displacement in the image domain. Instead, the position in the Fourier domain is the inverse of the wavelength of the Fourier transform of a delta function at that location.

 

[olivermsun]

I am currently computing the radial position of each point on my image in units of pixels, and then simply dividing each radial position by the width of the image (it's a square image).

Are you computing a 2-d FFT? If so, why are you using radial positions?

The 2-d Fourier transform will have 2 spatial frequency (wavenumber) axes. After sending the output of the FFT through fftshift(), then the zero frequencies will be near the center of the transformed image, but there is a slight offset because the frequencies are numbered from -N/2 to N/2-1 for an even-length transform, which is how you get the Matlab expression for the frequencies that I posted above.

 

[Dale]

Here is a little interactive widget I wrote in Wolfram's cdf format (http://www.wolfram.com/cdf-player/). It shows the image space result of a single "pixel" in the Fourier domain. Note that there is no "displacement" in the image domain since each "pixel" in the Fourier domain covers the entire image domain.

 

[lams]

@olivermsun - I'm applying a cutoff frequency. So I need to calculate the frequency associated with every point in my Fourier domain image and exclude pixels above the cutoff. Geometrically this gives me a circular unit step-function mask centred at the origin in the Fourier domain.

Are you computing a 2-d FFT? If so, why are you using radial positions?

The 2-d Fourier transform will have 2 spatial frequency (wavenumber) axes. After sending the output of the FFT through fftshift(), then the zero frequencies will be near the center of the transformed image, but there is a slight offset because the frequencies are numbered from -N/2 to N/2-1 for an even-length transform, which is how you get the Matlab expression for the frequencies that I posted above.

 

[lams]

No, you need two pixels minimum to fit one cycle. So 128 cycles is the maximum number of cycles that can fit in a 256 pixel image. So the maximum wave number in the image is 128 cycles/25.6 μm = 5E6 m^-1. You would never divide by 12.8 μm nor any other distance, you always divide by the FOV = 25.6 μm.

You mention displacement a couple of times. Displacement is irrelevant, don't think in terms of displacement. The frequency domain is in cycles per m or cycles per FOV, depending on if you are using physical wavelengths or image-based wavelength measures. The "position" in the Fourier domain is unrelated to any displacement in the image domain. Instead, the position in the Fourier domain is the inverse of the wavelength of the Fourier transform of a delta function at that location.

@DaleSpam - This is a very helpful clarification. Your CDF widget is very cool, I'll use it to help explain the frequency filtering aspect of the model to my group once I've got it all working. Thank you!

 

[olivermsun]

@olivermsun - I'm applying a cutoff frequency.

I understand what you are doing now. Cool.

 

[Dale]

@DaleSpam - This is a very helpful clarification. Your CDF widget is very cool, I'll use it to help explain the frequency filtering aspect of the model to my group once I've got it all working. Thank you!

You are welcome, I am glad you like it. If you need modifications, let me know. They are actually relatively easy to make.