CN111145179B

CN111145179B - 一种基于水平集的灰度不均图像分割方法

Info

Publication number: CN111145179B
Application number: CN201911145061.6A
Authority: CN
Inventors: 房巾莉; 吕毅斌; 王樱子
Original assignee: Kunming University of Science and Technology
Current assignee: Kunming University of Science and Technology
Priority date: 2019-11-20
Filing date: 2019-11-20
Publication date: 2023-07-25
Anticipated expiration: 2039-11-20
Also published as: CN111145179A

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、在原始灰度不均图像中设置圆形初始演化曲线的位置，定义水平集函数的形式为符号距离函数；其中，图像各像素点的初始水平集函数值/>为：计算各像素点到初始演化曲线的最短距离；初始演化曲线上的各像素点到自身的最短距离为0，称作零水平集；n表示迭代次数；

S3、利用原始灰度不均图像的各像素值、各像素点的水平集函数值，计算图像各像素点的全局像素灰度均值c₁、c₂和局部区域拟合均值a₁、a₂；

S4、利用原始灰度不均图像的各像素值、各像素点的全局像素灰度均值和全局像素灰度均值，计算图像各像素点的全局拟合方差和局部拟合方差/>

S5、利用图像各像素的水平集函数值，计算长度项Length；计算惩罚项S6、利用速度停止算子、全局拟合方差、局部拟合方差、长度项和惩罚项，构造能量泛函，得到水平集函数迭代公式，计算迭代公式得到更新的水平集函数；

S7、检查是否满足迭代停止准则，若满足，则迭代停止，输出分割结果图像；若不满足，则返回步骤S3。

所述步骤S1中，各像素点的速度停止算子值，公式为：

式中，|▽G_σ*I(x)|是高斯平滑后图像G_σ*I(x)各像素点的梯度值，G_σ代表均值为0且方差为σ的高斯核，*代表卷积算子，▽代表梯度算子；e是自然底数；0＜q≤1。

所述步骤S2，具体为：在原始灰度不均图像中设置圆心为x₀、半径为r的圆形初始演化曲线C的位置，定义水平集函数为符号距离函数；其中图像各像素点的初始水平集函数值/>为：计算各像素点到初始演化曲线的最短距离；初始演化曲线上的各像素点到自身的最短距离为0，称作零水平集；

其中，演化曲线C与水平集函数的关系表达为：

式中，Ω表示图像区域；inside(C)表示演化曲线C的内部区域；outside(C)表示演化曲线C的外部区域；|dist(x,x₀)|表示任意像素点x到圆心像素x₀的距离。

所述步骤S3，具体为：利用原始灰度不均图像的各像素值、各像素点的水平集函数值，计算图像各像素点的全局算术均值的公式为：

局部拟合均值的公式为：

式中，c₁、c₂分别表示演化曲线内部和外部的全局像素灰度均值；a₁、a₂分别表示平均卷积算子后的图像演化曲线内部和外部的局部区域拟合均值；I(x)是原始灰度不均图像；K_A表示均值核；ε＝1为正规化参数，H_ε(·)为阶跃函数：则则/>

所述步骤S4，具体为：利用原始灰度不均图像的各像素值、各像素点的全局算术均值和局部拟合均值，计算图像各像素点的全局拟合方差的公式为：

局部拟合方差的公式为：

式中，分别表示像素点x在演化曲线内部和外部的全局拟合方差，分别表示像素点x在演化曲线内部和外部的局部拟合方差。

所述步骤S5中：

长度项Length的公式为：

惩罚项的公式为：

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

所述步骤S6中，构造的能量泛函为：

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

计算水平集函数迭代演化公式：得到更新的水平集函数值/>

式中，λ₁、λ₂是可变参数；div是散度算子，▽是梯度算子，Δ是Laplacian算子；n是迭代次数，Δt是步长。

所述步骤S7中，迭代停止准则包括：迭代次数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所示，一种基于水平集的灰度不均图像分割方法，所述方法步骤如下：

S5、利用图像各像素的水平集函数值，计算长度项Length来惩罚零水平集的曲线弧长；计算惩罚项使水平集函数始终近似符号距离函数，避免了水平集函数重新初始化的复杂过程；

S6、利用速度停止算子、全局拟合方差、局部拟合方差、长度项和惩罚项，构造能量泛函，得到水平集函数迭代公式，计算迭代公式得到更新的水平集函数；

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

