A Simple Model for Sharpness in Digital Cameras – II

Now that we know from the introductory article that the spatial frequency response of a typical perfect digital camera and lens can be modeled simply as the product of the Modulation Transfer Function of the lens and pixel area, convolved with a Dirac delta grid at cycles-per-pixel spacing

(1)   \begin{equation*} MTF_{Sys2D} = \left|(\widehat{ PSF_{lens} }\cdot \widehat{PIX_{ap} })\ast\ast\: \delta\widehat{\delta_{pitch}}\right|_{pu} \end{equation*}

we can take a closer look at each of those components (pu here indicating normalization).   I used Matlab to generate the examples below but you can easily do the same in a spreadsheet.  Here is the code if you wish to follow along.

Diffraction’s MTF

The lens in our simple example is assumed to be diffraction limited with no aberrations, therefore just producing an Airy pattern on the sensing plane.  At f/16 the pattern would look as shown in Figure 1 of the last post, a tiny white dot at the center of a 1024×1024 linear unit black image.  It can be modeled according to the following continuous equation in the frequency domain:

(2)   \begin{equation*} MTF_{\text diff} = \frac{2}{\pi}[\arccos(s)-s\sqrt{1-s^2}] \end{equation*}

with s the linear spatial frequency f normalized for extinction: s = \frac{f}{\lambda N}.  The MTF is zero when s is greater than 1.  This is what it looks like as a solid and in profile.

Airy MTF 25D
Figure 1.

The graphs show MTF on the Z axis and the spatial resolution in cycles per pixel (c/p) on the  X and Y axes.  As mentioned earlier a pixel here is considered to be a square 5 linear units on the side (more on this later).

If we Fourier Transform the data of the f/16 Airy disc image we get exactly the same result, indicating that the formula is correct and that Matlab’s FFT routine works accurately.  Below is the horizontal central radial slice, parallel to the X-axis of the two dimensional diffraction MTF above, shown as the ‘Modeled’ solid curve; and the same slice off the solid obtained after the Airy image is Discrete Fourier Transformed, shown as ‘Measured’ dots.  It is similar to the graph in Figure 1b but showing positive frequencies only, zoomed into the area of detail that we are normally used to seeing with MTF curves, 0-1 c/p.  It represents the one dimensional MTF in that one direction only:

Airy SFR Model vs Measured
Figure 2.

They are identical.  Because the Airy disc is rotationally symmetric, the one dimensional MTF profile will be the same no matter what radial slice out of the 2D MTF we look at.  This graph shows numerically that the theoretical lens MTF matches that obtained by taking the Discrete Fourier Transform of the image of its Point Spread Function on the sensing plane.

Pixel Aperture MTF

We can perform the same exercise with the Pixel Aperture (aka footprint or area) component.   As discussed earlier a square can be modeled in the frequency domain as

(3)   \begin{equation*} MTF_{PIX} = \left|\frac{sin(\pi f_{x} w)}{\pi f_{x} w}\right|\left|\frac{sin(\pi f_{y} w)}{\pi f_{y} w}\right| \end{equation*}

with f_{x} and f_{y} the horizontal and vertical spatial frequencies better specified later and w the relative width of a pixel.  Equation (3) is shown below as a solid on the left; and on the right its X-Z axis profile, as before.

Square MTF 25D
Figure 3.

A horizontal radial slice from the origin with no vertical frequency component can be modeled in one dimension as follows:

(4)   \begin{equation*} MTF_{PIXh} = |sinc(fw)|= \left|\frac{sin(\pi f w)}{\pi f w}\right| \end{equation*}

with f* the linear spatial frequency in the horizontal direction.  It would of course look the same in the vertical direction, with no horizontal component.

We want to compare this theoretical model to the actual discrete Fourier Transform of the 5×5 linear unit square footprint of the Pixel.  To do so we generate an image of the same size as the Airy image (1024×1024 linear units), fill it with zeros everywhere except in the center where a 5×5 square area of ones is inserted, and take the two dimensional DFT of it.  The plot below shows a horizontal radial slice from the result (‘Measured’) and from the modeled solid(‘Modeled’)

SFR Square and ROI
Figure 4. The impact of the size and shape of the Region of Interest on discrete Fourier transformed image data

Hmm, not perfect.  It looks like the fit is pretty good up to about 1 c/p but after that the model and measured curves part ways.  That’s because, contrary to the Fourier Transform shown as the ‘Modeled’ blue line above, which works on unbounded continuous curves resulting in only a baseband signal, the DFT works on a finite set of discrete input and output values, which produce an infinite set of replicas of the baseband signal in the frequency domain.  If the signal is not band limited to start with, energy from these replicas seeps back into baseband in the form of aliasing.

A rect function is definitely not band limited, so the aliased energy from the replicas adds up quickly.  Its DFT is therefore not the sinc function shown in equation (4) of the continuous model, but an aliased sinc (asinc), periodic sinc, Dirichlet or diric function, which is separable and – when the rect is symmetric – can be written as follows for each axis:

(5)   \begin{equation*} DFT_{rect} = asinc(f w) = \frac{sin( \pi f w)}{sin(\pi f)} \end{equation*}

