Image Deblurring

Image deblurring recovers a sharp image from a blurred one. Blur appears whenever the light from one point of the scene spreads over several pixels: a lens out of focus, a camera or an object moving during the exposure, diffraction in a microscope objective. When the spread is the same everywhere, the blurred image is the convolution of the sharp image with the point-spread function (PSF), and deblurring means undoing that convolution, which is why it is also called deconvolution. The difficulty is that blur discards information. Fine detail is attenuated, some of it to nothing, and any method that restores it also amplifies the noise recorded at the same spatial frequencies. Every deblurring method is a compromise between sharpness and noise, set by a regularization parameter or by a number of iterations.

Microscopists use deconvolution routinely (Sage et al., 2017): in fluorescence microscopy the PSF can be measured and the noise follows photon statistics, and the Richardson–Lucy algorithm used there was derived for image restoration by Richardson (1972) and for statistical astronomy by Lucy (1974). In photography the blur is usually unknown and must be estimated with the image. Image processing places deblurring among the restoration problems in outline; removing noise alone is image denoising, and increasing resolution is super-resolution.

Problem Definition

The observation model writes the blurred image $y$ as the sharp image $x$ convolved with the PSF $h$, plus noise $n$:

$$y = h * x + n , \qquad \text{or, for images written as vectors,} \qquad y = Hx + n ,$$

where $H$ is the matrix of the discretized convolution (Sage et al., 2017). The PSF is the image of a single point of light, the system’s impulse response, as convolution explains. Szeliski (2022, Section 10.1.4) lists what shapes it in a camera: lens blur, anti-aliasing filters in front of the sensor, and the shape and extent of each pixel’s active area. In a microscope it is the diffraction pattern described in optical microscopy. The noise has a signal-dependent part, the Poisson-distributed number of photons arriving at a pixel, and a part usually modelled as additive Gaussian, which collects background, read-out and quantization noise (Sage et al., 2017); under Poisson statistics $n$ depends on the data, which the additive notation hides. Image noise gives the combined model.

The model assumes a linear, shift-invariant system, characterized by one PSF for the whole field. It also assumes that pixel values are proportional to light. Nah, Kim and Lee (2017) note that non-uniform deblurring becomes significantly more difficult when a nonlinear camera response is involved, and they invert an approximate gamma curve ($\gamma = 2.2$) before synthesizing blur; in our reading, the tone curve of a digital camera has to be undone in the same way before a photograph is deconvolved.

By the convolution theorem of the Fourier transform, the model is a product frequency by frequency,

$$\hat y(\omega) = \hat h(\omega)\, \hat x(\omega) + \hat n(\omega) ,$$

the form Levin et al. (2009) use, written here with hats for transforms. A blur is a low-pass filter: the transfer function $\hat h$ is close to 1 at low frequencies, small at high ones, and for some PSFs exactly zero at certain frequencies. Recovering $\hat x$ means dividing by $\hat h$ wherever it is small, which is where the noise $\hat n$ dominates. Lucy (1974) explained the consequence for direct inversion in general: the observed distribution, being blurred, is smoother than the true one, so solving for the true one greatly magnifies short-wavelength errors. In Szeliski’s (2022, Section 4.2) terms, deblurring is an inverse problem that is ill-posed, because many sharp images explain the same data, and ill-conditioned, because small changes in the data can produce large changes in the solution. Linear algebra works through the condition number of a small blur matrix and what quantization noise does to its exact inverse.

Input and Output

Non-blind deblurring takes the blurred image, the PSF and assumptions about the noise, and returns an estimate $\tilde x$ of the sharp image (Sage et al., 2017). Blind deblurring takes the image alone and must return both the sharp image and the PSF (Levin et al., 2009). The image can be a photograph, one channel of a multichannel image, or a 3D stack; microscopy deconvolves z-stacks with a 3D PSF (Sage et al., 2017).

The PSF comes from a measurement, a model or the image itself. An experimental microscope PSF is obtained by imaging a sparse field of very small beads (Sage et al. give 100 nm as an example diameter) and averaging several of them; it captures the microscope’s aberrations but is limited to the acquisition’s resolution. A theoretical PSF is computed from a model with the numerical aperture, wavelength and pixel size as parameters (Sage et al., 2017). For cameras, Szeliski (2022, Section 10.1.4) describes calibration targets of slanted edges, which give a one-dimensional profile of the PSF, and a pattern with edges at all orientations, from which the two-dimensional PSF is recovered by linear least squares against the known pattern; otherwise the PSF is inferred from knowledge of the motion, the defocus and the optics, or estimated from the image (Szeliski, 2022, Section 10.3).

Variants

• Non-blind and blind. With a known PSF the problem is linear in $x$. Blind deconvolution, the simultaneous inference of PSF and sharp image (Szeliski, 2022, Section 10.3), is nonlinear and has many more solutions; it is outlined below.
• Shift-invariant and spatially varying blur. The convolution model holds only when one PSF describes the whole image. In large micrographs the PSF varies across the field of view (Sage et al., 2017), and refractive-index mismatch between immersion medium and specimen makes it vary with depth, at which point the blur is no longer a convolution. Levin et al. (2009) found that even for a planar scene photographed with a hand-held camera, point traces in the four corners differed, which suggested that in-plane rotation is a significant component of hand shake. In dynamic scenes, Nah, Kim and Lee (2017) list multiple object motions, camera shake, depth variation and occlusions at motion boundaries as sources of non-uniform blur.
• Single and multiple images. Szeliski (2022, Section 10.3) lists ways of reducing blur with more data: combining a sharp but noisy image with a blurrier but cleaner one, taking many quick exposures, coded apertures that estimate depth and reduce blur together, inertial measurement units that help determine the blur kernel, and dual-pixel sensors that aid the deblurring of misfocused images. Coded exposure and coded apertures are described in computational photography.

