CN104463811A - 基于能量泛函的图像平滑与锐化算法 - Google Patents

Abstract

本发明涉及一种基于能量泛函的图像平滑与锐化算法，包括以下步骤：(1)将各向异性扩散方程转化为最小化能量泛函；(2)将动态方程以内在坐标形式表示；(3)建立能量泛函(4)建立图像的梯度阈值函数k＝e^-αt；(5)引入保真项

Description

基于能量泛函的图像平滑与锐化算法

技术领域

本发明涉及图像去噪算法领域，尤其是基于能量泛函的图像去噪算法。

背景技术

图像噪声的主要来源是图像在采集过程中的随机高斯噪声和图像传播过程中的椒盐噪声。传统的去噪方法有中值滤波、同态滤波，逆滤波等，这些方法在一定程度上可以达到去除噪声的目的，但它们有一个共同的弱点，在去噪的同时，也会使图像的边缘模糊化，甚至使图像的细节纹理信息丢失。相比于传统方法，基于偏微分方程的图像处理方法具有更强的局部自适应能力和更高的灵活性，在图像的去噪、分割、边缘检测、增强等方面都有重要应用。

基于偏微分方程的图像处理方法是在图像的连续数学模型基础上，令图像遵循某一指定的偏微分方程发生变化，基于变分模型最具代表性的算法是非线性偏微分方程去噪算法(TV算法)，该算法将图像归类为变分有界函数空间，采用一次范数全变分作为其“平滑性”的度量，并沿着梯度垂直的方向进行约束，很好地保护图像的边缘，但是由于它并不完全符合图像处理的形态学原则，迭代多次后通常会产生阶梯效应。

发明内容

针对上述现有技术中的缺陷，本发明提供一种基于能量泛函的图像平滑与锐化算法，其应用于图像去噪时效性高、复杂度低、精确度高。

本发明采用以下技术方案：基于能量泛函的图像平滑与锐化算法，其特征在于，包括以下步骤：

步骤一、建立噪声图像模型I₀＝KI，为使去噪后的图像与原始图像充分接近，建立所述噪声图像模型的最小化模型其动态方程为引入Neumann边界，该动态方程等价于各向异性扩散方程 $\left\{\begin{array}{c}\frac{\∂I}{\∂t}=\mathrm{div}(g\left(\right|\▿I\left|\right)\▿I)\\ \frac{\∂I}{\∂n}=0\\ I(x,y,0)=I_0\end{array}\right.,$ 即各向异性扩散方程转化为最小化能量泛函，其中I是无噪图像，I₀是降质图像，K是降质算子，Ω是图像区域；

步骤二、将步骤一中的动态方程以内在坐标形式表示为

步骤三、建立能量泛函式(1)表示为

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

取 $g\left(\right|\▿I\left|\right)=\frac{1+k^2}{{|\▿I|}^2+k^2},$ 代入式2)得

$\frac{\∂I}{\∂t}=\frac{1+k^2}{{|\▿I|}^2+k^2}{I}_{\mathrm{\ξ\ξ}}+\frac{(1+k^2)}{{|\▿I|}^2+k^2}\frac{(k^2-{|\▿I|}^2)}{{|\▿I|}^2+k^2}{I}_{\mathrm{\η\η}}---3)$

在图像边缘时，|▽I|＞＞k，沿梯度方向逆扩散；在图像平坦区域时，|▽I|<<k，沿着边缘方向和梯度方向同时进行平滑处理，其中f(·)≥0，k是图像I的梯度阈值；

步骤四、建立图像的梯度阈值函数k＝e^-αt，其为随扩散时间和扩散次数变化的函数，其中t＝n，t为扩散时间，n为扩散次数，max{W,H}为图像的最大灰度值；

步骤五、引入保真项建立基于能量泛函的图像平滑与锐化算法

$\left\{\begin{array}{c}\frac{\&PartialD;I}{\&PartialD;t}=\frac{1+k^2}{{|\&dtri;I|}^2+k^2}{I}_{\mathrm{\&xi;\&xi;}}+\frac{(1+k^2)}{{|\&dtri;I|}^2+k^2}\frac{(k^2-{|\&dtri;I|}^2)}{{|\&dtri;I|}^2+k^2}{I}_{\mathrm{\&eta;\&eta;}}-\mathrm{\&lambda;}(1-\frac{1}{{|\&dtri;I|}^2+1})(I-I_0)\\ I(x,y,0)=I_0\end{array}\right.---4)$

