CN105427262A - 基于双向增强扩散滤波的图像去噪方法 - Google Patents

Abstract

本发明公开了一种基于双向增强扩散滤波的图像去噪方法，简化了扩散方程，建立双向扩散系数，使模型在扩散过程中能够实现平滑与锐化的双向过程，为进一步加强平滑和锐化强度，对图像做增强处理，并用小波变换增强图像的整体轮廓，弱化图像的纹理细节，然后，对梯度阈值进行了自适应设计和改进，使其根据图像的最大灰度值和迭代次数自动控制梯度阈值，进一步保留图像边缘和细节特征，最后，对提出的模型进行仿真，用MATLAB软件对方法进行仿真验证，能够兼顾图像噪声的去除和边缘、纹理等细节信息的保护，峰值信噪比有了大幅提高，去噪性能较经典模型更具优越性，具有良好的应用前景。

Description

基于双向增强扩散滤波的图像去噪方法

技术领域

本发明涉及图像去噪技术领域，具体涉及一种基于双向增强扩散滤波的图像去噪方法。

背景技术

在图像处理与计算机领域中，图像去噪是最基本的问题之一。当今，许多科学领域都包含着数字图像处理技术，而基于偏微分方程的图像去噪方法中，各向异性扩散已成为当前研究的热点。Perona和Malik在1990年首次提出了PM模型，该模型是令传导系数依赖于图像梯度的二阶偏微分模型。1992年，Rudin，Osher，Fatemi等从能量泛函的角度提出了基于图像全变差的正则化TV模型，该模型能更好地保持图像的边缘纹理等细节特征，是偏微分方程在图像去噪中的重要体现。GuyGilboa等人在2002年提出了复模型，该模型能够同时实现前后向扩散，将PM模型从实数域扩展到了复数域，实现了图像去噪范围的扩大。自PM模型提出后，近些年，基于偏微分方程的图像去噪方法取得了很多有价值的成果。

最近，从另一些角度考虑，出现了许多较新的图像去噪算法，A.Buades等人提出了NLM滤波器，该滤波器较好的的解决了滤波中保持图像纹理和细节特征的问题，但是计算复杂，处理速度慢，对原图像的结构信息保护还不够好；BM3D方法去噪后的图像不仅有较高的信噪比，视觉效果也很好，但是时间复杂度相对较高；C.A.Deledalle等人提出了PPB(ProbablisticPatch‐based)滤波器，该滤波器对合成孔径雷达(syntheticapertureradar,SAR)图像滤波效果较好，但该滤波器在滤波过程中耗时过长，而其非迭代的滤波方法虽然解决了耗时过长的问题，但该滤波器不能有效地保持图像的纹理和细节特征。

如何设计一种图像去噪算法，便于实现，且能够兼顾图像噪声的去除和边缘、纹理等细节信息的保护，提高去噪性能，是当前急需解决的问题。

发明内容

本发明的目的是为了克服现有的图像去噪算法存在的不足。本发明的基于双向增强扩散滤波的图像去噪方法，能够兼顾图像噪声的去除和边缘、纹理等细节信息的保护，峰值信噪比有了大幅提高，去噪性能较经典模型更具优越性，具有良好的应用前景。

为了达到上述的目的，本发明所采用的技术方案是：

一种基于双向增强扩散滤波的图像去噪方法，其特征在于：包括以下步骤，

步骤(A)，建立双向扩散系数，根据PM模型的最小能量泛函，得到的扩散方程，如公式(1)所示，

