CN103049892A - 基于相似块矩阵秩最小化的非局部图像去噪方法 - Google Patents

Abstract

本发明公开了一种基于相似块矩阵秩最小化的非局部图像去噪方法。其实现步骤为：(1)输入一幅N行M列的含噪图像；（2）估计含噪图像的噪声标准差，并根据标准差设置参数；（3）逐块的利用图像块的DCT特征计算图像块之间的欧氏距离；（4）对欧式距离由小到大排序，选取前k个对应的样本组成相似性矩阵；（5）对相似性矩阵进行秩最小化逼近，得到低秩矩阵；（6）对去噪后的图像块样本集进行聚集，得到去噪后图像；(7)判断是否达到迭代次数，若未达到迭代次数，转到步骤2到步骤6，否则输出结果。本发明具有重建结果边缘纹理结构信息保持好的优点，可用于医学影像、天文学影像、视频多媒体等领域的数字图像预处理。

Description

基于相似块矩阵秩最小化的非局部图像去噪方法

技术领域

本发明属于数字图像处理技术领域，具体地说是一种对含噪相似块矩阵低秩逼近无噪矩阵的去噪方法，可用于医学影像、天文学影像、视频多媒体等领域的数字图像预处理。

背景技术

图像已成为人类活动中最常用的信息载体，它们包含着物体的大量信息，成为人们获取外界原始信息的主要途径。然而图像在拍摄、采样、传输和存贮等过程中，经常会受到外界噪声的干扰和影响而使图像降质而不能真实地反映景物，并且图像预处理算法的好坏又直接关系到后续图像处理的效果，如图像分割、特征提取、目标识别等，所以在图像处理中，图像去噪起着重要作用。它可以有效地抑制各种噪声、保持边缘等细节信息，从而改善后续处理工作的质量。

目前存在的图像去噪方法有很多，主要有空域滤波和变换域滤波方法。空域滤波是在原图像上直接进行数据运算，对像素的灰度值进行处理。常见的空域滤波方法有Lee滤波、Frost滤波、非局部均值滤波NLM，基于稀疏表示和KSVD字典学习滤波KSVD的方法、基于块的最优滤波PLOW等。变换域滤波方法是对图像进行某种变换，将图像从空域转换到变换域，再对变换域中的变换系数进行处理，最后进行反变换将图像从变换域转换到空域来达到去除噪声的目的。将图像从空域转换到变换域的变换方法很多，如小波变换方法、三维块匹配去噪滤波BM3D法等。

由Buades等人在文献“A non local algorithm for image denoising.IEEE CVPR，2005,vol.2,pp.60-65.”中提出非局部均值方法，通过利用图像中模式结构的自相似来恢复原始图像，具体实现是计算图像中两个像素点为中心的图像块的欧氏距离来确定两点的相似度，由得到的相似度来确定该点对所要求的点的信息补偿程度，并设计负指数函数将两点间的欧氏距离转化为相似度权值，中心点的灰度值为邻域内像素点灰度值的加权平均。相比其它方法，它能准确地表征像素的特征，并能很好地解决图像边缘和线性体等细节的保留问题，但该方法过于强调邻域内像素点的非局部信息而忽略了中心像素点自身的局部信息，导致图像边缘细节信息丢失严重，图像平滑区域也存在过平滑现象。图像中存在的大量分形模式作为图像的先验信息，在越来越多的算法中得以应用并取得不错的效果。

由Elad等人在文献“Image denoising via spare and redundant representations overlearned dictionaries.IEEE TIP，2006,vol.15no.12,pp.3736-3745.”中提出的稀疏表示下的图像去噪方法，采用的是图像块在冗余字典上的稀疏近似来实现图像去噪，图像块在冗余字典上越稀疏，则逼近原图像块的信息就越准确。但该方法仍然存在的不足是，只是对图像块自身进行稀疏逼近，利用的只是图像块自身的局部信息，而忽视了邻域图像块的非局部信息，使得图像的平滑区域出现伪纹理的现象并且纹理区域去噪效果也不够理想。

