CN111145142A - 一种基于水平集算法的灰度不均囊肿图像分割方法 - Google Patents

Abstract

本发明公开了一种基于水平集算法的灰度不均囊肿图像分割方法，本发明在保留C‑V模型的全局方差信息的同时，在能量泛函中引入局部信息，使之与全局灰度均值叠加使算法对灰度不均图像的边缘具有全局和局域化效果，避免灰度不均效应使轮廓曲线的膨胀力和收缩力在非边缘处就制约从而导致分割失败。提出一个新的速度停止算子并引入驱动力项中，在迭代过程调节水平集曲线演化速度，避免演化陷入局部极小值，从而得到更光滑的演化曲线，通过在驱动力项引入速度停止算子，可以得到更光滑准确的分割结果。

Description

一种基于水平集算法的灰度不均囊肿图像分割方法

技术领域

本发明涉及一种基于水平集算法的灰度不均囊肿图像分割方法，属于图像处理领域。

背景技术

目前，数字图像处理技术广泛应用于工程学、计算机科学、统计学、物理、化学、医学、遥感等领域，而图像分割是图像处理过渡到图像分析的关键步骤，因此研究图像分割尤为重要。经过多年的研究，国内外研究者提出大量的图像分割方法，活动轮廓模型是一种有良好性能的分割算法，得到了广泛应用，也是计算机视觉和图像处理领域的研究热点，它包括以Snake模型为代表的参数活动轮廓模型和基于水平集的几何活动轮廓模型。

1988年Kass等人提出原始的参数化活动轮廓模型并应用于图像分割，但不易处理分割目标的拓扑结构，如单位法矢和曲率等曲线参数的计算比较费力。1988年Sethian和Osher等人首次提出依赖时间的运动界面的水平集描述，可以很自然地处理目标对象拓扑结构的改变，并构造了水平集方程的高精度稳定数值解法。1993、1995年Caselles等人将水平集理论和主动轮廓模型结合提出几何活动轮廓模型，也被称为水平集方法，就是将轮廓曲线在分割过程中作为零水平集被隐式地包含在水平集函数中，从而很自然地处理演化曲线的拓扑结构变化，并且由于水平集函数演化过程中始终保持为函数，很容易实现其数值近似算法。

早期的水平集方法主要是边缘型(edge-based)方法，依赖于图像的局部边缘信息，因此对于没有明显梯度变化或者梯度无意义的弱边缘的逼近效果不理想。对此，提出了区域型(region-based)方法，即利用区域信息来引导曲线向目标轮廓进行逼近，不依赖于梯度信息，可以分割出没有明显边缘的目标。其中，Chan和Vese在2001年提出的Chan-Vese(C-V)模型最具有代表性，对噪声有一定的鲁棒性，但仍然存在缺点：水平集对初始轮廓敏感，不能分割一些灰度不均图像，复杂的重新初始化数值步骤。为了解决这些问题，结合区域与边界信息的改进的水平集方法被提出，成为图像分割的研究热点。

而灰度不均图像的分割，还要求分割方法对灰度不均效应具有鲁棒性。图像中出现灰度不均现象有两方面原因，一是硬件干扰因素，如不均匀光照；二是成像物体本身因素，如物体的形状和位置。而医学图像由于存在局部体积效应、人体组织器官相互重叠和其成像过程带来的噪声等，灰度不均现象更加常见。虽然学者们已经提出了很多灰度不均匀校正算法，但实际上消除灰度不均匀效应至今仍然是一个难以解决的问题，因此直接在灰度不均图像中研究水平集演化以逼近真实边缘尤为重要。

发明内容

本发明提供了一种基于水平集算法的灰度不均囊肿图像分割方法。

本发明的技术方案是：一种基于水平集算法的灰度不均囊肿图像分割方法，所述方法步骤如下：

S1、读取原始灰度不均囊肿图像的各像素值I(x)，对图像进行高斯平滑，计算平滑后的图像各像素点的梯度值，并根据梯度值计算各像素点x的速度停止算子g(x)；