$\left\{\begin{array}{c}\frac{\∂I}{\∂t}=div(\frac{f^{\′}\left(\right|\▿I\left|\right)}{|\▿I|}\·\▿I)\\ I(t=0)=I_0\end{array}\right.---\left(1\right)$

其中，I为图像，div、▽分别为图像I的散度算子和梯度算子，I₀为初始图像，f(|▽I|)为非增扩散函数、f'(|▽I|)为函数f(|▽I|)的一阶导数、I(t＝0)为图像在t＝0时的值，t为扩散时间尺度，扩散系数为关于梯度的函数，为经典PM模型下的扩散系数；

构建公式(1)内在的坐标形式，设图像中某点o(i,j)的内在坐标系为(η,ξ)，η为图像I的梯度方向；ξ为垂直于梯度的方向， $\mathrm{\η}=\frac{(I_x,I_y)}{\sqrt{I_x^2+I_y^2}},\mathrm{\ξ}=\frac{(-I_y,I_x)}{\sqrt{I_x^2+I_y^2}}$

则， $I_{\mathrm{\ξ}\mathrm{\ξ}}=\frac{I_y^2{I}_{xx}-2I_x{I}_y{I}_{xy}+I_x^2{I}_{yy}}{I_x^2+I_y^2},I_{\mathrm{\η}\mathrm{\η}}=\frac{I_x^2{I}_{xx}+2I_x{I}_y{I}_{xy}+I_y^2{I}_{yy}}{I_x^2+I_y^2}$

其中，x,y为两个正交方向，I_x为图像I在x方向分量，I_y为图像I在y方向分量，I_xx，I_xy分别为I_x在x,y方向的方向导数，I_yy为I_y在y方向的方向导数，I_ηη表示图像I沿着边缘梯度方向的方向导数，I_ξξ表示图像I沿着边缘切线方向的方向导数；

扩散方程的最终表达式，如公式(2)所示，

$\frac{\∂I}{\∂t}=\frac{f^{\′}\left(\right|\▿I\left|\right)}{2|\▿I|}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+\frac{f^{\′\′}\left(\right|\▿I\left|\right)}{2}{I}_{\mathrm{\η}\mathrm{\η}}---\left(2\right)$

为了使公式(1)有解，f(|▽I|)单调递增，f′(|▽I|)＞0；由公式(2)可知，扩散系数为沿图像的边缘方向，即沿着边缘方向进行正向扩散，使图像模糊；而扩散系数为沿图像的梯度方向，其正负号不确定，当f”(|▽I|)＞0时，沿梯度方向正扩散，使图像模糊，当f”(|▽I|)＜0时，沿梯度方向逆扩散，锐化图像的边缘，保护图像的边缘；

步骤(B)，突出图像的整体结构，增强边缘的锐化作用，对公式(2)做进一步化简，如公式(3)所示，

$\begin{array}{c}\frac{\∂I}{\∂t}=\frac{f^{\′}\left(\right|\▿I\left|\right)}{2|\▿I|}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+{\[|\▿I|\·\frac{f^{\′}\left(\right|\▿I\left|\right)}{2|\▿I|}\]}^{\′}{I}_{\mathrm{\η}\mathrm{\η}}\\ =g\left(\right|\▿I\left|\right)I_{\mathrm{\ξ}\mathrm{\ξ}}+{\[|\▿I|\·g\left(\right|\▿I\left|\right)\]}^{\′}{I}_{\mathrm{\η}\mathrm{\η}}\end{array}---\left(3\right)$

其中，为经典PM模型简化分析后的扩散系数；

令k是图像I的梯度阈值，当|▽I|→0时，g(|▽I|)→1；当|▽I|→∞时，g(|▽I|)→0，扩散系数的选取条件，将代入公式(3)，得到公式(4)，

$\begin{array}{c}\frac{\∂I}{\∂t}=\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+{\left(\frac{|\▿I|(1+k^2)}{|\▿I{|}^2+k^2}\right)}^{\′}{I}_{\mathrm{\η}\mathrm{\η}}\\ =\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+\left(\frac{(1+k^2)(k^2-|\▿I{|}^2)}{{\left(\right|\▿I{|}^2+k^2)}^2}\right)I_{\mathrm{\η}\mathrm{\η}}\\ =\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+\frac{(1+k^2)}{|\▿I{|}^2+k^2}\frac{(k^2-|\▿I{|}^2)}{|\▿I{|}^2+k^2}{I}_{\mathrm{\η}\mathrm{\η}}\end{array}---\left(4\right)$

当在图像边缘时，|▽I|＞＞k，得到公式(5)，

$\begin{array}{c}\frac{\∂I}{\∂t}\≈\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+\frac{(1+k^2)(0-|\▿I{|}^2)}{|\▿I{|}^2+0|\▿I{|}^2+k^2}{I}_{\mathrm{\η}\mathrm{\η}}\\ =\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}-\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\η}\mathrm{\η}}\\ =g\left(\right|\▿I\left|\right)I_{\mathrm{\ξ}\mathrm{\ξ}}-g\left(\right|\▿I\left|\right)I_{\mathrm{\η}\mathrm{\η}}\end{array}---\left(5\right)$

由公式(5)可知，图像沿梯度方向逆扩散，锐化图像的边缘，能够改善图像边缘模糊的现象，保护图像的边缘；

