CN103761710B - 基于边缘自适应的高效图像盲去模糊方法 - Google Patents

Abstract

本发明公开了一种基于边缘自适应的图像盲去模糊方法，针对现有全变差去模糊算法易模糊图像边缘和细节的问题，构造了去均值的梯度全变差正则模型，并利用图像梯度的局部方差自适应计算迭代加权系数，有效提升了去模糊算法恢复图像边缘和细节的能力。其实现步骤是：（1）输入模糊图像,对梯度域清晰图像及模糊核交替求解，得到模糊图像的初始模糊核；（2）使用初始模糊核对模糊图像进行一次非盲去模糊，得到初始清晰图像；（3）对初始清晰图像聚类，更新去均值正则模型中的均值和加权系数并重新求解模糊核；（4）使用新的模糊核进行二次非盲去模糊，得到清晰图像。实验结果表明，本发明比现有技术具有更好的去模糊效果，可用于图像恢复。

Description

基于边缘自适应的高效图像盲去模糊方法

技术领域

本发明涉及图像处理技术领域，具体涉及图像的盲去模糊方法，可用于图像恢复。

背景技术

在相机成像过程中，由于相机抖动、离焦、散焦，或者目标物体运动，人们拍摄的照片经常会存在模糊。从单幅模糊图像恢复出清晰的数字图像在数字媒体娱乐、国防安全等方面均具有重要的需求。通常情况下，相机相对于景物的运动参数等信息是未知的，为了去除图像的模糊，需要同时估计图像模糊核和清晰的数字图像，在模糊核未知情况下对图像去模糊即是图像盲去模糊，图像盲去模糊是一个较为困难的图像恢复问题。

图像盲去模糊问题可以表示为其中y表示模糊图像，k表示图像模糊核，x表示原始清晰图像，表示卷积算子，n表示加性高斯噪声。图像盲去模糊算法通常包含两个交替迭代的步骤：1)图像模糊核k估计，即利用初始恢复清晰图像估计得到图像的模糊核；2)清晰图像x估计，利用估计得到的模糊核进行非盲图像去模糊，得到去模糊图像x。这两个过程依次交替，最终得到清晰的图像。

图像盲去模糊是一个欠定的图像逆问题。为了获得清晰的数字图像，需要利用图像的先验知识，如何挖掘、利用准确的图像先验知识是图像盲去模糊的关键因素。传统的图像去模糊算法广泛使用图像的光滑性，即用L₂范数对图像的梯度进行正则约束，该方法虽能有效去除噪声，但会丢失图像的细节。为了保持图像细节，研究人员也提出用L_p，0.7≤p≤1范数对图像的梯度进行正则化，相比L₂，L_p范数能更好的保持了图像的边缘结构和局部细节。此外，人们也提出利用学习的方法来获得自适应的滤波器以代替梯度算子，从而能更好的保持图像的边缘和纹理结构。这些改进的图像先验模型，在图像模糊核已知的情况下能较好的恢复出清晰图像。但是，对于图像盲去模糊，利用上述现有的图像先验模型，通常难以得到一个清晰的图像。这是因为，上述先验图像模型对于模糊图像能获得更小的正则项能量，从而在目标函数能量最小化求解过程中得到原始模糊图像的局部极值。为了获得清晰图像，通常需要采用一些复杂的方法来对图像边缘进行锐化，但这种方法不仅复杂，而且从理论和效果上也难以保证其可行性。

发明内容

本发明的目的在于针对上述现有技术的不足，提出一种基于边缘自适应的高效图像盲去模糊方法，以提高恢复图像的清晰度，使其最大程度的逼近原始清晰图。

实现本发明目的技术方案是：构造新的梯度加权正则模型，利用梯度的方差的倒数作为加权系数，使得图像先验模型对清晰自然图像具有更小的能量，而对模糊图像具有较大的能量，从而避免原始模糊图像成为目标函数的解；通过构造去均值的梯度正则模型，提升图像边缘和细节。具体步骤包括如下：

(1)输入模糊图像y，其模糊核为k，待求解的清晰图像为x，设置清晰图像x的初始解为模糊图像y，设待求解清晰图像的梯度域图像为Dx；

(2)初始化待求解图像的梯度域图像Dx的均值μc为0，初始化非盲去模糊次数f＝0；

