Wavefront to PSF to MTF: Physical Units

In the last article we saw that the Point Spread Function and the Modulation Transfer Function of a lens could be easily obtained numerically by applying Discrete Fourier Transforms to its generalized exit pupil function P twice in sequence.[1]

Obtaining the 2D DFTs is easy: simply feed MxN numbers representing the two dimensional complex image of the pupil function in its uv space to a fast fourier transform routine and, presto, it produces MxN numbers that represent the amplitude of the PSF on the xy sensing plane.  Figure 1a shows a simple case where pupil function P is a uniform disk representing the circular aperture of a perfect lens with MxN = 1024×1024.  Figure 1b is the resulting PSF.

Figure 1. 1a Left: Array of numbers representing a circular aperture (zeros for black and ones for white).  1b Right: Array of numbers representing the PSF of image 1a (contrast slightly boosted).

Simple and fast.  Wonderful.  Below is a slice through the center, the 513th row, zoomed in.  Hmm….  What are the physical units on the axes of displayed data produced by the DFT?

Figure 2. A slice through the center of the PSFshown in figure 1b.

Less easy – and the subject of this article as seen from a photographic perspective.

A Numerical Perspective

For this discussion we take a numerical view of the simplified diagram representing the exit pupil and sensing plane (Figure 3 below, not drawn to scale).  Thus the exit pupil function P(u,v) refers to an array of MxN numbers indexed by integers u and v which range from 0 to M-1 and 0 to N-1 respectively.  Similarly for the relative Point Spread Function on the sensing plane I(x,y) and Modulation Transfer Function MTF(f_x,f_y).  The only physical parameters on the diagram are shown in red: aperture diameter D, effective focal length f and sampling pitch du and dx (not shown dv and dy which would however be the same  as du and dx  resp. in a square-pixeled sensor).

Figure 3.  Simplified representation of the air filled space between the exit pupil of the lens and the focal plane, not to scale.

The Physical Units of the Exit Pupil Function

To get the axes of the PSF we first have to determine the physical units of the exit pupil function, which fortunately we do have.  Using the same example as in the last article, with  the diameter of the disk D at 5mm, we can easily put axes on Figure 1a above since I counted that exactly 9 diameters fit within a width and a height of the exit pupil image[2].  That means that the uv pupil plane must correspond to a 45mmx45mm space.

If the plane is sampled n times in the horizontal direction, with n = N = 1024 per Figure 1,  that means that linear physical sample spacing du is 45mm/1024 = 0.0439mm on the physical pupil plane.

Denoting the number of times the diameter of the aperture fits horizontally within the pupil plane as Q, the horizontal axis of the pupil function P(u,v) in Figure 1a at index position u would then correspond to physical value

(1)   \begin{equation*} \begin{aligned} u \cdot \frac{QD}{n} \text{ ,  with u = 0, 1, 2, 3 ...n-1} , \\ u \in \mathbb{Z} \end{aligned} \end{equation*}

and the spacing between samples to du = \frac{QD}{n}, all in the same physical units as D since Q and n are dimensionless.  \mathbb{Z} in this case represents integers 0,1,2… to n-1.

That’s with the origin top left, as shown in Figure 1.  The convention is however to put the origin of the pupil function P(0,0) in the center of the image as shown in Figure 3.  If we shift the image halfway down and to the right[1] we need to let the units change accordingly by modifying u and v.  Shifted u then becomes

    \[ u = -n/2, -n/2+1, -n/2+2,  ...n/2-1}. \]

Of course du remains unaffected by shifting.  Assuming square pixels, dv is the same as du and physical values for v can be obtained similarly (in fact the axes are the same in the case of a square pixel and a square image, as in this example).

From Pupil Function to PSF Units

As described in the previous article, to generate the Point Spread Function I(x,y) on the xy sensing plane we take the Fourier Transform of the Pupil Function P(u,v) after performing coordinate substitution

(2)   \begin{equation*} f_u = \frac{x}{\lambda f} \text{  and } f_v = \frac{y}{\lambda f}, \end{equation*}

with \lambda the wavelength of light and f the effective focal length of the lens (the distance between parallel uv and xy planes in Figure 3 above).

Since our data is digitized with samples at regular intervals we can use the Discrete Fourier Transform to perform the conversion numerically.  The DFT produces an output data set of the same size as the input MxN data set and at a physical spacing equal to the inverse of the input sample spacing.  One over a spacing is sometimes referred to as a ‘frequency’. We can think of it like that if it helps us follow the math – though in this case the result of the Fourier Transform leads to a PSF, clearly a quantity in the spatial domain.

