理查森-露西算法是如何工作的？代码示例？

Question

我想弄清楚反褶积是如何工作的。我理解它背后的想法，但我想了解实现它的一些实际算法-算法，将带有点样本函数（模糊核）的模糊图像作为输入，并生成作为输出的潜在图像。

到目前为止，我发现了Richardson–Lucy算法，其中的数学似乎没有那么困难，但我无法理解实际的算法是如何工作的。维基百科上说：

This leads to an equation for which can be solved iteratively according...

但是它没有显示实际的循环。有谁能给我指一个解释实际算法的资源吗。在Google上，我只能找到使用Richardson-Lucy作为步骤之一的方法，而不是实际的Richardson-Lucy算法。

任何语言或伪代码中的算法都是不错的，但是如果Python中有算法，那就太棒了。

先走一步。

编辑

基本上我想知道的是模糊图像（nxm）：

x00 x01 x02 x03 .. x0n
x10 x11 x12 x13 .. x1n
...
xm0 xm1 xm2 xm3 .. xmn

以及用于获取模糊图像的内核（ixj）：

p00 p01 p02 .. p0i
p10 p11 p12 .. p1i
...
pj0 pj1 pj2 .. pji

为了找出原始图像，理查森-露西算法的具体步骤是什么。

Answer 1

下面是一个非常简单的Matlab实现：

function result = RL_deconv(image, PSF, iterations)
    % to utilise the conv2 function we must make sure the inputs are double
    image = double(image);
    PSF = double(PSF);
    latent_est = image; % initial estimate, or 0.5*ones(size(image)); 
    PSF_HAT = PSF(end:-1:1,end:-1:1); % spatially reversed psf
    % iterate towards ML estimate for the latent image
    for i= 1:iterations
        est_conv      = conv2(latent_est,PSF,'same');
        relative_blur = image./est_conv;
        error_est     = conv2(relative_blur,PSF_HAT,'same'); 
        latent_est    = latent_est.* error_est;
    end
    result = latent_est;

original = im2double(imread('lena256.png'));
figure; imshow(original); title('Original Image')

hsize=[9 9]; sigma=1;
PSF = fspecial('gaussian', hsize, sigma);
blr = imfilter(original, PSF);
figure; imshow(blr); title('Blurred Image')

res_RL = RL_deconv(blr, PSF, 1000); 
figure; imshow(res_RL); title('Recovered Image')

你也可以在频域工作，而不是像上面那样在空域工作。在这种情况下，代码是：

function result = RL_deconv(image, PSF, iterations)
fn = image; % at the first iteration
OTF = psf2otf(PSF,size(image)); 
for i=1:iterations
    ffn = fft2(fn); 
    Hfn = OTF.*ffn; 
    iHfn = ifft2(Hfn); 
    ratio = image./iHfn; 
    iratio = fft2(ratio); 
    res = OTF .* iratio; 
    ires = ifft2(res); 
    fn = ires.*fn; 
end
result = abs(fn); 

唯一我不太明白的是PSF的空间反转是如何工作的以及它的作用是什么。如果有人能解释给我听那就太酷了！我也在寻找一个简单的Matlab R-L实现空间变量psf（即空间非齐次点扩散函数）-如果有人会有一个请让我知道！

为了去掉边上的人工制品，可以在边上镜像输入图像，然后裁剪镜像的位，或者在调用RL_deconv之前使用Matlab的image = edgetaper(image, PSF)。

本地的Matlab实现deconvlucy.m有点复杂，顺便说一句，它的源代码可以找到here，并使用accelerated version of the basic algorithm。

Answer 2

下面是一个开源的Python实现： http://code.google.com/p/iocbio/wiki/IOCBioMicroscope

Answer 3

如果这有帮助的话，我写了一个实现，其中包括一些文档。。。。

https://github.com/bnorthan/projects/blob/master/truenorthJ/ImageJ2Plugins/functions/src/main/java/com/truenorth/functions/fft/filters/RichardsonLucyFilter.java

Richardson Lucy是许多其他反褶积算法的基础。例如，上面的iocbio示例修改了算法以更好地处理噪声。

这是一个相对简单的算法（在这些情况下），也是更复杂算法的起点，因此您可以找到许多不同的实现。