S2、在原始灰度不均囊肿图像中设置圆形的初始轮廓曲线C的位置，定义水平集函数

的形式为符号距离函数；其中图像各像素点的初始水平集函数值

为：各像素点到初始轮廓曲线的最短距离；初始轮廓曲线上的各像素点到自身的最短距离为0，称作零水平集；n表示迭代次数；

S3、利用原始灰度不均囊肿图像的各像素值、各像素点的水平集函数值，计算原始灰度不均囊肿图像各像素点的全局算术均值c₁、c₂；利用全局算术均值c₁、c₂，计算原始灰度不均囊肿图像各像素点的全局拟合方差

将速度停止算子g(x)加权到全局拟合方差

来构造并计算全局能量项

S4、利用原始灰度不均囊肿图像的各像素值、各像素点的水平集函数值，计算各像素点的局部拟合均值f₁(x)、f₂(x)；对局部拟合均值f₁(x)、f₂(x)进行加权，构造局部拟合图像，记F(x)为局部拟合图像各像素值；

S5、计算局部拟合图像各像素点的全局算术均值d₁、d₂；利用全局算术均值d₁、d₂，计算局部拟合图像各像素点的全局拟合方差

S6、计算局部拟合图像的各像素点的梯度值，并根据梯度值计算各像素点的速度停止算子V(x)，将速度停止算子V(x)加权到局部拟合图像的全局拟合方差

构造并计算局部能量项

S7、利用各像素点的水平集函数值，计算长度项Length；计算惩罚项

S8、利用

Length和

构造能量泛函，得出并计算水平集函数迭代更新公式，得到图像各像素点的更新的水平集函数值；

S9、迭代过程中检查是否满足迭代终止条件，若满足，则迭代停止，输出分割结果图像；若不满足，则返回步骤S3。

所述步骤S1中，速度停止算子值，公式为：

式中，

是高斯平滑后图像G_σ*I(x)各像素点的梯度值，G_σ代表均值为0且方差为σ的高斯核，*代表卷积算子，

代表梯度算子。

所述步骤S2，具体为：在原始灰度不均囊肿图像中设置圆心像素为x₀、半径为r的圆形初始轮廓曲线C的位置，定义水平集函数

为符号距离函数；其中图像各像素点的初始水平集函数值

为：计算各像素点到初始轮廓曲线的最短距离；初始轮廓曲线上的各像素点到自身的最短距离为0，称作零水平集；

其中，轮廓曲线C与水平集函数

的关系表达为：

式中，Ω表示图像区域；inside(C)表示轮廓曲线C的内部区域；outside(C)表示轮廓曲线C的外部区域；

|dist(x,x₀)|表示任意像素点x到圆心像素x₀的距离。

所述步骤S3中，利用原始灰度不均囊肿图像的各像素值、各像素点的水平集函数值，计算原始灰度不均囊肿图像各像素点的全局算术均值c₁、c₂：

式中，c₁、c₂分别表示原始灰度不均囊肿图像在轮廓曲线C内部和外部的像素点的全局算术均值，ε＝1为正规化参数，H_ε(·)为阶跃函数，表示为：

利用全局算术均值c₁、c₂，计算I(x)各像素点的全局拟合方差：

式中，

分别表示原始灰度不均囊肿图像在轮廓曲线C内部和外部的像素点的全局拟合方差；

将速度停止算子g(x)加权到I(x)的拟合方差

来构造全局能量项：

所述步骤S4中，利用原始灰度不均囊肿图像的各像素值、各像素点的水平集函数值，计算各像素点的局部拟合均值f₁(x)、f₂(x)：

式中，

表示大小为m₁×m₁的均值核，

表示大小为m₂×m₂的均值核；f₁(x)、f₂(x)分别表示原始灰度不均囊肿图像在轮廓曲线内部和外部的像素点的局部拟合均值；

利用局部拟合均值f₁(x)、f₂(x)构造局部拟合图像，记F(x)为局部拟合图像的各像素值：

所述步骤S5中，计算局部拟合图像各像素点的全局算术均值d₁、d₂：

