Point Spread Function

The point spread function (PSF) is the image that an imaging system forms of a single point of light. In a camera it is the combined blur of the lens and the sensor: diffraction at the aperture, defocus, aberrations, any optical low-pass filter and the light-sensitive area of each pixel. When one PSF describes the whole image, the recorded image is the scene convolved with it, so the PSF characterizes the system’s sharpness completely. Its Fourier transform, the optical transfer function (OTF), says what happens to each spatial frequency, and the magnitude of the OTF, the modulation transfer function (MTF), is the curve lens and camera tests report. The site writes the term both as “point spread function” and as “point-spread function”.

Image formation introduces the PSF and the MTF in a paragraph each, as one stage of the imaging chain. This article treats them in full: the components and how they combine, the diffraction-limited OTF of a circular pupil and its cutoff frequency, the pixel aperture’s transfer function, the sign changes that the MTF hides, and the slanted-edge method by which an MTF is measured from an image.

Definition

The PSF $h(x, y)$ of a linear imaging system is the irradiance it produces in the image plane from a point source of unit power. For mutually incoherent light, such as ordinary scenes, irradiances add, so the image of an extended scene is the sum of shifted and scaled copies of the PSF, one per scene point (Bankhead, 2022–2025, “Blur & the PSF”; Optical Microscopy gives the same statement for microscopes). If $h$ does not depend on where the point lies, the system is shift-invariant and the image is a convolution:

$$g(x, y) = (h * f)(x, y) = \iint h(x - x', y - y')\, f(x', y')\, dx'\, dy'.$$

Bankhead describes the PSF as a linear filter applied to the sharp data we would like to acquire but never can. Two derived functions are used for measurement. The line spread function (LSF) is the image of an infinitely thin line, the PSF integrated along the line’s direction; the edge spread function (ESF) is the image of a straight step edge, the integral of the LSF across the edge. Both follow from the definition by linearity (our statement).

Intuition

Each point of the scene becomes a small blob in the image, and blobs from neighbouring points overlap and add. Bankhead illustrates this with an image seen as a set of points, each replaced by a more diffuse blob; the larger the blob, the blurrier the image. Fine detail, such as two close points or a fine grating, survives only if the blobs are small compared with its spacing.

The frequency view says the same thing in another way. Any image is a sum of sinusoidal gratings at different frequencies and orientations, and a shift-invariant blur changes each grating only in its contrast and position, never its frequency. Levoy’s slides (Sampling and Pixels, slide 17) describe the MTF as “contrast at each frequency relative to original signal”, lost to misfocus, aberrations, diffraction and manufacturing defects among other causes. A wide PSF lowers the contrast of fine gratings strongly; a narrow one keeps it. Levoy (slide 18) contrasts two MTF curves: in one, contrast stays high but drops off at a relatively low resolution; in the other, finer detail is preserved but at low contrast throughout. In our reading, a single number such as a resolution limit therefore does not describe a lens.

Formal Definition

Optical and modulation transfer functions

The OTF is the two-dimensional Fourier transform of the PSF, normalized to 1 at zero frequency:

$$H(\nu_x, \nu_y) = \frac{\iint h(x, y)\, e^{-2\pi i (\nu_x x + \nu_y y)}\, dx\, dy}{\iint h(x, y)\, dx\, dy}, \qquad \mathrm{MTF}(\nu_x, \nu_y) = |H(\nu_x, \nu_y)|.$$

The normalization makes the OTF independent of the system’s overall transmission; in discrete code it means the sampled PSF sums to 1, or its transform is divided by its value at zero frequency. Frequencies are in cycles per unit length on the sensor: cycles/mm (equal to line pairs per mm, one cycle being a dark and a light line) or cycles per pixel. The Nyquist frequency of a sensor with pitch $p$ is $1/(2p)$, 0.5 cycles/pixel; sampling and aliasing explains why.

By the convolution theorem (Fourier transform), the image spectrum is the scene spectrum times $H$. A sinusoidal grating $1 + m \cos(2\pi \nu x)$ with modulation $m$, the Michelson contrast $(I_{\max} - I_{\min})/(I_{\max} + I_{\min})$ that Image Formation’s example measures, comes out as $1 + m\,|H(\nu)| \cos(2\pi \nu x + \phi(\nu))$ for a real PSF, with $\phi$ the phase of $H$ (our derivation). The MTF is therefore the ratio of output to input modulation of a sinusoid, and the phase is a shift of the grating. A bar target is not a sinusoid: it contains harmonics, so its contrast at a given bar frequency is not the MTF at that frequency.

For a PSF that is symmetric about its centre, $H$ is real, and its phase is 0 or $\pi$: a negative OTF means that a grating at that frequency comes out with its bright and dark bars exchanged, a contrast reversal. The MTF, being a magnitude, does not show it.

