CN104331869A - 梯度与曲率相结合的图像平滑方法 - Google Patents

Abstract

本发明首先提供了一种图像曲率平滑方法，将图像水平集曲率作为一个二阶微分量，用来描述图像形态学特征，改善了图像处理效果。由于当图像受到噪声污染时，曲率会发生显著变化，本发明进一步将图像的水平集曲率作为一个检测因子代入到PM方法中，提出了一种梯度与曲率相结合的图像平滑方法，该方法不但可以保持图像的特征，而且峰值信噪比较以往的各向异性方法大幅度提高，复杂度低，时效性高，图像清晰度增加，相较以往图像平滑算法更具有效性和准确性，图像处理效果更好。

Description

梯度与曲率相结合的图像平滑方法

技术领域

本发明属于图像处理技术领域，具体涉及一种将梯度与曲率相结合的图像平滑方法。

背景技术

在图像处理与计算机领域中，图像去噪是最基本的问题之一。图像平滑就是通过减少任何可能存在的退化而增强图像的技术。近年来，偏微分方程(PDE)方法在图像分析处理和计算机视觉领域的应用非常广泛。偏微分方程能反映未知变量关于时间的导数和关于空间变量的导数之间制约关系。基于偏微分方程的图像处理方法将离散图像表示成连续的数学模型，利用偏微分方程完善的数值分析理论对图像进行处理。相比于传统方法，该方法具有更强的局部自适应能力和更高的灵活性，因此在图像的去噪、分割、边缘检测、增强等方面都有重要应用。

20世纪90年代初期，Perona和Malik首先提出了各向异性扩(Anisotropic Diffusion)的理论(即PM方法)，将扩散系数由常数改为关于梯度模值的单调非增函数，使滤波模型在去除噪声的同时还能对边缘起到保护作用。

PM方法的扩散方程为：

