Richardson-Lucy算法如何工作？代码示例？ [英] How does Richardson–Lucy algorithm work? Code example?

问题描述

我试图找出解卷积的工作原理。我理解它背后的想法，但我想了解一些实现它的实际算法 - 算法将带有点样本函数（模糊内核）的模糊图像作为输入，并产生潜像的输出。

到目前为止，我找到了

  hsize = [9 9];西格玛= 1; 
 PSF = fspecial（'gaussian'，hsize，sigma）; 
 blr = imfilter（original，PSF）; 
数字; imshow（BLR）;标题（'模糊图像'）
  

  res_RL = RL_deconv（blr，PSF，1000）; TOC; 
数字; imshow（res_RL2）; title（'Recovered Image'）
  

您也可以在频域中工作，而不是如上所述在空间域中工作。在这种情况下，代码将是：

 函数结果= RL_deconv（图像，PSF，迭代）
 fn = image ;第一次迭代时的％
 OTF = psf2otf（PSF，size（image））; i = 1的
：迭代
 ffn = fft2（fn）; 
 Hfn = OTF。* ffn; 
 iHfn = ifft2（Hfn）; 
 ratio = image./iHfn; 
 iratio = fft2（比率）; 
 res = OTF。* iratio; 
 ires = ifft2（res）; 
 fn = ires。* fn; 
结束
结果= abs（fn）; 
  

我唯一不太了解的是PSF的这种空间逆转是如何工作的以及它的用途是什么。如果有人能为我解释那会很酷！我也在寻找一个简单的Matlab RL实现空间变异PSF（即空间非均匀点扩散函数） - 如果有人有一个请告诉我！

To摆脱边缘处的人工制品，你可以在边缘镜像输入图像，然后在之后裁掉镜像的位或使用Matlab的 image = edgetaper（图像，PSF）之前你打电话给 RL_deconv 。

原生Matlab实现deconvlucy.m有点复杂btw - 源代码可以找到此处并使用基本算法的加速版本。

I am trying to figure out how deconvolution works. I understand the idea behind it but I want to understand some of the actual algorithms which implement it - algorithms which take as input a blurred image with its point sample function (blur kernel) and produce as output the latent image.

So far I found Richardson–Lucy algorithm where the math does not seem to be that difficult however I can't figure how the actual algorithm works. At Wikipedia it says:

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

however it does not show the actual loop. Can anyone point me to a resource where the actual algorithm is explained. On Google I only manage to find methods which use Richardson–Lucy as one of its steps, but not the actual Richardson–Lucy algorithm.

Algorithm in any language or pseudo-code would be nice, however if one is available in Python, that would be amazing.

Thanx in advance.

Edit

Essentially what I want to figure out is given blurred image (nxm):

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

and the kernel (ixj) which was used in order to get the blurred image:

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

What are the exact steps in the Richardson–Lucy algorithm in order to figure out the original image.

解决方案

Here is a very simple Matlab implementation :

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); toc;
figure; imshow(res_RL2); title('Recovered Image')

You can also work in the frequency domain instead of in the spatial domain as above. In that case the code would be :

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); 

Only thing I don't quite understand is how this spatial reversal of the PSF works and what it's for. If anyone could explain that for me that would be cool! I'm also looking for a simple Matlab R-L implementation for spatially variant PSFs (ie spatially nonhomogeneous point spread functions) - if anyone would have one please let me know!

To get rid of the artefacts at the edges you could mirror the input image at the edges and then crop away the mirrored bits afterwards or use Matlab's image = edgetaper(image, PSF) before you call RL_deconv.

The native Matlab implementation deconvlucy.m is a bit more complicated btw - the source code of that one can be found here and uses an accelerated version of the basic algorithm.

这篇关于Richardson-Lucy算法如何工作？代码示例？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！