CN103198455A - 一种利用全变差最小化和灰度共生矩阵的图像去噪方法 - Google Patents

Abstract

本发明的目的在于提供一种利用全变差最小化和灰度共生矩阵的图像去噪方法，包括以下步骤：对原始含噪图像进行高斯滤波，使用检测窗遍历得到的图像，求出每个检测窗内子图像块的四个灰度共生矩阵，由得到的灰度共生矩阵求对比度图像，利用得到的对比度图像，并结合全变差最小化模型及各项扩散模型去除原始含噪图像中的噪声干扰。本发明提高了对边缘等纹理信息位置的检测精度，且使用对比度图像来自适应的在全变差最小化去噪方法和各项同性扩散去噪方法之间过渡，兼顾了二者在去噪和保护边缘的方面优点，并能有效地减少阶梯效应的影响。

Description

一种利用全变差最小化和灰度共生矩阵的图像去噪方法

技术领域

本发明涉及的是一种图像处理方法。

背景技术

在图像成像及传输过程中，难免会受到各种噪声的干扰，给图像的后续使用和解析带来很多不利影响。因此，图像去噪是数字图像处理的基本任务之一，是图像特征提取、分割和模式识别等图像操作的前提。

近年来，RUDIN等人提出的基于全变差最小化的图像去噪方法因其在去噪的同时能够有效地保护图像边缘而在图像去噪领域得到广泛的研究及应用，受到越来越多研究者的重视。通常含噪图像的全变差比不含噪图像的全变差明显大，全变差最小化去噪的基本思想就是通过减少图像的全变差以达到去噪的目的。其实质是先将图像去噪问题转化为一个泛函求极值问题，再使用变分法导出一个欧拉-拉格朗日方程。然后在给定初始条件和边界条件下，通过数值计算求解出去噪图像。但该去噪模型易将噪声所在位置当做图像边缘的来处理，在图像平滑区域产生阶梯效应。

为了解决全变差最小化去噪模型的不足，很多学者对该模型提出了许多改进方法。Chambolle提出一种自适应全变差去噪模型，通过比较梯度模值和给定阈值之间的大小关系来判断是使用全变差最小化模型还是各项同性扩散模型来对图像进行去噪处理。其实质是在利用全变差最小化去噪模型能够保护图像边缘性质的同时，在图像非边缘区域使用各项同性扩散去噪模型抑制噪声，取得了不错的去噪效果。但是，该模型对阈值的选取较为敏感，而且使用梯度模值来判断图像边缘的位置精度不高，容易误判。

因此，如何利用含噪图像自身的信息，找到一种切实有效的方法来检测图像边缘等纹理信息，将含噪图像的纹理区域和平滑区域区分开，并使用一定的技术手段在含噪图像上两种不同区域上自适应的选择全变差最小化去噪方法或各项同性扩散去噪方法来进行处理，是很多研究人员不断努力的目标。

发明内容

本发明的目的在于提供在尽可能不破坏图像原始有用信息和完整性的前提下，在去除噪声的同时很好地保护边缘，减少或者消除阶梯效应的影响的一种利用全变差最小化和灰度共生矩阵的图像去噪方法。

本发明的目的是这样实现的：

本发明一种利用全变差最小化和灰度共生矩阵的图像去噪方法，其特征是：

（1）对原始含噪图像进行高斯滤波：

设含噪图像为X，其大小为M×N，灰度级范围为[0,255]，用高斯滤波器对图像X进行预处理，去掉非边缘区域孤立的噪声点，其中高斯滤波器的窗口大小为G×G、方差为σ，经过高斯滤波后的得到的图像记为X′；

（2）使用检测窗遍历由步骤（1）得到的图像，求出每个检测窗内子图像块的四个灰度共生矩阵：

1）将图像X′的灰度级由256降为32：

则X′的灰度值范围变为[1,32]；

2）选取大小为M_x×M_y的检测窗口在图像X′上沿水平和垂直方向上移动，每次移动1个像素距离，截取出大小为M_x×M_y的子图像块，并记在图像X′上点(i,j)处截取的子图像块为X_i,j，(i＝1,?,M,j＝1,?,N)；

在使用检测窗口移动前，将图像在左侧、右侧、顶部和底部四个方向分别依次扩展

行和

列：

在顶部和底部：

在左侧和右侧：

其中

