CN101957984A - 基于非局部萎缩因子参数估计的图像去噪方法 - Google Patents

Abstract

本发明公开了一种基于非局部萎缩因子参数估计的图像去噪方法，主要克服自然图像原始自适应萎缩去噪中的边缘吉布斯效应明显、视觉质量不理想的问题。其实现过程是：(1)对输入的自然图像进行多尺度变换，得到需处理的子带系数；(2)对子带系数进行初始掩码估计；(3)依次计算子带系数的似然比、先验比和初始萎缩因子；(4)根据子带系数计算非局部滤波子带权重；(5)根据非局部滤波子带权重更新初始萎缩因子；(6)根据非局部权重萎缩因子更新子带系数；(7)用更新的子带系数进行多尺度反变换得到去噪结果。本发明能很好的减弱吉布斯效应，得到较高的PSNR值，可用于对自然图像的去噪处理中。

Description

基于非局部萎缩因子参数估计的图像去噪方法 

技术领域

本发明属于图像处理技术领域，具体地说是一种参数估计方法，可用于对自然图像自适应萎缩去噪处理时更新萎缩因子以提高去噪效果。 

背景技术

数字图像处理由于成像设备及成像条件的限制，使得图像在采集，转换，以及运输过程中不可避免受到噪声的污染。因此图像去噪在图像处理领域占据着举足轻重的地位，成为该领域最基本技术之一。图像处理中许多实际的噪声可以近似的认为是高斯白噪声，因此去除图像中的高斯白噪声成为图像去噪领域中一个重要的方向。 

传统的去噪方法大致可以分为两类，一类是基于空域的方法，一类是基于变换域的方法。 

空域去噪方法中比较经典的方法包括高斯滤波，中值滤波，双边滤波等。它们的共同特点就是利用局部窗口内像素灰度值的连续性来对当前像素进行灰度调整。这些方法的缺点是在去除噪声的同时模糊了图像的细节信息，例如图像的边缘，纹理等。在空域去噪方法中，非局部均值去噪方法以当前像素为中心取大小一定的窗口，在整幅图像内寻找与其具有相似结构的窗口，以窗口之间的相似度为权值对当前像素的灰度值进行调整。这一思想有效的结合了图像系数间的相关性。 

基于变换域的去噪方法比较成熟的是小波域的各种去噪方法，但由于它缺少方向选择性，不适宜表示图像边缘、轮廓等线性奇异性的结构特征，为此，一些新的具有多尺度多方向特性的变换应运而生，如：Brushlet变换、Curvelet变换、Contourlet变换和非下采样Contourlet变换等。 

自适应去噪算法是一种基于变换域的图像去噪方法。这种方法是利用图像在变换域产生的系数进行先验知识——似然比和先验比的计算来更新萎缩系数，从而在变换域对系数进行重新估计，最后进行重构来完成图像去噪。这类自适应去 噪方法能够有效的去除噪声，但它对于萎缩因子的计算方法较为简单，并没有考虑到变换域系数之间的有效相关性，因此去噪结果往往造成平滑效果较差。 

发明内容

本发明的目的在于克服自适应图像去噪的不足，提出了一种基于非局部萎缩因子参数估计的图像去噪方法，以有效地改善平滑力度，提高去噪效果。 

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

1)对输入的含噪图像c进行多尺度变换，将其分解为K＝4层子带，每层子带分解为L＝4个方向，子带系数为： 

，k＝1，…K；l＝1，…L；j＝1，2…512·512，最低频子带系数不作处理； 

2)对子带系数 

进行初始掩码估计，得到掩码估计值 

δ^l为Donoho提出的鲁棒中值阈值， 

为子带系数 的无斑系数，取 

3)计算子带系数 

的似然比 

和方向性先验比 

