基于L0稀疏先验的盲相机抖动去模糊方法

摘要

基于L0稀疏先验的盲相机抖动去模糊方法属于数字图像处理技术领域，是一种对相机抖动造成的模糊图像进行去模糊的方法，能够估计各种空间不变的相机抖动模糊核（点扩散函数）；该方法不仅避免了当前变分贝叶斯估计方法计算复杂度高的问题，而且解决了当前最大后验估计方法缺乏严格最优化理论支撑的问题；首先，引入基于L0范数的显著边缘稀疏先验，利用迭代硬阈值收缩法实现显著边缘特征的隐性自动预测；其次，引入基于L0范数的相机抖动模糊核稀疏先验，利用迭代重新加权最小二乘法实现模糊核的快速估计；最后，利用基于超拉普拉斯先验的图像非盲去模糊方法，获得高质量去模糊图像；该方法流程图如附图1所示。

说明书

技术领域

本发明属于数字图像处理技术领域，具体涉及对相机抖动造成的模糊图像进行去模糊的领域。

背景技术

盲相机抖动去模糊是指去除或减轻相机拍摄数字图像过程中发生的图像运动模糊现象，是近年来数字图像处理中非常重要又极具挑战性的研究内容，其核心是估计对应各种相机抖动的模糊核（点扩散函数）。当前，研究人员已经提出了多种不同方法。

最经典的盲相机抖动去模糊方法是的最大似然估计法，已集成到数值计算软件MATLAB中的图像处理工具箱，但是该方法在实际应用中具有明显的局限性，不仅要求模糊图像具有较高的信噪比，而且要求模糊核的尺寸较小。该方法的另一个缺点是，去模糊图像中的边缘等细节信息不能很好地恢复，常伴有明显的振铃效应。

目前，国际上提出了多种盲相机抖动去模糊的新方法。根据模糊核的估计准则，这些方法主要分为两大类：变分贝叶斯估计方法和最大后验估计方法。变分贝叶斯估计方法的计算复杂度较高，其中的两种代表性方法是：Fergus等人提出了基于高斯混合模型的变分贝叶斯盲相机抖动去模糊方法，参见文献《Removing camera shake from a single photograph》,ACM Trans.Graph.,2006,vol.25,no.3,pp.787–794；之后，Levin等人基于类似的图像建模思想降低了上述方法的计算复杂度，参见文献《Efficient marginallikelihood optimization in blind deconvolution》,IEEE Conf.Computer Vision and Pattern Recognition,2011,pp.2657–2664。相比变分贝叶斯估计方法，盲相机抖动去模糊最大后验估计方法的计算复杂度要低得多。代表性方法参见文献《Fast motion deblurring》,ACM Transactions on Graphics(SIGGRAPH ASIA),2009,vol.28,no.5,article no.145,《Psf estimation using sharp edge prediction》,IEEE Conf.Computer Vision and PatternRecognition,2008,pp.1–8。但是，当前的最大后验估计方法往往先要利用平滑滤波器与冲击滤波器预测图像边缘信息，然后再利用预测边缘图像估计模糊核，经过这两个过程的反复迭代最终估计获得模糊核，缺乏严格的最优化理论支撑，所以无法保证解的（局部）最优性。

发明内容

本发明的目的在于解决现有技术方法在理论和实现两个方面的不足，提出一种基于L0稀疏先验的最大后验模糊核估计新方法。

本发明技术方案是：首先，引入基于L0范数的显著边缘稀疏先验，利用迭代硬阈值收缩法实现显著边缘特征的隐性自动预测；其次，引入基于L0范数的相机抖动模糊核稀疏先验，利用迭代重新加权最小二乘法实现模糊核的快速估计；最后，利用基于超拉普拉斯先验的图像非盲去模糊方法，实现高质量图像去模糊。

1、本发明具体实施步骤：

(1)利用水平方向和垂直方向的一阶导数算子▽_h＝[1,-1;0,0],▽_v＝[1,0;-1,0]，得到相机抖动模糊图像y的梯度图像$>y_d={▿}_d\otimes y,d\in \Lambda ,\Lambda =\{h,v\};$>

(2)令待估计模糊核k的尺寸为Z×Z，为提高模糊核估计方法的收敛性，采用多尺度实现方式迭代估计模糊核；

(3)令初始模糊核k(1)的尺寸为3×3，且k⁽¹⁾＝[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(y_d,1)*ksize(s)/Z);