式中，x代表图像的各像素点；|▽G_σ*I(x)|是高斯平滑后图像G_σ*I(x)各像素点的梯度值，G_σ代表均值为0且方差为σ的高斯核，*代表卷积算子，▽代表梯度算子；e是自然底数，e＝2.718…；0＜q≤1；0＜q≤1主要调节g(x)曲线在边界区域的下降快慢，防止弱边缘区域出现边界泄露；速度停止算子的目的是：利用高斯平滑去除噪声的同时保留图像边缘信息，避免演化陷入局部极小值，从而得到更光滑的演化曲线；调节轮廓曲线演化速度，并使轮廓曲线在目标边缘演化停止。

进一步地，可以设置所述步骤S2，具体为：在原始灰度不均图像中设置圆心为x₀、半径为r的圆形初始演化曲线C的位置，定义水平集函数为符号距离函数；其中图像各像素点的初始水平集函数值/>为：计算各像素点到初始演化曲线的最短距离；初始演化曲线上的各像素点到自身的最短距离为0，称作零水平集；

其中，演化曲线C与水平集函数的关系表达为：

或者可以描述为：

对于演化曲线C对应的任意像素x：

对于演化曲线C内部区域inside(C)对应的任意像素x：

对于演化曲线C外部区域outside(C)对应的任意像素x：式中，Ω表示图像区域；inside(C)表示演化曲线C的内部区域；outside(C)表示演化曲线C的外部区域；|dist(x,x₀)|表示任意像素点x到圆心像素x₀的距离。

进一步地，可以设置所述步骤S3，具体为：利用原始灰度不均图像的各像素值、各像素点的水平集函数值，计算图像各像素点的全局算术均值的公式为：

局部拟合均值的公式为：

式中，c₁、c₂分别表示演化曲线内部和外部的全局像素灰度均值，考虑了水平集函数内外的平均像素，因此该模型对噪声有一定的鲁棒性；a₁、a₂分别表示平均卷积算子后的图像演化曲线内部和外部的局部区域拟合均值；I(x)是原始灰度不均图像；K_A表示均值核；ε＝1为正规化参数，H_ε(·)为阶跃函数：则/> 则因此利用H_ε(·)可以方便地表示演化曲线的内部和外部。将局部与全局信息结合，一方面可以对窗口大小的选择不过于敏感，另一方面也可得到良好的演化效果。

进一步地，可以设置所述步骤S4，具体为：利用原始灰度不均图像的各像素值、各像素点的全局算术均值和局部拟合均值，计算图像各像素点的全局拟合方差的公式为：

局部拟合方差的公式为：

式中，分别表示像素点x在演化曲线内部和外部的全局拟合方差，分别表示像素点x在演化曲线内部和外部的局部拟合方差。对于灰度不均图像，将图像每个像素在其邻域内的局部拟合方差结合目标和背景的全局拟合方差后，可以更准确地计算出内部能量和外部能量，表达为轮廓曲线上的膨胀力和收缩力的大小，从而促使轮廓曲线逼近目标真正边缘。

进一步地，可以设置所述步骤S5中：

长度项Length的公式为：

惩罚项的公式为：

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

进一步地，所述步骤S6中，由于灰度不均效应在图像局部区域是缓慢变化的，甚至在小区域中的灰度是定值，因此我们结合全局拟合方差、局部拟合方差，利用速度停止算子、长度项和惩罚项，构造的能量泛函为：

计算水平集函数迭代演化公式：得到更新的水平集函数值/>

式中，λ₁、λ₂是可变参数；div是散度算子，▽是梯度算子，Δ是Laplacian算子；n是迭代次数，Δt是步长；对于灰度均匀的图像，λ₁和λ₂的数值接近或者相等，即全局项和局部项所起的作用应是一样的；对于灰度不均的图像，λ₁的数值应该小于λ₂的数值，从而使得局部项的效果得到进一步加强。驱动力项由全局拟合方差和局部拟合方差构成，这里将速度停止算子加权到驱动力项，除了可以调节演化曲线速度，还可以避免演化曲线陷入局部最优，使演化曲线更光滑；

进一步地，可以设置所述步骤S7中，迭代停止准则包括：准则1：迭代次数n达到设置的最大迭代次数Numiter；或者准则2：在规定的连续迭代次数Num内，演化曲线内部在迭代更新前后的面积差均小于给定阈值ω(迭代停止准则1或2满足一个即停止，输出结果)。可以针对不同灰度不均图像设置不同的Numiter、Num和ω的值。迭代过程中，检查是否满足迭代停止准则2，若满足，则迭代停止，输出分割结果图像；若不满足，则返回步骤S3，直到水平集函数迭代次数达到Numiter时停止迭代，输出分割结果图像。

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