${\mathrm{\ξ}}_{k,j}^l=\left\{\begin{array}{c}\exp (-\mathrm{\α}),m_{k,j}^l>(1-\mathrm{\δ})T_{k,j}^l\\ \frac{m_{k,j}^l-(1-\mathrm{\α})T_{k,j}^l}{(1+\mathrm{\α})T_{k,j}^l-m_{k,j}^l},(1-\mathrm{\δ})T_{k,j}^l\≤m_{k,j}^l\≤(1+\mathrm{\δ})T_{k,j}^l\\ \exp \left(\mathrm{\α}\right),m_{k,j}^l>(1+\mathrm{\δ})T_{k,j}^l\end{array}\right.,$

其中， 

是 

对应的方向块阈值；δ是方向块控制参数，取值为0.5；α是似然比公式参数，取值为0.5； 是子带系数的模值； 

${{\mathrm{\η}}^l}_{k,j}=\exp \{\mathrm{\γ}\·\underset{i=\left\{\mathrm{1,2,3,4}\right\}}{\max }[\underset{k\∈{\∂}_i}{\mathrm{\Σ}}2{\hat{x}}_{k,j}^l-1\left]\right\},$

其中，γ是似然比控制参数，取值为0.5； 为四个方向各项异性模型； 

4)根据似然比 

和方向性先验比 

，计算子带系数 的初始萎缩因子 

${\mathrm{\ρ}}_{k,j}^l=\frac{{\mathrm{\η}}_{k,j}^l{\mathrm{\ξ}}_{k,j}^l}{1+{\mathrm{\η}}_{k,j}^l{\mathrm{\ξ}}_{k,j}^l};$

5)用非局部方法，在搜索窗Δ内对初始萎缩因子 

进行修正，得到修正后的非局部萎缩因子 

${{\mathrm{\ρ}}_{k,j}^l}^{\′}=\frac{\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{k,j}^l\·\mathrm{\ω}(x_{k,j}^l,x_{k,i}^l)}{\underset{i=1}{\overset{r}{\mathrm{\Σ}}}\mathrm{\ω}(x_{k,j}^l,x_{k,i}^l)},$ r＝1，2...W·W 

其中W是搜索窗Δ的尺度，W＝21； 

是子带系数 与 

的权重，Λ_u，Λ_v分别表示Δ内以子带系数 

为中心的大小为M×M的块，M取7；d(Λ_u，Λ_v)是Λ_u与Λ_v的相似性，通过高斯加权的欧氏距离来衡量，h是平滑参数，h＝0.5σ，其中σ是搜索窗Δ的标准差； 

6)用非局部萎缩因子 

更新子带系数 

，得到新的子带系数： 

k＝1，…K；l＝1，…L；j＝1，2…512·512， 

7)对新的子带系数 

进行多尺度反变换，得到去噪后的图像。 

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

1.本发明由于用非局部萎缩因子参数对初始萎缩因子进行更新，增加了萎缩因子的局部相关性，进而提高了萎缩因子的有效度； 

2.本发明由于用非局部萎缩因子参数对多尺度变换子带系数进行修正，不仅能最大程度的减弱吉布斯现象，而且能够较好抑制噪声，同时保持自然图像的边缘和纹理细节。 

附图说明

图1是本发明的流程图； 

图2是现有非下采样contourlet变换4方向各项异性模型图； 

图3本发明使用的第一幅测试图像； 

图4本发明使用的第二幅测试图像； 

图5是本发明使用的第一幅含噪图像； 

图6是本发明对图5进行去噪的结果图； 

图7是分别用本发明与三种现有方法对图5进行去噪的局部结果图； 

图8是本发明使用的第二幅含噪图像； 

图9是本发明对图8进行去噪的结果图； 

图10是分别用本发明与三种现有方法对图8进行去噪的局部结果图。 

图11是lena图像去噪结果PSNR曲线对比图。 

具体实施方式

参照附图1，本发明基于非局部萎缩因子参数估计的图像去噪方法，包括如下步骤： 

步骤1，对含噪图像进行多尺度变换。 

输入尺度为512×512的含噪自然图像，对该含噪自然图像进行多尺度变换：变换将其分解为K层子带，K＝4；每层子带分解为L个方向，L＝4；分解子带记为： 

，k＝1，…K；l＝1，…L； 

子带系数记为： 

，k＝1，…K；l＝1，…L；j＝1，2…512·512 

