CN102184526A - 基于字典学习和块匹配的自然图像去噪方法 - Google Patents

Abstract

本发明公开了一种基于字典学习和块匹配的自然图像去噪方法，主要克服现有自然图像去噪中纹理细节易丢失和同质区域不平滑的问题。其实现过程是：(1)设定去噪目标函数，输入含噪图像z(x)；(2)令原始图像等于含噪图像，即y(x)＝z(x)，令字典D为冗余DCT字典；(3)采用KSVD算法对字典D的原子和相应系数矩阵α_ij进行更新；(4)用BM3D算法对含噪图像z(x)去噪，获得初步去噪结果

(5)将更新后的D、α_ij和

代入原始图像的估计公式，得到含噪图像z(x)的去噪结果

本发明相对于现有的经典去噪方法能够更好地去除噪声，在平滑同质区域的同时兼顾保留图像的纹理、轮廓、边缘细节信息，可用于自然图像的去噪。

Description

基于字典学习和块匹配的自然图像去噪方法

技术领域

本发明属于图像处理技术领域，具体地说是一种稀疏表示和字典学习的去噪方法，可用于图像处理，模式识别和生物医学等领域。

背景技术

图像去噪的目的是在去除图像噪声的同时保留图像特征信息，像纹理，边缘，轮廓和点状目标等。自然图像的噪声具有加性背景，去噪可以从空域和变换域两方面来进行，目前效果较好的空域滤波方法有非局部均值滤波方法、稀疏表示下的图像去噪等，较好的变换域的滤波方法有三维块匹配方法BM3D等。

非局部均值方法，通过计算图像中两个像素点邻域的相似度来确定该点对所要求的点的信息补偿程度，中心点的灰度值为邻域内像素点灰度值的加权平均。它虽然能很好的解决图像边缘和线性体的保留问题，但仍存在平滑区域过平滑现象。

图像稀疏表示下的去噪方法，是新近提出的一种空域图像去噪新方法，它采用图像在冗余字典上的稀疏近似来实现噪声去除，如DCT字典去噪和KSVD字典学习去噪。然而，DCT字典不能有效逼近原图的边缘与细节信息，而KSVD字典学习的误差控制方法粗糙，易造成原图中部分纹理细节的丢失，而且这种方法忽视了图像的全局结构。

三维块匹配方法BM3D，是一种可有效滤除高斯噪声的去噪方法，该方法不仅利用图像的结构信息，而且结合了变换域的阈值方法，是一种非局部、多尺度的自适应非参数化滤波技术，它虽然能较好保留图像边缘与纹理等信息，但是对某些细小的条纹状纹理去噪时会出现失真，且在噪声较大的情况下，去噪结果中会出现划痕。

发明内容

本发明的目的在于克服上述已有技术的不足，提出一种基于字典学习和块匹配的自然图像去噪方法，以在提高同质区域平滑性的同时兼顾保留图像的纹理、轮廓、边缘细节信息，提高图像去噪效果。

实现本发明目的的技术思路是在稀疏表示和字典学习图像去噪的框架下，将空域字典学习和非局部、多尺度的BM3D方法结合，同时考虑图像的局部和全局信息，设计了基于稀疏表示误差、空域估计误差、非局部多尺度估计误差，以及稀疏系数零范数最小化的目标函数，实现对自然图像去噪，其具体步骤包括如下：

(1)设含噪图像z(x)＝y(x)+n(x)去噪后的原始图像的估计

为：

$\hat{y}\left(x\right)=\arg \underset{D,{\mathrm{\α}}_{\mathrm{ij}},y\left(x\right)}{\min }\mathrm{\λ}{\left|\right|y\left(x\right)-z\left(x\right)\left|\right|}_2^2+\underset{\mathrm{ij}}{\mathrm{\Σ}}\mathrm{\μ}{\left|\right|{\mathrm{\α}}_{\mathrm{ij}}\left|\right|}_0+\underset{\mathrm{ij}}{\mathrm{\Σ}}{\left|\right|{\mathrm{D\α}}_{\mathrm{ij}}-R_{\mathrm{ij}}y\left(x\right)\left|\right|}_2^2+v\left|\right|y\left(x\right)-{\hat{y}}_1\left(x\right){\left|\right|}_2^2$