由Dabov等人在文献“Image denoising by sparse3-D transform-domain collaborativefiltering.IEEE Transactions on Image Processiag，2007,vol.16no.8,pp.2080-2095.”中提出来的BM3D去噪方法。该方法结合了三维变换、非局部和维纳滤波，并采用了基于两步的去噪方式。该方法的主要步骤：第一步，对噪声图像进行基于三维变换的硬阈值；第二步，将硬阈值的去噪结果作为初始估计再对原噪声图像进行基于三维变换的维纳滤波。但该方法仍然存在的不足有：该方法中第一步的去噪方式采用的是硬阈值，通过这样的方式易造成图像边缘和纹理细节的丢失；该方法中第一步的去噪结果是第二步的初始估计，如果初始估计不准确，那么直接地影响了第二步的去噪结果。

由Priyam Chatterjee在文献“Patch-Based Near-Optimal Image Denoising.IEEETransactions on Image Processing，2012,vol.21,NO.4,pp.1635-1649”中提出的PLOW方法，结合了基于结构聚类的PCA变换、非局部和维纳滤波。该方法首先对图像进行逐像素分块，然后再对图像块进行聚类并计算每一类的均值和方差，最后在非局部搜索窗内找相似块，并结合所属类的均值和方差利用设计的维纳滤波器对图像块去噪。该方法的不足是：利用的是每一类的均值和方差对图像块进行去噪，由于有些相似块不够相似，使得均值和方差的估计都不够准确，导致去噪效果不够理想。

发明内容

本发明的目的在于针对上述去噪方法的不足，提出一种基于相似块矩阵秩最小化的非局部图像去噪方法，以在低秩的条件下利用块之间的相似性对相似性矩阵进行秩最小化约束逼近，减小相似块融合不准确的现象，提高去噪效果。

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

（1）输入一幅N行M列的含噪图像X^k，其中迭代次数k=0；

（2）对该含噪图像X^k进行小波分解，得到的第一层高频系数CI；

（3）估计含噪图像X^k的噪声标准差σ_n：

${\mathrm{\σ}}_n=\frac{\mathrm{median}\left(\right|\mathrm{CI}\left|\right)}{0.6745}$

其中|·|是取绝对值操作，median(·)是取中值操作；

（4）根据噪声标准差的大小设置迭代次数t、图像块的边长l、搜索窗的边长s、相似块总数k，并且逐像素地提取图像块样本，将所有的样本组成一个图像样本集Y；

（5）以像素点i∈[1，N×M]为中心选取一个边长为l的图像块y_i，对边长为s的搜索框中的任一点j∈[1，s²]，选取一个边长为l的图像块y_j，利用下式计算图像块之间的欧式距离d(y_i,y_j)

$d(y_i,y_j)={\left|\right|\mathrm{\γ}\left(T\left(y_i\right)\right)-\mathrm{\γ}\left(T\left(y_j\right)\right)\left|\right|}_2^2$

其中T是离散余弦线性变换，γ是阈值系数为1.9σ_n的阈值算子，

表示二范数；

（6）对得到的欧式距离由小到大排序，选取前k个对应的样本作为y_i的相似样本，得到图像块对应的相似性矩阵M_i；

（7）对样本集Y中每一个图像块获得的相似性矩阵M_i利用下式进行秩最小化逼近，得到已去噪的样本集

$D_i^{*}=\underset{D_i}{\arg }\min \frac{1}{2}{\left|\right|M_i-D_i\left|\right|}_2^2+{\left|\right|D_i\left|\right|}_{*}^W$

其中，D_i是待求的变量，表示加权核范数，

是等式右边的最优值。该等式的最优解是通过求解以下等式：

$D_i^{*}=T_w\left(M_i\right)=U_i{T}_w\left(S_i\right)V_i^T$

其中，T_w()表示阈值算子，U_i是左奇异矩阵，S_i是奇异值对角矩阵，V_i是右奇异矩阵；

（8）对自适应去噪后的图像样本集

进行聚集，并对图像块与图像块之间重叠的像素点取均值，得到去噪后图像

