CN101430788A - 基于第二代Curvelet变换的弧形窗局部Wiener滤波方法 - Google Patents

Abstract

本发明公开一种基于第二代Curvelet变换的弧形窗局部Wiener滤波方法，它涉及数字图像处理领域，主要克服现有基于第二代Curvelet变换的阈值法对系数的“过扼杀”及方形窗口方向选择性差的不足。其实现步骤为：(1)选取测试对象，加入高斯噪声，得到噪声图像；(2)对噪声图像进行第二代Curvelet变换；(3)在变换域进行各子带噪声方差的估计；(4)在变换后的各子带中，选取一个以当前要处理系数为中心的弧形区域；(5)估计当前系数的信号方差；(6)对每一个系数进行局部Wiener滤波；(7)对处理过的系数进行第二代Curvelet反变换，得到原始图像的估计。本发明具有系数萎缩灵活、方向选择性好、边缘及细节信息清晰和峰值信噪比高的优点，可用于滤除自然图像中的高斯噪声。

Description

基于第二代Curvelet变换的弧形窗局部Wiener滤波方法

技术领域

本发明属于数字图像处理领域，具体地说是一种基于第二代Curvelet变换的弧形窗Wiener滤波方法，该方法可用于滤除被污染的自然图像中的噪声。

背景技术

图像在生成和传输过程中常常因受到各种噪声的干扰和影响而降质，这对后续的图像处理，如分割、压缩等，将产生不利影响。所以在对图像进一步分析之前，要滤除噪声。图像去噪的主要目标是降低噪声水平的同时又保护图像的特征，而人们也根据实际图像的特点、噪声的统计特征和频谱分布的规律，发展了各式各样的去噪方法，其中最为直观的方法是根据噪声能量一般集中于高频，而图像频谱则分布于一个有限区间的这一特点，采用低通滤波方式来进行去噪的方法，例如滑动平均窗滤波器，还有Wiener线性滤噪器等。其他的去噪方法还有基于秩-阶滤波的方法、基于马尔可夫场模型和基于偏微分方程的方法和LP正则化方法等。

1992年，Donoho和Johnstone提出了小波阈值萎缩方法；与此同时，Krim等人运用Rissanen的MDL准则，也得到了相同的阈值公式，此后小波阈值萎缩方法被用到各种去噪应用中，并取得了很大的成功，对高斯噪声尤其如此。但是Donoho和Johnstone给出的通用阈值，由于有很严重的“过扼杀”系数的倾向，因此人们纷纷对阈值的选择进行了研究，并提出了多种不同的阈值确定方法。相对阈值萎缩方法来说，比例萎缩有更大的灵活性，从某种意义上说，可以认为阈值萎缩是比例萎缩的一种特例，比例萎缩的特点主要在于它具有对信号的某一局部的适应能力，2006年，Pizurica.A等人提出了UWT域空间自适应Bayesian缩减方法，该方法采用方形窗口来估计系数信号方差，考虑了一定的邻域信息，应用于图像去噪，得到了很好的效果，但是这类方形窗口的不足在于它是各向同性的，没有考虑方向信息，所以有其自身的局限性。

小波作为一种逼近的工具，在处理一维和二维的具有点状奇异性的目标中体现出了良好的性能。然而对图像而言，边缘的不连续性是按空间分布的，这种奇异性影响了小波展开级数中的许多项，所以，小波变换在处理图像边缘时的效果不是很好。Candes和Donoho发展了Curvelet变换，这种变换沿着曲线描绘边缘和其他异常点，比其他传统变换更有效。它可以用很少的系数便能实现相当精确的重构。Curvelet变换和小波变换一样，也是一种多尺度变换，它的结构元素也包括尺度和位置参数，但和小波变换不一样的是，它还包括方向参数，使得Curvelet变换具有良好的方向特性。由于Curvelet变换具有多尺度特性以及良好的方向特性，噪声信息和边缘信息能够很好的分开，在保持边缘的同时，使噪声抑制达到了一个很好的效果。Candes和Donoho等人提出的基于Curvelet变换的硬阈值去噪方法，就证明了这一点，但这种方法的最大不足在于对系数的“过扼杀”，即只要比阈值小的系数，一律被作为噪声处理，设置为零，这样就导致图像边缘出现模糊，影响图像质量。

发明内容

本发明的目的在于克服上述已有技术的不足，提供一种基于第二代Curvelet变换的弧形窗局部Wiener滤波方法，该方法以第二代Curvelet为变换工具，采用更能与子带系数的能量分布趋势相一致的弧形窗口来估计系数的信号方差，以达到更好的方向选择性，保持图像边缘清晰，提高图像质量。

实现本发明技术目的技术方案，包括如下步骤：

