发明内容
本发明的目的在于解决现有技术方法在理论和实现两个方面的不足,提出一种基于L0稀疏先验的最大后验模糊核估计新方法。
本发明技术方案是:首先,引入基于L0范数的显著边缘稀疏先验,利用迭代硬阈值收缩法实现显著边缘特征的隐性自动预测;其次,引入基于L0范数的相机抖动模糊核稀疏先验,利用迭代重新加权最小二乘法实现模糊核的快速估计;最后,利用基于超拉普拉斯先验的图像非盲去模糊方法,实现高质量图像去模糊。
1、本发明具体实施步骤:
(1)利用水平方向和垂直方向的一阶导数算子▽h=[1,-1;0,0],▽v=[1,0;-1,0],得到相机抖动模糊图像y的梯度图像
(2)令待估计模糊核k的尺寸为Z×Z,为提高模糊核估计方法的收敛性,采用多尺度实现方式迭代估计模糊核;
(3)令初始模糊核k(1)的尺寸为3×3,且k(1)=[0,0,0;1,1,1;0,0,0]3,根据(3.1)-(3.12)确定尺度总数scales_num及各尺度下模糊核k(s)的尺寸ksize(s)×ksize(s):
(3.1) s=1;
(3.2) tmp=3;
(3.3) while(tmp<Z)
(3.4) ksize(s)=tmp;
(3.5) s=s+1;
(3.6) tmp=ceil(tmp*sqrt(2));
(3.7) if(mod(tmp,2)==0)
(3.8) tmp=tmp+1;
(3.9) end;
(3.10) end;
(3.11) ksize(s)=Z;
(3.12) scales_num=s;
(4)相应于各尺度下的模糊核k
(s),根据(4.1)-(4.3)确定模糊梯度图像y
d在各尺度下的尺寸r
(s)×c
(s)以及对应的插值模糊图像
(4.1) r(s)=floor(size(yd,1)*ksize(s)/Z);
(4.2) c(s)=floor(size(yd,2)*ksize(s)/Z);
(4.3)
(5)令显著边缘图像与模糊核在各尺度下的循环迭代次数为out_iter_num,循环迭代初始次数为out_iter=1,初始尺度为s=1,显著边缘L0稀疏先验的正则化参数为η,模糊核L0稀疏先验的正则化参数为γ,K
old对应模糊核k
(1)的二维矩阵,
对应模糊插值图像
的一维向量,k
old对应模糊核k
(1)的一维向量;
(6)令
τ
old=1,利用迭代硬阈值收缩法估计第s个尺度下的显著边缘
主要利用步骤(6.1)-(6.4)进行m次循环估计:
(6.1)
(6.2)
(6.3)τnew=τold2;
(6.4)
其中,硬阈值算子ΦHARD(·,·)定义如下:
(7)令权矩阵
利用迭代重新加权最小二乘方法估计第s个尺度下的模糊核k
(s)=k
new,具体而言利用步骤(7.1)-(7.4)进行n次循环估计:
(7.1)估计
(7.2)将knew投影到约束集
(7.3)更新
(7.4)kold=knew;
(8)更新out_iter:out_iter=out_iter+1;
(9)如果out_iter<out_iter_num,转到(6),否则转到(10);
(10)更新第s+1个尺度下模糊核k的初始值:k(s+1)=imresize(k(s),[ksize(s),ksize(s)],'bilinear');
(11)更新s:s=s+1;
(12)如果s<scales_num,返回(6),否则转到(13);
(14)利用基于超拉普拉斯先验的图像非盲去模糊方法,最终获得去模糊图像
有益效果:
(1)发明方法的模糊核估计是个严格意义上的稀疏最优化问题;
(2)发明方法的实现方式简单,无需平滑滤波、冲击滤波等任何迭代预处理;
(3)发明方法的计算复杂度低;
(4)发明方法的模糊核估计准确度高;
估计出相机抖动模糊核之后,利用基于超拉普拉斯先验的图像非盲去模糊方法,从而获得更高质量的去模糊图像。
2、本发明模型推导过程:
不失一般性,相机抖动模糊可利用如下卷积型观察模型进行表述
其中,y是获取的相机抖动模糊图像,f是真实清晰图像,k是空间不变的相机抖动模糊核,n是服从高斯分布的加性随机噪声,
代表卷积;为了表述的方便,该观察模型还可写成如下矩阵-向量表达形式
y=Kf+n
其中,y,f,以及n分别是y,f,以及n的向量表达形式,K是模糊核k的矩阵表达形式。采取分而治之的策略解决相机抖动去模糊问题,分成两大步骤:(1)模糊核估计;(2)非盲图像恢复;
模糊核估计是在图像梯度域实现的,为此,建立如下梯度域观察模型:
其中, 类似的,上述梯度域卷积型观察模型可以写成如下矩阵-向量表达形式
yd=Kfd+nd=Fdk+nd
其中,yd,fd,nd,k分别是yd,fd,nd,k的向量表达形式,Fd是fd的矩阵表达形式。
图像中的显著边缘是精确估计模糊核的重要所在。为了便于自动预测显著边缘,提高模糊核估计的准确性,首先提出基于L0范数的显著边缘先验及其优化模型:
其中,
是对应f
d的最优化结果,η是正则化参数,Ξ({f
d}
d∈Λ)是基于L0范数的显著边缘先验,定义为:
根据迭代硬阈值收缩方法,上述最优化问题可以下步骤(1)-(3)进行迭代求解:
(1)
其中,K
old对应之前估计的模糊核,
对应之前估计的显著边缘,
对应当前更新的显著边缘,Φ
HARD(·,·)代表硬阈值算子,定义如下:
(2)τnew=τold2,τold对应当前的调整参数,τnew代表更新的调整参数;
(3)τold=τnew;
当显著边缘给定时,利用相机抖动模糊核的稀疏特性,从而能够建立基于L0范数的模糊核先验及其优化模型:
其中,Ξ(k)=k0是基于L0范数的模糊核先验,kest是对应k的最优化结果,γ是正则化参数;
根据迭代重新加权最小二乘方法,首先上述将最优化问题转化为如下近似形式:
其中,W=diag({1(kl+0.0001)2},然后利用以下步骤(4)-(6)进行迭代求解:
(4)利用共轭梯度法估计
其中,W为之前估计的权矩阵
为之前估计的显著边缘的矩阵表达形式,k
new为当前更新的模糊核;
(5)进一步将knew投影到约束集
(6)更新
利用上述方法原理,采用多尺度实现方式迭代估计显著边缘图像与相机抖动模糊核,从而得到最终的估计模糊核
。
具体实施方式
(1)利用水平方向和垂直方向的一阶导数算子▽h=[1,-1;0,0],▽v=[1,0;-1,0],利用MATLAB函数conv2获得相机抖动模糊图像y的梯度图像
yh=conv2(y,▽h,'valid');
yv=conv2(y,▽v,'valid');
(2)令待估计模糊核k的尺寸为Z×Z,为提高模糊核估计方法的收敛性,采用多尺度实现方式迭代估计模糊核;
(3)令初始模糊核k(1)=[0,0,0;1,1,1;0,0,0]3(尺寸为3×3),并且利用以下MATLAB代码确定尺度总
数scales_num及各尺度下模糊核k(s)的尺寸ksize(s)×ksize(s):
s=1;
tmp=3;
while(tmp<Z)
ksize(s)=tmp;
s=s+1;
tmp=ceil(tmp*sqrt(2));
if(mod(tmp,2)==0)
tmp=tmp+1;
end;
end;
ksize(s)=Z;
scales_num=s;
(4)相应于各尺度下的模糊核k
(s),利用以下MATLAB代码确定模糊梯度图像y
d在各尺度下的尺寸r
(s)×c
(s)以及对应的插值模糊图像
r(s)=floor(size(yd,1)*ksize(s)/Z);
c(s)=floor(size(yd,2)*ksize(s)/Z);
(5)令显著边缘图像与模糊核在各尺度下的循环迭代次数为out_iter_num,循环迭代初始次数为out_iter=1,初始尺度为s=1,显著边缘L0稀疏先验的正则化参数为η,模糊核L0稀疏先验的正则化参数为γ,K
old对应模糊核k(1)的二维矩阵,k
old对应模糊核k
(1)的一维向量,
对应模糊插值图像
的一维向量;
(6)令
τ
old=1,利用迭代硬阈值收缩法估计第s个尺度下的显著边缘
(迭代硬阈值收缩法具体参见文献《Iterative hard thresholding for compressed sensing》,(T.Blumensath,M.E.Davies),Appliedand Computational Harmonic Analysis,2009,vol.27,no.3,pp.265–274)。具体利用以下MATLAB代码进行m次循环估计:
τnew=τold/2;
τold=τnew;
(7)令权矩阵
利用迭代重新加权最小二乘方法估计第s个尺度下的模糊核k
(s)=k
new。迭代重新加权最小二乘方法可参见文献《Iteratively reweighted least squares minimization forsparse recovery》(I.Daubechies,et al.),Communications on Pure and Applied Mathematics,2010,vol.63,pp.1–38。具体利用以下步骤进行n次循环估计:
(7.1)利用MATLAB函数pcg(共轭梯度法)求解如下方程组,估计knew:
共轭梯度法可参见参考文献《Convex optimization》(S.Boyd,L.Vandenberghe),Cambridge University Press,2004。
(7.2)将knew投影到约束集
主要利用以下MATLAB代码实现:
knew(knew<0)=0;
sumk=sum(knew(:));
knew=knew./sumk;
(7.3)更新
(7.4)kold=knew;
其中,
对应显著边缘
的二维矩阵,k
new对应k
new的矩阵形式;
(8)更新out_iter:out_iter=out_iter+1;
(9)如果out_iter<out_iter_num,转到(6),否则转到(10);
(10)更新第s+1个尺度下模糊核k的初始值:k(s+1)=imresize(k(s),[ksize(s),ksize(s)],'bilinear');
(11)更新s:s=s+1;
(12)如果s<scales_num,返回(6),否则转到(13);
(14)利用基于超拉普拉斯先验的图像非盲去模糊方法恢复出去模糊图像
具体参见文献《Fast imagedeconvolution using hyper-laplacian priors》(D.Krishnan and R.Fergus),Advances in Neural InformationProcessing Systems,2009,vol.22,pp.1033-1041.
(15)为了验证基于L0稀疏先验的盲相机抖动去模糊方法的有效性,设计4组仿真实验,每组实验均提供了真实图像及真实模糊核(附图2、附图3、附图4、附图5),并与最大似然估计盲去模糊方法(简写为MaxL)以及基于规范化稀疏度量的盲去模糊方法(简写为NormSpar)进行比较,MaxL可调用MATLAB图像处理工具箱中的函数deconvblind进行实现,NormSpar具体参见文献《Blind deconvolution using a normalizedsparsity measure》(D.Krishnan,et al.),IEEE Conf.Computer Vision and Pattern Recognition,2011:233–240。为了衡量不同方法的去模糊效果,采用峰值信噪比(简称PSNR)和结构相似性度量(简称SSIM)这两个客观标准进行度量。所有方法的参数均采用人工方式选取,以取得最佳PSNR、SSIM、以及视觉效果为准。此外,具体的实验环境包括:软件MATLAB版本号为v7.0,电脑的配置为Pentium(R)Core-Duo T4200CPU(2GHZ),2GB RAM,Microsoft Windows XP(version2002,Service Pack3).
(16)PSNR的定义:
f为真实图像,
为去模糊图像,M为f的像素个数。
(17)SSIM的具体定义与实现参见文献《Image quality assessment:from error measurement to structuralsimilarity》(Z.Wang,et al.),IEEE Transactios on Image Processing,2004,vol.13,no.4,pp.600–612.
(18)表1提供了上述3种盲去模糊方法在4组仿真实验中的PSNR和SSIM。由表1可见,本发明在所有中均取得了最高PSNR,特别是对于图像Lena、House、Boat,本发明的优势明显;从SSIM的角度看,本发明在对应图像Lena、House、Boat的3组实验中的去模糊效果优于MaxL和NormS,对图像Cameraman本发明取得了与NormS基本相似的去模糊效果,且两者均优于MaxL。
表1.各种盲去模糊方法在4组仿真实验中的PSNR和SSIM
(19)从视觉感知的角度看,本发明也取得了最优的去模糊效果。附图6给出了对应图像Lena的相机抖动模糊图;附图7给出了基于规范化稀疏度量盲去模糊方法的恢复图像及估计模糊核;附图8给出了本发明基于L0稀疏先验的盲去模糊方法的恢复图像及估计模糊核。容易看出,本发明估计的模糊核更接近真实模糊核,且恢复的图像具有更清晰的边缘和较少的振铃效应;而基于规范化稀疏度量盲去模糊方法的恢复图像与估计模糊核与真实情况相去甚远。此外,附图9给出了基于L0稀疏先验的盲去模糊方法的显著边缘图像。
以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应了解,本发明不受上述实施例的限制,上述实施例和说明书中描述的只是说明本发明的原理,在不脱离本发明精神和范围的前提下,本发明还会有各种变化和改进,这些变化和改进都落入要求保护的本发明范围内,本发明要求保护范围由所附的权利要求书其等效物界定。