CN105184766A - 一种频域边界能量模型的水平集图像分割方法 - Google Patents

Abstract

本发明提供了一种频域边界能量模型的水平集图像分割方法，包括：在初始化轮廓下，对输入图像的轮廓内外区域分别进行傅里叶变换；设计频域滤波器对轮廓内外区域进行滤波从而提取边界能量项；将水平集函数正则项、长度约束项和边界能量项作为整体能量项；结合阶跃函数与冲击函数，运用水平集方法对能量项进行最小化，从而得到理想的分割结果。本发明所述的频域边界能量模型的水平集分割算法，选择带通滤波器作为最优滤波器提取目标边界信息，而排除其他干扰信息，如强噪声、偏移场、规则纹理等，能够对分割结果边界的平滑度进行控制，对同时存在强噪声和偏移场的医学图像有高精度的分割结果。

Description

一种频域边界能量模型的水平集图像分割方法

技术领域

本发明涉及图像处理技术领域，具体涉及一种能够有效分割图像的频域边界能量模型的水平集图像分割方法。

背景技术

图像分割一直是图像处理领域的研究热点和难点，其目的是将图像分割为具有相同特性的若干子区域，如亮度、纹理、密度等。图像分割是进行高级图像如图像分析、目标识别、图像理解和计算机视觉的关键步骤。图像分割在医学领域、安防系统、军事领域等诸多方面都有广泛应用。虽然经过多年研究，许多学者已经提出大量有效的图像分割算法，但图像分割没有一种通用的有效算法，针对不同的应用场景，往往需要进行研究改进现有算法从而达到最好效果。在众多图像分割算法中，基于活动轮廓模型的方法能够获取亚像素精度结果以及平滑的轮廓而被广泛应用于医学图像分割。然而，在医学图像中，由于成像技术水平的限制，获取的图像常常包含了强噪声或偏移场。而且，两种干扰信息常常同时出现，例如超声波图像，光学相干层析图像(OpticalCoherenceTomography，OCT)图像等。这些干扰信息严重影响了医学图像的精细分割，因此实现一种能够消除强噪声和偏移场干扰的分割方法是十分必要的。

发明内容

发明目的：本发明所要解决的技术问题是针对现有技术的不足，提供一种频域边界能量模型的水平集图像分割方法。

为了解决上述技术问题，本发明公开了一种频域边界能量模型的水平集图像分割方法，在频域内，将图像分割问题归结为如何设计滤波器有效地提取边界能量信息，包括以下步骤：

步骤1，在图像区域Ω范围内，对初始化轮廓外区域和内区域进行傅里叶变换。

步骤2，利用设计的滤波器，如带通滤波器，提取图像的边界能量项。

步骤3，融合水平集函数正则化项，轮廓约束项及边界能量项，获得整体的能量项。

步骤4，实例化阶跃函数H(·)函数和冲击函数δ(·)，计算使得能量函数变小的水平集函数的变化值。对水平集函数进行更新。若能量函数不收敛，则返回步骤3；否则输出分割结果。

其中，步骤1包括以下步骤：

步骤1-1，对于定义在Ω上的输入图像u₀，输入图像u₀上的每个像素对应一个Ω上的坐标向量x。在Ω上设置初始化轮廓C₀，例如可以采用通过交互式方法。设t时刻，在Ω上的初始化轮廓演化为C_t，将C_t用水平集函数Φ表示，即C_t＝{x|Φ(x)＝0}。C_t将Ω划分为两个子区域，即C_t外区域Ω₁及内区域Ω₂。利用阶跃函数H(·)及水平集函数Φ(·)，Ω₁及Ω₂可以分别表示为u₁和u₂，如下式所示：

u₁(x)＝u₀(x)H(φ(x)),u₂(x)＝u₀(x)(1-H(φ(x)))。

步骤1-2,在Ω范围内分别对u₁及u₂进行傅里叶变换，得到变换结果：

其中，s为频域坐标向量，在Ω上对u_i进行傅里叶变换，F_i(s)为变换结果。

步骤2包括以下步骤：

步骤2-1，设置带通滤波器作为最优的边界滤波器。使用DoG算子作为带通滤波器，公式如下：