式中，d₁、d₂分别表示局部拟合图像在轮廓曲线C内部和外部的像素点的全局算术均值，再计算局部拟合图像各像素点的全局拟合方差：

式中，

分别表示局部拟合图像在轮廓曲线C内部和外部的像素点的全局拟合方差。

所述步骤S6中，利用局部拟合图像，计算各像素点的梯度值

并根据梯度值计算各像素点的速度停止算子V(x)：

将速度停止算子V(x)加权到局部拟合图像的拟合方差

来构造并计算局部能量项：

所述步骤S7中，利用各像素点的水平集函数值，计算长度项Length：

式中，μ为长度项的系数，取μ＝τ×255²，τ∈[0,1]；ε＝1为正规化参数，δ_ε(·)为狄克拉函数的近似表达，

代表梯度算子；

计算惩罚项

所述步骤S8中，利用全局能量项、局部能量项、长度项和惩罚项，构造的能量泛函为：

对能量泛函利用变分法的梯度下降流，得到水平集函数的迭代演化的Euler-Lagrange方程：

计算水平集函数更新演化的迭代计算公式为：

得到更新的水平集函数值

式中，λ₁、λ₂是可变参数；div是散度算子，

是梯度算子，Δ是Laplacian算子；n是迭代次数，Δt是步长。

所述步骤S9中，迭代终止条件为：迭代次数n达到设置的最大迭代次数Numiter；或者在规定的连续迭代次数Num内，均满足：

其中ω表示阈值。

本发明的有益效果是：本发明在保留C-V模型的全局方差信息的同时，在能量泛函中引入局部信息，使之与全局灰度均值叠加使算法对灰度不均图像的边缘具有全局和局域化效果，避免灰度不均效应使轮廓曲线的膨胀力和收缩力在非边缘处就制约从而导致分割失败。提出一个新的速度停止算子并引入驱动力项中，在迭代过程调节水平集曲线演化速度，避免演化陷入局部极小值，从而得到更光滑的演化曲线，通过在驱动力项引入速度停止算子，可以得到更光滑准确的分割结果。

附图说明

图1为本发明的流程图；

图2为实验1待分割的灰度不均图像的示意图；

图3为实验1待分割的灰度不均图像及其初始演化曲线的示意图；

图4为实验1灰度不均图像的分割结果及其最终演化曲线的示意图；

图5为实验1初始三维水平集函数的示意图；

图6为实验1最终的三维水平集函数的示意图；

图7为实验2待分割的灰度不均图像的示意图；

图8为实验2待分割的灰度不均图像及其初始演化曲线的示意图；

图9为实验2灰度不均图像的分割结果及其最终演化曲线的示意图；

图10为实验2初始三维水平集函数的示意图；

图11为实验2最终的三维水平集函数的示意图。

具体实施方式

实施例1：如图1-11所示，一种基于水平集算法的灰度不均囊肿图像分割方法，所述方法步骤如下：

S1、读取原始灰度不均囊肿图像的各像素值I(x)，对图像进行高斯平滑，计算平滑后的图像各像素点的梯度值，并根据梯度值计算各像素点x的速度停止算子g(x)；

S2、在原始灰度不均囊肿图像中设置圆形的初始轮廓曲线C的位置，定义水平集函数

的形式为符号距离函数；其中图像各像素点的初始水平集函数值

为：各像素点到初始轮廓曲线的最短距离；初始轮廓曲线上的各像素点到自身的最短距离为0，称作零水平集；n表示迭代次数；

S3、利用原始灰度不均囊肿图像的各像素值、各像素点的水平集函数值，计算原始灰度不均囊肿图像各像素点的全局算术均值c₁、c₂；利用全局算术均值c₁、c₂，计算原始灰度不均囊肿图像各像素点的全局拟合方差

将速度停止算子g(x)加权到全局拟合方差

来构造并计算全局能量项

S4、利用原始灰度不均囊肿图像的各像素值、各像素点的水平集函数值，计算各像素点的局部拟合均值f₁(x)、f₂(x)；对局部拟合均值f₁(x)、f₂(x)进行加权，构造局部拟合图像，记F(x)为局部拟合图像各像素值；

