基于TV和小波正则化的空间变化模糊图像复原方法

摘要

本发明公开了一种基于TV和小波正则化的空间变化模糊图像复原方法，包括：(1)对模糊图像进行灰度化；(2)根据灰度化的模糊图像构建模糊核分解模型，在所述模糊核分解模型中，采用奇异值分解方法将模糊核分解为基滤波矩阵和系数矩阵；(3)应用所述模糊核分解模型，并结合TV正则项和小波正则项构建去模糊模型；(4)将所述去模糊模型转化为增广拉格朗日形式后，在对增广拉格朗日形式的去模糊模型进行改进，得到新去模糊模型；(5)采用ADMM算法对新去模糊模型进行数值迭代求解，获得复原图像。该方法解决了TV正则化算法在复原过程中产生的细节信息丢失问题。

说明书

技术领域

本发明属于图像处理技术领域，具体涉及一种基于TV和小波正则化的空间变化模糊图像复原方法。

背景技术

图像去模糊旨在通过算法对被噪声污染的图像进行某种处理，以降低噪声对原始有用信息的影响，根据退化图像尽可能地求解出清晰图像，具体可以分为三大类，分别为图像增强、图像复原和超分辨率重构。

图像复原是以获取视觉质量某种程度的改善为目的，根据某些特定的图像退化模型进行估计计算，以实现对退化图像的复原。

对于许多成像装置来说，他们的图像退化模型虽然可以被认为是线性的，但在图像获取传输过程中，由成像系统、传输介质方面的原因造成图像退化模型并不是空间不变(Space-invariant，SI)的，而是空间变化的 (Space-variant，SV)。空间变化退化模型与空间不变退化模型差异主要表现在：空间不变的模糊核(Point spread function，PSF)在图像退化模型中是固定不变，而空间变化的模糊核在图像的不同区域是不同的。

申请公布号为CN105741243A的专利申请公开了一种模糊图像复原方法，应用于由相机和被拍摄目标相对运动形成的模糊图像，该方法包括：通过图像得到图像的平滑区域；通过图像的平滑区域计算平滑区域标记矩阵；根据图像生成初始模糊核；通过对初始模糊核优化得到模糊核的估计值；通过模糊核的估计值对模糊图像进行复原。该方法通过对图像的平滑区域进行标记和充分利用图像中的平滑区域信息对模糊核估计过程施加约束，保证了模糊核估计的准确性，并利用估计的模糊核对模糊图像进行复原。该方法采用了空间变化的模糊核对模糊图像进行复原。

在图像复原过程中，图像上的一点点噪声可能就会对复原的结果产生非常大的影响，因为很多复原算法都会放大噪声。这时候需要在最优化问题模型中添加一些正则项来保持图像的光滑性，全变分(Total Variation， TV)是常用的一种正则项。TV项用在图像复原和去噪中的作用就是保持图像的光滑性，消除图像复原可能带来的伪影。

由于TV正则项通常在平坦区域内和边缘采取不同平滑速度的非线性策略，让其在平坦区域加速平滑，在边缘处抑制平滑，以保护图像边缘，导致了在去模糊过程中阶梯效应的产生，导致细节信息的丢失，使复原图像中出现虚假边缘现象。

奇异值分解是一种基于特征向量的矩阵变换方法，在信号处理、模式识别、数字水印技术等方面都得到了应用。

发明内容

本发明的目的是提供一种基于TV和小波正则化的空间变化模糊图像复原方法。通过添加小波正则项，利用小波良好的重建能力来补充TV正则化算法在复原过程中产生的细节信息丢失问题，同时TV正则项可以解决小波分解产生的边缘模糊问题。

为实现上述发明目的，本发明提供以下技术方案：

一种基于TV和小波正则化的空间变化模糊图像复原方法，包括以下步骤：

(1)对模糊图像进行灰度化；

(2)根据灰度化的模糊图像构建模糊核分解模型，在所述模糊核分解模型中，采用奇异值分解方法将模糊核分解为基滤波矩阵和系数矩阵；

(3)应用所述模糊核分解模型，并结合TV正则项和小波正则项构建去模糊模型；

(4)将所述去模糊模型转化为增广拉格朗日形式的去模糊模型后，在对增广拉格朗日形式的去模糊模型进行改进，得到新去模糊模型；

(5)采用ADMM算法对新去模糊模型进行数值迭代求解，获得复原图像。

本发明具有的有益效果为：

建立基于奇异值分解技术的模糊核分解模型，该模型解决了常用的加性卷积模型不适用于空间变化模糊图像复原的问题。采用TV函数和小波框架为正则项，将两者的线性组合作为先验信息用于图像复原算法中，使阶梯效应得到较好的抑制和纹理细节的保持。最后利用ADMM算法对图像去模糊图像进行求解，实现模型快速收敛。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图做简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动前提下，还可以根据这些附图获得其他附图。

