CN102622731B - 基于二维Otsu的轮廓波域维纳滤波图像去噪方法 - Google Patents

Abstract

本发明公开了一种基于二维Otsu的轮廓波域维纳滤波图像去噪方法，主要解决现有去噪方法去噪效果不佳的问题。其实现步骤是：(1)对含噪图像TI进行轮廓波分解；(2)分别对分解出的各个高频子带进行二维Otsu分割，得到重要系数和非重要系数；(3)分别计算高频子带的椭圆窗口，根据椭圆窗口估计高频子带的信号方差，对重要系数和非重要系数分别进行维纳滤波；(4)对去噪后的高频子带进行轮廓波逆变换，得到去噪图像FI；(5)对FI进行非局部均值滤波，得到去噪输出。本发明能有效去除含有高斯白噪声的自然图像中的噪声，可用于变化检测，目标识别时对图像的预处理。

Description

基于二维Otsu的轮廓波域维纳滤波图像去噪方法

技术领域

本发明属于图像处理技术领域，涉及含噪情况下自然图像的去噪，可用于森林资源调查、环境变化评估、医学影像等领域的数字图像预处理。

背景技术

由于受各种条件的限制，图像在获取、编码、传输过程中会受到各种各样的噪声影响，图像去噪解决了由噪声干扰导致的图像质量下降问题，抑制了噪声影响，提高了图像质量，因此图像去噪技术在很多领域中都占据着重要位置。

图像去噪方法中应用广泛的是基于小波变换和轮廓波变换的滤波方法，轮廓波变换采用可调的长方形基函数，能准确地将图像中的边缘捕获到不同尺度、不同方向的子带中，实现了对图像信号的稀疏分离，因此十分适合去噪。

X Li等在文献“Image Denoising via Doubly Wiener Filtering with AdaptiveDirectional Windows and Mean Shift Algorithm in Wavelet Domain.ICMA，2010：114-118.”中提出了一种基于自适应窗口的小波域双维纳滤波去噪方法。该方法首先对含噪图像进行一次维纳滤波得到参考图像，再利用均值漂移对参考图像进行分割，然后变换到小波域估计信号方差，对含噪图进行第二次维纳滤波。该方法的缺点是第一次维纳滤波时图像信号损失过多，导致估计的信号方差误差过大，去噪图像噪声残留过多。

Z.F.Zhou等在文献“Contourlet-based Image Denoising Algorithm using AdaptiveWindows.ICIEA，2009：3654-3657.”中提出了一种基于自适应窗口的轮廓波域维纳滤波去噪方法。该方法首先在高频子带内选择各向异性窗口，然后对高频系数进行该窗口的维纳滤波。该方法的缺点是轮廓波不具有平移不变性且维纳滤波的萎缩因子选取不够准确，导致去噪图像含有伪Gibbs现象且噪声有残留。

Q Zhao等在文献“Q Zhao，B Ye，J.L.Cao.Image Denoising Based on ImprovedNon-local Means and Nonsubsampled Contourlet Transform Wiener Filtering.Journal ofComputational Information Systems，2010，6(2)：601-610.”中提出一种基于非局部均值的NSCT域维纳滤波去噪方法去噪方法，该方法首先对图像进行非局部均值预滤波去除大噪声，然后将滤波后图像变换到非下采样轮廓波(Non-Subsampled Contourlet，NSCT)域，对NSCT高频子带进行方形窗口的维纳滤波。该方法的缺点是第一步非局部均值预滤波时图像信号损失过多，导致估计的维纳萎缩因子过小，去噪图像仍含有很多噪声。

发明内容

本发明的目的在于针对上述去噪方法的不足，提出一种基于二维Otsu的轮廓波域维纳滤波图像去噪方法，以保持图像细节，减少噪声残留，提高去噪图像的清晰度。

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

(1)对大小为A×B的含噪图像TI先进行轮廓波分解，分解层数为3，第一层分解得到16个高频子带，第二层分解得到8个高频子带，第三层分解得到8个高频子带，再采用绝对值中值法估计含噪图像TI的噪声noise_TI；