其中，y(x)为原始图像，n(x)为标准差为σ的高斯白噪声，x∈X，X为图像像素的坐标集合，λ为langrage乘子，且λ＝30/σ，代表取二范数平方，μ为残差控制因子，α_ij为图像块R_ijy(x)的稀疏表示系数，||||₀代表取零范数，D为过完备稀疏表示字典，R_ij代表从图像中取出

大小的图像块的取块操作，i，j为取出的图像块在原图像中的坐标位置，v为权重参数，它的取值范围是(78/σ)²～(1049/σ)²，

为BM3D方法对含噪图像z(x)的初步去噪结果；

(2)令原始图像等于含噪图像，即y(x)＝z(x)，令过完备稀疏表示字典D为冗余DCT字典，且D∈R^n×K，R^n×K代表n×K大小的实数矩阵，取n＝64，K＝256；

(3)采用KSVD字典学习方法更新过完备稀疏表示字典D和稀疏表示系数α_ij，此步骤迭代J次，当σ≤5时，J＝5，当σ＞5时，J＝10；

(4)用BM3D方法对含噪图像z(x)去噪，得到初步去噪结果

(5)将更新后的过完备稀疏表示字典D、更新后的稀疏表示系数α_ij和初步去噪结果

代入步骤(1)设定的原始图像的估计公式，得到原始图像的估计

$\hat{y}\left(x\right)=[(\mathrm{\λ}+v)I+\underset{\mathrm{ij}}{\mathrm{\Σ}}{R_{\mathrm{ij}}^T{R}_{\mathrm{ij}}]}^{-1}[\mathrm{\λz}\left(x\right)+\underset{\mathrm{ij}}{\mathrm{\Σ}}{R}_{\mathrm{ij}}^T{\mathrm{D\α}}_{\mathrm{ij}}+{v\hat{y}}_1\left(x\right)]$

其中，I为图像大小的单位矩阵，

为R_ij的转置。

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

1、本发明由于利用图像的多尺度变换域和空域信息，并基于图像在变换域和空域中过完备冗余字典下的稀疏表示，设计了含噪图像的去噪目标函数，因而在有效去除噪声的同时，兼顾保留了图像的纹理、轮廓、边缘细节信息，提高了图像去噪效果。

2、由于本模型也可以处理非高斯加性噪声，乘性斑点噪声也可以转化成加性高斯白噪声进行去噪，故本发明可以方便的推广到非高斯加性和斑点噪声图像去噪。

附图说明

图1是本发明的流程图；

图2是本发明实验输入的原始图像；

图3是对图2加了噪声标准差为20的高斯白噪声后的含噪图像；

图4是用现有KSVD字典学习方法对图3去噪的结果图；

图5是用现有BM3D方法对图3去噪的结果图；

图6是用本发明对图3去噪的结果图。

具体实施方式

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

步骤1：设计含噪图像z(x)＝y(x)+n(x)的去噪目标函数，最小化该目标函数得到原始图像的估计。

自然图像的噪声具有加性背景，常见的加性噪声为高斯白噪声，它服从正态分布N～(0，σ²)，其中σ为噪声标准差，根据高斯白噪声的性质，设原始图像的估计

为：

$\hat{y}\left(x\right)=\arg \underset{D,{\mathrm{\α}}_{\mathrm{ij}},y\left(x\right)}{\min }\mathrm{\λ}{\left|\right|y\left(x\right)-z\left(x\right)\left|\right|}_2^2+\underset{\mathrm{ij}}{\mathrm{\Σ}}\mathrm{\μ}{\left|\right|{\mathrm{\α}}_{\mathrm{ij}}\left|\right|}_0+\underset{\mathrm{ij}}{\mathrm{\Σ}}{\left|\right|{\mathrm{D\α}}_{\mathrm{ij}}-R_{\mathrm{ij}}y\left(x\right)\left|\right|}_2^2+v\left|\right|y\left(x\right)-{\hat{y}}_1\left(x\right){\left|\right|}_2^2$