(1)选取测试对象，加入高斯噪声，得到噪声图像；

(2)对所要处理的噪声图像进行第二代Curvelet变换；

(3)对变换后的噪声图像进行各子带噪声方差

的估计；

(4)根据Curvelet的曲线特征，在变换后的各子带中，选取一个以当前要处理系数为中心的弧形区域；

(5)利用选取的弧形区域中的所有系数估计当前系数的信号方差

(6)把当前系数的信号方差

和所在子带的噪声方差

代入到局部Wiener滤波公式，对该当前系数进行局部Wiener滤波；

(7)对滤波后的系数进行第二代Curvelet反变换，得到原始图像的估计。

上述的滤波方法步骤(4)所说的在变换后的各子带中，选取一个以当前要处理系数为中心的弧形区域，按如下过程进行：

1)设定弧形区域的公式：

W＝{(m，n)||m|≤rsinθ，|n|≤r(1-cosθ)}

式中，r是半径，r越大，窗内包含的系数个数越多；θ是弧形窗口对应中心角的一半，通过调节θ的值能够使弧形窗有更好的方向选择性；

2)选取合适的r值：根据多次试验结果证明，r取2时效果最好；

3)选取合适的θ值：

${\mathrm{\θ}}_0=(\frac{3}{4}-\frac{1}{2^l})\mathrm{\π},$ $\mathrm{\θ}={\mathrm{\θ}}_0-\frac{(k-1)\mathrm{\π}}{{2M}_l}$

式中，l是当前子带的分解尺度，l＝1，...L，L为分解的最大尺度数，M_l为每一个分解尺度l上的子带总数，k表示当前子带。

上述的滤波方法步骤(5)所说的系数信号方差

的估计过程如下：

1)对于当前待处理系数C(l，k，i，j)，确定所选弧形区域中的所有系数，C(l，k，i+p，j+q)，(p，q)∈W，W表示弧形区域，l是当前子带的分解尺度，k表示当前子带，(i，j)是当前系数在子带中的坐标；

2)利用下式估计当前系数的信号方差

