# Richardson-Lucy Deconvolution

## Richard-Lucy as an EM algorithm

This paper presents an EM algorithm for reconstruction from emission
tomography data that is exactly Richardson-Lucy:

> Green, Peter J.  "Bayesian Reconstructions from Emission Tomography
  Data Using a Modified EM Algorithm".  *IEEE Transactions on Medical
  Imaging*, vol. 9, pp. 84-93, 1990.

Suppose we wish to reconstruct a vector `%x%` given a vector `%y%` in
the matrix equation `%y=Ax%`.  Written in terms of elements, this is:

     %% y_i = \sum_j A_{ij} x_j

In tomography, the `%\{x_j\}%` are cells in space that emit radiation,
and the `%\{y_i\}%` are "bins" of photons collected at particular
angles and locations.  In deblurring, `%x%` is the original image and
`%y%` is the blurry image.

The premise behind Richardson-Lucy is that each contribution from an
element of `%x%` to an element of `%y%` has independent Poisson noise:

     %% y_i = \sum_j B_{ij}, \qquad B_{ij} \sim \text{Poisson}(A_{ij}x_j)

For the "M" step, we find the most likely value of `%x%` based on
`%B%`.  The most likely value of `%x_j%` depends on the `%j%`th column
of `%B%`, which in the absence of noise would be `%x_j%` times the
corresponding column of `%A%`:

     %% \hat{x}_j = \frac{\sum_i B_{ij}}{\sum_i A_{ij}}

For the "E" step, we find the expected value of `%B_{ij}%` given
`%x%`.  The calculation amounts to distributing the value of each
`%y_i = \sum_j B_{ij}%` amongst the `%i%`th row of `%B_{ij}%` in
proportions according to the known values of `%A_{ij}x_j%`:

     %% B_{ij} = y_i \frac{A_{ij}x_j}{\sum_k A_{ik}x_k}
     %%   = \frac{y_i}{(Ax)_i} A_{ij}x_j

So, here is a way to envision the R-L algorithm at work.  Start with a
matrix whose elements are `%A_{ij}x_j%` for some initial guess of
`%x%`.  Then independently scale each row so that it sums to `%y_i%`
(turning the matrix into an estimate of `%B_{ij}%`).  Then, replace
each column with a multiple of the corresponding column of `%A%`
having the same sum, and repeat the process.  In other words, we
alternate between enforcing two constraints: summing across the rows
of our matrix should yield `%y%`, and each column should be a multiple
of the same column in `%A%`.  When the process converges, our matrix
contains valid values of `%A_{ij}x_j%`.

To show that this actually is R-L, we plug `%(4)%` into `%(3)%`.  We
assume that `%A%` is energy-preserving, so that `%\sum_i A_{ij}%` is `%1%`.

     %% x'_j = \sum_i \frac{y_i}{(Ax)_i}A_{ij}x_j
     %% = x_j\sum_i(A^T)_{ji}\frac{y_i}{(Ax)_i}

     %% x' = x \ast A^T\frac{y}{Ax}

In the last equation, for notational convenience, the `%\ast%`
represents element-wise multiplication, and the fraction bar
represents element-wise division.

**Ringing:** Given the "matrix-update" formulation of R-L given above,
we can see how it is easy for ringing to come about.  The ringing
arises in the reconstructed image `%x%`, each element of which
corresponds to the scale of a column of the matrix.  Basically,
neither update step hinders ringing at all; operating on columns of
the matrix does not relate the `%\{x_j\}%` to each other, and when
operating on the rows, only the sum of each row is considered, which
ignores the oscillation of ringing.

## Acceleration

Performing basic R-L as described above converges exceedingly slowly
and will not match the quality of Matlab's deconvlucy.  Matlab uses an
acceleration technique which essentially extends each step by a small
factor.  The Matlab source code references the following papers:

    "Acceleration of iterative image restoration algorithms, by D.S.C. Biggs 
    and M. Andrews, Applied Optics, Vol. 36, No. 8, 1997.

    "Deconvolutions of Hubble Space Telescope Images and Spectra",
    R.J. Hanisch, R.L. White, and R.L. Gilliland. in "Deconvolution of Images 
    and Spectra", Ed. P.A. Jansson, 2nd ed., Academic Press, CA, 1997.

The acceleration technique is described in the first paper.  Matlab
implements the algorithm from sections 2C-D of this paper, and it is
not a hard algorithm to implement yourself.

-----

[Back](index.html)

[View texdown source](deconvlucy.text)