The 1-D Fourier transform of the LSF along a direction equals the OTF along the line through the origin in that direction (the projection-slice property of the Fourier transform; our statement). This is why an edge or a line, which gives the LSF, measures one slice of the OTF at a time.

Components, which convolve

A camera’s PSF is the convolution of the PSFs of its stages, so its OTF is the product of their OTFs (convolution theorem). Image Formation lists the stages that make up a pixel’s PSF, the lens blur convolved with the light-sensitive area of the photosite, and Image Deblurring gives Szeliski’s list of lens blur, anti-aliasing filter and pixel shape. Three components have closed forms.

Diffraction. For an aberration-free lens with a circular pupil in incoherent light, ISETCam’s dlCore.m computes the OTF as

$$H_{\mathrm{diff}}(\nu) = \frac{2}{\pi}\left( \arccos s - s\sqrt{1 - s^2} \right), \quad s = \frac{\nu}{\nu_c} < 1, \qquad H_{\mathrm{diff}} = 0 \text{ for } s \ge 1,$$

with the incoherent cutoff frequency $\nu_c = A/(\lambda\, d_i)$, where $A$ is the aperture diameter and $d_i$ the distance from aperture to detector; dlMTF.m notes that for a distant object $d_i$ is the focal length, so $\nu_c = 1/(\lambda N)$ with $N$ the f-number. (The comment in dlCore.m writes the square root as $\sqrt{(1 - s)^2}$; the code computes $\sqrt{1 - s^2}$, which is the form above.) The comments cite a 1990 Applied Optics paper by Subbarao for the formula, which was not read here. $H_{\mathrm{diff}}$ is real and non-negative and falls to zero at $\nu_c$: no lens of f-number $N$ passes any contrast above $1/(\lambda N)$, whatever its quality. Optical Microscopy describes this cutoff as one way to define the diffraction limit, beyond which detail is absent rather than blurred. For green light, $\lambda = 0.55$ µm, $\nu_c$ is 649 cycles/mm at f/2.8, 227 at f/8 and 114 at f/16 (our arithmetic); a sensor with a 3 µm pitch has its Nyquist frequency at 167 cycles/mm, so at f/16 the optics cut off below it.

The corresponding PSF is the Airy pattern, $h(r) \propto \left( 2 J_1(v)/v \right)^2$ with $v = \pi r/(\lambda N)$ on the sensor, which is Optical Microscopy’s form $v = 2\pi\,\mathrm{NA}\,r/\lambda$ with $\mathrm{NA} \approx 1/(2N)$, the approximation Levoy gives (Optics II, slide 39). Its first zero at $v = 3.83$ gives the radius $1.22\,\lambda N$ that ISETCam’s airyDisk.m returns and that Image Formation evaluates as 1.9, 5.4 and 10.7 µm at f/2.8, f/8 and f/16. Levoy (slides 36 and 40) summarizes the practical consequence: diffraction blur grows with $N$, and well-corrected lenses are called diffraction-limited. The example below checks that the Fourier transform of this pattern is $H_{\mathrm{diff}}$.

Pixel aperture. A pixel integrates the irradiance over its light-sensitive area; Levoy (Sampling and Pixels, slide 54) describes this prefilter as a pixel-sized rect. For a square aperture of side $a$ (equal to the pitch at 100% fill factor), the OTF is the transform of a box (our derivation):

$$H_{\mathrm{pix}}(\nu_x, \nu_y) = \operatorname{sinc}(a \nu_x)\, \operatorname{sinc}(a \nu_y), \qquad \operatorname{sinc}(u) = \frac{\sin \pi u}{\pi u}.$$

At the Nyquist frequency of a 100%-fill sensor, $a\nu = 0.5$, this is $2/\pi \approx 0.64$: the pixel alone removes about a third of the contrast there, and it still passes contrast above Nyquist, which aliases (Image Formation; sampling and aliasing).

Defocus. In geometric optics an out-of-focus point becomes a uniform disc, the circle of confusion, whose diameter $c$ the thin lens model derives. The OTF of a uniform disc of diameter $c$ is (our derivation)

$$H_{\mathrm{disc}}(\nu) = \frac{2 J_1(\pi c \nu)}{\pi c \nu},$$

which is negative between its first zero, at $\nu = 1.22/c$, and its second, at $2.23/c$, and changes sign again at each later zero. Image Deblurring notes that some PSFs have transfer functions that are exactly zero at certain frequencies; the disc is one, and between its zeros it reverses contrast. Multiplying $H_{\mathrm{diff}}$ by $H_{\mathrm{disc}}$ treats diffraction and geometric defocus as two independent blurs. Wave optics gives a different defocused OTF, which ISETCam’s opticsDefocusedMTF.m computes from a formula its comments attribute to Hopkins and to Marimont and Wandell; the product is a simpler model, and it is the one the example uses.