其中，y(x)为原始图像，n(x)为标准差为σ的高斯白噪声，x∈X，X为图像像素的坐标集合，λ为langrage乘子，且λ＝30/σ，代表取二范数平方，μ为残差控制因子，α_ij为图像块R_ijy(x)的稀疏表示系数，||||₀代表取零范数，D为过完备稀疏表示字典，R_ij代表从图像中取出

大小的图像块的取块操作，i，j为取出的图像块在原图像中的坐标位置，v为权重参数，它的取值范围是(78/σ)²～(1049/σ)²，在实验中根据不同噪声方差选择最佳参数，且v随σ的增大而减小，为BM3D方法对含噪图像z(x)的初步去噪结果。

步骤2：令原始图像等于含噪图像，即y(x)＝z(x)，令过完备稀疏表示字典D为冗余DCT字典，且D∈R^n×K，R^n×K代表n×K大小的实数矩阵，取n＝64，K＝256。由于原始图像的估计公式中有四个未知量原始图像y(x)、过完备稀疏表示字典D、稀疏表示系数α_ij和初步去噪结果

因此要先初始化原始图像y(x)和过完备稀疏表示字典D再计算稀疏表示系数α_ij。

步骤3：在已知原始图像y(x)、过完备稀疏表示字典D和稀疏表示系数式_ij的情况下，采用KSVD字典学习方法更新过完备稀疏表示字典D和稀疏表示系数α_ij，此步骤迭代J次，当σ≤5时，J＝5，当σ＞5时，J＝10。

3.1)将原始图像y(x)分解成为

(本文实验中取为8)大小的图像块y_ij(x)＝R_ijy(x)，用追踪算法计算图像块y_ij(x)的稀疏表示系数α_ij，即

${\mathrm{\α}}_{\mathrm{ij}}=\arg \underset{{\mathrm{\α}}_{\mathrm{ij}}}{\min }{\left|\right|{\mathrm{\α}}_{\mathrm{ij}}\left|\right|}_{0}s.t{\left|\right|{\mathrm{D\α}}_{\mathrm{ij}}-y_{\mathrm{ij}}\left(x\right)\left|\right|}_2^2\≤{\left(1.15\mathrm{\σ}\right)}^2.$

3.2)对于字典D的任一原子d_l，l＝1，2，...，K，找出用到原子d_l的那些图像块的坐标位置集合，即w_l＝{(i，j)|α_ij(l)≠0}，其中α_ij(l)表示稀疏系数α_ij的第l个元素；

3.3)对于

计算残差

其中d_m为字典D的第m个原子，α_ij(m)表示稀疏系数α_ij的第m个元素；

3.4)令组成的集合为E_l，即

3.5)对所述E_l进行奇异值分解，得到E_l＝UΔV^T，其中U为左奇异矩阵，Δ为奇异值矩阵，V^T右奇异矩阵的转置，用U中第一列更新原子d_l，用V^T的第一行乘以Δ(1，1)以对α_ij(l)进行更新，其中

步骤4：用BM3D方法对含噪图像z(x)去噪，得到初步去噪结果

4.1)将含噪图像z(x)，分解成8×8大小的图像块，依次取出一块图像块作为当前参考块，设为令所有图像块为侯选块，设为Z_x；

4.2)用块匹配法计算当前参考块

和候选块Z_x的距离

$d(Z_{x_R},Z_x)=\frac{{\left|\right|\mathrm{\γ}(T_{2D}\left(Z_{x_R}\right)-\mathrm{\γ}(T_{2D}\left(Z_x\right)\left|\right|}_2^2}{8^2}$

其中，T_2D表示二维小波线性变换，γ代表硬阈值算子，设硬阈值为0，

代表取二范数平方，取出满足

的侯选块作为与参考块

相似的块，将这些相似块堆积成三维形式，表示为其中

表示这些相似块的坐标集合；

4.3)在三维域对所述的

作硬阈值联合滤波，得到

的去噪结果

并将去噪结果

拆分二维形式，表示为

x₁∈S₁，其中，T_3D表示三维线性变换，γ′表示硬阈值算子，取硬阈值为2.7σ，表示逆三维线性变换；

4.4)重复步骤5.2)-5.3)对所有参考块进行处理，得到所有参考块对应相似块的估计