表示向下取整，M、N表示原始图像的行和列，M_x、M_y表示检测窗口的行和列，得到新的X′，有X′＝X′′′，X′的大小变为

3）计算上述每个子图像块四个方向上的灰度共生矩阵：子图像块X_i,j生成的灰度共生矩阵记为表示从灰度值为m的点到灰度值为n的点的概率，其中灰度值为m，n两点间距离为d，两点连线与x轴的夹角为θ，取θ=0°,45°,90°,135°四个方向，则得到子图像块X_i,j的四个灰度共生矩阵；

（3）由步骤（2）中得到的灰度共生矩阵求对比度图像：

定义一个大小为M×N的零矩阵CON，对以上求出的每个图像子块X_i,j(i＝1,...,M,j＝1,...,N)的四个方向上的灰度共生矩阵

θ=0°,45°,90°,135°，分别计算其对比度，记为Con1(i,j),Con2(i,j),Con3(i,j),Con4(i,j)，其中对比度定义为：

$\mathrm{Con}=\underset{m=1}{\overset{L-1}{\mathrm{\Σ}}}\underset{n=1}{\overset{L-1}{\mathrm{\Σ}}}{(m-n)}^2{P}_d^{\mathrm{\θ}}(m,n)$

其中L＝32为子图像块的灰度级，然后求出四个对比度的平均值并赋值给CON矩阵中相应的位置，即：

CON(i,j)＝(Con1(i,j)+Con2(i,j)+Con3(i,j)+Con4(i,j))/4

求出的矩阵CON即为由灰度共生矩阵得到对比度特征图像；

（4）利用步骤（3）得到的对比度图像，并结合全变差最小化模型及各项扩散模型去除原始含噪图像中的噪声干扰：

去噪模型为：

$\underset{u\∈\mathrm{BV}\left(\mathrm{\Ω}\right)}{\min }(1-\mathrm{\φ}\left(X\right)){\∫\∫}_{\mathrm{\Ω}}|\▿u|\mathrm{dxdy}+\mathrm{\φ}\left(X\right){\∫\∫}_{\mathrm{\Ω}}{|\▿u|}^2\mathrm{dxdy}+\frac{\mathrm{\λ}}{2}{\∫\∫}_{\mathrm{\Ω}}{(u-X)}^2\mathrm{dxdy}$

其中u为清晰的原始信号，Ω为信号的积分区域，BV代表有界变差空间，λ表示用以平衡去噪模型前两项和第三项权重的拉格朗日乘子，φ(X)＝1/1+CON，▽u＝(u_x,u_y)，表示u的梯度模值；

去噪模型的欧拉-拉格朗日方程可表示为：

$-(1-\mathrm{\φ}\left(X\right))\mathrm{div}\left(\frac{\▿u}{|\▿u|}\right)-2\mathrm{\φ}\left(X\right)\▿u+\mathrm{\λ}(u-X)=0$

其中

表示散度算子，

表示拉普拉斯算子；

采用人工时间演化方法对欧拉-拉格朗日方程进行求解，其数值计算形式为：

$u^{n+1}=u^n+d_t(1-\mathrm{\φ}\left(X\right))\frac{u_{\mathrm{xx}}^n{\left(u_y^n\right)}^2+u_{\mathrm{yy}}^n{\left(u_x^n\right)}^2-u_x^n{u}_y^n{u}_{\mathrm{xy}}^n-u_x^n{u}_y^n{u}_{\mathrm{yx}}^n}{{({\left(u_x^n\right)}^2+{\left(u_y^n\right)}^2)}^{\frac{3}{2}}}+2\mathrm{\φ}\left(X\right)(u_{\mathrm{xx}}^n-u_{\mathrm{yy}}^n)-\mathrm{\λ}(u^n-X)$

其中d_t表示时间迭代步长；

将上式中的u⁰初始化为X并按上式迭代，即可得到去噪后的图像。

本发明的优势在于：

（1）本发明使用由灰度共生矩阵得到的对比度图像来区分含噪图像的纹理区域和光滑区域，充分利用了灰度共生矩阵的对比度特征图像灰度值在反映图像平滑性方面的优势，提高了对边缘和光滑区域的区分能力；

