Richardson-Lucy Deconvolution: Iteration Count, Regularisation, and Reading the Artefacts

The setting that decides whether a Richardson–Lucy (RL) result is publishable is the iteration count, and the most common mistake is treating it as a quality dial to be turned up until the image looks crisp. It is not a quality dial. It is the regularisation parameter — scikit-image says exactly that in one line of its API reference, describing num_iter as the parameter that “plays the role of regularisation” — and the algorithm it controls has no natural stopping point to converge to.

That is not a folk belief about a fussy algorithm. Liu, Panezai, Wang and Stallinga showed in Nature Communications in 2025 that RL is necessarily ill-convergent: an analytical Cramér–Rao Lower Bound for the estimate diverges for spatial frequencies approaching the diffraction limit from below, so a regular optimum of the likelihood does not exist. Their conclusion, by reductio ad absurdum, is that “RL-deconvolution cannot converge.” Running more iterations does not walk you towards a right answer that is waiting at infinity. It walks you towards noise.

Everything else on this page follows from that: pick a stopping rule and defend it, get the point spread function (PSF) right because the algorithm cannot fix a wrong one, and learn to recognise the four or five artefacts RL reliably produces so you do not report them as biology.

What one RL iteration actually does

RL is a multiplicative, iterative maximum-likelihood estimator for an image formed by convolution with a known PSF and corrupted by Poisson (photon-counting) noise. Each iteration forward-convolves the current object estimate with the PSF, divides the measured image by that prediction to form a ratio image, correlates the ratio back through the PSF, and multiplies the current estimate by the result. Where the prediction matches the data the ratio is one and nothing changes; where the prediction is too dim the estimate is scaled up.

It was published twice independently — by William Richardson in JOSA in 1972 as a Bayesian iterative restoration method, and by Leon Lucy in The Astronomical Journal in 1974 as a technique for rectifying observed distributions. Shepp and Vardi later showed in IEEE Transactions on Medical Imaging (1982) that the same update is the expectation–maximisation algorithm for maximum-likelihood emission tomography, which is why the identical iteration shows up in PET reconstruction, astronomy and fluorescence microscopy under different names.

Two structural properties matter more than the derivation, and both come from the multiplicative form. Prato and colleagues state them plainly in Astronomy & Astrophysics (2012): the method “preserves the number of counts and the non-negativity of the original object.”

• Non-negativity is built in, not imposed. A strictly positive estimate multiplied by a positive ratio stays positive. Linear inverse or Wiener-type filters have no such guarantee and routinely return negative intensities in dark regions, which then have to be clipped.
• Total flux is conserved. RL redistributes photons; it does not create or destroy them. This is what makes integrated-intensity measurements survive deconvolution in a way that peak-intensity measurements do not.

Why more iterations is not more resolution

The behaviour practitioners describe as “it gets sharper, then it gets grainy” is semiconvergence: the estimate approaches the true object for a while, then diverges from it while the likelihood keeps improving, because the extra likelihood is being bought by fitting noise.

Liu et al. characterised what that noise looks like. In their microscopy experiments they ran 300 iterations and stored intermediate estimates at iterations 1, 20, 100 and 300. The amplified noise has “a spatial frequency content that strongly peaks for the spatial frequencies just below the cutoff,” and that peak migrates monotonically closer to the diffraction cutoff as iteration proceeds. In practice you see it as a fine, high-frequency granularity at roughly the resolution scale of the system — which is exactly the scale at which you were hoping to find real structure, and exactly why it is so easy to mistake for signal.

The same paper deals with the super-resolution question directly. RL does extrapolate some signal beyond the diffraction limit, but the spectral signal-to-noise ratio gain for those frequencies “drops to zero when the number of iterations goes to infinity.” Any apparent super-resolution is a transient of early iterations, not a limit the algorithm reaches. If you need a genuine, physically grounded resolution gain rather than a contrast improvement, that is a hardware or sample-preparation decision — see our guides to the optical microscopy resolution limit and to expansion microscopy.

A 2024 preprint by Hendrix and colleagues frames the practical consequence well: RL “is known to prefer sparse solutions and overfit noise, leading to high-frequency artifacts,” and “the structure of these artifacts is sensitive to the number of RL iterations.” The artefact pattern is not fixed — it is a function of a number you chose.