当在图像平坦区域时，|▽I|＜＜k，得到公式(6)，

$\begin{array}{c}\frac{\∂I}{\∂t}\≈\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+\frac{(1+k^2)}{0+k^2}\frac{(k^2-0)}{|\▿I{|}^2+k^2}{I}_{\mathrm{\η}\mathrm{\η}}\\ =\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\η}\mathrm{\η}}\\ =g\left(\right|\▿I\left|\right)I_{\mathrm{\ξ}\mathrm{\ξ}}+g\left(\right|\▿I\left|\right)I_{\mathrm{\η}\mathrm{\η}}\end{array}---\left(6\right)$

由公式(6)可知，图像的平坦区域，沿着图像的边缘方向和梯度方向同时平滑，具有去噪能力；

步骤(C)，根据图像的最大灰度值和迭代次数自动控制梯度阈值k，进一步保留图像边缘和细节特征，梯度阈值k为一个随扩散时间和扩散次数变化的函数，如公式(7)所示，

k(t)＝e^-αt(7)

其中，t为扩散时间，扩散次数n代表扩散时间，则t＝n，max{W,H}为图像的最大灰度值，随着扩散次数增加或者扩散时间增加，梯度阈值k不断减小，准确的缩小边缘判断的范围，保留图像边缘和细节；

在滤波过程中为了不让图像失真，给图像添加一项保真项，保真项，如公式(8)所示，

$D={\mathrm{\λ}}^{*}(1-\frac{1}{|\▿I{|}^2+1})(I-I_0)---\left(8\right)$

其中，λ^*是调整参数，在图像的边缘，此时|▽I|→∞，滤波结果会接近原始图像，加强了保边缘的效果；在图像的平坦区域，此时|▽I|→0，则扩散程度达到最大，去除噪声；

步骤(D)，根据公式(4)和公式(8)，建立双向增强的扩散模型，如公式(9)所示，