S5、计算局部拟合图像各像素点的全局算术均值d₁、d₂；利用全局算术均值d₁、d₂，计算局部拟合图像各像素点的全局拟合方差

S6、计算局部拟合图像的各像素点的梯度值，并根据梯度值计算各像素点的速度停止算子V(x)，将速度停止算子V(x)加权到局部拟合图像的全局拟合方差

构造并计算局部能量项

S7、利用各像素点的水平集函数值，计算长度项Length来惩罚轮廓曲线的弧长；计算惩罚项

使水平集函数始终近似符号距离函数，避免重新初始化的复杂过程；

S8、利用

Length和

构造能量泛函，得出并计算水平集函数迭代更新公式，得到图像各像素点的更新的水平集函数值；

S9、迭代过程中检查是否满足迭代终止条件，若满足，则迭代停止，输出分割结果图像；若不满足，则返回步骤S3。

进一步地，可以设置所述步骤S1中，速度停止算子值，公式为：

式中，

是高斯平滑后图像G_σ*I(x)各像素点的梯度值，G_σ代表均值为0且方差为σ的高斯核，*代表卷积算子，

代表梯度算子；本算子利用高斯平滑去除噪声的同时保留图像边缘信息，避免演化陷入局部极小值，从而得到更光滑的轮廓曲线；调节轮廓曲线演化速度，并使轮廓曲线在目标边缘演化停止；

进一步地，可以设置所述步骤S2，具体为：在原始灰度不均囊肿图像中设置圆心像素为x₀、半径为r的圆形初始轮廓曲线C的位置，定义水平集函数

为符号距离函数；其中图像各像素点的初始水平集函数值

为：计算各像素点到初始轮廓曲线的最短距离；初始轮廓曲线上的各像素点到自身的最短距离为0，称作零水平集；

其中，轮廓曲线C与水平集函数

的关系表达为：

式中，Ω表示图像区域；inside(C)表示轮廓曲线C的内部区域；outside(C)表示轮廓曲线C的外部区域；

|dist(x,x₀)|表示任意像素点x到圆心像素x₀的距离。

可以看出，曲线C上的各像素点到自身的最短距离为0，因此曲线C被称作零水平集；曲线C内部像素的初始水平集函数值为正，曲线C外部像素的初始水平集函数值为负。

进一步地，可以设置所述步骤S3中，利用原始灰度不均囊肿图像的各像素值、各像素点的水平集函数值，计算原始灰度不均囊肿图像各像素点的全局算术均值c₁、c2：

式中，c₁、c₂分别表示原始灰度不均囊肿图像在轮廓曲线C内部和外部的像素点的全局算术均值，ε＝1为正规化参数，H_ε(·)为阶跃函数，表示为：

因此利用H_ε(·)可以方便地表示曲线C的内部和外部。

利用全局算术均值c₁、c₂，计算I(x)各像素点的全局拟合方差：

式中，

分别表示原始灰度不均囊肿图像在轮廓曲线C内部和外部的像素点的全局拟合方差；

将速度停止算子g(x)加权到I(x)的拟合方差

来构造全局能量项：

进一步地，可以设置所述步骤S4中，利用原始灰度不均囊肿图像的各像素值、各像素点的水平集函数值，计算各像素点的局部拟合均值f₁(x)、f₂(x)：

式中，

表示大小为m₁×m₁的均值核，

表示大小为m₂×m₂的均值核；f₁(x)、f₂(x)分别表示原始灰度不均囊肿图像在轮廓曲线内部和外部的像素点的局部拟合均值；可以针对囊肿图像的灰度不均的不同程度的设置不同的m₁和m₂的值。

利用局部拟合均值f₁(x)、f₂(x)构造局部拟合图像，记F(x)为局部拟合图像的各像素值：

进一步地，可以设置所述步骤S5中，计算局部拟合图像各像素点的全局算术均值d₁、d₂：

式中，d₁、d₂分别表示局部拟合图像在轮廓曲线C内部和外部的像素点的全局算术均值，再计算局部拟合图像各像素点的全局拟合方差：