Aberrations (Levoy, Optics II, slides 3–23, and Image Formation) and motion during exposure (Image Deblurring) add further factors. Image Formation notes that distortion moves image points without blurring them, and Optical Microscopy that the PSF of a real objective is never better than the aberration-free ideal; uniform motion multiplies the OTF by a sinc with zeros.

Properties

• Products, not sums. The system MTF is the product of the component MTFs, so the weakest component limits it; doubling the pixel count does not help when the optics already cut off below the new Nyquist frequency.
• A hard cutoff. Diffraction sets $H = 0$ above $1/(\lambda N)$. A sensor whose Nyquist frequency lies above that cutoff cannot alias, and deconvolution, or super-resolution computed from the recorded images, cannot restore frequencies the optics removed (Optical Microscopy; Image Deblurring; Super-Resolution). Below the cutoff, deconvolution divides by $H$ and amplifies noise where it is small, which Wiener filtering regularizes.
• Normalization at zero frequency. The MTF is 1 at $\nu = 0$ by definition, so it carries no information about transmission or exposure.
• Phase. A symmetric PSF has a real OTF whose sign can change, as with defocus. An asymmetric PSF, such as coma, gives a phase that is not linear in frequency, shifting gratings by different amounts and making edges asymmetric; a phase linear in frequency only shifts the whole image (our statement).
• Shift variance. Real PSFs change across the field: Image Formation calls the single-PSF model a simplification, and Image Deblurring cites Levin et al. (2009) and Sage et al. (2017) on PSFs that vary over the image.
• Orientation. A rotationally symmetric PSF has the same MTF in every direction. The pixel aperture is not rotationally symmetric: its MTF along a diagonal differs from that along a row.
• Units. Values in cycles/mm, cycles/pixel and line pairs/mm are related by the pitch: with $p = 3$ µm, 1 cycle/pixel is 333 cycles/mm.

Measuring a PSF or an MTF

Point sources. A small enough source images as the PSF itself. In microscopy, a sub-resolution fluorescent bead serves, and a measured PSF differs from a theoretical one partly through the bead’s finite size and partly through the objective’s imperfections (Optical Microscopy, after Bankhead). In cameras, in our reading, a point source gives little light and its image covers few pixels, so its PSF is undersampled.

Edges. A straight edge gives far more light and samples the ESF along a whole line. The slanted-edge method, implemented in Peter Burns’s sfrmat programs and adapted in ISETCam’s ISO12233.m, measures one slice of the MTF from an image of an edge tilted a few degrees from the pixel columns. Both programs, read at the commit cited here, take the same steps:

1. Locate the edge: differentiate each row, take the centroid of each row’s derivative, and fit a line through the centroids. Burns’s sfrmat4 allows a polynomial of order 1 to 5 and notes that order 1, a straight line, is what the current ISO 12233 standard specifies (the standard itself was not read here).
2. Bin every pixel by its distance from the fitted edge along its row and average the pixel values in bins a quarter of a pixel wide (nbin = 4). Because the edge is slanted, successive rows cross it at different sub-pixel offsets, so the bins fill with samples at four times the pixel rate: a supersampled ESF.
3. Differentiate the ESF with a three-tap difference filter to obtain the LSF, and apply a Hamming window centred on the LSF.
4. Take the magnitude of the discrete Fourier transform and divide by its value at zero frequency.

The two programs differ in a detail that changes the result. sfrmat4 multiplies the MTF by a correction for the response of the difference filter (its fir2fix, limited to a factor of 10) and scales the sampling interval by the cosine of the edge angle, so that frequencies refer to the edge normal. ISETCam’s ISO12233.m defines fir2fix but does not call it and applies no cosine factor. A three-tap central difference with step $\delta$ multiplies a frequency $\nu$ by $\sin(2\pi\nu\delta)$ instead of the derivative’s $2\pi\nu\delta$ (our derivation), so without the correction the MTF reads low by the factor $\sin(2\pi\nu\delta)/(2\pi\nu\delta)$, 0.90 at the Nyquist frequency for $\delta = 1/4$ pixel. ISO12233.m also warns when the edge contrast is below 20%, which it says can lead to high error in the result.

The method measures the whole system as the image records it: optics, any anti-aliasing filter, the pixel aperture and, if it is applied, the camera’s processing. Sharpening, for example, raises the measured curve above the optics’ MTF (our reasoning: sharpening amplifies mid and high frequencies). It measures the MTF along one direction, near the edge’s position in the field.

Examples

The example checks the formulas above and measures an MTF by the slanted-edge method on a synthetic image whose true MTF is known exactly. It uses NumPy 2.3.5 and SciPy 1.16.3, has no random element, and runs in a few seconds.