待分割的灰度不均图像如图2、图7所示(实验1为灰度不均的合成图、实验2为灰度不均囊肿图)，则基于水平集的灰度不均图像分割方法包括：

S1、先计算高斯平滑后的图像G_σ*I(x)各像素点的梯度值|▽G_σ*I(x)|，再计算各像素点的速度停止算子：实验1、实验2均取σ＝0.8。

S2、在原始图像I(x)中选定圆心为x₀、半径为r的圆形初始演化曲线C的位置，计算初始水平集函数|dist(x,x₀)|表示任意像素点x到圆心像素x₀的距离。

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

S3、利用水平集函数，计算全局算术均值c₁、c₂，局部拟合均值a₁、a₂：

S4、利用c₁、c₂、a₁、a₂，计算图像各像素点的全局拟合方差局部拟合方差

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

S6、利用g(x)、以及长度项和惩罚项，计算水平集函数的迭代更新公式/>得到更新的图像各像素点的水平集函数值对于实验1，取λ₁＝0.1，λ₂＝7；对于实验2，取λ₁＝0.1，λ₂＝1。

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

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

Claims

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

S4、利用原始灰度不均图像的各像素值、各像素点的全局像素灰度均值和局部区域拟合均值，计算图像各像素点的全局拟合方差和局部拟合方差/>

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

S7、检查是否满足迭代停止准则，若满足，则迭代停止，输出分割结果图像；若不满足，则返回步骤S3；

所述步骤S6中，构造的能量泛函为：

计算水平集函数迭代演化公式：得到更新的水平集函数值/>式中，λ₁、λ₂是可变参数；div是散度算子，/>是梯度算子，Δ是Laplacian算子；n是迭代次数，Δt是步长；Ω表示图像区域；/>是高斯平滑后图像G_σ*I(x)各像素点的梯度值，G_σ代表均值为0且方差为σ的高斯核，*代表卷积算子；I(x)是原始灰度不均图像；μ为长度项的系数；ε为正规化参数，δ_ε(·)为狄克拉函数的近似表达。

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

式中，代表梯度算子；e是自然底数；0<q≤1。

3.根据权利要求1所述的基于水平集的灰度不均图像分割方法，其特征在于：所述步骤S2，具体为：在原始灰度不均图像中设置圆心为x₀、半径为r的圆形初始演化曲线C的位置，定义水平集函数为符号距离函数；其中图像各像素点的初始水平集函数值/>为：计算各像素点到初始演化曲线的最短距离；初始演化曲线上的各像素点到自身的最短距离为0，称作零水平集；

其中，演化曲线C与水平集函数的关系表达为：

式中，inside(C)表示演化曲线C的内部区域；outside(C)表示演化曲线C的外部区域；|dist(x,x₀)|表示任意像素点x到圆心像素x₀的距离。

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

局部区域拟合均值的公式为：

式中，c₁、c₂分别表示演化曲线内部和外部的全局像素灰度均值；a₁、a₂分别表示平均卷积算子后的图像演化曲线内部和外部的局部区域拟合均值；K_A表示均值核；ε＝1为正规化参数，H_ε(·)为阶跃函数：则/>则/>

5.根据权利要求4所述的基于水平集的灰度不均图像分割方法，其特征在于：所述步骤S4，具体为：利用原始灰度不均图像的各像素值、各像素点的全局像素灰度均值和局部区域拟合均值，计算图像各像素点的全局拟合方差的公式为：

局部拟合方差的公式为：

式中，分别表示像素点x在演化曲线内部和外部的全局拟合方差，/>分别表示像素点x在演化曲线内部和外部的局部拟合方差。

6.根据权利要求5所述的基于水平集的灰度不均图像分割方法，其特征在于：所述步骤S5中：

长度项Length的公式为：

惩罚项的公式为：

式中，取μ＝τ×255²，τ∈[0,1]；ε＝1；代表梯度算子。

7.根据权利要求1所述的基于水平集的灰度不均图像分割方法，其特征在于：所述步骤S7中，迭代停止准则包括：迭代次数n达到设置的最大迭代次数Numiter；或者在规定的连续迭代次数Num内，演化曲线内部在迭代更新前后的面积差均小于给定阈值ω。