CN106934815A - 基于混合区域的活动轮廓模型图像分割方法 - Google Patents

Abstract

本发明公开了一种基于混合区域的活动轮廓模型图像分割算法，该方法包括：通过图像增强算法对原始图像进行增强处理；同时使用原始图像的局部灰度拟合项和增强图像的全局灰度拟合项，构建基于混合区域的活动轮廓模型的能量泛函；通过变分法求解构建能量泛函，使用变分水平集的求解框架进行数学计算，即轮廓曲线以隐式水平集的形式表示，然后通过梯度下降流和欧拉一拉格朗日方程将构建的能量泛函转化为一个偏微分方程，并通过迭代逼近的方式求出所述偏微分方程的最优解，进而获得最终的目标边界轮廓。本发明对初始轮廓曲线的位置以及图像的灰度不均匀性不敏感，分割效率和精度都比较高，且受到噪声的影响较小，得到的分割结果鲁棒性较强。

Description

基于混合区域的活动轮廓模型图像分割方法

技术领域

本发明涉及一种图像处理领域的图像分割方法，具体涉及一种基于混合区域的活动轮廓模型图像分割方法。

技术背景

人类大脑从外界获取的各种信息中，绝大部分信息是通过眼睛获得的，物体或者场景发出的光或者反射的光通过眼睛中的晶状体折射到视网膜上，在视网膜上形成其倒立的图像，再由视觉感知神经传达给大脑，经大脑处理后即可看到该物体或场景的图像，而图像中包含其所表达物体或场景丰富的描述信息，因此图像信息是人类日常生活中获取信息的主要途径。

随着科学技术的飞速发展，特别是成像技术、计算机技术和信号处理理论的发展，人们将从各种成像器件中获得的图像信息，转变成数字信号，并在计算机中对其进行各种形式的加工和处理，如降噪、压缩、分割、增强、复原、特征抽取、分析和理解等，这些加工和处理即为人们常说的数字图像处理过程。图像分割是指按照图像中各个区域特有的性质，把图像分成几个区域，然后再把感兴趣的区域从这些区域中提取出来的技术。在数字图像处理的过程中，许多加工和处理都是建立在好的图像分割的基础上，只有先从图像中准确分割出目标物体，才能更好得进行下一步的图像处理及分析的工作。因此，图像分割是数字图像处理中一个非常基本、非常关键的技术。虽然图像分割技术很重要，而且到目前为止相关的科研人员也提出了种类多样的图像分割算法，但是由于成像系统的多样性和复杂性，造成了所成图像的多样性和复杂性，所以即使到目前为止科研人员已经提出了许多优秀的图像分割方法，但是还没有一个通用的图像分割算法，各种分割方法也只能用于针对某一类图像的分割问题，也因为这个原因，使得图像分割技术一直是数字图像处理技术当中研究的重点和难点。

传统的图像分割方法包括阈值法、边缘检测法以及基于区域提取的方法。而随着待解决问题复杂程度的提高，需要越来越多的高层知识融入到图像分割当中。此时出现了基于活动轮廓模型的图像分割方法，该方法依据轮廓曲线的表达形式的不同而分为参数活动轮廓模型和几何活动轮廓模型。而几何活动轮廓模型具有能有效处理曲线拓扑结构变化的优点，使得研究人员对其的研究越来越多。几何活动轮廓模型中，利用图像全局灰度信息的模型会对图像灰度不均匀性敏感；利用图像局部灰度信息的模型会对初始轮廓位置敏感；而结合图像全局和局部灰度信息的模型会存在一个全局项权重系数设置的问题。

发明内容

本发明的目的在于提供一种基于混合区域的活动轮廓模型图像分割方法。

实现本发明目的的技术方案为：一种基于混合区域的活动轮廓模型图像分割方法，包括以下步骤：

第一步，原始图像的增强处理：通过图像增强算法对原始图像进行增强处理；

第二步，构建能量泛函：同时使用原始图像的局部灰度拟合项和增强图像的全局灰度拟合项，构建基于混合区域的活动轮廓模型的能量泛函；

第三步，通过变分法求解构建能量泛函，使用变分水平集的求解框架进行数学计算，即轮廓曲线以隐式水平集的形式表示，然后通过梯度下降流和欧拉一拉格朗日方程将构建的能量泛函转化为一个偏微分方程，并通过迭代逼近的方式求出所述偏微分方程的最优解，进而获得最终的目标边界轮廓。

与现有技术相比，本发明的显著优点为：