Motion blur

Raskar, Agrawal and Tumblin (2006) describe the exposure time as a temporal box filter that smears a moving object across the image by convolution; in the frequency domain the signal is multiplied by a sinc function with zeros and significant attenuation at most other frequencies, so that deblurring becomes ill-posed. For uniform motion of $L$ pixels in a straight line, the PSF is a line segment of $L$ equal weights, and in our derivation its transfer function along the direction of motion has zeros at the non-zero multiples of $1/L$ cycles per pixel up to the Nyquist frequency; frequencies at those zeros are not attenuated but erased. Raskar et al.’s fluttered shutter, which computational photography describes, changes the box filter into a broadband filter and makes the deconvolution well-posed. Camera shake is rarely a straight line: in Levin et al.’s (2009) recordings, made with a partly loosened tripod (see Datasets and Benchmarks), an 85 mm lens and 0.3 s exposures, the support of the measured kernels varied from 10 to 25 pixels.

Defocus blur

In geometric optics a point outside the focal plane is imaged as a patch whose size, the circle of confusion, depends on how far the image plane is from focus relative to the aperture diameter (Szeliski, 2022, Section 2.2.3). Szeliski notes that with a non-circular aperture, such as a hexagonal diaphragm, the shape can sometimes be seen in the actual blur function. A defocus PSF is therefore approximately a scaled image of the aperture, a uniform disc for a circular one; in our derivation, as for the line segment, its transfer function has zeros, here on rings, so defocus also erases frequencies. Because the circle of confusion grows with distance from the focal plane, defocus in a scene with depth is spatially varying. In a microscope, defocus is the axial extent of the 3D PSF, which optical microscopy describes.

Challenges

Besides noise amplification and the model violations listed under Variants, two difficulties recur.

• An inexact PSF. Nah, Kim and Lee (2017) point out that kernel estimation is sensitive to noise and saturation and that incorrectly estimated kernels produce artifacts; a measured PSF is noisy and limited to the acquisition’s resolution, and a modelled one omits aberrations (Sage et al., 2017).
• Choosing the regularization. The parameter that balances sharpness against noise, or the number of iterations, decides the result, and without ground truth it is set by a rule or by eye (Sage et al., 2017).

Boundaries and ringing

Fourier-based implementations compute circular convolutions, so the data are implicitly periodic: structures near one border spill over to the opposite border after processing, as ghosts (Sage et al., 2017). The data are also treated as band-limited, so sharp edges incur local overshoots and undershoots, the ringing artifacts that Sage et al. relate to the Gibbs phenomenon. A non-negativity constraint can cancel only the undershoots that would otherwise be negative. Sage et al. implement two countermeasures: apodization, which multiplies the data by a window that falls gradually to zero near the borders and cancels the ghost but also reduces the intensity of the data; and zero-padding, which keeps the data intact and extends the image with zeros, at a higher computational cost.

A real photograph does not wrap around: its border pixels contain light from outside the frame, which no deconvolution of the frame alone can account for. Benchmarks that synthesize blur can sidestep this: Zhang et al. (2017) assume circular boundary conditions when blurring their test images, so that their quadratic step can be computed with the FFT. scikit-image’s wiener also relies on Fourier diagonalization and therefore on a circulant blur matrix, while its richardson_lucy uses SciPy’s convolve in same mode, a linear convolution that treats pixels outside the image as zero (scikit-image 0.25.2 source). The second example below compares periodic processing, apodization and mirror padding on a blur that does not wrap.

Approaches

The methods below are classical except the last. Sage et al. (2017) cite Wiener’s 1949 monograph for the Wiener filter, Tikhonov’s 1963 paper for regularization and Landweber’s 1951 paper for his iteration; Richardson (1972) and Lucy (1974) published the iteration that bears their names. Deconvolution of micrographs began in the 1980s, as fluorescence microscopy recounts. Before 2009, blind deconvolution of photographs had progressed by modelling the statistics of natural images, as Levin et al. (2009) review, and learned methods followed in the 2010s (Szeliski, 2022, Section 10.3).

Inverse filtering

The simplest approach minimizes the least-squares cost $\lVert y - Hx \rVert^2$, which is the maximum-likelihood estimate under Gaussian noise. The solution is a division in the Fourier domain, $\hat x = \hat y / \max(\hat h, \epsilon)$ with a small constant $\epsilon$ against division by zero; it is parameter-free and fast but amplifies noise into spurious high-frequency oscillations (Sage et al., 2017). On their synthetic test volume, Sage et al. describe its results as non-exploitable.

Tikhonov regularization and the Wiener filter

Tikhonov regularization adds a penalty on large values of the solution to the least-squares cost, and a regularized inverse filter penalizes the energy of a derivative instead (Sage et al., 2017):

