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 systematically 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 Aberrations there is no longer a single plane of focus so SA and defocus interact to produce a plane of minimum blur (presumably one that results in the Circle of Least Confusion). 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 generalized pupil 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 Differences from the reference sphere at the Exit Pupil – of third order SA and Defocus respectively, in the same units as light wavelength
$J_0$ , the zeroth order Bessel function of the first kind.

Mathematically inclined readers will recognize the heart of Equation (1) as a Fourier Transform of complex pupil function $\gamma$ , as better described in the next section. A slice through the center of the PSF looks as follows when varying defocus coefficient $W_{020}$ only (no SA)

Figure 2. Center slice of PSF of ideal unaberrated lens with circular aperture and varying amounts of defocus.

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 only third order Spherical Aberration coefficient $W_{040}$ is varied:

Figure 3. Center slice of PSF of Ideal circular aperture lens, with varying amounts of third order Spherical Aberrations only.

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 the origin.

SA and Defocus Modulation Transfer Functions

We are after a radial slice through the center of the 2D MTF, which yields the spectral performance 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, as opposed to 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(f) =\left|\int_{0}^{\infty} PSF_{lens}(r)\cdot J_0(2\pi f \cdot r)\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 third order SA with zero defocus:

Figure 4. Radial slice of 2D MTF of ideal circular aperture lens with third order Spherical Aberrations only – no defocus,.

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:

Figure 5. Radial slice of 2D MTF of ideal unaberrated circular aperture lens with Defocus only.

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

Figure 6. Comparison of Defocus MTF produced by Hopkins’ Equation 1 in the last post and from Equation (3) in this one.

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 sharpest ‘focus’ as intended by photographers is the one where the extent of the PSF on the sensing plane is smallest.

Geometrically, the resulting disk is denoted the circle of least confusion (CLC) and Wyant suggests that its diameter 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 the PSF has infinite support per Equation (1) so its energy is no longer neatly encircled in a disk. Definitions of ‘extent’ can therefore vary significantly depending on the application, for instance percentage of encircled energy (EE) or full width at half maximum (FWHM). One of the better ones for our purposes is minimum RMS wavefront error, which results in the ‘smallest’ value when $W_{020}$ = – $W_{040}$ :

MTF of W040 vs W020 — Figure 7. Modulation Transfer Functions with SA3 W040 held constant at half a wavelength of p-v Optical Path Difference and defocus W020 is varied at shown. The MTF curve closest to pure diffraction occurs when W020 = -W040.

Many other criteria are possible, for instance in the frequency domain where one could envision using the area under the MTF curve, or perceptually weighted SQF and CPIQ Acutance.

In terms of the MTF50 metric used extensively in these pages, performing an empirical simulation using the code described in the post scriptum article yields a $\frac{W_{020}}{W_{040}}$ factor of about -0.9/-0.95 (‘a’ in the Figure below) . That is peak MTF50 is achieved when $W_{020}$ is about equal to -0.95* $W_{040}$ , close to minimum RMS.

Figure 8. MTF50 when SA3 W040 is held constant at the shown p-v Optical Path Difference and defocus W020 is allowed to vary proportionately (their ratio indicated by factor a). The maximum occurs approximately when W020 = -W040.

Bringing It All Together

So Equation (3) allows us to calculate the combined two dimensional Modulation Transfer Function of diffraction, defocus and third order Spherical Aberrations for an otherwise perfect lens with a circular aperture ( $MTF_{lens}$ ). Together with the Transfer Functions of AA, pixel aperture and sampling pitch 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 $**$ 2-dimensional 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, as explored in a couple of following posts.

PS. This article allows the simulation of more complex lens imperfections and aberrations by taking a wholly numerical approach.

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 .}

6 thoughts on “A Simple Model for Sharpness in Digital Cameras – Spherical Aberrations”

Sheng says:

August 20, 2023 at 8:26 am

Thanks a lot for the amazing series of posts. I learnt a lot.
I am struggling with using the MTF or the OTF (seems you are mixing them in the figure, but it is OK) formula (including the Matlab code generated numerical values) to perform image convolution:
1. Assume I don’t use the PSF, using PSF should be OK since FFT(PSF) will have both real and imaginary parts.
2. Using OTF formula (your Matlab code), spin it (assume rotational symmetric) to get 2D OTF. But I only have the real part, no phase info. Would this create a problem if I directly use it to * fft(image)? If there is issue, anyway to solve it?
Thank you.

1. Jack says:
  
  August 20, 2023 at 1:49 pm
  
  Hi Sheng, the code calculates the radial MTF, which is just a normalized version of the radial |OTF|. So to get the radial OTF simply take away the abs() function in the relative mtf line, and normalize appropriately.
  
  Or you could take a numerical approach:
  https://www.strollswithmydog.com/wavefront-to-psf-to-mtf-physical-units/
  
  Jack
  
Sheng says:

August 23, 2023 at 6:26 am

Thanks Jack. Probably I didn’t ask the correct question. I understand the differences between MTF and OTF. My question is more about if I only know OTF’s real part, how can I get PSF or convolve with an image (by directly multiply the fft(image), then ifft). The reason I am asking is because:
1. The 2nd eq. in your previous chapter (https://www.strollswithmydog.com/a-simple-model-for-sharpness-in-digital-cameras-defocus/), only has absolute value.
2. I calculated the value of 3rd eq. in this chapter above, without using abs in front, the resulting value (should be OTF) is still a vector of real values, there are no imaginary part.

I thought the OTF should have both real and imaginary parts as if I directly perform fft(PSF). Thanks.

1. Jack says:
  
  August 23, 2023 at 8:39 am
  
  Hello Sheng,
  
  That’s because the Fourier Transform of a fully real, rotationally symmetric, 2D function (like those PSFs) is also only real – so I think you are good to go. More detail on this in the section on the OTF in the following article
  https://www.strollswithmydog.com/fourier-optics-complex-pupil-function/#OTF
  
  Jack
  
  1. Sheng LIu says:
    
    August 24, 2023 at 5:21 am
    
    Thanks, Jack! I see above in eq(3) you are using the Bessel function for calculating MTF/OTF, instead of using fft(PSF), is this because the PSF you got in eq.(1) is 1D PSF instead of 2D? And the fft(1D PSF) ≠ cross-section(fft(2D PSF))? Moreover, is the fft(2D PSF) equivalent of performing 2D conversion of the 1D OTF get from eq.(3) (by rotating the 1D curve) ?
    
    Thank you.
    
    1. Jack says:
      
      August 25, 2023 at 3:07 pm
      
      Hi Sheng,
      
      Equation (1) produces a radial slice of the relative 2D rotationally symmetric PSF. Equation (3) is the FT of such a PSF in just one direction, that of the slice.
      
      A slice is different than a projection (such as a Line Spread Function, LSF) – the latter is made up of components of every pixel in the 2D PSF, the former is not.
      
      You can obtain the 2D MTF by taking the 2D FFT of the PSF slice from Equation (1) rotated around 360 degrees. Careful though with slices of the 2D MTF. The Fourier-Slice Theorem says that a slice of the 2D MTF corresponds to a 1D projection of the 2D PSF, in the direction perpendicular to that of the slice.
      
      Jack

Strolls with my Dog