基于非局部小波系数收缩的图像去噪方法

摘要

本发明提出了一种基于非局部小波系数收缩的图像去噪方法，主要解决现有图像去噪方法去噪声时丢失图像细节的问题。其实现步骤为：(1)构建含噪图像的相似组，并对相似组中相似块作二维小波变换,计算相似组小波系数的非局部均值；(2)使用双L1范数模型收缩小波系数，再进行小波反变换，得到相似块估计值，并对其进行整合，得到一次估计图像；(3)对一次估计图像进行残差回补，并执行步骤(1)-(2)，得到基础估计图像；(4)构建基础估计图像的相似组，进而得到含噪图像的相似组；(5)对含噪图像的相似组进行维纳协同滤波，得到去噪图像。本发明能在平滑噪声的同时，更好地保持图像的边缘纹理，可用于对自然图像的去噪处理。

说明书

技术领域

本发明属于图像处理技术领域，具体地说是基于非局部小波系数收缩的图像去噪方 法，可用于对自然图像的去噪处理。

背景技术

图像是人们获取信息的重要来源，但图像在生成和传输过程中常常会受到各种噪声 的干扰，这不仅影响到图像的视觉效果，也阻碍了特征提取、目标识别等后续工作的进 行。因此，图像去噪是图像处理领域至关重要的一部分。

图像去噪的目的就是从含噪图像中恢复出高质量清晰的图像，去噪的同时尽可能地 保持图像的固有特征信息。目前，大量的去噪方法已经被提出，其中正则化方法得到了 广泛的研究，该方法将观测图像和先验知识两者结合于变分公式，强调找到合适的图像 先验模型是至关重要的。传统的正则化方法，比如二次Tikhonov方法、TV方法由于其分 段常数的假设，在去噪的同时会过平滑图像，因此，很快被基于稀疏性的正则化方法所 取代。该方法利用了图像的局部稀疏性，用字典中的若干原子的线性组合来表示图像块， 每个图像块是单独进行稀疏表示的，并没有考虑到与其他图像块之间的相关性。

2005年，Buades等人提出了非局部的图像去噪方法，该方法突破了传统局部滤波的 思想，充分利用图像的非局部结构相似性，获得了显著的去噪效果。随后，出现了很多 基于非局部思想的去噪方法，如BM3D、CSR等方法，BM3D方法将结构相似的二维图像 块聚集在一起形成三维数组，通过对这些三维数组联合滤波，聚合图像块的估计值，进 一步提高了去噪效果。但该方法由于对相似组进行块间一维变换，使得图像的部分细节 被弱化，导致边缘区域变模糊。CSR方法将字典学习和结构聚类结合起来，使图像的稀 疏编码噪声足够小，从而提高去噪效果，但该方法实现起来比较耗时，而且部分边缘的 去噪效果仍不理想。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于非局部小波系数收缩的 图像去噪的实现方法，以实现对自然图像去噪中边缘和平滑区域的兼顾，提高图像去噪 效果。

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

(1)在输入含噪图像中以步长3取参考图像块，根据欧式距离公式计算该参考块与其 邻域内所有图像块的距离d(Z_i,Z_i,j)，选取距离最小的N₂个图像块构成该参考块的相似 组：S_i＝{Z_i,j:min(d(Z_i,Z_i,j)),j＝1...N₂,i表示相似组序号}，其中，Z_i为参考图像块，Z_i,j为Z_i邻域中的图像块，N₂为相似组中相似块个数；

(2)对相似组中的相似块进行二维小波变换,得到相似块的小波系数：

α_i,j＝Τ_2D(Z_i,j),Z_i,j∈S_i

其中，Τ_2D表示二维小波变换，α_i,j为第i个相似组中第j个相似块的小波系数；

(3)根据非局部均值公式，计算每个相似组小波系数的非局部均值μ_i：

${\mu }_i=\underset{j\in S_i}{\Sigma }{\omega }_{i,j}{\alpha }_{i,j},j=1...N_2$

ω_i,j＝exp(-d(Z_i,Z_i,j))/h)/W

其中，ω_i,j为相似块对应的权值，h为12σ_n，σ_n为噪声标准差，W表示归一化操作；

(4)使用双L₁范数模型的收缩函数，计算相似块的小波系数估计值