(3)用迭代优化求解算法更新模糊核k：

(3a)设置迭代次数L，初始化迭代序号l为1，初始化迭代前的模糊核初始化迭代前的待求解图像的梯度域图像

(3b)按下式优化求解待求解图像的梯度域图像：

$D{\hat{x}}^{\left(l\right)}=\arg \underset{Dx}{min}\left|\right|Dy-{\hat{k}}^{(l-1)}\⊗D{\hat{x}}^{(l-1)}|{|}_2^2+\underset{c=1}{\overset{X}{\Σ}}\underset{i\∈S_c}{\Σ}{\mathrm{\λ}}_c^{(l-1)}||{{\mathrm{Dx}}_i}^{(l-1)}-{\mathrm{\μ}}_c|{|}_1,$

其中，表示第l次解优化式后获得的待求解清晰图像的梯度域图像，表示目标函数取最小值时Dx所取的值，Dy表示梯度域模糊图像，表示第l-1次迭代后模糊核的取值，表示第l-1次迭代后的待求解清晰图像的梯度域图像，表示卷积运算，表示2范数的平方，X表示像素分类数，c表示当前像素的像素分类序号，1≤c≤X，表示对各类的正则项求和，S_c表示像素类别序号是c的清晰图像x中的像素集合，表示对像素属于集合S_c的正则项求和，表示本次迭代的正则项加权系数，其值为梯度的方差的倒数。||Dx_i ^(l-1)-μ_c||₁表示去均值的梯度正则模型，Dx_i ^(l-1)表示第l-1次迭代后的待求解清晰图像的梯度域图像的第i个像素点值，||·||₁表示1范数；

(3c)根据待求解图像的梯度域图像，按下式优化求解模糊核：

${\hat{k}}^{\left(l\right)}=\arg \underset{k}{min}\left|\right|Dy-{\hat{k}}^{(l-1)}\⊗D{\hat{x}}^{\left(l\right)}|{|}_2^2+\mathrm{\η}|\left|{\hat{k}}^{(l-1)}\right|{|}_1,s.t.\left\{\begin{array}{c}\left|\right|{\hat{k}}^{\left(l\right)}|{|}_1=1\\ {\hat{k}}^{\left(l\right)}\≥0\end{array}\right.,$

其中，表示第l次解优化式后获得的模糊核，表示目标函数取最小值时模糊核k所取的值，η表示正则系数，表示中的所有元素都大于0；

(3d)令l＝l+1，重复步骤(3b)～(3c)共L次，得到模糊核

(4)根据计算的模糊核k非盲去模糊，按下式优化求解待求解的清晰图像:

$x^{\left(f\right)}=\arg \underset{x}{min}\underset{i}{\Σ}\left|\right|W^{\left(i\right)}|{|}_2+\frac{\mathrm{\β}}{2}\underset{i}{\Σ}||w^{\left(i\right)}-(D^{\left(i\right)}x-{\mathrm{\μ}}_c)|{|}_2^2+\frac{u}{2}|\[y-k\⊗x|{|}_2^2,$

其中x^(f)表示解优化式后更新的待求解清晰图像，表示目标函数取最小值时x所取的值，i取值为1和2，w⁽ⁱ⁾表示逼近正则项变量，||·||₂表示2范数，D⁽ⁱ⁾表示方向梯度算子，D⁽ⁱ⁾＝[D⁽¹⁾,D⁽²⁾]，D⁽¹⁾＝[-11],β表示惩罚因子，u为范围在0到1之间的数据项系数；

(5)更新非盲去模糊次数f＝1，对更新后待求解的清晰图像x⁽⁰⁾进行K均值算法聚类，按像素分类计算待求解的清晰图像x⁽⁰⁾梯度域像素图像的均值μ_c'，用所得均值更新步骤(3b)公式中待求解图像的梯度域图像的均值，即μ_c＝μ_c'；

(6)重复步骤(3)～(4)，得到清晰图像x＝x⁽¹⁾。

本发明与现有技术相比具有以下优点：

第一，本发明由于构造了梯度域图像优化式，且在该优化式中设置了去均值的梯度正则模型，使得正则项更稀疏，极大地提升了清晰图像的边缘和细节。

第二，本发明由于在梯度域图像优化式中，使用正则项加权系数对去均值的梯度正则模型进行梯度加权，利用梯度的方差的倒数作为正则项加权系数，使得优化式对清晰自然图像具有更小的能量，而对模糊图像具有较大的能量，从而避免原始模糊图像成为目标函数的解，克服了传统先验正则模型偏向选择模糊图像的问题。

附图说明

图1为本发明的实现流程图；

图2为待处理的模糊图；

图3为用本发明和现有四种方法对图2进行去模糊后的结果对比图。

具体实施方式

下面结合附图对本发明做进一步详细描述。

参照图1，本发明的实施步骤如下：

步骤1：输入模糊图像，设置待求解清晰图像的梯度域图像。

输入一副模糊图像y，该模糊图像的模糊核为k；设待求解的清晰图像为x，待求解清晰图像x的梯度域图像为Dx。

步骤2：初始化参数：

初始化待求解图像的梯度域图像Dx的均值μ_c为0；初始化去模糊次数f＝0；将清晰图像x的初始解设为模糊图像y。

步骤3：更新模糊核。

3a)设置迭代次数L，初始化迭代序号l为1，初始化迭代前的模糊核为值全1的矩阵，该矩阵维数根据具体模糊图像设置，初始化迭代前的待求解图像的梯度域图像式中D⁽ⁱ⁾＝[D⁽¹⁾,D⁽²⁾]，D⁽¹⁾＝[-11], 表示卷积运算；

3b)计算加性高斯噪声标准差和当前像素梯度的标准差