(2)采用二维Otsu分别对以上分解出的各个高频子带进行分割，将高频子带中的系数分为重要系数和非重要系数；

(3)用含噪图像TI的噪声noise_TI分别估计各个高频子带的噪声方差并计算各个高频子带的椭圆窗口P；

(4)对各个高频子带进行维纳滤波：

(4a)根据椭圆窗口P计算高频子带中点(i，j)处的信号方差：

${\mathrm{\σ}}_k^2(i,j)=\max (0,\frac{1}{M_s}\underset{(u,v)\∈N_s}{\mathrm{\Σ}}{(W(u,v)-m_k)}^2-{\mathrm{\σ}}_n^2)$

其中N_s为以点(i，j)为中心，形状与椭圆窗口P相同的窗口，M_s为窗口N_s内重要系数的个数，m_k为窗口N_s内所有重要系数的均值，(u，v)是窗口N_s内任意点的坐标，W(u，v)是高频子带在点(u，v)处的系数值，为点(i，j)处的信号方差，为高频子带的噪声方差；

(4b)根据步骤(2)中高频子带系数的分割，若高频子带中的点(i，j)为重要系数点，按下式进行维纳滤波，否则进入步骤(4c)：

$x(i,j)=m_k+\frac{{\mathrm{\σ}}_k^2(i,j)}{{\mathrm{\σ}}_k^2(i,j)+{\mathrm{\σ}}_n^2}(W(i,j)-m_k)$

其中W(i，j)为高频子带在点(i，j)处的系数值，m_k为窗口N_s内所有重要系数的均值，为点(i，j)处的信号方差，为高频子带的噪声方差，x(i，j)为点(i，j)处滤波后的系数值；

(4c)若高频子带中的点(i，j)为非重要系数点，则按下式进行维纳滤波：

$x(i,j)=\frac{{\mathrm{\σ}}_k^2(i,j)}{{\mathrm{\σ}}_k^2(i,j)+{\mathrm{\σ}}_n^2}W(i,j)$

其中为点(i，j)处的信号方差，为高频子带的噪声方差，W(i，j)为高频子带在点(i，j)处的系数值，x(i，j)为点(i，j)处滤波后的系数值；

(5)对32个经过维纳滤波的高频子带进行轮廓波重构，得到大小为A×B的去噪图像FI；

(6)采用绝对值中值法估计去噪图像FI的噪声noise_FI，然后对大小为A×B的去噪图像FI进行非局部均值滤波，得到最终的去噪图像SI。

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

(a)本发明采用抗噪性良好的二维Otsu分割方法将高频系数分为重要系数和非重要系数，用重要系数估计信号方差，得到了准确的滤波因子，减少了去噪图像的噪声残留；

(b)本发明对重要系数和非重要系数进行了不同的维纳滤波，更好的保留了图像的细节信息，提高了去噪图像的清晰度。

附图说明

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

图2是本发明仿真使用的未加噪图像和加噪图像；

图3是采用三种对比方法以及本发明方法去噪后的图像。

具体实施方式

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

步骤1，对大小为A×B的含噪图像TI先进行轮廓波分解，分解层数为3，第一层分解得到16个高频子带，第二层分解得到8个高频子带，第三层分解得到8个高频子带，再采用绝对值中值法估计含噪图像TI的噪声noise_TI：

$\mathrm{noise}\_\mathrm{TI}=\frac{\mathrm{median}\left(\right|I\left|\right)}{0.6745}$

其中I为含噪图像TI进行轮廓波分解得到的第一层轮廓波高频系数，||是取绝对值操作，median是取中值操作。

步骤2，采用二维Otsu分别对以上分解出的各个高频子带进行分割，将高频子带中的系数分为重要系数和非重要系数：

(2a)对高频子带的系数取幅值得到幅值矩阵，利用二维Otsu对该幅值矩阵进行分割，得到两个阈值Th1，Th2，其中Th1≥Th2；

