CN103700072A - 一种基于自适应小波阈值和双边滤波器的图像去噪方法 - Google Patents

Abstract

本发明提出一种基于自适应小波阈值和双边滤波器的图像去噪方法，目的是为了改善小波阈值去噪算法的效果以及更好地保护图像的边缘和细节信息。该算法包括以下步骤：采用离散小波对图像进行分解，获得多个子带及其小波系数；选取随着小波分解尺度和子带的变化而自适应变化的阈值，并采用软阈值函数进行阈值量化处理；进行小波逆变换，获得重构图像；采用双边滤波器对重构图像进行滤波，以便获得清晰的图像。本发明通过利用自适应于小波分解尺度和子带的阈值进行小波阈值去噪，并结合双边滤波器进行滤波，设计出的去噪算法不但能够有效地去除高斯白噪声，而且能很好地保留图像的边缘和细节信息。

Description

一种基于自适应小波阈值和双边滤波器的图像去噪方法

技术领域

本发明涉及图像处理技术领域，特别涉及一种基于自适应小波阈值和双边滤波器的去噪方法。

背景技术

近年来，小波变换由于具有多分辨率分析的特性、在时频域内良好地表征信号的能力以及大小固定形状可变的窗口等特点，被广泛应用于图像去噪领域中。常用的小波域去噪方法是对小波系数进行萎缩处理，即小波阈值去噪方法。其中，经典的阈值是Donoho和Johnstone提出的依赖于噪声标准差和信号长度的统一阈值。

现有的方法存在的不足：一方面，传统的小波阈值去噪方法没有考虑选取的阈值应该随着小波分解尺度和子带的变化而自适应变化这一因素；另一方面，小波阈值去噪方法在去除噪声的同时，可能也会去除高频细节，使图像的边缘变得模糊，而且易于产生显著的伪影，如低频噪声和边缘振荡。

发明内容

本发明的目的旨在解决上述技术缺陷，改善小波阈值去噪算法的效果以及更好地保护图像的边缘和细节信息。

为了达到上述目的，本发明提出一种基于自适应小波阈值和双边滤波器的去噪方法，包括以下几个步骤：

S1：采用离散小波对含有噪声的图像g进行分解，获得多个子带及其对应的小波系数矩阵；

S2：选取每个高频子带对应的随着小波分解尺度和子带变化而自适应变化的阈值，并对每个高频子带的小波系数矩阵进行阈值量化处理，具体如下：

S2.1：根据获得的子带中高频子带的小波系数矩阵，计算所有高频子带的自适应阈值，其中，第k层高频子带的自适应阈值T^k的计算公式如下：

$T^k=\frac{\mathrm{\α}}{k}\mathrm{\σ}\sqrt{2\ln N}$

其中，k为高频子带的小波分解尺度，α为与第k层高频子带有关的可调参数，α的取值范围是(0,1]，ln(·)是以数学常数e为底的对数函数，N是第k层高频子带的尺寸，σ为噪声标准差，且由高频子带HH₁的鲁棒中值来估计；

S2.2：根据获得的子带中高频子带的小波系数矩阵和选取的自适应阈值，采用软阈值函数，对所有高频子带的小波系数矩阵进行阈值量化处理，获得处理后的小波系数矩阵，其中，第k层高频子带在(i,j)处的小波系数通过下面的软阈值函数公式进行阈值量化处理：