$\begin{array}{cc}{\mathrm{\σ}}_n^{\left(l-1\right)}=\sqrt{\frac{1}{M\×N-1}\underset{i=1}{\overset{M\×N}{\Σ}}{\left({\mathrm{Dy}}_i-{\hat{k}}^{\left(l-1\right)}\⊗{{\mathrm{Dx}}_i}^{\left(l-1\right)}\right)}^2}& {\mathrm{\σ}}_c^{\left(l-1\right)}=\sqrt{\frac{1}{M\×N-1}\underset{i=1}{\overset{M\×N}{\Σ}}{\left({{\mathrm{Dx}}_i}^{\left(l-1\right)}-{\mathrm{\μ}}_c\right)}^2}\end{array}$

其中M×N表示清晰图像的大小，Dy_i表示模糊图像y的梯度域图像第i个像素点值，Dx_i ^(l-1)表示第l-1次迭代后待求解清晰图像x的梯度域图像的第i个像素点值，表示第l-1次迭代后模糊核的取值；

3c)计算正则项加权系数 其中λ为0到1之间的调整参数；

3d)按下式优化求解待求解清晰图像x的梯度域图像：

$D{\hat{x}}^{\left(l\right)}=\arg \underset{Dx}{min}\left|\right|Dy-{\hat{k}}^{(l-1)}\⊗D{\hat{x}}^{(l-1)}|{|}_2^2+\underset{c=1}{\overset{X}{\Σ}}\underset{i\∈S_c}{\Σ}{\mathrm{\λ}}_c^{(l-1)}||{{\mathrm{Dx}}_i}^{(l-1)}-{\mathrm{\μ}}_c|{|}_1,$

其中，表示第l次解优化式后获得的待求解清晰图像x的梯度域图像，表示目标函数取最小值时Dx所取的值，Dy表示模糊图像的梯度域图像，表示2范数的平方，X表示清晰图像x的像素分类数，c表示当前像素的像素分类序号，1≤c≤X，表示对各类的正则项求和，S_c表示像素类别序号是c的清晰图像x中的像素集合，表示对像素属于集合S_c的正则项求和，||Dx_i ^(l-1)-μ_c||₁表示去均值的梯度正则模型，||·||₁表示1范数；

3e)根据待求解图像的梯度域图像，按如下优化式求解第l次解优化式后获得的模糊核