其中λ是调整参数，在图像的边缘，|▽I|→∞，在图像的平坦区域，|▽I|→0，则扩散程度达到最大；

步骤六、用中心差分数值算法对步骤五的结果进一步处理。

所述步骤六的中心差分数值算法具体为，令

$T=\mathrm{div}\left(\frac{\&dtri;I}{|\&dtri;I|}\right)=\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}}}$

$P=\mathrm{div}\left(g\right|\&dtri;I|\&dtri;I)=\frac{k^2{I}_{\mathrm{xx}}+k^2{I}_{\mathrm{yy}}-I_x^2{I}_{\mathrm{yy}}+I_y^2{I}_{\mathrm{xx}}-4I_x{I}_y{I}_{\mathrm{xy}}+I_x^2{I}_{\mathrm{yy}}-I_y^2{I}_{\mathrm{yy}}}{{\left(\frac{k^2+I_x^2+I_y^2}{k}\right)}^2}$

式4)的离散化形式为

$\frac{I^{n+1}-I^n}{\mathrm{\&Delta;I}}=\frac{1+k^2}{{\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n+k^2}{\left(I_{\mathrm{\&xi;\&xi;}}\right)}_{\mathrm{ij}}^n+\frac{(1+k^2)(k^2-{\left(I_x^2\right)}_{\mathrm{ij}}^n-{\left(I_y^2\right)}_{\mathrm{ij}}^n)}{{({\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n+k^2)}^2}{\left(I_{\mathrm{\&eta;\&eta;}}\right)}_{\mathrm{ij}}^n-\mathrm{\&lambda;}(1-\frac{1}{{\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n+1})(I-I_0)$

式中，n＝0,1,2,……,表示分解尺度，则式4)的离散形式为

$I^{n+1}-I^n=(\frac{1+k^2}{{\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n+k^2}{\left(I_{\mathrm{\&xi;\&xi;}}\right)}_{\mathrm{ij}}^n+\frac{(1+k^2)(k^2-{\left(I_x^2\right)}_{\mathrm{ij}}^n-{\left(I_y^2\right)}_{\mathrm{ij}}^n)}{{({\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n+k^2)}^2}{\left(I_{\mathrm{\&eta;\&eta;}}\right)}_{\mathrm{ij}}^n-\mathrm{\&lambda;}(1-\frac{1}{{\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n+1})(I-I_0))\mathrm{\&Delta;t}.$

本发明达到的有益效果：在复杂度方面，利用内在坐标形式进行性能分析，得到扩散滤波算法，算法需要的信息量少，方法简单，实现图像去噪，使受污染图像更接近原始图像；在时效性方面，因为需要的信息量少，实施的复杂度低，从而降低了算法的处理时间；在去噪性能方面，图像的峰值信噪比大幅提高，受噪声污染的图像经本算法处理后更加接近原始图像。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

第一步：建立噪声图像模型I₀＝KI，为了使去噪后的图像与原始图像充分接近，现建立最小化模型求该式的拉格朗日方程，建立动态方程引入Neumann边界，从而动态方程就等价于各向异性扩散方程 $\left\{\begin{array}{c}\frac{\&PartialD;I}{\&PartialD;t}=\mathrm{div}(g\left(\right|\&dtri;I\left|\right)\&dtri;I)\\ \frac{\&PartialD;I}{\&PartialD;n}=0\\ I(x,y,0)=I_0\end{array}\right.,$ 这样从最小化问题出发得到了各向异性扩散方程，即各向异性扩散方程可以归结为最小化能量泛函问题，其中I是无噪图像，I₀是降质图像，K表示降质算子，；

第二步：建立第一步中的内在坐标形式由于第一步中的最小化模型有解，则故沿边缘方向的扩散系数在扩散过程中沿边缘方向正扩散，模糊图像；而梯度方向的扩散系数其正负号不确定，当时，沿梯度方向正扩散，模糊图像，当时，沿梯度方向逆扩散，锐化图像的边缘，能很好的保护图像的边缘。其中(η,ξ)为内在坐标系，η为图像的梯度方向，即垂直于图像特征(边缘)的方向；ξ为垂直于梯度的方向，即沿图像特征(边缘)的方向， $\mathrm{\&eta;}=\frac{(I_x,I_y)}{\sqrt{I_x^2+I_y^2}},\mathrm{\&xi;}=\frac{(-I_y,I_x)}{\sqrt{I_x^2+I_y^2}},I_{\mathrm{\&xi;\&xi;}}=\frac{I_x^2{I}_{\mathrm{xx}}-2I_x{I}_y{I}_{\mathrm{xy}}\text{+}{I}_x^2{I}_{\mathrm{yy}}}{I_x^2+I_y^2},$ $I_{\mathrm{\&eta;\&eta;}}=\frac{I_x^2{I}_{\mathrm{xx}}+2I_x{I}_y{I}_{\mathrm{xy}}\text{+}{I}_y^2{I}_{\mathrm{yy}}}{I_x^2+I_y^2};$

第三步：建立能量泛函由第一步解该式的最小能量泛函，建立数学模型由第二步建立数学模型的内在坐标形式 $\frac{\&PartialD;I}{\&PartialD;t}=g\left(\right|\&dtri;I\left|\right)I_{\mathrm{\&xi;\&xi;}}+{\left(\right|\&dtri;I\left|g\left(\right|\&dtri;I\left|\right)\right)}^{\&prime;}{I}_{\mathrm{\&eta;\&eta;}},$ 现取 $g\left(\right|\&dtri;I\left|\right)=\frac{1+k^2}{{|\&dtri;I|}^2+k^2},$ 代入数学模型中 $\frac{\&PartialD;I}{\&PartialD;t}=\frac{1+k^2}{{|\&dtri;I|}^2+k^2}{I}_{\mathrm{\&xi;\&xi;}}+\frac{(1+k^2)}{{|\&dtri;I|}^2+k^2}\frac{(k^2-{|\&dtri;I|}^2)}{{|\&dtri;I|}^2+k^2}{I}_{\mathrm{\&eta;\&eta;}},$ 根据第二步，当在图像边缘时，|▽I|＞＞k，沿梯度方向逆扩散，锐化图像的边缘，改善了边缘模糊的现象，能很好的保护图像的边缘；当在图像平坦区域时，|▽I|<<k，在平坦区域，沿着边缘方向和切线方向同时进行平滑，转化为各向同性的扩散方程，因此具有较强的去噪能力。其中Ω为图像区域，f(·)≥0，f′(·)＞0，k是图像I的梯度阈值，下同；

第四步：建立图像的梯度阈值函数k。根据第三步本发明的整个扩散过程，实质上是图像的梯度与阈值大小进行比较实现的，故将梯度阈值k设计为一个随扩散时间和扩散次数变化的函数k＝e^-αt，其中扩散次数代表扩散时间，即t＝n，max{W,H}为图像的最大灰度值。随着扩散次数增加，即扩散时间增加，k不断减小，即梯度阈值不断减小，这样就使得下一次的边缘判断更加准确，从而保留了更多的边缘信息；

第五步：建立保真项使图像不失真。其中λ是调整参数，一般取较小的值，有重要作用，在图像的边缘，此时|▽I|→∞， 这样滤波结果会尽可能地接近原始图像，进一步加强了保边缘的效果；在图像的平坦区域，此时|▽I|→0，则扩散程度达到最大，尽可能的去除噪声。在此基础上，建立基于能量泛函的图像平滑与锐化算法

$\left\{\begin{array}{c}\frac{\&PartialD;I}{\&PartialD;t}=\frac{1+k^2}{{|\&dtri;I|}^2+k^2}{I}_{\mathrm{\&xi;\&xi;}}+\frac{(1+k^2)}{{|\&dtri;I|}^2+k^2}\frac{(k^2-{|\&dtri;I|}^2)}{{|\&dtri;I|}^2+k^2}{I}_{\mathrm{\&eta;\&eta;}}-\mathrm{\&lambda;}(1-\frac{1}{{|\&dtri;I|}^2+1})(I-I_0)\\ I(x,y,0)=I_0\end{array}\right.;$

第六步：根据第五步得到的基于能量泛函的图像平滑与锐化算法，用中心差分数值算法进行数值计算。

为了实现本发明算法,我们采用有限差分法，令

$T=\mathrm{div}\left(\frac{\&dtri;I}{|\&dtri;I|}\right)=\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}}}$

$P=\mathrm{div}\left(g\right|\&dtri;I|\&dtri;I)=\frac{k^2{I}_{\mathrm{xx}}+k^2{I}_{\mathrm{yy}}-I_x^2{I}_{\mathrm{yy}}+I_y^2{I}_{\mathrm{xx}}-4I_x{I}_y{I}_{\mathrm{xy}}+I_x^2{I}_{\mathrm{yy}}-I_y^2{I}_{\mathrm{yy}}}{{\left(\frac{k^2+I_x^2+I_y^2}{k}\right)}^2}$

式中，k为图像的梯度阈值。现采用中心差分法，则本发明算法的离散化形式为

$\frac{I^{n+1}-I^n}{\mathrm{\&Delta;I}}=\frac{1+k^2}{{\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n+k^2}{\left(I_{\mathrm{\&xi;\&xi;}}\right)}_{\mathrm{ij}}^n+\frac{(1+k^2)(k^2-{\left(I_x^2\right)}_{\mathrm{ij}}^n-{\left(I_y^2\right)}_{\mathrm{ij}}^n)}{{({\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n+k^2)}^2}{\left(I_{\mathrm{\&eta;\&eta;}}\right)}_{\mathrm{ij}}^n-\mathrm{\&lambda;}(1-\frac{1}{{\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n+1})(I-I_0)$

式中，n＝0,1,2,……,表示时间水平。其中

$\frac{I^{n+1}-I^n}{\mathrm{\&Delta;t}}=\frac{I_{\mathrm{ij}}^{n+1}-I_{\mathrm{ij}}^n}{\mathrm{\&Delta;t}}$

${\left(I_{\mathrm{\&xi;\&xi;}}\right)}_{\mathrm{ij}}^n=\frac{{\left(I_x^2\right)}_{\mathrm{ij}}^n{\left(I_{\mathrm{yy}}\right)}_{\mathrm{ij}}^n-2{\left(I_x\right)}_{\mathrm{ij}}^n{\left(I_y\right)}_{\mathrm{ij}}^n{\left(I_{\mathrm{xy}}\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n{\left(I_{\mathrm{xx}}\right)}_{\mathrm{ij}}^n}{({\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n)}$

${\left(I_{\mathrm{\&eta;\&eta;}}\right)}_{\mathrm{ij}}^n=\frac{{\left(I_x^2\right)}_{\mathrm{ij}}^n{\left(I_{\mathrm{xx}}\right)}_{\mathrm{ij}}^n-2{\left(I_x\right)}_{\mathrm{ij}}^n{\left(I_y\right)}_{\mathrm{ij}}^n{\left(I_{\mathrm{xy}}\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n{\left(I_{\mathrm{yy}}\right)}_{\mathrm{ij}}^n}{({\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n)}$

$T_{\mathrm{ij}}^n=\frac{{\left(I_x^2\right)}_{\mathrm{ij}}^n{\left(I_{\mathrm{yy}}\right)}_{\mathrm{ij}}^n-2{\left(I_x\right)}_{\mathrm{ij}}^n{\left(I_y\right)}_{\mathrm{ij}}^n{\left(I_{\mathrm{xy}}\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n{\left(I_{\mathrm{xx}}\right)}_{\mathrm{ij}}^n}{{({\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n)}^{\frac{3}{2}}}$

$P_{\mathrm{ij}}^n=\frac{k^2{\left(I_{\mathrm{xy}}\right)}_{\mathrm{ij}}^n+k^2{\left(I_{\mathrm{yy}}\right)}_{\mathrm{ij}}^n-{\left(I_x^2\right)}_{\mathrm{ij}}^n{\left(I_{\mathrm{xx}}\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n{\left(I_{\mathrm{xx}}\right)}_{\mathrm{ij}}^n-4{\left(I_x\right)}_{\mathrm{ij}}^n{\left(I_y\right)}_{\mathrm{ij}}^n{\left(I_{\mathrm{xy}}\right)}_{\mathrm{ij}}^n+{\left(I_x^2\right)}_{\mathrm{ij}}^n{\left(I_{\mathrm{yy}}\right)}_{\mathrm{ij}}^n-{\left(I_y^2\right)}_{\mathrm{ij}}^n{\left(I_{\mathrm{yy}}\right)}_{\mathrm{ij}}^n}{{\left(\frac{k^2+{\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n}{k}\right)}^2}$

最终本发明算法的离散形式为

$I^{n+1}-I^n=(\frac{1+k^2}{{\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n+k^2}{\left(I_{\mathrm{\&xi;\&xi;}}\right)}_{\mathrm{ij}}^n+\frac{(1+k^2)(k^2-{\left(I_x^2\right)}_{\mathrm{ij}}^n-{\left(I_y^2\right)}_{\mathrm{ij}}^n)}{{({\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n+k^2)}^2}{\left(I_{\mathrm{\&eta;\&eta;}}\right)}_{\mathrm{ij}}^n-\mathrm{\&lambda;}(1-\frac{1}{{\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n+1})(I-I_0))\mathrm{\&Delta;t}.$

综上所述，本发明专利所述基于能量泛函的图像平滑与锐化算法复杂度低，时效性高，去噪后的图像峰值信噪比与PM算法相比提高了15个dB左右，与TV算法相比提高了8个dB左右，能够增强与锐化图像中的边缘纹理等重要信息。

用受高斯随机噪声(σ＝20)污染的图进行仿真实验，把本发明算法与经典的PM算法与TV算法进行比较，由表1和表2，可以看出，本发明算法较PM算法的峰值信噪比提高了15个dB左右，较TV算法提高了8个dB左右，进一步证实了本发明算法的有效性与合理性。

表1各去噪算法的MSE和PSNR比较

表2各去噪算法不同噪声方差的MSE和PSNR比较

以上是本发明的较佳实施方式，但本发明的保护范围不限于此。任何熟悉本领域的技术人员在本发明所揭露的技术范围内，未经创造性劳动想到的变换或替换，都应涵盖在本发明的保护范围之内。因此本发明的保护范围应以权利要求所限定的保护范围为准。

Claims (2)

1.基于能量泛函的图像平滑与锐化算法，其特征在于，包括以下步骤：

步骤一、建立噪声图像模型I₀＝KI，为使去噪后的图像与原始图像充分接近，建立所述噪声图像模型的最小化模型其动态方程为引入Neumann边界，该动态方程等价于各向异性扩散方程 $\left\{\begin{array}{c}\frac{\&PartialD;I}{\&PartialD;t}=\mathrm{div}(g\left(\right|\&dtri;I\left|\right)\&dtri;I)\\ \frac{\&PartialD;I}{\&PartialD;n}=0\\ I(x,y,0)=I_0\end{array}\right.,$ 即各向异性扩散方程转化为最小化能量泛函，其中I是无噪图像，I₀是降质图像，K是降质算子，Ω是图像区域；

步骤二、将步骤一中的动态方程以内在坐标形式表示为

步骤三、建立能量泛函式1)表示为

$\frac{\&PartialD;I}{\&PartialD;t}=g\left(\right|\&dtri;I\left|\right)I_{\mathrm{\&xi;\&xi;}}+\left(\right|\&dtri;I\left|g\left(\right|\&dtri;I\left|\right)\right)\&prime;I_{\mathrm{\&eta;\&eta;}}---2)$

取代入式2)得

$\frac{\&PartialD;I}{\&PartialD;t}=\frac{1+k^2}{{|\&dtri;I|}^2+k^2}{I}_{\mathrm{\&xi;\&xi;}}+\frac{(1+k^2)}{{|\&dtri;I|}^2+k^2}\frac{(k^2-{|\&dtri;I|}^2)}{{|\&dtri;I|}^2+k^2}{I}_{\mathrm{\&eta;\&eta;}}---3)$

在图像边缘时，|▽I|>>k，沿梯度方向逆扩散；在图像平坦区域时，|▽I|<<k，沿着边缘方向和梯度方向同时进行平滑处理，其中f(·)≥0，k是图像I的梯度阈值；

步骤四、建立图像的梯度阈值函数k＝e^-αt，其为随扩散时间和扩散次数变化的函数，其中t＝n，t为扩散时间，n为扩散次数，max{W,H}为图像的最大灰度值；

步骤五、引入保真项建立基于能量泛函的图像平滑与锐化算法

$\left[\begin{array}{c}\frac{\&PartialD;I}{\&PartialD;t}=\frac{1+k^2}{{|\&dtri;I|}^2+k^2}{I}_{\mathrm{\&xi;\&xi;}}+\frac{(1+k^2)}{{|\&dtri;I|}^2+k^2}\frac{(k^2-{|\&dtri;I|}^2)}{{|\&dtri;I|}^2+k^2}{I}_{\mathrm{\&eta;\&eta;}}-\mathrm{\&lambda;}(1-\frac{1}{{|\&dtri;I|}^2+1})(I-I_0)---4)\\ I(x,y,0)=I_0\end{array}\right]$

其中λ是调整参数，在图像的边缘， $|\&dtri;I|\&RightArrow;\&infin;,\frac{1}{{|\&dtri;I|}^2+1}\&RightArrow;0,(1-\frac{1}{{|\&dtri;I|}^2+1})\&ap;1;$ 在图像的平坦区域， $|\&dtri;I|\&RightArrow;0,\frac{1}{{|\&dtri;I|}^2+1}\&RightArrow;1,(1-\frac{1}{{|\&dtri;I|}^2+1})\&ap;0;$ 则扩散程度达到最大；

步骤六、用中心差分数值算法对步骤五的结果进一步处理。

2.根据权利要求1所述的基于能量泛函的图像平滑与锐化算法，其特征在于，所述步骤六的中心差分数值算法具体为，令

$T=\mathrm{div}\left(\frac{\&dtri;I}{|\&dtri;I|}\right)=\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}}}$

$P=\mathrm{div}\left(g\right|\&dtri;I|\&dtri;I)=\frac{k^2{I}_{\mathrm{xx}}+k^2{I}_{\mathrm{yy}}-I_x^2{I}_{\mathrm{xx}}+I_y^2{I}_{\mathrm{xx}}-4I_x{I}_y{I}_{\mathrm{xy}}+I_x^2{I}_{\mathrm{yy}}-I_y^2{I}_{\mathrm{yy}}}{{\left(\frac{k^2+I_x^2+I_y^2}{k}\right)}^2}$

式4)的离散化形式为

$\frac{I^{n+1}-I^n}{\mathrm{\&Delta;I}}=\frac{1+k^2}{{\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n+k^2}{\left(I_{\mathrm{\&xi;\&xi;}}\right)}_{\mathrm{ij}}^n+\frac{(1+k^2)(k^2-{\left(I_x^2\right)}_{\mathrm{ij}}^n-{\left(I_y^2\right)}_{\mathrm{ij}}^n)}{{({\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n+k^2)}^2}{\left(I_{\mathrm{\&eta;\&eta;}}\right)}_{\mathrm{ij}}^n-\mathrm{\&lambda;}(1-\frac{1}{{\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n+1})(I-I_0)$

式中，n＝0,1,2,······,表示分解尺度，则式4)的离散形式为

$I^{n+1}-I^n=(\frac{1+k^2}{{\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n+k^2}{\left(I_{\mathrm{\&xi;\&xi;}}\right)}_{\mathrm{ij}}^n+\frac{(1+k^2)(k^2-{\left(I_x^2\right)}_{\mathrm{ij}}^n-{\left(I_y^2\right)}_{\mathrm{ij}}^n)}{{({\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n+k^2)}^2}{\left(I_{\mathrm{\&eta;\&eta;}}\right)}_{\mathrm{ij}}^n-\mathrm{\&lambda;}(1-\frac{1}{{\left(I_x^2\right)}_{\mathrm{ij}}^n+{\left(I_y^2\right)}_{\mathrm{ij}}^n+1})(I-I_0))\mathrm{\&Delta;t}.$