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 be absorbed by a 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 as it arrives at its virtual back end, the Exit Pupil, we’ll call this blocking/distorting function $P$ . Other than in very simple cases, the Exit Pupil does not necessarily coincide with a specific physical element or Principal surface.^[iv] It is a convenient mathematical construct which condenses all of the light transforming properties of a lens into a single plane.

The complex light field at the Exit Pupil’s two dimensional $uv$ plane is then $U\cdot P$ as shown below (not to scale, the product of the two arrays is element-by-element):

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

Goodman shows that – with a linear, space invariant imaging system, the Fresnel, paraxial and a few other approximations typically valid for well corrected photographic lenses – the intensity of image $I$ projected by the Exit Pupil 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*}$

with

- - $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 distance between the Exit Pupil and Sensing planes, equal to nominal focal length if the lens has unity pupil magnification and is focused at infinity
  - $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 system non-idealities at location (u,v) on the Exit Pupil plane
  - array products are element-by-element
  - $\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 element-wise product of the complex light field from the scene times the generalized Exit Pupil function. In the absence of aberrations the pupil function $P$ 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 two-dimensional phase term will also be present in the pupil function $P$ – representing the net effect of setup and lens imperfections on the light field, as suggested in earlier articles (you can also see an example of defocus and third order spherical aberration at the bottom of this one).

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 uniform 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 uniform field to be normalized to 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. Such an impulse response is referred to as the system’s intensity Point Spread Function $PSF(x,y)$ .

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 take the magnitude of the result squared and, after normalization, obtain an approximation to $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 shown as a solid (left) and as a slice through its 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 circularly 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 its Optical Transfer Function, which just happens to be 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 effective 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 Exit Pupil function $P(u,v)$ , just three lines of code in Octave or Matlab can approximate this apparently complicated feat.

$Numerical Fourier Optics: Diffracted Amplitude Point Spread Function; intensity PSF; MTF$

The substance of the first line is represented by the two dimensional Fast Fourier Transform (fft2) of the sampled Exit Pupil function $P$ , resulting in an estimate of its diffracted pattern on the image plane otherwise known as the Amplitude PSF (the shift is just there so that the PSF will be centered). The second line is a better and faster way to compute the magnitude squared of the amplitude, which is an approximation of the system Intensity PSF. The third line produces the magnitude of the second Discrete Fourier Transform, which is referred to as the system MTF. The values that result above are unnormalized and can then be normalized as desired.

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 Aberration only. Using the rotationally symmetric generalized pupil function described in Equation 2 of the earlier article on defocus and spherical aberrations

$P(\rho) = e^{i\frac{2\pi}{\lambda}(W_{040}\rho^4+W_{020}\rho^2)}$

with 2.5mm pupil radius $\rho$ rotated through 180 degrees to populate the $P(u,v)$ plane for the occasion, 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 the earlier article when analytical formulas were applied instead. They also match similar curves in Wyant^[4]. Here are slices of the relative MTFs then resulting from the third line of code:

Looking good, as expected.

This approach is quite useful because in optics departures from lens perfection are often specified in terms of polynomial fits to 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, the harder part: 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. Basic Wavefront Aberration Theory for Optical Metrology, James C. Wyant and Katherine Creath.}_{5. Thanks to Alan Robinson for his mentoring on this subject.

6. If you would like the Matlab routine used to generate these plots let me know via the form available in the ‘About’ page top right. [Edit: I got enough requests for the code that I decided to make it available despite the fact that it is not designed for wide distribution, see note 5 and the relative disclaimer in the following article.]}

Strolls with my Dog