(2b)将高频子带中幅值大于阈值Th1的系数分为重要系数，其他系数分为非重要系数。

步骤3，用含噪图像TI的噪声noise_TI分别估计各个高频子带的噪声方差并计算各个高频子带的椭圆窗口P：

(3a)估计高频子带的噪声方差

${\mathrm{\σ}}_n^2={(\mathrm{noise}\_\mathrm{TI}*e^{1-n^{1/2}})}^2,$

其中noise_TI为含噪图像TI的噪声，n为子带所在层，n＝1，2，3；

(3b)计算高频子带的椭圆窗口P：

P＝{(a，b)∈Z|(acosθ+bsinθ)²+s⁴(-asinθ+bcosθ)²≤s²r²}

其中它为高频子带的窗口角度，n为高频子带所在分解层，n＝1，2，3，k表示高频子带在当前层中的子带标号，Z表示中心点为(0，0)，其尺寸为4r+1的正方形窗口，(a，b)为正方形窗口Z内任意点的坐标，(r，s)为窗口参数，其中r是尺寸参数，s是形状参数，第一层高频子带的窗口参数(r，s)选为(3，2)，第二层高频子带的窗口参数(r，s)选为(2，1.8)，第三层高频子带的窗口参数(r，s)选为(3，1.5)。

步骤4，对各个高频子带进行维纳滤波：

(4a)根据椭圆窗口P计算高频子带中点(i，j)处的信号方差：

${\mathrm{\σ}}_k^2(i,j)=\max (0,\frac{1}{M_s}\underset{(u,v)\∈N_s}{\mathrm{\Σ}}{(W(u,v)-m_k)}^2-{\mathrm{\σ}}_n^2),$

其中N_s为以点(i，j)为中心，且形状与椭圆窗口P相同的窗口，M_s为窗口N_s内重要系数的个数，m_k为窗口N_s内所有重要系数的均值，(u，v)是窗口N_s内任意点的坐标，W(u，v)是高频子带在点(u，v)处的系数值，为点(i，j)处的信号方差，为高频子带的噪声方差；

(4b)根据步骤(2)中高频子带系数的分割，若高频子带内点(i，j)为重要系数点，按下式进行维纳滤波，否则进入步骤(4c)：

$x(i,j)=m_k+\frac{{\mathrm{\σ}}_k^2(i,j)}{{\mathrm{\σ}}_k^2(i,j)+{\mathrm{\σ}}_n^2}(W(i,j)-m_k),$

其中W(i，j)为高频子带在点(i，j)处的系数值，m_k为窗口N_s内所有重要系数的均值，为点(i，j)处的信号方差，为高频子带的噪声方差，x(i，j)为点(i，j)处滤波后的系数值；

(4c)若高频子带中的点(i，j)为非重要系数点，则按下式进行维纳滤波：

$x(i,j)=\frac{{\mathrm{\σ}}_k^2(i,j)}{{\mathrm{\σ}}_k^2(i,j)+{\mathrm{\σ}}_n^2}W(i,j),$

其中为点(i，j)处的信号方差，为高频子带的噪声方差，W(i，j)为高频子带在点(i，j)处的系数值，x(i，j)为点(i，j)处滤波后的系数值。

步骤5，对32个经过维纳滤波的高频子带进行轮廓波逆变换，得到大小为A×B的去噪图像FI。

步骤6，采用绝对值中值法估计去噪图像FI的噪声noise_FI，然后对大小为A×B的去噪图像FI进行非局部均值滤波，得到最终的去噪图像SI：

(6a)图像FI的噪声noise_FI估计如下：

$\mathrm{noise}\_\mathrm{FI}=\frac{\mathrm{median}\left(\right|Y\left|\right)}{0.6745},$

Y为去噪图像FI小波分解三层时第一层的高频系数，||是取绝对值操作，median是取中值操作；

(6b)将大小为A×B的去噪图像FI分为若干个7×7的图像块，组成图像块集合{F_s}， 表示向下取整；

