基于分裂Bregman迭代的快速鲁棒图像运动去模糊方法

摘要

本发明属于数字图像处理领域，提供了一种基于分裂Bregman迭代的快速鲁棒图像运动去模糊方法，通过直接利用图像梯度与运动模糊核的L

说明书

技术领域

本发明属于数字图像处理领域，具体涉及利用相机捕获的单模糊图像自动估计对应相机各种随机抖动的点扩散函数(point spread function)或称为运动模糊核(motion blurkernel)的方法。

背景技术

在相机拍摄过程中，由于某些不可控因素造成的相机随机抖动往往导致捕获图像呈现运动模糊的现象，最典型的情形是相机在弱光环境中的长时间曝光发生的随机抖动。

处理运动模糊图像的技术核心是自动估计对应相机各种随机抖动的点扩散函数。目前，绝大多数点扩散函数估计方法都是基于贝叶斯统计框架，根据推断准则的不同，主要分为两大类：均值场变分近似估计方法(mean field variational approximation estimation)和最大后验估计方法(maximum a posterior estimation)。最近，Krishnan等人指出：上述两类估计方法实质上遵循了相似的最优化准则，核心都是要为自然图像和点扩散函数赋予适当的稀疏先验模型，参见文献《Blind deconvolution with re-weighted sparsitypromotion》,arXiv:1311.4029,2013。

均值场变分近似方法能够较精确地估计运动模糊核，但是计算复杂度较高。相比而言，运动模糊核函数的最大后验估计方法的计算复杂度要低得多。此外，该类方法更容易理解、包容性更强。

然而，当前运动模糊核估计方法中的图像先验模型往往都是高度非凸的，且其具有一个共同的特点，即：本质上都是利用不同的技巧实现对L₀范数的一种近似逼近。

发明内容

本发明针对上述现有技术存在的问题作出改进，即本发明要解决的技术问题是提供一种基于分裂Bregman迭代的快速鲁棒图像运动去模糊方法，这种方法旨在避免现有最大后验估计方法在建模方面和均值场变分近似方法在实现方面的不足，从而使得图像运动去模糊技术具有更强的实用性。

为了解决上述技术问题，本发明提供了如下的技术方案：

一种基于分裂Bregman迭代的快速鲁棒图像运动去模糊方法，包括如下步骤：

首先，通过直接利用图像梯度与运动模糊核的L₀范数并结合其各自的L₂范数，构建运动模糊核估计的非凸非光滑能量泛函；其次，通过耦合算子分裂和增广拉格朗日方法，设计运动模糊核的分裂Bregman迭代求解格式；最后，利用基于全变差先验的图像非盲去模糊方法，实现图像的快速去模糊。

本发明具体实施步骤：

(1)给定待运动去模糊图像y，给定待估计运动模糊核h的尺寸为Z×Z；

(2)采用多尺度实现方式迭代估计模糊核，设定尺度总数为S＝4；

(3)令y⁽⁴⁾＝y，利用下列MATLAB代码计算其它尺度下的运动模糊图像y^(s)(1≤s≤3)：

(3.1)for s＝3:-1:1

(3.2)y^(s)＝imresize(y^(s+1),0.5)；

(3.3)end

(4)利用下列MATLAB代码初始化运动模糊核h⁽⁰⁾：

(4.1)hsize＝ceil(Z/2^(3))；

(4.2)cen＝floor((hsize+1)/2)；

(4.3)h⁽⁰⁾＝zeros(hsize)；

(4.4)h⁽⁰⁾(cen(1),cen(2))＝1；

(5)设定参数λ,β_u,β_k,τ_u,τ_k的取值，其中，λ为保证项的参数，β_u为图像梯度L₀范数的参数，β_k为运动模糊核L₀范数的参数，τ_u为图像梯度L₂范数的参数，τ_k为运动模糊核L₂范数的参数；

(6)设定每个尺度下的内、外部循环迭代次数分别为10、10，内、外部循环迭代初始次数l、i分别取为0，且初始尺度s取为1；

(7)令o＝y^(s)，且全部设定为0向量，利用下列分裂Bregman迭代方法估计对应每个尺度s的运动模糊核

(7.1)for i＝0:1:9

(7.2)for l＝0:1:9

$>\left(7.3\right)---w_i^{l+1}={\Theta }_{\mathrm{Hard}}(▿u_i^l+\frac{1}{{\gamma }_u}{\mu }_u^l,{\left(\frac{2{\beta }_u}{{\gamma }_u}\right)}^{\frac{1}{2}}),$>