H_bp(s)＝DoG(s)＝G(s,σ_s1)-G(s,σ_s2)，

$G(s,{\mathrm{\σ}}_{si})=\frac{1}{2{{\mathrm{\π\σ}}_{si}}^2}\exp (-\frac{s^{T}s}{2{\mathrm{\σ}}_{si}^2}),$

G(s,σ_si)为变量为s标准差为σ_si的高斯滤波器。σ_s1和σ_s2分别为两个高斯滤波器的标准差，满足：σ_s1＞0，σ_s2＞0；

步骤2-2，使用频域滤波器H_bp(s)对F_i(s)进行滤波，并进行傅里叶逆变换得到目标的边界信息。

其中，f_b1(x)为图像区域上C_t外部的边界信息，f_b2(x)为图像区域上C_t内部的边界信息；为傅里叶逆变换。

步骤2-3将C_t表示为水平集函数，计算边界能量项E^FBE(φ,f_bi(x))：

$E^{FBE}(\mathrm{\φ},f_{bi}(x\left)\right)=\underset{u_1,\mathrm{\Ω}}{\∫}{f_{b1}}^2\left(x\right)H\left(\mathrm{\φ}\right)dx+\underset{u_2,\mathrm{\Ω}}{\∫}{f_{b2}}^2\left(x\right)(1-H(\mathrm{\φ}\left)\right)dx,$

其中，Φ为C_t对应的水平集函数。；H(·)为阶跃函数。

步骤3包括以下步骤：

步骤3-1，引入水平集函数正则化项，水平集函数的正则化项如下式所示：

$E^{LSR}={\▿}^2\mathrm{\φ}-div\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right),$

其中▽为梯度运算，div为散度运算。

步骤3-2，引入轮廓的长度约束项，长度约束项如下式所示：

$E^{Len}=\mathrm{\δ}\left(\mathrm{\φ}\right)div\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right),$

其中，δ(φ)为冲击函数；

步骤3-3，融合水平集函数正则项、轮廓长度约束项以及定义的频域边界能量项，则整体的能量函数为：

E(φ,f_b1(x))＝λE^FBE(φ,f_b1(x))+μE^LSR(φ)+νE^Len(φ)，

其中，λ，μ，ν分别为三个能量项的权值，满足：0<λ，0<μ，0<ν。

步骤4包括以下步骤：

步骤4-1：将阶跃函数和冲击函数分别实例化为，当ε趋于无限大时，两个函数都趋于理想的阶跃函数和冲击函数。

$H_{\mathrm{\&epsiv;}}\left(z\right)=\frac{1}{2}\&lsqb;1+\frac{2}{\mathrm{\&pi;}}arctan\left(\frac{z}{\mathrm{\&epsiv;}}\right)\&rsqb;,$

${\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(z\right)=\frac{1}{\mathrm{\&pi;}}\frac{\mathrm{\&epsiv;}}{{\mathrm{\&epsiv;}}^2+z^2};$

其中，z与ε是表示上述函数的辅助变量。当ε趋于无限大时，两个函数都趋于理想的阶跃函数H(z)和冲击函数δ(z)；步骤4-2，计算频域边界能量项下，水平集函数的更新方程，对水平集函数引入辅助变量α及η，则：固定f_b1(x)、f_b2(x)，对Φ求偏导，使α趋近于0：

$\begin{array}{c}\frac{\&part;E^{FBE}}{\&part;\mathrm{\&phi;}}=\underset{\mathrm{\&alpha;}\&RightArrow;0}{lim}\frac{d}{d\mathrm{\&alpha;}}{\left(\underset{\mathrm{\&Omega;}}{\&Integral;}{f}_{b1}\right(x\left)H\right(\widetilde{\mathrm{\&phi;}}\left)\right)}^2+{\left(f_{b2}\right(x\left)\right(1-H\left(\widetilde{\mathrm{\&phi;}}\right)\left)\right)}^{2}dx)\\ =\underset{\mathrm{\&alpha;}\&RightArrow;0}{lim}(2\underset{\mathrm{\&Omega;}}{\&Integral;}\&lsqb;{f_{b1}}^{2}H(\widetilde{\mathrm{\&phi;}})-{f_{b2}}^2\left(x\right)(1-H(\widetilde{\mathrm{\&phi;}}))\&rsqb;\mathrm{\&delta;}\left(\widetilde{\mathrm{\&phi;}}\right)\mathrm{\&eta;}dx\\ =2\underset{\mathrm{\&Omega;}}{\&Integral;}\&lsqb;{f_{b1}}^2\left(x\right)H\left(\mathrm{\&phi;}\right)-{f_{b2}}^2\left(x\right)(1-H(\mathrm{\&phi;}))\&rsqb;\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)\mathrm{\&eta;}dx\end{array}$