(6c)以图像块F_s为中心，取大小为21×21的搜索窗，由搜索窗的最上角开始，每隔3个像素取一个7×7的图像块，共得到49个相似图像块集合{G_q}，1≤q≤49；

(6d)利用相似图像块集合{G_q}对图像块F_s进行滤波：

$F_s^{\′}=\frac{1}{\underset{q\∈\mathrm{\Λ}}{\mathrm{\Σ}}{\mathrm{weight}}_q}\underset{q\∈\mathrm{\Λ}}{\mathrm{\Σ}}{\mathrm{weight}}_q{F}_s,$

其中F′_s为滤波后的图像块，Λ为集合{1，…，49}，weight_q为相似图像块G_q与图像块F_s之间的相似度：

${\mathrm{weight}}_q=e^{-\frac{1}{h}}{(F_s-G_q)}^2{e}^{-\frac{{(F_s-G_q)}^2}{2*\mathrm{noise}\_{\mathrm{FI}}^2}},$

其中G_q为相似图像块，1≤q≤49，noise_FI为去噪图像FI的噪声，h＝10noise_FI，是非局部均值滤波的平滑因子；

(6e)重复步骤(6c)到步骤(6d)，对图像块集合{F_s}中的图像块逐一进行滤波，得到最终的去噪图像SI。

本发明的效果可以进一步通过以下内容进行说明：

1.实验数据

实验所使用的输入图像如图2所示，图2(a)为Lena图像，图2(b)为Barara图像，图2(c)为Boats图像，图2(d)为加入噪声水平为20的高斯随机白噪声的lena图像，四幅图像均为灰度图像，大小均为512×512，灰度级为256。

2.实验评价指标

图像去噪效果的评价分为主观和客观两个方面。在主观上评价一幅图像去噪效果的优劣主要是通过人眼的视觉特性来衡量，图像质量好，感觉清晰则去噪效果好，反之效果则差。在客观上评价一幅图像的去噪效果，本发明采用峰值信噪比PSNR来衡量。令真实图像为V，去噪图像为V′，V(t)为真实图像V中第t个像素的幅值，V′(t)为去噪图像V′中第t个像素的幅值，Ω为图像中的像素集合，令V_max＝max{V(t)，t∈Ω}。则峰值信噪比按如下公式计算：

$\mathrm{PSNR}=10\log \left[\frac{V_{\max }^2}{\frac{1}{\left|\mathrm{\Ω}\right|}\underset{t\∈\mathrm{\Ω}}{\mathrm{\Σ}}{[V^{\′}\left(t\right)-V\left(t\right)]}^2}\right].$

3.实验内容

对比方法1：ZF Zhou等在文献“Contourlet-based Image Denoising Algorithm usingAdaptive Windows.ICIEA，2009：3654-3657.”中提出的轮廓波域自适应窗口维纳滤波去噪方法，简称C_W，利用该方法对图2(d)所示噪声水平为20的加噪lena图进行去噪，结果如图3(a)；

对比方法2：X Li，X.Q.Su等在文献“Image Denoising via Doubly Wiener Filteringwith Adaptive Directional Windows and Mean Shift Algorithm in Wavelet Domain.ICMA，2010：114-118.”中提出的小波域自适应窗口双维纳滤波去噪方法，简称W_DW，利用该方法对图2(d)所示噪声水平为20的加噪lena图进行去噪，结果如图3(b)；

对比方法3：Q Zhao等在文献“Q Zhao，B Ye，J.L.Cao.Image Denoising Based onImproved Non-local Means and Nonsubsampled Contourlet Transform Wiener Filtering.Journal of Computational Information Systems，2010，6(2)：601-610.”中提出的基于非局部均值的非下采样轮廓波域维纳滤波去噪方法，简称NSCT_W，利用该方法对图2(d)所示噪声水平为20的加噪lena图进行去噪，结果如图3(c)；