The sequence of values output by the DFT corresponds to an integer multiple of this frequency, from 0 to n-1 times, by definition.  Here is the Discrete Fourier Transform expressed with the variables at hand[3]

(3)   \begin{equation*} \begin{aligned} I(f_u,f_v) = \sum\limits_{u=0}^{\text{N}-1} \sum\limits_{v=0}^{\text{M}-1} e^{-2\pi i(f_u u/\text{N}+ f_v v/\text{M})} P(u,v), \\ f_u,f_v \in \mathbb{Z} \end{aligned} \end{equation*}

\mathbb{Z} in this case represents integers 0,1,2… to N-1 and M-1 for f_u and f_v respectively.  So with  physical units of du = \frac{QD}{n} the horizontal axis of the PSF in Figure 1b at integer index position f_u would correspond to physical value

    \[ f_u \cdot \frac{1}{QD} \text{ , with $f_u$ = 0, 1, 2, 3 ... n-1}. \]

To convert to units of x on the sensor we simply multiply f_u by \lambda f according to (2) so that for instance the fourth pixel on the sensing plane now represents 3\frac{\lambda f}{QD} in physical units.

Substituting  effective f-number N for f/D in that fraction we obtain the following physical quantity at integer index position x on the horizontal axis of the xy sensing plane

(4)   \begin{equation*} \begin{aligned} x \cdot \frac{\lambda N}{Q} \text{ , with $x$ = 0, 1, 2, 3 ... n-1}, \\ x \in \mathbb{Z} \end{aligned} \end{equation*}

and dx = \frac{\lambda N}{Q}, all in the same physical units as \lambda because N and Q are dimensionless.  That’s with the origin top left.  If we want the PSF to be centered we shift the origin half an image both ways as before:

    \[ x = -n/2, -n/2+1, -n/2+2,  ... n/2-1. \]

With the origin shifted to the middle of the image, the 513th sample in Figure 2 becomes the origin.  It is then easy to see from its central slice that the Airy pattern hits first zero at x = 11, corresponding to 11/9*0.5*10 = 6.11 um in physical units (from the last article Q = 9, \lambda = 0.5 um and N = 10).  Correct, because theory predicts that the Airy’s first ring radius should be 1.22 \lambda N.

PSF to MTF Units

To get Modulation Transfer Function MTF(f_x,f_y) from Point Spread Function I(x,y) we take a second Discrete Fourier Transform of the data, as shown by the third line of code up top.  This time there is no coordinate substitution so the physical units corresponding to the output of the DFT are the bona fide linear spatial frequencies present in the PSF.  Since dx = \frac{\lambda N}{Q} the physical units corresponding to integer index f_x on the horizontal axis of the relative MTF are

(5)   \begin{equation*} \begin{aligned} f_x \cdot \frac{Q}{\lambda N} \text{ ,  with $f_x$ = 0, 1, 2, 3 ...n-1}, \\ f_x \in \mathbb{Z} \end{aligned} \end{equation*}

and df_x = \frac{Q}{\lambda N}, all in cycles per the same units as \lambda.  Sometimes it is useful to shift the origin of the frequency diagram to the center of the image.  If so proceed as before to also shift f_x.

Bringing Them All Together

In Octave or Matlab, the three lines of code up top are therefore complemented by the three lines of code below to generate axes for the three planes discussed in this post:

Plots of the pupil, image and frequency spectrum can therefore now be shown with axes in physical units:

Figure 4. Figure 1 and 2 with axes representing physical units.

y, dy, f_y and df_y in the vertical xy direction can be obtained using the same procedure starting from dv.

On Q and \lambda N

You may have noticed that the only variables required to determine the physical units of exit pupil, PSF and MTF are Q, \lambda, f and D.   A Q value of about one order of magnitude is mainly necessary in order to obtain relatively unaliased data from the discrete fast fourier transform routine but for continuous photographic thought experiments it can be considered to be equal to one.  In that case  the physical units of PSFs and MTFs are seen to depend only on a single figure, the product \lambda N.  This realization can provide some neat insights into the basic properties of light on the sensing plane after it has gone through a lens.


Notes and References

1. The substance of the first line is represented by the Fast Fourier Transform of the Pupil Function. The shifts are there to ensure that the resulting PSF is centered. Thanks to Alan Robinson for this particular shifting sequence. The second line is a better and faster way to compute the magnitude of the intensity spectrum squared, thus producing the PSF. The third line takes the magnitude of the Discrete Fourier Transform of the PSF, which produces the MTF.
2. Simplified Figure 3 is not drawn to scale.
3. The definition of the Discrete Fourier Transform can be found here.