具体来说：低频子带系数为 

，次低频子带系数为 

，中频子带系数为 

，高频子带系数为 ，本发明将低频子带 

看做无噪声系数，进行保留，对其余三层子带进行如下步骤2-步骤8的操作。 

步骤2，初始掩码估计。 

对子带系数 

进行掩码估计时，定义该子带系数 

的无噪声系数为 

是比子带系数 

低频一层的相应子带系数； 

初始掩码估计的公式如下： 

${\hat{x}}_{k,j}^l=\left\{\begin{array}{c}0,|c_{k,j}^l\&CenterDot;{\hat{y}}_{k,j}^l|<{\left({\mathrm{\&delta;}}^l\right)}^2\\ 1,|c_{k,j}^l\&CenterDot;{\hat{y}}_{k,j}^l|\&GreaterEqual;{\left({\mathrm{\&delta;}}^l\right)}^2\end{array}\right.,$

其中δ^l为Donoho提出的鲁棒中值阈值： 

Median是指取中值； 即为相应 

层对应的初始掩码估计值。 

步骤3，似然比计算。 

3.1)定义子带系数 的方向性阈值为 

$T_{k,j}^l=3\&CenterDot;\sqrt{\left(\frac{\mathrm{\&Sigma;B}(:)}{7\&CenterDot;7}\right)/0.6745},$

式中，B是指子带系数 

的方向块，如图2四方向各项异性方向模板所示，对子带系数 

计算方向块B，就是取子带系数 为中心的7×7块与子带系数 

所属方向模板的乘积；B(：)是指取方向块B中的所有元素； 

3.2)计算似然比 

${\mathrm{\&xi;}}_{k,j}^l=\left\{\begin{array}{c}\exp (-\mathrm{\&alpha;}),m_{k,j}^l>(1-\mathrm{\&delta;})T_{k,j}^l\\ \frac{m_{k,j}^l-(1-\mathrm{\&alpha;})T_j^l}{(1+\mathrm{\&alpha;})T_{k,j}^l-m_{k,j}^l},(1-\mathrm{\&delta;})T_{k,j}^l\&le;m_{k,j}^l\&le;(1+\mathrm{\&delta;})T_{k,j}^l\\ \exp \left(\mathrm{\&alpha;}\right),m_{k,j}^l>(1+\mathrm{\&delta;})T_{k,j}^l\end{array}\right.,$

δ是方向块控制参数，取值为0.5；α是似然比公式参数，取值为0.5； 

为 

的模值，取 

步骤4，先验比计算。 

根据初始掩码估计值 

，计算先验比 

${{\mathrm{\&eta;}}^l}_{k,j}=\exp \{\mathrm{\&gamma;}\&CenterDot;\underset{i=\left\{\mathrm{1,2,3,4}\right\}}{\max }[\underset{k\&Element;{\&PartialD;}_i}{\mathrm{\&Sigma;}}2{\hat{x}}_{k,j}^l-1\left]\right\},,$

为图2的四个方向各项异性模板，且大小都延拓为9×9，γ是似然比控制参数，取值为0.5。 

步骤5，初始萎缩因子计算。 

用步骤3和步骤4得到的似然比 和先验比 

，计算子带系数 

的初 始萎缩因子 

${\mathrm{\&rho;}}_{k,j}^l=\frac{{\mathrm{\&eta;}}_{k,j}^l{\mathrm{\&xi;}}_{k,j}^l}{1+{\mathrm{\&eta;}}_{k,j}^l{\mathrm{\&xi;}}_{k,j}^l};$

步骤6，非局部子带权重计算。 

6.1)定义Δ是尺度为W×W的搜索窗，W取21； 

6.2)定义子带 

内待估计的子带系数是 

，定义搜索窗Δ中的异于 

的子带系数是 

；Λ_u是以待估计子带系数 

为中心的系数值矩阵，Λ_v是搜索窗Δ中以子带系数 为中心的系数值矩阵，且Λ_u和Λ_v的尺度均为M×M，取M＝7； 

6.3)通过高斯加权的欧氏距离计算Λ_u与Λ_v的相似性： 