根据欧拉-拉格朗日公式，得到：

[f_b1 ²(x)H(φ)-f_b2 ²(x)(1-H(φ))]δ(φ)＝0，

运用梯度下降法则得到：

$\frac{\&part;\mathrm{\&phi;}}{\&part;t}=\&lsqb;{f_{b1}}^2\left(x\right)H\left(\mathrm{\&phi;}\right)-{f_{b2}}^2\left(x\right)(1-H(\mathrm{\&phi;}\left)\right)\&rsqb;\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)$

步骤4-3，结合水平集函数正则项和长度约束项的水平集函数更新方式以及实际使用的阶跃函数和冲击函数，得到了水平集函数的最终更新值：

上式为边界能量项； $\mathrm{\&mu;}({\&dtri;}^2\mathrm{\&phi;}-\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right))$ 为水平函数正则化项；为轮廓长度约束项；

步骤4-4，利用上述最终更新值，对水平集函数进行更新，并判断水平集函数的变化是否小于特定值，若小于，则收敛，输出分割结果，该特定值取值范围为自然数；否则返回步骤4-3重新计算水平集函数的更新值，利用更新后水平集函数重新计算新的能量函数。

本发明根据图像边界的高频能量，定义高频能量项；运用水平集函数，结合FBE-DoG模型，(FBE，FrequencyBoundaryEnergy，频域边界能量)计算理想分割边界。本发明提供了一种新的图像分割思路，选择带通滤波器作为最优滤波器提取目标边界信息，而排除其他干扰信息，如强噪声，偏移场，纹理等，能够对分割结果边界的平滑度进行控制。高频能量项的提出，结合水平集对图像分割，对同时存在强噪声和偏移场的医学图像有高精度的分割结果。仿真实验结果表明，该方法不仅能有效地对图像进行分割，较传统的图像分割方法，该方法也有较好的效果。

附图说明

下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述和/或其他方面的优点将会变得更加清楚。

图1为输入的待分割图像。

图2为初始化轮廓示意；其中红白虚线为初始化轮廓。

图3a为迭代过程中的C_t；

图3b为C_t对应的水平集函数。

图4a为C_t取图3a所示时的u₁；

图4b为C_t取图3a所示时的u₂；

图5为带通滤波器特性。

图6为迭代过程中，能量值的变化。

图7为图1的最终分割结果。

图8为不同分割算法对比图。

图9为本发明流程图。

具体实施方式

实施例

由于本发明的特殊性，本实施例中图1～图8不可避免的使用到灰度图。但是本领域技术人员根据本发明文字的描述不依赖实施例的附图内容就可以理解本发明内容。

如图9所示，本发明公开了一种频域边界能量模型的水平集图像分割方法，包括以下步骤：步骤1，在图像区域Ω范围内，对初始化轮廓内外区域进行傅里叶变换。详细步骤如下：

步骤1-1，对于定义在Ω上的输入图像u₀(如图1所示)，图像u₀上的每个像素对应一个Ω上的坐标向量x。通过交互式方法在Ω上设置初始化轮廓C₀，如图2所示。设t时刻，在Ω上的初始化轮廓演化为C_t，如图3a所示。将C_t用水平集函数Φ表示，即C_t＝{x|Φ(x)＝0}，如图3b所示。C_t将Ω划分为两个子区域，即C_t外区域Ω₁及内区域Ω₂。利用阶跃函数H(·)及水平集函数Φ(·)，Ω₁及Ω₂可以分别表示为u₁和u₂，分别如图4a，图4b所示，公式如下所示：