Choosing an iteration count you can defend

There is no universal number, and the spread of defaults across mature implementations is itself the evidence:

• MATLAB deconvlucy: default NUMIT = 10.
• scikit-image richardson_lucy: default num_iter = 50.
• ImageJ/Fiji guidance: “Start with the default values and set iterations to 10 initially.”
• Astronomical practice: Prato et al. note RL is inefficient and that “in most cases, acceptable solutions are obtained at the cost of hundreds or thousands of iterations.”

Those differ by two orders of magnitude because the right number depends on what you are imaging and how noisy it is. The controlling variable is object sparsity. Prato et al. give the useful boundary case: for “point-wise objects such as binaries or open star clusters, iterations can be pushed to convergence,” because a sparse, spiky solution is what the object actually is. Dense, extended, textured objects — most cells, most tissue — are the opposite case, and are where early stopping matters most, because RL’s preference for sparse solutions is actively wrong about them.

A workable protocol:

1. Sweep, do not guess. Save the estimate at a geometric ladder of iterations (for example 5, 10, 20, 50, 100) on one representative field and inspect them side by side at full zoom.
2. Look for the onset of granularity, not for maximum sharpness. The stopping point is the last iteration before background regions start to break up into speckle. Maximum apparent sharpness is always past it.
3. Use an objective criterion where you can. Liu et al. propose splitting the observed image into “two noise independent halves” and using the agreement between independently deconvolved halves to determine the iteration count without eyeballing it — a noise-split approach rather than an ad-hoc cutoff.
4. Freeze it across the comparison set. Every image you will compare quantitatively must be deconvolved with the same PSF, the same algorithm and the same iteration count. Different iteration counts produce differently distributed intensities, so a per-image “optimised” count silently destroys comparability.

Because RL needs many iterations, acceleration is a standard part of production implementations rather than an exotic option. MATLAB’s deconvlucy cites Biggs and Andrews’ vector-extrapolation acceleration (Applied Optics, 1997); Prato et al. report speed-ups from a factor of about 4 to more than 30 using a scaled gradient projection method with RL-style scaling, and up to two orders of magnitude on GPU. Acceleration changes how fast you traverse the iteration path, not where the path goes — an accelerated iteration is not equivalent to one plain iteration, so an iteration count is only meaningful alongside the implementation that produced it.

Adding a real regulariser

Early stopping is regularisation by truncation. If you need to run further — or want a stopping point that is less sensitive to exactly where you stop — add an explicit prior.

Total variation (RL-TV)

The best-established option in fluorescence microscopy is total-variation regularised RL. Dey and colleagues (Microscopy Research and Technique, 2006) combined RL with a TV constraint for 3D confocal microscopy, on the grounds that confocal images are “degraded by residual out-of-focus light and by Poisson noise due to photon-limited detection” and that RL “tends to amplify noise.” TV “suppresses unstable oscillations while preserving object edges,” and they report improvement over unregularised RL on both simulated and real images, visually and quantitatively. The trade-off TV brings is its own: it favours piecewise-constant solutions, so smooth intensity gradients can be flattened into staircase-like patches.

Damping

MATLAB’s deconvlucy exposes a DAMPAR threshold (default 0, i.e. off) that suppresses further iteration for pixels whose deviation between iterations is small relative to the noise. The documented purpose is to “control noise amplification by suppressing iterations for pixels that deviate a small amount compared to the noise” — in effect, letting the algorithm keep working on structure while leaving flat background alone.

Noise-model regularisation

Commercial packages usually expose regularisation through a noise parameter rather than an iteration count. Scientific Volume Imaging’s Huygens documentation frames its signal-to-noise ratio setting as the control that determines “which fluctuations in your image intensities will be interpreted as noise or real structure,” noting that “the more noise is taken into account in deconvolution the smoother the resulting image will become.” Their published typical ranges: confocal 5–10 (low quality), 15–20 (average), above 30 (high quality); widefield 5–15 (low) and above 40 (good); cooled CCD cameras 50–100. Treat these as orientation for a vendor-specific parameter, not as universal constants — the mapping from a nominal SNR number to an amount of smoothing is implementation-specific.

Other classical regularisers

