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

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

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

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(即空间非齐次点扩散函数)——如果有人想要，请告诉我!

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

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