${\hat{\alpha }}_{i,j}=\left(\begin{array}{cc}{S}_{{\tau }_1,{\tau }_2,{\mu }_i}\left({\alpha }_{i,j}\right)& {\mu }_i\ge 0\\ -S_{{\tau }_1,{\tau }_2,-{\mu }_i}\left({\alpha }_{i,j}\right)& {\mu }_i<0\end{array}\right)$

其中，为收缩函数，其定义为：

$S_{{\tau }_1,{\tau }_2,b}\left(t\right)=\left(\begin{array}{cc}t+{\tau }_1+{\tau }_2& t<-{\tau }_1-{\tau }_2\\ 0& -{\tau }_1-{\tau }_2\le t\le {\tau }_1-{\tau }_2\\ t-{\tau }_1+{\tau }_2& {\tau }_1-{\tau }_2\le t\le {\tau }_1-{\tau }_2+b\\ b& {\tau }_1-{\tau }_2+b\le t\le {\tau }_1+{\tau }_2+b\\ t-{\tau }_1-{\tau }_2& t>{\tau }_1+{\tau }_2+b\end{array}\right)$

其中，t为待收缩系数，b是待收缩系数对应的非局部均值，τ₁和τ₂为两个不同的收 缩阈值，c₁，c₂为两个不同数值的常数，c₁取值为0.1，c₂取值为0.9，σ_n为噪声标准差，σ_i是由α_i估计得到的标准差，α_i＝{α_i,1...α_i,j}，δ_i由α_i-μ_i 估计得到的标准差；

(5)对小波系数估计值进行小波反变换，得到相似块估计值，整合所有估计值， 得到一次估计图像Y：

${\hat{Y}}_{i,j}=T_{2D}^{-1}\left({\hat{\alpha }}_{i,j}\right)$

其中，X为含噪图像，x为X中的像素点，表示二维小波反变换，χ_j是图像块的特征函数，时，χ_j取值为1，否则为0；

(6)对得到的一次估计图像Y进行残差补回，在残差补回图上重复步骤(1)-(5)，得到 基础估计图像Y'；当σ_n＞50时，对基础估计图像Y'再执行一次相同的残差补回操作。

(7)对得到的基础估计图像Y'，以步长N_s选取参考块，根据欧式距离公式构建其相 似组并记录相似组中相似块的坐标信息：

$S_i^1=\{{\hat{Y}}_{i,j}^{\text{'}}:\frac{{\left|\right|{\hat{Y}}_i^{\text{'}}-{\hat{Y}}_{i,j}^{\text{'}}\left|\right|}_2^2}{{\left(N_1^{\text{'}}\right)}^2}<\tau \}$

其中，和分别是基础估计图像中的参考块和候选块，表示的 L₂范数平方，N₁'×N₁'是图像块大小，τ是判断两个图像块是否相似的阈值，σ_n≤40时， τ取值为400，σ_n＞40时，τ取值为3500；

(8)根据相似块坐标信息，从含噪图像中提取相应的图像块，构成相似组对基 础估计图像的相似组和含噪图像的相似组均进行三维变换，得到各自的变换系数， 根据相似组的变换系数计算维纳收缩系数W_i，再根据相似组的变换系数与维纳收缩 系数W_i，得到相似组中相似块估计值：

$W_i=\frac{{\left|T_{3D}\left(S_i^1\right)\right|}^2}{{\left|T_{3D}\left(S_i^1\right)\right|}^2+{{\sigma }_n}^2}$