式中，

分别表示局部拟合图像在轮廓曲线C内部和外部的像素点的全局拟合方差。

进一步地，可以设置所述步骤S6中，利用局部拟合图像，计算各像素点的梯度值

并根据梯度值计算各像素点的速度停止算子V(x)：

将速度停止算子V(x)加权到局部拟合图像的拟合方差

来构造并计算局部能量项：

进一步地，可以设置所述步骤S7中，利用各像素点的水平集函数值，计算长度项Length：

式中，μ为长度项的系数，取μ＝τ×255²，τ∈[0,1]；一般在分割小目标时，τ取较小的数值，分割大目标时，τ取较大的数值；ε＝1为正规化参数，δ_ε(·)为狄克拉函数的近似表达，

代表梯度算子。

计算惩罚项

进一步地，可以设置所述步骤S8中，利用全局能量项、局部能量项、长度项和惩罚项，构造的能量泛函为：

对于灰度较均匀的囊肿图像，λ₁和λ₂的数值接近或者相等，即全局能量项和局部能量项所起的作用应是一样的；对于灰度不均的囊肿图像，λ₁的数值应该小于λ₂的数值，从而使得局部能量项的效果得到进一步加强。驱动力项由全局能量项和局部能量项构成，这里将两种速度停止算子分别加权到驱动力项，除了可以调节轮廓曲线速度，还可以避免轮廓曲线陷入局部最优，使轮廓曲线更光滑。

对能量泛函利用变分法的梯度下降流，得到水平集函数的迭代演化的Euler-Lagrange方程：

计算水平集函数更新演化的迭代计算公式为：

得到更新的水平集函数值

式中，λ₁、λ₂是可变参数；div是散度算子，

是梯度算子，Δ是Laplacian算子；n是迭代次数，Δt是步长。

进一步地，可以设置所述步骤S9中，迭代终止条件为：迭代次数n达到设置的最大迭代次数Numiter；或者在规定的连续迭代次数Num内，均满足：

其中ω表示阈值。

可以针对囊肿图像的灰度不均的不同程度的设置不同的Numiter、Num和ω的值。迭代时，检查是否满足迭代终止条件，若满足，则迭代停止，输出分割结果图像；若不满足，则返回步骤S3。

下面通过代入具体的灰度不均囊肿图像对实施例1作进一步的说明：

待分割的原始灰度不均囊肿图像如图2、图7所示(实验1和实验2分别为两幅不同的子宫囊肿图)，则基于水平集算法的灰度不均囊肿图像分割方法包括：

S1、读入原始的灰度不均囊肿图像。计算高斯平滑后图像G_σ*I(x)的各像素点的梯度值

计算各像素点的速度停止算子：

实验1、实验2均取高斯核窗口大小为15，σ＝0.8。

S2、在原始图像中设置圆心像素为x₀、半径为r的圆形初始轮廓曲线C，定义水平集函数为符号距离函数；计算图像各像素点的初始水平集函数值

即计算各像素点到曲线C的最短距离为：

|dist(x,x₀)|表示任意像素点x到圆心像素x₀的距离。

二维灰度不均囊肿图像及其初始演化轮廓曲线如图3、图8所示，三维初始水平集函数

如图5、图10所示。

S3、利用原始图像各像素值、水平集函数值，计算图像在曲线C内部和外部的像素全局灰度算术均值

利用c₁、c₂，计算原始图像在轮廓曲线C内部和外部的像素点的全局拟合方差：

将速度停止算子g(x)加权到拟合方差

来构造全局能量项：

S4、利用原始图像的各像素值I(x)、各像素点的水平集函数值

计算平均卷积算子后的图像像素点x的两种局部拟合均值f₁(x)、f₂(x)：

表示大小为m₁×m₁的均值核，

表示大小为m₂×m₂的均值核。利用f₁(x)、f₂(x)构造局部拟合图像，记F(x)为局部拟合图像的各像素值：

对于实验1，取m₁＝5，m₂＝7；对于实验2，取m₁＝m₂＝5。

