# Adapting Richardson Lucy (RL) Deconvolution for Shot Noise Limited Coherent Imaging

I am an experimental physicist who is collecting a series of coherent imaging of trapped gas. If you are familiar with phase contrast imaging, you may understand what I mean by coherent imaging. The information about the sample is retrieved as interference between a phase-shifted version illumination and light scattered by the sample. The light I am using is monochromatic. Since the imaging technique is coherent, the point-spread function (PSF) I need to use is $$J_1(x)/x$$ instead of the usual $$(J_1(x)/x)^2$$ for a circular aperture. This means that PSF has both positive and negative values. (It may be even complex in real-life situations.) I will describe the details of what I actually observe in an image below.

Also, since I am working with a low light-level condition, I thought it is best to use the Richardson-Lucy (RL) algorithm, which seems tailored for deconvolving images with Poisson noise. My understanding of the RL algorithm is based on Gordon Wetzstein - Stanford EE367/CS448I: Computational Imaging and Display - Lecture 10 (slide numbers 26-30). The algorithm assumes that you have an input $$\vec{x}$$, which goes through a linear operator $$\hat{A}$$ (convolution with PSF), and the output is a Poisson random variable $$\vec{b}$$ with average $$\langle\vec{b}\rangle = \hat{A}\vec{x}$$.

The assumptions used for the RL algorithm are consistent with typical incoherent imaging, where the PSF and the image are everywhere non-negative, but they are not in general for coherent imaging which may have negative values (because negative values wouldn't make much sense for a Poisson distribution).

Question: would there be a way to adopt the RL algorithm for coherent imaging, or more in general, PSF with negative values?

More detailed description of a typical image. $$I$$ denotes intensity. $$E$$ denotes complex electric field, and $$I = |E|^2/2$$.

\begin{align*} I_{\text{image}} &= \frac{|E_{\text{scattered}} + i E_{\text{illumination}}|^2}{2} - I_{\text{illumination}}\\ &= I_{\text{scattered}} + \mathrm{Im}(E_{\text{scattered}} E_{\text{illumination}}^*)\\ E_{\text{illumination}} &= \text{plane wave}\\ E_{\text{scattered}} &\sim\text{spherical wave focused to a spot}\\ |E_{\text{scattered}}| &<< |E_{\text{illumination}}| \end{align*} the last condition allows us to ignore $$I_{\text{scattered}}$$ and focus on the $$\mathrm{Im}(E_{\text{scattered}} E_{\text{illumination}}^*)$$. Since $$E_{\text{illumination}}$$ is more-or-less flat, you could say we are imaging a magnified version of $$E_{\text{scattered}}$$ times its relative phase shift $$\phi$$ (small) with respect to the illumination.

\begin{align*} E_{\text{scattered}} &= E_{\text{illumination}} (e^{i\phi_{gas}}-1)\\ &\approx E_{\text{illumination}}i\phi_{gas}\\ I_{\text{image}}[x,y] &\approx |E_{\text{illumination}}|^2\cdot\phi_{gas}[x,y] \end{align*}

Here is a sample image that shows fringes. • What about adding large value to the output so all values are non negative, apply the algorithm both to the image and to the added value and then remove it?
– Royi
Oct 31 '20 at 13:30
• @Royi, I think the problem is that it is unclear what to do with the normalization factor in the RL algorithm. See the term $\mathrm{diag}(A^T\,1)$ in bigwww.epfl.ch/deconvolution/challenge2013/…. EDIT: never mind, I think as long as the integral of kernel over all space is positive, it should be okay.
– wcc
Nov 2 '20 at 18:01