${\hat{Y}}_{S_i^2}^{\text{'}}=T_{3D}^{-1}\left(W_i{T}_{3D}\left(S_i^2\right)\right)$

其中，Τ_3D表示由二维DCT变换和块间一维小波变换组成的三维变换，表示三 维反变换，表示对相似组的变换系数取绝对值平方，W_i为相似组对应的维 纳收缩系数；

(9)对得到的相似组中相似块估计值加权平均，得到去噪图像Y''：

其中，ω_i是相似组对应的权值，表示W_i的L₂范数平方。

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

1.本发明由于是在小波域中进行的，实现过程简单、速度快。

2.本发明充分利用了图像的非局部相似性，使相似块的小波系数α_i,j接近相似组小 波系数的非局部均值μ_i，能够有效地收缩相似块小波系数，得到准确的相似块估计值。

3.本发明由于使用了双L₁范数模型收缩函数，能够得到准确的基础估计图像，进而 得到准确的含噪图像相似组和维纳收缩系数W_i，从而准确计算出像素点估计值，能够 平滑噪声的同时保持和恢复图像的边缘和纹理细节。

附图说明

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

图2是本发明实验使用的三幅测试图像；

图3是用现有的三种方法和本发明对House图像的去噪结果；

图4是用现有的三种方法和和本发明对Cameraman图像的去噪结果；

图5是用现有的三种方法和本发明对Barbara图像的去噪结果。

具体实施方式

参照附图1，本发明基于非局部小波系数收缩的图像去噪方法，包括如下步骤：

步骤1，对含噪图像X构建其相似组。

1.1)在无噪测试图像中添加标准差为σ_n的噪声，得到含噪图像X：

X=U+σ_n*randn(N)，

其中，U为无噪测试图，N为U中的像素总数，randn()为matlab语言中产生随机 数的函数；

1.2)在含噪图像X中以步长3取参考块Z_i，根据距离公式计算该参考块与其大小为 w×w邻域中候选块的距离d(Z_i,Z_i,j)及相应的权值ω_i,j：

$d(Z_i,Z_{i,j})=\frac{{\left|\right|Z_i-Z_{i,j}\left|\right|}_2^2}{{\left(N_1\right)}^2}$

ω_i,j＝exp(-d(Z_i,Z_i,j))/h)/W，

其中，Z_i为参考块，Z_i,j为Z_i邻域中的图像块，表示Z_i-Z_i,j的L₂范数平 方，N₁×N₁为图像块大小，N₁取值为8，w取值为39，exp表示指数运算，h为12σ_n，σ_n为噪声标准差，W表示归一化操作；

1.3)取距离最小的N₂个图像块，组成参考块Z_i的相似组S_i：

S_i＝{Z_i,j:min(d(Z_i,Z_i,j)),j＝1...N₂,i表示相似组序号}，

其中，min表示取N₂个最小值，N₂为相似组中相似块个数，σ_n≤30时，N₂取值为 16；σ_n＞30时，N₂取值为20。

步骤2，对相似组中相似块进行小波变换，得到相似块的小波系数α_i,j：

α_i,j＝Τ2D(Z_i,j),Z_i,j∈S_i，

其中，Τ_2D表示二维bior1.5小波变换，且为三层分解，α_i,j为第i个相似组中第j个 相似块的小波系数。

步骤3，根据权值ω_i,j和相似块小波系数α_i,j，利用非局部均值公式计算每个相似组 小波系数的非局部均值μ_i：

${\mu }_i=\underset{j\in S_i}{\Sigma }{\omega }_{i,j}{\alpha }_{i,j},j=1...N_2.$

步骤4，使用双L₁范数模型的收缩函数对相似块小波系数α_i,j进行收缩。

4.1)由相似块小波系数α_i,j组成相似组的小波系数集合α_i，α_i＝{α_i,1...α_i,j}，用α_i的 每列减去非局部均值μ_i得到新集合α_i-μ_i，α_i-μ_i＝{(α_i,1-μ_i)...(α_i,j-μ_i)}；

4.2)计算小波系数集合α_i的标准差σ_i和新集合α_i-μ_i的标准差δ_i：

${\sigma }_i=\sqrt{\frac{1}{N_2}\Sigma {\left({\alpha }_i\right)}^2},{\delta }_i=\sqrt{\frac{1}{N_2}\Sigma {({\alpha }_i-{\mu }_i)}^2},$

4.3)根据小波系数集合α_i的标准差σ_i计算第一收缩阈值τ₁，根据新集合α_i-μ_i的标 准差δ_i计算第二收缩阈值τ₂：

${\tau }_1=\frac{c_1\sqrt{2}{{\sigma }_n}^2}{{\sigma }_i},{\tau }_2=\frac{c_2\sqrt{2}{{\sigma }_n}^2}{{\delta }_i},$

其中，σ_n为噪声标准差，c₁，c₂为两个不同数值的常数，c₁取值为0.1，c₂取值为 0.9；

4.4)定义双L₁范数模型的收缩函数为：

$S_{{\tau }_1,{\tau }_2,b}\left(t\right)=\left(\begin{array}{cc}t+{\tau }_1+{\tau }_2& t<-{\tau }_1-{\tau }_2\\ 0& -{\tau }_1-{\tau }_2\le t\le {\tau }_1-{\tau }_2\\ t-{\tau }_1+{\tau }_2& {\tau }_1-{\tau }_2\le t\le {\tau }_1-{\tau }_2+b\\ b& {\tau }_1-{\tau }_2+b\le t\le {\tau }_1+{\tau }_2+b\\ t-{\tau }_1-{\tau }_2& t>{\tau }_1+{\tau }_2+b\end{array}\right),$