$\widehat{{\mathrm{\σ}}_s^2}(l,k,i,j)={(\frac{1}{\#W}\underset{(p,q)\∈W}{\mathrm{\Σ}}{C}^2(l,k,i+p,j+q)-\widehat{{\mathrm{\σ}}^2}(l,k))}_{+}$

式中，W表示弧形区域，#W表示弧形区域内系数的个数，

${(\frac{1}{\#W}\underset{(p,q)\∈W}{\mathrm{\Σ}}{C}^2(l,k,i+p,j+q)-\widehat{{\mathrm{\σ}}^2}(l,k))}_{+}=\max \{0,\frac{1}{\#W}\underset{(p,q)\∈W}{\mathrm{\Σ}}{C}^2(l,k,i+p,j+q)-\widehat{{\mathrm{\σ}}^2}(l,k)\}.$

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

1、本发明的方差估计具有空间局部适应性

由于信号方差

的值不仅随尺度、方向的变化而变化，而且在同一子带内，不同的空间位置的相差也很大，因此，在同一子带内采用相同的

估计是不合理的，保证方差估计具有空间局部自适应特点的方法是Mihcak提出的基于局部自适应窗口的极大似然去噪方法和局部自适应窗口的最大后验概率估计法，当前系数的

由该系数邻域的系数来估计，选择以该系数为中心的正方形窗口。

由于邻域大小的选择在

的估计和整个去噪效果上的影响较大，窗口增大时，将导致

估计偏小，适应低信噪比的系数；而窗口减小时可能导致估计的准确性降低，适应高信噪比的系数，因此，在不同子带内选用不同的邻域大小，原则是细尺度用大窗口，粗尺度用小窗口，这样虽然可以减少去噪误差，但计算复杂性也增加了，而且，方形窗本身就缺乏方向选择性，不能更好的使

的估计具有局部空间适应性，本发明提出的弧形窗克服了这一不足，把这个弧形窗口用于第二代Curvelet变换域中，使窗口的方向与子带系数的能量分布趋势相一致，所估计的方差具有更好的空间局部适应性。

2、本发明的系数萎缩方法具有较大的灵活性

Candes和Donoho等人提出的基于Curvelet变换的硬阈值萎缩去噪方法，对每一个子带中所有的系数采用同一个阈值进行萎缩处理，只要比阈值小的系数，一律被作为噪声处理，设置为零，很明显，该方法灵活性很小，相对阈值萎缩方法来说，比例萎缩有更大的灵活性，从某种意义上说，可以认为阈值萎缩是比例萎缩的一种特例，比例萎缩的特点主要在于它具有对信号的某一局部的适应能力，本发明采用的局部Wiener滤波就是一种比例萎缩方法，所以该方法具有更大的灵活性。

仿真实验结果表明。本发明在对加有高斯噪声的自然图像去噪中，得到了更高的峰值信噪比。

附图说明

图1是本发明的主要操作过程示意图；

图2是用本发明与已有方法对测试图像Lena512×512的去噪效果对比图；

图3是用本发明与已有方法对测试图像Barbara512×512的去噪效果对比图；

图4是用本发明与已有方法对测试图像Peppers256×256的去噪效果对比图。

具体实施方式

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

步骤1：选取测试对象，加入高斯噪声，得到噪声图像，该噪声图像的数学表达式为：

f＝g+n

其中g＝{g(i，j)|i，j＝1，2，...N}表示原始图像，n＝{n(i，j)|i，j＝1，2，...N}表示零均值和方差等于σ²的高斯噪声，含噪图像记为f＝{f(i，j)|i，j＝1，2，...N}，N表示图像大小。

步骤2：对所要处理的噪声图像进行第二代Curvelet变换，得到含噪图像在第l层的第k个子带位置(i，j)处的Curvelet系数：

C(l，k，i，j)＝s(l，k，i，j)+ε(l，k，i，j)

其中s(l，k，i，j)和ε(l，k，i，j)分别是原始图像和噪声在第l层的第k个子带位置(i，j)处的Curvelet系数。

步骤3：对变换后的噪声图像进行各子带噪声方差

的估计，即采用J-L.Starck等人给出的Monte-Carlo方法对该噪声方差进行求解，具体过程通过程序仿真实现。

步骤4：根据Curvelet的曲线特征，在变换后的各子带中，选取一个以当前要处理系数为中心的弧形区域。在Curvelet变换域中，曲线特性更为显著，据此，以当前待处理系数为中心，选取弧形区域，该弧形所对应的中心角由当前子带的尺度信息和方向信息共同决定，从而使弧形区域有更好的方向选择性，对于弧形窗所对应圆的半径，选取多次实验后效果最好的那个值。定义弧形窗表达式为：

W＝{(m，n)||m|≤rsinθ，|n|≤r(1-cosθ)}

式中，r是半径，r越大，窗内包含的系数个数越多；θ是弧形窗口对应中心角的一半，通过调节θ的值能够使弧形窗有更好的方向选择性，可以看出θ＝90°时，该弧形窗就变为方形窗口，所以弧形窗比方窗更具有通用意义，这里，我们选取：

${\mathrm{\θ}}_0=(\frac{3}{4}-\frac{1}{2^l})\mathrm{\π},$ $\mathrm{\θ}={\mathrm{\θ}}_0-\frac{(k-1)\mathrm{\π}}{{2M}_l}$

式中，l是当前子带的分解尺度，l＝1，...L，L为分解的最大尺度数，M_l为每一个分解尺度l上的子带总数，k表示当前子带。

对于r的选取，根据多次实验，通过比较发现，r取2时效果最好。

步骤5：利用选取的弧形区域中的所有系数估计当前系数的信号方差

当前系数的信号方差由当前系数C(l，k，i，j)的弧形邻域中的系数来估计，具体求解公式如下：

$\widehat{{\mathrm{\σ}}_s^2}(l,k,i,j)={(\frac{1}{\#W}\underset{(p,q)\∈W}{\mathrm{\Σ}}{C}^2(l,k,i+p,j+q)-\widehat{{\mathrm{\σ}}^2}(l,k))}_{+}$

式中，W表示弧形区域，#W表示弧形区域内系数的个数，

${(\frac{1}{\#W}\underset{(p,q)\∈W}{\mathrm{\Σ}}{C}^2(l,k,i+p,j+q)-\widehat{{\mathrm{\σ}}^2}(l,k))}_{+}=\max \{0,\frac{1}{\#W}\underset{(p,q)\∈W}{\mathrm{\Σ}}{C}^2(l,k,i+p,j+q)-\widehat{{\mathrm{\σ}}^2}(l,k)\}.$

步骤6：把当前系数的信号方差

和所在子带的噪声方差

代入到局部Wiener滤波公式，得到如下滤波公式为：

$\hat{s}(l,k,i,j)=\frac{\widehat{{\mathrm{\σ}}_s^2}(l,k,i,j)}{\widehat{{\mathrm{\σ}}_s^2}(l,k,i,j)+\widehat{{\mathrm{\σ}}^2}(l,k)}C(l,k,i,j)$

利用该滤波公式对当前系数进行局部Wiener滤波。

步骤7：对滤波后各子带的Curvelet系数采用EJ.Candes等人给出的Wrap算法进行第二代Curvelet反变换，得到原始图像的估计，从而实现滤除噪声的目的。

本发明的效果通过以下仿真进一步说明。

1、仿真条件

采用图像去噪中常用的标准Lena 512×512、Barbara 512×512和peppers256×256图像，并对这三幅图像分别加入标准差σ为10，15，20和25的零均值加性高斯白噪声，共产生12幅含噪图像作为测试数据。分别用基于第二代Curvelet变换弧形窗局部Wiener滤波方法、基于第二代Curvelet变换的硬阈值去噪方法、小波BayesShrink阈值方法和小波BivaShrink阈值方法进行去噪处理。

2、仿真结果分析

仿真的结果如图2、图3和图4。

从图2中可以看出，本发明去噪后不仅有较高的峰值信噪比，而且在边缘和细节信息的保持方面也优于其它三种去噪方法，比如经本发明去噪后Lena的帽沿及帽子上的纹理比较清晰，而用基于第二代Curvelet变换的硬阈值方法去噪后，帽沿比较模糊，帽子上的纹理平滑得也很严重，用基于小波的两种方法去噪后虽然帽沿也比较清晰，但帽子上的纹理有些模糊。

从图3中可以看出，本发明在滤除噪声的同时更多的保持了Barbara裤子上的边缘和纹理细节；可见，本发明在细节保持方面有相当大的优势。

从图4中可以看出，本发明去噪后峰值信噪比提高相对不是很大，但各种蔬菜的边缘很清晰，这一点明显比其他三种方法去噪后的效果要好。

总之，本发明不仅视觉效果好而且峰值信噪比也比较高，尤其在边缘保持方面优势更大，这一点在图像数据信息比较多，即图像比较大时更明显一些。

Claims (3)

1、一种基于第二代Curvelet变换的弧形窗局部Wiener滤波方法，包括如下步骤：

(1)选取测试对象，加入高斯噪声，得到噪声图像；

(2)对所要处理的噪声图像进行第二代Curvelet变换；

(3)对变换后的噪声图像进行各子带噪声方差

(l，k)的估计；

(4)根据Curvelet的曲线特征，在变换后的各子带中，选取一个以当前要处理系数为中心的弧形区域；

(5)利用选取的弧形区域中的所有系数估计当前系数的信号方差(l，k，i，j)；

(6)将当前系数的信号方差

(l，k，i，j)和所在子带的噪声方差

(l，k)代入到局部Wiener滤波公式，对该当前系数进行局部Wiener滤波；

(7)对滤波后的系数进行第二代Curvelet反变换，得到原始图像的估计。

2、根据权利要求1所述的滤波方法，其中步骤(4)所说的在变换后的各子带中，选取一个以当前要处理系数为中心的弧形区域，按如下步骤进行：

(2a)设定弧形区域的公式：

W＝{(m，n)‖m|≤rsinθ，|n|≤r(1-cosθ)}

式中，r是半径，θ是弧形窗口对应中心角的一半；

(2b)选取合适的r值：根据多次试验结果证明，r取2时效果最好；

(2c)选取合适的θ值：

${\mathrm{\θ}}_0=(\frac{3}{4}-\frac{1}{2^l})\mathrm{\π},$   $\mathrm{\θ}={\mathrm{\θ}}_0-\frac{(k-1)\mathrm{\π}}{2M_l}$

式中，l是当前子带的分解尺度，l＝1，...L，L为分解的最大尺度数，M_l为每一个分解尺度l上的子带总数，k表示当前子带。

3、根据权利要求1所述的滤波方法，其中步骤(5)所说的系数信号方差

的估计，按如下步骤进行：

(3a)对于当前待处理系数C(l，k，i，j)，确定所选弧形区域中的所有系数为：C(l，k，i+p，j+q)，(p，q)∈W，W表示弧形区域，l是当前子带的分解尺度，k表示当前子带，(i，j)是当前系数在子带中的坐标；

(3b)利用下式估计当前系数的信号方差

$\widehat{{\mathrm{\σ}}_s^2}(l,k,i,j)={(\frac{1}{\#W}\underset{(p,q)\∈W}{\mathrm{\Σ}}{C}^2(l,k,i+p,j+q)-\widehat{{\mathrm{\σ}}^2}(l,k))}_{+}$

式中，W表示弧形区域，#W表示弧形区域内系数的个数，

${(\frac{1}{\#W}\underset{(p,q)\∈W}{\mathrm{\Σ}}{C}^2(l,k,i+p,j+q)-\widehat{{\mathrm{\σ}}^2}(l,k))}_{+}=\max \{0,\frac{1}{\#W}\underset{(p,q)\∈W}{\mathrm{\Σ}}{C}^2(l,k,i+p,j+q)-\widehat{{\mathrm{\σ}}^2}(l,k)\}.$

Images