$$C(x) = \lVert y - Hx \rVert^2 + \lambda \lVert Lx \rVert^2 , \qquad \tilde x = \left( H^{\mathsf{T}} H + \lambda L^{\mathsf{T}} L \right)^{-1} H^{\mathsf{T}} y ,$$

with $L = I$ for Tikhonov regularization and a discrete differential operator, the Laplacian in DeconvolutionLab2, for the regularized inverse filter. The parameter $\lambda$ balances data fidelity against regularization; $\lambda = 0$ returns the inverse filter, and larger values give smoother images. Sage et al. also read the penalty as prior information about the signal, which makes the estimate a maximum a posteriori one, and the least squares article treats the same penalty as generalized Tikhonov regularization. With circular boundaries both $H$ and $L$ are diagonalized by the discrete Fourier transform, and the solution becomes the filter

$$\hat{\tilde x}(\omega) = \frac{\hat h^{*}(\omega)}{|\hat h(\omega)|^2 + \lambda\, |\hat l(\omega)|^2}\, \hat y(\omega) ,$$

which is how scikit-image’s wiener computes it, with a Laplacian $\hat l$ by default (scikit-image 0.25.2 documentation). Where $|\hat h|$ is large the filter divides by $\hat h$; where $|\hat h|^2$ falls below $\lambda |\hat l|^2$ it attenuates instead of amplifying.

The Wiener filter is the special case in which the weights come from statistics: Sage et al. note that when $L^{\mathsf{T}} L$ has a whitening effect on $x$ and $\lambda$ is set by the noise variance, the regularized inverse filter amounts to Wiener filtering. Szeliski (2022, Section 3.4.1) describes it as the optimum restoration filter that follows from assuming that an image is a sample of a correlated Gaussian random field, combined with a statistical model of the measurement, and explains that he dropped its derivation from his second edition because the filter is almost never used in practice any more, having been replaced by better-performing nonlinear filters; its formula in terms of the image and noise power spectra is in image processing. In Sage et al.’s experiments, Tikhonov regularization and the regularized inverse filter gave decent results, but ringing remained; on an image of a bead, too much regularization blurred the bead and overestimated its size, and too little produced ringing.

Quadratic penalties favour smooth images. Sage et al. also implement penalties on the $\ell_1$ norm of wavelet coefficients, solved by fast iterative soft-thresholding, which better preserve details and discontinuities, and total variation, the $\ell_1$ norm of the gradient, which yields piecewise-constant results; the same penalties are used for image denoising.

Landweber iteration

The Landweber algorithm minimizes the same least-squares cost as the inverse filter, but by iterative gradient descent instead of direct inversion, and DeconvolutionLab2 uses the iteration to impose non-negativity at every step (Sage et al., 2017):

$$x^{(k+1)} = \max\!\left( x^{(k)} + c\, H^{\mathsf{T}} \big( y - H x^{(k)} \big),\; 0 \right) ,$$

with a step size $c$ and the maximum taken pixel by pixel. At convergence the solution overfits the noise, but stopping early gives a satisfactory trade-off between deconvolution and noise amplification, so the number of iterations acts as a pseudo-regularization parameter; Sage et al. note that this happens for all maximum-likelihood algorithms. Iterating in the same way on the regularized cost gives the Tikhonov–Miller algorithm.

In our derivation, without the projection and starting from $x^{(0)} = 0$, the $k$-th iterate is the filter

$$\hat x^{(k)}(\omega) = \frac{1 - \left( 1 - c\, |\hat h(\omega)|^2 \right)^k}{\hat h(\omega)}\, \hat y(\omega) ,$$

which equals the inverse filter at frequencies where $c|\hat h|^2$ is close to 1 and, where it is small, approximately $k c\, \hat h^{*} \hat y$, a damped response that grows with every iteration. Each iteration therefore releases a little more of the poorly transmitted frequencies, and $1/(ck)$ plays roughly the role of $\lambda$ in Tikhonov regularization with $L = I$.

Richardson–Lucy

The Richardson–Lucy (RL) algorithm, named after Richardson (1972) and Lucy (1974), is the maximum-likelihood method for Poisson noise. Sage et al. (2017) write its cost as

$$C(x) = \mathbf{1}^{\mathsf{T}} Hx - y^{\mathsf{T}} \log(Hx)$$

and its iteration as a multiplicative update

$$x^{(k+1)} = x^{(k)} \odot H^{\mathsf{T}} \left( \frac{y}{H x^{(k)}} \right),$$

where the product $\odot$, the division and the logarithm act element by element and $H^{\mathsf{T}}$ is the adjoint of the blur, a convolution with the mirrored PSF. In our derivation, setting the gradient of $C$ to zero gives the fixed point $x = x \odot H^{\mathsf{T}}(y / Hx) / H^{\mathsf{T}}\mathbf{1}$; the update as written assumes a PSF that sums to one, so that $H^{\mathsf{T}}\mathbf{1} = \mathbf{1}$ away from the borders. Because the update is multiplicative, a non-negative start stays non-negative (Sage et al., 2017). Lucy (1974), deriving the method from Bayes’ theorem for frequency distributions in astronomy, showed that it conserves normalization and non-negativity and increases the likelihood of the observed sample at every iteration. Since the cost is the Poisson negative log-likelihood up to terms that do not depend on $x$, minimizing it is maximum-likelihood estimation under the noise model of photon-counting images.