（9）判断是否达到迭代次数t，如果t>k，则k＝k+1，且利用公式(3)的残差补偿得到下一次迭代的噪声图像X^k，把噪声图像X^k作为输入图像重复步骤（2）到步骤（8）；否则终止迭代，输出最终的去噪结果

$X^k={\hat{X}}^k+\mathrm{\δ}(X^0-{\hat{X}}^k),$

其中，X⁰是原始的含噪图像，δ是比例因子。

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

1．本发明是一种结合了局部、非局部和低秩逼近技术，在去噪方式上由于采用了一种噪声重新估计的迭代技术，在去除噪声的同时更能实现抑制伪纹理，并能很好的保留待去噪图像中的纹理细节信息。

2．本发明由于对搜寻区域内的的图像块利用提取的特征进行相似块的获取，增强了相似块的准确性。

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

4．本发明基于块的结构信息调节相似性矩阵秩的范围，最大限度的平滑光滑区域，并保持纹理信息；

5．本发明由于是在空域中进行，实现过程简单，且达到较好的效果。

附图说明

图1是本发明的实现流程图；

图2是本发明实验采用的四幅无噪自然图像；

图3是本发明实验输入的Fingerprint、Monarch、Barbara、Lena四幅有噪图像：

图4是用现有的三种方法和本发明方法对Fingerprint图像的去噪效果比较图；

图5是用现有的三种方法和本发明方法对Barbara图像的去噪效果比较图；

图6是图5的局部放大图。

具体实施方式

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

步骤1，输入一幅N行M列的含噪图像X^k，其中迭代次数k=0，并对该含噪图像X^k进行小波分解，得到的第一层高频系数CI。

步骤2，估计含噪图像X^k的噪声标准差σ_n：

${\mathrm{\σ}}_n=\frac{\mathrm{median}\left(\right|\mathrm{CI}\left|\right)}{0.6745}$

其中|·|是取绝对值操作，median(·)是取中值操作。

步骤3，根据噪声标准差的大小设置参数：

参数包括迭代次数t、图像块的边长l、搜索窗的边长s、相似块总数k，和比例因子δ：

如果噪声标准差σ_n≤30，则l＝7，s＝21，k=66，t=9，δ=0.031；

如果噪声标准差σ_n>30，则l＝9，s＝31，k=88，t＝10，δ=0.031。

步骤4，构建图像块样本集Y

首先从含噪图像X^k中逐像素地提取边长为l的图像块；然后将图像块转换成向量样本；最后将所有的样本组成一个图像块样本集Y。

步骤5，以像素点i∈[1，N×M]为中心选取一个边长为l的图像块y_i，对边长为s的搜索框中的任一点j∈[1，s²]，选取一个边长为l的图像块y_j，利用下式计算图像块之间的欧式距离d(y_i,y_j)：

$d(y_i,y_j)={\left|\right|\mathrm{\γ}\left(T\left(y_i\right)\right)-\mathrm{\γ}\left(T\left(y_j\right)\right)\left|\right|}_2^2$

其中，T是离散余弦线性变换，γ是阈值系数为1.9σ_n的阈值算子，

表示二范数。

步骤6，对得到的欧式距离d(y_i,y_j)由小到大排序，选取前k个对应的样本作为y_i的相似样本，得到图像块对应的相似性矩阵M_i。

步骤7，对样本集Y中每一个图像块获得的相似性矩阵M_i采用以下等式进行秩最小化逼近，得到已去噪的样本集

$D_i^{*}=\underset{D_i}{\arg }\min \frac{1}{2}{\left|\right|M_i-D_i\left|\right|}_2^2+{\left|\right|D_i\left|\right|}_{*}^W$

其中，D_i是待求的变量，表示加权核范数，是等式右边的最优值。

该等式的最优解是通过求解以下等式：

$D_i^{*}=T_w\left(M_i\right)=U_i{T}_w\left(S_i\right)V_i^T$

其中，T_w()表示阈值算子，U_i是左奇异矩阵，S_i是奇异值对角矩阵，V_i是右奇异矩阵；

