CN105260998A - 基于mcmc采样和阈值低秩逼近的图像去噪方法 - Google Patents

Abstract

本发明公开了一种基于MCMC和阈值低秩逼近的图像去噪方法。本方法在去噪过程中，利用蒙特卡罗采样产生图像块，并根据直方图多个统计特征，得到满足马尔科夫链的相似度判决函数；对各类图像块簇进行奇异值分解，并依据图像相应的先验信息，对奇异值进行自适应阈值估计；基于分解后的低秩结构，借助低秩逼近算法进行图像重构，从而达到去噪目的。本方法充分利用图像非局部几何结构信息相似的特性以及低秩矩阵能较好地处理高维数据的优势，克服了以往非局部均值遍历搜索方法选块复杂度高的缺点，增强了选块的鲁棒性，并在一定程度上减少了算法的运行时间。

Description

基于MCMC采样和阈值低秩逼近的图像去噪方法

技术领域

本发明涉及一种基于MCMC采样和阈值低秩逼近的图像去噪方法，属于视频图像数据处理的技术领域。

背景技术

图像是人类活动中最常用的信息载体，已成为人们获取外界原始信息的主要途径。然而，复杂的监测场景、恶劣的天气条件、光照变化以及采集设备的精度等易使图像受噪声污染严重、视觉效果模糊、目标特征难以分辨，从而对图像的后续处理产生影响，无法满足应用环境的需求。因此，作为图像预处理中不可或缺的环节，图像去噪对后续图像特征提取、分割、压缩等起到很大作用，能够使图像更加真实地被还原，并再现目标场景。

近年来，稀疏表示得到广泛关注，Rubinstein等人在文献“Doublesparsity:learningsparseddictionariesforsparsesignalapproximation”提出了稀疏K-SVD字典学习方法。该方法同时兼顾了Wavelet、Curvelet、DCT等字典的简单性，以及K-SVD、在线学习等方法训练获得的字典具有的灵活性，解决了K-SVD字典结构性差、高维数据处理难的问题。然而，此方法未利用各图像块之间存在的几何结构相似特性。

非局部均值(Non-LocalMeans,NLM)去噪方法是Buades等在2005年提出的，此方法提出了非局部思想的概念，充分利用了图像内部的自相似性。该方法在保持纹理和边缘信息方面有不错的效果，但是在图像的匀质区域会出现伪纹理的现象。同时，该方法是从像素点到像素点逐个进行处理的，导致计算复杂度很高，且未考虑图像块的能量集中性。

Dabov等人在非局部相似的理论基础上，于2007年提出了BM3D算法。BM3D算法是一种空域和频域相结合的去噪方式，不仅包含了非局部去噪的思想，而且利用了变换域滤波方法。该方法去噪后的图像具有较高的信噪比，且纹理信息保留较好，特别是对变化比较缓慢的边缘保持效果较好，是目前最理想的图像去噪算法之一。BM3D算法利用了图像变换的能量集中性，但相似块群的构造仅仅考虑了图像块的整体信息而并未利用图像的结构信息。

为克服BM3D算法的不足，基于图像局部结构相似特性，一些结构聚类的图像去噪算法被相继提出。典型代表是基于聚类的稀疏表示算法(Clustering-basedSparseRepresentation,CSR)，CSR结合了NLM和BM3D的思想，通过减小稀疏表示误差(SparseCodingNoise,SCN)以提高图像的重构质量。该方法的去噪效果优于BM3D，具有较高的信噪比，对纹理和边缘的保持也取得了非常好的效果。然而，由于通过不断的迭代更新完备主成分分析(PrincipalComponentAnalysis,PCA)字典和稀疏表示系数，对图像块进行稀疏表示，导致算法的时间复杂度大大增加。

发明内容

本发明目的在于结合上述去噪方法的优势并针对其不足，提出一种基于马尔科夫蒙特卡罗采样和自适应阈值低秩逼近的图像去噪方法，结合非局部相似的特性以及低秩逼近处理在处理高维数据时的优势，克服了字典训练的复杂度高的问题。本发明旨在通过蒙特卡罗采样产生图像块，并根据直方图多个统计特征，采用满足马尔科夫链的相似度判决函数选择图像相似块，以充分利用图像的结构信息。此外，在低秩条件下利用块之间的相似性对相似性矩阵进行秩最小化约束逼近，以减小相似块融合不准确的现象，有效地保留图像的细节、纹理等结构信息，从而提高去噪效果。