$>\left(7.4\right)---u_i^{l+1}={(\lambda K_i^{*}{K}_i+({\tau }_u+\frac{{\gamma }_u}{2})▿*▿)}^{-1}(\lambda K_i^{*}o+\frac{{\gamma }_u}{2}▿*(w_i^l-\frac{1}{{\gamma }_u}{\mu }_u^l)),$>

$>\left(7.5\right)---g_i^{l+1}={\Theta }_{\mathrm{Hard}}(k_i^l+\frac{1}{{\gamma }_k}{\mu }_k^l,{\left(\frac{2{\beta }_k}{{\gamma }_k}\right)}^{\frac{1}{2}}),$>

$>\left(7.6\right)---k_i^{l+1}={(\lambda U_i^{*}{U}_i+({\tau }_k+\frac{{\gamma }_k}{2})I)}^{-1}(\lambda U_i^{*}o+\frac{{\gamma }_k}{2}(g_i^l-\frac{1}{{\gamma }_k}{\mu }_k^l)),$>

$>\left(7.7\right)---{\mu }_u^{l+1}={\mu }_u^l+{\gamma }_u(▿u_i^{l+1}-w_i^{l+1}),$>

$>\left(7.8\right)---{\mu }_u^{l+1}={\mu }_k^l+{\gamma }_k(k_i^{l+1}-g_i^{l+1}).$>

(7.9)end

(7.10)end

其中，▽为对应水平方向和垂直方向一阶导数算子▽_h＝[1,-1；0,0],▽_v＝[1,0；-1,0]的卷积算子，U_i,K_i为对应卡通图像和运动模糊核的卷积算子，γ_u,γ_k为设定的增广拉格朗日惩罚算子，硬阈值算子Θ_Hard(·,·)定义为Θ_Hard(α,T)＝α·[α≥T]，且在傅立叶变换域计算式(7.2)、(7.4)，

式(7.2)、(7.4)的傅立叶变换域计算公式如下：

其中，代表a的二维傅立叶变换，ifft2为MATLAB中的二维逆傅立叶变换函数；

(8)利用下列MATLAB代码，将h^(s)投影到约束集{h|h≥0,Σ_rΣ_th(r,t)＝1}：

(8.1)h^(s)(h^(s)<0)＝0；

(8.2)sumh＝sum(h^(s)(:))；

(8.3)h^(s)＝h^(s)./sumh；

(9)输出最终估计的运动模糊核

(10)利用基于全变差先验的图像非盲去模糊方法，最终获得去模糊图像

本发明的有益效果：

(1)本发明方法的运动模糊核估计基于混合L₀、L₂范数的严格稀疏最优化问题；

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

(3)本发明方法的实时性强；

(4)本发明方法的准确度高；

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：

图1为本发明方法流程图；

图2为不同运动去模糊方法的PSNR互补累计直方图；

图3为本发明方法和Xu&Jia估计的32个运动模糊核；

图4为不同运动去模糊方法恢复标准测试集第16幅运动模糊图像的实验结果；

图5为本发明方法在最大尺度下每次外部迭代的图像与模糊核混合L₀、L₂范数的能量曲线图。

具体实施方式

如图1-5所示，本发明公开一种基于分裂Bregman迭代的快速鲁棒图像运动去模糊方法，具体步骤为：

首先，通过直接利用图像梯度与运动模糊核的L₀范数并结合所述图像梯度与运动模糊核各自的L₂范数，构建运动模糊核估计的非凸非光滑能量泛函；

其次，通过耦合算子分裂和增广拉格朗日方法，设计运动模糊核的分裂Bregman迭代求解格式；

最后，利用基于全变差先验的图像非盲去模糊方法，实现图像的快速去模糊。

本发明具体实施步骤：

(1)给定待运动去模糊图像y，给定待估计运动模糊核h的尺寸为Z×Z；

(2)为避免模糊核估计陷入无效的局部极小点，采用多尺度实现方式迭代估计模糊核，设定尺度总数为S；

(3)令y^(S)＝y，计算其它尺度下的运动模糊图像y^(s)(1≤s≤S-1)；

(4)初始化运动模糊核h⁽⁰⁾；

(5)设定参数λ,β_u,β_k,τ_u,τ_k的取值，其中，λ为保证项的参数，β_u为图像梯度L₀范数的参数，β_k为运动模糊核L₀范数的参数，τ_u为图像梯度L₂范数的参数，τ_k为运动模糊核L₂范数的参数；

(6)设定每个尺度下的内、外部循环迭代次数分别为L、I，内、外部循环迭代初始次数l、i分别取为0，且初始尺度s取为1；

(7)令o＝y^(s)，且全部设定为0向量，利用分裂Bregman迭代方法估计对应每个尺度s的运动模糊核具体步骤如(7.1)-(7.6)所示：

for i＝0:1:I-1

for l＝0:1:L-1