S5、计算局部拟合图像F(x)在轮廓曲线C内部和外部的各像素点的全局灰度均值：

利用d₁、d₂计算图像F(x)在轮廓曲线C内部和外部的各像素点的全局拟合方差:

S6、利用局部拟合图像，计算各像素点的梯度值

并根据梯度值计算各像素点的速度停止算子:

将速度停止算子V(x)加权到局部拟合图像的拟合方差

计算局部能量项：

S7、利用水平集函数，计算两个规则化项：长度项和惩罚项。实验1、实验2取μ＝0.01×255²。

S8、结合全局能量项、局部能量项，采用长度项和惩罚项来规则化能量项，迭代计算水平集函数更新演化的公式：

得到更新的水平集函数值

对于实验1、实验2，取步长Δt＝0.1，λ₁＝λ₂＝1。

S9、迭代过程中判断是否满足迭代终止条件，若满足，则结束迭代，输出分割结果；不满足，则返回步骤S3继续迭代。最终得到的轮廓曲线如图4、图9所示，最终的三维水平集函数如图6、图11所示。对于实验1，取Numiter＝165，Num＝15，ω＝5；对于实验2，取Numiter＝50，Num＝15，ω＝5。

上面结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。

Claims (10)

1.一种基于水平集算法的灰度不均囊肿图像分割方法，其特征在于：所述方法步骤如下：

S1、读取原始灰度不均囊肿图像的各像素值I(x)，对图像进行高斯平滑，计算平滑后的图像各像素点的梯度值，并根据梯度值计算各像素点x的速度停止算子g(x)；

S2、在原始灰度不均囊肿图像中设置圆形的初始轮廓曲线C的位置，定义水平集函数

的形式为符号距离函数；其中图像各像素点的初始水平集函数值

为：各像素点到初始轮廓曲线的最短距离；初始轮廓曲线上的各像素点到自身的最短距离为0，称作零水平集；n表示迭代次数；

S3、利用原始灰度不均囊肿图像的各像素值、各像素点的水平集函数值，计算原始灰度不均囊肿图像各像素点的全局算术均值c₁、c₂；利用全局算术均值c₁、c₂，计算原始灰度不均囊肿图像各像素点的全局拟合方差

将速度停止算子g(x)加权到全局拟合方差

来构造并计算全局能量项

S4、利用原始灰度不均囊肿图像的各像素值、各像素点的水平集函数值，计算各像素点的局部拟合均值f₁(x)、f₂(x)；对局部拟合均值f₁(x)、f₂(x)进行加权，构造局部拟合图像，记F(x)为局部拟合图像各像素值；

S5、计算局部拟合图像各像素点的全局算术均值d₁、d₂；利用全局算术均值d₁、d₂，计算局部拟合图像各像素点的全局拟合方差

S6、计算局部拟合图像的各像素点的梯度值，并根据梯度值计算各像素点的速度停止算子V(x)，将速度停止算子V(x)加权到局部拟合图像的全局拟合方差

构造并计算局部能量项

S7、利用各像素点的水平集函数值，计算长度项Length；计算惩罚项

S8、利用

Length和

构造能量泛函，得出并计算水平集函数迭代更新公式，得到图像各像素点的更新的水平集函数值；

S9、迭代过程中检查是否满足迭代终止条件，若满足，则迭代停止，输出分割结果图像；若不满足，则返回步骤S3。

2.根据权利要求1所述的基于水平集算法的灰度不均囊肿图像分割方法，其特征在于：所述步骤S1中，速度停止算子值，公式为：

式中，

是高斯平滑后图像G_σ*I(x)各像素点的梯度值，G_σ代表均值为0且方差为σ的高斯核，*代表卷积算子，

代表梯度算子。

3.根据权利要求1所述的基于水平集算法的灰度不均囊肿图像分割方法，其特征在于：所述步骤S2，具体为：在原始灰度不均囊肿图像中设置圆心像素为x₀、半径为r的圆形初始轮廓曲线C的位置，定义水平集函数