$d({\mathrm{\&Lambda;}}_u,{\mathrm{\&Lambda;}}_v)=\underset{r\&Element;\mathrm{\&Lambda;}}{\mathrm{\&Sigma;}}{|{\mathrm{\&Lambda;}}_{u,r}-{\mathrm{\&Lambda;}}_{v,r}|}^2,$

Λ_u，r表示矩阵Λ_u的第r个元素，Λ_v，r表示矩阵Λ_v的第r个元素； 

6.4)计算待估计子带系数 与搜索窗Δ中异于 

的子带系数 

的非局部权重为： 

$\mathrm{\&omega;}(x_{k,j}^l,x_{k,i}^l)=\exp (-\frac{1}{h^2}d({\mathrm{\&Lambda;}}_u,{\mathrm{\&Lambda;}}_v)),$

其中，h是一个平滑参数，h＝0.5σ，σ是搜索窗Δ的标准差。 

步骤7，更新非局部权重萎缩因子。 

根据步骤6，得到以子带系数 

为中心的W×W邻域内共有W²个权重 

，根据权重计算非局部权重萎缩因子 

${{\mathrm{\&rho;}}_{k,j}^l}^{\&prime;}=\frac{\underset{i=1}{\overset{r}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_{k,j}^l\&CenterDot;\mathrm{\&omega;}(x_{k,j}^l,x_{k,i}^l)}{\underset{i=1}{\overset{r}{\mathrm{\&Sigma;}}}\mathrm{\&omega;}(x_{k,j}^l,x_{k,i}^l)},$ r＝1，2...W·W。 

步骤8，子带系数更新。 

用非局部权重萎缩因子对多尺度分解的子带系数 进行更新，得到更新的多尺度变换子带： 

${c_{k,j}^l}^{\&prime;}=c_{k,j}^l\&CenterDot;{{\mathrm{\&rho;}}_{k,j}^l}^{\&prime;},$ k＝2，3，4； 

步骤9，多尺度反变换得到去噪图像。 

用最低频子带 

与更新的多尺度变换子带 进行多尺度反变换，得到去噪后图像。 

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

一.实验条件和内容 

实验条件：实验所使用的输入图像如图3和4所示，其中，图3是测试图像lena，图4是测试图像peppers，图5是图3加入噪声标准差为20的含噪图像，图8是图4加入噪声标准差为20的含噪图像，实验过程中，具体实施方式以现有的NSCT多尺度变换为例； 

实验内容：在上述实验条件下，选用当前去噪领域内几种典型的算法和本发明方法进行实验，它们是：(1)NSCT域自适应萎缩去噪算法NSCT-Shrink；(2)SWT域BayesShrink去噪算法，简称为SWT-BayesShrink；(3)非局部均值滤波原始方法NLM；(4)本发明简称为NLM-Shrink，去噪结果的客观评价指标用峰值信噪比PSNR衡量。 

二.实验结果 

用NSCT-Shrink算法对图5做四层分解的条件下得到的去噪结果，如图7(b)所示。从图7(b)可以看出，lena的头顶部和眼部明显的表现出图像褶皱对去噪结果的影响，说明该方法有一定的噪声抑制能力，但在抑制噪声的同时欠平滑，造成去噪图像有一定褶皱，影响去噪结果。 

用SWT-BayesShrink算法对图5做五层分解的条件下，得到的去噪结果如图7(c)所示。从图7(c)可以看出，，该方法对噪声抑制能力有限，主观效果不好，整幅图像存在小型条状干扰严重。 

用NLM算法对图5进行去噪，其中搜寻窗大小为21×21，相似窗大小为7×7，平滑参数h＝15σ，σ为图像所含噪声标准差，得到的去噪结果如图7(d)所示。从图7(d)可以看出，此方法具备一定噪声抑制能力，但如lena帽子上的纹路明显减少，说明该方法造成图像细节信息丢失的状况严重，不能很好的保持图像的边缘和纹理信息。 

本发明方法对图5的去噪结果如图7(a)所示。从图7(a)可以看出，它的去噪效果要优于上面所提到的所有现有方法，同质区域较平滑，图像的亮度保 持效果较好，而且图像的边缘和细节也得到了很好的保持。 

将本发明和现有三种方法进行去噪的结果用PSNR作为去噪效果的评价指标，做曲线图对比得到图11。从图11可直观的看出本发明得到的PSNR值明显高于现有三种方法得到PSNR值。 

同时，为了更好的说明本发明的有效性，实验中分别用本发明和三种现有方法对图8进行了第二组实验，实验结果如图9和图10所示。 

用NSCT-Shrink算法对图8做四层分解的条件下得到的去噪结果，如图10(b)所示。从图10(b)可以看出，pepper的轮廓出现褶皱，说明该去噪结果产生褶皱，影响了去噪效果。 

用SWT-BayesShrink算法对图8做五层分解的条件下，得到的去噪结果如图10(c)所示。从图10(c)可以看出，该方法对噪声抑制能力有限，整幅图像存在小型条状干扰严重。 

用NLM算法对图8进行去噪，其中搜寻窗大小为21×21，相似窗大小为7×7，平滑参数h＝15σ，σ为图像所含噪声标准差，得到的去噪结果如图10(d)所示。从图10(d)可以看出，此方法具备一定噪声抑制能力，但不能很好的保持图像的边缘和纹理信息。 

本发明方法对图8的去噪结果如图10(a)所示。从图10(a)可以看出，它的去噪效果要优于上面所提到的所有现有方法，图像边缘和细节得到了很好的保持。 

对图3和图4中的测试图像分别加入噪声标准差为10，20，30，40的高斯加性白噪声，用PSNR作为去噪效果的评价指标，将上述三种现有的去噪方和本发明的方法进行比较，各种方法的去噪效果PSNR值列在表1中。 

表1各种去噪算法的PSNR(dB)比较 

表1中的结果均为10次平均后的结果，从表1中可以看出，本发明方法的去噪效果比NSCT-Shrink算法，SWT-BayesShrink算法，NLM算法在PSNR值上都有很大的提高。 

表2是本发明与NSCT域自适应萎缩去噪算法NSCT-Shrink对比： 

表2本发明改进PSNR(dB)提高 

从表2看出，本发明与NSCT-Shrink相比，当噪声标准差σ_n越大，本发明的PSNR值提高越多。说明本发明适用在噪声含量较高的图像去噪应用中。 

综合以上，本发明无论在客观指标还是主观效果上，都表现出了较好的性能，更好地平滑噪声的同时保持自然图像的边缘和纹理等细节，同时消除了吉布斯现象。 

Claims (2)

1.一种基于非局部萎缩因子参数估计的图像去噪方法，包括如下步骤：

1)对输入的含噪图像c进行多尺度变换，将其分解为K＝4层子带，每层子带分解为L＝4个方向，子带系数为：