$Y_{\mathrm{ij}}^{\&prime;k}=\left\{\begin{array}{cc}\left[\mathrm{sign}\left(Y_{\mathrm{ij}}^k\right)\right]\left(\right|Y_{\mathrm{ij}}^k|-T^k),& \left|Y_{\mathrm{ij}}^k\right|\&GreaterEqual;T^k\\ 0,& \left|Y_{\mathrm{ij}}^k\right|<T^k\end{array}\right.$

其中，

和

分别是上述含噪图像g的第k层高频子带在(i,j)处的小波系数和阈值量化处理后的第k层高频子带在(i,j)处的小波系数，sign(·)是符号函数，T^k为上述选取的第k层高频子带的自适应阈值；

S3：利用上述含噪图像g的低频子带LL_J的小波系数和阈值量化处理后的所有高频子带的小波系数进行小波逆变换，获得重构图像u；

S4：采用双边滤波器对重构图像u进行滤波，获得结果图像v。

步骤S1中所述的多个子带包括一个低频子带LL_J和多个高频子带，其中第k层高频子带包括水平子带HL_k，垂直子带LH_k和对角子带HH_k，k是高频子带的小波分解尺度，k＝1,2,…,J，J为小波分解的最大尺度。

步骤S4中所述的结果图像v中像素点(i,j)的像素值v(i,j)通过重构图像u中像素点(i,j)的邻域像素点(k,l)的像素值u(k,l)的加权组合计算得到，具体公式如下：

$v(i,j)=\frac{{\mathrm{\&Sigma;}}_{(k,l)\&Element;{\mathrm{\&Omega;}}_M(i,j)}w(i,j,k,l)u(k,l)}{{\mathrm{\&Sigma;}}_{(k,l)\&Element;{\mathrm{\&Omega;}}_M(i,j)}w(i,j,k,l)}$

其中，Ω_M(i,j)表示重构图像u中像素点(i,j)处的(2M+1)×(2M+1)的邻域窗，M为邻域半径，像素点(k,l)为邻域窗内的像素，w(i,j,k,l)为权重系数，它取决于空域核函数d(i,j,k,l)和值域核函数r(i,j,k,l)的乘积，其表达式为：

$d(i,j,k,l)=\exp (-\frac{{(i-k)}^2+{(j-l)}^2}{2{\mathrm{\&sigma;}}_d^2})$

$r(i,j,k,l)=\exp (-\frac{{[u(i,j)-u(k,l)]}^2}{2{\mathrm{\&sigma;}}_r^2})$

w(i,j,k,l)＝d(i,j,k,l)*r(i,j,k,l)

其中，exp(·)是以数学常数e为底数的指数函数，参数σ_d和σ_r分别为空域核函数d(i,j,k,l)和值域核函数r(i,j,k,l)的方差。

有益效果

本发明通过利用自适应于小波分解尺度和子带的阈值进行小波阈值去噪，并结合双边滤波器进行滤波，设计出的去噪方法不但能够有效地去除高斯白噪声，而且能很好地保留图像的边缘和细节信息。

附图说明

图1为本发明流程图；

图2为本发明一个实施例的两层小波分解的子带结构示意图，其中，图2(a)为小波分解后的各子带结构，图2(b)为lena图像小波分解后的各子带图像。

具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能解释为对本发明的限制。

如图1所示，本发明的实施例的基于自适应小波阈值和双边滤波器的去噪算法，包括以下几个步骤：

步骤S1，采用离散小波对图像进行分解，获得多个子带及其小波系数。

在本发明的一个实施例中，首先，读入含有噪声的图像信号，即含噪图像g。其中，含噪图像g在像素点(i,j)处的像素值可以表示为g(i,j)＝f(i,j)+n(i,j)，i,j＝1,2,…,m。其中，g(i,j)为含噪图像g在像素点(i,j)处的像素值，f(i,j)为不含噪声的原始图像f在像素点(i,j)处的像素值，n(i,j)为噪声矩阵n在像素点(i,j)处的噪声值，它是独立同分布的高斯白噪声，服从N(0,σ²)分布，m是含噪图像g的宽度。然后，采用离散小波对含噪图像g进行多层小波分解，获得一个低频子带LL_J和多个高频子带及各子带的小波系数。第k层高频子带包括水平子带HL_k，垂直子带LH_k和对角子带HH_k。其中，k是高频子带的小波分解尺度，k＝1,2,…,J，J为小波分解的最大尺度。如图2所示为一个示例的两层小波分解得到的子带结构示意图。其中，图2(a)为小波分解后的各子带结构，图2(b)为lena图像小波分解后的各子带图像。含噪图像g经小波分解后有Y＝X+V，其中，Y、X和V分别为g、f和n的小波系数矩阵。得到的小波系数矩阵Y在(i,j)处的小波系数可以表示为Y_ij。小波阈值去噪算法就是利用一种阈值函数，处理含噪图像的所有高频子带的每一个小波系数Y_ij，来获得阈值量化处理后的小波系数，再经过小波逆变换，得到重构图像的估计量f'＝W^-1X'，其中W^-1为小波逆变换算子，X'为含噪图像g的低频子带LL_J的小波系数与阈值量化处理后的高频子带的小波系数组成的矩阵。

步骤S2，选取随着小波分解尺度和子带的变化而自适应变化的阈值，并采用软阈值函数进行阈值量化处理。

具体地，可以分为以下2个小步骤：

步骤S2.1：选取随着小波分解尺度和子带的变化而自适应变化的阈值。

由信号的奇异性理论可知，利用小波对图像进行多层分解的过程中，噪声和图像信号会表现出不同的性质，即噪声具有负的奇异性，随着分解尺度的增大，噪声的幅度和稠密度会快速减小，噪声的小波系数也会相应地减小，而图像信号则相反，它的小波系数会越来越大。所以在对图像进行小波阈值去噪时，为了最大程度地把每一层的噪声和图像有用信息分离出来，不同分解尺度的阈值的选取应该不同，而且选取的阈值应该随着分解尺度的增大而减小。

此外，对于同一分解尺度，不同方向上的高频子带的小波系数分布也是不同的。在分解后的水平、垂直和对角三个高频子带中，对角子带所含的系数值和取值的数量以及所占的能量百分比相对于其他两个高频子带都要小一些，这说明它所含的信息量更少，噪声更多。因而，不同子带上选取的阈值也应该有所不同，而且对角子带上选取的阈值应该比其他两个高频子带选取的阈值大一些。

但是传统的阈值选取方式并没有随着小波分解尺度和子带的变化而自适应的变化，因此，本发明提出一种随着小波分解尺度和子带变化而变化的自适应阈值。根据上述获得的子带中高频子带的小波系数，计算所有高频子带的自适应阈值。其中，第k层高频子带的自适应阈值T^k可以通过以下公式进行计算：

$T^k=\frac{\mathrm{\&alpha;}}{k}\mathrm{\&sigma;}\sqrt{2\ln N}$

其中，k为高频子带的小波分解尺度，α为与第k层高频子带有关的可调参数，α的取值范围是(0,1]，随着噪声方差的增大，α的取值也会逐渐地增大，且对于同一分解尺度k，对角子带HH_k上选取的参数α的值要比垂直子带LH_k和水平子带HL_k上选取的α值都大一些。ln(·)是以数学常数e为底的对数函数，N是第k层高频子带的尺寸，σ为噪声标准差，且由高频子带HH₁的鲁棒中值来估计：

σ＝median(Y₁)/0.6745

其中，median(·)为中值函数，功能是求取数组的中值，假设高频子带HH₁的小波系数矩阵为

它的行数和列数分别为a和b，对矩阵

中的每个小波系数取绝对值就生成一个新的尺寸为a*b的矩阵

把矩阵

中所有的元素排成一列即得到Y₁。

应理解，上述示例仅为示意性的实施例，并不用于限制本发明，本领域的技术人员还可使用其他方法计算图像的噪声标准差，这些改变和变化均应包含在本发明的保护范围内。

在本发明中，选取的这个自适应阈值是以统一阈值为基础，在其分子上添加与子带有关的可调参数α，在分母上添加小波分解尺度k，使得选取的阈值能够自适应地随着小波分解尺度和子带的变化而变化。该阈值具有以下优点：符合经过小波分解后不同分解尺度和不同子带上的小波系数在对应于噪声和图像信号的比例分布上有所不同的事实，减小了阈值处理后的小波系数和原系数之间的偏差，增强了阈值的实用性，减小了由阈值误断引起的偏差。

步骤S2.2：根据上述获得的子带中高频子带的小波系数和上述选取的自适应阈值，采用软阈值函数，对所有高频子带的小波系数进行阈值量化处理，获得处理后的小波系数。其中，第k层高频子带在(i,j)处的小波系数可以通过下面的软阈值函数公式进行阈值量化处理：

$Y_{\mathrm{ij}}^{\&prime;k}=\left\{\begin{array}{cc}\left[\mathrm{sign}\left(Y_{\mathrm{ij}}^k\right)\right]\left(\right|Y_{\mathrm{ij}}^k|-T^k),& \left|Y_{\mathrm{ij}}^k\right|\&GreaterEqual;T^k\\ 0,& \left|Y_{\mathrm{ij}}^k\right|<T^k\end{array}\right.$

其中，

和

分别是上述含噪图像g的第k层高频子带在(i,j)处的小波系数和阈值量化处理后的第k层高频子带在(i,j)处的小波系数，sign(·)是符号函数，T^k为上述选取的第k层高频子带的自适应阈值。

步骤S3，进行小波逆变换，获得重构图像。

本发明利用上述含噪图像g的低频子带LL_J的小波系数和上述阈值量化处理后的所有高频子带的小波系数进行小波逆变换，获得重构图像u。

步骤S4，采用双边滤波器对重构图像进行滤波，以便获得清晰的图像。

由于小波阈值去噪方法在去除噪声的同时，可能也会去除高频细节，使图像的边缘变得模糊，而且易于产生显著的伪影，如低频噪声和边缘振荡。而双边滤波器则能够通过邻域像素值的非线性加权，来达到去除噪声的同时保留图像边缘的目的。因而可以用它来弥补小波阈值去噪会损失图像边缘的局限。双边滤波器是结合图像的空间邻近度和像素值相似度的一种折衷处理，它包含了空间域的滤波及像素值域的滤波。其中前者考虑的是像素的邻近关系；后者则考虑了图像像素值的变化。作为这两种滤波的非线形组合，双边滤波器的权值不仅与像素的空间距离有关，而且与像素值有关。

在本发明中，利用双边滤波器对上述获得的重构图像u进行滤波，就可以得到本发明的去噪算法的结果图像v。其中，像素点(i,j)处经过双边滤波后的输出像素值v(i,j)能够通过它的邻域像素点(k,l)的像素值u(k,l)的加权组合直接计算出来，它的公式如下：

$v(i,j)=\frac{{\mathrm{\&Sigma;}}_{(k,l)\&Element;{\mathrm{\&Omega;}}_M(i,j)}w(i,j,k,l)u(k,l)}{{\mathrm{\&Sigma;}}_{(k,l)\&Element;{\mathrm{\&Omega;}}_M(i,j)}w(i,j,k,l)}$

其中，Ω_M(i,j)表示像素点(i,j)处的(2M+1)×(2M+1)的邻域窗，M为邻域半径，M越大，平滑效果越强，但M太大则会导致图像模糊。像素点(k,l)为邻域窗内的像素，w(i,j,k,l)为权重系数，它取决于空域核函数d(i,j,k,l)和值域核函数r(i,j,k,l)的乘积，其表达式为：

$d(i,j,k,l)=\exp (-\frac{{(i-k)}^2+{(j-l)}^2}{2{\mathrm{\&sigma;}}_d^2})$

$r(i,j,k,l)=\exp (-\frac{{[u(i,j)-u(k,l)]}^2}{2{\mathrm{\&sigma;}}_r^2})$

w(i,j,k,l)＝d(i,j,k,l)*r(i,j,k,l)

其中，exp(·)是以数学常数e为底数的指数函数，参数σ_d和σ_r分别为空域核函数d(i,j,k,l)和值域核函数r(i,j,k,l)的方差，而且它们分别控制空域和值域中权重因子的衰减程度，若它们的取值太大（趋于无穷），双边滤波器则退化为均值滤波器，取值太小则会削弱平滑效果。

最后，双边滤波器输出的图像v就是本发明的去噪算法的结果图像。

根据本发明实施例的基于自适应小波阈值和双边滤波器的去噪算法，通过利用自适应于小波分解尺度和子带的阈值进行小波阈值去噪，并结合双边滤波器进行滤波，设计出的去噪算法不但能够有效地去除高斯白噪声，而且能很好地保留图像的边缘和细节信息。

尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同限定。

Claims (3)

1.一种基于自适应小波阈值和双边滤波器的图像去噪方法，其特征在于包括以下步骤：

S1：采用离散小波对含有噪声的图像g进行分解，获得多个子带及其对应的小波系数矩阵；

S2：选取每个高频子带对应的随着小波分解尺度和子带变化而自适应变化的阈值，并对每个高频子带的小波系数矩阵进行阈值量化处理，具体如下：

S2.1：根据获得的子带中高频子带的小波系数矩阵，计算所有高频子带的自适应阈值，其中，第k层高频子带的自适应阈值T^k的计算公式如下：

$T^k=\frac{\mathrm{\&alpha;}}{k}\mathrm{\&sigma;}\sqrt{2\ln N}$

其中，k为高频子带的小波分解尺度，α为与第k层高频子带有关的可调参数，α的取值范围是(0,1]，ln(·)是以数学常数e为底的对数函数，N是第k层高频子带的尺寸，σ为噪声标准差，且由高频子带HH₁的鲁棒中值来估计；

S2.2：根据获得的子带中高频子带的小波系数矩阵和选取的自适应阈值，采用软阈值函数，对所有高频子带的小波系数矩阵进行阈值量化处理，获得处理后的小波系数矩阵，其中，第k层高频子带在(i,j)处的小波系数通过下面的软阈值函数公式进行阈值量化处理：

$Y_{\mathrm{ij}}^{\&prime;k}=\left\{\begin{array}{cc}\left[\mathrm{sign}\left(Y_{\mathrm{ij}}^k\right)\right]\left(\right|Y_{\mathrm{ij}}^k|-T^k),& \left|Y_{\mathrm{ij}}^k\right|\&GreaterEqual;T^k\\ 0,& \left|Y_{\mathrm{ij}}^k\right|<T^k\end{array}\right.$

其中，和

分别是上述含噪图像g的第k层高频子带在(i,j)处的小波系数和阈值量化处理后的第k层高频子带在(i,j)处的小波系数，sign(·)是符号函数，T^k为上述选取的第k层高频子带的自适应阈值；

S3：利用上述含噪图像g的低频子带LL_J的小波系数和阈值量化处理后的所有高频子带的小波系数进行小波逆变换，获得重构图像u；

S4：采用双边滤波器对重构图像u进行滤波，获得结果图像v。

2.根据权利要求1所述的一种基于自适应小波阈值和双边滤波器的图像去噪方法，其特征在于，步骤S1中所述的多个子带包括一个低频子带LL_J和多个高频子带，其中第k层高频子带包括水平子带HL_k，垂直子带LH_k和对角子带HH_k，k是高频子带的小波分解尺度，k＝1,2,…,J，J为小波分解的最大尺度。

3.根据权利要求1所述的一种基于自适应小波阈值和双边滤波器的图像去噪方法，其特征在于，步骤S4中所述的结果图像v中像素点(i,j)的像素值v(i,j)通过重构图像u中像素点(i,j)的邻域像素点(k,l)的像素值u(k,l)的加权组合计算得到，具体公式如下：

$v(i,j)=\frac{{\mathrm{\&Sigma;}}_{(k,l)\&Element;{\mathrm{\&Omega;}}_M(i,j)}w(i,j,k,l)u(k,l)}{{\mathrm{\&Sigma;}}_{(k,l)\&Element;{\mathrm{\&Omega;}}_M(i,j)}w(i,j,k,l)}$

其中，Ω_M(i,j)表示重构图像u中像素点(i,j)处的(2M+1)×(2M+1)的邻域窗，M为邻域半径，像素点(k,l)为邻域窗内的像素，w(i,j,k,l)为权重系数，它取决于空域核函数d(i,j,k,l)和值域核函数r(i,j,k,l)的乘积，其表达式为：

$d(i,j,k,l)=\exp (-\frac{{(i-k)}^2+{(j-l)}^2}{2{\mathrm{\&sigma;}}_d^2})$

$r(i,j,k,l)=\exp (-\frac{{[u(i,j)-u(k,l)]}^2}{2{\mathrm{\&sigma;}}_r^2})$

w(i,j,k,l)＝d(i,j,k,l)*r(i,j,k,l)

其中，exp(·)是以数学常数e为底数的指数函数，参数σ_d和σ_r分别为空域核函数d(i,j,k,l)和值域核函数r(i,j,k,l)的方差。

Images