$>\left(7.1\right)---w_i^{l+1}={\Theta }_{\mathrm{Hard}}(▿u_i^l+\frac{1}{{\gamma }_u}{\mu }_u^l,{\left(\frac{2{\beta }_u}{{\gamma }_u}\right)}^{\frac{1}{2}}),$>

$>\left(7.2\right)---u_i^{l+1}={(\lambda K_i^{*}{K}_i+({\tau }_u+\frac{{\gamma }_u}{2})▿*▿)}^{-1}(\lambda K_i^{*}o+\frac{{\gamma }_u}{2}▿*(w_i^l-\frac{1}{{\gamma }_u}{\mu }_u^l)),$>

$>\left(7.3\right)---w_i^{l+1}={\Theta }_{\mathrm{Hard}}(▿u_i^l+\frac{1}{{\gamma }_u}{\mu }_u^l,{\left(\frac{2{\beta }_u}{{\gamma }_u}\right)}^{\frac{1}{2}}),$>

$>\left(7.4\right)---u_i^{l+1}={(\lambda K_i^{*}{K}_i+({\tau }_u+\frac{{\gamma }_u}{2})▿*▿)}^{-1}(\lambda K_i^{*}o+\frac{{\gamma }_u}{2}▿*(w_i^l-\frac{1}{{\gamma }_u}{\mu }_u^l)),$>

$>\left(7.5\right)---g_i^{l+1}={\Theta }_{\mathrm{Hard}}(k_i^l+\frac{1}{{\gamma }_k}{\mu }_k^l,{\left(\frac{2{\beta }_k}{{\gamma }_k}\right)}^{\frac{1}{2}}),$>

$>\left(7.6\right)---{\mu }_k^{l+1}={\mu }_k^l+{\gamma }_k(k_i^{l+1}-g_i^{l+1}).$>

end

end

其中，▽为对应水平方向和垂直方向一阶导数算子▽_h＝[1,-1；0,0],▽_v＝[1,0；-1,0]的卷积算子，U_i,K_i为对应卡通图像和运动模糊核的卷积算子，γ_u,γ_k为设定的增广拉格朗日惩罚算子，硬阈值算子Θ_Hard(·,·)定义为Θ_Hard(α,T)＝α·[α≥T]，且在傅立叶变换域计算式(7.2)、(7.4)；

(8)将h^(s)投影到约束集{h|h≥0,Σ_rΣ_th(r,t)＝1}；

(9)输出最终估计的运动模糊核

(10)利用基于全变差先验的图像非盲去模糊方法，最终获得去模糊图像

其中，本发明模型推导过程：

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

y＝h*x+n

其中，y是相机随机抖动后获取的模糊图像，x是原始清晰图像，h是空间不变的对应相机随机抖动的运动模糊核，n是服从高斯分布的加性随机噪声，*代表卷积算子；与当前方法类似，本发明采取分而治之的策略解决相机随机抖动去模糊问题，主要分为两大步：(1)运动模糊核估计；(2)非盲图像去模糊；

图像运动去模糊是个严重不适定的数学反问题，为了实现运动模糊核的稳定有效的估计，需要为自然图像和点扩散函数赋予适当的稀疏先验模型。

卡通图像或者说图像中的显著边缘是精确估计运动模糊核的重要所在。为了有效区分模糊图像和清晰图像，提高运动模糊核估计的准确性和稳定性，提出基于混合L₀、L₂范数的图像先验模型：

其中，▽为对应水平方向和垂直方向一阶导数算子的卷积算子，β_u,τ_u是正则化参数。

此外，考虑到运动模糊核自身的稀疏性物理特性以及计算稳定性，提出基于混合L₀、L₂范数的运动模糊核先验模型：

其中，β_k,τ_k是正则化参数。

基于新提出的图像和运动模糊核先验模型，构建运动模糊核估计的能量最小化泛函：

$>{\min }_{x,h}\{\lambda {\left|\right|h*x-y\left|\right|}_2^2+{\beta }_u{\left|\right|▿x\left|\right|}_0+{\tau }_u{\left|\right|▿x\left|\right|}_2^2+{\beta }_k{\left|\right|h\left|\right|}_0+{\tau }_k{\left|\right|h\left|\right|}_2^2\}.$>

以下介绍对应尺度为s的能量最小化泛函求解方法。令o＝y^(s)，k＝h^(s)，u＝x^(s)，则能量最小化问题转化为：

$>{\min }_{u,k}\{\lambda {\left|\right|k*u-o\left|\right|}_2^2+{\beta }_u{\left|\right|▿u\left|\right|}_0+{\tau }_u{\left|\right|▿u\left|\right|}_2^2+{\beta }_k{\left|\right|k\left|\right|}_0+{\tau }_k{\left|\right|k\left|\right|}_2^2\}.$>