Sage et al.’s DeconvolutionLab2 paper (Methods, 2017) benchmarks RL against regularised inverse filtering, Tikhonov regularisation, Landweber, Tikhonov–Miller and fast iterative shrinkage-thresholding on 3D microscopy data, comparing image quality, performance, usability and computational cost. If you want to justify choosing RL over the alternatives for a specific dataset rather than by habit, that benchmark and its accompanying synthetic and real test cases are the place to start.

The five artefacts, and what each one is telling you

1. High-frequency granularity or speckle

Fine grain appearing first in dim, flat regions, at a scale near the resolution limit, intensifying with iteration. This is the semiconvergence signature described above: noise fitted as signal, peaking just below the diffraction cutoff. Fix: fewer iterations, or an explicit regulariser. It is not resolved structure, however plausible the spacing looks.

2. Ringing at strong edges

Bright and dark fringes flanking high-contrast boundaries — the Gibbs phenomenon, arising from representing a sharp step with a band-limited Fourier basis. MATLAB’s documentation is explicit that the output “could exhibit ringing introduced by the discrete Fourier transform used in the algorithm” and recommends applying edgetaper to the input before deconvolving. Fix: apodise the input; consider a TV regulariser, which is designed to preserve edges without oscillating around them.

3. Border artefacts

A bright or dark frame, or ripples running along the edges of the field. FFT-based convolution is circular: the algorithm assumes the image wraps around, so structure at the left edge is treated as adjacent to structure at the right edge. Bertero and Boccacci addressed this in Astronomy & Astrophysics (2005) by reconstructing the object over a domain broader than the detected image — in effect letting RL choose its own boundary conditions — and showed the reconstruction inside the domain is “considerably improved.” Fix: use an implementation with non-circulant or extended-domain boundary handling, or pad the stack generously and crop after. Never quantify inside the outer margin.

4. PSF-mismatch artefacts

Asymmetric flaring, elongation along one axis, or structure that gets worse rather than better with iteration. RL assumes the PSF is exact; a wrong PSF is a wrong forward model, and iteration drives the estimate harder into that wrong model. The most common cause in 3D microscopy is depth dependence: SVI’s documentation notes that under strong refractive-index mismatch “the PSF will then be different for different imaging depths, which is not accounted for when using a measured PSF,” and recommends optimising the setup and then deconvolving with a depth-dependent theoretical PSF.

5. Artefacts from a broken noise model

RL’s likelihood assumes Poisson photon statistics on the raw signal. Saturated (clipped) pixels are not Poisson — their true value is unknown and larger than recorded — and RL will spread the deficit into the neighbourhood as it tries to reconcile the mismatch. An uncorrected camera offset or background pedestal similarly violates the model. Fix: exclude saturated fields rather than deconvolving them; subtract the detector offset before deconvolution. MATLAB’s WEIGHT array exists partly for this, letting you exclude specific pixels from the restoration by assigning them a weight of 0.

Get the PSF right before you touch the iteration count

The PSF is the input that determines whether deconvolution helps at all. It is set by the optics — principally numerical aperture, wavelength and immersion medium, which is why objective selection is upstream of every deconvolution decision.

• One PSF per configuration. The ImageJ documentation is blunt about this: “A PSF should be obtained for every combination of z-step intervals, wavelengths, objectives, and immersion medium.” A PSF borrowed from a different channel or a different z-step is a mismatched forward model.
• Measured PSFs need correction, not just acquisition. SVI list three reasons a raw bead image is not a PSF: beads have finite size, so “it is more accurate to first correct for the bead size”; “bead images contain photon noise while true PSFs do not,” so noise correction is desirable; and a measured PSF cannot represent depth-dependent aberration.
• Theoretical PSFs are not the inferior option. They are noise-free and can be made depth-dependent. SVI note that their software “automatically corrects for spherical aberration when deconvolving with a theoretical PSF” by using different PSFs at different depths — something a single measured PSF cannot do.
• Measure the PSF anyway, as quality control. Even where you deconvolve with a theoretical PSF, a bead measurement tells you whether the system matches the model. Cole, Jinadasa and Brown’s Nature Protocols (2011) protocol is the standard reference for measuring and interpreting PSFs for confocal resolution and QC.
• Sample at Nyquist first. Deconvolution cannot recover frequencies your pixel grid never captured; undersampled data yields aliased artefacts that iteration then sharpens. The required pixel and z-step size depends on NA and wavelength, so compute it for your configuration rather than reusing a number.

