testing fft residual

FFT_img.py

#! /usr/bin/env python

'''

Reads in image data from a FITS file. It transforms the image using FFT.

Plot the original image, the real part and the imaginary part of the transform, and the inverse transform -- real part, which should be equal to the original image, and a imaginary part.

The last panel = residual, and shows the quality of the FFT = original image - real part of the inverse FFT image.

'''


import pyfits
import sys
from numpy.fft import fft2 as fft2d
from numpy.fft import ifft2 as ifft2d

# for future testing fftw performance
# import pyfftw       # pyfftw package
# import fftw3        # pyfftw3 package

import numpy as np
import pylab

file = 'HST_9744_75_ACS_WFC_F555W_drzcroplinearShift.fits'
# Open the FITS file for reading the data
hdr = pyfits.open(file)
print hdr.info()                               
naxis = hdr[0].header['naxis']
print 'Number of axes in this image = %d ' % (naxis)

img = hdr[0].data                               
print "Dimensions of image data = %s.  Type of data is: %s" % (img.shape, type(img))

if (naxis == 2):
    pylab.subplot(231)
    pylab.imshow(np.log(abs(img) + 1))
    (l, m) = img.shape
    mi = img.min()
    ma = img.max()

C = fft2d(img)    
pylab.subplot(232)
pylab.title('FFT and inverse for a FITS image')
pylab.imshow(C.real)
pylab.subplot(233)
pylab.imshow(C.imag)

D = ifft2d(C, C.shape)
pylab.subplot(234)
pylab.imshow(D.imag)

InvFFT = D.real
pylab.subplot(235)
pylab.imshow(np.log(abs(InvFFT) + 1))


diff = InvFFT - np.array(img)
pylab.subplot(236)
pylab.imshow(np.log(abs(diff) + 1), vmin=mi / 10., vmax=ma / 10.)
pylab.xlabel('Residual original image - inverse FFT', fontsize=8)

hdr.close()         
pylab.show()