u₁(x)＝u₀(x)H(φ(x)),u₂(x)＝u₀(x)(1-H(φ(x)))。

步骤1-2,在Ω范围内分别对u₁及u₂进行傅里叶变换，得到变换结果：

其中，s为频域坐标向量，在Ω上对u_i进行傅里叶变换，F_i(s)为变换结果。

步骤2，利用设计的滤波器，如带通滤波器，提取图像的边界能量项。详细步骤如下：

步骤2-1，设置带通滤波器作为最优的边界滤波器。使用DoG算子作为带通滤波器，如图5所示，公式如下：

H_bp(s)＝DoG(s)＝G(s,σ_s1)-G(s,σ_s2)，

$G(s,{\mathrm{\&sigma;}}_{si})=\frac{1}{2{{\mathrm{\&pi;\&sigma;}}_{si}}^2}\exp (-\frac{s^{T}s}{2{\mathrm{\&sigma;}}_{si}^2}),$

G(s,σ_si)为变量为s标准差为σ_si的高斯滤波器。σ_s1和σ_s2分别为两个高斯滤波器的标准差，满足：σ_s1＞0，σ_s2＞0

步骤2-2，使用频域滤波器H_bp(s)对F_i(s)进行滤波，并进行傅里叶逆变换得到目标的边界信息。

其中，f_b1(x)为图像区域上C_t外部的边界信息，f_b2(x)为图像区域上C_t内部的边界信息；为傅里叶逆变换。

步骤2-3将C_t表示为水平集函数，计算边界能量项E^FBE(φ,f_bi(x))：

$E^{FBE}(\mathrm{\&phi;},f_{bi}(x\left)\right)=\underset{u_1,\mathrm{\&Omega;}}{\&Integral;}{f_{b1}}^2\left(x\right)H\left(\mathrm{\&phi;}\right)dx+\underset{u_2,\mathrm{\&Omega;}}{\&Integral;}{f_{b2}}^2\left(x\right)(1-H(\mathrm{\&phi;}\left)\right)dx,$

其中，Φ为C_t对应的水平集函数；H(·)为阶跃函数。

步骤3，融合水平集函数正则化项，轮廓约束项及边界能量项，获得整体的能量项。详细步骤如下：

步骤3-1，引入水平集函数正则化项，水平集函数的正则化项如下式所示：

$E^{LSR}={\&dtri;}^2\mathrm{\&phi;}-div\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right),$

其中▽为梯度运算，div为散度运算。

步骤3-2，引入轮廓的长度约束项，长度约束项如下式所示：

$E^{Len}=\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)div\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right),$

其中，δ(φ)为冲击函数；

步骤3-3，融合水平集函数正则项、轮廓长度约束项以及定义的频域边界能量项，则整体的能量函数为：

E(φ,f_b1(x))＝λE^FBE(φ,f_b1(x))+μE^LSR(φ)+νE^Len(φ)，

其中，λ，μ，ν分别为三个能量项的权值，满足：0<λ，0<μ，0<ν。

步骤4，实例化阶跃函数H(·)函数和冲击函数δ(·)，计算使得能量函数变小的水平集函数的变化值。对水平集函数进行更新。若能量函数不收敛，则返回步骤3；否则输出分割结果,迭代过程中能量的变化，如图6所示。详细步骤如下：

步骤4-1：将阶跃函数和冲击函数分别实例化为，当ε趋于无限大时，两个函数都趋于理想的阶跃函数和冲击函数。

$H_{\mathrm{\&epsiv;}}\left(z\right)=\frac{1}{2}\&lsqb;1+\frac{2}{\mathrm{\&pi;}}arctan\left(\frac{z}{\mathrm{\&epsiv;}}\right)\&rsqb;,$