为实现上述目的，本发明的实现方法，包括如下步骤：

(1)输入一幅N行M列的含噪图像Y：

(2)设置参数：

观测图像Y∈R^m×n，样本点s_j，s_j的集合S_j，样本点数m₀，空间方差σ_s，期望图像加性噪声N，相似度评价函数Ψ；

(3)采用块直方图多个统计特征相似度判据进行马尔科夫蒙特卡罗采样，即对观测图像Y中的每个点s_j进行迭代操作，以获得相似匹配块簇；

3a)初始化参数：S_j＝φ，s_j＝s_j0，k＝1；

3b)对整幅图像进行采点：

按照高斯建议性分布 $Q({s^{\′}}_{jk}/s_{jk-1})=\frac{1}{2{\mathrm{\π\σ}}_s}\exp (-({\left({s^{\′}}_{jk}-s_{jk-1}\right)}^2)/\left(2{\mathrm{\σ}}_s\right)),$ 从位置s_j0开始，随机取一点s'_jk；

3c)通过s'_jk的接受概率决定是否接受s'_jk为序列集Ω中的点，接受概率为：

$\mathrm{\α}({s^{\′}}_{jk}/s_{jk-1})=min\{1,\frac{\mathrm{\Ψ}({s^{\′}}_{jk}/s_{j0})}{\mathrm{\Ψ}(s_{jk-1}/s_{j0})}\}$

其中，Ψ(s'_jk/s_j0)是用来度量采样点s'_jk与估计点s_j0近似程度的一个客观评价函数；

3d)本方法采用直方图多个统计特征的组合来度量图像块之间的相似度，其结构相似性评价函数为：

$\mathrm{\Ψ}({s^{\′}}_{jk}/s_{j0})=\underset{{\mathrm{\α}}_i\∈\mathrm{\ρ}}{\Π}\exp \[-{\left(F_{\mathrm{\ρ}}^{{s^{\′}}_{jk}}\right({\mathrm{\α}}_i)-{F_{\mathrm{\ρ}}}^{s_{j0}}({\mathrm{\α}}_i\left)\right)}^2\]$

3e)若Ψ(s'_jk/s_j0)值越大，则表示s'_jk与s_j0结构信息越相似，可通过判定来决定是否放入序列集Ω中；

3f)若k＜m₀+1，则转到3b)，否则转到(4)；

(4)基于序列集Ω，获得每一个图像块的相似性矩阵Y(s_j)，c_j(Y(s_j))＝{s_j0,s_j1,…,s_jk}，令c_j(Y(s_j))＝Y_j；

(5)依据图像的先验信息估计阈值矢量w，对序列集Ω的每一个图像块相似矩阵Y_j进行SVD分解(U_j,S_j,V_j)＝svd(Y_j)，其中U_j是左奇异矩阵，S_j是奇异值对角矩阵，V_j是右奇异矩阵；从而获得估计值

(6)对分解后的低秩结构进行秩最小化逼近，按如下公式进行秩最小化逼近：

${\hat{X}}_j=\underset{X_j}{\mathrm{argmin}}\frac{1}{{\mathrm{\σ}}_n^2}\left|\right|Y_j-X_j|{|}_F^2+|\left|X_j\right|{|}_{w,*}$

其中，X_j是待求的变量，||.||_w,*表示加权核范数，是等式右边的最优值。该等式的最优解是通过求解以下等式：

$X_j=S_w\left(Y_j\right)=U_j{S}_w\left({\mathrm{\Σ}}_j\right)V_j^T$

其中，S_w(.)表示阈值算子。因此，经秩最小化逼近后得到去噪后的样本集

(7)若含噪图像中的每个像素点s_j都处理完毕，转到(8)；否则，转到(3)；

(8)聚集去噪后的图像样本集以形成整幅干净图像X。

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

1.本发明结合了局部、非局部和低秩逼近技术，在去除噪声的同时更能实现和抑制为纹理，并能很好的保留待去噪图像中的纹理细节信息。

2.本发明由于对搜寻区域内的图像块利用直方图统计特性评价函数，并通过马尔科夫蒙特卡罗采样获取相似块簇，增强了相似块的准确性。

3.本发明由于使用相似性矩阵秩最小化的逼近问题，避免了非局部均值中权值趋同的问题，以及字典学习稀疏表示不够准确的问题。