${\hat{k}}^{\left(l\right)}=\arg \underset{k}{min}\left|\right|Dy-{\hat{k}}^{(1-1)}\⊗D{\hat{x}}^{\left(l\right)}|{|}_2^2+\mathrm{\η}|\left|{\hat{k}}^{(l-1)}\right|{|}_1,s.t.\left\{\begin{array}{c}\left|\right|{\hat{k}}^{\left(l\right)}|{|}_1=1\\ {\hat{k}}^{\left(l\right)}\≥0\end{array}\right.,$

其中，表示目标函数取最小值时模糊核所取的值，η表示值为0到1之间的正则系数，s.t.表示模糊核的约束范围，表示中的所有元素都大于0；

3f)令l＝l+1，重复步骤(3b)～(3e)共L次，得到模糊图像y的模糊核

步骤4：求解待求解的清晰图像。

4a)计算逼近正则项变量w⁽ⁱ⁾：

w⁽ⁱ⁾＝max(|D⁽ⁱ⁾x-μ_c|-(λ/β),0)×sign(D⁽ⁱ⁾x-μ_c)+μ_c，

其中max(|D⁽ⁱ⁾x-μ_c|-(λ/β),0)表示对矩阵|D⁽ⁱ⁾x-μ_c|-(λ/β)的每个元素与0比较取最大值，β表示初值为15的惩罚因子，|·|表示取绝对值，sign(·)表示符号函数；

4b)根据上述步骤3f)得到的模糊图像y的模糊核k和步骤4a)得到的正则项变量w⁽ⁱ⁾，并更新β＝β×1.09，按下式优化求解待求解的清晰图像：

$x^{\left(t\right)}=\arg \underset{x}{min}\underset{i}{\Σ}\left|\right|W^{\left(i\right)}|{|}_2+\frac{\mathrm{\β}}{2}\underset{i}{\Σ}||W^{\left(i\right)}-(D^{\left(i\right)}x-{\mathrm{\μ}}_c)|{|}_2^2+\frac{u}{2}|\[y-k\⊗x|{|}_2^2,$

其中x^(f)表示解优化式更新后的待求解清晰图像，表示目标函数取最小值时x所取的值，u为范围在0到1之间的数据项系数。

4c)w⁽ⁱ⁾和x^(f)交替优化求解V次，得到更新后的待求解清晰图像。V值根据模糊图像的模糊程度设置。

步骤5：更新非盲去模糊次数为f＝1，求解初次更新的待求解图像的梯度域图像的均值μ_c'。

计算初次更新的待求解清晰图像x⁽⁰⁾梯度域图像的均值μ_c'，用所得均值更新步骤(3b)公式中待求解图像的梯度域图像的均值，即μ_c＝μ_c'，μ_c'的求解步骤如下：

5a)对初次更新的待求解清晰图像x⁽⁰⁾进行K均值聚类：

5a1)对初次更新的待求解清晰图像x⁽⁰⁾，以第i个像素为中心，1≤i≤M×N，取大小为n的块按每个块与周围块的能量大小，在块周围找到与其能量最相似的m个相似块，生成相似块矩阵索引集合G_i＝[j_i1,j_i2...,j_is,...,j_im]^T和权重矩阵W_ij＝[w_i1,w_i2...,w_it,...,w_ij]，1≤s≤m，1≤t≤j，1≤j≤M×N；

5a2)将每个像素点的相似块矩阵索引集合G_i按列排列，得到最终的相似块索引集合G，G是一个m行M×N列的矩阵，第i列表示与图像中第i个像素块相似的m个块的索引；

5a3)将每个像素点的权重矩阵W_ij按行排列，得到最终的权重矩阵W，权重矩阵W是M×N的方阵；

5a4)将相似块矩阵索引集合G中第i列的索引值记为z，将初次更新的待求解清晰图像x⁽⁰⁾中第i个像素块相似的m个块的对应权重存放在权重矩阵W的第i行第z列，第i行的其余M×N-m个位置存放0值；

5b)计算初次更新的待求解清晰图像x⁽⁰⁾的梯度域图像的均值。

设初次更新的待求解清晰图像x⁽⁰⁾的梯度域列向量是则x⁽⁰⁾的梯度均值图像的列向量将用matlab中的reshape函数还原成M×N的矩阵，得到初次更新的待求解清晰图像x⁽⁰⁾的梯度域图像均值μ_c'。

步骤6：重复步骤3～步骤4，得到清晰图像x＝x⁽¹⁾，x⁽¹⁾是第二次更新的待求解清晰图像，也是最终恢复的清晰图像。

下面结合仿真实验对本发明的去模糊效果做进一步的描述。