$\left\{\begin{array}{c}\frac{\∂I}{\∂t}=\frac{1+k{\left(t\right)}^2}{|\▿\hat{I}{|}^2+k{\left(t\right)}^2}{\hat{I}}_{\mathrm{\ξ}\mathrm{\ξ}}+\frac{(1+k{\left(t\right)}^2)(k{\left(t\right)}^2-|\▿\hat{I}{|}^2)}{|\▿\hat{I}{|}^2+k{\left(t\right)}^2|\▿\hat{I}{|}^2+k{\left(t\right)}^2}{\hat{I}}_{\mathrm{\η}\mathrm{\η}}\\ -{\mathrm{\λ}}^{*}(1-\frac{1}{|\▿I{|}^2+1})(I-I_0)\\ I(x,y,0)=I_0\end{array}\right.---\left(9\right)$

其中，分别表示为图像I是经过小波增强处理后沿着边缘切线方向的方向导数和沿着边缘梯度方向的方向导数，I(x,y,0)＝I₀表示I(x,y,t)在t＝0时，其即为初始图像I₀；

步骤(E)，通过Matlab软件进行仿真，验证步骤(D)建立的双向增强的扩散模型有效，峰值信噪比高，去噪性能强。

本发明的有益效果是：本发明的基于双向增强扩散滤波的图像去噪方法，简化了扩散方程，建立双向扩散系数，使模型在扩散过程中能够实现平滑与锐化的双向过程，为进一步加强平滑和锐化强度，对图像做增强处理，并用小波变换增强图像的整体轮廓，弱化图像的纹理细节，然后，对梯度阈值进行了自适应设计和改进，使其根据图像的最大灰度值和迭代次数自动控制梯度阈值，进一步保留图像边缘和细节特征，最后，对提出的模型进行仿真，用MATLAB软件对方法进行仿真验证，能够兼顾图像噪声的去除和边缘、纹理等细节信息的保护，峰值信噪比有了大幅提高，去噪性能较经典模型更具优越性，具有良好的应用前景。

附图说明

图1是本发明的基于双向增强扩散滤波的图像去噪方法的流程图。

图2是本发明的图像中某点局部坐标的示意图。

具体实施方式

下面将结合说明书附图，对本发明作进一步的说明。

本发明的基于双向增强扩散滤波的图像去噪方法，简化了扩散方程，建立双向扩散系数，使模型在扩散过程中能够实现平滑与锐化的双向过程，为进一步加强平滑和锐化强度，对图像做增强处理，并用小波变换增强图像的整体轮廓，弱化图像的纹理细节，然后，对梯度阈值进行了自适应设计和改进，使其根据图像的最大灰度值和迭代次数自动控制梯度阈值，进一步保留图像边缘和细节特征，最后，对提出的模型进行仿真，用MATLAB软件对方法进行仿真验证，能够兼顾图像噪声的去除和边缘、纹理等细节信息的保护，峰值信噪比有了大幅提高，去噪性能较经典模型更具优越性，如图1所示，包括以下步骤，

步骤(A)，建立双向扩散系数，根据PM模型的最小能量泛函，得到的扩散方程，如公式(1)所示，

$\left\{\begin{array}{c}\frac{\∂I}{\∂t}=div(\frac{f^{\′}\left(\right|\▿I\left|\right)}{|\▿I|}\·\▿I)\\ I(t=0)=I_0\end{array}\right.---\left(1\right)$

其中，I为图像，div、▽分别为图像I的散度算子和梯度算子，I₀为初始图像，f(|▽I|)为非增扩散函数、f'(|▽I|)为函数f(|▽I|)的一阶导数、I(t＝0)为图像在t＝0时的值、t为扩散时间尺度，扩散系数为关于梯度的函数，为经典PM模型的扩散系数，PM模型的扩散方程为，

$\frac{\∂I}{\∂t}=div\left(g\right(|\▿I|)\·\▿I)$

其中，div、▽分别为散度算子和梯度算子，I₀代表初始图像，它是通过原图像与高斯核卷积获得，即I(x,y,t)＝I₀*G(x,y,t)，扩散系数g(|▽I|)满足是关于梯度的函数，即在梯度低的区域平滑效果很强，而在接近边缘梯度较高的区域，PM方法尽可能地减少平滑，实现了保边缘，Perona和Mailik提出的g(|▽I|)有以下两种形式：

$g_1\left(\right|\▿I\left|\right)=\frac{1}{1+{\left(\right|\▿I|/k)}^2}$ 或者 $g_2\left(\right|\▿I\left|\right)=\exp \[-{\left(\frac{|\▿I|}{k}\right)}^2\];$

构建公式(1)内在的坐标形式，如图2所示，设图像中某点

o(i,j)的内在坐标系为(η,ξ)，η为图像I的梯度方向；ξ为垂直于梯度的方向， $\mathrm{\η}=\frac{(I_x,I_y)}{\sqrt{I_x^2+I_y^2}},\mathrm{\ξ}=\frac{(-I_y,I_x)}{\sqrt{I_x^2+I_y^2}}$

则， $I_{\mathrm{\ξ}\mathrm{\ξ}}=\frac{I_y^2{I}_{xx}-2I_x{I}_y{I}_{xy}+I_x^2{I}_{yy}}{I_x^2+I_y^2}{I}_{\mathrm{\η}\mathrm{\η}}=\frac{I_x^2{I}_{xx}+2I_x{I}_y{I}_{xy}+I_y^2{I}_{yy}}{I_x^2+I_y^2}$

其中，x,y为两个正交方向，I_x为I在x方向分量，I_y为I在y方向分量，I_xx，I_xy分别为I_x在x,y方向的方向导数，I_yy为I_y在y方向的方向导数，I_ηη表示I沿着边缘梯度方向的方向导数，I_ξξ表示I沿着边缘切线方向的方向导数；

扩散方程的最终表达式，如公式(2)所示，

$\frac{\∂I}{\∂t}=\frac{f^{\′}\left(\right|\▿I\left|\right)}{2|\▿I|}{z}_{\mathrm{\ξ}\mathrm{\ξ}}+\frac{f^{\′\′}\left(\right|\▿I\left|\right)}{2}{I}_{\mathrm{\η}\mathrm{\η}}---\left(2\right)$

为了使公式(1)有解，f(|▽I|)单调递增，f′(|▽I|)＞0；由公式(2)可知，扩散系数为沿图像的边缘方向，即沿着边缘方向进行正向扩散，使图像模糊；而扩散系数为沿图像的梯度方向，其正负号不确定，当f”(|▽I|)＞0时，沿梯度方向正扩散，使图像模糊，当f”(|▽I|)＜0时，沿梯度方向逆扩散，锐化图像的边缘，保护图像的边缘；

步骤(B)，突出图像的整体结构，增强边缘的锐化作用，对公式(2)做进一步化简，如公式(3)所示，

$\begin{array}{c}\frac{\∂I}{\∂t}=\frac{f^{\′}\left(\right|\▿I\left|\right)}{2|\▿I|}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+{\[|\▿I|\·\frac{f^{\′}\left(\right|\▿I\left|\right)}{2|\▿I|}\]}^{\′}{I}_{\mathrm{\η}\mathrm{\η}}\\ =g\left(\right|\▿I\left|\right)I_{\mathrm{\ξ}\mathrm{\ξ}}+{\[|\▿I|\·g\left(\right|\▿I\left|\right)\]}^{\′}{I}_{\mathrm{\η}\mathrm{\η}}\end{array}---\left(3\right)$

其中，为经典模型简化分析后的扩散系数；

令k是图像I的梯度阈值，当|▽I|→0时，g(|▽I|)→1；当|▽I|→∞时，g(|▽I|)→0，扩散系数的选取条件，将代入公式(3)，得到公式(4)，

$\begin{array}{c}\frac{\∂I}{\∂t}=\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+{\left(\frac{|\▿I|(1+k^2)}{|\▿I{|}^2+k^2}\right)}^{\′}{I}_{\mathrm{\η}\mathrm{\η}}\\ =\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+\left(\frac{(1+k^2)(k^2-|\▿I{|}^2)}{{\left(\right|\▿I{|}^2+k^2)}^2}\right)I_{\mathrm{\η}\mathrm{\η}}\\ =\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+\frac{(1+k^2)}{|\▿I{|}^2+k^2}\frac{(k^2-|\▿I{|}^2)}{|\▿I{|}^2+k^2}{I}_{\mathrm{\η}\mathrm{\η}}\end{array}---\left(4\right)$

当在图像边缘时，|▽I|＞＞k，得到公式(5)，

$\begin{array}{c}\frac{\∂I}{\∂t}\≈\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+\frac{(1+k^2)(0-|\▿I{|}^2)}{|\▿I{|}^2+0|\▿I{|}^2+k^2}{I}_{\mathrm{\η}\mathrm{\η}}\\ =\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}-\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\η}\mathrm{\η}}\\ =g\left(\right|\▿I\left|\right)I_{\mathrm{\ξ}\mathrm{\ξ}}-g\left(\right|\▿I\left|\right)I_{\mathrm{\η}\mathrm{\η}}\end{array}---\left(5\right)$