Lucy also described how the iteration proceeds. Deviations between the data and the re-blurred estimate on length scales large compared with the PSF are removed in essentially one iteration, while short-wavelength deviations are largely averaged out by the blur and produce only small corrections; he regarded this as desirable, since short-wavelength deviations are the most likely to be noise, and expected that after a few iterations further ones would mainly fit the statistical fluctuations. Raskar et al. (2006) describe RL as a nonlinear, ratio-based method that produces non-negative grey levels. Levin et al. (2009), evaluating blind deconvolution code that used RL as its final non-blind step, called it a source of ringing artifacts and replaced it with a sparse deconvolution.

Adding a total-variation penalty to the cost (RL-TV) counters the noise amplification and favours piecewise-constant results that preserve edges (Sage et al., 2017). On their synthetic microtubule data, RL-TV gave results similar to plain RL, which Sage et al. attribute to thin filaments being hard to recover with a regularizer that promotes piecewise-constant surfaces.

Regularization and stopping

The catch of the iterative maximum-likelihood methods is the same as for inverse filtering. A maximum-likelihood solution fits the noise as well as the signal, so RL and Landweber amplify noise as they converge, and the number of iterations acts as the regularization parameter: the algorithm is stopped before noise dominates, at a count that Sage et al. (2017) describe as set heuristically. They give two rules for the iteration cap: methods known to reach the desired solution at convergence need it large enough, while methods that amplify noise during convergence need it small enough, and then its choice has a crucial impact on the result.

Lucy (1974) documented the effect in a test on a one-dimensional distribution sampled 150 times. The likelihood increased and the misfit to the data decreased monotonically with the iterations, but the misfit to the known true distribution was smallest at three iterations and grew afterwards; its minimum coincided with the re-blurred estimate fitting the data more closely than the true distribution itself would, so the later iterations were merely fitting the statistical fluctuations. Repeated on forty independent samples, stopping at three iterations gave good fits in general. Because the truth is unknown in practice, Lucy proposed stopping on the goodness of fit to the data: with a $\chi^2$ test, as soon as the probability that random sampling would give a larger $\chi^2$ is not small (not below 0.05, for example). In our reading, this is a discrepancy rule: stop when the residual is as small as the noise alone would leave it. The second example below applies such a rule to Landweber iteration and to the regularized filter.

Blind deconvolution

Blind deconvolution attempts to estimate the object and the PSF jointly from the data alone. It is strongly ill-posed and nonlinear, and among other degeneracies it must resolve the trade of scale, since $(a h) * (x / a) = h * x$ for any non-zero $a$ (Sage et al., 2017). Sage et al. report that blind algorithms for fluorescence microscopy alternate between deconvolution and PSF estimation, some with an optically motivated parameterization of the PSF.

Levin et al. (2009) analysed blind deconvolution of photographs, where one undesirable solution that satisfies the model perfectly is the no-blur explanation: a delta kernel and $x = y$. They showed that the simplest approach, maximizing the posterior probability of the image and kernel together under a sparse prior on image derivatives, mostly favours that explanation, because under such a prior a blurred image is usually more probable than a sharp one. Estimating the kernel alone, by marginalizing over the image, is well constrained because the kernel has far fewer unknowns than the image, and with that estimator even a weak Gaussian prior recovered the blur; they concluded that the choice of estimator matters more than the choice of prior. In their evaluation on recorded camera shake, the variational Bayes method of Fergus et al. (2006) outperformed the alternatives they tested.

Learned deblurring

Szeliski (2022, Section 10.3) notes a large number of learning-based deblurring algorithms introduced over the past decade. Two designs illustrate the range. Nah, Kim and Lee (2017) trained a multi-scale convolutional network to map blurred to sharp images end to end, without estimating a kernel, arguing that kernel convolution cannot model occluded regions or depth variation and that finding a spatially varying kernel for every pixel is costly; their training pairs average successive frames of high-speed video, as described under Datasets and Benchmarks. Zhang et al. (2017) keep the observation model and learn only the prior. Half-quadratic splitting alternates between a quadratic step, $x_{k+1} = (H^{\mathsf{T}} H + \mu I)^{-1} (H^{\mathsf{T}} y + \mu z_k)$, which has the Tikhonov form above, and a step that, by their Bayesian reading, is Gaussian denoising of $x_{k+1}$ at noise level $\sqrt{\lambda / \mu}$; they plug in a set of 25 convolutional denoisers trained for noise levels from 0 to 50, decrease the denoiser’s noise level over 30 iterations, and so need far fewer networks than training one model for each degradation.

In microscopy, learned restoration recovers images from far fewer photons, but in our reading it infers structure from training data rather than measuring it, and restoration methods applied improperly can give misleading images, as bioimage analysis notes.

Practical Example: Richardson–Lucy at Two Photon Levels

