NEWS FAQS Overview Registration Downloads UPLOADS Documentation Contact

RICHARDSON-LUCY ALGORITHM

Summary

Here we clarify a practical implication of the fact that our forward model relies on Toeplitz boundary conditions. To illustrate our point, we rederive the standard Richardson-Lucy algorithm. The same type of modification may be required for your own algorithm if you originally developed it for circulant boundary conditions.

Derivation of the algorithm

The Richardson-Lucy algorithm is a maximum-likelihood algorithm that is based on a Poisson noise model. Using matrix-vector notations, we have

$$\mathbf{y} \sim \mathcal{P}(\mathbf{Ax} + \mathbf{b})$$ where $\mathbf{x}$ is a vector corresponding to the ground-truth image, $\mathbf{A}$ is a matrix representing the forward model, $\mathbf{b}$ is a constant vector representing the background signal and $\mathbf{y}$ is a random vector modeling the measurements. The corresponding Poisson likelihood is $$p(\mathbf{y}|\mathbf{x}) = \exp\big({-}(\mathbf{Ax} + \mathbf{b})^T\mathbf{1}\big) \times \exp\big(\log(\mathbf{Ax} + \mathbf{b})^T\mathbf{y}\big) \times \prod_{n=1}^N 1/\mathbf{y}_n!,$$ where the logarithm function is applied component-wise and $N$ is the dimension of the measurement vector $\mathbf{y}$.

Maximizing $p(\mathbf{y}|\mathbf{x})$ with respect to $\mathbf{x}$ is equivalent to maximizing the log-likelihood function

$$L(\mathbf{x}) = \log\big(p(\mathbf{y}|\mathbf{x})\big).$$ The gradient of this function is $$\nabla L(\mathbf{x}) = \mathbf{A}^T \mathrm{diag}(\mathbf{Ax}+\mathbf{b})^{-1} \mathbf{y} - \mathbf{A}^T\mathbf{1}.$$ Imposing that this quantity vanishes leads to the following multiplicative-update algorithm: $$\mathbf{x}_{k+1} = \mathrm{diag}(\mathbf{A}^T\mathbf{1})^{-1} \mathrm{diag}\big[\mathbf{A}^T \mathrm{diag}(\mathbf{Ax}_k+\mathbf{b})^{-1} \mathbf{y}\big] \mathbf{x}_k.$$

Take-home message

Usually the matrix $\mathbf{A}$ is taken to be block-circulant and normalized such that $\mathbf{A}^T\mathbf{1} = \mathbf{1}$; in this case the Richardson-Lucy iteration reduces to $\mathbf{x}_{k+1} = \mathrm{diag}\big[\mathbf{A}^T \mathrm{diag}(\mathbf{Ax}_k+\mathbf{b})^{-1} \mathbf{y}\big] \mathbf{x}_k$. However, in the framework of this challenge the matrix $\mathbf{A}$ has a block-Toeplitz structure and therefore the normalization factor $\mathrm{diag}(\mathbf{A}^T\mathbf{1})^{-1}$ is non-trivial. We refer to our reference implementation of the Richardson-Lucy algorithm (RLdeblur3D.m) for more details.

Important Dates