由公式(5)可知，图像沿梯度方向逆扩散，锐化图像的边缘，能够改善图像边缘模糊的现象，保护图像的边缘；

当在图像平坦区域时，|▽I|＜＜k，得到公式(6)，

$\begin{array}{c}\frac{\∂I}{\∂t}\≈\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+\frac{(1+k^2)}{0+k^2}\frac{(k^2-0)}{|\▿I{|}^2+k^2}{I}_{\mathrm{\η}\mathrm{\η}}\\ =\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\η}\mathrm{\η}}\\ =g\left(\right|\▿I\left|\right)I_{\mathrm{\ξ}\mathrm{\ξ}}+g\left(\right|\▿I\left|\right)I_{\mathrm{\η}\mathrm{\η}}\end{array}---\left(6\right)$

由公式(6)可知，图像的平坦区域，沿着图像的边缘方向和梯度方向同时平滑，具有去噪能力；

随着梯度的增大，当梯度达到某一值时，扩散系数出现了负值，符合以上的分析过程，本发明实质上是通过图像梯度与阈值的大小进行比较实现的，通常在扩散过程中，会预先设定一个梯度阈值k，k一般是个常量，即不论像素点所在区域与相邻区域的灰度值存在怎样的关系，都采用这个预先设定好的常量梯度阈值，这样不利于某些区域图像的细节信息和图像边缘的保留。随着扩散时间和扩散次数的增加，像素点所在区域的梯度阈值不断减小，为了正确区分图像边缘与噪声，梯度阈值k也应当随着扩散时间的增加而减小，这样就能够更准确地缩小边缘判断的范围，保留更多的图像边缘和细节；

步骤(C)，根据图像的最大灰度值和迭代次数自动控制梯度阈值k，进一步保留图像边缘和细节特征，梯度阈值k为一个随扩散时间和扩散次数变化的函数，如公式(7)所示，

k(t)＝e^-αt(7)