7a)对步骤6获得每个图像块的相似性矩阵M_i进行SVD分解，即M_i=U_iS_iV_i，其中，U_i是左奇异矩阵，S_i是奇异值对角矩阵，V_i是右奇异矩阵；

7b）对奇异值对角矩阵S_i进行矢量操作，即diag(S_i)=[λ_i1,...,λ_im,...,λ_ir]，其中，λ_im是M_i相应的奇异值，m∈[1,r]，r为奇异值总数，diag()表示将对角化矩阵转成矢量的运算符号；

7c）计算奇异值的阈值：

${\mathrm{\τ}}_{i,m}=c\frac{{\left({\mathrm{\σ}}_n\right)}^2}{{\mathrm{\λ}}_{i,m}},$

其中c是阈值系数，取值为2.1σ_n，σ_n是含噪图像X^k的噪声标准差，λ_i，m是第i个像素点对应的相似性矩阵M_i的第m个奇异值，τ_i，m表示M_i的第m个奇异值阈值；

7d）利用下式获得相似性矩阵M_i对应的无噪矩阵的近似奇异值矩阵

$T_{S_i}=\mathrm{diag}({({\mathrm{\λ}}_{i,1}-{\mathrm{\τ}}_{i,1})}_{+},...,{({\mathrm{\λ}}_{i,m}-{\mathrm{\τ}}_{i,m})}_{+},...,{({\mathrm{\λ}}_{i,r}-{\mathrm{\τ}}_{i,r})}_{+}),$

其中，diag(·)为将矢量转成对角化矩阵的操作，(·)₊为取正运算符，当矩阵中的元素数值小于0时，对其元素取值为0；

7e）根据步骤(7d)的结果重构出低秩的相似矩阵

${\hat{M}}_i=U_i{T}_{S_i}{V}_i,$

7f)对图像样本集Y中的所有图像块样本y_i，i＝1,2,…,N×M，逐个进行步骤（7a）到（7e）,得到已去噪的样本集

步骤8，对自适应去噪后的图像样本集

进行聚集，并对图像块与图像块之间重叠的像素点取均值，得到去噪后图像

步骤9，判断是否达到迭代次数t，如果t>k，则k＝k+1，且利用下式残差补偿得到下一次迭代的噪声图像X^k，把噪声图像X^k作为输入图像重复步骤2到步骤8；否则终止迭代，输出最终的去噪结果

$X^k={\hat{X}}^k+\mathrm{\δ}(X^0-{\hat{X}}^k),$

其中，X⁰是原始的含噪图像，δ是比例因子。

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

1．仿真条件：

在CPU为pentium(R)4处理器：主频1.86GHZ，内存2G，操作***：WINDOWSXP SP3，仿真平台：Matlab7.10平台上进行。

2．评价指标：

图像去噪效果的评价分为主观和客观两个方面。在主观上评价一幅图像去噪效果的优劣主要是通过人眼的视觉特性来衡量，图像质量好，感觉清晰则去噪效果好，反之效果则差；在客观上评价一幅图像的去噪效果，本发明采用峰值信噪比PSNR来衡量。令真实图像为Y，去噪图像为Y′，Y(t)为真实图像Y中第t个像素的幅值，Y′(t)为去噪图像Y′中第t个像素的幅值，Ω为图像中的像素集合，令Y_max=max{Y(t),t∈Ω}，则峰值信噪比按如下公式计算：

$\mathrm{PSNR}=10\log \left[\frac{Y_{\max }^2}{\frac{1}{\left|\mathrm{\Ω}\right|}\underset{t\∈\mathrm{\Ω}}{\mathrm{\Σ}}{[Y^{\′}\left(t\right)-Y\left(t\right)]}^2}\right]$

3．仿真图像：