（2）将本发明与Chambolle的去噪方法相比较，使用对比度图像来自适应的在全变差最小化去噪方法和各项同性扩散去噪方法之间过渡，而不是根据给定的阈值使用梯度模值来硬性决定使用全变差最小化去噪方法和各项同性扩散去噪方法之一，提高了鲁棒性。

（3）本发明所述去噪模型兼顾了全变差最小化去噪方法和各项同性扩散去噪方法在去噪和保护边缘方面的优点，有效地减少了阶梯效应的影响。

附图说明

图1为本发明的流程图；

图2为本发明的由检测窗获取子图像块示意图；

图3为本发明的灰度共生矩阵像素对关系示意图；

图4a为加入噪声后的lena图像；图4b为使用Chambolle的方法对含噪lena图像的去噪效果；图4c为使用本发明方法对含噪lena图像的去噪效果；

图5a为加入噪声后的cameraman图像；图5b为使用Chambolle的方法对含噪cameraman图像的去噪效果；图5c为使用本发明方法对含噪cameraman图像的去噪效果。

具体实施方式

下面结合附图举例对本发明做更详细地描述：

结合图1～5，本发明的利用全变差最小化和灰度共生矩阵的图像去噪方法实现流程如图1所示，包括以下几个步骤：

步骤一、对原始含噪图像进行高斯滤波。

设含噪图像为X，其大小为M×N，灰度级范围为[0,255]。用高斯滤波器对图像X进行预处理，去掉非边缘区域一些孤立的噪声点，减少将这些噪声当做边缘的可能性。其中，高斯滤波器的窗口大小为G×G，方差为σ。经过高斯滤波后的得到的图像记为X′。

步骤二、使用检测窗遍历由步骤一得到的图像，求出每个检测窗内子图像块的四个灰度共生矩阵；

1）将图像X′的灰度级由256降为32，具体操作如下式所示。

（3）

则X′的灰度值范围变为[1,32]。

2）选取大小为M_x×M_y的检测窗口在图像X′上沿水平和垂直方向上移动。每次移动1个像素距离，截取出大小为M_x×M_y的子图像块，如图2所示，并记在图像X′上点(i,j)处截取的子图像块为X_i,j，(i＝1,...,M,j＝1,...，N)。由于检测窗口在X′的左侧、右侧、顶部和底部超出了图像X′的范围，因此本发明给出的补偿办法是在使用检测窗口移动前，将图像在左侧、右侧、顶部和底部四个方向分别依次扩展

行和列，具体做法如下：

在顶部和底部：

在左侧和右侧：

（5）

上述

表示向下取整，M、N表示原始图像的行和列，M_x、M_y表示检测窗口的行和列。由以上方法得到新的X′，有X′＝X′′′。且经过式（4）和式（5）的操作，X′的大小变为

3）计算上述每个子图像块四个方向上的灰度共生矩阵。子图像块X_i,j生成的灰度共生矩阵记为表示从灰度值为m的点到灰度值为n的点的概率，如图3所示。其中灰度值为m，n两点间距离为d，两点连线与x轴的夹角为θ。取θ=0°,45°,90°,135°四个方向，则得到子图像块X_i,j的四个灰度共生矩阵。

步骤三、由步骤二中得到的灰度共生矩阵求对比度图像。

定义一个大小为M×N的零矩阵CON。对以上求出的每个图像子块X_i,j(i＝1,...,M,j＝1,...,N)的四个方向上的灰度共生矩阵

θ=0°,45°,90°,135°，分别计算其对比度，记为Con1(i,j),Con2(i,j),Con3(i,j),Con4(i,j)。其中对比度定义为：

$\mathrm{Con}=\underset{m=1}{\overset{L-1}{\mathrm{\Σ}}}\underset{n=1}{\overset{L-1}{\mathrm{\Σ}}}{(m-n)}^2{P}_d^{\mathrm{\θ}}(m,n)$ （6）

其中，L＝32为子图像块的灰度级。然后求出四个对比度的平均值并赋值给CON矩阵中相应的位置，即：

CON(i,j)＝(Con1(i,j)+Con2(i,j)+Con3(i,j)+Con4(i,j))/4（7）

求出的矩阵CON即为由灰度共生矩阵得到对比度特征图像。

步骤四、利用步骤三得到的对比度图像，并结合全变差最小化模型及各项扩散模型去除原始含噪图像中的噪声干扰。

可采用人工时间演化方法对式(2)进行求解，其数值计算形式为：