(1)本发明对初始轮廓曲线的位置以及图像的灰度不均匀性不敏感，分割效率和精度都比较高，且受到噪声的影响较小，得到的分割结果鲁棒性较强；

(2)本发明能够克服现有技术的不足，将图像中感兴趣的目标从背景中清除、完整地分割出来，并且降低了算法的复杂程度。

附图说明

图1是本发明基于混合区域的活动轮廓模型的图像分割方法的流程图。

图2(a)和图2(b)分别是本发明仿真实验中原始图像和增强图像示意图。

图3(a)～图3(c)是本发明仿真实验中选取不同初始轮廓位置时的分割结果图。

图4(a)～图4(c)是本发明仿真试验中对待分割图像加入不同噪声后的分割结果图。

图5(a)～图5(c)是仿真实验中CV模型的分割结果示意图。

图6(a)～图6(c)是仿真实验中LRCV模型的分割结果示意图。

图7(a)～图7(c)是本发明仿真实验中分割结果示意图。

图8(a)～图8(c)是本发明仿真实验中与手工分割结果的对比图。

具体实施方式

结合图1，本发明的一种基于混合区域的活动轮廓模型图像分割方法，包括以下步骤：

第一步，原始图像的增强处理：通过图像增强算法对原始图像进行增强处理，即使用各项同性非抽取小波变换算法改善原始图像中存在的灰度不均匀现象，同时凸显图像中的目标边界；

第二步，构建能量泛函：同时使用原始图像的局部灰度拟合项和增强图像的全局灰度拟合项，构建基于混合区域的活动轮廓模型的能量泛函；

第三步：通过变分法求解构建能量泛函，使用变分水平集的求解框架进行数学计算，即轮廓曲线以隐式水平集的形式表示，然后通过梯度下降流和欧拉一拉格朗日方程将构建的能量泛函转化为一个偏微分方程，并通过迭代逼近的方式求出所述偏微分方程的最优解，进而获得最终的目标边界轮廓。

进一步的，所述的图像增强算法包括以下步骤：

S11，通过低通滤波器处理原始图像I₀，得到一个尺度系数a_j，其中a₀＝I₀；

S12，将相邻两个尺度系数a_j和a_j-1相减获得小波系数w_j；a_j＝h^j↑*a_j-1，h^j↑为低通滤波器；

S13，选择小波系数进行相加得到增强图像I_e＝w₁+w₂+w₃。

进一步的，第二步中构建的能量泛函为：

其中，λ₁、λ₂、η₁和η₂皆为非负常数，c_i(x)为图像局部灰度拟合项，是与图像大小一致的矢量，用于近似表示轮廓曲线内外局部区域的灰度；c_i为图像全局灰度拟合项，为轮廓曲线内外的灰度均值，R为水平集正则项，其定义如下：

其中μ和ν分别为非负常数，为梯度运算符；R的前半部分为距离惩罚项，可以保证φ在演化过程中与符号距离函数的一致性；后半部分为长度约束项，可以保证轮廓曲线C的光滑性。

由于轮廓曲线难以显示表示，因此在实际运算中常用水平集函数φ表示轮廓曲线，从而得到：

其中，M₁＝H(φ)和M₂＝1-H(φ)分别用于描述曲线C的内部和外部区域，H(φ)为正则化Heaviside函数，ε为大于零的常数。

进一步的，根据变分法将能量泛函转化为偏微分方程的形式为：

F₁＝-λ₁(I₀(x)-c₁(x))²+λ₂(I₀(x)-c₂(x))²

F₂＝-η₁(I_e(x)-c₁)²+η₂(I_e(x)-c₂)²

其中，div(·)表示曲率运算符，为Dirac函数；

根据变分法，得到常数c₁和c₂满足如下欧拉-拉格朗日方程：

∫(I_e(x)-c_i)M_idx＝0

经求解上式得到：

根据LRCV(Local region-based Chan-Vese)模型中局部灰度拟合项的求解得到：

式中，K_σ(x-y)为高斯核函数，x为图像中的某一点，y为x邻域范围内的某一点。

下面将结合附图，对本发明的实施例进行详细描述。

实施例

如图1所示，本发明的基于混合区域的活动轮廓模型的图像分割方法，通过以下实施步骤实现：

步骤1，对输入的待分割图像I₀进行增强处理：首先对图像I₀进行低通滤波，滤波器为上采样滤波算子h^j↑，是通过对h₀中相邻系数间插入2^j-1个零值获得，而h₀＝[1,4,6,4,1]/16则由三次B样条曲线衍生获得。经低通滤波获得各级尺度系数c_j，然后将相邻的两个尺度系数c_j和c_j-1相减可计算获得小波系数w_j，再将1-3级小波系数相加得到增强图像，即I_e＝w₁+w₂+w₃。