Units are micrometres on the sensor and cycles/pixel or cycles/mm in frequency. The wavelength is 0.55 µm and the pixel pitch 3 µm, with 100% fill. Image coordinates are $x$ (column) and $y$ (row, downwards), with pixel centres at integer positions. The edge passes through the centre of a 64 × 64 pixel region, its normal at $\theta = 5°$ from the $x$ axis, so the edge is tilted 5° from the columns, and the bright side, value 1, is where the signed distance $d = (x - 31.5)\cos\theta + (y - 31.5)\sin\theta$ is positive; the dark side is 0. Values are linear in light, with no noise, clipping or quantization.

The script has three parts.

1. It samples the Airy pattern on a 2048 × 2048 grid at 0.5 µm spacing, takes the magnitude of its discrete Fourier transform and compares it along one axis with $H_{\mathrm{diff}}$, below and above the cutoff, at f/2.8, f/8 and f/16. The Airy pattern contains no frequencies above $\nu_c$, so 0.5 µm samples, finer than $1/(2\nu_c) = 0.77$ µm at f/2.8, lose nothing to aliasing (our reasoning); what remains is the truncation of its slowly decaying rings at the edge of the grid.
2. It checks the disc OTF formula against a numerical transform of the disc’s projection, a semicircle, for $c = 12$ µm.
3. It renders the slanted edge for four systems, f/2.8, f/8 and f/16 in focus and f/2.8 with a 12 µm (4 pixel) defocus disc, and measures their MTFs. The optics’ ESF is computed from its OTF, $\mathrm{ESF}(x) = \tfrac12 + \int_0^{\nu_c} H(\nu) \sin(2\pi\nu x)/(\pi\nu)\, d\nu$ (our derivation for a real, symmetric OTF), and each pixel value is the mean of that ESF over 16 × 16 sub-sample positions inside the pixel, so the pixel aperture enters by integration, not by its formula. The edge method then follows the steps above, except that the edge’s position and angle are taken from the construction rather than estimated, and pixels are binned by their distance along the edge normal, as the batch plan describes, so that no cosine factor is needed. It reports the MTF both without the difference-filter correction, as ISO12233.m computes it, and with a correction $t/\sin t$, $t = 2\pi\nu\delta$, which plays the role of sfrmat4’s fir2fix. It also reports a signed estimate: because the edge lies at the centre of the bin array, removing the phase of that half-array shift leaves the real, signed OTF. The reference is the analytic product $H_{\mathrm{diff}}\, H_{\mathrm{disc}}\, \operatorname{sinc}(\nu\cos\theta)\operatorname{sinc}(\nu\sin\theta)$ along the edge normal, with $\nu$ in cycles/pixel.

The code re-implements the slanted-edge steps from their description; it does not port the sfrmat code.

import os
for v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
    os.environ[v] = "1"
import numpy as np
from scipy.special import j1
from scipy.interpolate import CubicSpline

lam, p = 0.55, 3.0                 # wavelength and pixel pitch (um)
theta = np.deg2rad(5.0)            # edge normal, 5 degrees from the x (column) axis
n, ss, nbin = 64, 16, 4            # ROI size (px), sub-samples per pixel side, bins per px

def h_diff(nu, nu_c):              # diffraction-limited OTF of a circular pupil
    s = np.clip(np.abs(nu) / nu_c, 0, 1)
    return (2 / np.pi) * (np.arccos(s) - s * np.sqrt(1 - s**2))

def disc_otf(nu, c):               # uniform disc of diameter c: 2 J1(pi c nu) / (pi c nu)
    a = np.pi * c * np.asarray(nu, float)
    out = np.ones_like(a)
    out[a != 0] = 2 * j1(a[a != 0]) / a[a != 0]
    return out

def optics_esf(N, c):              # ESF(x) = 1/2 + int OTF(nu) sin(2 pi nu x) / (pi nu) dnu
    nu_c = p / (lam * N)           # cutoff 1/(lambda N), in cycles/pixel
    dnu = 0.02 * p / 1000          # 0.02 cycles/mm
    nu = np.arange(0, nu_c + dnu / 2, dnu)
    otf = h_diff(nu, nu_c) * (disc_otf(nu, c) if c else 1.0)
    w = np.full(nu.size, dnu); w[[0, -1]] = dnu / 2
    x = np.round(np.arange(-4000, 4001) * 0.01, 2)          # -40 ... 40 px
    nz = np.where(nu > 0, nu, 1.0)
    esf = np.empty_like(x)
    for i in range(0, x.size, 400):
        xx = x[i:i + 400, None]
        k = np.where(nu > 0, np.sin(2 * np.pi * nu * xx) / (np.pi * nz), 2 * xx)
        esf[i:i + 400] = 0.5 + k @ (otf * w)
    return x, esf

yy, xx = np.mgrid[0:n, 0:n].astype(float)
def dist(ox=0.0, oy=0.0):          # signed distance (px) to the edge through the ROI centre
    return ((xx + ox - (n - 1) / 2) * np.cos(theta)
            + (yy + oy - (n - 1) / 2) * np.sin(theta))

