Optical Processing

Digital image processing lets us perform operations such as background subtraction, noise reduction, contrast enhancement, filtering, etc. on a digitized image, consisting of a rectangular array of numbers.  One of the operations we can perform is Fourier transform filtering.

A digital image records light intensity as a function of position.  According to Fourier's theorem any reasonably continuous function defined over some distance L can be synthesized by a sum of harmonic (sine and cosine) functions whose wavelengths are integral submultiples of L, (such as L/2, L/3, ...).  Let f(x) be such a periodic function.  Then we may write

,

or f(x) = S-¥+¥Cnexp(iknx)

where the coefficients A0, An, and Bn or Cn can be calculated efficiently using computer algorithms. 
The kn = n2
p/L are equal to 2p times "frequencies" of the component harmonic waves.  We actually do not sum over infinitely many frequencies, but only over frequencies up to the sampling frequency.

Links:

Making Waves  Click on "Run Now" to explore how to synthesize various functions using only sine and cosine waves. The Fourier transform of the amplitude versus frequency plot, which tells us the amplitude (how much) of each sine and cosine function we have to use to  reproduce the chosen function.  This is called a plot of the function in the "frequency domain".

Fourier Series 1  Click the "Applet" button and practice generating various functions using superpositions of sine and cosine waves.

Fourier Series 2  Explore how adding harmonic waves of higher and higher frequency (shorter and shorter wavelength) helps you reproduce sharper features of the function.

1d fast Fourier transform  Choose an actual image, examine the brightness of the pixels along a line through the image (function, upper plot) and the Fourier transform of this function of brightness versus position (function, lower plot).  The lower plot is a plot of the amplitude (how much) versus the frequency of each sine and cosine function we have to use to  reproduce the upper plot.

The Fourier transform of a two-dimensional image is a two-dimensional array of complex numbers giving the coefficients Cn.  It is an important image processing tool.  It represents the image in the frequency plane, i.e. each point represents a particular frequency contained in the image plane.  We can now manipulate the transform in the frequency plane and observe how this affects the synthesized image in the image plane.  Image filtering, smoothing, and edge enhancement are some of the effects we can achieve by manipulation the Fourier transform.

Examine the simple image of a square and its Fourier transform below.

 

Image

Fourier transform

The center of the Fourier transform plot represents the amplitudes of the low frequency sine and cosine waves that make up the image, while the outer regions represent the high frequency waves.

A low-pass filter eliminates all the frequencies above a cutoff frequency.  If we synthesize the image only using only waves with low frequencies (long wavelengths), then sharp features are no longer resolved.   We have smoothed the image.

Applying a low-pass filter to the Fourier transform

The resulting smoothed image
 

A high-pass filter eliminates all the frequencies below a cutoff frequency.  If we synthesize the image only using only waves with high frequencies (short wavelengths), then smooth features of the image are largely eliminated and the sharp edges dominate.   We have edge-enhanced the image.

Applying a high-pass filter to the Fourier transform

The resulting edge-enhanced image
 

This interactive tutorial lets you manipulate the Fourier transform.  You can construct low-pass and high pass filters to smooth the image or produce edge enhancement effects.

An optical processor lets us analyze an image and synthesize the image with various modification by optical means.  A standard optical processor consists of two equal converging lenses with focal lengths f separated be a distance 2f.  The image to be analyzed is placed in the object focal plane of the first lens and illuminated with coherent light incident parallel to the optical axis.  An inverted image with unit magnification appears in the image focal plane of the second lens.  This is the synthesized image and it can be manipulated and modified by placing masks and filters in the image focal plane of the first lens.  This plane is called the frequency plane.  It is also the object focal plane of the second lens.

Why does this work?

Assume the original image in the object focal plane of the first lens is a single slit.  

We have already calculated the Fraunhofer diffraction patter of a single slit.  It is the far-field diffraction pattern of the slit, which can be projected by a converging lens into the focal plane of the lens.  Fresnel's principle guaranties that the relative phases in the focal plane of the lens are equal to the relative phases for a focus at infinity.  We found for the field amplitude a distance x'' from the optical axis

.

Converting to exponential notation we have

.

Let us define kx'' = ksinq @ kq.   We have x'' = f tanq @ fq and therefore x'' µ kx''.  We therefore write

.

Here g(x) = 1, -a/2 < x < a/2, and g(x) = 0 everywhere else. 

The electric field amplitude E(kx'') in the frequency plane (the image focal plane of the first lens) is proportional to the Fourier transform of of the field amplitude g(x) in the object focal plane of the first lens. 

Review of the mathematical properties of the Fourier series and the Fourier transform.

We have found that for the single slit the electric field amplitude E(kx'') in the frequency plane (the image focal plane of the first lens) is proportional to the Fourier transform of of the field amplitude g(x) in the object focal plane of the first lens.  This not only holds for the single slit, but for any (real or complex) amplitude function g(x).  If the field amplitude g(x,y) in the object plane of the first lens depends on two coordinates x and y, then the field amplitude in the frequency plane E(kx'',ky'') is proportional to the the two-dimensional Fourrier transform of g(x,y).

For the single slit we have found the intensity pattern in the frequency plane.

The pattern <I(x'')> has the shape of the square of a sinc function.

For the minima we have sinq = ml/a.  Assume the diameter of the lenses of the optical processor is d.  Because d is finite, the lens can only project M minima into the frequency plane. Light scattered at angles q > qmax misses the lens.  If the numerical aperture of the lens is NA = d/2f = sinqmax, then M = sinqmaxa/l = da/(2fl).

Let us now look at the second lens of the optical processor.  From the field amplitude in the frequency plane it synthesizes an image in its image focal plane.  The amplitude at a point a distance x from the optical axis is found by adding the contributions from all points in the frequency plane.  If the first lens projects the complete diffraction pattern we have

If we let x' = -x then

Therefore E(x') µ g(x),  E(x') is the scaled Fourier transform of E(kx''), and kx'' is the wave number ("spatial frequency") of the component harmonic waves that add up to produce g(x).  The synthesized image is inverted, but otherwise identical to the original image.

E(x'') extends to x'' = ±¥.  But the magnitude of E(x'') decreases rapidly as x'' becomes large.  If the diameter of the lens is large enough so that it intercepts a large number M of minima of the intensity distribution, then the synthesized image will be true to the original.  We need M = (d/2f)(a/l) >> 1 for the synthesized image to be true to the original.

We can modify the synthesized image by modifying the amplitude and phase of E(x'') in the frequency plane.

Assume we put an aperture in the frequency plane with diameter d' such that (d'/2f)(a/l) @ 1.  We now restrict the range of "frequencies" kx'' which can contribute to the synthesized amplitude E(x').  (We only transmit the central maximum of the diffraction pattern.)  We have |kx''| < kd'/2f = kl/a = 2p/a.

Let us choose our units so that a = 1.  The E(kx'')= C sin(kx''/2)/(kx''/2) and

Using Microsoft Excels FFT function we find <I(x')> as shown below.

The image now has smoothed edges and the intensity distribution across the image is no longer uniform.