Richardson-Lucy 算法是如何工作的？代码示例？答案

【问题标题】：How does Richardson–Lucy algorithm work? Code example?Richardson-Lucy 算法是如何工作的？代码示例？
【发布时间】：2012-03-24 18:27:32
【问题描述】：

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

到目前为止，我发现 Richardson–Lucy 算法的数学似乎并不那么困难，但是我无法弄清楚实际算法是如何工作的。在维基百科它说：

这导致一个方程可以根据...迭代求解...

但是它没有显示实际的循环。谁能指出我解释实际算法的资源。在 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

为了找出原始图像，Richardson-Lucy 算法的具体步骤是什么。

【问题讨论】：

标签： python algorithm image-processing convolution

【解决方案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。

【讨论】：

当您考虑像 [0,0,1] 这样的内核时，空间反转的需求变得显而易见，它将整个图像向右移动一个像素而不会使它模糊。然后必须将误差估计反馈给左邻居，这是通过反转内核来实现的。

【解决方案2】：

Wikipedia 上的方程给出了迭代 t+1 的迭代函数 t。您可以通过以下方式实现这种类型的迭代算法：

def iter_step(prev):
  updated_value = <function from Wikipedia>
  return updated_value

def iterate(initial_guess):
  cur = initial_guess
  while True:
    prev, cur = cur, iter_step(cur)
    if difference(prev, cur) <= tolerance:
      break
  return cur

当然，您必须实现自己的difference 函数，该函数对于您正在处理的任何类型的数据都是正确的。您还需要处理从未达到收敛的情况（例如限制迭代次数）。

【讨论】：

【解决方案3】：

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

【讨论】：

【解决方案4】：

如果有帮助的话，我编写的一个实现包含一些文档......

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

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

这是一个相对简单的算法（就像这些事情一样），并且是更复杂算法的起点，因此您可以找到许多不同的实现。

【讨论】：