Fourier transform of an image

Question

I am doing a project where I need to take the FFT and IFFT of a photo in MATLAB. There is a principle in optics where the intensity distribution from a lens is equal to the Fourier transform of the aperture that the light initially goes through. I'm having issues finding the FFT and IFFT of an ideal source.

A circular aperture should have an intensity pattern called Airy disk, which should result from the Fourier Transform of a circle. Here is a photo of the Airy disk that I'm using in my code:

Taking the inverse Fourier transform of the Airy disk should result in an image of a circular aperture, but all I'm seeing is black when I convert to uint8. If I leave it as type double then I see things that also look like Airy disks in the corners. Here is the code I'm using:

pic=imread('airy.jpg');  
pic=rgb2gray(pic); %Convert to grayscale  
imshow(pic)  
title('Initial');  

pic=double(pic); %Convert to double before transforms  
pic=ifft2(pic); %Take IFFT of Airy disk  
pic=real(pic); %Take real part  

figure()  
pic=uint8(pic); %If commented out, black screen has Airy disks in corners  
imshow(pic);  

Does anyone have any ideas on why the IFFT of the Airy disk doesn't return a circle? Thank you

Answer 1

You have the gist of it, but there are a couple of problems. I'm sorry, I don't have MATLAB, but I'll do my best to give you MATLAB-like pseudo-code to work with.

The main problem is that your Airy disk is an image, i.e. amplitude only (the result of taking abs). To reconstruct the circular aperture, you need phase as well. Here's code to first generate your complex Airy disk (PSF), from which you can then reconstruct the circular aperture (pupil function). Starting with this file, circ.jpg:

pupil = imread('circ.jpg');
pupil = rgb2gray(pupil);
% use fftshift to block-swap the fft output (move DC to the center)
complex_airy = fftshift(fft2(pupil));

% now, starting with this complex airy pattern, you were close:
pupil_reconstructed = ifft2(complex_airy);
pupil_reconstructed = abs(pupil_reconstructed);

Below are the amplitude and phase of the complex_airy pattern, obtained using abs(complex_airy) and angle(complex_airy), respectively. (I've included only the centers of the images, to show detail.) The problem with your initial approach is that it's an ill-formed inverse problem: a single Airy amplitude pattern can correspond to an infinite number of pupil functions. When you inverse transformed the amplitude only, you assume the phase is everywhere zero, which--as is clear from the phase image below--it's not. There are statistical approaches for estimating the pupil function from the Airy amplitude alone, using lots of other constraints gleaned from the light source properties and optical system, but I don't think that's what you were after.