Measure a point spread function from beads and deconvolve with it (Python)

The problem

You have an image of sub-resolution fluorescent beads and one of a sample, taken with the same objective, filters, and camera. You want a point spread function (PSF) from the beads, the sample deconvolved with it, and a check that it fits. A real PSF is close to an Airy pattern: a core with its first dark ring at 0.61 λ/NA, 0.23 µm at 520 nm and NA 1.4, and faint rings beyond. The peaks and crops come from Counting cells with scipy.ndimage, the residual and the stopping rule from Richardson-Lucy deconvolution with scikit-image.

The example is a generated bead slide and a sample with two beads 0.15 µm apart. Swap in your own images, pixel size, wavelength, and NA; every length in the code follows from the last three.

The code

The section marked # ---- images generates both pictures and is the only one you replace. Your images must be in detected photons, because a camera stores an offset plus a gain times the photoelectrons. OFFSET is in counts and GAIN in counts per photoelectron, one over the electrons per count a data sheet gives. FLOOR is the residual a right PSF reaches, 1 when photon noise is all the noise there is.

import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage as ndi, special
from scipy.signal import fftconvolve
from skimage.restoration import richardson_lucy

# ---- your numbers
PIXEL, WAVELENGTH, NA = 0.05, 0.52, 1.4          # µm per pixel, µm, numerical aperture
OFFSET, GAIN = 0.0, 1.0                           # counts = OFFSET + GAIN * photoelectrons; GAIN = 1 / (e⁻ per count)
FLOOR, COUNTS = 1.0, (25, 50, 100, 200, 400)      # residual of a right PSF (1: photon noise alone); iterations
RING = 0.61 * WAVELENGTH / NA / PIXEL             # first dark ring of the Airy pattern, in pixels
HALF = round(4 * RING)                            # crop half-width: holds the first three bright rings

# ---- images: generated; replace with calib = ..., sample = ... (counts) and VIEW = the region to show
rng = np.random.default_rng(342)
v = lambda r: 2 * np.pi * NA * PIXEL / WAVELENGTH * np.maximum(r, 1e-9)
airy = lambda r: (2 * special.j1(v(r)) / v(r)) ** 2 * v(1) ** 2 / (4 * np.pi)   # photons per pixel, 1 in total
row, col = np.indices((320, 320))
calib = 10 + 10 * col / 319                       # uneven illumination, 10 to 20 photons per pixel
for r0, c0 in rng.uniform(0, 320, (30, 2)):
    calib = calib + 50_000 * airy(np.hypot(row - r0, col - c0))
truth = np.zeros((96, 96))
truth[48, [47, 50]] = truth[[20, 20, 76, 76], [20, 76, 20, 76]] = 20_000
blurred = fftconvolve(truth, airy(np.hypot(*np.indices((121, 121)) - 60)), mode="same")
calib, sample, VIEW = rng.poisson(calib), rng.poisson(np.clip(blurred, 0, None) + 10), np.s_[44:53, 43:55]
calib, sample = (calib - OFFSET) / GAIN, (sample - OFFSET) / GAIN   # <- in photons; the generated images already are

# ---- find isolated beads and average them
def edge_median(a, w):                            # median of the strip w pixels wide along the border
    return np.median(np.concatenate([a[:w], a[-w:], a[:, :w].T, a[:, -w:].T], axis=None))
smooth = ndi.gaussian_filter(calib, RING / 4)
rows, cols = np.nonzero((smooth == ndi.maximum_filter(smooth, size=2 * round(RING) + 1)) & (smooth > 0.2 * smooth.max()))
def bead_crops(half=HALF):                        # a neighbor's core within half + RING enters the crop
    full = (rows >= half) & (cols >= half) & (rows < calib.shape[0] - half) & (cols < calib.shape[1] - half)
    alone = np.sum(np.maximum(abs(rows - rows[:, None]), abs(cols - cols[:, None])) <= half + round(RING), axis=1) == 1
    crops = [calib[r - half:r + half + 1, c - half:c + half + 1] for r, c in zip(rows[full & alone], cols[full & alone])]
    return crops, np.sum(~full), np.sum(full & ~alone)