，k＝1，…K；l＝1，…L；j＝1，2…512·512，最低频子带系数不作处理；

2)对子带系数

进行初始掩码估计，得到掩码估计值

δ^l为Donoho提出的鲁棒中值阈值，为子带系数

的无斑系数，取

3)计算子带系数

的似然比

和方向性先验比

${\mathrm{\&xi;}}_{k,j}^l=\left\{\begin{array}{c}\exp (-\mathrm{\&alpha;}),m_{k,j}^l>(1-\mathrm{\&delta;})T_{k,j}^l\\ \frac{m_{k,j}^l-(1-\mathrm{\&alpha;})T_{k,j}^l}{(1+\mathrm{\&alpha;})T_{k,j}^l-m_{k,j}^l},(1-\mathrm{\&delta;})T_{k,j}^l\&le;m_{k,j}^l\&le;(1+\mathrm{\&delta;})T_{k,j}^l\\ \exp \left(\mathrm{\&alpha;}\right),m_{k,j}^l>(1+\mathrm{\&delta;})T_{k,j}^l\end{array}\right.,$

其中，

是

对应的方向块阈值；δ是方向块控制参数，取值为0.5；α是似然比公式参数，取值为0.5；

是子带系数的模值；

${{\mathrm{\&eta;}}^l}_{k,j}=\exp \{\mathrm{\&gamma;}\&CenterDot;\underset{i=\left\{\mathrm{1,2,3,4}\right\}}{\max }[\underset{k\&Element;{\&PartialD;}_i}{\mathrm{\&Sigma;}}2{\hat{x}}_{k,j}^l-1\left]\right\},$

其中，γ是似然比控制参数，取值为0.5；为四个方向各项异性模型；

4)根据似然比

