本发明涉及机器视觉领域，特别提出了一种基于超拉普拉斯先验的图像盲去模糊方法。
背景技术：
：
图像作为人类获取信息的重要来源，在现代社会中占据着重要的地位。然而对于手持拍摄设备而言，在曝光的过程中经常会发生抖动现象，导致获取的图像变得模糊，这严重影响图像的使用及其后续的处理。图像去模糊的目的就是仅从退化了的模糊图像中恢复出清晰的、包含丰富信息的图像。模糊核是指在曝光的过程中拍摄设备因发生抖动而产生的运动轨迹。根据模糊核是否已知可将图像去模糊过程分为盲去模糊和非盲去模糊。其中，图像盲去模糊是在模糊核无法获知情形下的去模糊过程，是一个严重的不适定问题，具有很大的挑战性。图像盲去模糊最主要的任务是估计出模糊核，之后的问题就转化为了非盲去模糊问题，最终基于非盲去模糊方法实现图像的复原。图像去模糊作为数字图像处理和计算机视觉领域的重要研究课题，在摄影学、光学、天文学、医疗图像、监控、遥感探测和军事研究等领域都具有广泛的应用价值，也是近年来学术界普遍关注的一个热点课题，具有非常重要的理论意义和现实意义。
近年来，针对盲去模糊严重不适定性的问题，国内外学者在基于最大后验概率框架下不断提出新的方法。然而在该框架下提出的去模糊方法中，往往存在两个问题。首先是不能选取有效的图像先验正则项，这会导致在模糊核估计过程中，降低中间图像的估计质量，进而导致模糊核估计不准。为了求得模糊核，通常采用由粗尺度到精尺度的方法，在每个尺度下的图像与模糊核交替迭代更新，最终得到模糊核，这一过程中的图像称作中间图像。因此中间图像的求解直接影响到模糊核估计的质量。其次是选取的正则项可能将求解过程变为非凸稀疏优化的问题，导致难以求得全局最优解，因此在这种情况下选用何种方法来求解也变成了一个难题。
本发明针对中间图像估计不准的问题，引入了超拉普拉斯先验来正则化自然图像梯度的重尾分布特性，并通过引入广义软阈值算子解决了非凸稀疏求解最优化的问题，最终提出一种基于超拉普拉斯先验的图像盲去模糊方法，该方法能够有效地估计模糊核，提高模糊图像复原的质量。
技术实现要素：
：
本发明的主要目的是提出一种基于超拉普拉斯先验的图像盲去模糊方法，利用超拉普拉斯先验模拟自然图像梯度的重尾分布特性，有效求解中间图像，进而有利于估计模糊核，最终提高对模糊图像复原的质量。
为了实现上述目的，本发明提供如下技术方案：
步骤一、读入单幅模糊图像B，判断B是否为灰度图像，若不是灰度图像则对B进行灰度化处理；
步骤二、对B构建图像金字塔n为图像金字塔的总层数，i为图像金字塔的当前所在层，且i初始为1；将作为第一层初始的中间图像I1；
步骤三、由中间图像Ii求得显著性图像边缘
步骤四、由显著性图像边缘更新模糊核k，通过进行求最优解k，表示第i层金字塔图像的梯度，γ为权值；
步骤五、由模糊核k更新中间图像Ii，通过进行求最优解Ii，表示第i层中间图像的梯度，p表示超拉普拉斯分布函数的参数，λ为权值；
步骤六、判断第i层内模糊核与中间图像迭代更新的次数j是否达到预设最大迭代次数m，若达到则执行步骤七，否则迭代次数j加1且返回步骤三，j的取值为1到m；
步骤七、判断层数i是否达到金字塔第n层，若达到第n层，则执行步骤八，否则层数i加1并返回步骤三，且将第i+1层的中间图像Ii+1初始化为Ii；
步骤八、将第n层求得的k作为模糊图像B对应的最优模糊核K，选择合适的非盲去模糊算法复原出最终的图像。
与现有的技术相比，本发明具有以下有益效果：
(1)通过步骤五在求解中间图像时，本发明采用超拉普拉斯先验作为中间图像梯度的正则化约束项，相比较拉普拉斯或高斯等分布而言，超拉普拉斯分布可以更好地模拟中间图像梯度的重尾分布特性，因此能够有效地求解中间图像，进而有效求解模糊核，最终提高图像复原的质量。
(2)针对步骤五的中间图像求解方法，本发明使用了广义软阈值算子，很好地解决了因引入超拉普拉斯先验而造成的非凸稀疏最优值求解的问题，避免了局部最优的问题。
因此，本发明在摄影学、光学、天文学、医疗图像、监控、遥感探测和军事研究等领域将具有广泛的应用前景。
附图说明：
图1本发明的算法流程图；
图2构建模糊图像的金字塔示意图；
图3自然图像梯度分布及其模拟示意图；
图4图例的模糊核估计示意图；
图5图例的模糊核在估计过程中平方误差和的曲线图；
图6 Levin数据库中图像im01的模糊核估计结果图；
图7算法在Levin数据库的平均结构相似性指标分析图；
图8算法在Levin数据库的平均峰值信噪比指标分析图；
图9算法在Levin数据库的误差率曲线分析图；
图10算法对于自然模糊图像的复原效果图。
具体实施方式：
为了更好的说明本发明的目的、具体步骤以及特点，下面结合附图对本发明作进一步详细的说明：
参考图1，本发明提出的一种基于超拉普拉斯先验的图像盲去模糊方法，主要包含以下步骤：
步骤一、读入单幅模糊图像B，判断B是否为灰度图像，若不是灰度图像则对B进行灰度化处理；
步骤二、对B构建图像金字塔n为图像金字塔的总层数，i为图像金字塔的当前所在层，且i初始为1；将作为第一层初始的中间图像I1；
步骤三、由中间图像Ii求得显著性图像边缘
步骤四、由显著性图像边缘更新模糊核k，通过进行求最优解k，表示第i层金字塔图像的梯度，γ为权值；
步骤五、由模糊核k更新中间图像Ii，通过进行求最优解Ii，表示第i层中间图像的梯度，p表示超拉普拉斯分布函数的参数，λ为权值；
步骤六、判断第i层内模糊核与中间图像迭代更新的次数j是否达到预设最大迭代次数m，若达到则执行步骤七，否则迭代次数j加1且返回步骤三，j的取值为1到m；
步骤七、判断层数i是否达到金字塔第n层，若达到第n层，则执行步骤八，否则层数i加1并返回步骤三，且将第i+1层的中间图像Ii+1初始化为Ii；
步骤八、将第n层求得的k作为模糊图像B对应的最优模糊核K，选择合适的非盲去模糊算法复原出最终的图像。
上述技术方案中，步骤二中对模糊图像B构建n层图像金字塔，且n由给定模糊核k的尺寸决定。以Levin数据库[1]中图像im01-ker08为例说明构建金字塔的过程。该图像B的尺寸为255×255，模糊核k的尺寸为23×23。根据下列公式可计算出n：
其中k1和k2分别表示k的高和宽，floor表示向下取整操作符。这里k1＝k2＝23，代入公式求解得n等于5。最终构建如图2所示的五层分辨率递增的图像金字塔，分辨率分别为63×63，89×89，127×127，179×179和255×255，且每一层对应的模糊核尺寸分别为7×7，9×9，13×13，17×17和23×23。
上述技术方案中，步骤三中通过对中间图像Ii先求解显著性结构然后求得显著性边缘，具体求解步骤包括：
(1)首先从中间图像Ii中提取主要结构Is：
其中x表示像素点；Ii(x)表示x处的像素值，Is(x)表示Ii(x)的结构部分；Ii(x)-Is(x)表示Ii(x)的细节部分；表示水平和垂直方向的梯度微分算子；θ变量初值设定为1，且每当上式执行1次之后θ按照θ/1.1更新，即随着中间图像和模糊核之间迭代次数的递增θ值逐渐减小，使得图像细节部分的权值逐渐增加，因此保证之后的运算过程包含更多的细节部分。θ变量由金字塔各层模糊核求解过程共享，因此仅需要在金字塔第一层的第一次迭代过程中初始化一次即可。ω(x)为自适应变量，且
ω(x)＝exp(-|r(x)|0.8)，
其中
Nh(x)是以中心像素点x的h×h窗口，选择h等于5，y表示窗口区域内的像素点，表示窗口区域内像素点梯度的矢量和，表示窗口区域内像素点梯度的绝对值之和。在图像平坦的区域梯度值一般比较接近，所以矢量和较小；而在结构显著的区域梯度值一般相差较大，所以矢量和较大。因此相对而言，r在图像平坦的区域取较小值，而在图像结构显著的区域取较大值。可以通过自适应调整权值的方式，实现对图像结构的有效选取，即对估计模糊核而言提高有效的结构并减少对无效结构的选取。
(2)其次通过冲击滤波器求得显著性结构
其中Δ表示拉普拉斯算子，表示梯度算子，sign表示符号函数。ΔI＝Ix2Ixx+2IxIyIxy+Iy2Iyy，其中Ix和Iy分别表示水平与垂直方向的一阶微分，Ixx，Ixy和Iyy表示二阶微分。
(3)然后通过显著性结构求解显著性边缘
表示单位二值化掩膜函数，定义如下：
其中表示梯度幅值，t表示梯度幅值的阈值，初值设定为0.5，且每当内循环执行1次之后t按照t/1.1的方式更新，即阈值逐渐减小，以使得之后的运算可以包含更多的结构。t变量由金字塔各层模糊核求解过程共享，因此仅需要在金字塔第一层的第一次迭代过程中初始化一次即可。
上述技术方案中，步骤四中通过显著性边缘求解模糊核k，求解方法为：
其中γ为权值，γ＝10-2，上式可用迭代加权最小二乘法[3]求解。
上述技术方案中，步骤五中通过步骤四求得的模糊核k来求解中间图像Ii。如图3所示，自然图像的梯度分布呈重尾分布，由图可见超拉普拉斯分布相对于高斯分布和拉普拉斯分布能更好的模拟自然图像的梯度分布情况。因此可以更有效地求解中间图像Ii，进而在与模糊核迭代求解的过程中，可以求得更佳的模糊核。中间图像Ii可通过下式求解：
其中p通常取值范围是0.5到0.8之间，这里根据经验选择为0.8，λ为权值，λ＝4×10-3，对于该式求解可以针对梯度引入辅助变量G＝(Gh,Gv)T[4]，则上式转化为：
其中β为变化的权值，初始值定为2λ。对上式分别迭代求解Ii和G，满足迭代终止条件后，最终得到Ii，具体迭代步骤如下：
(1)首先根据下式，先固定Ii，求解最优G：
这里求解G的过程是求非凸稀疏最优问题，通过引入广义软阈值算子(GST)[5]求解：
其中表示G的最小值，表示阈值，sgn表示符号函数。
(2)然后根据下式，固定G，求解最优Ii：
可用最小二乘的方法求解：
其中F(·)和F-1(·)分别表示快速傅立叶变换和反变换；表示快速傅立叶变换的复共轭运算；和分别表示Ii的水平和垂直梯度。若β<βmax则继续迭代，且将β变为2β。直到β≥βmax停止迭代运算，得到最终的中间图像Ii，这里βmax设定为105。
上述技术方案中，步骤六用于判断金字塔第i层内，模糊核与中间图像迭代更新的次数j是否达到预设最大迭代次数m，若达到则执行步骤七，否则j加1返回步骤三。通过大量的实验权衡了准确率和时间两方面因素，最终设定迭代次数m为5。
上述技术方案中，步骤八将第n层求得的k作为模糊图像B对应的最优模糊核K，选择合适的非盲去模糊算法复原出最终的图像。本发明选用全变差图像重构的方法[6]来得到最终的结果，最终复原的图像X可以通过下式求解：
表示图像X的梯度，λf为权值，λf＝10-3。
上述技术方案中，图4给出了算法对Levin数据库im01-ker08实例的模糊核估计过程的示意图。图5给出了该实例的模糊核在估计过程中平方误差和的曲线图，通过本发明估计出的模糊核与真实模糊核的误差和进行比较，可以发现随着迭代次数的增加误差逐渐减少，直至近似收敛。图6给出了Levin数据库中图(a)的模糊核估计结果对比图，(b)为原始的模糊核，(c)为本发明估计出的模糊核，从定性角度可见本发明估计出的模糊核与原始模糊核基本一致。图7给出了算法在Levin数据库中的平均结构相似性指标图，结构相似性指标越接近于是1表示复原的图像在结构方面与原始清晰图像相似度越高。图8给出了算法在Levin数据库中的平均峰值信噪比图，峰值性噪比越大，表示复原的质量越好。图9给出了算法在Levin数据库中的误差率与成功率曲线图，通常认为误差率在小于等于2的区域时图像在视觉上是可以接受的，可见在数据库所有图像中，本发明在此区域已经占到了接近90％。图10给出了算法在自然模糊图像实例上的实验结果，可见模糊核轨迹清晰，且复原图像结果从视觉上看清晰度显著提高，图像细节还原度好。通过实验测试，实验结果以定量和定性进行说明，验证算法的有效性。
本专利运用超拉普拉斯先验来模拟自然图像的梯度分布，并结合提取图像显著性结构的方法来得到显著性边缘，有效地求解模糊核。本专利在利用模糊核求解中间图像的过程中，采用了广义软阈值算子有效地解决了非凸稀疏优化的问题。通过提高中间图像的估计效果，来增加模糊核估计的准确度，最终提高图像复原的质量。
上面结合附图对本发明的具体实施方式做了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化。
[1]http://www.wisdom.weizmann.ac.il/～levina/papers/LevinEtalCVPR09Data.rar
[2]Jinshan Pan,Risheng Liu,Zhixun Su,and Xianfeng Gu,Kernel estimation from salient structure for robust motion deblurring.IEEE Signal Processing Letters,vol.20,pp.841-844,2013.
[3]Anat Levin,Rob Fergus,Fredo Durand,and Bill Freeman,Image and depth from a conventional camera with a coded aperture.ACM Transactions on Graphics,vol.26,no.3,pp.70-78,2007.
[4]Jinshan Pan,Zhe Hu,Zhixun Su,and Ming-Hsuan Yang,Deblurring text images via L0-regularized intensity and gradient prior.IEEE Conference on Computer Vision and Pattern Recognition,pp.2901-2908,2014.
[5]Wangmeng Zuo,Deyu Meng,Lei Zhang,Xiangchu Feng,and David Zhang,A gen-eralized iterated shrinkage algorithm for non-convex sparse coding.IEEE International Conference on Computer Vision,pp.217-224,2013.
[6]Yilun Wang,Junfeng Yang,Wotao Yin,and Yin Zhang,A new alternating minimi-zation algorithm for total variation image reconstruction.SIAM J.Image,vol.1,no.3,pp.248-272,2008.