x₁∈S₁，加权平均

得到原始图像的基本估计

$\∀x\∈X$

其中，

为位于坐标x₁∈X的图像块的特征函数矩阵，为

的权值，

式中N′为硬阈值去噪后

中非零系数的个数；

4.5)将求得的

分解成8×8大小的图像块，依次取出一块图像块作为当前参考块，设为

令所有图像块为侯选块，设为

4.6)用块匹配法计算当前参考块

和候选块的距离

取出满足的侯选块作为与参考块

相似的块，将这些相似块堆积成三维形式，表示为

其中

表示这些相似块的坐标集合；从含噪图像z(x)中找出位于坐标集合S₂的那些图像块，并堆积成三维形式，表示为

4.7)在三维域对所述的

进行Wiener联合滤波，得到

的去噪结果

其中，T′_3D表示在三维域作2D-DCT线性变换；

代表Wiener滤波，即Wiener经验收缩系数，且

||²为求模的平方，

表示逆三维线性变换；将去噪结果

拆分成二维形式，表示为

x₂∈S₂；

4.8)重复步骤5.6)-5.7)对所有参考块进行处理，得到所有参考块对应相似块的估计

x₂∈S₂，

加权平均

得含噪图像z(x)的初步去噪结果 ${\hat{y}}_1\left(x\right):$

$\∀x\∈X$

其中，

为位于坐标x₂∈X的图像块的特征函数矩阵，为

的权值，表示为

式中

表示对二范数平方取逆。

步骤5：将更新后的过完备稀疏表示字典D、更新后的稀疏表示系数α_ij和初步去噪结果

代入步骤(1)设定的原始图像的估计公式，得到原始图像的估计

$\hat{y}\left(x\right)=[(\mathrm{\λ}+v)I+\underset{\mathrm{ij}}{\mathrm{\Σ}}{R_{\mathrm{ij}}^T{R}_{\mathrm{ij}}]}^{-1}[\mathrm{\λz}\left(x\right)+\underset{\mathrm{ij}}{\mathrm{\Σ}}{R}_{\mathrm{ij}}^T{\mathrm{D\α}}_{\mathrm{ij}}+{v\hat{y}}_1\left(x\right)]$

其中，I为图像大小的单位矩阵，

为R_ij的转置。

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

一.实验条件和内容

实验条件：实验所使用的输入图像如图2所示，其中图2(a)是大小为512×512，格式为raw的barbara图像，图2(b)是图2(a)的局部放大图，实验中，各种去噪方法都是使用matlab语言编程实现。

实验内容：在上述实验条件下，分别用KSVD字典学习方法，BM3D方法和本发明方法进行实验。

二.实验结果

(1)使用KSVD字典学习方法对图3所示的含噪图像进行去噪，其结果如图4所示，其中图4(a)是去噪后的图，图4(b)是局部放大图。从图4可以看出，KSVD字典学习方法可以较好的保持同质区域的平滑性，边缘的连续性也相对较好，但有部分纹理细节丢失。

(2)使用BM3D方法对图3所示的含噪图像的进行去噪，其结果如图5所示，其中图5(a)是去噪后的图，图5(b)是局部放大图。从图5可以看出，该滤波方法能够较好保留图像的细节和纹理信息，边缘轮廓明确，但是对某些细小的条纹状纹理的去噪结果失真。

(3)使用本发明对图3所示的含噪图像的进行去噪，其结果如图6所示，其中图6(a)是去噪后的图，图6(b)是局部放大图，从图6可以看出，对比于上述各种去噪方法，本发明的去噪图像清晰，可以在保持同质区域的平滑性的同时保留同质区域的纹理信息，还能很好的保留线目标及一些小的细节信息，边缘的连续性也保持较好。

(4)用峰值信噪比PSNR作为去噪结果的定量评价指标。将KSVD字典学习方法、BM3D方法和本发明方法对大小为512×512的自然图像去噪，其去噪结果的PSNR值对比情况分别如表1，表2和表3所示。其中表1是对barbara图像去噪的PSNR值对比结果；表2是对goldhill图像去噪的PSNR值对比结果；表3是对finger图像去噪的PSNR值对比结果。