The first example simulates widefield fluorescence images of a synthetic specimen whose ground truth is known and deconvolves a bright and a dim frame with the known PSF. The specimen, on a 128 × 128 grid, contains a uniform disc standing in for a nucleus, three thin lines standing in for filaments, two isolated point sources, and a pair of point sources 4 pixels apart; it emits a stated number of photons per pixel per frame, above a constant background of 2 photons per pixel. The PSF is a Gaussian with $\sigma = 2$ pixels, the in-focus approximation that optical microscopy mentions; the pair’s separation, $2\sigma$, is where the sum of two such Gaussians no longer has a dip between them. The specimen bleaches mono-exponentially with a time constant of 20 frames over 60 frames, so the last frame carries about 5% of the first frame’s photons, and each frame is a Poisson sample of its expected counts, with no read noise. The simulation is that of the example in fluorescence microscopy, which measures the bleaching itself; its set-up lines are repeated here so that the block runs on its own.

The code runs scikit-image’s richardson_lucy with the known PSF for increasing numbers of iterations on three inputs: the noise-free first frame, the noisy first frame and the noisy last frame. The implementation in 0.25.2 starts from a constant image of 0.5, uses zero-padded convolution with the PSF and its mirror, and by default clips the result to $[-1, 1]$, so raw counts are passed with clip=False (scikit-image 0.25.2 source and API reference). The error is the root-mean-square difference from the unblurred specimen plus background, relative to the norm of that reference, away from the borders; the pair’s dip ratio is the value at its midpoint divided by the lower of the values at the two point positions, so a value below 1 means the pair is seen as two. The last lines repeat the search for the best iteration count with random seeds 0 to 4. It uses NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2; all parameters were fixed before the first run, and no change was made after it.

import numpy as np
from scipy.signal import fftconvolve
from skimage.draw import disk, line
from skimage.restoration import richardson_lucy

# Ground truth: emission in photons per pixel per frame at t = 0
obj = np.zeros((128, 128))
obj[disk((40, 40), 12)] = 100.0                       # uniform "nucleus"
for r0, c0, r1, c1 in [(90, 10, 120, 60), (70, 70, 70, 120), (10, 80, 50, 120)]:
    obj[line(r0, c0, r1, c1)] = 150.0                 # thin "filaments"
for r, c in [(100, 90), (40, 95), (100, 100), (100, 104)]:
    obj[r, c] = 5000.0                                # point sources; last two are a pair 4 px apart
background = 2.0                                      # photons/px, does not bleach

y, x = np.mgrid[-12:13, -12:13]
psf = np.exp(-(x**2 + y**2) / (2 * 2.0**2))
psf /= psf.sum()                                      # Gaussian PSF, sigma = 2 px

blurred = fftconvolve(obj, psf, mode="same")
tau, T = 20.0, 60                                     # bleaching time constant (frames), frame count

def expected(t):
    return blurred * np.exp(-t / tau) + background

def series(seed):
    rng = np.random.default_rng(seed)
    return np.stack([rng.poisson(expected(t)) for t in range(T)]).astype(float)

frames = series(0)

# Richardson-Lucy with the known PSF
inner = (slice(12, 116), slice(12, 116))
iters = [0, 5, 10, 20, 50, 100, 200, 500]

def rl_curve(img, truth):
    errs, dips = [], []
    for n in iters:
        est = img if n == 0 else richardson_lucy(img, psf, num_iter=n, clip=False)
        d = est[inner] - truth[inner]
        errs.append(np.linalg.norm(d) / np.linalg.norm(truth[inner]))
        dips.append(est[100, 102] / min(est[100, 100], est[100, 104]))
    return errs, dips

cases = [("noise-free t=0", expected(0), obj + background),
         ("noisy t=0", frames[0], obj + background),
         ("noisy t=59", frames[59], obj * np.exp(-59 / tau) + background)]
print("iterations:      " + " ".join(f"{n:>6d}" for n in iters))
for name, img, truth in cases:
    errs, dips = rl_curve(img, truth)
    print(f"{name:<15}  rel. RMSE " + " ".join(f"{e:6.3f}" for e in errs))
    print(f"{'':<15}  pair dip  " + " ".join(f"{d:6.3f}" for d in dips))

print("best iteration count and its rel. RMSE over seeds 0-4:")
for t in (0, 59):
    truth = obj * np.exp(-t / tau) + background
    best = []
    for seed in range(5):
        rng = np.random.default_rng(seed)
        img = np.stack([rng.poisson(expected(s)) for s in range(T)])[t].astype(float)
        errs, _ = rl_curve(img, truth)
        k = int(np.argmin(errs))
        best.append((iters[k], errs[k]))
    n_best = [b[0] for b in best]
    e_best = [b[1] for b in best]
    print(f"  t={t:2d}: iterations {min(n_best)}-{max(n_best)}, rel. RMSE {min(e_best):.3f}-{max(e_best):.3f}")

iterations:           0      5     10     20     50    100    200    500
noise-free t=0   rel. RMSE  0.947  0.929  0.904  0.869  0.800  0.726  0.630  0.479
                 pair dip   1.068  1.315  1.247  1.058  0.716  0.487  0.304  0.135
noisy t=0        rel. RMSE  0.945  0.929  0.905  0.871  0.804  0.733  0.639  0.480
                 pair dip   1.211  1.386  1.330  1.122  0.782  0.547  0.346  0.147
noisy t=59       rel. RMSE  0.913  0.858  0.850  0.844  0.848  0.873  0.936  1.102
                 pair dip   2.111  1.568  1.698  1.681  1.417  1.146  0.870  0.499