Quantifying after deconvolution

Because RL preserves total counts, integrated intensity over a region is the measurement most likely to survive deconvolution intact. Peak and mean intensity are not: concentrating the same photons into a smaller volume raises peak values and changes means over fixed-size regions, by an amount that depends on the iteration count.

The practical rules that follow are unglamorous but decide whether reviewers accept the figure:

• Deconvolve every image in a comparison with identical settings, or deconvolve none of them.
• Do not derive a ratio, colocalisation coefficient or intensity threshold from images deconvolved with different parameters.
• Keep the raw stacks. Deconvolution is not reversible, and the raw data is what a reanalysis needs.
• Report it. The community checklists published by Schmied et al. in Nature Methods (2024) state that “advanced image processing steps (e.g., deconvolution, Noise2Void, CARE) may need to be indicated in the figure/figure legend, in addition to the material and methods sections,” and that “when images are quantitatively compared in an experiment, the same adjustments and processing steps must be applied.” A complete report names the software and version, the algorithm, the PSF source (theoretical or measured, and how measured), the iteration count and any regularisation parameter.

Implementations

• DeconvolutionLab2 (Sage et al., Methods 2017) — open-source, runs inside ImageJ/Fiji or standalone, hosts RL alongside several classical algorithms, and ships synthetic and real benchmark cases. The reference choice when you need a documented, citable pipeline.
• ImageJ/Fiji plugins — Iterative Deconvolve 3D and the parallel iterative deconvolution plugins, for constrained iterative deconvolution inside an existing Fiji workflow.
• MATLAB deconvlucy — RL with damping (DAMPAR), per-pixel weights (WEIGHT), sub-sampled PSF support (SUBSMPL) and Biggs–Andrews acceleration.
• scikit-image richardson_lucy — the minimal reference implementation for scripted pipelines; note clip=True by default, which thresholds output values outside [-1, 1], and filter_epsilon for guarding the division by small numbers.
• Commercial packages (Huygens and equivalents) — parameterised by noise level and quality threshold rather than raw iteration count, with depth-dependent theoretical PSFs and automated parameter selection. Convenient, but you still have to record what the automation chose.

Frequently asked questions

How many Richardson-Lucy iterations should I use?

There is no transferable number, and the defaults in mature software span an order of magnitude — 10 in MATLAB’s deconvlucy, 50 in scikit-image, hundreds to thousands in astronomical practice. Sweep a ladder of counts on one representative field, stop at the last iteration before flat background begins to break into speckle, then fix that count for the whole dataset. Sparse, point-like objects tolerate far more iterations than dense extended structures.

Does Richardson-Lucy deconvolution break the diffraction limit?

Not durably. RL does extrapolate some signal beyond the diffraction cutoff at low iteration counts, but Liu et al. (2025) show the spectral SNR gain for those frequencies drops to zero as iterations go to infinity. Treat RL as a contrast and out-of-focus-light restoration method. If you need a real resolution gain, that comes from the optics or the sample.

Why does my deconvolved image look grainy or speckled?

You have run past the semiconvergence point. RL prefers sparse solutions and begins fitting noise, producing high-frequency structure that peaks just below the diffraction cutoff and moves closer to it with each iteration. Reduce the iteration count or add a regulariser such as total variation.

Why is there a bright or dark frame around my deconvolved image?

FFT-based deconvolution treats the image as circularly wrapped, so opposite edges are convolved into each other. Use an implementation with non-circulant or extended-domain boundary handling — the approach Bertero and Boccacci (2005) formalised — or pad the stack before deconvolving and crop afterwards. Never make measurements inside the outer margin.

Do I need a measured PSF, or is a theoretical one good enough?

Both are defensible; the failure mode is a PSF that does not match the acquisition. Theoretical PSFs are noise-free and can be made depth-dependent, which matters under refractive-index mismatch. Measured PSFs capture the real system but carry photon noise, are convolved with a finite bead, and are fixed at one depth. Many labs deconvolve with a theoretical PSF and measure a bead PSF periodically as quality control.

Can I still measure fluorescence intensity after deconvolution?

