Goodman, in his excellent Introduction to Fourier Optics^{[1]}, describes how an image is formed on a camera sensing plane starting from first principles, that is electromagnetic propagation according to Maxwell’s wave equation. If you want the play by play account I highly recommend his math intensive book. But for the budding photographer it is sufficient to know what happens at the exit pupil of the lens because after that the transformations to Point Spread and Modulation Transfer Functions are straightforward, as we will show in this article.

The following diagram exemplifies the last few millimeters of the journey that light from the scene has to travel in order to smash itself against our camera’s sensing medium. Light from the scene in the form of field arrives at the front of the lens. It goes through the lens being partly blocked and distorted by it (we’ll call this blocking/distorting function ) and finally arrives at its back end, the exit pupil. The complex light field at the exit pupil’s two dimensional plane is now as shown below:

Goodman shows that – with a linear, space invariant imaging system, the paraxial and a few other approximations typically valid for well corrected photographic lenses – the intensity of image projected by exit pupil onto the focal plane is^{[5-16]}

with

- the light intensity of the image on the sensing plane in two dimensions at location (x,y)
- the wavelength of light
- the effective focal length of the lens
- the light field from the scene at the front of the lens in the two dimensional (u,v) coordinates of the exit pupil plane
- the generalized pupil function, condensing the net effects of aperture and lens aberrations at location (u,v) on the exit pupil plane
- indicating time averaging to deal with incoherent light
^{[2]}and weighted averaging to deal with polychromatic light.^{[3]}

#### Simpler Than It Seems

As formidable as it appears this equation is actually not that complicated. The term is just the product of light from the scene times the generalized pupil function. In the absence of aberrations the pupil function of a lens is real and simply the shape and transmittance of the aperture projected onto the exit pupil plane. In photography we typically assume that it is a circular disk opening that lets all light inside it through and stops all light outside of it. On the other hand in the presence of aberrations a phase term will also be present in the pupil function – representing the net effect of lens imperfections on the light field, as suggested in earlier articles.

With a simple change of coordinates

the rest of Equation 1, the integrals and exponential term, can be seen to be the two dimensional Fourier Transform of the above aberrated light field at the exit pupil plane.

#### From Exit Pupil Wavefront to Point Spread Function

An impulse from a single on-axis point-source at infinity will produce an instantaneous coherent field at the front of the lens, which means Equation 1 can be applied as-is for monochromatic light^{[3]}. If we further assume this field to be of unit amplitude, that is , then it is clear from Equation 1 that the intensity on the sensing plane is just the normalized magnitude of the two dimensional Fourier Transform of the exit pupil function squared. The intensity on the sensing plane resulting from the impulse response of a lens is called its intensity Point Spread Function.

Let’s try the theory out. Assume that we have such an impulse of light going through a perfect, unaberrated lens of circular aperture. Since there are no imperfections the exit pupil function will just be a disk, letting through light inside it and blocking light outside of it. Mathematically we can represent this transmittance as ones inside the disk and zeros outside of it, as shown below for an aperture with diameter of 5mm (white representing ones and black zeros).

If we have this image of ones and zeros digitized by, say, a 1024×1024 sample grid we can feed it to a two dimensional Discrete and Fast Fourier Transform – then square, normalize and take the magnitude of the result to obtain , the intensity Point Spread Function that such a perfect lens would project onto the sensing plane.

That’s pretty easy to accomplish in Matlab, Octave or even Excel. Assuming mean light of wavelength = 0.5um, aperture diameter = 5mm and focal length = 50mm, the resulting PSF in 3D (left) and as a slice through the center (right) appear as follows:

Hey, it’s our old friend the Airy Disk, the diffraction pattern expected from a lens of circular aperture without aberrations. Effective f-number is about 10 () so first zero hits at per theory, so far so good. What if the aperture is not circular but, as often happens with photographic cameras, it is a polygon – say an octagon? No problem, let’s draw up an octagonal exit pupil and fill it with ones