表1对barbara图像去噪的PSNR(dB)对比结果

表2对goldhill图像去噪的PSNR(dB)对比结果

表3对finger图像去噪的PSNR(dB)对比结果

从以上表格中可以发现，KSVD字典学习去噪方法和BM3D去噪方法在各种噪声方差上都能保持一定的去噪能力，而本发明相对前面两种方法，其PSNR评价指标和细节显示图在不同噪声方差上都占有绝对的优势。

以上实验结果表明，本发明相对于已有的去噪方法具有更好的性能，平滑同质区域的同时能更好的保持自然图像的边缘和纹理细节。

Claims (3)

1.一种基于字典学习和块匹配的自然图像去噪方法，包括如下步骤：

(1)设含噪图像z(x)＝y(x)+n(x)去噪后的原始图像的估计为：

$\hat{y}\left(x\right)=\arg \underset{D,{\mathrm{\α}}_{\mathrm{ij}},y\left(x\right)}{\min }\mathrm{\λ}{\left|\right|y\left(x\right)-z\left(x\right)\left|\right|}_2^2+\underset{\mathrm{ij}}{\mathrm{\Σ}}\mathrm{\μ}{\left|\right|{\mathrm{\α}}_{\mathrm{ij}}\left|\right|}_0+\underset{\mathrm{ij}}{\mathrm{\Σ}}{\left|\right|{\mathrm{D\α}}_{\mathrm{ij}}-R_{\mathrm{ij}}y\left(x\right)\left|\right|}_2^2+v\left|\right|y\left(x\right)-{\hat{y}}_1\left(x\right){\left|\right|}_2^2$

其中，y(x)为原始图像，n(x)为标准差为σ的高斯白噪声，x∈X，X为图像像素的坐标集合，λ为langrage乘子，且λ＝30/σ，

代表取二范数平方，μ为残差控制因子，α_ij为图像块R_ijy(x)的稀疏表示系数，||||₀代表取零范数，D为过完备稀疏表示字典，R_ij代表从图像中取出

大小的图像块的取块操作，i，j为取出的图像块在原图像中的坐标位置，v为权重参数，它的取值范围是(78/σ)²～(1049/σ)²，

为BM3D方法对含噪图像z(x)的初步去噪结果；

(2)令原始图像等于含噪图像，即y(x)＝z(x)，令过完备稀疏表示字典D为冗余DCT字典，且D∈R^n×K，R^n×K代表n×K大小的实数矩阵，取n＝64，K＝256；

(3)采用KSVD字典学习方法更新过完备稀疏表示字典D和稀疏表示系数α_ij，此步骤迭代J次，当σ≤5时，J＝5，当σ＞5时，J＝10；

(4)用BM3D方法对含噪图像z(x)去噪，得到初步去噪结果

(5)将更新后的过完备稀疏表示字典D、更新后的稀疏表示系数α_ij和初步去噪结果

代入步骤(1)设定的原始图像的估计公式，得到原始图像的估计

$\hat{y}\left(x\right)=[(\mathrm{\λ}+v)I+\underset{\mathrm{ij}}{\mathrm{\Σ}}{R_{\mathrm{ij}}^T{R}_{\mathrm{ij}}]}^{-1}[\mathrm{\λz}\left(x\right)+\underset{\mathrm{ij}}{\mathrm{\Σ}}{R}_{\mathrm{ij}}^T{\mathrm{D\α}}_{\mathrm{ij}}+{v\hat{y}}_1\left(x\right)]$

其中，I为图像大小的单位矩阵，

为R_ij的转置。

2.根据权利要求1所述的基于字典学习和块匹配的自然图像去噪方法，其中步骤(3)所述的采用KSVD字典学习方法更新过完备稀疏表示字典D和稀疏表示系数α_ij，按如下步骤进行：

3a)将原始图像y(x)分解成为

大小的图像块y_ij(x)＝R_ijy(x)，用追踪算法计算图像块y_ij(x)的稀疏表示系数α_ij，即