def render(N, c):                  # pixel = mean of the optics ESF over 16 x 16 sub-samples
    x, esf = optics_esf(N, c)
    spl = CubicSpline(x, esf)
    off = (np.arange(ss) + 0.5) / ss - 0.5
    img = sum(spl(dist(ox, oy)) for oy in off for ox in off) / ss**2
    return img, esf

def edge_mtf(img):                 # edge method with the edge known by construction
    nb, delta = n * nbin, 1 / nbin
    k = np.floor((dist().ravel() + n / 2) / delta).astype(int)
    keep = (k >= 0) & (k < nb)
    cnt = np.bincount(k[keep], minlength=nb)
    esf = np.bincount(k[keep], img.ravel()[keep], minlength=nb) / cnt
    lsf = np.empty(nb)
    lsf[1:-1] = (esf[2:] - esf[:-2]) / 2                     # 3-tap derivative
    lsf[[0, -1]] = lsf[[1, -2]]
    F = np.fft.rfft(lsf * np.hamming(nb))
    nu = np.arange(F.size) / (nb * delta)                    # cycles/pixel
    t = 2 * np.pi * nu * delta
    corr = np.ones_like(nu); corr[1:] = t[1:] / np.sin(t[1:])   # derivative response
    shift = np.exp(2j * np.pi * np.arange(F.size) * (nb - 1) / 2 / nb)
    mtf_u = np.abs(F) / abs(F[0])
    signed = np.real(F * shift) / np.real(F[0]) * corr
    return nu, mtf_u, mtf_u * corr, signed, cnt

def system_otf(nu, N, c):          # analytic product along the edge normal (signed)
    D = disc_otf(nu, c) if c else 1.0
    pix = np.sinc(nu * np.cos(theta)) * np.sinc(nu * np.sin(theta))
    return h_diff(nu, p / (lam * N)) * D * pix

