I'm working on a feature for a Free Software astronomical image processing program. I want to implement a damped Richardson-Lucy deconvolution as described by Richard L White (Journal of the Space Telescope Science Institute, 1994). He offers the following term for each successive iteration:
$$O_{new}(j) = O(j)\frac{\left(\sum(i)P(i|j)(U_i^N-1(N-(N-1)U_i)\frac{D(i)-I(i)}{I(i)}\right)}{\sum(i)P(i|j)}$$
where: $O$ is the un-blurred object, $I$ is the blurry image, $D$ is the observed image, $N$ is a constant, $U_i$ is a formula depending only on per-element values $D(i)$ and $I(i)$.
I want to carry out this algorithm using FFTs rather than naive convolutions, both for speed and as I already have a perfectly functioning standard R-L implementation using FFTs, so I should have all the building blocks I need.
I can make the link between the basic R-L formulation expressed as $$O_{new}(j) = O(j)\frac{\sum(i)P(i|j)\frac{D(i)}{I(i)}}{\sum(i)P(i|j)}$$ and its expression using using FFTs as $$O_{new} = O\cdot \texttt{IFFT}\left(\frac{D}{\texttt{IFFT}(\texttt{FFT}(O)\cdot\texttt{FFT}(P))}\cdot\texttt{FFT}(P^\star)\right)$$
But I'm completely failing to understand how to rewrite the more complicated formula for White's damped version in terms of FFTs and inverse FFTs (and element-wise operations on the image pixels). Sorry if it's a bit of a basic question: I've read up as much as I can find and have asked elsewhere but without success, so please could someone show me how to do this?