$\left\{\begin{array}{c}\frac{\∂I}{\∂t}\text{=div}(g\left(\right|\▿I\left|\right)\·\▿I)\\ I(x,y,0)=I_0\end{array}\right.---\left(1\right)$

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

$g_1\left(\right|\▿I\left|\right)=\frac{1}{1+{\left(\right|\▿I|/k)}^2}---\left(2\right)$

$g_2\left(\right|\▿I\left|\right)=\exp [-{\left(\frac{|\▿I|}{k}\right)}^2]---\left(3\right)$

自PM方法提出后，各向异性扩散技术有了长足的发展。随着对该技术研究的不断深入，许多实验结果表明，PM方法中存在缺陷，这些缺点使得处理后的图像降噪不稳定，有明显的“阶梯”效应，并且保边缘性不是很好。

由于方法的稳定性直接影响了滤波的效果，许多学者对PM方法中扩散系数进行了改进,试图建立起更有效的保边缘平滑滤波器,并取得了一些进展。然而大部分改进方法都没有从根本上解决PM方法的缺陷。综上，现有的方法中主要存在如下问题：当图像中存在噪声时，由于噪声点的梯度比较大，扩散系数较小，噪声不能完全去除,且带有噪声的边界也会被破坏掉；方法本身在数学上是病态的，解的存在性和唯一性不能得到保证，也就导致了图像处理效果并不稳定。

发明内容

我们通过对现有方法的研究发现，仅用一阶微分量(梯度)来表征图像局部特征是不够的，二阶微分量中含有更丰富的信息。因此，为了解决现有技术中的缺陷，本发明将图像水平集曲率作为一个二阶微分量，用来描述图像形态学特征，提出了一种图像曲率平滑方法。

为了达到上述目的，本发明提供如下技术方案：

图像曲率平滑方法，包括如下步骤：

步骤A，将噪声图像进行滤波处理；

步骤B，由PM模型输出二值图像；

步骤C，将图像水平集的曲率作为一个二阶微分量描述图像形态学特征，建立曲率平滑模型：

$\left\{\begin{array}{c}\frac{\∂I}{\∂t}=\mathrm{div}(\frac{1}{1+{\left(\right|\mathrm{\κ}|/l)}^2}\▿I)\\ I(t=0)=I_0\end{array}\right.$

其中为水平集曲率，|κ|为曲率模值，l为阈值，I₀代表初始图像，I是通过原图像与高斯核卷积获得，即I(x,y,t)＝I₀*G(x,y,t)；

步骤D，对步骤C中提供的图像平滑模型用半隐式加性算子分裂数值方法进行计算。

进一步的，所述步骤D具体包括如下过程：

将图像I尺度化分解在[0,1]区间，并进行简化,完成Iⁿ后进行下列步骤：

步骤D-1，令

步骤D-2，计算f_σ＝f*G_σ，|▽f_σ|_ij，|κ|_ij，

步骤D-3，当i＝1,…,M时，计算的三个对角线上的元素： $({\mathrm{\β}}_k^{\left(i\right)},k=1,...,N-1),({\mathrm{\γ}}_k^{\left(i\right)},k=2,...,N),$ 并采用追赶法求解 $(I-2\mathrm{\τ}{A}_{x,i}^n)I_{1i}^{n+1}=I_{1i}^n,$ 得到

步骤D-4，当j＝1,…,N时，同样计算的三个对角线上的元素，并采用追赶法求解 $(I-2\mathrm{\τ}{A}_{y,j}^n)I_{2j}^{n+1}=I_{2j}^n,$ 得到

步骤D-5，计算 $I^{n+1}=\frac{1}{2}(I_1^{n+1}+I_2^{n+1}).$

进一步的，所述步骤B中PM模型如下：

$\frac{\∂I}{\∂t}\text{=div}(\frac{1}{1+\left(\right|\▿I|/k)}\▿I).$

当图像受到噪声污染时，曲率会发生显著变化，因此我们将图像的水平集曲率作为一个检测因子代入到PM方法中，提出了一种梯度与曲率相结合的图像平滑方法。

梯度与曲率相结合的图像平滑方法，包括如下步骤：

步骤A，将噪声图像进行滤波处理；

步骤B，由PM模型输出二值图像；

步骤C，将图像的水平集曲率作为一个检测因子代入到PM方法中，建立一个梯度与曲率相结合的图像平滑模型：

$\frac{\∂I}{\∂t}\text{=div}(\frac{1}{1+(\left(\right|\▿I|+|\mathrm{\κ}\left|\right)/k)}\▿I)$

其中，水平集曲率为 $\mathrm{\κ}=\frac{I_x^2{I}_{\mathrm{yy}}-2I_x{I}_y{I}_{\mathrm{xy}}+I_y^2{I}_{\mathrm{xx}}}{{(I_x^2+I_y^2)}^{\frac{3}{2}}},$

▽I为水平集，div、▽分别为散度算子和梯度算子，k为阈值；

步骤D，对步骤E中提供的图像平滑模型用半隐式加性算子分裂数值方法进行计算。

进一步的，所述步骤D具体包括如下过程：

将图像I尺度化分解在[0,1]区间，并进行简化,完成Iⁿ后进行下列步骤：

步骤D-1，令

步骤D-2，计算f_σ＝f*G_σ，|▽f_σ|_ij，|κ|_ij，

步骤D-3，当i＝1,…,M时，计算的三个对角线上的元素： $({\mathrm{\β}}_k^{\left(i\right)},k=1,...,N-1),({\mathrm{\γ}}_k^{\left(i\right)},k=2,...,N),$ 并采用追赶法求解 $(I-2\mathrm{\τ}{A}_{x,i}^n)I_{1i}^{n+1}=I_{1i}^n,$ 得到

步骤D-4，当j＝1,…,N时，同样计算的三个对角线上的元素，并采用追赶法求解 $(I-2\mathrm{\τ}{A}_{y,j}^n)I_{2j}^{n+1}=I_{2j}^n,$ 得到

步骤D-5，计算 $I^{n+1}=\frac{1}{2}(I_1^{n+1}+I_2^{n+1}).$

进一步的，所述步骤B中PM模型如下：

$\frac{\∂I}{\∂t}\text{=div}(\frac{1}{1+\left(\right|\▿I|/k)}\▿I).$

与现有技术相比，本发明具有如下优点和有益效果：

本发明首先将图像水平集曲率作为一个二阶微分量描述图像形态学特征，改善了图像处理效果，又进一步将图像的水平集曲率作为一个检测因子代入到PM方法中去，不但可以保持图像的特征，而且峰值信噪比较以往的各向异性方法大幅度提高，复杂度低，时效性高，图像清晰度增加，相较以往图像平滑算法更具有效性和准确性，图像处理效果更好。

附图说明

图1为本发明步骤流程图；

图2为受噪声污染的图像；

图3为PM方法处理后的仿真图；

图4为基于曲率方法处理后的仿真图；

图5为本发明提供的梯度与曲率相结合的图像平滑方法处理后的仿真图；

图6为PM方法处理后的局部仿真图；

图7为梯度与曲率相结合的图像平滑方法处理后的局部仿真图

具体实施方式

以下将结合具体实施例对本发明提供的技术方案进行详细说明，应理解下述具体实施方式仅用于说明本发明而不用于限制本发明的范围。

由于二阶微分量中含有更丰富的信息，本发明将形态学中的水平集曲率作为表征图像的局部特征，提供了一种图像曲率平滑方法，观察平滑后的图像处理结果后，发现当图像受到噪声污染时，曲率会发生显著变化，因此本发明进一步将图像的水平集曲率作为一个检测因子代入到PM方法中，建立一个梯度与曲率相结合的图像平滑方法，最后用MATLAB软件对方法进行仿真，得到去噪图像和数值结果用来评价方法。本发明整体流程图如图1所示，其中图像曲率平滑方法具体包括如下步骤：

步骤A，将噪声图像进行Gauss滤波，去除较大噪声；

步骤B，首先采用公式(2)表达的g(x)，由经典的PM方法输出二值图像，其中div、▽分别为散度算子和梯度算子，k为阈值，下同。观察该二值图像，得出该传统的图像去噪复原方法及其后来的改进方法均以|▽I|作为图像的边缘检测算子，而这类方法存在一个不足：它不完全符合图像处理的形态学原则，在图像扩散过程中，它的变化不仅取决于水平集(由▽I表征)，同时还取决于灰度值(I)，故上述方法在边缘和图像灰度渐变区，及去除平坦区的孤立噪声时就显得无能为力，会产生阶梯效应。

步骤C，由上述步骤分析可知，仅用一阶微分量(梯度)来表征图像局部特征是不够的，二阶微分量中含有更丰富的信息。由上面曲率的定义及其表达式可知，曲率是几何体不平坦程度的一种度量，能够区分边缘等细节信息，本发明将图像水平集的曲率作为一个二阶微分量，用来描述图像形态学特征，由此建立曲率平滑模型：

$\left\{\begin{array}{c}\frac{\∂I}{\∂t}=\mathrm{div}(\frac{1}{1+{\left(\right|\mathrm{\κ}|/l)}^2}\▿I)\\ I(t=0)=I_0\end{array}\right.$

其中为水平集曲率，|κ|为曲率模值，l为阈值，I₀代表经过步骤B处理后的二值图像，I是通过原图像I₀与高斯核卷积获得，即I(x,y,t)＝I₀*G(x,y,t)。

前述的水平集曲率通过如下过程获得：

设平面曲线的切矢量为T，且为单位矢量，记为

T＝C_s,|T|＝1          (4)

由于C_s为单位矢量，它与自身的内积为1，即<C_s,C_s>＝||C_s||²＝1，该式两边求导可得<C_s,C_ss>＝0，可见矢量C_ss与单位切矢量C_s正交。现定义与T构成右手坐标系的单位矢量为法矢量N。可见矢量C_ss与N共线，故可表达为

C_ss＝κN          (5)

式中，比例系数κ称为曲率。

由式(4)、式(5)可得

T_s＝κN          (6)

式(6)中，曲率κ(s)的几何意义是切矢量T随弧长的变化率。现假定曲线上某一点s的切矢量和法矢量分别为T(s)＝(cosθ,sinθ)，N(s)＝(-sinθ,cosθ)，式中θ表示T与x轴的夹角。当沿曲线移动到s+Δs点时， $\frac{T(s+\mathrm{\&Delta;s})-T\left(s\right)}{\mathrm{\&Delta;s}}\&cong;(-\sin \mathrm{\&theta;}\frac{\mathrm{\&Delta;\&theta;}}{\mathrm{\&Delta;s}},\cos \mathrm{\&theta;}\frac{\mathrm{\&Delta;\&theta;}}{\mathrm{\&Delta;s}})=\frac{\mathrm{\&Delta;\&theta;}}{\mathrm{\&Delta;s}}N\left(s\right)$ 取Δs→0的极限，并与式(6)比较得

$\mathrm{\&kappa;}=\frac{\mathrm{d\&theta;}}{\mathrm{ds}}---\left(7\right)$

这说明曲率是切矢量的旋转角速度，同时也是法矢量的旋转角速度。

又因为单位法矢量N(s)＝(-sinθ,cosθ)，所以

$\frac{\&PartialD;n1}{\&PartialD;x}=\frac{\&PartialD;n1}{\&PartialD;\mathrm{\&theta;}}\frac{\&PartialD;\mathrm{\&theta;}}{\&PartialD;s}\frac{\&PartialD;s}{\&PartialD;x}=-\cos \mathrm{\&theta;}\&CenterDot;\mathrm{\&kappa;}\&CenterDot;\frac{\&PartialD;s}{\&PartialD;x},\frac{\&PartialD;n2}{\&PartialD;x}=-\sin \mathrm{\&theta;}\&CenterDot;\mathrm{\&kappa;}\&CenterDot;\frac{\&PartialD;s}{\&PartialD;y}---\left(8\right)$

但ds＝dxcosθ+dysinθ，将它与相比较，可得将它们带入式(6)，可得 $\frac{\&PartialD;n1}{\&PartialD;x}=-\mathrm{\&kappa;}{\cos }^2\mathrm{\&theta;},\frac{\&PartialD;n2}{\&PartialD;y}=-\mathrm{\&kappa;}{\sin }^2\mathrm{\&theta;},$ 从而可得

$\mathrm{\&kappa;}=-(\frac{\&PartialD;n1}{\&PartialD;x}+\frac{\&PartialD;n2}{\&PartialD;y})=-\mathrm{div}\left(N\right)---\left(9\right)$

对于平面封闭曲线，

C＝{(x,y),I(x,y)＝0}        (10)

式中，I(x,y)为某一个二维函数。也就是说曲线C满足I(x,y)＝c的点集，称为函数I(x,y)的一个水平(线)集，称I(x,y)是曲线C的嵌入函数。如果在水平集的某一点p沿水平集的切线方向对I(x,y)求方向导数，由于I(x,y)沿水平集保持不变，则式中的θ表示切矢量t与x轴的夹角。可见，I(x,y)的梯度矢量为

$\&dtri;I=(\frac{\&PartialD;I}{\&PartialD;x},\frac{\&PartialD;I}{\&PartialD;y})---\left(11\right)$

它与水平集的切矢量T(s)＝(cosθ,sinθ)相垂直，即与水平集的法矢量平行。另一方面，根据式(9)，梯度矢量总是指向I值增大的方向，可见水平集的单位法矢量可表示为

$N=\&PlusMinus;\frac{\&dtri;I}{|\&dtri;I|}---\left(12\right)$

一般约定式(12)取负号，把式(12)代入式(9)中，便可求得嵌入函数I(x,y)水平集曲率为

$\begin{array}{c}\mathrm{\&kappa;}=\mathrm{div}\left(\frac{\&dtri;I}{|\&dtri;I|}\right)=\mathrm{grad}\left(\frac{1}{|\&dtri;I|}\mathrm{div}(\&dtri;I)\right)\\ =\left(\begin{array}{c}-\frac{1}{2}{(I_x^2+I_y^2)}^{-\frac{3}{2}}(2I_x{I}_{\mathrm{xx}}+2I_y{I}_{\mathrm{yy}})\\ -\frac{1}{2}{(I_x^2+I_y^2)}^{-\frac{3}{2}}(2I_x{I}_{\mathrm{xy}}+2I_y{I}_{\mathrm{yy}})\end{array}\right)\left(\begin{array}{c}{I}_x\\ I_y\end{array}\right)+\frac{I_{\mathrm{xx}}+I_{\mathrm{yy}}}{{(I_x^2+I_y^2)}^{\frac{1}{2}}}\\ =-\frac{I_x^2{I}_{\mathrm{xx}}+2I_x{I}_y{I}_{\mathrm{xy}}+I_y^2{I}_{\mathrm{yy}}}{{(I_x^2+I_y^2)}^{\frac{3}{2}}}+\frac{(I_{\mathrm{xx}}+I_{\mathrm{yy}})(I_x^2+I_y^2)}{{(I_x^2+I_y^2)}^{\frac{3}{2}}}\\ =\frac{I_x^2{I}_{\mathrm{yy}}-2I_x{I}_y{I}_{\mathrm{xy}}+I_y^2{I}_{\mathrm{xx}}}{{(I_x^2+I_y^2)}^{\frac{3}{2}}}\end{array}---\left(13\right)$

步骤D，根据步骤C得到的平滑模型，用半隐式加性算子分裂(AOS)数值方法进行数值计算，运算的具体过程如下：

将图像I尺度化分解在[0,1]区间，并进行简化，其简化过程如下：

当用一维矩阵向量表示法表示时，其迭代方案如下：

$I^{n+1}={[1-\mathrm{\&tau;A}\left(I^n\right)]}^{-1}{I}^n---\left(14\right)$

其中，τ是时间步长，A(Iⁿ)＝[a_ij(Iⁿ)]，且

$a_{\mathrm{ij}}\left(I^n\right)=\left\{\begin{array}{cc}\frac{{\mathrm{\&gamma;}}_i^n+{\mathrm{\&gamma;}}_j^n}{2h^2}& j\&Element;N_i\\ -\underset{k\&Element;N_i}{\mathrm{\&Sigma;}}\frac{{\mathrm{\&gamma;}}_i^n+{\mathrm{\&gamma;}}_k^n}{2h^2}& j=i\\ 0& \mathrm{otherwise}\end{array}\right.---\left(15\right)$

式中，γ_i＝a_ig_i，h是离散化步长。以此类推，当用N维矩阵向量表示法表示时，其迭代方案如下：

$I^{n+1}={[1-\mathrm{\&tau;}\underset{l=1}{\overset{n}{\mathrm{\&Sigma;}}}{A}_l\left(I^n\right)]}^{-1}{I}^n---\left(16\right)$

式中，矩阵A_l＝(a_ijl)_ij。

当完成Iⁿ后，通过下述过程进行迭代：

步骤D-1，令

步骤D-2，计算f_σ＝f*G_σ，|▽f_σ|_ij，|κ|_ij，

步骤D-3，当i＝1,…,M时，计算的三个对角线上的元素： $({\mathrm{\&beta;}}_k^{\left(i\right)},k=1,...,N-1),({\mathrm{\&gamma;}}_k^{\left(i\right)},k=2,...,N),$ 并采用追赶法求解 $(I-2\mathrm{\&tau;}{A}_{x,i}^n)I_{1i}^{n+1}=I_{1i}^n,$ 得到

步骤D-4，当j＝1,…,N时，同样计算的三个对角线上的元素，并采用追赶法求解 $(I-2\mathrm{\&tau;}{A}_{y,j}^n)I_{2j}^{n+1}=I_{2j}^n,$ 得到

步骤D-5，计算 $I^{n+1}=\frac{1}{2}(I_1^{n+1}+I_2^{n+1}).$

通过上述步骤C-1～C-5便完成了一次迭代，经过多次迭代操作便可得到一幅很清晰的图像。

由于曲率是几何体不平坦程度的一种度量，对于一幅图像来说，其水平集应该是光滑的，而当图像受到噪声污染时，其曲率会发生显著变化，因此我们还可以对本发明方法进一步改进，提出一种梯度与曲率相结合的图像平滑方法，首先同样要进行图像去噪步骤(前述步骤A)和由PM方法输出二值图像的步骤(前述步骤B)，随后采用梯度与曲率相结合的模型取代前述曲率平滑模型，即：

步骤C，将图像的水平集曲率作为一个检测因子代入到PM方法中，建立一个梯度与曲率相结合的图像平滑模型其中div、▽分别为散度算子和梯度算子，k为阈值，下同。

步骤D，建立梯度与曲率相结合的图像平滑模型之后，采用半隐式加性算子分裂(AOS)数值方法进行数值计算，计算过程与曲率平滑模型的计算过程相同(前述步骤D)。

为了验证本发明方法的有效性，用受高斯随机噪声(σ＝20)污染的图Nuist(600×600)进行仿真实验，在Matlab环境下编程实现，将本发明提供的曲率平滑方法和梯度与曲率相结合的图像平滑方法与经典的PM方法进行比较，采用上述几种方法同时对受噪声污染的图2进行图像平滑处理，各方法处理后图像如图3～图7所示。其中，图3为PM方法处理后的仿真图，图4为基于曲率的方法处理后的仿真图，由图3、图4对比可见，用基于曲率的方法处理图像，图像平滑后的图4较图3更加清晰。图5为梯度与曲率相结合的图像平滑方法处理后的仿真图，图像平滑后较图3与图4更清晰。图6、图7分别为PM方法处理后的局部仿真图及本发明梯度与曲率相结合的图像平滑方法处理后的局部仿真图，与图6相比，梯度与曲率相结合的方法保持图像细微结构的性能更好，这样不但可以保持图像的特征，如：直线、曲线的边缘，角点，斜坡和小尺度特征等。

上述几种方法得到的数值结果如表1所示，数值对比过程是将均方差(MSE)和峰值信噪比(PSNR)作为方法性能评价标准，它们的定义为

$\mathrm{MSE}=\frac{1}{W\&times;H}\underset{i=1}{\overset{W}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{H}{\mathrm{\&Sigma;}}}{[I(i,j)-I_0(i,j)]}^2---\left(17\right)$

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

$\mathrm{Definition}=\frac{1}{W\&times;H}\underset{i=1}{\overset{W}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{H}{\mathrm{\&Sigma;}}}\sqrt{{(I(i,j)-(i-1,j))}^2+{(I(i,j)-I(i,j-1))}^2}---\left(19\right)$

图像的分辨率为W×H，I和I₀表示原始图像和去噪后的图像，均方差(MSE)越小越好，峰值信噪比(PSNR)越大越好。

表1各方法处理后数值结果比较

由上表可知，曲率平滑方法对噪声抑制的效果要优于PM模型，其峰值信噪较PM模型提高了4个db左右，而梯度与曲率相结合的平滑方法峰值信噪比跟PM方法相比提高了近44个dB，增强了尖锐的边缘，是一种更为理想的方法，降低了方法复杂度，提高了方法时效性，从而简化了图像平滑问题。

本发明方案所公开的技术手段不仅限于上述实施方式所公开的技术手段，还包括由以上技术特征任意组合所组成的技术方案。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本发明的保护范围。

Claims (6)

1.图像曲率平滑方法，其特征在于，包括如下步骤：

步骤A，将噪声图像进行滤波处理；

步骤B，由PM模型输出二值图像；

步骤C，将图像水平集的曲率作为一个二阶微分量描述图像形态学特征，建立曲率平滑模型：

$\left\{\begin{array}{c}\frac{\&PartialD;I}{\&PartialD;t}=\mathrm{div}(\frac{1}{1+{\left(\right|\mathrm{\&kappa;}|/l)}^2}\&dtri;I)\\ I(t=0)=I_0\end{array}\right.$

其中为水平集曲率，|κ|为曲率模值，l为阈值，I₀代表初始图像，I是通过原图像与高斯核卷积获得，即I(x,y,t)＝I₀*G(x,y,t)；

步骤D，对步骤C中提供的图像平滑模型用半隐式加性算子分裂数值方法进行计算。

2.根据权利要求1所述的图像曲率平滑方法，其特征在于：所述步骤D具体包括如下过程：将图像I尺度化分解在[0,1]区间，并进行简化，完成Iⁿ后进行下列步骤：

步骤D-1，令

步骤D-2，计算f_σ＝f*G_σ，|κ|_ij，

步骤D-3，当i＝1,…,M时，计算的三个对角线上的元素： $({\mathrm{\&beta;}}_k^{\left(i\right)},k=1,\&CenterDot;\&CenterDot;\&CenterDot;,N-1),({\mathrm{\&gamma;}}_k^{\left(i\right)},k=2,\&CenterDot;\&CenterDot;\&CenterDot;,N),$ 并采用追赶法求解 $(I-2\mathrm{\&tau;}{A}_{x,i}^n)I_{1i}^{n+1}=I_{1i}^n,$ 得到

步骤D-4，当j＝1,…,N时，同样计算的三个对角线上的元素，并采用追赶法求解 $(I-2\mathrm{\&tau;}{A}_{y,j}^n)I_{2j}^{n+1}=I_{2j}^n,$ 得到

步骤D-5，计算 $I^{n+1}=\frac{1}{2}(I_1^{n+1}+I_2^{n+1}).$

3.根据权利要求1或2所述的图像曲率平滑方法，其特征在于，所述步骤B中PM模型如下：

$\frac{\&PartialD;I}{\&PartialD;t}=\mathrm{div}(\frac{1}{1+\left(\right|\&dtri;I|/k)}\&dtri;I).$

4.梯度与曲率相结合的图像平滑方法，其特征在于，包括如下步骤：

步骤A，将噪声图像进行滤波处理；

步骤B，由PM模型输出二值图像；

步骤C，将图像的水平集曲率作为一个检测因子代入到PM方法中，建立一个梯度与曲率相结合的图像平滑模型：

$\frac{\&PartialD;I}{\&PartialD;t}=\mathrm{div}(\frac{1}{1+(\left(\right|\&dtri;I|+|\mathrm{\&kappa;}\left|\right)/k)}\&dtri;I).$

其中，水平集曲率为 $\mathrm{\&kappa;}=\frac{I_x^2{I}_{\mathrm{yy}}-2I_x{I}_y{I}_{\mathrm{xy}}+I_y^2{I}_{\mathrm{xx}}}{{(I_x^2+I_y^2)}^{\frac{3}{2}}},$

为水平集，div、分别为散度算子和梯度算子，k为阈值；

步骤D，对步骤E中提供的图像平滑模型用半隐式加性算子分裂数值方法进行计算。

5.根据权利要求4所述的梯度与曲率相结合的图像平滑方法，其特征在于，所述步骤D具体包括如下过程：

将图像I尺度化分解在[0,1]区间，并进行简化，完成Iⁿ后进行下列步骤：

步骤D-1，令

步骤D-2，计算f_σ＝f*G_σ，|κ|_ij，

步骤D-3，当i＝1,…,M时，计算的三个对角线上的元素： $({\mathrm{\&beta;}}_k^{\left(i\right)},k=1,\&CenterDot;\&CenterDot;\&CenterDot;,N-1),({\mathrm{\&gamma;}}_k^{\left(i\right)},k=2,\&CenterDot;\&CenterDot;\&CenterDot;,N),$ 并采用追赶法求解 $(I-2\mathrm{\&tau;}{A}_{x,i}^n)I_{1i}^{n+1}=I_{1i}^n,$ 得到

步骤D-4，当j＝1,…,N时，同样计算的三个对角线上的元素，并采用追赶法求解 $(I-2\mathrm{\&tau;}{A}_{y,j}^n)I_{2j}^{n+1}=I_{2j}^n,$ 得到

步骤D-5，计算 $I^{n+1}=\frac{1}{2}(I_1^{n+1}+I_2^{n+1}).$

6.根据权利要求1或2梯度与曲率相结合的图像平滑方法，其特征在于，所述步骤B中PM模型如下：

$\frac{\&PartialD;I}{\&PartialD;t}=\mathrm{div}(\frac{1}{1+\left(\right|\&dtri;I|/k)}\&dtri;I).$