${\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(z\right)=\frac{1}{\mathrm{\&pi;}}\frac{\mathrm{\&epsiv;}}{{\mathrm{\&epsiv;}}^2+z^2};$

其中，z与ε是表示上述函数的辅助变量。当ε趋于无限大时，两个函数都趋于理想的阶跃函数H(z)和冲击函数δ(z)；

步骤4-2，计算频域边界能量项下，水平集函数的更新方程，对水平集函数引入辅助变量α及η，则：固定f_b1(x)、f_b2(x)，对Φ求偏导，使α趋近于0：

$\begin{array}{c}\frac{\&part;E^{FBE}}{\&part;\mathrm{\&phi;}}=\underset{\mathrm{\&alpha;}\&RightArrow;0}{lim}\frac{d}{d\mathrm{\&alpha;}}{\left(\underset{\mathrm{\&Omega;}}{\&Integral;}{f}_{b1}\right(x\left)H\right(\widetilde{\mathrm{\&phi;}}\left)\right)}^2+{\left(f_{b2}\right(x\left)\right(1-H\left(\widetilde{\mathrm{\&phi;}}\right)\left)\right)}^{2}dx)\\ =\underset{\mathrm{\&alpha;}\&RightArrow;0}{lim}(2\underset{\mathrm{\&Omega;}}{\&Integral;}\&lsqb;{f_{b1}}^{2}H(\widetilde{\mathrm{\&phi;}})-{f_{b2}}^2\left(x\right)(1-H(\widetilde{\mathrm{\&phi;}}))\&rsqb;\mathrm{\&delta;}\left(\widetilde{\mathrm{\&phi;}}\right)\mathrm{\&eta;}dx\\ =2\underset{\mathrm{\&Omega;}}{\&Integral;}\&lsqb;{f_{b1}}^2\left(x\right)H\left(\mathrm{\&phi;}\right)-{f_{b2}}^2\left(x\right)(1-H(\mathrm{\&phi;}))\&rsqb;\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)\mathrm{\&eta;}dx\end{array}$

根据欧拉-拉格朗日公式，得到：

[f_b1 ²(x)H(φ)-f_b2 ²(x)(1-H(φ))]δ(φ)＝0，

运用梯度下降法则得到：

$\frac{\&part;\mathrm{\&phi;}}{\&part;t}=\&lsqb;{f_{b1}}^2\left(x\right)H\left(\mathrm{\&phi;}\right)-{f_{b2}}^2\left(x\right)(1-H(\mathrm{\&phi;}\left)\right)\&rsqb;\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)$

步骤4-3，结合水平集函数正则项和长度约束项的水平集函数更新方式以及实际使用的阶跃函数和冲击函数，得到了水平集函数的最终更新值：

上式为边界能量项；为水平函数正则化项；为轮廓长度约束项；