best iteration count and its rel. RMSE over seeds 0-4:
  t= 0: iterations 500-500, rel. RMSE 0.469-0.513
  t=59: iterations 20-20, rel. RMSE 0.844-0.855

At the start of the series, with about 200 photons at the peak of each point source, the noisy curves follow the noise-free ones closely. The error keeps falling up to the largest tested count, 500 iterations, for all five seeds, and the pair, a single peak in the input, shows a dip from 50 iterations on. In the first few iterations the dip ratio rises; in our reading, the algorithm first sharpens the blended light into one peak before it separates the two sources. Even after 500 iterations the error is still about one half, in our reading because single-pixel points and lines carry spatial frequencies that the Gaussian blur attenuates almost completely, which RL recovers only slowly, as Lucy’s description of short-wavelength corrections leads one to expect.

At the end of the series, with about 5% of the original photons, the error is lowest at 20 iterations for every seed; for seed 0 it then grows, exceeding that of the blurred input by 200 iterations, as the iterations fit the Poisson noise. The pair splits only from 200 iterations on, where the image as a whole is already worse than the input. At low counts, stopping early limits noise but leaves close structures merged, and the best number of iterations differs between frames of the same experiment: the stopping problem that Sage et al. (2017) describe, made concrete by bleaching. The set-up is deliberately favourable, with a known, shift-invariant Gaussian PSF, no read noise and no out-of-focus light; we would expect real data, whose PSF is only measured or modelled, to fare worse.

Practical Example: Motion and Defocus Blur

The second example blurs scikit-image’s bundled camera image with two PSFs, a 21-pixel horizontal motion and a uniform disc of radius 6 pixels standing in for defocus, and compares a regularized filter under three boundary treatments with Landweber iteration. The blur is applied to the whole 512 × 512 image and only the central 384 × 384 pixels are kept, so that, as in a photograph, the borders contain light from outside the field and the data do not wrap around. Gaussian noise with a standard deviation of 0.01, on intensities from 0 to 1, is added, one realization with seed 0. The filter is scikit-image’s wiener, the regularized inverse filter above with its default Laplacian penalty, for five values of its balance parameter $\lambda$. It is applied to the observation as it is (periodic), after multiplying it by a Tukey window that tapers the outer 12.5% of each side to zero (apodized, one of Sage et al.’s countermeasures), and after mirroring 32 pixels of the image outward on every side and cropping the result back (mirror padding, our choice). The PSNR is measured against the sharp crop; the residual is the root-mean-square difference between the observation and the re-blurred estimate, to be compared with the noise level of 0.01. For the best balance, the PSNR is also given separately for the 32-pixel band along the borders and for the interior. Landweber iteration with non-negativity runs on the mirror-padded image with step size 1, starting from the observation. The example uses NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2. One setting was changed after the first run: with the planned balance values, $10^{-4}$ to $10^{-1}$, the mirror-padded PSNR was still rising at the largest value for both blurs, so the grid did not contain the best balance; a second run with $10^{-3}$ to $10$, and a diagnostic run, still put the best balance for motion at the end of the grid, so the grid was moved to $10^{-2}$–$10^{2}$. Nothing else was changed, and the rows for $10^{-2}$ and $10^{-1}$ and the Landweber rows are the same as in the first run.

import numpy as np
from scipy.signal import fftconvolve
from scipy.signal.windows import tukey
from skimage import data
from skimage.draw import disk
from skimage.metrics import peak_signal_noise_ratio
from skimage.restoration import wiener

sharp = data.camera() / 255.0                         # bundled 512 x 512 image
field = (slice(64, 448), slice(64, 448))              # 384 x 384 field of view
truth = sharp[field]
sigma, P = 0.01, 32                                   # noise level, padding width

motion = np.zeros((21, 21)); motion[10, :] = 1.0      # 21-px horizontal motion
defocus = np.zeros((13, 13)); defocus[disk((6, 6), 6)] = 1.0   # disc, radius 6 px
psfs = {"motion": motion / motion.sum(), "defocus": defocus / defocus.sum()}
rng = np.random.default_rng(0)

def psnr(x, ref=truth):
    return peak_signal_noise_ratio(ref, x, data_range=1.0)