仿真图像选择图2所示的四幅原始测试图像，其中，图2（a）为Fingerprint图像，图2（b）为Monarch图像，图2（c）为Barbara图像，图2（d）为Lena图像，图3为图2所示的四幅原始测试图像加入标准差为30的高斯随机白噪声的噪声图像，其中图3（a）为含噪Fingerprint图像，图3（b）为含噪Monarch图像，图3(c)为含噪Barbara图像，图3(d)为含噪Lena图像。四幅图像均为灰度图像，灰度级均为256。

4.仿真对比方法：

对比方法1：Buades等人在文献“A non local algorithm for image denoising.IEEECVPR，2005,vol.2,pp.60-65.”中提出非局部均值方法，简称NLM。

对比方法2：Elad等人在文献“Image denoising via spare and redundantrepresentations over learned dictionaries.IEEE TIP，2006,vol.15no.12,pp.3736-3745.”中提出的稀疏表示下的图像去噪方法，简称KSVD。

对比方法3：Dabov等在文献“Image denoising by sparse3-D transform-domaincollaborative filtering.IEEE TIP，2007,vol.16no.8,pp.2080-2095.”中提出的基于三维块匹配去噪方法，简称BM3D。

对比方法4：J.Mairal,和F.Bach在文献“Non-local sparse models for imagerestoration.IEEE 12th ICCV，2009，pp.2272-2279.”中提出的同步稀疏方法，简称LSSC。

对比方法5：Priyam Chatterjee和Peyman Milanfar在文献“Clustering-baseddenoising with locally learned dictionaries.IEEE TIP，2012,vol.21,no.4,pp.1635-1649.”中提出的基于聚类的局部字典学习去噪方法，简称KLLD。

5.仿真内容

仿真1，利用现有的五种对比方法和本发明方法对图2中的四幅图分别在噪声方差为5、10、20、25和50五种情况下进行了去噪仿真，仿真结果如表1。

表1 四幅图像在六种情况下的峰值信噪比结果 其中sigma为噪声方差。

从表1的四幅图像在五种噪声水平下的去噪后的PSNR可以看出，本发明方法在五种噪声水平下的PSNR大都高于其他五种方法，尤其是对于纹理内容较多的Monarch、Fingerprint和Barbara图像的去噪性能较好。从单幅图像的相同噪声水平下处理，本发明方法的PSNR大部分比其它方法高，所以可以充分地验证本发明方法较其它五种方法对自然图像去噪更具有效性。

仿真2，利用现有的三种方法和本发明方法对图3(a)进行去噪仿真，结果如图4所示，其中，图4(a)、图4(b)、图4(c)和图4(d)分别是利用BM3D方法、KSVD方法、NLM方法和本发明方法对图3(a)进行去噪的结果；

仿真3，利用现有的三种方法和本发明方法对图3(c)进行去噪仿真，结果如图5所示，其中，图5(a)、图5(b)、图5(c)和图5(d)分别是利用BM3D方法、KSVD方法、NLM方法和本发明方法对图3(c)进行去噪的结果图。

图6是图5的局部放大图，可以看出本发明的噪声边缘和细节保持能力优于NLM方法和KSVD方法，轮廓纹理保持效果优于BM3D方法。

综上所述，本发明的去噪效果优于现有的三种去噪方法。

Claims (2)

1.一种基于相似块矩阵秩最小化的非局部图像去噪方法，包括如下步骤：

（1）输入一幅N行M列的含噪图像X^k，其中k=0；

（2）对该含噪图像X^k进行小波分解，得到的第一层高频系数CI；

（3）估计含噪图像X^k的噪声标准差σ_n：

${\mathrm{\σ}}_n=\frac{\mathrm{median}\left(\right|\mathrm{CI}\left|\right)}{0.6745}$

其中|·|是取绝对值操作，median(·)是取中值操作；

（4）根据噪声标准差的大小设置迭代次数t、图像块的边长l、搜索窗的边长s、相似块总数k，其中本发明实例中迭代次数t取值为9，比例因子δ取值为0.031，如果噪声标准差σ_n≤30，则l＝7，s＝21，k=66；如果σ>30，则l＝9，s=31，k=88，并且逐像素地提取图像块样本，将所有的样本组成一个图像样本集Y；