其中，t为扩散时间，扩散次数n代表扩散时间，则t＝n，max{W,H}为图像的最大灰度值，随着扩散次数增加或者扩散时间增加，梯度阈值k不断减小，准确的缩小边缘判断的范围，保留图像边缘和细节；

在滤波过程中为了不让图像失真，给图像添加一项保真项，保真项，如公式(8)所示，

$D={\mathrm{\λ}}^{*}(1-\frac{1}{|\▿I{|}^2+1})(I-I_0)---\left(8\right)$

其中，_λ*是调整参数，在图像的边缘，此时|▽I|→∞，滤波结果会接近原始图像，加强了保边缘的效果；在图像的平坦区域，此时|▽I|→0，则扩散程度达到最大，去除噪声；

步骤(D)，根据公式(4)和公式(8)，建立双向增强的扩散模型，如公式(9)所示，

$\left\{\begin{array}{c}\frac{\∂I}{\∂t}=\frac{1+k{\left(t\right)}^2}{|\▿\hat{I}{|}^2+k{\left(t\right)}^2}{\hat{I}}_{\mathrm{\ξ}\mathrm{\ξ}}+\frac{(1+k{\left(t\right)}^2)(k{\left(t\right)}^2-|\▿\hat{I}{|}^2)}{|\▿\hat{I}{|}^2+k{\left(t\right)}^2|\▿\hat{I}{|}^2+k{\left(t\right)}^2}{\hat{I}}_{\mathrm{\η}\mathrm{\η}}\\ -{\mathrm{\λ}}^{*}(1-\frac{1}{|\▿I{|}^2+1})(I-I_0)\\ I(x,y,0)=I_0\end{array}\right.---\left(9\right)$

其中，分别表示为I是经过小波增强处理后沿着边缘切线方向的方向导数和沿着边缘梯度方向的方向导数。I(x,y,0)＝I₀表示I(x,y,t)在t＝0时，其即为初始图像I₀；

步骤(E)，通过Matlab软件进行仿真，验证步骤(D)建立的双向增强的扩散模型有效，峰值信噪比高，去噪性能强，本发明用小波sym4对Nuist图像进行处理，设定图像的高频系数为350，若大于该高频系数，则使高频系数增大为原来的2倍，否则缩小为原来的一半，以此来突出图像的轮廓与整体结构，弱化细节，双向增强的扩散模型实现了平滑和锐化的双向过程，比较与现有的算法的均方差(MSE)和峰值信噪比(PSNR)来评价算法的有效性，其中，

$MSE=\frac{1}{W\×H}\underset{i=1}{\overset{W}{\Σ}}\underset{j=1}{\overset{H}{\Σ}}{\[I(i,j)-I_0(i,j)\]}^2$

$PSNR=10{\log }_{10}\left(\frac{{255}^2}{MSE}\right)$

图像的分辨率为W×H，I₀和I分别上表示原始图像和去噪后的图像，均方差(MSE)越小越好，峰值信噪比(PSNR)越大越好，具体结果如下表1，

表1各模型计算得到的MSE和PSNR比较

由上表可知，PM模型去噪后的图像噪声基本没有被去除，且图像中产生了“块状”效应；TV模型去噪后的图像中出现了一些孔洞，这是由于该模型不符合形态学原则，破坏了图像，但该模型有一定的去噪效果，图像的边缘保持较好；NLM算法取得了较好的去噪效果，但是对原图像结构信息保护不够，NLM处理后的图像丢失了很多边缘纹理等细节信息；BM3D去算法不仅有较高的信噪比，视觉效果也较好，但不能有效的保护图像细节信息。IterativePPB算法对图像的处理效果较好，但是该方法耗时过长，不利于实际应用，而Non‐iterativePPB算法虽然解决了IterativePPB算法耗时过长的问题，但是，破坏了图像的边缘纹理等细节信息；本发明飞双向增强扩散模型的可视性最好，主要由于该模型建立平滑与锐化的双向扩散过程，更好地考虑图像的局部特征，兼顾了去噪与边缘纹理等细节信息的保护。因此，得到双向增强扩散模型效果最好，与滤波结果的可视性相一致。

以上显示和描述了本发明的基本原理、主要特征及优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (1)

1.基于双向增强扩散滤波的图像去噪方法，其特征在于：包括以下步骤，

步骤(A)，建立双向扩散系数，根据PM模型的最小能量泛函，得到的扩散方程，如公式(1)所示，

$\left\{\begin{array}{c}\frac{\∂I}{\∂t}=div(\frac{f^{\′}\left(\right|\▿I\left|\right)}{|\▿I|}\·\▿I)\\ I(t=0)=I_0\end{array}\right.---\left(1\right)$

其中，I为图像，div、▽分别为图像I的散度算子和梯度算子，I₀为初始图像，f(|▽I|)为非增扩散函数、f'(|▽I|)为函数f(|▽I|)的一阶导数、I(t＝0)为图像在t＝0时的值，t为扩散时间尺度，扩散系数为关于梯度的函数，为经典PM模型下的扩散系数；

构建公式(1)内在的坐标形式，设图像中某点o(i,j)的内在坐标系为(η,ξ)，η为图像I的梯度方向；ξ为垂直于梯度的方向， $\mathrm{\η}=\frac{(I_x,I_y)}{\sqrt{I_x^2+I_y^2}},\mathrm{\ξ}=\frac{(-I_y,I_x)}{\sqrt{I_x^2+I_y^2}}$

则， $I_{\mathrm{\ξ}\mathrm{\ξ}}=\frac{I_y^2{I}_{xx}-2I_x{I}_y{I}_{xy}+I_x^2{I}_{yy}}{I_x^2+I_y^2},I_{\mathrm{\η}\mathrm{\η}}=\frac{I_x^2{I}_{xx}+2I_x{I}_y{I}_{xy}+I_y^2{I}_{yy}}{I_x^2+I_y^2},$

其中，x,y为两个正交方向，I_x为图像I在x方向分量，I_y为图像I在y方向分量，I_xx，I_xy分别为I_x在x,y方向的方向导数，I_yy为I_y在y方向的方向导数，I_ηη表示图像I沿着边缘梯度方向的方向导数，I_ξξ表示图像I沿着边缘切线方向的方向导数；

扩散方程的最终表达式，如公式(2)所示，

$\frac{\∂I}{\∂t}=\frac{f^{\′}\left(\right|\▿I\left|\right)}{2|\▿I|}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+\frac{f^{\′\′}\left(\right|\▿I\left|\right)}{2}{I}_{\mathrm{\η}\mathrm{\η}}---\left(2\right)$

为了使公式(1)有解，f(|▽I|)单调递增，f′(|▽I|)＞0；由公式(2)可知，扩散系数为沿图像的边缘方向，即沿着边缘方向进行正向扩散，使图像模糊；而扩散系数为沿图像的梯度方向，其正负号不确定，当f”(|▽I|)＞0时，沿梯度方向正扩散，使图像模糊，当f”(|▽I|)＜0时，沿梯度方向逆扩散，锐化图像的边缘，保护图像的边缘；

步骤(B)，突出图像的整体结构，增强边缘的锐化作用，对公式(2)做进一步化简，如公式(3)所示，

$\begin{array}{c}\frac{\∂I}{\∂t}=\frac{f^{\′}\left(\right|\▿I\left|\right)}{2|\▿I|}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+{\[|\▿I|\·\frac{f^{\′}\left(\right|\▿I\left|\right)}{2|\▿I|}\]}^{\′}{I}_{\mathrm{\η}\mathrm{\η}}\\ =g\left(\right|\▿I\left|\right)I_{\mathrm{\ξ}\mathrm{\ξ}}+{\[|\▿I|\·g\left(\right|\▿I\left|\right)\]}^{\′}{I}_{\mathrm{\η}\mathrm{\η}}\end{array}---\left(3\right)$

其中，为经典PM模型简化分析后的扩散系数；

令k是图像I的梯度阈值，当|▽I|→0时，g(|▽I|)→1；当|▽I|→∞时，g(|▽I|)→0，扩散系数的选取条件，将代入公式(3)，得到公式(4)，

$\begin{array}{c}\frac{\∂I}{\∂t}=\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+{\left(\frac{|\▿I|(1+k^2)}{|\▿I{|}^2+k^2}\right)}^{\′}{I}_{\mathrm{\η}\mathrm{\η}}\\ =\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+\left(\frac{(1+k^2)(k^2-|\▿I{|}^2)}{{\left(\right|\▿I{|}^2+k^2)}^2}\right)I_{\mathrm{\η}\mathrm{\η}}\\ =\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+\frac{1+k^2}{|\▿I{|}^2+k^2}\frac{(k^2-|\▿I{|}^2)}{|\▿I{|}^2+k^2}{I}_{\mathrm{\η}\mathrm{\η}}\end{array}---\left(4\right)$

当在图像边缘时，|▽I|＞＞k，得到公式(5)，

$\begin{array}{c}\frac{\∂I}{\∂t}\≈\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+\frac{(1+k^2)}{|\▿I{|}^2+0}\frac{(0-|\▿I{|}^2)}{|\▿I{|}^2+k^2}{I}_{\mathrm{\η}\mathrm{\η}}\\ =\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}-\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\η}\mathrm{\η}}\\ =g\left(\right|\▿I\left|\right)I_{\mathrm{\ξ}\mathrm{\ξ}}-g\left(\right|\▿I\left|\right)I_{\mathrm{\η}\mathrm{\η}}\end{array}---\left(5\right)$

