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 <code>x</code> given a vector <code>y</code> in the matrix equation <code>y=Ax</code>. Written in terms of elements, this is:

<code>y_i = sum_j A_ij x_j</code>

In tomography, the <code>x_j</code> are cells in space that emit radiation, and the <code>y_i</code> are "bins" of photons collected at particular angles and locations. In deblurring, <code>x</code> is the original image and <code>y</code> is the blurry image.

The premise behind Richardson-Lucy is that each contribution from an element of <code>x</code> to an element of <code>y</code> has independent Poisson noise:

<code>y_i = sum_j B_ij qquad B_ij sim textPoisson(A_ijx_j)</code>

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

<code>hatx_j = fracsum_i B_ijsum_i A_ij</code>

For the "E" step, we find the expected value of <code>B_ij</code> given <code>x</code>. The calculation amounts to distributing the value of each <code>y_i = sum_j B_ij</code> amongst the <code>i</code>th row of <code>B_ij</code> in proportions according to the known values of <code>A_ijx_j</code>:

<code>B_ij = y_i fracA_ijx_jsum_k A_ikx_k = fracy_i(Ax)_i A_ijx_j</code>

So, here is a way to envision the R-L algorithm at work. Start with a matrix whose elements are <code>A_ijx_j</code> for some initial guess of <code>x</code>. Then independently scale each row so that it sums to <code>y_i</code> (turning the matrix into an estimate of <code>B_ij</code>). Then, replace each column with a multiple of the corresponding column of <code>A</code> 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 <code>y</code>, and each column should be a multiple of the same column in <code>A</code>. When the process converges, our matrix contains valid values of <code>A_ijx_j</code>.

To show that this actually is R-L, we plug <code>(4)</code> into <code>(3)</code>. We assume that <code>A</code> is energy-preserving, so that <code>sum_i A_ij</code> is <code>1</code>.

<code>x_j = sum_i fracy_i(Ax)_iA_ijx_j = x_jsum_i(A^T)_jifracy_i(</code>

<code>x = x ast A^TfracyAx</code>

In the last equation, for notational convenience, the <code>ast</code> 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 <code>x</code>, 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 <code>x_j</code> 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

View texdown source