${\mathrm{\α}}_{\mathrm{ij}}=\arg \underset{{\mathrm{\α}}_{\mathrm{ij}}}{\min }{\left|\right|{\mathrm{\α}}_{\mathrm{ij}}\left|\right|}_{0}s.t{\left|\right|{\mathrm{D\α}}_{\mathrm{ij}}-y_{\mathrm{ij}}\left(x\right)\left|\right|}_2^2\≤{\left(1.15\mathrm{\σ}\right)}^2;$

3b)对于字典D的任一原子d_l，l＝1，2，...，K，找出用到原子d_l的那些图像块的坐标位置集合，即w_l＝{(i，j)|α_ij(l)≠0}，其中式α_ij(l)表示稀疏系数α_ij的第l个元素；

3c)对于

计算残差

其中d_m为字典D的第m个原子，α_ij(m)表示稀疏系数α_ij的第m个元素；

3d)令

组成的集合为E_l，即

3e)对所述的E_l进行奇异值分解，得E_l＝UΔV^T，其中U为左奇异矩阵，Δ为奇异值矩阵，V^T右奇异矩阵的转置，用U中第一列更新原子d_l，用V^T的第一行乘以Δ(1，1)以对α_ij(l)进行更新，其中

3.根据权利要求1所述的基于字典学习和块匹配的自然图像去噪方法，其中步骤(4)所述的用BM3D方法对含噪图z(x)去噪，得到初步去噪结果

按如下步骤进行：

4a)将含噪图像z(x)，分解成8×8大小的图像块，依次取出一块图像块作为当前参考块，设为

令所有图像块为侯选块，设为Z_x；

4b)用块匹配法计算当前参考块

和候选块Z_x的距离

$d(Z_{x_R},Z_x)=\frac{{\left|\right|\mathrm{\γ}(T_{2D}\left(Z_{x_R}\right)-\mathrm{\γ}(T_{2D}\left(Z_x\right)\left|\right|}_2^2}{8^2}$

其中，T_2D表示二维小波线性变换，γ代表硬阈值算子，设硬阈值为0，代表取二范数平方，取出满足

的侯选块作为与参考块相似的块，将这些相似块堆积成三维形式，表示为

其中

表示这些相似块的坐标集合；

4c)在三维域对所述的

作硬阈值联合滤波，得到

的去噪结果

并将去噪结果

拆分二维形式，表示为

x₁∈S₁，其中，T_3D表示三维线性变换，γ′表示硬阈值算子，取硬阈值为2.7σ，

表示逆三维线性变换；

4d)重复步骤4b)-4c)对所有参考块进行处理，得到所有参考块对应相似块的估计

x₁∈S₁，加权平均得到原始图像的基本估计

$\∀x\∈X$

其中，

为位于坐标x₁∈X的图像块的特征函数矩阵，

为

的权值，式中N′为硬阈值去噪后

中非零系数的个数；

4e)将求得的

分解成8×8大小的图像块，依次取出一块图像块作为当前参考块，设为

令所有图像块为侯选块，设为

4f)用块匹配法计算当前参考块

和候选块

的距离

取出满足

的侯选块作为与参考块

相似的块，将这些相似块堆积成三维形式，表示为

其中

表示这些相似块的坐标集合；从含噪图像z(x)中找出位于坐标集合S₂的那些图像块，并堆积成三维形式，表示为

4g)在三维域对所述的进行Wiener联合滤波，得到

的估计

${\hat{Y}}_{S_2}=T_{3D}^{\′-1}\left(W_{S_2}{T}_{3D}^{\′}\left(Z_{S_2}\right)\right)$

其中，T′_3D表示在三维域作2D-DCT线性变换；

代表Wiener滤波，即Wiener经验收缩系数，且

||²为求模的平方，表示逆三维线性变换；将去噪结果

拆分成二维形式，表示为

x₂∈S₂；

4h)重复步骤4f)-4g)对所有参考块进行处理，得到所有参考块对应相似块的估计

X₂∈S₂，加权平均

得含噪图像z(x)的初步去噪结果

$\∀x\∈X$

其中，

为位于坐标x₂∈X的图像块的特征函数矩阵，

为

的权值，表示为

式中

表示对二范数平方取逆。

Images