Aberrated Wave to Image Intensity to MTF

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  U 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 P) and finally arrives at its back end, the exit pupil.   The complex light field at the exit pupil’s two dimensional uv plane is now  U\cdot P as shown below:

Figure 1. Simplified schematic diagram of the space between the exit of a camera lens and its sensing plane. The space is filled with air.

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 I projected by exit pupil P onto the xy focal plane is[5-16]

    \begin{equation*} \begin{split} I(x,y) = \langle\frac{1}{\lambda^2 f^2}\left|\iint\displaylimits_{-\infty}^{\infty}U(u,v)P(u,v)e^{-i2\pi\frac{xu+yv}{\lambda f}}dudv\right|^2\rangle \\ \text{Equation (1)} \end{split} \end{equation*}


  • I(x,y) the light intensity of the image on the sensing plane in two dimensions at location (x,y)
  • \lambda the wavelength of light
  • f the effective focal length of the lens
  • U(u,v) the light field from the scene at the front of the lens in the two dimensional (u,v) coordinates of the exit pupil plane
  • P(u,v) the generalized pupil function, condensing the net effects of aperture and lens aberrations at location (u,v) on the exit pupil plane
  • \langle\rangle  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  U\cdot P 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 uv 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

    \[ f_u = \frac{x}{\lambda f} \text{   , and   } f_v = \frac{y}{\lambda f}. \]

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 uv 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  U=1, then it is clear from Equation 1 that the intensity I(x,y) on the sensing plane is just the normalized magnitude of the two dimensional Fourier Transform of the exit pupil function P(u,v) 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 P(u,v) 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 D of 5mm (white representing ones and black zeros).

Figure 2.  Pupil function of a lens without aberrations and effective circular aperture with diameter 5mm.

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 I(x,y), 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  \lambda = 0.5um, aperture diameter D = 5mm and focal length f = 50mm, the resulting PSF in 3D (left) and as a slice through the center (right) appear as follows:

Figure 3.  3D  and central slice of PSF produced by applying Equation (1) to the pupil function of Figure 2.

Hey, it’s our old friend the Airy Disk, the diffraction pattern expected from a lens of circular aperture without aberrations.  Effective f-number N is about 10 (=f/D) so first zero hits at 1.22 \lambda N = 6.1 um 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

Figure 4.  Pupil Function of a perfect lens with octagonal aperture of 5mm diameter.

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)

Figure 5. 3D rendition and central slice of PSF obtained by applying Equation (1) to the pupil function shown in Figure 4.

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:

Figure 6.  MTF of the PSFs in Figures 3 and 5, representing the linear spatial resolution performance of perfect lenses with 5mm circular and octagonal apertures respectively.

As expected extinction occurs at 1/\lambda N or 200 lp/mm  with f-number N equal to 10 and mean light wavelength \lambda 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 P(u,v), 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:

Figure 7.  Seidel defocus and 3rd order Spherical Aberration with W020 and W040 of 0, 0.25, 0.5, .075, 1 wavelength optical path difference with otherwise perfect circular aperture lenses.

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:

Figure 8. MTFs of the PSFs in Figure 7.

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 P(u,v) – in Seidel, Zernike or other notation.  And we have seen in this article that if we know P(u,v) 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.