$u^{n+1}=u^n+d_t(1-\mathrm{\φ}\left(X\right))\frac{u_{\mathrm{xx}}^n{\left(u_y^n\right)}^2+u_{\mathrm{yy}}^n{\left(u_x^n\right)}^2-u_x^n{u}_y^n{u}_{\mathrm{xy}}^n-u_x^n{u}_y^n{u}_{\mathrm{yx}}^n}{{({\left(u_x^n\right)}^2+{\left(u_y^n\right)}^2)}^{\frac{3}{2}}}+2\mathrm{\φ}\left(X\right)(u_{\mathrm{xx}}^n-u_{\mathrm{yy}}^n)-\mathrm{\λ}(u^n-X)$

（8）

其中d_t表示时间迭代步长，φ(X)＝1/1+CON，CON为由步骤三中得到对比度特征图像。为了防止由于CON的值过小使得φ(X)难以逼近0，本发明具体实验中将CON的值扩大2倍。且为了避免(8)式分母中不出现为零的情况，引入一个小的正参数ε，使得(8)式分母变为

将式(8)中的u⁰初始化为X并按式(8)迭代N次，即可得到去噪后的图像。

本发明的效果可以通过以下实验进行验证，并与Chambolle方法的去噪效果进行比较。

（1）实验条件

实验仿真环境为：Matlab R2012a CPU Pentium(R)2.5GHz，2G内存；实验采用图像为：lena和cameraman灰度图像，大小均为256×256，灰度值的变化范围为[0,255]；实验参数为：对两幅图像分别加入方差为20的高斯白噪声，高斯滤波器的窗口大小为3×3，方差为1，检测窗口大小为5×5（M_x×M_y），迭代次数N=80，λ=0.2，时间迭代步长d_t=0.01，小的正参数ε=1。

（2）实验内容

（2.1）用本发明和Chambolle的方法对lena图像进行实验。图4(a)为加入噪声后的lena图像，图4(b)为使用Chambolle的方法对含噪lena图像的去噪效果，图4(c)为使用本发明方法对含噪lena图像的去噪效果；

（2.2）用本发明和Chambolle的方法对cameraman图像进行实验。图5(a)为加入噪声后的cameraman图像，图5(b)为使用Chambolle的方法对含噪cameraman图像的去噪效果，图5(c)为使用本发明方法对含噪cameraman图像的去噪效果。

（3）实验结果分析

从图4和图5所示去噪的视觉效果图可以看到，Chambolle方法去噪效果不是很好，在光滑区域会产生阶梯效应，而本发明在边缘上得到的效果比较好，且有效地抑制了阶梯效应的产生。

表1给出了本发明方法和Chambolle方法的去噪效果统计分析。表1中的数据为去噪图像与原始清晰图像对比所计算出来的峰值信噪比，表示图像质量的改善情况。

表1本发明和Chambolle方法去噪结果的峰值信噪比比较

从表1中可知，本发明去噪结果与Chambolle方法去噪结果相比较，图像质量明显改善。

结合对图4和图5的视觉效果以及对表1的统计分析，可知本发明的利用全变差最小化和灰度共生矩阵的图像去噪方法，通过使用由灰度共生矩阵得到的对比度图像来区分含噪图像的纹理区域和光滑区域，相对于根据给定的阈值使用梯度模值来区分这两种区域，提高了对边缘等纹理信息位置的检测精度。且使用对比度图像来自适应的在全变差最小化去噪方法和各项同性扩散去噪方法之间过渡，兼顾了二者在去噪和保护边缘的方面优点，并能有效地减少阶梯效应的影响。

Claims (1)

1.一种利用全变差最小化和灰度共生矩阵的图像去噪方法，其特征是：

（1）对原始含噪图像进行高斯滤波：

设含噪图像为X，其大小为M×N，灰度级范围为[0,255]，用高斯滤波器对图像X进行预处理，去掉非边缘区域孤立的噪声点，其中高斯滤波器的窗口大小为G×G、方差为σ，经过高斯滤波后的得到的图像记为X′；

（2）使用检测窗遍历由步骤（1）得到的图像，求出每个检测窗内子图像块的四个灰度共生矩阵：

1）将图像X′的灰度级由256降为32：

则X′的灰度值范围变为[1,32]；