then take the normalized magnitude of its 2D DFT squared to obtain the relative intensity PSF on the sensing plane. Note that the PSF from the octagonal pupil function is no longer rotationally symmetric (below left), the octagonal pattern is visible in the ripples and there are some differences in the intensity distribution seen in slices through the center of the PSF on a logarithmic scale (below right)

Good enough for most purposes at f/10, as we all know from happily using cameras with octagonal stops.

By the way, the aperture’s transmission does not necessarily have to be all or nothing. Some lens designers attempt to vary transmission attenuation near the aperture’s edges smoothly by using an apodization element in order to induce desirable qualities in their product, for example softer bokeh. The Minolta/Sony 135mm STF, FE 100mm F2.8 STF, upcoming Canon EF 135mm f/2 and others come to mind. To simulate those effects all we would need to do is draw a grayscale exit pupil of the chosen aperture shape with transitions varying smoothly between zero and one as desired.

#### From PSF to Modulation Transfer Function

The Modulation Transfer Function of a linear, space invariant imaging system is just the normalized magnitude of the Fourier Transform of the relative intensity PSF. Performing that operation on the PSFs calculated above we obtain the following two dimensional spatial frequency response for the lenses, of which the central slice in the vertical direction is shown below:

As expected extinction occurs at 1/ or 200 lp/mm with f-number equal to 10 and mean light wavelength 0.5um. A perfect lens with a square aperture would produce a straight line MTF from 1 to the extinction spatial frequency. Intuitively the performance of the octagonal aperture should fall in between circular and square, as indeed seen above.

#### Easy Numerical Application

So in the end apparently complex Equation 1 turns out to have quick and easy numerical application. To go from the two dimensional aberrated wavefront at the exit pupil to intensity PSF on the sensing plane to linear spatial resolution MTF in the frequency domain all we need to do is take the Fourier Transform of the complex pupil function twice, plus a bit of massaging in between. Given a complex pupil function , just three lines of code in Octave or Matlab accomplish this apparently complicated feat.

The substance of the first line is represented by the two dimensional fast fourier transform (fft2) of the pupil function, the shifts are just there so that the PSF will be centered.^{[4]} The second line is a better and faster way to compute the magnitude squared of the amplitude spectrum, which is the system PSF. The third line produces the magnitude of the second Discrete Fourier Transform, which is the system MTF.

#### Checking Numerical Results Against Theory

Here are a few further checks with circular aperture pupil functions that include varying degrees of Seidel defocus and third order Spherical Aberrations only. The first two lines of code result in the following central slices of the relative PSFs:

If these graphs look familiar it’s because they reflect the same parameters and format used in a previous article when analytical formulas were applied instead. They also match similar curves in Wyant^{[5]}. Here are slices of the relative MTFs resulting from the third line of code:

They seem to work fine to me but if you spot any inconsistencies please let me know.

This approach is quite useful because in optics departures from lens perfection are often specified in terms of the parameters of the generalized exit pupil function – in Seidel, Zernike or other notation. And we have seen in this article that if we know we can easily estimate numerically the performance of the relative lens on the sensing plane in the spatial and frequency domains.

Next: how the units work.

#### Notes and References

_{1. Introduction to Fourier Optics, 3rd Edition. Joseph W Goodman.
2. Equation 1 can be used without time averaging with coherent light or with incoherent light when the coherence area on the object is small compared with a resolution cell size in object space (Introduction to Fourier Optics section 6.1.3).
3. For polychromatic light Equation 1 needs to be evaluated at the various wavelengths present and the results combined in a weighted average according to each wavelength’s relative energy. We did something similar here.
4. Thanks to Alan Robinson for the precise shifting sequence in the first line of code.
5. Basic Wavefront Aberration Theory for Optical Metrology, James C. Wyant and Katherine Creath.
6. If you would like the Matlab scripts used to generate these plots let me know via the form available in the ‘About’ page top right.
}