def average(crops):                               # each crop minus the background where its bead sits
    psf = np.clip(np.mean([c - edge_median(c, 3) for c in crops], axis=0), 0, None)
    return psf / psf.sum()
crops, border, close = bead_crops()
psf = average(crops)

# ---- deconvolve with the PSF and with a Gaussian of the Airy core's FWHM, 0.51 λ/NA; stop by the discrepancy rule
fwhm, x = 0.51 * WAVELENGTH / NA, np.arange(-HALF, HALF + 1)
gauss = np.exp(-(x[:, None] ** 2 + x[None, :] ** 2) * (2.355 * PIXEL / fwhm) ** 2 / 2)
gauss /= gauss.sum()
background = edge_median(sample, 4)
image = np.clip(sample - background, 0, None)
def deconvolve(p, n):
    return richardson_lucy(np.pad(image, h := p.shape[0] // 2, mode="symmetric"), p, num_iter=n, clip=False)[h:-h, h:-h]
def residual(est, p):
    model = fftconvolve(np.pad(est, p.shape[0] // 2, mode="symmetric"), p, mode="valid") + background
    return np.mean((sample - model) ** 2 / model)
def ring_light(p):                                # fraction of the light outside the first dark ring
    return p[np.hypot(*(np.indices(p.shape) - p.shape[0] // 2)) >= RING].sum() / p.sum()
def stop(p):                                      # first count with residual <= FLOOR (else the last), all residuals
    res = {n: residual(deconvolve(p, n), p) for n in COUNTS}
    return next((n for n in COUNTS if res[n] <= FLOOR), COUNTS[-1]), res
(n_stop, res_psf), (_, res_gauss) = stop(psf), stop(gauss)
estimates = {"Gaussian": deconvolve(gauss, n_stop), "measured PSF": deconvolve(psf, n_stop)}

# ---- report and plot
print(f"{len(rows)} peaks, {len(crops)} beads kept, {border} rejected at the border, {close} for a neighbor; stop at {n_stop}")
for name, p, res in [("measured PSF", psf, res_psf), (f"Gaussian, FWHM {fwhm:.2f} µm", gauss, res_gauss)]:
    print(f"{name:26s} ring light {100 * ring_light(p):4.1f} %   residual {res[n_stop]:.2f}, at 400 {res[400]:.2f}")
fig, axs = plt.subplots(1, 3, figsize=(8, 2.9), layout="constrained", width_ratios=[1.5, 1, 1])
for p, color, marker, name, xy in [(psf, "#c8553d", "o", "measured", (-0.62, 3e-2)), (gauss, "#2a7f9e", "", "Gaussian", (0, 3e-4))]:
    axs[0].plot(x * PIXEL, np.where(p[HALF] > 0, p[HALF] / p.max(), np.nan), marker=marker, color=color, ms=3, lw=1.2)
    axs[0].text(*xy, name, color=color, ha="center")
axs[0].vlines([-RING * PIXEL, RING * PIXEL], 1e-4, 1.2, color="#8a8f98", ls="--", lw=1)   # first dark ring
axs[0].set(yscale="log", ylim=(1e-4, 1.2), xlabel="x / µm", ylabel="intensity / peak")
axs[0].spines[["top", "right"]].set_visible(False)
for ax, (name, est), res in zip(axs[1:], estimates.items(), [res_gauss, res_psf]):
    ax.imshow(est[VIEW], cmap="gray", vmin=0, vmax=max(e[VIEW].max() for e in estimates.values()))
    ax.text(0, 1.03, f"{name}, residual {res[n_stop]:.2f}", transform=ax.transAxes, color="#1f2a44")
    ax.set_axis_off()
axs[2].plot([0.5, 0.5 + 0.25 / PIXEL], [7.6, 7.6], color="white", lw=3)   # scale bar
axs[2].text(0.5 + 0.125 / PIXEL, 7.1, "0.25 µm", color="white", ha="center")
plt.show()

30 peaks, 20 beads kept, 8 rejected at the border, 2 for a neighbor; stop at 100
measured PSF               ring light 11.7 %   residual 0.98, at 400 0.95
Gaussian, FWHM 0.19 µm     ring light  1.3 %   residual 1.58, at 400 1.54

Left: center row of the measured PSF (dots) and a Gaussian of the same core width (line), intensity relative to the peak on a log axis against x in µm; only the measured PSF has rings past the dashed first dark ring at 0.23 µm. Middle: the pair deconvolved with the Gaussian, residual 1.58. Right: with the measured PSF, residual 0.98.

Of the 30 peaks, 20 are isolated beads; 8 sat too close to the edge for a full crop and 2 had a neighbor. Ring light, the fraction of a PSF's light outside its first dark ring, is 11.7 % for the measured PSF and 1.3 % for the Gaussian. The residual tells the two apart without a truth. With the measured PSF it falls to 0.98 at 100 iterations; with the Gaussian it is 1.58 there and levels off at 1.54 by 400. The Gaussian does not have the shape of this PSF.

The next cell checks both against the generated truth, turns three knobs, each line a new PSF judged by the residual alone, and compares the two cores:

for name, est in estimates.items():               # generated data only: dip between the pair, error against the truth
    dip, err = 1 - est[48, 48:50].min() / est[48, [47, 50]].min(), np.sqrt(np.mean((est - truth) ** 2)) / truth.max()
    print(f"{name:12s} dip {100 * dip:3.0f} %   RMS error {100 * err:.2f} % of the peak")

def report(label, crops):
    p = average(crops)
    n, res = stop(p)
    print(f"{label:22s} {len(crops):3d} beads   ring light {100 * ring_light(p):4.1f} %   "
          f"residual at 400 {res[400]:.2f}   stop {n if res[n] <= FLOOR else 'never'}")
for m in (1, 2, 3, 4, 6):
    report(f"half-width {m} RING", bead_crops(round(m * RING))[0])
for n_beads in (1, 3, 10):
    report(f"first {n_beads} of {len(crops)}", crops[:n_beads])
report("neighbors kept", [calib[r - HALF:r + HALF + 1, c - HALF:c + HALF + 1] for r, c in zip(rows, cols)
                          if HALF <= min(r, c) and max(r, c) < calib.shape[0] - HALF])

def core(p):                                      # share of the core's light that lies within half the ring radius
    r = np.hypot(*(np.indices(p.shape) - p.shape[0] // 2))
    return p[r < RING / 2].sum() / p[r < RING].sum()
print(f"core light within half a ring radius: measured PSF {100 * core(psf):.0f} %, Gaussian {100 * core(gauss):.0f} %")

Gaussian     dip  43 %   RMS error 0.99 % of the peak
measured PSF dip  73 %   RMS error 0.70 % of the peak
half-width 1 RING       30 beads   ring light  0.4 %   residual at 400 1.04   stop never
half-width 2 RING       25 beads   ring light  6.9 %   residual at 400 0.95   stop 100
half-width 3 RING       24 beads   ring light 10.1 %   residual at 400 0.95   stop 100
half-width 4 RING       20 beads   ring light 11.7 %   residual at 400 0.95   stop 100
half-width 6 RING        7 beads   ring light 14.9 %   residual at 400 1.06   stop never
first 1 of 20            1 beads   ring light 12.8 %   residual at 400 1.03   stop never
first 3 of 20            3 beads   ring light 12.2 %   residual at 400 1.00   stop 200
first 10 of 20          10 beads   ring light 11.8 %   residual at 400 0.99   stop 200
neighbors kept          22 beads   ring light 12.9 %   residual at 400 1.04   stop never
core light within half a ring radius: measured PSF 80 %, Gaussian 74 %

The truth agrees: the pair dips 73 % between the beads with the measured PSF against 43 % with the Gaussian, and the RMS error is 0.70 % against 0.99 % of the peak.

The knobs

The crop size decides how much of the rings the PSF keeps. A half-width of one radius cuts every ring off: ring light drops to 0.4 %, and the residual levels at 1.04, not at the Gaussian's 1.54. The rings are the lesser part of the Gaussian's error. The rest is the core: the measured one holds 80 % of its light within half the ring radius, the Gaussian 74 %, whose skirt still glows where the Airy core has gone dark.

Six radii leave only 7 isolated beads, and the residual climbs back to 1.06 through noise: a noisy PSF is a wrong PSF too. Take three to four ring radii, as long as a dozen isolated beads still fit. The noise of the average falls as one over the square root of the number of beads: one bead never reaches 1, three and ten reach it at 200 iterations, and all twenty at 100.

Kept beads with a neighbor hold the residual at 1.04: the neighbor's core enters the average as a false ring. Two beads closer than the core pass as one bead twice as bright, so drop a kept bead far brighter than the rest. The threshold, 0.2 of the brightest smoothed peak, misses single beads once a clump is five times brighter; lower it then. Centering on the smoothed peak pixel shifts each crop by up to half a pixel, 0.025 µm against a 0.19 µm core.

The averaged PSF belongs to this objective and wavelength at the beads' depth and focus, convolved with the bead and the half-pixel jitter. Beads of 100 nm or less are small against the 0.19 µm core; larger beads widen it. It is a 2D in-focus PSF: a sample deep in a medium of another refractive index has a wider, lopsided one, and a z-stack needs a 3D stack of beads.

A residual at FLOOR says the PSF explains the sample down to the noise. Read noise σ_r, in electrons on the data sheet, lifts the floor to about 1 + σ_r²/b, with b the background in photons per pixel: 1.11 for 1.5 e⁻ over 20 photons. An EMCCD doubles the variance: multiply GAIN by its multiplication gain and set FLOOR to 2. A residual above the floor at every count says the PSF is wrong. The discrepancy stop is the safe, early one, not the best count.

Pitfalls

The numbers below come from this cell:

est = estimates["measured PSF"]
print(f"PSF summing to 1.2: residual {residual(est, 1.2 * psf):.2f}; "
      f"PSF times 3: total {deconvolve(3 * psf, n_stop).sum() / est.sum():.3f} of the run with the PSF as is")
raw = np.mean(crops, axis=0)                      # the same crops, background left in
raw /= raw.sum()
print(f"crops not background-subtracted: ring light {100 * ring_light(raw):.1f} %, residual at 400 {stop(raw)[1][400]:.2f}")
print("crops near the edge:", calib[8 - HALF:8 + HALF + 1, 200 - HALF:200 + HALF + 1].shape,
      calib[200 - HALF:200 + HALF + 1, 312 - HALF:312 + HALF + 1].shape)

PSF summing to 1.2: residual 1.34; PSF times 3: total 1.000 of the run with the PSF as is
crops not background-subtracted: ring light 38.9 %, residual at 400 3.15
crops near the edge: (0, 37) (37, 26)

A PSF that does not sum to 1. richardson_lucy does not care: the PSF times 3 gives the same estimate, with a total of 1.000 of the run with the PSF as is. Everything that blurs the estimate again does care. A PSF summing to 1.2 puts the residual at 1.34 instead of 0.98, and the discrepancy rule never stops. Divide by the sum after averaging, every time.

Background left in the crops. Average the raw crops and the 10 to 20 photons per pixel of background become a flat pedestal under the PSF. Ring light jumps from 11.7 % to 38.9 %, and the residual levels at 3.15, nowhere near 1: Richardson-Lucy reads the pedestal as a wide blur. Subtract each crop's border median, and check that the averaged PSF's edge sits at about zero.

Beads near the border. A negative slice start does not wrap around, it counts from the end. A bead 8 pixels from the top edge gives calib[-10:27], an empty crop of shape (0, 37), and a bead 8 pixels from the right edge a short crop of shape (37, 26), because a slice past the end stops at the edge. Averaging crops of unequal shape then fails with ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape. The code skips such beads, and the same holds for a bead cut by the edge of a sample tile.