步骤2，基于混合区域的活动轮廓模型的能量泛函：

将原始图像I₀的局部灰度拟合项与增强图像I₀的全局灰度拟合项一起用于构建该能量泛函，其中局部灰度拟合项为选定像素点领域内外的灰度均值，而全局灰度拟合项为轮廓曲线内外的灰度均值。而只使用全局项和局部项构建能量泛函不足以获得满意的结果，故在实际应用中会加入R为水平集正则项用于保持轮廓曲线的平滑性、连续性以及保持轮廓的局部细节。

步骤3，变分求解框架：

对构建的完整的能量泛函求极小值，该过程使用变分法将其转化为偏微分方程的求解，即使用梯度下降流和欧拉一拉格朗日方程获得所需的一个关于时间的偏微分方程，并以迭代逼近的方式获取最终的目标边界。具体求解微分方程为:

F₁＝-λ₁(I₀(x)-c₁(x))²+λ₂(I₀(x)-c₂(x))²

F₂＝-η₁(I_e(x)-c₁)²+η₂(I_e(x)-c₂)²

步骤4，图像分割：

迭代逼近的求解偏微分方程需要设置初始的轮廓曲线，并设定曲线内表示目标所在的区域，曲线外围为背景区域。在微分方程的迭代计算下，轮廓曲线逐渐逼近所求目标的边界，从而将所需的目标分割出来。

下面介绍本发明的仿真实验：

(1)仿真条件：

本发明的仿真在Win7-64Intel(R)Core(TM)i3-2100 CPU@3.10GHz 3.10GHz RAM4GB平台上的MATLAB 2012a软件上进行仿真模拟的，仿真数据选用荷兰一个视网膜成像项目组公开发表视网膜影像数据(Digital Retinal Images for VesselExtraction,DRIVE)，该数据分为训练数据和测试数据两部分，每部分有20幅图像，图像大小为768×584；此外，该套数据还提供两组人工提取的血管作为实际所需的血管组织。

(2)仿真结果：

图2(a)和图2(b)显示原始图像存在很严重的灰度不均匀性，尤其是其中心位置处，背景信息与血管之间的灰度差异很小，导致二者难于分辨很困难，而增强图像中灰度不均匀性得到很明显的改善，血管和背景之间的灰度差异增大，用肉眼即可分辨出图像中的大部分血管。

图3(a)、图3(b)、图3(c)显示本发明在设置不同初始轮廓位置时的分割结果，图中方框表水初始轮廓位置，从图中可以看到分割结果基本稳定一致，均可以分割出视网膜图像中的大部分血管，而且对于图像中某些细小血管的分割也能取得较好的效果。

图4(a)、图4(b)、图4(c)显示分别给待分割图像加入不同强度的高斯白噪声时的本发明的分割结果，从图中可以看出，当噪声强度增加时，分割结果基本一致，效果都较好。

图5(a)～图5(c)，图6(a)～图6(c)，图7(a)～图7(c)，图8(a)～图8(c)分别显示CV模型、LRCV模型、本发明以及手工分割结果对比，从图中可以看出，本发明分割效果较好且准确率较高。

使用测试数据进行的仿真实验中，实验结果通过准确率ACC以及相似率DSC两个个指标进行量化。它们分别表示如下：

其中，tp表示正确定义的血管像素数，tn表示正确定义背景像素总数，fp表示错误定义的血管像素总数，fn表示错误背景像素总数，tp、tn、fp、fn的值皆为算法分割和手工分割比较后得到的数据。S₁和S₂分别表示算法分割和手工分割的结果，N()表示给定的集合中总的像素个数，∩表示交集符号。这两个量化标准的值越大，说明分割结果越接近手工分割结果，而且DSC值大于0.7便可视为优秀的分割结果。

对比实验的统计分析结果如下的表1所示：

表1

从表1可以看出：本文算法在定量对比时平均值最高，要优于其他两种算法。

基于本发明的公开的视网膜图像数据(Digital Retinal Images for VesselExtraction,DRIVE)的仿真实验表明：本发明能够从各种灰度不均匀的视网膜图像中较为准确完整的分割出其中的血管部分，且分割结果几乎不随初始轮廓位置变化而变化，同时有较强的抗噪声能力，因此解决上述的技术问题。

以上所述的仅是本发明的优选实施方式，应当指出，对于本领域的普通技术人员来说，在不脱离本发明创造构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。