def otf(psf, shape):                                  # PSF centred at the origin
    k = np.zeros(shape)
    k[:psf.shape[0], :psf.shape[1]] = psf
    k = np.roll(k, (-(psf.shape[0] // 2), -(psf.shape[1] // 2)), axis=(0, 1))
    return np.fft.rfft2(k)

def mirror(y):
    return np.pad(y, P, mode="symmetric")

def residual(y, x_pad, H):                            # RMS of y - h*x over the field
    hx = np.fft.irfft2(H * np.fft.rfft2(x_pad), s=x_pad.shape)[P:-P, P:-P]
    return np.sqrt(np.mean((y - hx) ** 2))

window = np.outer(tukey(384, 0.25), tukey(384, 0.25))
band = np.ones((384, 384), bool); band[P:-P, P:-P] = False

for name, psf in psfs.items():
    y = fftconvolve(sharp, psf, mode="same")[field]  # light from outside the field
    y = y + sigma * rng.standard_normal(y.shape)
    H = otf(psf, (384 + 2 * P, 384 + 2 * P))
    print(f"{name}: blurred input PSNR {psnr(y):.2f} dB")
    print("  balance   periodic  apodized  mirror   residual RMS (mirror)")
    results = {}
    for bal in [1e-2, 1e-1, 1, 10, 100]:
        x_per = wiener(y, psf, bal, clip=False)
        x_apo = wiener(y * window, psf, bal, clip=False)
        x_mir_pad = wiener(mirror(y), psf, bal, clip=False)
        x_mir = x_mir_pad[P:-P, P:-P]
        results[bal] = (x_per, x_apo, x_mir)
        print(f"  {bal:7.0e}  {psnr(x_per):8.2f}  {psnr(x_apo):8.2f}  {psnr(x_mir):6.2f}"
              f"   {residual(y, x_mir_pad, H):.4f}")
    best = max(results, key=lambda b: psnr(results[b][2]))
    print(f"  balance {best:.0e}, PSNR border band / interior:")
    for label, x in zip(["periodic", "apodized", "mirror"], results[best]):
        print(f"    {label:<9} {psnr(x[band], truth[band]):6.2f} / {psnr(x[~band], truth[~band]):6.2f}")

    # Landweber with non-negativity, mirror padding, step 1
    yp = mirror(y)
    x = yp.copy()
    HtY = np.conj(H) * np.fft.rfft2(yp)
    HtH = np.abs(H) ** 2
    checkpoints = [1, 3, 10, 30, 100, 300, 1000, 3000, 10000]
    rows = []
    for k in range(1, checkpoints[-1] + 1):
        X = np.fft.rfft2(x)
        x = np.maximum(np.fft.irfft2(X + (HtY - HtH * X), s=yp.shape), 0.0)
        if k in checkpoints:
            rows.append((k, psnr(x[P:-P, P:-P]), residual(y, x, H)))
    print("  Landweber iterations " + " ".join(f"{r[0]:>6d}" for r in rows))
    print("  PSNR (dB)            " + " ".join(f"{r[1]:6.2f}" for r in rows))
    print("  residual RMS         " + " ".join(f"{r[2]:6.4f}" for r in rows))

motion: blurred input PSNR 21.26 dB
  balance   periodic  apodized  mirror   residual RMS (mirror)
    1e-02      8.46     11.50    9.16   0.0097
    1e-01     13.75     11.46   15.35   0.0120
    1e+00     18.50     11.39   21.41   0.0158
    1e+01     20.75     11.31   21.94   0.0214
    1e+02     20.47     11.25   21.11   0.0304
  balance 1e+01, PSNR border band / interior:
    periodic   19.86 /  21.21
    apodized    6.52 /  20.82
    mirror     23.89 /  21.30
  Landweber iterations      1      3     10     30    100    300   1000   3000  10000
  PSNR (dB)             21.61  22.02  22.90  23.92  22.05  16.97  11.80   8.23   5.04
  residual RMS         0.0171 0.0147 0.0124 0.0104 0.0096 0.0093 0.0089 0.0088 0.0087
defocus: blurred input PSNR 22.31 dB
  balance   periodic  apodized  mirror   residual RMS (mirror)
    1e-02     15.84     11.56   21.56   0.0100
    1e-01     20.87     11.51   24.63   0.0107
    1e+00     22.55     11.45   23.53   0.0120
    1e+01     22.11     11.39   22.63   0.0163
    1e+02     21.00     11.32   21.54   0.0278
  balance 1e-01, PSNR border band / interior:
    periodic   17.50 /  24.03
    apodized    6.57 /  23.18
    mirror     26.07 /  24.12
  Landweber iterations      1      3     10     30    100    300   1000   3000  10000
  PSNR (dB)             22.68  23.04  23.69  24.64  25.01  21.21  14.63  10.22   7.53
  residual RMS         0.0155 0.0132 0.0115 0.0104 0.0093 0.0088 0.0082 0.0076 0.0073

Too little regularization is catastrophic: at $\lambda = 10^{-2}$ the motion-blurred image comes out at 9 dB, far worse than the blurred input, as the filter nearly divides by a transfer function that is close to zero along whole lines of frequencies. The best balance improves on the input by 2.3 dB for defocus but by only 0.7 dB for motion, where the residual at the best balance is twice the noise level. In our reading, that residual shows the filter smoothing detail the blur let through: a Laplacian penalty is weakest at low frequencies, where the motion blur has its first zeros (at 1/21 cycles per pixel), so a balance large enough to control the noise near those zeros oversmooths the rest of the image.

The boundary rows show where periodic processing fails. At the best balance its interior is almost as good as with mirror padding, but its border band is 4 to 9 dB worse, the ghosts and ringing of a periodic extension that joins unrelated image content. Apodization removes them but darkens the border band to about 6.5 dB, which caps its full-field PSNR near 11 dB whatever the balance, the reduction of intensity Sage et al. (2017) mention. Mirror padding, which makes the periodic extension continuous, gives the best border band in both cases.

Landweber iteration shows the behaviour Lucy described. The PSNR rises to 23.9 dB for motion at 30 iterations and to 25.0 dB for defocus at 100, better than the best filter in both cases, then falls far below the input while the residual keeps decreasing: the later iterations fit the noise. In our reading, Landweber does better on the motion blur because its damping depends only on $|\hat h|$, so it holds back every poorly transmitted frequency, low or high. A discrepancy rule that stops at the first checkpoint with a residual below the noise level picks 100 iterations in both cases, the best checkpoint for defocus and one checkpoint late for motion, at 22.1 dB. Applied to the filter, scanning from strong to weak regularization, the same rule would not stop before the weakest balance in the grid, $\lambda = 10^{-2}$, where the residual reaches the noise level but the image is worse than the blurred input for both blurs; in our reading the residual cannot see errors at frequencies the blur nearly removes, since they barely change the re-blurred estimate.

Datasets and Benchmarks

• Levin et al. (2009). To compare blind deconvolution algorithms on real camera shake with ground truth, Levin et al. photographed a planar scene with the camera on a tripod whose rotation about the optical axis was locked while the other two axes were loose, so that the blur was approximately shift-invariant. They calibrated the kernels of 8 blurred images from a sharp reference and cropped four 255 × 255 windows from each, giving 32 test images.
• GOPRO (Nah, Kim & Lee, 2017). Blurred images synthesized by averaging 7 to 13 successive frames of 240 frames-per-second GoPro video after inverting an approximate camera response, with the middle frame as the sharp reference: 3,214 pairs at 1280 × 720, of which 2,103 were used for training and 1,111 for testing. The text also gives averaging 15 frames as an example, which lies outside the stated range. The blur is spatially varying, as in their example of a moving person in front of a static background.
• Köhler et al. As Nah, Kim and Lee (2017) describe it, this dataset has 4 sharp images and 12 blurred versions of each, created by replaying recorded camera motion with six degrees of freedom; it comes with its own evaluation code.
• Microscopy. Sage et al. (2017) evaluated DeconvolutionLab2’s algorithms on a synthetic microtubule volume blurred with a theoretical PSF and corrupted by mixed Gaussian and Poisson noise, on a bead of known diameter and on a C. elegans embryo. They report a signal-to-noise ratio after RL of 3.66 dB on the synthetic volume, the highest of the classical algorithms in DeconvolutionLab2, against 7.27 dB for the learned “Learn 2D, Apply 3D” method that won the 2014 3D deconvolution microscopy challenge (Table 2); the parameters of each algorithm were set by visual assessment, and the authors caution that the ranking depends on the data.

Evaluation Metrics

On synthetic data the result is compared with the known sharp image by PSNR, SSIM or a relative error, as in the examples above; perceptual measures and no-reference quality assessment are in image quality assessment. Levin et al. (2009) observed that wider kernels give larger errors even with the true kernel, so they score a blind method by the ratio of the error after deconvolution with its estimated kernel to the error with the true kernel; empirically, ratios above 2 were already visually implausible. On GOPRO, Nah, Kim and Lee (2017, Table 2) report a PSNR of 29.23 dB and an SSIM of 0.9162 for their network with two scales, against 24.64 dB and 23.64 dB for the methods of Sun et al. and Kim and Lee. On the Köhler dataset, whose evaluation code reports multi-scale SSIM, their Table 3 gives the highest PSNR (26.48 dB) to the three-scale network and the highest multi-scale SSIM (0.8116) to the two-scale one, while the text states the reverse. When blurred and sharp images are not aligned, as in Lai et al.’s dataset, Nah et al. note that PSNR and SSIM correlate less with perceived quality. In microscopy, Sage et al. (2017) also report the full width at half maximum of a deconvolved bead of known diameter.

Applications

• Microscopy. Widefield, confocal and light-sheet stacks are deconvolved routinely, as fluorescence microscopy describes. Light-sheet data need PSFs of their own: Becker et al. (2019) deconvolved a light-sheet stack of 667 slices recorded with a 2.5× objective (NA 0.12) in DeconvolutionLab with two modelled confocal PSFs and found the results poor or even worse than the original stack. With RL and a PSF modelled for light-sheet image formation, which matched measured bead PSFs for 10× and 20× objectives, they needed no regularization for the cleared-tissue stacks they show, which they attribute to the high signal-to-noise ratio of light-sheet recordings; after Gaussian noise was added to a stack, RL without damping amplified it and 5% damping suppressed it.
• Astronomy. Lucy (1974) applied his iteration to rectifying distributions of projected stellar rotational velocities and reported calculations for correcting radio-astronomical observations for beam smoothing.
• Photography. Motion blur from camera shake and moving objects, and defocus, are the main targets, with the blur usually unknown; coded exposure and coded apertures change the capture so that the blur can be inverted, as computational photography describes. Szeliski (2022, Section 10.1.4) notes that the colour-dependent blur of chromatic aberration can also be removed with deblurring techniques.

Open Problems

• Spatially varying blur. Real camera shake includes rotation, scenes have depth, and microscope PSFs vary across the field and with depth, while most algorithms assume one PSF (Levin et al., 2009; Sage et al., 2017; Nah, Kim & Lee, 2017).
• Choosing the regularization without ground truth. Iteration counts and regularization weights are still set heuristically or by eye (Sage et al., 2017), and the first example shows the best count changing within one experiment as the signal fades.
• Learned methods outside their training data. A network trained on synthesized blur learns that blur’s statistics; in our reading, structure it restores in other data is inferred rather than measured, which matters most where images are measurements, as in microscopy.