Re-thinking Richardson-Lucy without Iteration Cutoffs: Physically Motivated Bayesian Deconvolution

TL;DR

Bayesian deconvolution is presented, a rigorous deconvolution framework that combines a physically accurate image formation model avoiding the challenges inherent to the RL approach and is amenable to fast, parallelizable computation.

Abstract

Richardson-Lucy deconvolution is widely used to restore images from degradation caused by the broadening effects of a point spread function and corruption by photon shot noise, in order to recover an underlying object. In practice, this is achieved by iteratively maximizing a Poisson emission likelihood. However, the RL algorithm is known to prefer sparse solutions and overfit noise, leading to high-frequency artifacts. The structure of these artifacts is sensitive to the number of RL iterations, and this parameter is typically hand-tuned to achieve reasonable perceptual quality of the inferred object. Overfitting can be mitigated by introducing tunable regularizers or other ad hoc iteration cutoffs in the optimization as otherwise incorporating fully realistic models can introduce computational bottlenecks. To resolve these problems, we present Bayesian deconvolution, a rigorous deconvolution framework that combines a physically accurate image formation model avoiding the challenges inherent to the RL approach. Our approach achieves deconvolution while satisfying the following desiderata: I deconvolution is performed in the spatial domain (as opposed to the frequency domain) where all known noise sources are accurately modeled and integrated in the spirit of providing full probability distributions over the density of the putative object recovered; II the probability distribution is estimated without making assumptions on the sparsity or continuity of the underlying object; III unsupervised inference is performed and converges to a stable solution with no user-dependent parameter tuning or iteration cutoff; IV deconvolution produces strictly positive solutions; and V implementation is amenable to fast, parallelizable computation.

Figures (5)

• Figure 1: Bayesian deconvolution alternative to Richardson-Lucy. a. Image formation model where: 1) the imaged object $\rho$ convolved with the PSF to produce an expected convolved image $\mu$; 2) Poisson distributed photon counts $\varphi$ were replicated by corrupting the image with shot noise; and finally 3) normally distributed measurements were reproduced by further degrading with camera noise and pixelization. b. Inverse strategy where samples are generated from a posterior probability distribution using a Monte Carlo scheme. c. Mean of N Monte Carlo samples drawn from the posterior along with their Fourier spectra. Red circles represent the diffraction limit. d. Increasing peak sign-to-noise ratio (PSNR) and decreasing root mean square error (RMSE) as a function of the number of collected samples is desired and achieved by the method proposed herein.
• Figure 2: Deconvolution of a simulated image.a. Ground truth (GT) and simulated widefield raw image (WF) assuming Poisson noise only with peak pixel photon count of 150 approximately. Yellow circles represent the diffraction limit in terms of inter-spoke spacing. The scale bar represents $0.5\upmu \mathrm{m}$. b. RL deconvolution (RL) with increasing iterations (I) from left to right. c. Deconvolution by our Bayesian algorithm (B) as the number of samples (N) used to compute the mean increases. d-e. PSNR and RMSE for RL (left) and Bayesian deconvolutions (right). We note optimal RL performance at around 10 iterations fo this particular sample before both PSNR and RMSE degrade motivating efforts to seek an iteration cutoff in traditional RL.
• Figure 3: Deconvolution of mitochondrial network in HeLa cells labeled with Mitotracker Deep Red dye.a. Raw image after subtracting average sCMOS camera offset of 100 ADU and then dividing by the average gain of 2 ADU/e$^{-}$, showing approximate photon counts. Scale bar $5.0\upmu$m. b. The leftmost panel in the upper row shows the widefield raw image (WF) from the $150\times150$ pixel inset in a. The inset lies near the edge of illumination. Scale bar $1.5\upmu$m. This image has been brightened compared to a to show structures more clearly. The remaining images in the upper row show deconvolution by RL algorithms (RL) as the number of iterations (I) increases. The panels in the row below show deconvolution by our Bayesian algorithm (B) as the number of MCMC samples (N) used to compute the mean increases. c. PSNR computed for the Bayesian deconvolved image at $N=100$ with RL deconvolved images by iteration I as reference, demonstrating similarity of converged Bayesian and optimal RL results.
• Figure 4: Coefficient of Variation.a. Restored image from Figure \ref{['fig:siemens']}c obtained using Bayesian deconvolution with 100 MCMC samples. Scale bar $0.5\upmu$m. b. Ratio of standard deviation and mean for the Bayesian deconvolved image in the inset, providing a relative measure of uncertainty. As expected, the uncertainty is higher where information content (photon counts) is low.
• Figure 5: Parallelized Gibbs sampling.a. PSF and Raw Image as input for deconvolution. b. Padding of the size of the PSF is added to the image to perform convolutions conveniently. c. Padded image is chunked. Padding in the image interior overlaps with neighboring chunks and exchanged after each MCMC iteration. d. At each MCMC iteration, Gibbs sampling requires sequentially updating each pixel using convolutions of the object in the red box with PSF. Each updated pixel is immediately used for updating remaining pixels during the iteration. If Gibbs sampling is initiated in a corner chunk, it can be parallelized by exchanging pixels updated at iteration $I=n$ with edge-sharing neighbor chunks and immediately moving on to the next iteration $I=n+1$ while the neighboring chunks are now being updated at iteration $I=n$. This procedure results in wavefront parallelization where all chunks are being updated simultaneously but are at different MCMC iterations.