4.本发明基于先验信息的特点，针对不同的奇异值采用不同的奇异值阈值计算方式，使得相似性矩阵秩逼近更准确，最大限度的平滑光滑区域，并保持纹理信息。

5.本发明由于是在空域中进行的，实现过程简单，且去噪效果良好。

附图说明

图1为本发明的总体流程图。

图2为不同去噪方法对图像Road的去噪效果图。

图3为不同去噪方法对图像Rail的去噪效果图。

具体的实施方式

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

步骤1，输入一幅N行M列的含噪图像Y：

步骤2，设置参数：

观测图像Y∈R^m×n，s_j表示样本点，S_j表示s_j的集合，m₀表示样本点数，σ_s表示空间方差，表示期望图像，N表示加性噪声，Ψ表示相似度评价函数；

步骤3，采用块直方图多个统计特征相似度判据进行马尔科夫蒙特卡罗采样，即对观测图像Y中的每个点s_j进行迭代操作，以获得相似匹配块簇；

3a)初始化参数：S_j＝φ，s_j＝s_j0，k＝1；

3b)对整幅图像进行采点：

按照高斯建议性分布 $Q({s^{\′}}_{jk}/s_{jk-1})=\frac{1}{2{\mathrm{\π\σ}}_s}\exp (-({\left({s^{\′}}_{jk}-s_{jk-1}\right)}^2)/(2{\mathrm{\σ}}_s\left)\right),$ 从位置s_j0开始，随机取一点s'_jk；

3c)通过s'_jk的接受概率决定是否接受s'_jk为序列集Ω中的点，接受概率为：

$\mathrm{\α}({s^{\′}}_{jk}/s_{jk-1})=min\{1,\frac{\mathrm{\Ψ}({s^{\′}}_{jk}/s_{j0})}{\mathrm{\Ψ}(s_{jk-1}/s_{j0})}\}$

其中，Ψ(s'_jk/s_j0)是用来度量采样点s'_jk与估计点s_j0近似程度的一个客观评价函数；

3d)本方法采用直方图多个统计特征的组合来度量图像块之间的相似度，其结构相似性评价函数为：

$\mathrm{\Ψ}({s^{\′}}_{jk}/s_{j0})=\underset{{\mathrm{\α}}_i\∈\mathrm{\ρ}}{\Π}\exp \[-{\left(F_{\mathrm{\ρ}}^{{s^{\′}}_{jk}}\right({\mathrm{\α}}_i)-{F_{\mathrm{\ρ}}}^{s_{j0}}({\mathrm{\α}}_i\left)\right)}^2\]$

3e)若Ψ(s'_jk/s_j0)值越大，则表示s'_jk与s_j0结构信息越相似，可通过判定来决定是否放入序列集Ω中；

第一步：随机产生服从[0,1]均匀分布的数l；

第二步：判决A为计算接受概率

第三步：判决B为根据判决条件判决是否接受s'_jk；若α＜l，转到第四步，否则转到第五步；

第四步：若s'_jk被接受，令S_j＝S_j∪{s'_jk}，k＝k+1；

第五步：若s'_jk被拒绝，令S_j＝S_j，s_jk＝s_jk-1，k＝k；

3f)若k＜m₀+1，则转到3b)，否则转到(4)；

步骤4，基于序列集Ω，获得每一个图像块的相似性矩阵Y(s_j)，c_j(Y(s_j))＝{s_j0,s_j1,…,s_jk}，令c_j(Y(s_j))＝Y_j；

步骤5，依据图像的先验信息估计阈值矢量w，对序列集Ω的每一个图像块相似矩阵Y_j进行SVD分解(U_j,S_j,V_j)＝svd(Y_j)，其中U_j是左奇异矩阵，S_j是奇异值对角矩阵，V_j是右奇异矩阵；从而获得估计值

步骤6，对分解后的低秩结构进行秩最小化逼近：

第一步：对奇异值对角矩阵S_j进行矢量操作，即diag(S_j)＝[λ_j1,...,λ_jm,...,λ_jr]，其中，λ_jm是Y_j相应的奇异值，m∈[1,r]，r为奇异值总数，diag(.)表示将对角化矩阵转成矢量的运算符号；

第二步：阈值w_jm被假定为与λ_jm(Y_j的第m个奇异值)有关，且阈值与λ_jm成反比，则奇异值的阈值计算如下：