2）选取大小为M_x×M_y的检测窗口在图像X′上沿水平和垂直方向上移动，每次移动1个像素距离，截取出大小为M_x×M_y的子图像块，并记在图像X′上点(i，j)处截取的子图像块为X_i,j，(i＝1,...,M,j＝1,...,N)；

在使用检测窗口移动前，将图像在左侧、右侧、顶部和底部四个方向分别依次扩展

行和

列：

在顶部和底部：

在左侧和右侧：

其中

表示向下取整，M、N表示原始图像的行和列，M_x、M_y表示检测窗口的行和列，得到新的X′，有X′＝X′′′，X′的大小变为

3）计算上述每个子图像块四个方向上的灰度共生矩阵：子图像块X_i,j生成的灰度共生矩阵记为

表示从灰度值为m的点到灰度值为n的点的概率，其中灰度值为m，n两点间距离为d，两点连线与x轴的夹角为θ，取θ=0°,45°,90°,135°四个方向，则得到子图像块X_i,j的四个灰度共生矩阵；

（3）由步骤（2）中得到的灰度共生矩阵求对比度图像：

定义一个大小为M×N的零矩阵CON，对以上求出的每个图像子块X_i,j(i＝1,...,M,j＝1,...,N)的四个方向上的灰度共生矩阵

θ=0°,45°,90°,135°，分别计算其对比度，记为Con1(i,j),Con2(i,j),Con3(i,j),Con4(i,j)，其中对比度定义为：

$\mathrm{Con}=\underset{m=1}{\overset{L-1}{\mathrm{\Σ}}}\underset{n=1}{\overset{L-1}{\mathrm{\Σ}}}{(m-n)}^2{P}_d^{\mathrm{\θ}}(m,n)$

其中L＝32为子图像块的灰度级，然后求出四个对比度的平均值并赋值给CON矩阵中相应的位置，即：

CON(i,j)＝(Con1(i,j)+Con2(i,j)+Con3(i,j)+Con4(i,j))/4

求出的矩阵CON即为由灰度共生矩阵得到对比度特征图像；

（4）利用步骤（3）得到的对比度图像，并结合全变差最小化模型及各项扩散模型去除原始含噪图像中的噪声干扰：

去噪模型为：

$\underset{u\∈\mathrm{BV}\left(\mathrm{\Ω}\right)}{\min }(1-\mathrm{\φ}\left(X\right)){\∫\∫}_{\mathrm{\Ω}}|\▿u|\mathrm{dxdy}+\mathrm{\φ}\left(X\right){\∫\∫}_{\mathrm{\Ω}}{|\▿u|}^2\mathrm{dxdy}+\frac{\mathrm{\λ}}{2}{\∫\∫}_{\mathrm{\Ω}}{(u-X)}^2\mathrm{dxdy}$

其中u为清晰的原始信号，Ω为信号的积分区域，BV代表有界变差空间，λ表示用以平衡去噪模型前两项和第三项权重的拉格朗日乘子，φ(X)＝11+CON，▽u＝(u_x,u_y)，

表示u的梯度模值；

去噪模型的欧拉-拉格朗日方程可表示为：

$-(1-\mathrm{\φ}\left(X\right))\mathrm{div}\left(\frac{\▿u}{|\▿u|}\right)-2\mathrm{\φ}\left(X\right)\▿u+\mathrm{\λ}(u-X)=0$

其中

表示散度算子，

表示拉普拉斯算子；

采用人工时间演化方法对欧拉-拉格朗日方程进行求解，其数值计算形式为：

$u^{n+1}=u^n+d_t(1-\mathrm{\φ}\left(X\right))\frac{u_{\mathrm{xx}}^n{\left(u_y^n\right)}^2+u_{\mathrm{yy}}^n{\left(u_x^n\right)}^2-u_x^n{u}_y^n{u}_{\mathrm{xy}}^n-u_x^n{u}_y^n{u}_{\mathrm{yx}}^n}{{({\left(u_x^n\right)}^2+{\left(u_y^n\right)}^2)}^{\frac{3}{2}}}+2\mathrm{\φ}\left(X\right)(u_{\mathrm{xx}}^n-u_{\mathrm{yy}}^n)-\mathrm{\λ}(u^n-X)$

其中d_t表示时间迭代步长；

将上式中的u⁰初始化为X并按上式迭代，即可得到去噪后的图像。

Images