图1是实施例提供的空间变化模糊图像复原方法的流程图；

图2是实验图像，其中，(a)为清晰图像，(b)～(g)为不同模糊程度的模糊图像；

图3是对比算法对模糊图像的复原结果图，其中，(a)～(f)是对图2中的(b)～(g)模糊图像的复原结果图；

图4是本发明方法对模糊图像的复原结果图，其中，(a)～(f)是对图2中的(b)～(g)模糊图像的复原结果图。

具体实施方式

为使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例对本发明进行进一步的详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不限定本发明的保护范围。

图1为空间变化模糊图像复原方法的流程图。如图1所示，该空间变化模糊图像复原方法包括以下步骤：

S101，输入模糊图像g，并对模糊图像g进行灰度化；

S102，设置相关参数。

在采用ADMM算法对新去模糊模型进行数值迭代求解前，设置相关参数，具体包括分解基数k，保真项参数μ，正则项参数α和θ，ADMM 迭代辅助变量ω，v，H，q，p，t的初值和相应的增广拉格朗日项参数γ，β和ρ；ADMM的最大迭代次数MAX，并设其迭代次数i的初值为0，其中MAX的取值为1，2，3，……，N，N为正整数。

S103，根据灰度化的模糊图像构建模糊核分解模型。

具体地，在所述模糊核分解模型中，采用奇异值分解方法将模糊图像的模糊核分解为：

其中，A为待分解的模糊核矩阵，Y，Z分别为模糊核矩阵A的左奇异矩阵和右奇异矩阵，S为模糊核矩阵A的奇异值矩阵，分解出的奇异值个数为分解基数k，模糊核分解后的基滤波矩阵B＝Y，模糊核分解后的系数矩阵M＝S×Z^T。

S104，应用模糊核分解模型，并结合TV正则项和小波正则项构建去模糊模型。

具体地，构建的去模糊模型为：

E(u,g)＝Q(u,g)+λJ(u)

其中，Q(u,g)为保真项，λJ(u)为正则项，λ为正则化参数，用于平衡正则化项所占的权重，保真项Q(u,g)为：

其中，u表示清晰图像的矩阵表示；g表示模糊图像的矩阵表示；

正则项用来恢复图像的轮廓信息，正则项λJ(u)具体为：

λJ(u)＝α||Cu||₁+θ||Wu||₁

其中，C为差分算子；W为用来恢复图像的细节信息的小波框架变换，为简化后续计算，取W为紧框架变换，即WTW＝I(I表示单位矩阵)，α，θ分别为非负正则化参数；

将保真项和正则项带入到去模糊模型中，则将求解去模糊模型转化为求最小化问题，即

S105，将去模糊模型转化为增广拉格朗日形式后，在对增广拉格朗日形式的去模糊模型进行改进，得到新去模糊模型。

具体地，在去模糊模型中，引入辅助变量ω，v，H，转化为约束最优化问题：

将约束项ω＝Mu，v＝Cu，H＝Wu代入到去模糊模型中，获得增广拉格朗日形式的去模糊模型：

其中，q，p，t为增广拉格朗日形式引入的辅助变量。

采用ADMM算法(可以参见文献Wahlberg B,Boyd S,Annergren M,et al.An ADMMalgorithm for a class of total variation regularized estimation problems[J].IFAC Proceedings Volumes,2012,45(16):83-88)对增广拉格朗日形式的去模糊模型进行数值迭代求解，迭代算法如下：

求解式(2)的变量u，ω，v，H，只需在对应项令关于该变量的导数为零即可，即求解极小化问题，函数uⁱ⁺¹是关于u的极小化问题，可得：

(γM^TM+βC^TC+ρW^TW)u＝γM^T(ω+q)+βC^T(v+p)+ρW^T(H+t)(3)

若要求解式子(3)，需要求出(γM^TM+βC^TC+ρW^TW)^-1，其中W为紧小波框架分解算子，所以有W^TW＝I，I表示单位矩阵；C为差分算子，C^TC为拉普拉斯算子矩阵，因此，βC^TC和ρW^TW都可以在频域中快速计算。

然而，由于M^TM无法在频域中表示，(M^TM)^-1无法在频域中求解，所以在式子(1)中引入广义逆矩阵M⁺＝(M^TM)^-1M^T，将转化为即得到新去模糊模型：

S106，采用ADMM算法对新去模糊模型进行数值迭代求解，获得复原图像。

具体地，迭代过程如下式(5)所示：