步骤4-4，利用上述最终更新值，对水平集函数进行更新，并判断水平集函数的变化是否趋于稳定。这里将连续5次迭代变化的方差大小作为是否稳定的标准。若变化的方差小于特定值(比如0.01倍方差)，则收敛，输出分割结果，如图7所示；否则返回步骤4-3重新计算水平集函数的更新值，利用更新后水平集函数重新计算新的能量函数。与其他6种方法的对比结果如图8所示。图8为不同分割算法对比图，包括BSLS(B-SplineLevel-set，B样条水平集)、SBGFR(SelectiveBinaryandGaussianFilteringRegularizedLevelSet，选择性二元高斯滤波正规化水平集、LRCAM(Localizingregion-basedactivecontoursmodel，局部活动轮廓模型)、LGDF(localgaussiandistributionfittingengery，局部高斯分布拟合能量)、LBFM(localbinaryfittingmodel，局部二元拟合模型)、L2S(LegendreLevelSet，Legendre水平集)六种算法的对比效果。

本发明提供了一种频域边界能量模型的水平集图像分割方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。

Claims (5)

1.一种频域边界能量模型的水平集图像分割方法，其特征在于，包括以下步骤：

步骤1，在图像区域Ω范围内，对初始化轮廓的内区域和外区域进行傅里叶变换；

步骤2，利用滤波器提取图像的边界能量项；

步骤3，融合水平集函数正则化项、轮廓约束项及边界能量项，获得整体的能量项；

步骤4，实例化阶跃函数H(·)函数和冲击函数δ(·)，计算使得能量函数变小的水平集函数的变化值；对水平集函数进行更新，若能量函数不收敛，则返回步骤3；否则输出分割结果。

2.如权利要求1所述的一种频域边界能量模型的水平集图像分割方法，其特征在于，步骤1包括以下步骤：

步骤1-1，对于图像区域Ω上的输入图像u₀，输入图像u₀上每个像素对应一个Ω上的坐标向量x，在图像区域Ω上设置初始化轮廓C₀，设t时刻，在Ω上的轮廓为C_t，将C_t用水平集函数Φ表示，即C_t＝{x|Φ(x)＝0}，C_t将Ω划分为两个子区域，即C_t外区域Ω₁及内区域Ω₂，利用阶跃函数H(·)及水平集函数Φ(·)，将Ω₁及Ω₂分别表示为u₁和u₂，如下式所示：

u₁(x)＝u₀(x)H(φ(x)),u₂(x)＝u₀(x)(1-H(φ(x)))；

步骤1-2,在图像区域Ω范围内分别对u₁及u₂进行傅里叶变换：

其中，s为频域坐标向量，在Ω上对u_i进行傅里叶变换，F_i(s)为变换结果。

3.如权利要求2所述的一种频域边界能量模型的水平集图像分割方法，其特征在于，步骤2包括以下步骤：

步骤2-1，设置带通滤波器作为最优的边界滤波器，使用DoG算子作为带通滤波器，频域滤波器H_bp(s)的计算公式如下：

H_bp(s)＝DoG(s)＝G(s,σ_s1)-G(s,σ_s2)，

$G(s,{\mathrm{\&sigma;}}_{si})=\frac{1}{2{{\mathrm{\&pi;\&sigma;}}_{si}}^2}\exp (-\frac{s^{T}s}{2{\mathrm{\&sigma;}}_{si}^2}),i=1,2,$

G(s,σ_si)为变量为s标准差为σ_si的高斯滤波器，σ_s1和σ_s2分别为两个高斯滤波器的标准差，满足：σ_s1＞0，σ_s2＞0；

步骤2-2，使用频域滤波器H_bp(s)对F_i(s)进行滤波，并进行傅里叶逆变换得到目标的边界信息：

其中，f_b1(x)为图像区域上C_t外部的边界信息，f_b2(x)为图像区域上C_t内部的边界信息；为傅里叶逆变换；

步骤2-3，将C_t表示为水平集函数φ，计算边界能量项E^FBE(φ,f_bi(x))：

$E^{FBE}(\mathrm{\&phi;},f_{bi}(x\left)\right)=\underset{u_1,\mathrm{\&Omega;}}{\&Integral;}{f_{b1}}^2\left(x\right)H\left(\mathrm{\&phi;}\right)dx+\underset{u_2,\mathrm{\&Omega;}}{\&Integral;}{f_{b2}}^2\left(x\right)(1-H(\mathrm{\&phi;}\left)\right)dx,$

其中，φ为C_t对应的水平集函数；H(·)为阶跃函数。

4.如权利要求3所述的一种频域边界能量模型的水平集图像分割方法，其特征在于，步骤3包括以下步骤：

步骤3-1，采用以下公式计算水平集函数正则化项E^LSR：

$E^{LSR}={\&dtri;}^2\mathrm{\&phi;}-div\left(\frac{\&dtri;\mathrm{\&phi;}}{\left|\&dtri;\mathrm{\&phi;}\right|}\right),$

其中▽为梯度运算，div为散度运算；

步骤3-2，采用以下公式计算轮廓长度约束项E^Len：

$E^{Len}=\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)div\left(\frac{\&dtri;\mathrm{\&phi;}}{\left|\&dtri;\mathrm{\&phi;}\right|}\right),$

其中，δ(φ)为冲击函数；

步骤3-3，融合水平集函数正则项E^LSR、轮廓长度约束项E^Len以及频域边界能量项，采用如下公式计算整体的能量函数E(φ,f_b1(x))：

E(φ,f_b1(x))＝λE^FBE(φ,f_b1(x))+μE^LSR(φ)+νE^Len(φ)，

其中，λ，μ，ν分别为三个能量项的权值，满足：0<λ，0<μ，0<ν。