(4.2)   c^(s)=floor(size(y_d,2)*ksize(s)/Z);

(4.3)$>y_d^{\left(s\right)}=\mathrm{imresize}(y_d,[r^{\left(s\right)},c^{\left(s\right)}],\mathrm{bilinear}_{\prime }^{\prime });$>

(5)令显著边缘图像与模糊核在各尺度下的循环迭代次数为out_iter_num，循环迭代初始次数为out_iter=1，初始尺度为s=1，显著边缘L0稀疏先验的正则化参数为η，模糊核L0稀疏先验的正则化参数为γ，K^old对应模糊核k⁽¹⁾的二维矩阵，对应模糊插值图像的一维向量，k^old对应模糊核k⁽¹⁾的一维向量；

(6)令τ^old＝1，利用迭代硬阈值收缩法估计第s个尺度下的显著边缘主要利用步骤(6.1)-(6.4)进行m次循环估计：

(6.1)$>{\tilde{f}}_d^{\mathrm{new}}=f_d^{\mathrm{old}}-{\tau }^{\mathrm{old}}·{\left(K^{\mathrm{old}}\right)}^T(K^{\mathrm{old}}{f}_d^{\mathrm{old}}-y_d^{\left(s\right)});$>

(6.2)$>{\left(f_d^{\mathrm{new}}\right)}_j={\Phi }_{\mathrm{HARD}}({\left({\tilde{f}}_d^{\mathrm{new}}\right)}_j,{(\eta /2)}^{-1/2});$>

(6.3)τ^new＝τ^old2;

(6.4)$>f_d^{\mathrm{old}}=f_d^{\mathrm{new}};$>

其中，硬阈值算子Φ_HARD(·,·)定义如下：

$>{\Phi }_{\mathrm{HARD}}(a,b)=\left(\begin{array}{cc}a,& \mathrm{if}\left|a\right|\ge b\\ 0,& \mathrm{if}\left|a\right|<b\end{array}\right);$>

(7)令权矩阵利用迭代重新加权最小二乘方法估计第s个尺度下的模糊核k^(s)=k^new，具体而言利用步骤(7.1)-(7.4)进行n次循环估计：

(7.1)估计$>k^{\mathrm{new}}=\gamma ·{(W+\gamma ·{\Sigma }_{d\in \Lambda }{\left(F_d^{\mathrm{old}}\right)}^T\left(F_d^{\mathrm{old}}\right))}^{-1}{\Sigma }_{d\in \Lambda }{\left(F_d^{\mathrm{old}}\right)}^T{y}_d^{\left(s\right)};$>

(7.2)将k^new投影到约束集$>\left\{k\right|k\ge 0,{\Sigma }_{l=1}^{L}k\left(l\right)=1\};$>

(7.3)更新$>W=\mathrm{diag}\left(\right\{1/{(k_l^{\mathrm{new}}+0.0001)}^2\left\}\right);$>

(7.4)k^old=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)；

(13)输出最终估计的模糊核

(14)利用基于超拉普拉斯先验的图像非盲去模糊方法，最终获得去模糊图像

有益效果：

（1）发明方法的模糊核估计是个严格意义上的稀疏最优化问题；

（2）发明方法的实现方式简单，无需平滑滤波、冲击滤波等任何迭代预处理；

（3）发明方法的计算复杂度低；

（4）发明方法的模糊核估计准确度高；

估计出相机抖动模糊核之后，利用基于超拉普拉斯先验的图像非盲去模糊方法，从而获得更高质量的去模糊图像。

2、本发明模型推导过程：

不失一般性，相机抖动模糊可利用如下卷积型观察模型进行表述

$>y=k\otimes f+n$>

其中，y是获取的相机抖动模糊图像，f是真实清晰图像，k是空间不变的相机抖动模糊核，n是服从高斯分布的加性随机噪声，代表卷积；为了表述的方便，该观察模型还可写成如下矩阵-向量表达形式

y=Kf+n

其中，y，f，以及n分别是y，f，以及n的向量表达形式，K是模糊核k的矩阵表达形式。采取分而治之的策略解决相机抖动去模糊问题，分成两大步骤：(1)模糊核估计；(2)非盲图像恢复；

模糊核估计是在图像梯度域实现的，为此，建立如下梯度域观察模型：

$>y_d=k\otimes f_d+n_d$>

