图像去模糊（二）——Richardson–Lucy算法

1、基本模型假设
Y = Λ ∗ P
其中：
Y：退化图像，
Λ : 原始图像，
P : 点扩散函数，
∗ : 卷积操作。
　　这里的点扩散函数（Point Spread Function，PSF）描述了成像系统对点光源或点状物的响应。

2、算法思想
在假设模型

Y = Λ ∗ P 中，我们有：

简单起见，可以假设d=n，则观测图像和真实图像包含数量相同的像素点。
根据yj的物理意义知其满足泊松分布。事实上，

其中

3、求解

3.1 泊松分布
泊松分布描述了一段时间内某件事发生一定次数的概率。事件需要满足以下条件：
1）每次事件之间相互独立；
2）事件发生的概率在一段时间内是稳定的。
泊松分布的概率函数为：

参数λ是单位时间(或单位面积)内随机事件的平均发生率，也是泊松分布的期望和方差。
进一步的理解和应用实例可以参看博文：
http://www.ruanyifeng.com/blog/2013/01/poisson_distribution.html

3.2 迭代式推导
采用高斯型点扩展函数（实际问题中扩散函数形式是多种多样的，此处以该函数模型为例）：

D(i,j)表示点i和j的距离。
由前面的分析知yj满足泊松分布，得：

该式表示了点i在点j上的映射贡献，故yj可由下式表示：

当图像Λ已知时，z(i, j)可由下式估计：

而当z(i,j)已知时，λ可由下式估计：

因此最终的迭代式为：

更详细的求解过程用到了极大似然法和EM算法，这两个概念可以参看博文：
http://blog.csdn.net/zouxy09/article/details/8537620

4、去卷积效果实例
原图：

高斯模糊后：

5次迭代恢复：

20次迭代恢复：
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-4EtsvP2m-1595927090226)(https://img-blog.csdn.net/20161017223930211)]
可以看到，迭代次数越多，图像边界区域的振铃效应越明显，但图像细节恢复越清晰。
从整体上定量来看，迭代次数在某个值附近达到最佳效果，这个值因图而异。本例定量计算信噪比(SNR)及峰值信噪比(PSNR)结果如下图所示，大约在14次迭代时结果最优。图中横坐标为迭代次数，纵坐标单位为db：

附上实验用到的matlab代码：

clc;

I = imread('original.JPG');
figure(1);
imshow(I);
title('Original Image');

PSF = fspecial('gaussian',5,5);
Blurred = imfilter(I,PSF,'symmetric','conv');
figure(2);
imshow(Blurred);
title('Gaussian Blurred');
imwrite(Blurred,'高斯模糊.png');

lucy = deconvlucy(Blurred,PSF,5);
figure(3);
imshow(lucy);
title('Restored Image, NUMIT = 5');
imwrite(lucy,'5次迭代.png');

lucy = deconvlucy(Blurred,PSF,20);
figure(4);
imshow(lucy);
title('Restored Image, NUMIT = 20');
imwrite(lucy,'20次迭代.png');

function BestRL
    clc;
    generateImg();%
    calcPSNR();
end

function generateImg
    I = imread('original.JPG');

    PSF = fspecial('gaussian',5,5);
    Blurred = imfilter(I,PSF,'symmetric','conv');
    %imwrite(Blurred,'高斯模糊.png');

    for i=1:1:20
        rI = deconvlucy(Blurred,PSF,i);
        imwrite(rI,sprintf('%d.png',i));
    end
end

function calcPSNR
    PSNR = 1:1:20;
    SNR = 1:1:20;
    I = imread('original.JPG');
    for i=1:1:20
        Blurred = imread(sprintf('%d.png',i));
        [PSNR(i),SNR(i)] = psnr(Blurred,I);
    end
    figure(1);
    plot(PSNR);
    title('PSNR');
    
    figure(2);
    plot(SNR);
    title('SNR');
end