A Simple Model for Sharpness in Digital Cameras – Spherical Aberrations

Spherical Aberration (SA) is one key component missing from our MTF toolkit for modeling an ideal imaging system’s ‘sharpness’ in the center of the field of view in the frequency domain.  In this article formulas will be presented to compute the two dimensional Point Spread and Modulation Transfer Functions of the combination of diffraction, defocus and third order Spherical Aberration for an otherwise perfect lens with a circular aperture.

Spherical Aberrations result because most photographic lenses are designed with quasi spherical surfaces that do not necessarily behave ideally in all situations.  For instance, they may focus light on slightly different planes depending on whether the respective ray goes through the exit pupil closer or farther from the optical axis, as shown below:

371px-spherical_aberration_2
Figure 1. Top: an ideal spherical lens focuses all rays on the same focal point. Bottom: a practical lens with Spherical Aberration focuses rays that go through the exit pupil based on their radial distance from the optical axis. Image courtesy Andrei Stroe.

If one imagines the sensing plane being inserted parallel to the lens at various distances from it, it can be seen that in the presence of SA  light from a distant star projected by the lens onto the sensing plane would not converge to a single point but to a disk – and that the size and position along the optical axis of the smallest achievable disk (the circle of least confusion) would change as the lens is stopped down, resulting in focus shift.  This effect interacts with axial chromatic aberrations  (spherochromatism) and tends to be the key determinant of lens ‘sharpness’ loss in the center of most in-spec, well corrected lenses at wide apertures[1].

SA and Defocus Point Spread Function

As figure 1 shows, with Spherical Aberration there is no longer a single plane of focus so SA and defocus interact.  Wyant and Creath[2] suggest that the irradiance on the sensing plane (PSF) of the diffraction pattern of a circular aperture with a rotationally symmetric pupil function illuminated by a uniform light source can be written as

(1)   \begin{equation*} PSF_{lens}(r) = \frac{P\pi}{(\lambda N)^2}\left|\int_{0}^{1}\gamma(\rho)\cdot J_0 \left[\frac{\pi r}{\lambda N}\rho\right]\cdot \rho \cdot d\rho\right| ^2 \end{equation*}

with \gamma(\rho) representing the following function if defocus and third order SA are present

(2)   \begin{equation*} \gamma(\rho) = e^{i\frac{2\pi}{\lambda}(W_{040}\rho^4+W_{020}\rho^2)} \end{equation*}

where:

  • r, the radial distance from the center of the circularly symmetric PSF on the sensing plane
  • P, the power of light contained within the pupil
  • \lambda, the wavelength of light
  • N, the effective f-number
  • \rho, the radial distance on the exit pupil plane from its center, normalized to one
  • W_{040} and W_{020}, the peak Optical Path Difference from the principal plane of third order SA and Defocus respectively, in wavelengths.

A slice through the center of the PSF looks like this when varying defocus coefficient W_{020} only (no SA)

Figure 2. PSF of ideal lens with circular aperture and varying amounts of defocus.
Figure 2. PSF of ideal lens with circular aperture and varying amounts of defocus only. The units of the x-axis are incorrectly labeled, they should read ‘lambda.N’.

When the image is in focus the PSF has the shape of an Airy disk, as expected.  In that case a slice of the PSF would look as follows when Spherical Aberration coefficient W_{040} is varied:

Figure 3. Ideal Lens with circular aperture and varying amounts of third order SA aberration.
Figure 3. Ideal Lens with circular aperture, no defocus and varying amounts of third order Spherical Aberration only.  The units of the x-axis are incorrectly labeled, they should read ‘lambda.N’.

These graphs thankfully look like the ones in Wyant so apparently my Octave/Matlab implementation of the formula works[3].  They represent a slice through the center of the 2D PSFs, so to obtain the full 2D PSF all one needs to do is rotate the relative radial slice around zero through 180 degrees.

SA and Defocus Modulation Transfer Functions

We are after a radial slice through the center of the 2D MTF, which yields the Transfer Function of the system in the direction of the slice.  Taking the discrete Fourier Transform of a one dimensional slice as depicted above will not produce the desired result, even though the PSF is circularly symmetric.    We have two options: either expand the radial PSF above to two dimensions by rotating it and taking the normalized modulus of the 2D DFT of that; or obtain the MTF directly as the normalized modulus of the Fourier-Bessel Transform of equation (1) as is.  I opted for the latter in order to stay in FT, instead of DFT, math.