由公式(5)可知，图像沿梯度方向逆扩散，锐化图像的边缘，能够改善图像边缘模糊的现象，保护图像的边缘；

当在图像平坦区域时，|▽I|＜＜k，得到公式(6)，

$\begin{array}{c}\frac{\∂I}{\∂t}\≈\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}+\frac{(1+k^2)}{0+k^2}\frac{(k^2-0)}{\▿I{|}^2+k^2}{I}_{\mathrm{\η}\mathrm{\η}}\\ =\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\ξ}\mathrm{\ξ}}-\frac{1+k^2}{|\▿I{|}^2+k^2}{I}_{\mathrm{\η}\mathrm{\η}}\\ =g\left(\right|\▿I\left|\right)I_{\mathrm{\ξ}\mathrm{\ξ}}-g\left(\right|\▿I\left|\right)I_{\mathrm{\η}\mathrm{\η}}\end{array}---\left(6\right)$

由公式(6)可知，图像的平坦区域，沿着图像的边缘方向和梯度方向同时平滑，具有去噪能力；

步骤(C)，根据图像的最大灰度值和迭代次数自动控制梯度阈值k，进一步保留图像边缘和细节特征，梯度阈值k为一个随扩散时间和扩散次数变化的函数，如公式(7)所示，

k(t)＝e^-αt(7)