其中，$>y_d={▿}_d\otimes y,f_d={▿}_d\otimes f,n_d={▿}_d\otimes n,d\in \Lambda ,\Lambda =\{h,v\},\mathrm{and}{▿}_h=[1,-1;\mathrm{0,0}],{▿}_v=[\mathrm{1,0};-\mathrm{1,0}].$>类似的，上述梯度域卷积型观察模型可以写成如下矩阵-向量表达形式

y_d＝Kf_d+n_d＝F_dk+n_d

其中，y_d,f_d,n_d,k分别是y_d,f_d,n_d,k的向量表达形式，F_d是f_d的矩阵表达形式。

图像中的显著边缘是精确估计模糊核的重要所在。为了便于自动预测显著边缘，提高模糊核估计的准确性，首先提出基于L0范数的显著边缘先验及其优化模型：

$>{\left\{f_d^{\mathrm{est}}\right\}}_{d\in \Lambda }=\arg \underset{{\left\{f_d\right\}}_{d\in \Lambda }}{\min }\{\Xi \left({\left\{f_d\right\}}_{d\in \Lambda }\right)+\frac{\eta }{2}·{\Sigma }_{d\in \Lambda }{\left|\right|y_d-{\mathrm{Kf}}_d\left|\right|}_2^2\}$>

其中，是对应f_d的最优化结果，η是正则化参数，Ξ({f_d}_d∈Λ)是基于L0范数的显著边缘先验，定义为：

$>\Xi \left({\left\{f_d\right\}}_{d\in \Lambda }\right)={\Sigma }_{d\in \Lambda }{\left|\right|f_d\left|\right|}_0$>

根据迭代硬阈值收缩方法，上述最优化问题可以下步骤(1)-(3)进行迭代求解：

(1)$>f_d^{\mathrm{new}}={\Phi }_{\mathrm{HARD}}(f_d^{\mathrm{old}}-{\tau }^{\mathrm{old}}·{\left(K^{\mathrm{old}}\right)}^T(K^{\mathrm{old}}{f}_d^{\mathrm{old}}-y_d),{(\eta /2)}^{-1/2}),$>其中，K^old对应之前估计的模糊核，对应之前估计的显著边缘，对应当前更新的显著边缘，Φ_HARD(·,·)代表硬阈值算子，定义如下：

$>{\Phi }_{\mathrm{HARD}}(a,b)=\left(\begin{array}{cc}a,& \mathrm{if}\left|a\right|\ge b\\ 0,& \mathrm{if}\left|a\right|<b\end{array}\right);$>

(2)τ^new＝τ^old2，τ^old对应当前的调整参数，τ^new代表更新的调整参数；

(3)τ^old＝τ^new；

当显著边缘给定时，利用相机抖动模糊核的稀疏特性，从而能够建立基于L0范数的模糊核先验及其优化模型：

$>k^{\mathrm{est}}=\arg \underset{k}{\min }\{\Xi \left(k\right)+\frac{\gamma }{2}·{\Sigma }_{d\in \Lambda }{\left|\right|y_d-F_{d}k\left|\right|}_2^2\}$>

其中，Ξ(k)＝k₀是基于L0范数的模糊核先验，k^est是对应k的最优化结果，γ是正则化参数；

根据迭代重新加权最小二乘方法，首先上述将最优化问题转化为如下近似形式：

$>k^{\mathrm{est}}=\arg \underset{k}{\min }\{{\left|\right|k\left|\right|}_W^2+\frac{\gamma }{2}·{\Sigma }_{d\in \Lambda }{\left|\right|y_d-F_{d}k\left|\right|}_2^2\}$>