和方向性先验比

，计算子带系数

的初始萎缩因子

${\mathrm{\&rho;}}_{k,j}^l=\frac{{\mathrm{\&eta;}}_{k,j}^l{\mathrm{\&xi;}}_{k,j}^l}{1+{\mathrm{\&eta;}}_{k,j}^l{\mathrm{\&xi;}}_{k,j}^l};$

5)用非局部方法，在搜索窗Δ内对初始萎缩因子

进行修正，得到修正后的非局部萎缩因子

${{\mathrm{\&rho;}}_{k,j}^l}^{\&prime;}=\frac{\underset{i=1}{\overset{r}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_{k,j}^l\&CenterDot;\mathrm{\&omega;}(x_{k,j}^l,x_{k,i}^l)}{\underset{i=1}{\overset{r}{\mathrm{\&Sigma;}}}\mathrm{\&omega;}(x_{k,j}^l,x_{k,i}^l)},$ r＝1，2...W·W，

其中W是搜索窗Δ的尺度，W＝21；

是子带系数

与

的权重，Λ_u，Λ_v分别表示Δ内以子带系数为中心的大小为M×M的块，M取7；d(Λ_u，Λ_v)是Λ_u与Λ_v的相似性，通过高斯加权的欧氏距离来衡量，h是平滑参数，h＝0.5σ，其中σ是搜索窗Δ的标准差；

6)用非局部萎缩因子

更新子带系数，得到新的子带系数：

k＝1，…K；l＝1，…L；j＝1，2…512·512，

7)对新的子带系数

进行多尺度反变换，得到去噪后的图像。

2.根据权利要求1所述的基于非局部萎缩因子参数估计的图像去噪方法，其特征在于对步骤(5)得到的萎缩系数加入非局部权重进行修正，按如下步骤进行：

2a.对多尺度变换得到的分解子带系数

进行非局部权重计算：

设

是预进行非局部权重修正的子带系数，

是搜寻窗Δ内，预与

进行相似性比较的子带系数，子带系数

与

的权重为：

$\mathrm{\&omega;}(x_{k,j}^l,x_{k,i}^l)=\exp (-\frac{1}{h^2}d({\mathrm{\&Lambda;}}_u,{\mathrm{\&Lambda;}}_v)),$

其中，h为平滑控制参数，h＝0.5σ，σ是搜索窗Δ的标准差；

Λ_u，Λ_v分别表示Δ中以

为中心的大小为M×M的块，M取7；Λ_u，r，Λ_v，r分别是Λ_u和Λ_v的第r个系数值；Λ_u与Λ_v的相似性通过高斯加权的欧氏距离来衡量，即：

2b.计算非局部权重萎缩因子参数

以子带系数

为中心的W×W邻域内共有W²个权重

，根据权重计算非局部权重萎缩因子

${{\mathrm{\&rho;}}_{k,j}^l}^{\&prime;}=\frac{\underset{i=1}{\overset{r}{\mathrm{\&Sigma;}}}{\mathrm{\&rho;}}_{k,j}^l\&CenterDot;\mathrm{\&omega;}(x_{k,j}^l,x_{k,i}^l)}{\underset{i=1}{\overset{r}{\mathrm{\&Sigma;}}}\mathrm{\&omega;}(x_{k,j}^l,x_{k,i}^l)},$ r＝1，2...W·W。

2c.用非局部权重萎缩因子对多尺度分解的子带系数

进行更新，得到更新的多尺度变换子带：

k＝2，3，4，l＝1，…L；j＝1，2…512·512；

Images