$w_{jm}=c\frac{{\left({\mathrm{\σ}}_n\right)}^2}{{\mathrm{\λ}}_{jm}+\mathrm{\ϵ}}$

其中，c是阈值系数，c＞0且为常数，n是Y_j中相似块的数量，ε＝10^-16是避免被0分解的困扰，σ_n是含噪图像Y的噪声标准差。

第三步：利用下式获得相似性矩阵Y_j对应的去噪声矩阵的近似奇异值矩阵

$S_{{\mathrm{\Σ}}_j}=diag({\left({\mathrm{\λ}}_{j1}-{\mathrm{\τ}}_{j1}\right)}_{+},...,{\left({\mathrm{\λ}}_{jm}-{\mathrm{\τ}}_{jm}\right)}_{+},...,{\left({\mathrm{\λ}}_{jr}-{\mathrm{\τ}}_{jr}\right)}_{+})$

其中，(.)₊为取正运算符，当矩阵中的元素值小于0时，其元素取值为0；

第四步：根据第二步的结果重构出低秩的相似矩阵按如下公式进行秩最小化逼近：

${\hat{X}}_j=\underset{X_j}{\mathrm{argmin}}\frac{1}{{\mathrm{\σ}}_n^2}\left|\right|Y_j-X_j|{|}_F^2+|\left|X_j\right|{|}_{w,*}$

其中，X_j是待求的变量，||.||_w,*表示加权核范数，是等式右边的最优值。该等式的最优解是通过求解以下等式：

$X_j=S_w\left(Y_j\right)=U_j{S}_w\left({\mathrm{\Σ}}_j\right)V_j^T$

其中，S_w(.)表示阈值算子，经秩最小化逼近后得到去噪后的样本集

${\hat{X}}_j=U_j{S}_{{\mathrm{\Σ}}_j}{V}_j$

步骤7，若含噪图像中的每个像素点s_j都处理完毕，转到步骤8；否则，转到步骤3；

步骤8，聚集去噪后的图像样本集以形成整幅干净图像X。

本发明的效果可以通过以下实验进一步证实：

1、仿真条件

本方法所采用PC机的主频为2GHzIntel2处理器、8GBRAM，并在2011版Matlab编程环境下进行仿真。

2、评价指标

图像重构质量的评价可以从主观度量和客观度量两个方面进行。

3、仿真图像

仿真图像选择图2(a)Rail图像，图3(a)Road图像。

4、仿真对比方法：

对比方法1：Dabov等在文献“Imagedenoisingbysparse3-Dtransform-domaincollaborativefiltering”中提出BM3D方法。

对比方法2：WeishengDong在文献“ImageDeblurringwithLow-rankApproximationStructuredSparseRepresentation”中提出的同步稀疏低秩逼近方法，简称LASSC。

对比方法3：WeishengDong在文献“Sparsity-basedImageDenoisingviaDictionaryLearningandStructuralClustering”中提出的CSR方法。

5、仿真内容

利用现有的三种对比方法和本发明方法对Rail和Road两幅原始测试图像分别在噪声方差为10、20、30、40、50、80、100七种情况下进行了仿真实验，仿真结果如下表，其中sigma为噪声方差，SSIM为结构相似度。

上表显示了四幅图像在七种噪声水平下去噪后图像的PSNR，从中可以看出，本发明方法在七种噪声水平下的PSNR以及SSIM都高于其他三种方法，并且对于图像的纹理、细节结构信息保留更多，显示出了较好的去噪效果。从单幅图像的相同噪声水平下处理，本发明的PSNR及SSIM比其他方法都高，这足以充分地验证本发明方法较其他三种方法对自然图像去噪更具有效性。

Claims (3)

1.一种基于马尔科夫蒙特卡罗(Markov-chainMonteCarlo,MCMC)采样和阈值低秩逼近的图像去噪方法，包括如下步骤：

(1)输入一幅含噪图像Y：

Y＝X+N，其中，Y＝{Y(s),s∈S}表示含噪图像，X＝{X(s),s∈S}表示不含噪图像，N＝{N(s),s∈S}表示加性噪声，S是离散格点s的集合；

(2)设置参数：

观测图像Y∈R^m×n，s_j表示样本点，S_j表示s_j的集合，m₀表示样本点数，σ_s表示空间方差，表示期望图像，N表示加性噪声，Ψ表示相似度评价函数；