（5）以像素点i∈[1，N×M]为中心选取一个边长为l的图像块y_i，对边长为s的搜索框中的任一点j∈[1，s²]，选取一个边长为l的图像块y_j，利用下式计算图像块之间的欧式距离d(y_i,y_j)：

$d(y_i,y_j)={\left|\right|\mathrm{\γ}\left(T\left(y_i\right)\right)-\mathrm{\γ}\left(T\left(y_j\right)\right)\left|\right|}_2^2$

其中，T是离散余弦线性变换，γ是阈值系数为1.9σ_n的阈值算子，表示二范数；

（6）对得到的欧式距离由小到大排序，选取前k个对应的样本作为y_i的相似样本，得到图像块对应的相似性矩阵M_i；

（7）对样本集Y中每一个图像块获得的相似性矩阵M_i利用下式进行秩最小化逼近，得到已去噪的样本集

$D_i^{*}=\underset{D_i}{\arg }\min \frac{1}{2}{\left|\right|M_i-D_i\left|\right|}_2^2+{\left|\right|D_i\left|\right|}_{*}^W$

其中，D_i是待求的变量，

表示加权核范数，

是等式右边的最优值。该等式的最优解是通过求解以下等式：

$D_i^{*}=T_w\left(M_i\right)=U_i{T}_w\left(S_i\right)V_i^T$

其中，T_w(．)表示阈值算子，U_i是左奇异矩阵，S_i是奇异值对角矩阵，V_i是右奇异矩阵；

（8）对自适应去噪后的图像样本集

进行聚集，并对图像块与图像块之间重叠的像素点取均值，得到去噪后图像

（9）判断是否达到迭代次数t，如果t>k，则k＝k+1，且利用下式的残差补偿得到下一次迭代的噪声图像X^k，把噪声图像X^k作为输入图像重复步骤（2）到步骤（8）；否则终止迭代，输出最终的去噪结果

$X^k={\hat{X}}^k+\mathrm{\δ}(X^0-{\hat{X}}^k)$

其中，X⁰是原始的含噪图像，δ是比例因子。

2.根据权利要求1所述的去噪方法,其中步骤（7）所述的对样本集Y中每一个图像块获得的相似性矩阵M_i进行秩最小化逼近，按如下步骤进行：

7a)对步骤6获得每个图像块的相似性矩阵M_i进行SVD分解，即M_i=U_iS_iV_i，其中，U_i是左奇异矩阵，S_i是奇异值对角矩阵，V_i是右奇异矩阵；

7b）对奇异值对角矩阵S_i进行矢量操作，即diag(S_i)=[λ_i1,...,λ_im,...,λ_ir]，其中，λ_im是M_i相应的奇异值，m∈[1,r]，r为奇异值总数，diag()表示将对角化矩阵转成矢量的运算符号；

7c）计算奇异值的阈值：

${\mathrm{\τ}}_{i,m}=c\frac{{\left({\mathrm{\σ}}_n\right)}^2}{{\mathrm{\λ}}_{i,m}}$

其中，c是阈值系数，取值为2.1σ_n，σ_n是含噪图像X^k的噪声标准差，λ_im是第i个像素点对应的相似性矩阵M_i的第m个奇异值，τ_i，m表示M_i的第m个奇异值阈值；

7d）利用下式获得相似性矩阵M_i对应的无噪矩阵的近似奇异值矩阵

$T_{S_i}=\mathrm{diag}({({\mathrm{\λ}}_{i,1}-{\mathrm{\τ}}_{i,1})}_{+},...,{({\mathrm{\λ}}_{i,m}-{\mathrm{\τ}}_{i,m})}_{+},...,{({\mathrm{\λ}}_{i,r}-{\mathrm{\τ}}_{i,r})}_{+}),$

其中，diag(·)为将矢量转成对角化矩阵的操作，(·)+为取正运算符，当矩阵中的元素数值小于0时，对其元素取值为0；

7e）根据步骤(7d)的结果重构出低秩的相似矩阵

${\hat{M}}_i=U_i{T}_{S_i}{V}_i.$

Images