根据分裂Bregman迭代思想，上述最优化问题可利用以下方法进行迭代求解。

首先，利用算子分裂方法，将最优化问题等价地转化为如下的约束L₀最小化问题：

min_w,u,g,kΨ(w,u,g,k)s.t.w＝▽u,g＝k，

其中，Ψ(w,u,g,k)定义为

$>\lambda {\left|\right|k*u-o\left|\right|}_2^2+{\beta }_u{\left|\right|w\left|\right|}_0+{\beta }_k{\left|\right|g\left|\right|}_0+{\tau }_u{\left|\right|▿u\left|\right|}_2^2+{\tau }_u{\left|\right|k\left|\right|}_2^2.$>

此时，可以通过交替迭代方式估计w,u以及g,k：

(w_i+1,u_i+1)＝arg>w,uΨ(w,u,g_i,k_i)s.t.w＝▽u,

(g_i+1,k_i+1)＝arg>g,kΨ(w_i,u_i,g,k)s.t.g＝k,

其中，0≤i≤I-1。

然后，利用增广拉格朗日方法，通过求解如下非约束最小化问题迭代估计w_i+1,u_i+1,g_i+1,k_i+1：

$>w_i^{l+1},u_i^{l+1}=\arg {\min }_{w,u}{\bar{\Psi }}^l(w,u,g_i,k_i),$>

$>g_i^{l+1},k_i^{l+1}=\arg {\min }_{g,k}{\bar{\Psi }}^l(w_i,u_i,g,k),$>

其中，l∈[0,L-1]，且定义为

$>\Psi (w,u,g,k)+{\mu }_u^{l*}(▿u-w)+\frac{{\gamma }_u}{2}{\left|\right|▿u-w\left|\right|}_2^2+{\mu }_k^{l*}(k-g)+\frac{{\gamma }_k}{2}{\left|\right|k-g\left|\right|}_2^2,$>

其中，γ_u,γ_k为增广拉格朗日惩罚算子，为约束w＝▽u和g＝k的拉格朗日乘子，且利用如下格式进行迭代更新：

$>{\mu }_u^{l+1}={\mu }_u^l+{\gamma }_u(▿u_i^{l+1}-w_i^{l+1}),$>

$>{\mu }_k^{l+1}={\mu }_k^l+{\gamma }_k(k_i^{l+1}-g_i^{l+1}).$>

最后，通过简单直接的计算，便可利用以下各式估计

$>w_i^{l+1}={\arg \min }_w{\bar{\Psi }}^l(w,u_i^l,g_i,k_i)={\Theta }_{\mathrm{Hard}}(▿u_i^l+\frac{1}{{\gamma }_u}{\mu }_u^l,{\left(\frac{2{\beta }_u}{{\gamma }_u}\right)}^{\frac{1}{2}}),$>

$>u_i^{l+1}={\arg \min }_u{\overline{\text{Ψ}}}^l(w_i^l,u,g_i,k_i)={(\lambda K_i^{*}{K}_i+({\tau }_u+\frac{{\gamma }_u}{2})▿*▿)}^{-1}(\lambda K_i^{*}o+\frac{{\gamma }_u}{2}▿*(w_i^l-\frac{1}{{\gamma }_u}{\mu }_u^l)),$>

$>g_i^{l+1}={\arg \min }_g{\bar{\Psi }}^l(w_i,u_i,g,k_i^l)={\Theta }_{\mathrm{Hard}}(k_i^l+\frac{1}{{\gamma }_k}{\mu }_k^l,{\left(\frac{2{\beta }_k}{{\gamma }_k}\right)}^{\frac{1}{2}}),$>

$>k_i^{l+1}={\arg \min }_k{\bar{\Psi }}^l(w_i,u_i,g_i^l,k)={(\lambda U_i^{*}{U}_i+({\tau }_k+\frac{{\gamma }_k}{2})I)}^{-1}(\lambda U_i^{*}o+\frac{{\gamma }_k}{2}(g_i^l-\frac{1}{{\gamma }_k}{\mu }_k^l)),$>

其中，为给定初始值，U_i,K_i为对应卡通图像和运动模糊核的卷积算子，且硬阈值算子Θ_Hard(·,·)定义为Θ_Hard(α,T)＝α·[α≥T]。

利用上述方法原理，采用多尺度实现方式迭代估计运动模糊核，从而能够得到最终的估计模糊核

综上，本发明对原有模型进行了提升，通过引入支撑连续性先验，提高了运动模糊核的估计精度；设计了一种新的基于分裂Bregman迭代的快速解法，大幅提高了运动模糊核的估计效率。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。