1.实验条件

本发明的实验运行系统为Intel(R)Core(TM)i7-2600CPU6503.40GHz，32位Windows操作系统，仿真软件采用MATLAB(R2012b)。实验使用的模糊图像来源于RolfKohler,MichaelHirsch,BettyMohler等人发表的文章“Recordingandplaybackofcamerashake:benchmarkingblinddeconvolutionwithareal-worlddatabase”中的标准数据库，如图2所示。

2.实验内容

用本发明和现有的四种去模糊方法对图2进行去模糊处理，结果如图3，其中

图3(a)是用Cho,S.,Lee等人在文章“Fastmotiondeblurring，”(In:ACMTransactionsonGraphicsSIGGRAPHASIA，2009)中公开的去模糊方法的对图2的去模糊结果；

图3(b)是用Xu,L.,Jia等人在文章“Two-phasekernelestimationforrobustmotiondeblurring，”(In:ProceedingsoftheEuropeanConferenceonComputerVision，ECCV，2010)中公开的去模糊方法对图2的去模糊结果；

图3(c)是用Shan,Q.,Jia,J.,Agarwala,A.等人在文章“High-qualitymotiondeblurringfromasingleimage，”(In:ACMTransactionsonGraphicsSIGGRAPH.2008)中公开的去模糊方法对图2的去模糊结果；

图3(d)是用Hirsch,M.,Schuler,C.J.等人在文章“Fastremovalofnon-uniformcamera-shake，”(In:ProceedingsoftheIEEEInternationalConferenceonComputerVision，ICCV，2011)中公开的去模糊方法对图2的去模糊结果；

图3(e)是用本发明方法对图2的去模糊结果。

3.仿真结果分析

仿真实验中，采用峰值信噪比PSNR指标来评价压缩感知实验结果，峰值信噪比PSNR定义为：

$PSNR=10{\log }_{10}\left(\frac{{255}^2}{\mathrm{\Σ}\left|\right|y-x|{|}_2^2}\right)$

其中，y为原图像，x为恢复出来的去模糊图像，峰值信噪比PSNR的值越大则去模糊性能越好。

本发明和现有的四种去模糊方法对图2进行去模糊处理所得去模糊图像的PSNR如下：

图3(a)的PSNR值为33.948；

图3(b)的PSNR值为33.236；

图3(c)的PSNR值为31.366；

图3(d)的PSNR值为33.161；

图3(e)的PSNR值为34.942；

从所得数据可以看出，本发明方法的去模糊结果的PSNR值高于其他方法的去模糊结果的PSNR值，即本发明比现有技术具有更好的去模糊效果。

Claims (3)

1.一种基于边缘自适应的高效图像盲去模糊方法，包括如下步骤：

(1)输入模糊图像y，其模糊核为k，待求解的清晰图像为x，设置清晰图像x的初始解为模糊图像y，设待求解清晰图像的梯度域图像为Dx；

(2)初始化待求解图像的梯度域图像Dx的均值μ_c为0，初始化去模糊次数f＝0；

(3)用迭代优化求解算法更新模糊核k：

(3a)设置迭代次数L，初始化迭代序号l为1，初始化迭代前的模糊核初始化迭代前的待求解图像的梯度域图像

(3b)按下式优化求解待求解图像的梯度域图像：

$D{\hat{x}}^{\left(l\right)}=\arg \underset{Dx}{min}\left|\right|Dy-{\hat{k}}^{(l-1)}\⊗D{\hat{x}}^{(l-1)}|{|}_2^2+\underset{c=1}{\overset{X}{\Σ}}\underset{i\∈S_c}{\Σ}{\mathrm{\λ}}_c^{(l-1)}||{\mathrm{Dx}}_i^{(l-1)}-{\mathrm{\μ}}_c|{|}_1,$

其中，表示第l次解优化式后获得的待求解清晰图像的梯度域图像，表示目标函数取最小值时Dx所取的值，Dy表示模糊图像的梯度域图像，表示第l-1次迭代后模糊核的取值，表示第l-1次迭代后的待求解清晰图像的梯度域图像，表示卷积运算，表示2范数的平方，X表示像素分类数，c表示当前像素的像素分类序号，1≤c≤X，表示对各类的正则项求和，S_c表示像素类别序号是c的清晰图像x中的像素集合，表示对像素属于集合S_c的正则项求和，表示本次迭代的梯度正则加权系数，其值为梯度的方差的倒数；表示去均值的梯度正则模型，Dx_i ^(l-1)表示第l-1次迭代后的待求解清晰图像的梯度域图像的第i个像素点值，||·||₁表示1范数；