式(5)中第一行所示的关于u的极小化问题的解如式(6)中第一行所示；式(5)中第二行为变量ω的极小化问题，则同样可以利用FFT 快速求解，如式(6)中第二行所示；式(5)中第三行、第四行分别是变量v和H的极小化问题，分别可以使用二维收缩算法和一维收缩算法快速求解，结果如式(6)中第三行和第四行所示。其中，二维收缩算法和一维收缩算法参见文献WANG Y,YANG J,YIN W,et al.A new alternating minimization algorithmfor total variation image reconstruction[J].SIAM journal on imaging sciences,2008,1(3):248-272.

式(6)中，F和F^-1分别代表傅里叶变换和傅里叶反变换。

具体地，求解的具体过程为：

(a)按照公式(3)的迭代过程求出uⁱ⁺¹，并计算uⁱ⁺¹的峰值信噪比(PSNR)>2和uⁱ的峰值信噪比P₁；

(b)令迭代次数i加1，即i＝i+1后，若迭代次数i大于等于最大迭代次数MAX，执行步骤(d)，若迭代次数i小于最大迭代次数MAX，执行步骤(c)；

(c)若峰值信噪比P₂大于峰值信噪比P₁，执行步骤(a)，若峰值信噪比P₂小于等于峰值信噪比P₁，执行步骤(d)；

(d)输出复原图像u_final＝u^i-1。

下面利用具体实验表明本发明提供的空间变化模糊图像复原方法的去噪效果。

实验条件：

实验原始图像分别选取分辨率为256×256的“Lena”图像如图2中(a)所示。实验编程环境为Matlab2016b。实验需要添加空间变化模糊，采用与对比算法(TV正则化算法)相同散焦模糊图像来模拟模糊图像，如图2所示，具体地，(b)～(g)为PSNR分别为18.5190、18.6569、19.0230、 19.1244、19.2189以及19.3087的模糊图像。

本发明的参数为：分解基数k＝10，保真项参数μ＝8e⁵，正则项参数α＝1、θ＝1、增广拉格朗日项参数β＝γ＝800、ρ＝6e^-3，和ADMM迭代辅助变量ω、v、H、q、p、t初始值均设为0，迭代次数i的初值设为0。ADMM>

每迭代一次，计算出迭代结果的评价值，即uⁱ⁺¹的峰值信噪比(PSNR)>

用来对比的空间变化TV正则化方法(即参考文献SROUBEK F, KAMENICKY J,LU YM.Decomposition of space-variant blur in image deconvolution.IEEE signalprocessing letters,2016,23(3)：346-350中提出的方法)参数：分解基数k＝10，保真项参数μ＝1e⁵，正则项参数α＝1，增广拉格朗日项参数ADMM迭代辅助变量ω、v、q、p 初始值均设为0；迭代次数i的初值设为0。ADMM的最大迭代次数 MAX＝500。每迭代一次，计算出迭代结果的评价值，即u₁ⁱ⁺¹的峰值信噪比>1ⁱ作为最终去模糊结果图像u_1final输出，如果迭代次数没有达到最大迭代次数MAX，则通过比较本次迭代后的输出图像u₁ⁱ⁺¹的峰值信噪比是否大于上一次迭代输出图像u₁ⁱ的峰值信噪比，来决定是否继续迭代：如果本次迭代后的输出图像u₁ⁱ⁺¹的峰值信噪比大于上一次迭代输出图像u₁ⁱ的峰值信噪比，则继续迭代；如果本次迭代结果u₁ⁱ⁺¹的峰值信噪比小于等于上一次迭代输出图像u₁ⁱ的峰值信噪比则将上一次迭代输出图像u₁ⁱ作为最优值u_1final输出。

实验内容：

按照上面所述的实验步骤进行实验仿真，并将本发明方法与空间变化TV正则化方法进行比较。

实验结果见表1和图3和图4，其中，图3为对比算法(空间变化TV 正则化方法)对图2中的(b)～(g)所示的模糊图像的复原结果图，具体地，(a)～(f)分别为PSNR为29.6647、29.5755、30.0782、29.8154、 29.6329以及30.0227的复原图像。图4为本发明方法对图2中的(b)～ (g)所示的模糊图像的复原结果图，具体地，(a)～(f)分别为PSNR 为32.0919、31.4006、31.1590、30.6378、29.5806以及30.0313的复原图像。

表1本发明方法和空间变化TV正则化方法去模糊后图像的峰值信噪比和结构相似度

经分析，实验结果表明本发明去空间变化模糊方法与空间变化TV正则化去模糊方法相比，本发明去模糊后图像的峰值信噪比和结构相似度更高，去模糊效果更好。

以上所述的具体实施方式对本发明的技术方案和有益效果进行了详细说明，应理解的是以上所述仅为本发明的最优选实施例，并不用于限制本发明，凡在本发明的原则范围内所做的任何修改、补充和等同替换等，均应包含在本发明的保护范围之内。