Claims (4)

1.一种基于混合区域的活动轮廓模型图像分割方法，其特征在于，包括以下步骤：

第一步，原始图像的增强处理：通过图像增强算法对原始图像进行增强处理；

第二步，构建能量泛函：同时使用原始图像的局部灰度拟合项和增强图像的全局灰度拟合项，构建基于混合区域的活动轮廓模型的能量泛函；

第三步，通过变分法求解构建能量泛函，使用变分水平集的求解框架进行数学计算，即轮廓曲线以隐式水平集的形式表示，然后通过梯度下降流和欧拉一拉格朗日方程将构建的能量泛函转化为一个偏微分方程，并通过迭代逼近的方式求出所述偏微分方程的最优解，进而获得最终的目标边界轮廓。

2.根据权利要求1所述的基于混合区域的活动轮廓模型图像分割方法，其特征在于，所述的图像增强算法包括以下步骤：

S11，通过低通滤波器处理原始图像I₀，得到一个尺度系数a_j，其中a₀＝I₀；

S12，将相邻两个尺度系数a_j和a_j-1相减获得小波系数w_j，a_j＝h^j↑*a_j-1，h^j↑为低通滤波器；

S13，选择小波系数进行相加得到增强图像I_e＝w₁+w₂+w₃。

3.根据权利要求1所述的基于混合区域的活动轮廓模型图像分割方法，其特征在于，第二步中构建的能量泛函为：

$E=\underset{i=1}{\overset{2}{\Σ}}\∫\[{\mathrm{\λ}}_i{|I_0\left(x\right)-c_i\left(x\right)|}^2+{\mathrm{\η}}_i{|I_e\left(x\right)-c_i|}^2\]dx+R$

其中，λ₁、λ₂、η₁和η₂皆为非负常数，c_i(x)为图像局部灰度拟合项，是与图像大小一致的矢量，用于近似表示轮廓曲线内外局部区域的灰度；c_i为图像全局灰度拟合项，为轮廓曲线内外的灰度均值，R为水平集正则项，其定义如下：

$R=\mathrm{\μ}\∫\frac{1}{2}{\left(\right|\▿\mathrm{\φ}|-1)}^{2}dx+v\∫|\▿H\left(\mathrm{\φ}\right(x\left)\right)|dx$

其中μ和ν分别为非负常数，为梯度运算符。

采用水平集函数φ表示轮廓曲线，从而得到：

$E=\underset{i=1}{\overset{2}{\Σ}}\∫M_i\[{\mathrm{\λ}}_i{|I_0\left(x\right)-c_i\left(x\right)|}^2+{\mathrm{\η}}_i{|I_e\left(x\right)-c_i|}^2\]dx+R$

其中，M₁＝H(φ)和M₂＝1-H(φ)分别用于描述曲线C的内部和外部区域，H(φ)为正则化Heaviside函数，ε为大于零的常数。

4.根据权利要求1所述的基于混合区域的活动轮廓模型图像分割方法，其特征在于，根据变分法将能量泛函转化为偏微分方程的形式为：

$\frac{\∂\mathrm{\φ}}{\∂t}=\mathrm{\δ}\left(\mathrm{\φ}\right)\[F_1+F_2\]+\mathrm{\μ}\mathrm{\δ}\left(\mathrm{\φ}\right)div\left(\frac{\▿\mathrm{\φ}}{\left|\▿\mathrm{\φ}\right|}\right)+v(\▿{\mathrm{\φ}}^2-div(\frac{\▿\mathrm{\φ}}{\left|\▿\mathrm{\φ}\right|}\left)\right)$

F₁＝-λ₁(I₀(x)-c₁(x))²+λ₂(I₀(x)-c₂(x))²

F₂＝-η₁(I_e(x)-c₁)²+η₂(I_e(x)-c₂)²

其中，div(·)表示曲率运算符，为Dirac函数；

根据变分法，得到常数c₁和c₂满足如下欧拉-拉格朗日方程：

∫(I_e(x)-c_i)M_idx＝0

经求解上式得到：

$c_i=\frac{\∫I_e\left(x\right)M_{i}dx}{\∫M_{i}dx}$

根据LRCV模型中局部灰度拟合项的求解得到：

$c_i\left(x\right)=\frac{\∫K_{\mathrm{\σ}}(x-y)\[M_i{I}_0\left(x\right)\]dy}{\∫K_{\mathrm{\σ}}(x-y)M_{i}dy}$

式中，K_σ(x-y)为高斯核函数，x为图像中的某一点，y为x邻域范围内的某一点。