其中，t为待收缩系数，b是待收缩系数对应的非局部均值；

4.5)根据两个收缩阈值τ₁、τ₂和相似组小波系数的非局部均值μ_i，利用收缩函数 收缩相似块小波系数α_i,j，得到相似块小波系数估计值

${\hat{\alpha }}_{i,j}=\left(\begin{array}{cc}{S}_{{\tau }_1,{\tau }_2,{\mu }_i}\left({\alpha }_{i,j}\right)& {\mu }_i\ge 0\\ -S_{{\tau }_1,{\tau }_2,-{\mu }_i}\left({\alpha }_{i,j}\right)& {\mu }_i<0\end{array}\right).$

步骤5，根据相似块小波系数估计值得到一次估计图像Y。

5.1)对相似块小波系数估计值进行小波反变换，得到相似块估计值

${\hat{Y}}_{i,j}=T_{2D}^{-1}\left({\hat{\alpha }}_{i,j}\right),$

其中，表示二维小波反变换；

5.2)整合所有相似块估计值得到一次估计图像Y：

其中，X为含噪图像，x为X中的像素点，χ_j是图像块的特征函数，时， χ_j取值为1，否则为0。

步骤6，对一次估计图像Y进行残差补回，得到残差补回图像

$\bar{Y}=Y+\delta (X-Y),$

其中，δ为残差补回参数，取值为0.02。

步骤7，对残差补回图像执行步骤1-5，得到基础估计图像Y'。

7.1)计算残差补回图像的噪声标准差

${\hat{\sigma }}_n=\lambda \sqrt{{{\sigma }_n}^2-{{\delta }_n}^2},$

其中，λ为常数，取值为0.23，为差图的方差，

7.2)根据残差补回图像的噪声标准差对残差补回图像执行步骤1-5，得到基 础估计图像Y'；

7.3)当σ_n＞50时，对基础估计图像Y'再执行一次步骤6、步骤7.1)和步骤7.2)；否 则，直接执行步骤8。

步骤8，对得到的基础估计图像Y'，以步长N_s选取参考块，根据距离公式构建其相 似组

$S_i^1=\{{\hat{Y}}_{i,j}^{\text{'}}:\frac{{\left|\right|{\hat{Y}}_i^{\text{'}}-{\hat{Y}}_{i,j}^{\text{'}}\left|\right|}_2^2}{{\left(N_1^{\text{'}}\right)}^2}<\tau \},$

其中，和分别是基础估计图像中的参考块和候选块，表示的 L₂范数平方，N₁'×N₁'是图像块大小，τ是判断两个图像块是否相似的阈值，σ_n≤40时， N_s取值为3，N₁'取值为8，τ取值为400；σ_n＞40时，N_s取值为6，N₁'取值为11，τ取 值为3500。

步骤9，根据相似组中相似块的位置信息，从含噪图像X中提取相应的图像块， 并用提取出的图像块组成含噪图像X的相似组

步骤10，对含噪图像X的相似组进行维纳协同滤波，得到相似组中相似块估 计值

10.1)对基础估计图像Y'的相似组进行三维变换，根据其变换系数计算维纳收缩 系数W_i：

$W_i=\frac{{\left|T_{3D}\left(S_i^1\right)\right|}^2}{{\left|T_{3D}\left(S_i^1\right)\right|}^2+{{\sigma }_n}^2},$

其中，Τ_3D表示由二维DCT变换和块间一维haar小波变换组成的三维变换，表示对相似组的变换系数取绝对值平方；

10.2)对含噪图像X的相似组进行三维变换，根据其变换系数与维纳收缩系数Wi 计算相似组中相似块的估计值