5.如权利要求4所述的一种频域边界能量模型的水平集图像分割方法，其特征在于，步骤4包括以下步骤：

步骤4-1：将阶跃函数和冲击函数分别实例化为H_ε(z)和δ_ε(z)，

$H_{\mathrm{\&epsiv;}}\left(z\right)=\frac{1}{2}\&lsqb;1+\frac{2}{\mathrm{\&pi;}}arctan\left(\frac{z}{\mathrm{\&epsiv;}}\right)\&rsqb;,$

${\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(z\right)=\frac{1}{\mathrm{\&pi;}}\frac{\mathrm{\&epsiv;}}{{\mathrm{\&epsiv;}}^2+z^2};$

其中，z与ε是表示上述函数的辅助变量；当ε趋于无限大时，H_ε(z)和δ_ε(z)两个函数都趋于理想的阶跃函数H(z)和冲击函数δ(z)；

步骤4-2，计算频域边界能量项下，水平集函数的更新方程，对水平集函数引入辅助变量α及辅助变量η，则：固定f_b1(x)、f_b2(x)，对Φ求偏导，使α趋近于0：

$\begin{array}{c}\frac{\&part;E^{FBE}}{\&part;\mathrm{\&phi;}}=\underset{\mathrm{\&alpha;}\&RightArrow;0}{\lim }\frac{d}{d\mathrm{\&alpha;}}{\left(\underset{\mathrm{\&Omega;}}{\&Integral;}{f}_{b1}\left(x\right)H\left(\widetilde{\mathrm{\&phi;}}\right)\right)}^2+{\left(f_{b2}\right(x\left)\right(1-H\left(\widetilde{\mathrm{\&phi;}}\right)\left)\right)}^{2}dx)\\ =\underset{\mathrm{\&alpha;}\&RightArrow;0}{\lim }(2\underset{\mathrm{\&Omega;}}{\&Integral;}\&lsqb;{f_{b1}}^2\left(x\right)H\left(\widetilde{\mathrm{\&phi;}}\right)-{f_{b2}}^2(x\left)\right(1-H\left)\right(\widetilde{\mathrm{\&phi;}}\left)\right)\&rsqb;\mathrm{\&delta;}\left(\widetilde{\mathrm{\&phi;}}\right)\mathrm{\&eta;}dx\\ =2\underset{\mathrm{\&Omega;}}{\&Integral;}\&lsqb;{f_{b1}}^2\left(x\right)H\left(\mathrm{\&phi;}\right)-{f_{b2}}^2\left(x\right)\left(1-H\right)\left(\mathrm{\&phi;}\right))\&rsqb;\mathrm{\&delta;}(\mathrm{\&phi;})\mathrm{\&eta;}dx\end{array}$

根据欧拉-拉格朗日公式，得到：

[f_b1 ²(x)H(φ)-f_b2 ²(x)(1-H(φ))]δ(φ)＝0，

运用梯度下降法则得到：

$\frac{\&part;\mathrm{\&phi;}}{\&part;t}=\&lsqb;{f_{b1}}^2\left(x\right)H\left(\mathrm{\&phi;}\right)-{f_{b2}}^2\left(x\right)(1-H(\mathrm{\&phi;}\left)\right)\&rsqb;\mathrm{\&delta;}\left(\mathrm{\&phi;}\right);$

步骤4-3，结合水平集函数正则项和长度约束项的水平集函数更新方式以及实际使用的阶跃函数和冲击函数，得到了水平集函数的最终更新值：

上式中为边界能量项； $\mathrm{\&mu;}({\&dtri;}^2\mathrm{\&phi;}-div(\frac{\&dtri;\mathrm{\&phi;}}{\left|\&dtri;\mathrm{\&phi;}\right|}\left)\right)$ 为水平函数正则化项；为轮廓长度约束项；

步骤4-4，利用上述最终更新值，对水平集函数进行更新，并判断水平集函数的变化是否小于特定值，若小于，则收敛，输出分割结果，该特定值取值范围为自然数；否则返回步骤4-3重新计算水平集函数的更新值，利用更新后水平集函数重新计算新的能量函数。