本发明方法：用本发明中提出的基于二维Otsu的轮廓波域维纳滤波图像去噪方法对图2(d)所示噪声水平为20的加噪lena图进行去噪，结果如图3(d)。

4.实验结果分析

从图3可以看出，图3(a)所示对比方法一的去噪结果噪声残留非常多；图3(b)所示对比方法二的去噪结果含有明显的划痕，伪Gibbs现象严重且噪声有残留；图3(c)所示对比方法三的去噪结果含有轻微的划痕，同质区域噪声有残留；图3(d)所示本发明方法的去噪结果没有划痕，噪声没有残留并且细节保持的也很好。可见，本发明的主观去噪效果优于现有三种对比方法。

表1列出了噪声水平从15到50时三种对比方法和本发明方法这四种去噪方法的峰值信噪比PSNR，表中sigma表示噪声水平，粗体部分为4组数据最大值。

表1  四种去噪方法对三幅图像在各种噪声水平下的峰值信噪比PSNR

由表1可以看出，所有噪声水平下本发明方法的PSNR都比三种对比方法的PSNR高，在客观上充分验证了本发明对自然图像去噪的有效性。

Claims (6)

1.一种基于二维Otsu的轮廓波域维纳滤波图像去噪方法，包括如下步骤：

（1）对大小为A×B的含噪图像TI先进行轮廓波分解，分解层数为3，第一层分解得到16个高频子带，第二层分解得到8个高频子带，第三层分解得到8个高频子带，再采用绝对值中值法估计含噪图像TI的噪声noise_TI；

（2）采用二维Otsu分别对以上分解出的各个高频子带进行分割，将高频子带中的系数分为重要系数和非重要系数；

（3）用含噪图像TI的噪声noise_TI分别估计各个高频子带的噪声方差并计算各个高频子带的椭圆窗口P；

（4）对各个高频子带进行维纳滤波：

（4a）根据椭圆窗口P计算高频子带中点(i,j)处的信号方差：

${\mathrm{\σ}}_k^2(i,j)=\max (0,\frac{1}{M_s}\underset{(u,v)\∈N_s}{\mathrm{\Σ}}{(W(u,v)-m_k)}^2-{\mathrm{\σ}}_n^2)$

其中N_s为以点(i,j)为中心，形状与椭圆窗口P相同的窗口，M_s为窗口N_s内重要系数的个数，m_k为窗口N_s内所有重要系数的均值，(u,v)是窗口N_s内任意点的坐标，W(u,v)是高频子带在点(u,v)处的系数值，为点(i,j)处的信号方差，为高频子带的噪声方差；

（4b）根据步骤（2）中高频子带系数的分割，若高频子带中的点(i,j)为重要系数点，按下式进行维纳滤波，否则进入步骤（4c）：

$x(i,j)=m_k+\frac{{\mathrm{\σ}}_k^2(i,j)}{{\mathrm{\σ}}_k^2(i,j)+{\mathrm{\σ}}_n^2}(W(i,j)-m_k)$

其中W(i,j)为高频子带在点(i,j)处的系数值，m_k为窗口N_s内所有重要系数的均值，为点(i,j)处的信号方差，为高频子带的噪声方差，x(i,j)为点(i,j)处滤波后的系数值；

（4c）若高频子带中的点(i,j)为非重要系数点，则按下式进行维纳滤波：

$x(i,j)=\frac{{\mathrm{\σ}}_k^2(i,j)}{{\mathrm{\σ}}_k^2(i,j)+{\mathrm{\σ}}_n^2}W(i,j)$

其中为点(i,j)处的信号方差，为高频子带的噪声方差，W(i,j)为高频子带在点(i,j)处的系数值，x(i,j)为点(i,j)处滤波后的系数值；

（5）对32个经过维纳滤波的高频子带进行轮廓波逆变换，得到大小为A×B的去噪图像FI；

（6）采用绝对值中值法估计去噪图像FI的噪声noise_FI，然后对大小为A×B的去噪图像FI进行非局部均值滤波，得到最终的去噪图像SI。