其中，W＝diag({1(k_l+0.0001)²}，然后利用以下步骤(4)-(6)进行迭代求解：

(4)利用共轭梯度法估计$>k^{\mathrm{new}}=\gamma ·{(W+\gamma ·{\Sigma }_{d\in \Lambda }{\left(F_d^{\mathrm{old}}\right)}^T\left(F_d^{\mathrm{old}}\right))}^T{\Sigma }_{d\in \Lambda }{\left(F_d^{\mathrm{old}}\right)}^T{y}_d$>其中，W为之前估计的权矩阵为之前估计的显著边缘的矩阵表达形式，k^new为当前更新的模糊核；

(5)进一步将k^new投影到约束集$>\left\{k\right|k\ge 0,{\Sigma }_{l=1}^{L}k\left(l\right)=1\};$>

(6)更新$>W=\mathrm{diag}\left(\right\{1/{(k_l^{\mathrm{new}}+0.0001)}^2\left\}\right).$>

利用上述方法原理，采用多尺度实现方式迭代估计显著边缘图像与相机抖动模糊核，从而得到最终的估计模糊核。

附图说明

附图1.发明方法流程图；

附图2.仿真实验的真实图像Lena及模糊核；

附图3.仿真实验的真实图像Cameraman及模糊核；

附图4.仿真实验的真实图像House及模糊核；

附图5.仿真实验的真实图像Boat及模糊核；

附图6.对应图像Lena的相机抖动模糊图；

附图7.基于规范化稀疏度量盲去模糊方法的恢复图像及估计模糊核；

附图8.基于L0稀疏先验的盲去模糊方法的恢复图像及估计模糊核；

附图9.基于L0稀疏先验的盲去模糊方法的显著边缘图像。

具体实施方式

(1)利用水平方向和垂直方向的一阶导数算子▽_h＝[1,-1;0,0],▽_v＝[1,0;-1,0]，利用MATLAB函数conv2获得相机抖动模糊图像y的梯度图像$>y_d={▿}_d\otimes y,d\in \Lambda ,\Lambda =\{h,v\}:$>

y_h=conv2(y,▽_h,'valid');

y_v=conv2(y,▽_v,'valid');

(2)令待估计模糊核k的尺寸为Z×Z，为提高模糊核估计方法的收敛性，采用多尺度实现方式迭代估计模糊核；

(3)令初始模糊核k⁽¹⁾＝[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(y_d,1)*ksize(s)/Z);

c^(s)=floor(size(y_d,2)*ksize(s)/Z);

$>y_d^{\left(s\right)}=\mathrm{imresize}(y_d,[r^{\left(s\right)},c^{\left(s\right)}],\mathrm{bilinear}_{\prime }^{\prime });$>

(5)令显著边缘图像与模糊核在各尺度下的循环迭代次数为out_iter_num，循环迭代初始次数为out_iter=1，初始尺度为s=1，显著边缘L0稀疏先验的正则化参数为η，模糊核L0稀疏先验的正则化参数为γ，K^old对应模糊核k(1)的二维矩阵，k^old对应模糊核k⁽¹⁾的一维向量，对应模糊插值图像的一维向量；

(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次循环估计：

$>{\tilde{f}}_d^{\mathrm{new}}=f_d^{\mathrm{old}}-{\tau }^{\mathrm{old}}·{\left(K^{\mathrm{old}}\right)}^T(K^{\mathrm{old}}{f}_d^{\mathrm{old}}-y_d^{\left(s\right)});$>

$>{\left(f_d^{\mathrm{new}}\right)}_j={\Phi }_{\mathrm{HARD}}({\left({\tilde{f}}_d^{\mathrm{new}}\right)}_j,{(\eta /2)}^{-1/2})$>

$>=\left(\begin{array}{cc}{\left({\tilde{f}}_d^{\mathrm{new}}\right)}_j,& \mathrm{if}\left|{\left({\tilde{f}}_d^{\mathrm{new}}\right)}_j\right|\ge {(\eta /2)}^{-1/2}\\ 0,& \mathrm{if}\left|{\left({\tilde{f}}_d^{\mathrm{new}}\right)}_j\right|<{(\eta /2)}^{-1/2}\end{array}\right);$>

τ^new＝τ^old/2;

$>f_d^{\mathrm{old}}=f_d^{\mathrm{new}};$>

τ^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（共轭梯度法）求解如下方程组，估计k^new:

$>(W+\gamma ·{\Sigma }_{d\in \Lambda }{\left(F_d^{\mathrm{old}}\right)}^T\left(F_d^{\mathrm{old}}\right))k^{\mathrm{new}}=\gamma ·{\Sigma }_{d\in \Lambda }{\left(F_d^{\mathrm{old}}\right)}^T{y}_d^{\left(s\right)},$>

共轭梯度法可参见参考文献《Convex optimization》(S.Boyd,L.Vandenberghe)，Cambridge University Press，2004。

(7.2)将knew投影到约束集主要利用以下MATLAB代码实现：

k^new(k^new<0)=0;

sumk=sum(k^new(:));

k^new=k^new./sumk;

(7.3)更新$>W=\mathrm{diag}\left(\right\{1/{(k_l^{\mathrm{new}}+0.0001)}^2\left\}\right);$>

(7.4)k^old=k^new;

其中，对应显著边缘的二维矩阵，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)；

(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稀疏先验的盲去模糊方法的显著边缘图像。

以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内，本发明要求保护范围由所附的权利要求书其等效物界定。