The Fourier Transform of a circularly symmetric function is itself circularly symmetric and can be obtained by applying the Fourier-Bessel Transformation (also known as the Hankel transform of zero order) to equation (1)[4]. The result is a radial slice through the 2D MTF of the 2D PSF of the combination of diffraction, third order Spherical Aberrations and defocus:

(3)   \begin{equation*} MTF_{lens} = \left|2\pi \int_{0}^{\infty} PSF_{lens}(r)\cdot J_0(2\pi rf)\cdot r\cdot dr\right|_{pu} \end{equation*}

with f[5] the spatial frequencies of interest, pu indicating normalization to one at the origin, PSF_{lens}(r) and r as defined earlier.   I don’t know how to solve equation (3) analytically but it’s easy enough to do it numerically in Octave/Matlab[3].  Doing so produces the following graph when varying SA with zero defocus:

Figure 4. Radial slice of 2D MTF of third order Spherical Aberration, zero defocus, N = 5.6, 0.535um light and pitch = 5.
Figure 4. Radial slice of 2D MTF of third order Spherical Aberration, zero defocus, N = 5.6, 0.535um light and pitch = 5um.

The first curve with zero SA and defocus is just the MTF of the diffraction pattern caused by a circular aperture.  When varying defocus without SA these are the resulting curves instead:

mtf-of-defocus-vs-w020
Figure 5. Radial slice of 2D MTF of Defocus, no Aberrations, N = 5.6, 0.535um light and pitch = 5um.

What do you know, this last graph matches MTFs generated by Hopkins’ dedicated defocus equation explored in the last post:

hopkins-vs-wyant
Figure 6. Comparison of Defocus MTF produced by Hopkins’ formula from the last post and from Equation (3) in this one. F/5.6, 0.535um light,

SA and Defocus Interact to Produce the CLC

The signs of the W_{040} and W_{020} coefficients refer to whether the wavefront is ahead or behind the reference sphere.  They interact, so when SA is present the plane of best focus is not necessarily with W_{020} at zero, as can be gleaned from figure 1.  The plane of best focus is the one where the diameter of the PSF on the sensing plane is smallest, the resulting ‘disk’ denoted the circle of least confusion (CLC).  Wyant suggests that, geometrically, the CLC is smallest when W_{020} is equal to -1.5 * W_{040}. Indeed a curve with the coefficients in those proportions is the best one below:

mtf-of-sa-and-defocus
Figure 6. Radial slice of 2D MTF of the interaction of Defocus and third order Spherical Aberrations, N = 5.6, 0.535um light and pitch = 5um.

Taking diffraction into account instead the smallest CLC should be produced when the combined optical path differences sum to zero.

So equation (3) allows us to calculate the combined two dimensional Modulation Transfer Function of diffraction, defocus and third order Spherical Aberration for an otherwise perfect lens with a circular aperture.

Bringing It All Together

Together with the Transfer Functions of pixel aperture, sampling pitch and AA discussed in previous articles we should now have all the principal ingredients necessary to roughly model the linear spatial resolution (‘sharpness’) of a photographic digital camera matched to an excellent lens in the center of the field of view:

(4)   \begin{equation*} MTF_{sys} = MTF_{lens} (\cdot MTF_{AA}) \cdot MTF_{pixel} **\delta_{pitch} \end{equation*}

with \cdot indicating element-wise multiplication and ** convolution.  We could for instance simulate the effect of varying f-number, pixel size and shape while taking into account these basic lens imperfections.  And by analyzing the raw data of an image captured with good technique we should be able to learn a fair amount about the imaging system that took it.

PS.  This article allows the simulation of complex lens imperfections and aberrations from a numerical perspective.

Notes and References

1. See this excellent paper for for a more in-depth treatment of SA.
2. Basic Wavefront Aberration Theory for Optical Metrology, James C. Wyant and Katherine Creath.
3. The Matlab/Octave scripts used to generate these plots can be downloaded from here.
4. Introduction to Fourier Optics, 3rd Edition. Joseph W Goodman, equation 2-31.
5. See here for a description of linear spatial frequency f.