2.根据权利要求1所述的图像去噪方法，步骤（1）所述的采用绝对值中值法估计含噪图像TI的噪声noise_TI，按如下公式进行：

$\mathrm{noise}\_\mathrm{TI}=\frac{\mathrm{median}\left(\right|I\left|\right)}{0.6745}$

其中，I为含噪图像TI进行轮廓波分解得到的第一层轮廓波高频系数，||是取绝对值操作，median是取中值操作。

3.根据权利要求1所述的图像去噪方法，步骤（2）中所述的采用二维Otsu分别对以上分解出的各个高频子带进行分割，按如下步骤进行：

（2a）对高频子带的系数取幅值得到幅值矩阵，利用二维Otsu对该幅值矩阵进行分割，得到两个阈值Th1，Th2，其中Th1≥Th2；

（2b）将高频子带中幅值大于阈值Th1的系数作为重要系数，其他系数作为非重要系数。

4.根据权利要求1所述的图像去噪方法，步骤（3）中所述的用含噪图像TI的噪声noise_TI分别估计各个高频子带的噪声方差是通过如下公式进行：

${\mathrm{\σ}}_n^2={(\mathrm{noise}\_\mathrm{TI}*e^{{1-n}^{1/2}})}^2$

其中noise_TI为含噪图像TI的噪声，n为高频子带所在分解层，n=1,2,3。

5.根据权利要求1所述的图像去噪方法，步骤（3）中所述的计算各个高频子带的椭圆窗口P，是通过如下公式进行：

P={(a,b)∈Z|(acosθ+bsinθ)²+s⁴(-asinθ+bcosθ)²≤s²r²}

其中θ为高频子带的窗口， $\mathrm{\θ}=(\frac{3}{4}-\frac{1}{2^{\mathrm{\ν}+1}}-\frac{(\mathrm{\κ}-1)}{2^{\mathrm{\ν}}})$

n为高频子带所在分解层，n=1,2,3；

k表示高频子带在所在分解层中的子带标号；

Z表示中心点为（0,0），其尺寸为4r+1的正方形窗口；

（a,b）为正方形窗口Z内任意点的坐标；

（r,s）为窗口参数，其中r是尺寸参数，s是形状参数，第一层高频子带的窗口参数（r,s）为（3,2），第二层高频子带的窗口参数（r,s）为（2,1.8），第三层高频子带的窗口参数（r,s）为（3,1.5）。

6.根据权利要求1所述的图像去噪方法，步骤（6）中所述的对大小为A×B的去噪图像FI进行非局部均值滤波，按如下步骤进行：

（6_a）将大小为A×B的去噪图像FI分为若干个7×7的图像块集合{F_s}，1≤s≤∟表示向下取整；

（6b）以图像块F_s为中心，取大小为21×21的搜索窗，由搜索窗的最上角开始，每隔3个像素取一个7×7的图像块，共得到49个相似图像块集合{G_q}，1≤_q≤49；

（6c）利用相似图像块集合{G_q}对图像块F_s进行滤波：

$F_s^{\′}=\frac{1}{\underset{q\∈\mathrm{\Λ}}{\mathrm{\Σ}}{\mathrm{weight}}_q}\underset{q\∈\mathrm{\Λ}}{\mathrm{\Σ}}{\mathrm{weight}}_q{F}_s$

其中F_s'为滤波后的图像块，Λ为集合{1,…,49}，weight_q为相似图像块Gq与图像块F_s之间的相似度：

${\mathrm{weight}}_q=e^{-\frac{1}{h}}{(F_s-G_q)}^2{e}^{-\frac{{(F_s-G_q)}^2}{{2*\mathrm{noise}\_\mathrm{FI}}^2}}$

其中G_q为相似图像块，1≤q≤49，noise_FI为去噪图像FI的噪声，h=10noise_FI，是非局部均值滤波的平滑因子；

（6d）重复步骤（6b）到步骤（6c），对图像块集合{F_s}中的图像块逐一进行滤波，得到最终的去噪图像SI。