为符号距离函数；其中图像各像素点的初始水平集函数值

为：计算各像素点到初始轮廓曲线的最短距离；初始轮廓曲线上的各像素点到自身的最短距离为0，称作零水平集；

其中，轮廓曲线C与水平集函数

的关系表达为：

式中，Ω表示图像区域；inside(C)表示轮廓曲线C的内部区域；outside(C)表示轮廓曲线C的外部区域；

|dist(x,x₀)|表示任意像素点x到圆心像素x₀的距离。

4.根据权利要求3所述的基于水平集算法的灰度不均囊肿图像分割方法，其特征在于：所述步骤S3中，利用原始灰度不均囊肿图像的各像素值、各像素点的水平集函数值，计算原始灰度不均囊肿图像各像素点的全局算术均值c₁、c₂：

式中，c₁、c₂分别表示原始灰度不均囊肿图像在轮廓曲线C内部和外部的像素点的全局算术均值，ε＝1为正规化参数，H_ε(·)为阶跃函数，表示为：

利用全局算术均值c₁、c₂，计算I(x)各像素点的全局拟合方差：

式中，

分别表示原始灰度不均囊肿图像在轮廓曲线C内部和外部的像素点的全局拟合方差；

将速度停止算子g(x)加权到I(x)的拟合方差

来构造全局能量项：

5.根据权利要求4所述的基于水平集算法的灰度不均囊肿图像分割方法，其特征在于：所述步骤S4中，利用原始灰度不均囊肿图像的各像素值、各像素点的水平集函数值，计算各像素点的局部拟合均值f₁(x)、f₂(x)：

式中，

表示大小为m₁×m₁的均值核，

表示大小为m₂×m₂的均值核；f₁(x)、f₂(x)分别表示原始灰度不均囊肿图像在轮廓曲线内部和外部的像素点的局部拟合均值；

利用局部拟合均值f₁(x)、f₂(x)构造局部拟合图像，记F(x)为局部拟合图像的各像素值：

6.根据权利要求5所述的基于水平集算法的灰度不均囊肿图像分割方法，其特征在于：所述步骤S5中，计算局部拟合图像各像素点的全局算术均值d₁、d₂：

式中，d₁、d₂分别表示局部拟合图像在轮廓曲线C内部和外部的像素点的全局算术均值，再计算局部拟合图像各像素点的全局拟合方差：

式中，

分别表示局部拟合图像在轮廓曲线C内部和外部的像素点的全局拟合方差。

7.根据权利要求6所述的基于水平集算法的灰度不均囊肿图像分割方法，其特征在于：所述步骤S6中，利用局部拟合图像，计算各像素点的梯度值

并根据梯度值计算各像素点的速度停止算子V(x)：

将速度停止算子V(x)加权到局部拟合图像的拟合方差

来构造并计算局部能量项：

8.根据权利要求7所述的基于水平集算法的灰度不均囊肿图像分割方法，其特征在于：所述步骤S7中，利用各像素点的水平集函数值，计算长度项Length：

式中，μ为长度项的系数，取μ＝τ×255²，τ∈[0,1]；ε＝1为正规化参数，δ_ε(·)为狄克拉函数的近似表达，

代表梯度算子；

计算惩罚项

9.根据权利要求8所述的基于水平集算法的灰度不均囊肿图像分割方法，其特征在于：所述步骤S8中，利用全局能量项、局部能量项、长度项和惩罚项，构造的能量泛函为：

对能量泛函利用变分法的梯度下降流，得到水平集函数的迭代演化的Euler-Lagrange方程：

计算水平集函数更新演化的迭代计算公式为：

得到更新的水平集函数值

式中，λ₁、λ₂是可变参数；div是散度算子，

是梯度算子，Δ是Laplacian算子；n是迭代次数，Δt是步长。

10.根据权利要求9所述的基于水平集算法的灰度不均囊肿图像分割方法，其特征在于：所述步骤S9中，迭代终止条件为：迭代次数n达到设置的最大迭代次数Numiter；或者在规定的连续迭代次数Num内，均满足：

其中ω表示阈值。

Images