Richardson-Lucy算法如何工作？代码示例？

Question

我想弄清楚解卷积是如何工作的.我理解它背后的想法,但我想了解一些实现它的实际算法 - 算法将带有点样本函数(模糊内核)的模糊图像作为输入,并产生潜像的输出.

到目前为止,我发现了Richardson-Lucy算法,其中数学似乎并不那么困难,但我无法确定实际算法是如何工作的.在维基百科,它说:

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

但是它没有显示实际的循环.任何人都可以指向我解释实际算法的资源.在Google上,我只设法找到使用Richardson-Lucy作为其中一个步骤的方法,而不是实际的Richardson-Lucy算法.

任何语言或伪代码中的算法都会很好,但是如果有一个在Python中可用,这将是惊人的.

Thanx提前.

编辑

基本上我想弄清楚的是模糊图像(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); toc;
figure; imshow(res_RL2); 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 RL实现方案，用于空间变化的PSF（即空间非均匀点扩散函数）-如果有人可以，请告诉我！

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

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

Answer 2

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