(3c)根据待求解图像的梯度域图像，按下式优化求解模糊核：

${\hat{k}}^{\left(l\right)}=\arg \underset{k}{min}\left|\right|Dy-{\hat{k}}^{(l-1)}\⊗D{\hat{x}}^{\left(l\right)}|{|}_2^2+\mathrm{\η}|\left|{\hat{k}}^{(l-1)}\right|{|}_1,s.t.\left\{\begin{array}{c}\left|\right|{\hat{k}}^{\left(l\right)}|{|}_1=1\\ {\hat{k}}^{\left(l\right)}\≥0\end{array}\right.,$

其中，表示第l次解优化式后获得的模糊核，表示目标函数取最小值时模糊核k所取的值，η表示值为0到1之间的正则系数，表示中的所有元素都大于或等于0；

(3d)令l＝l+1，重复步骤(3b)～(3c)共L次，得到模糊核

(4)根据计算的模糊核k非盲去模糊，按下式优化求解待求解的清晰图像:

$x^{\left(f\right)}=\arg \underset{x}{min}\underset{i}{\Σ}\left|\right|w^{\left(i\right)}|{|}_2+\frac{\mathrm{\β}}{2}\underset{i}{\Σ}||w^{\left(i\right)}-(D^{\left(i\right)}x-{\mathrm{\μ}}_c)|{|}_2^2+\frac{u}{2}\left|\right|y-k\⊗x|{|}_2^2,$

其中x^(f)表示解优化式后更新的待求解清晰图像，表示目标函数取最小值时x所取的值，i取值为1和2，w⁽ⁱ⁾表示逼近正则项变量，||·||₂表示2范数，D⁽ⁱ⁾表示方向梯度算子，D⁽ⁱ⁾＝[D⁽¹⁾,D⁽²⁾]，D⁽¹⁾＝[-11],β表示惩罚因子，u为范围在0到1之间的数据项系数；

(5)更新非盲去模糊次数f＝1，对更新后待求解的清晰图像x⁽⁰⁾进行K均值算法聚类，按像素分类计算待求解的清晰图像x⁽⁰⁾梯度域像素图像的均值μ_c'，用所得均值更新步骤(3b)公式中待求解图像的梯度域图像的均值，即μ_c＝μ_c'；

(6)重复步骤(3)～(4)，得到清晰图像x＝x⁽¹⁾。

2.根据权利要求1所述的一种基于边缘自适应的高效图像盲去模糊方法，其特征在于，步骤(3b)中本次迭代的梯度正则加权系数通过下式计算：

其中λ为0到1之间的常数，表示加性高斯噪声方差，表示当前像素梯度的标准差。

3.根据权利要求1所述的一种基于边缘自适应的高效图像盲去模糊方法，其特征在于，步骤(5)中按像素分类计算待求解的清晰图像x⁽⁰⁾梯度域像素图像的均值μ_c'，按如下步骤进行：

5a)对初次更新的待求解清晰图像x⁽⁰⁾进行K均值聚类，形成X个像素类别；

5b)按下式计算初次更新的待求解清晰图像x⁽⁰⁾的梯度域图像的均值：

${{\mathrm{\μ}}_c}^{\′}=\frac{1}{\left|S_c\right|}\underset{i\∈S_c}{\Σ}{\mathrm{Dx}}_i^{\left(0\right)},$

其中c表示当前像素的像素分类序号，1≤c≤X，S_c表示像素类别序号为c的清晰图像x⁽⁰⁾中的像素集合，|S_c|表示初次更新的清晰图像x⁽⁰⁾中的像素类别属于S_c的像素个数，表示初次更新的待求解清晰图像x⁽⁰⁾的第i个像素点值，表示初次更新的待求解清晰图像x⁽⁰⁾的梯度域图像的第i个像素点值，1≤i≤M×N，M×N是初次更新的待求解清晰图像x⁽⁰⁾的大小。