# 1. |FFT| of a finely sampled Airy PSF against the analytic MTF
M, s = 2048, 0.5                   # grid size and spacing (um)
r = np.hypot(*np.meshgrid((np.arange(M) - M // 2) * s, (np.arange(M) - M // 2) * s))
print("Airy PSF |FFT| vs analytic MTF   max err below cutoff   max |MTF| above")
for N in (2.8, 8, 16):
    v = np.pi * r / (lam * N)
    I = np.ones_like(v); I[v > 0] = (2 * j1(v[v > 0]) / v[v > 0])**2
    F = np.abs(np.fft.fft2(np.fft.ifftshift(I)))[0, :M // 2]
    mtf, nu = F / F[0], np.arange(M // 2) / (M * s)          # cycles/um
    nu_c = 1 / (lam * N)
    lo, hi = nu < nu_c, nu > nu_c
    print(f"  f/{N:<4} cutoff {1000 * nu_c:6.1f} cy/mm       "
          f"{np.max(np.abs(mtf[lo] - h_diff(nu[lo], nu_c))):.2e}"
          f"              {np.max(mtf[hi]):.2e}")

# 2. Disc OTF formula: 1-D transform of the disc's projection (a semicircle)
R = 6.0                            # disc radius (um), c = 12 um
xs = np.arange(-6000, 6001) * 0.001
lsf_d = 2 * np.sqrt(np.clip(R**2 - xs**2, 0, None)) / (np.pi * R**2)
nus = np.arange(334) / 1000.0      # 0 ... 0.333 cycles/um (0 ... 1 cy/px)
num = np.trapezoid(lsf_d[None, :] * np.cos(2 * np.pi * nus[:, None] * xs[None, :]), xs, axis=1)
print(f"disc OTF formula: max |numerical - 2J1(x)/x| = "
      f"{np.max(np.abs(num - disc_otf(nus, 2 * R))):.2e}")

# 3. Slanted-edge MTF against the analytic system MTF
conds = [(2.8, 0), (8, 0), (16, 0), (2.8, 4)]
at_nyq = {}
for N, c in conds:
    img, esf_tab = render(N, c)
    nu, mu, mc, sg, cnt = edge_mtf(img)
    A = system_otf(nu, N, c)
    b = nu <= 0.5
    eu, ec = np.abs(mu - np.abs(A))[b], np.abs(mc - np.abs(A))[b]
    print(f"\nf/{N}" + (f" + {c} px defocus disc" if c else " in focus")
          + f"   bins: {cnt.min()}-{cnt.max()} px each, ESF(0) = {esf_tab[4000]:.9f}, "
          f"ESF(-40) = {esf_tab[0]:.4f}, ESF(40) = {esf_tab[-1]:.4f}, "
          f"pixels in [{img.min():.10f}, {img.max():.10f}]")
    print("  nu(cy/px) cy/mm  analytic  corrected  uncorr" + ("  signed: analytic  edge" if c else ""))
    for j in range(0, 33, 4):
        line = f"   {nu[j]:.4f} {nu[j] * 1000 / p:6.1f}   {abs(A[j]):.4f}    {mc[j]:.4f}   {mu[j]:.4f}"
        print(line + (f"           {A[j]:+.4f}  {sg[j]:+.4f}" if c else ""))
    print(f"  max error 0-0.5 cy/px: corrected {ec.max():.2e} (at {nu[b][ec.argmax()]:.4f} cy/px), "
          f"uncorrected {eu.max():.2e}")
    if c:
        neg = b & (A < -0.038)
        print(f"  frequencies with analytic OTF < -0.038: {neg.sum()}, signed edge < 0: "
              f"{(sg[neg] < 0).sum()}, |edge| > 0: {(mc[neg] > 0).sum()}")
    else:
        at_nyq[N] = mc[32]
print("\ncorrected MTF at Nyquist:", ", ".join(f"f/{N} {v:.4f}" for N, v in at_nyq.items()))

Airy PSF |FFT| vs analytic MTF   max err below cutoff   max |MTF| above
  f/2.8  cutoff  649.4 cy/mm       5.75e-04              1.13e-06
  f/8    cutoff  227.3 cy/mm       1.64e-03              7.08e-06
  f/16   cutoff  113.6 cy/mm       3.27e-03              2.05e-05
disc OTF formula: max |numerical - 2J1(x)/x| = 8.05e-07

f/2.8 in focus   bins: 8-18 px each, ESF(0) = 0.500000000, ESF(-40) = 0.0008, ESF(40) = 0.9992, pixels in [0.0009703410, 0.9990296590]
  nu(cy/px) cy/mm  analytic  corrected  uncorr
   0.0000    0.0   1.0000    1.0000   1.0000
   0.0625   20.8   0.9530    0.9580   0.9565
   0.1250   41.7   0.8949    0.8983   0.8925
   0.1875   62.5   0.8278    0.8285   0.8166
   0.2500   83.3   0.7536    0.7513   0.7321
   0.3125  104.2   0.6747    0.6690   0.6425
   0.3750  125.0   0.5933    0.5844   0.5512
   0.4375  145.8   0.5115    0.4998   0.4613
   0.5000  166.7   0.4312    0.4171   0.3755
  max error 0-0.5 cy/px: corrected 1.41e-02 (at 0.5000 cy/px), uncorrected 5.57e-02

f/8 in focus   bins: 8-18 px each, ESF(0) = 0.500000000, ESF(-40) = 0.0024, ESF(40) = 0.9976, pixels in [0.0027730306, 0.9972269694]
  nu(cy/px) cy/mm  analytic  corrected  uncorr
   0.0000    0.0   1.0000    1.0000   1.0000
   0.0625   20.8   0.8778    0.8920   0.8906
   0.1250   41.7   0.7483    0.7592   0.7543
   0.1875   62.5   0.6171    0.6241   0.6151
   0.2500   83.3   0.4896    0.4927   0.4801
   0.3125  104.2   0.3706    0.3702   0.3556
   0.3750  125.0   0.2642    0.2620   0.2471
   0.4375  145.8   0.1736    0.1706   0.1575
   0.5000  166.7   0.1010    0.0980   0.0883
  max error 0-0.5 cy/px: corrected 1.56e-02 (at 0.0312 cy/px), uncorrected 1.71e-02

f/16 in focus   bins: 8-18 px each, ESF(0) = 0.500000000, ESF(-40) = 0.0047, ESF(40) = 0.9953, pixels in [0.0055450580, 0.9944549420]
  nu(cy/px) cy/mm  analytic  corrected  uncorr
   0.0000    0.0   1.0000    1.0000   1.0000
   0.0625   20.8   0.7630    0.7888   0.7875
   0.1250   41.7   0.5300    0.5473   0.5438
   0.1875   62.5   0.3177    0.3268   0.3221
   0.2500   83.3   0.1427    0.1450   0.1413
   0.3125  104.2   0.0242    0.0220   0.0212
   0.3750  125.0   0.0000    0.0038   0.0036
   0.4375  145.8   0.0000    0.0022   0.0021
   0.5000  166.7   0.0000    0.0000   0.0000
  max error 0-0.5 cy/px: corrected 3.00e-02 (at 0.0312 cy/px), uncorrected 2.96e-02

f/2.8 + 4 px defocus disc   bins: 8-18 px each, ESF(0) = 0.500000000, ESF(-40) = 0.0008, ESF(40) = 0.9992, pixels in [0.0009711799, 0.9990288201]
  nu(cy/px) cy/mm  analytic  corrected  uncorr  signed: analytic  edge
   0.0000    0.0   1.0000    1.0000   1.0000           +1.0000  +1.0000
   0.0625   20.8   0.8814    0.8865   0.8851           +0.8814  +0.8865
   0.1250   41.7   0.6459    0.6492   0.6450           +0.6459  +0.6492
   0.1875   62.5   0.3719    0.3729   0.3675           +0.3719  +0.3729
   0.2500   83.3   0.1366    0.1353   0.1319           +0.1366  +0.1353
   0.3125  104.2   0.0130    0.0157   0.0151           -0.0130  -0.0157
   0.3750  125.0   0.0709    0.0738   0.0696           -0.0709  -0.0738
   0.4375  145.8   0.0635    0.0652   0.0602           -0.0635  -0.0652
   0.5000  166.7   0.0291    0.0287   0.0258           -0.0291  -0.0287
  max error 0-0.5 cy/px: corrected 5.60e-03 (at 0.0312 cy/px), uncorrected 5.22e-03
  frequencies with analytic OTF < -0.038: 10, signed edge < 0: 10, |edge| > 0: 10

corrected MTF at Nyquist: f/2.8 0.4171, f/8 0.0980, f/16 0.0000

The numerical transform of the Airy pattern agrees with ISETCam’s formula to within 0.0033 below the cutoff at every f-number, and is below $3 \times 10^{-5}$ above it, so the closed-form OTF and the Airy PSF describe the same optics, and the cutoff at $1/(\lambda N)$ is real: 649, 227 and 114 cycles/mm. The error grows with $N$ because, in our reading, the same grid holds fewer rings of the wider pattern, so more of its light lies outside. The disc formula agrees with the transform of its projection to $8 \times 10^{-7}$.

The slanted-edge MTF with the difference-filter correction stays within about 0.03 of the true system MTF up to the Nyquist frequency in all four systems. Where its largest error falls says what limits it. At f/2.8 the largest error, 0.014, is at the Nyquist frequency, where the corrected estimate reads 0.4171 against 0.4312. In our reading this is mostly the quarter-pixel binning, which averages the ESF over each bin and so multiplies the measured MTF by $\operatorname{sinc}(\nu\delta) = 0.974$ at Nyquist (our derivation): $0.4312 \times 0.974 = 0.420$. At f/8, f/16 and with defocus, the largest error is at the lowest non-zero frequency, 0.0312 cycles/pixel, and the printed rows at 0.0625 cycles/pixel show the estimate reading high there: 0.8920 against 0.8778 at f/8 and 0.7888 against 0.7630 at f/16. In our reading this is the Airy pattern’s long tail. The measurement sees the ESF only over ±32 pixels, and the printed ESF table values at ±40 pixels, 0.0047 and 0.9953 at f/16, show that part of the light is spread farther still. The light missing from the measured tail would lower the MTF at the lowest frequencies, so cutting it off raises them. At f/16 the estimate also shows 0.0038 and 0.0022 at 125 and 146 cycles/mm, above the optics’ cutoff, where the true MTF is zero; in our reading this is leakage from the window and the truncated LSF.

Without the correction, the f/2.8 estimate falls to 0.3755 at Nyquist, 13% low, and its largest error, 0.056, is four times that of the corrected version; at f/8 the uncorrected error is 0.017 against 0.016. At f/16 and with defocus the uncorrected error is slightly the smaller, 0.0296 against 0.0300 and 0.0052 against 0.0056, as printed. In our reading, the corrected estimate’s largest error there lies at 0.0312 cycles/pixel, where the correction factor is close to 1, and where the estimate reads high, as the rows at 0.0625 cycles/pixel show, the slightly lower uncorrected value is the closer. Of the four systems, the choice between the two codes matters most at f/2.8, where the MTF is still high near Nyquist.

The MTF at the Nyquist frequency falls with the f-number, from 0.4171 at f/2.8 to 0.0980 at f/8 and 0.0000 at f/16. At f/16 the diffraction cutoff, 114 cycles/mm, lies below Nyquist, so the optics pass no contrast that could alias; at f/2.8 the image keeps 42% of the contrast of a Nyquist-frequency grating, and frequencies just above Nyquist alias. Image Formation’s arithmetic makes the same point in the spatial domain: the Airy radius is 1.9 µm at f/2.8 and 10.7 µm at f/16, against a 3 µm pitch.

With the 12 µm defocus disc, the true OTF passes through zero at 0.305 cycles/pixel (102 cycles/mm) and is negative from there to Nyquist. The MTF column rises again after the zero, from 0.0157 at 104 cycles/mm to 0.0738 at 125 cycles/mm, as if contrast were recovering; the signed estimate shows that it is reversed, −0.0738, in agreement with the analytic −0.0709. At all 10 frequencies where the true OTF is below −0.038, the signed estimate is negative while the MTF is positive. A defocused photograph of a fine grating at such a frequency shows bars of the right spacing in the wrong places, and an MTF curve that bounces after a zero is a sign of this reversal, not of recovered detail.

Planned outcomes. All values are from the printed output, and every hypothesis and check was written in the plan before any script ran.

• Hypothesis H1, that the transform of the sampled Airy pattern matches $H_{\mathrm{diff}}$ below the cutoff and is near zero above it within a tolerance of 0.0041, held at f/2.8, f/8 and f/16.
• Hypothesis H2, that the corrected slanted-edge MTF matches the analytic system MTF from 0 to 0.5 cycles/pixel within 0.038, held in all four systems.
• Hypothesis H2u, that the uncorrected estimate has a larger maximum error than the corrected one at f/2.8 and f/8, held at both. At f/16 and with defocus, outside its scope, the uncorrected error was the smaller, as printed.
• Hypothesis H3, that the MTF at Nyquist falls from f/2.8 to f/8 to f/16 and is within 0.038 of zero at f/16, held.
• Hypothesis H4, that with the defocus disc the signed estimate is negative and the MTF positive wherever the true OTF is below −0.038, held at all 10 such frequencies.
• The check of the disc OTF formula (difference at most $10^{-4}$) held.
• The data checks held: every quarter-pixel bin received 8 to 18 pixels, the ESF is 0.5 at the edge, and every pixel value lies between 0 and 1.

Process. The tolerances 0.0041 and 0.038 were set from a planning run of the same noise-free computation, made after the hypotheses were written and before the pilot; the plan records that “Because the tolerances below come from the same noise-free computation that the pilot repeats, H1 and H2 are illustrations of the size of each method’s error, not tests”. That planning run also printed the errors that decide H2u, so H2u is an illustration as well. H3 and H4 were written before any run, but use the same 0.038 as their threshold; since H3’s analytic values at Nyquist (about 0.43, 0.10 and 0, computed by hand in the plan) are far apart, the planning run’s largest errors already implied that H3 would hold, so H3 is an illustration as well. H4’s signs and the checks rest on predictions written before any run. The example departs from the batch plan’s description in six ways, each decided before the pilot:

• the edge’s position and angle come from the construction instead of being estimated;
• the MTF is reported both with and without the difference-filter correction, because sfrmat4 and ISETCam’s ISO12233.m differ;
• defocus is a geometric disc multiplied with the diffraction OTF, not the wave-optics defocused OTF of opticsDefocusedMTF.m;
• a signed estimate is added alongside the MTF to show the contrast reversal;
• the disc-formula check and the data checks were added;
• the scene is rendered from the analytic OTF through the ESF integral, not by spatial convolution with a sampled Airy pattern, whose slowly decaying rings a finite kernel would cut off; the Airy check in the first part ties the two descriptions together.

One edit was made to the script before its first run: the printed range of pixel values was given more decimals so that the data check could be read from the output. The pilot was approved without change, and the run shown is identical to it.

Common Misconceptions

• “The MTF is the PSF’s width in another form.” The MTF depends on the whole shape of the PSF, not only on its width (our statement from the definition). Levoy’s two MTF curves (slide 18) differ in shape, not in a single resolution figure, and a defocus disc has an MTF that falls to zero and rises again.
• “An MTF that rises again after a zero means detail comes back.” For a symmetric PSF it means the OTF has changed sign: the detail is there with its contrast reversed, as the defocus case above shows.
• “Higher-resolution sensors always give sharper images.” The system MTF is a product. When the optics’ cutoff $1/(\lambda N)$ lies below the sensor’s Nyquist frequency, smaller pixels add samples but no contrast (our reasoning from the product rule; Levoy, Optics II, slides 38–39, discusses diffraction blur against pixel size).
• “The edge method measures the lens.” It measures everything between the edge and the pixel values, including the pixel aperture, demosaicing and sharpening, along one direction at one place in the field.
• “The contrast of a bar chart is the MTF.” Bars are square waves with harmonics; the MTF is defined for sinusoids. Image Formation’s example measures bar contrast, a related but different quantity.
• “One PSF describes the whole image.” Shift invariance is a model. Real PSFs vary across the field, and in microscopy with depth (Image Deblurring).

Where It Is Used

• Camera and lens characterization. MTF curves describe lens sharpness in lens reviews (Levoy, Sampling and Pixels, slides 2 and 17). Multispectral satellite imaging reports Sentinel-2’s MTF requirement at the Nyquist frequency.
• System design. Choosing pixel pitch, f-number and anti-aliasing filters is a matter of balancing the component MTFs against the Nyquist frequency; vision system design and digital camera cover the practical side.
• Restoration. Image deblurring, Wiener filtering and super-resolution model the image as the scene convolved with a PSF and invert the model where the transfer function is not too small.
• Microscopy. Optical microscopy and fluorescence microscopy derive resolution limits from the PSF and measure it with beads; the PSF of a 3D microscope, an ultrasound pulse-echo PSF (ultrasound imaging) and the PSF of an MRI sampling mask (medical image reconstruction) are PSFs of other systems in the same sense.
• Simulation. Rendering training or test images with a known PSF, as the examples of image formation and image deblurring do, gives ground truth for sharpness and restoration.