(3)采用块直方图多个统计特征相似度判据进行马尔科夫蒙特卡罗采样，即对观测图像Y中的每个点s_j进行迭代操作，以获得相似匹配块簇；

3a)初始化参数：S_j＝φ，s_j＝s_j0，k＝1；

3b)对整幅图像进行采点：

按照高斯建议性分布即第k个候选点必须在以第k-1个样本点为中心点的高斯分布意义下取得，从位置s_j0开始，随机取一点s'_jk；

3c)通过s'_jk的接受概率决定是否接受s'_jk为序列集Ω中的点，接受概率为：

其中，Ψ(s'_jk/s_j0)是度量采样点s'_jk与估计点s_j0近似程度的客观评价函数；

3d)本方法采用直方图多个统计特征的组合来度量图像块之间的相似度，其结构相似性评价函数为：

上式中，ρ＝{m,σ,R}表示像素点所对应图像块的直方图多个统计特征集合，其中：均值m表示纹理平均亮度，标准差s表示纹理平均对比度，平滑度R表示纹理亮度的相对平滑度；对于灰度一致的区域，其平滑度R为1；对于灰度级的值有较大差异的区域，R为0，起到保护边缘和细节信息的作用；m、s、R都对应上式中的α_i；表示s'_jk、s_j0点对应图像块在直方图统计特征ρ下，统计特征分量α_i对应的值；

3e)若Ψ(s'_jk/s_j0)值越大，则表示s'_jk与s_j0结构信息越相似，可通过判定来决定是否放入序列集Ω中；

3f)若k＜m₀+1，则转到3b)，否则转到(4)；

(4)基于序列集Ω，获得每一个图像块的相似性矩阵Y(s_j)，c_j(Y(s_j))＝{s_j0,s_j1,…,s_jk}，令c_j(Y(s_j))＝Y_j；

(5)依据图像的先验信息估计阈值矢量w，对序列集Ω的每一个图像块相似矩阵Y_j进行SVD分解(U_j,S_j,V_j)＝svd(Y_j)，其中U_j是左奇异矩阵，S_j是奇异值对角矩阵，V_j是右奇异矩阵；从而获得估计值

(6)对分解后的低秩结构进行秩最小化逼近，按如下公式进行秩最小化逼近：

其中，X_j是待求的变量，表示加权核范数，是等式右边的最优值；该等式的最优解是通过求解以下等式：

X_j＝S_w(Y_j)＝U_jS_w(Σ_j)V_j ^T

其中，S_w(.)表示阈值算子，因此，经秩最小化逼近后得到去噪后的样本集

(7)若含噪图像中的每个像素点s_j都处理完毕，转到(8)；否则，转到(3)；

(8)聚集去噪后的图像样本集以形成整幅干净图像X。

2.根据权利要求1所述的去噪方法，其中步骤3e)中所述的通过判定决定是否将所选点放入序列集Ω中，s'_jk的判定按照如下步骤进行：

(1)随机产生服从[0,1]均匀分布的数l；

(2)判决A为计算接受概率

(3)判决B为根据判决条件判决是否接受s'_jk；若α＜l，转到第四步，否则转到第五步；

(4)若s'_jk被接受，令S_j＝S_j∪{s'_jk}，k＝k+1；

(5)若s'_jk被拒绝，令S_j＝S_j，s_jk＝s_jk-1，k＝k。

3.根据权利要求1所述的去噪方法，其中步骤(6)中所述的对分解后的低秩结构进行秩最小化逼近，按如下步骤进行：

(1)对奇异值对角矩阵S_j进行矢量操作，即diag(S_j)＝[λ_j1,...,λ_jm,...,λ_jr]，其中，λ_jm是Y_j相应的奇异值，m∈[1,r]，r为奇异值总数，diag(.)表示将对角化矩阵转成矢量的运算符号；

(2)阈值w_jm被假定为与λ_jm(Y_j的第m个奇异值)有关，且阈值与λ_jm成反比，则奇异值的阈值计算如下：

其中，c是阈值系数，c＞0且为常数，n是Y_j中相似块的数量，ε＝10^-16是避免被0分解的困扰，σ_n是含噪图像Y的噪声标准差；

(3)利用下式获得相似性矩阵Y_j对应的去噪声矩阵的近似奇异值矩阵

其中，(.)₊为取正运算符，当矩阵中的元素值小于0时，其元素取值为0；

(4)根据第二步的结果重构出低秩的相似矩阵

。