${\hat{Y}}_{S_i^2}^{\text{'}}=T_{3D}^{-1}\left(W_i{T}_{3D}\left(S_i^2\right)\right),$

其中，W_i为维纳收缩系数，表示三维反变换。

步骤11，对相似组中相似块估计值加权平均，得到去噪图像Y''：

$Y^{\text{'}}=\{{\hat{y}}^{\text{'}}\left(x\right)=\frac{\underset{i\in X}{\Sigma }\underset{j\in S_i^2}{\Sigma }{\omega }_i{\hat{Y}}_{i,j}^{\text{'}}\left(x\right)}{\underset{i\in X}{\Sigma }\underset{j\in S_i^2}{\Sigma }{\omega }_i{\chi }_j\left(x\right)},\forall x\in X\},$

其中，ω_i是相似组对应的权值，表示W_i的L₂范数平方。

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

一.实验条件

实验条件：实验中，使用六幅测试图像进行实验，各种去噪方法都是使用matlab语 言编程实现，且均使用randn('seed',0)添加噪声。

其中三幅测试图像如图2所示，其中，图2(a)是256×256大小的House图像，图2(b) 是对图2(a)图添加了噪声标准差为100噪声的含噪图像，图2(c)是256×256大小的 Cameraman图像，图2(d)是对图2(c)图添加了噪声标准差为75噪声的含噪图像，图2(e) 是512×512大小的Barbara图像，图2(f)是对图2(e)图添加了噪声标准差为50噪声的含 噪图像。

二.实验内容与结果分析：

2.1）在上述实验条件下，分别使用BM3D方法、CSR方法、LSSC方法和本发明进 行实验对图2进行去噪仿真实验。

实验1，使用BM3D方法、CSR方法、LSSC方法和本发明对图2(a)的House图像进 行去噪，结果如图3所示，其中，图3(a)是BM3D方法得到的去噪图像，图3(b)是CSR 方法得到的去噪图像，图3(c)是LSSC方法得到的去噪图像，图3(d)是本发明得到的去噪 图像；

实验2，使用BM3D方法、CSR方法、LSSC方法和本发明对图2(c)的Cameraman 图像进行去噪，结果如图4所示，其中，图4(a)是BM3D方法得到的去噪图像，图4(b) 是CSR方法得到的去噪图像，图4(c)是LSSC方法得到的去噪图像，图4(d)是本发明得 到的去噪图像；

实验3，使用BM3D方法、CSR方法、LSSC方法和本发明对图2(e)的Barbara图像 进行去噪，结果如图5所示，其中，图5(a)是BM3D方法得到的去噪图像，图5(b)是CSR 方法得到的去噪图像，图5(c)是LSSC方法得到的去噪图像，图5(d)是本发明得到的去噪 图像。

从图3、图4和图5可以看出，现有的BM3D方法能够较好的去除噪声，但是会使 图像的边缘模糊且不能保持边缘的连续性，去噪图像中产生了明显的吉布斯效应；相比 BM3D方法，现有CSR方法的去噪效果有所改善，可以较好地恢复边缘纹理区域，对应 的去噪图像中吉布斯效应有所减少，但在去噪的同时仍然会使部分边缘模糊，且均匀区 域不够平滑；现有LSSC方法对图像的平滑区域具有很强的噪声抑制能力，但是会“过 平滑”边缘和纹理区域，甚至无法恢复出部分边缘区域；相比上述现有的三种去噪方法， 本发明在平滑噪声时，能够更好地保留图像细节信息，并引入较少的假信号，尤其在噪 声强度较大时，去噪效果会显著提高。

2.2）计算六幅测试图像去噪后的峰值信噪比PSNR，用PSNR作为去噪的定量评价 指标，其计算方法为：

$\mathrm{PSNR}=101g\left[\frac{u_{\max }^2}{\frac{1}{\left|N\right|}\underset{i\in N}{\Sigma }{[\hat{v}\left(i\right)-u\left(i\right)]}^2}\right]$

其中，u(i)为原始无噪图，为去噪后的结果图，u_max＝max{u(i),i∈N}，N表 示图像大小。

在噪声标准差σ_n为25、35、50、75、100时，对六幅测试图像分别去噪，计算去噪后 的PSNR值，如表1所示。

表1中，每个单元格列出了现有的三种去噪方法和本发明对应的六幅测试图像的 PSNR值，其中每个单元格的左上方为BM3D方法的PSNR值，右上方为CSR方法的PSNR 值，左下方为LSSC方法的PSNR值，右下方为本发明的PSNR值。

表1去噪结果对比

从表1可以看出，本发明相对于现有的三种去噪方法的PSNR评价指标均有所提高。 本发明能够平滑噪声的同时，更好地保持图像的边缘和纹理细节。