其中，t为扩散时间，扩散次数n代表扩散时间，则t＝n， $\mathrm{\α}=\frac{1}{max\{W,H\}},$ max{W,H}为图像的最大灰度值，随着扩散次数增加或者扩散时间增加，梯度阈值k不断减小，准确的缩小边缘判断的范围，保留图像边缘和细节；

在滤波过程中为了不让图像失真，给图像添加一项保真项，保真项，如公式(8)所示，

$D={\mathrm{\λ}}^{*}(1-\frac{1}{|\▿I{|}^2+1})(I-I_0)---\left(8\right)$

其中，λ^*是调整参数，在图像的边缘，此时|▽I|→∞，滤波结果会接近原始图像，加强了保边缘的效果；在图像的平坦区域，此时|▽I|→0， $\frac{1}{|\▿I{|}^2+1}\→1,(1-\frac{1}{|\▿I{|}^2+1})\≈0,$ 则扩散程度达到最大，去除噪声；

步骤(D)，根据公式(4)和公式(8)，建立双向增强的扩散模型，如公式(9)所示，

$\left\{\begin{array}{c}\frac{\∂I}{\∂t}=\frac{1+k{\left(t\right)}^2}{|\▿\hat{I}{|}^2+k{\left(t\right)}^2}{\hat{I}}_{\mathrm{\ξ}\mathrm{\ξ}}+\frac{(1+k{\left(t\right)}^2)}{|\▿\hat{I}{|}^2+k{\left(t\right)}^2}\frac{(k{\left(t\right)}^2-|\▿\hat{I}{|}^2)}{|\▿\hat{I}{|}^2+k{\left(t\right)}^2}{\hat{I}}_{\mathrm{\η}\mathrm{\η}}\\ -{\mathrm{\λ}}^{*}(1-\frac{1}{|\▿I{|}^2+1})(I-I_0)\\ I(x,y,0)=I_0\end{array}\right.---\left(9\right)$

其中，分别表示为图像I是经过小波增强处理后沿着边缘切线方向的方向导数和沿着边缘梯度方向的方向导数，I(x,y,0)＝I₀表示I(x,y,t)在t＝0时，其即为初始图像I₀；

步骤(E)，通过Matlab软件进行仿真，验证步骤(D)建立的双向增强的扩散模型有效，峰值信噪比高，去噪性能强。