There are many applications for taking fourier transforms of images (noise filtering, searching for small structures in diffuse galaxies, etc.). I wanted to point out some of the python capabilities that I have found useful in my particular application, which is to calculate the power spectrum of an image (for later separation of the distribution of stars from the PSF and the noise; see Sheehy et al. 2006).
Below I have posted an example snippet of code, which you can also find on the wiki. But first, some figures to show what we are doing — an astronomical image (left), the 2D power spectrum of the image (middle), and the azimuthally averaged 1D power spectrum (right).
The bulk of the heavy lifting can be done using SciPy’s fftpack. I also wrote radialProfile.py and it is very crude at the moment.
from scipy import fftpack
import numpy as np
import pylab as py
image = pyfits.getdata(‘myimage.fits’)
# Take the fourier transform of the image.
F1 = fftpack.fft2(image)
# Now shift the quadrants around so that low spatial frequencies are in
# the center of the 2D fourier transformed image.
F2 = fftpack.fftshift( F1 )
# Calculate a 2D power spectrum
psd2D = np.abs( F2 )**2
# Calculate the azimuthally averaged 1D power spectrum
psd1D = radialProfile.azimuthalAverage(psd2D)
# Now plot up both
py.imshow( np.log10( image ), cmap=py.cm.Greys)
py.imshow( np.log10( psf2D ))
py.semilogy( psf1D )
If you have suggestions for the above examples or other related code snippets or helper functions, let me know or comment below. I have posted this example on the wiki and will keep it updated with your suggestions.