Integrated intensity over a region is the safest measurement, because RL preserves the total number of counts. Peak and mean intensities change with iteration count and are not comparable across differently deconvolved images. Deconvolve every image in a comparison identically, keep the raw data, and report the parameters.

How is Richardson-Lucy different from Wiener or inverse filtering?

Wiener and regularised inverse filters are single linear steps in the Fourier domain; RL is a non-linear, iterative maximum-likelihood estimator built on a Poisson noise model. The practical differences are that RL enforces non-negativity intrinsically and conserves total flux, where a linear filter can return negative intensities that then have to be clipped — and that RL costs orders of magnitude more computation. DeconvolutionLab2 benchmarks both families on the same 3D microscopy data if you need to justify the choice for a specific dataset.

Should I use total variation regularisation by default?

It is the best-supported regulariser for fluorescence microscopy — Dey et al. (2006) showed RL-TV beats unregularised RL visually and quantitatively on confocal data. But TV assumes piecewise-constant structure, so it can flatten genuine smooth gradients into staircase patches. If your signal is dominated by graded intensity rather than discrete objects, compare TV against plain early stopping before committing.

References

• Richardson, W. H. (1972). Bayesian-based iterative method of image restoration. Journal of the Optical Society of America, 62(1), 55–59. doi:10.1364/JOSA.62.000055
• Lucy, L. B. (1974). An iterative technique for the rectification of observed distributions. The Astronomical Journal, 79, 745. doi:10.1086/111605
• Shepp, L. A., & Vardi, Y. (1982). Maximum likelihood reconstruction for emission tomography. IEEE Transactions on Medical Imaging, 1(2), 113–122. doi:10.1109/TMI.1982.4307558
• Liu, Y., Panezai, S., Wang, Y., & Stallinga, S. (2025). Noise amplification and ill-convergence of Richardson-Lucy deconvolution. Nature Communications, 16, 911. PMC11751374
• Dey, N., Blanc-Féraud, L., Zimmer, C., Roux, P., Kam, Z., Olivo-Marin, J.-C., & Zerubia, J. (2006). Richardson-Lucy algorithm with total variation regularization for 3D confocal microscope deconvolution. Microscopy Research and Technique, 69(4), 260–266. PubMed 16586486
• Bertero, M., & Boccacci, P. (2005). A simple method for the reduction of boundary effects in the Richardson-Lucy approach to image deconvolution. Astronomy & Astrophysics, 437(1), 369–374. doi:10.1051/0004-6361:20052717
• Prato, M., Cavicchioli, R., Zanni, L., Boccacci, P., & Bertero, M. (2012). Efficient deconvolution methods for astronomical imaging: algorithms and IDL-GPU codes. Astronomy & Astrophysics, 539, A133. arXiv:1210.2258
• Biggs, D. S. C., & Andrews, M. (1997). Acceleration of iterative image restoration algorithms. Applied Optics, 36(8), 1766–1775. doi:10.1364/AO.36.001766
• Sage, D., Donati, L., Soulez, F., Fortun, D., Schmit, G., Seitz, A., Guiet, R., Vonesch, C., & Unser, M. (2017). DeconvolutionLab2: An open-source software for deconvolution microscopy. Methods, 115, 28–41. doi:10.1016/j.ymeth.2016.12.015
• Cole, R. W., Jinadasa, T., & Brown, C. M. (2011). Measuring and interpreting point spread functions to determine confocal microscope resolution and ensure quality control. Nature Protocols, 6, 1929–1941. doi:10.1038/nprot.2011.407
• Schmied, C., Nelson, M. S., Avilov, S., et al. (2024). Community-developed checklists for publishing images and image analyses. Nature Methods, 21, 170–181. PMC10922596
• Hendrix, Z., Brown, A., Flanagan, T., Shepherd, D. P., Saurabh, A., & Pressé, S. (2024). Re-thinking Richardson-Lucy without iteration cutoffs: physically motivated Bayesian deconvolution (preprint). arXiv:2411.00991
• MathWorks. deconvlucy — Image Processing Toolbox reference. mathworks.com
• scikit-image. skimage.restoration.richardson_lucy API reference. scikit-image.org
• ImageJ. Deconvolution. imagej.net
• Scientific Volume Imaging. Signal to Noise Ratio. svi.nl
• Scientific Volume Imaging. Point spread function: theoretical vs. experimental. svi.nl