with f* in cycles per sample and w the width of the square aperture as before.  This can be thought of as the ratio of two sincs, equation (4) divided by sinc(f):

(6)   \begin{equation*} DFT_{rect} = asinc(f w) = \frac{sinc(f w)}{sinc(f)} \end{equation*}

So in this case the output of the DFT differs from the continuous transform by this additional factor in the denominator (its effect shown as the yellow dashed curve above).  To make the output of the DFT match the continuous transform model we therefore need to multiply it by sinc(f), which brings us back to having ‘measured’ results accurately match theory for a square aperture, as shown below just in the horizontal direction:

Square SFR Model vs Measured
Figure 5. Horizontal Radial Slice of 2D MTF of perfect square pixel aperture (also referred to as footprint).

The fact that the DFT’s output was not readily apparent in the earlier Airy plots is because diffraction extinction happens around 0.11cycles/sample (0.55 c/p), much below Nyquist at 0.5 cycles/sample, so that signal is band limited.  That is typically not the case in photographs captured with good technique, where there is effectively one sample per pixel.

Asymmetric Squares

Contrary to the case with diffraction the response of a square is not rotationally symmetric, as is obvious from Figure 3a above.  Its edgy pixels affect image detail differently depending on how the detail is laid out on the sensor.  For instance if we had taken a 45 degree radial slice this would have been the result:

Square SFR Model vs Measured 45 deg
Figure 6.  Diagonal 45 degree radial slice of 2D MTF of perfect square pixel aperture.

The blue line of the model is once again invisible under the measured dots.  This 45 degree diagonal radial slice from the origin is affected equally by both horizontal and vertical spatial frequencies. It can therefore be modeled linearly in one dimension as follows:

(7)   \begin{equation*} MTF_{PIXdiag} =\left|\frac{sin(\pi f w)}{\pi f w}\right|^2 \end{equation*}

This is the MTF of a triangle, because that’s the intensity profile that arises when a square at 45 degrees to an axis is projected onto it.  The linear spatial frequency f* in cycles per picture here is the same as earlier, but care must be taken when plotting these values in c/p that pixels are \sqrt{2} c/p long in this direction.

Intermediate slices will give linear curves somewhere in between these two extremes per equation (3).

Modeling Lens and Pixel Aperture Together

The model predicts the ‘measurements’ obtained from each component in isolation.   We can check how well it works on the components together by applying equation (1).  We can generate a modeled 2D MTF by multiplying together equations (2) and (3) and compare the result to the 2D Discrete Fourier Transform of the image of the Airy disc multiplied point by point by the 2D DFT of the image of the 5×5 linear unit pixel aperture square.

Because diffraction at f/16 is dominant in this example the result appears to be rotationally symmetric, but it is not exactly so as can be gleaned by looking at the minimal differences below extinction (about 0.55 c/p) in Figures 5 and 6.  At larger relative apertures the differences would be more marked.

SFR Product 25D
Figure 7.

We can check theory against practice as before by plotting a radial slice of the Lens x Pixel Aperture solid and its X-Z profile, shown in Figure 7.  Since we now know that the model works well we can also work in one dimension by multiplying diffraction from equation (2) by a single sinc(f*pitch) function for the horizontal slice ‘Model’; and since it’s a square aperture by sinc(f*pitch)^2 for the diagonal slice, as shown below.  The curves are all on top of each other.

Product Slice
Figure 8.

Sampling and System MTF

All right, so now that we have a good model for the MTF of individual components and their interaction all that is left to try is to put it all together and actually ‘sample’ the original Airy image.  The Airy pattern is at the center of a 1024×1024 linear unit sized image –  if it helps pretend that a unit is something like a micrometer but ignore absolute dimensions for now.

Let’s assume for the sake of the exercise that the image of the Airy pattern is continuous, and that the pixelation seen in 3200% zoomed Figure 9a below to the left is simply a conveniently dense number of figures to feed to the DFT algorithm.  The brightness of the image represents the relative number of photons incident on the sensing plane in that small area.  The sensor sits below it and has perfect square pixels with no filters and 100% fill factor.  Each pixel’s area corresponds to 5×5 linear units.  Therefore the resulting ‘sampled’ image is 204×204 pixels.

Figure 9. The left image (9a) represents ‘continuous’ lens PSF intensity corresponding to the number of photons arriving during exposure on the sensing plane. The brightness of the right image (9b) represents the relative number of photons collected by the corresponding pixel on the underlying sensor. Zoomed 3200% and 16000% respectively (the latter 5x as much because one of its pixels corresponds to 5 linear units on the side).

The extent of the 2D MTF of the sampled image will also be one fifth the size, so it shrinks from 2.5 to 0.5 c/p .  Below is the sampled system 2D MTF rendition as a solid to the left; and to the right the customary horizontal slice.  Looks like the model worked perfectly once again, accurately predicting the response of the image as sampled by the square-pixeled sensor.  But wait, what is that little lip at the bottom near 0.5 c/p?

Sampled Diff+Square with Slice
Figure 10.

Oh, right.  Aliasing.  Next.


*Spatial frequency f is explained in this article.