0 引言

图像分割是图像处理过程中关键的一部分，目的是将目标对象与背景分离，为下一步的识别、配准等操作奠定基础，所以一个好的分割起着至关重要的作用。活动轮廓模型是一种基于水平集和曲线演变理论的有效的分割方法，是目前应用广泛且具有挑战性的研究。活动轮廓模型主要分为基于边缘^[1-3]的和基于区域^[4-6]的模型：前者依赖于目标对象的边缘信息^[7]，对弱边界目标处理较为困难；而后者依赖于图像区域的统计信息，对图像灰度分布的均匀性^[5-6]较为敏感。

测地活动轮廓(Geodesic Active Contour, GAC)模型^[1]是一种经典的基于边缘的活动轮廓模型，具有实施简单快速等特点，在图像分割领域有着较好的分割结果。但同时也存在许多缺点，其主要问题包括能量函数依赖图像边缘信息，对于弱边界图像或边界梯度无定义的图像处理困难，边缘停止函数的好坏直接影响最终的分割效果；水平集函数在迭代过程中容易产生震荡，偏离符号距离函数特性，需要进行定期重新初始化，导致演化过程中迭代次数较多、耗时较长等。当前基于边缘的活动轮廓模型主要的改进方法是加入一些先验信息^[8]，或与基于区域的模型相结合^[9-10]等。而对于水平集重新初始化问题，当前较为广泛的是距离正则化水平集演变(Distance Regularized Level Set Evolution, DRLSE)方法^[11-12]。近年文献[13]提出了一种无重新初始化模型——反应扩散(Reaction Diffusion, RD)模型，解决了DRLSE模型对弱边界效果差、对噪声敏感等问题。

本文提出了一种基于贝塞尔滤波的边缘停止函数，充分利用贝塞尔边缘检测滤波的优良性能来构建新的边缘停止函数，达到了较好的分割效果，同时采用文献[13]提出反应扩散方程，对水平集函数进行正则化处理，提出的模型不仅在时间上提高了效率，也达到了更好的分割效果。

1 研究背景 1.1 测地活动轮廓模型

活动轮廓模型的基本思想是构建关于曲线的能量泛函，把图像分割的过程转化为能量泛函最小化的过程。基于边缘的活动轮廓模型能得到更光滑的边界而被广泛采用，其主要解决方案是：首先通过边缘检测函数检测图像的边缘信息，通常这些信息是离散的，然后利用曲线演变理论，使轮廓线停止在图像的边缘，从而将稀疏的、不连续的边缘信息转化为一条闭合曲线。其能量泛函主要包含两个部分:一是控制曲线运动的力项，用来缩小或扩大活动轮廓；另一个是基于边缘的停止项，使曲线停止在所需的对象边界。停止项通常设计为基于图像梯度信息的边缘检测器，边缘函数的好坏直接影响活动轮廓演变的最终结果。

测地活动轮廓模型是一种应用广泛的基于边缘的活动轮廓模型，其能量泛函定义如下：

$ {{L}_{R}}(\mathit{\boldsymbol{C}})=\int_{0}^{L(\mathit{\boldsymbol{C}})}{g(|\nabla I[\mathit{\boldsymbol{C}}(s)]|)}\rm{d}\mathit{s} $

(1)

其中:C为闭合曲线，L(C)为C的弧长，L_R(C)为加权弧长; ▽I为图像的梯度值; g(·)为边缘检测器。通常设置为：

$ g=1/\left( 1+|\nabla {{G}_{\sigma }}*I{{|}^{2}} \right) $

(2)

其中：G_σ为高斯滤波，与图像卷积使图像光滑降低噪声。在图像灰度值均匀区域(即平坦区域)|▽I|≈0，于是g≈1；在图像灰度值急剧变化的区域(即边缘区域)|▽I|取得极大值，g取得极小值，同时能量函数达到最小值。于是求图像边缘轮廓的过程转变为求该能量泛函极小值的过程，通过梯度下降方法可以得到最小化能量方程：

$ \frac{\partial \mathit{\boldsymbol{C}}}{\partial t}=g(\mathit{\boldsymbol{C}})\kappa N-(\nabla g\cdot N)N $

(3)

其中：κ为曲线C的曲率；N为法线。法线方向总指向曲线内部。

在实际操作过程中，当图像中目标对象存在有较深的凹陷区域时，能量函数会陷入局部最小值，而提前停止演变，不能到达真正的边界。一种简单的解决方案是加入一个受到边缘停止函数控制的收缩力，力的方向总是指向轮廓线的内部，于是GAC模型的方程为：

$ \frac{\partial \mathit{\boldsymbol{C}}}{\partial t}=g(c+\kappa )N-(\nabla g\cdot N)N $

(4)

式中:c一般为常数，当曲线的曲率κ＜0时，在c的作用下曲线仍能够继续演变，达到真正的边界。如果收缩力过大，曲线可能会穿过目标边界；力过小又可能提前停止演变而达不到目标边界：所以需要在具体操作中选择合适的力度大小。

1.2 水平集方法

当前基于水平集的活动轮廓模型应用广泛，由于其可以对演化中的曲线直接进行数值计算而不必参数化，并且可以有效地处理拓扑变化。在水平集演变过程中，由于数值稳定性问题，需要进行正则化处理。传统解决办法是需要定期地将水平集重新初始化为符号距离函数以保持曲线演变的稳定，但这一过程时间复杂度较高。

当前许多模型都提出了无重新初始化方法来避免这一过程。应用最广泛的是距离正则化^[11-12]方法，其基本思想在能量函数中加入符号距离惩罚能量函数：

$ P(\mathit{\boldsymbol{\phi}} )=\int{\frac{1}{2}}{{(|\nabla {\mathit{\boldsymbol{\phi}}} (x)-1|)}^{2}}\rm{d}\mathit{x} $

(5)

其中:ϕ为水平集函数，▽为梯度算子。通过梯度下降法最小化该能量函数可得热扩散方程：

$ \frac{\partial \mathit{\boldsymbol{C}}}{\partial t}=\rm{div}[\mathit{r}(\mathit{\boldsymbol{\phi}} )\nabla \mathit{\boldsymbol{\phi}}] $

(6)

其中r为扩散率，且在文献[11]中有：

$ r({\mathit{\boldsymbol{\phi}}} ) = 1 - 1/|\nabla {\mathit{\boldsymbol{\phi}}} | $

(7)

由于当▽ϕ趋近于1时，r(ϕ)趋近于0，所以会导致水平集函数的震荡，于是文献[12]提出了改进的的惩罚函数：

$ P({\mathit{\boldsymbol{\phi}}} )=\left\{ \begin{align} & \frac{1}{{{(2\rm{ \mathsf{ π} })}^{2}}}(1-\cos (2\rm{ \mathsf{ π}}|\nabla {\mathit{\boldsymbol{\phi}}} |)), \;\;|\nabla {\mathit{\boldsymbol{\phi}}} |\le 1 \\ & \frac{1}{2}{{(|\nabla {\mathit{\boldsymbol{\phi}}} |-1)}^{2}}, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\rm{ }|\nabla {\mathit{\boldsymbol{\phi}}} |>1 \\ \end{align} \right. $

(8)

其对应的扩散率为：

$ r({\mathit{\boldsymbol{\phi}}} ) = \left\{ \begin{array}{l} \frac{{\sin (2{\rm{ \mathsf{ π} }}|\nabla {\mathit{\boldsymbol{\phi}}} |)}}{{2{\rm{ \mathsf{ π} }}|\nabla {\mathit{\boldsymbol{\phi}}} |}}, \;\;\;|\nabla {\mathit{\boldsymbol{\phi}}} | \le 1\\ 1 - \frac{1}{{|\nabla {\mathit{\boldsymbol{\phi}}} |}}, \;\;\;\;\;\;|\nabla {\mathit{\boldsymbol{\phi}}} | > 1 \end{array} \right. $

(9)

将惩罚项嵌入到活动轮廓能量方程中，使水平集在迭代过程中能够自动保持|▽ϕ|=1这一符号距离函数的特性，完全避免了重新初始化，提高了效率。但该方法仍存在一些缺点，如对弱边界防漏能力有限、对噪声敏感等。

近年来，Zhang等^[13]提出了新的无重新初始化方法——反应扩散模型，其主要思想是在能量方程中加入了扩散项来正则化水平集函数：

$ \frac{{\partial {\mathit{\boldsymbol{\phi}}} }}{{\partial t}} = \varepsilon \Delta {\mathit{\boldsymbol{\phi}}} - \frac{1}{\varepsilon }L({\mathit{\boldsymbol{\phi}}} ) $

(10)

其中：等号右边第一项为扩散项，用于惩罚复杂的区域边界使水平集保持平滑，ε＞0且为较小的常量；Δ为拉普拉斯算子，避免分割过程中产生孤立的小区域；L(ϕ)为反应项，控制曲线的演变。该方法对弱边界图像以及噪声图像均有良好的效果，但在实际运用中，该方法的时间效率相对DRLSE方法较低。

2 提出的活动轮廓模型 2.1 贝塞尔滤波

贝塞尔滤波器是一种线性滤波器，与高斯滤波器相似，随着滤波器阶数的增加趋向于相同的形状，与有限阶近似的高斯滤波器相比，贝塞尔滤波比相同阶数的高斯有更好的塑造因子、平坦的相位延迟和平坦的群延迟。文献[14]提出了一种基于M-S函数的贝塞尔滤波：

$ B(x, y;v, w) = {K_0}(\sqrt {({{(x - v)}^2} + {{(y - w)}^2})/\mu } ) $

(11)

其中:(x, y), (v, w)∈R²; μ为常数; K₀为第二类修正的贝塞尔函数。K₀积分表达式为：

$ {K_0}(x) = \int_0^{ + \infty } {\cos (x\sinh (t)){\rm{d}}t = } \int_0^{ + \infty } {\frac{{\cos (xt)}}{{\sqrt {1 + {t^2}} }}{\rm{d}}t}; x > 0 $

(12)

该贝塞尔边缘检测滤波具有尺度不变性，与其他滤波相比有更好的整体性能，实验结果显示其算法的数值实现具有较低的复杂度。于是文献[15]在此基础上提出了贝塞尔滤波的正则化和标准化函数为：

$ B(x) = \left\{ \begin{array}{l} \frac{{{K_0}(|x|/\sqrt \mu )}}{{{K_0}(\varepsilon /2\sqrt \mu )}}, |x| > \varepsilon /2\\ 1, {\rm{ }}\;\;\;\;\;\;\;\;\;\;\;|x| \le \varepsilon /2 \end{array} \right. $

(13)

而截断的贝塞尔滤波窗口函数为：

$ B(r, d) = \left\{ \begin{array}{l} 1, \;\;\;\;\;\;\;\;\;\;0 \le r < \varepsilon /2\\ \frac{{{K_0}(r/\sqrt \mu )}}{{{K_0}(\varepsilon /\sqrt \mu )}}, \varepsilon /2 \le r < d\\ 0, \;\;\;\;\;\;\;\;\;\;r > d \end{array} \right.\ $

(14)

其中:r为自变量表示点到滤波中心的距离；d为窗口大小；ε、μ为大于0的参数，分别控制水平和竖直方向上的滤波形状。在本文应用中统一设置ε=2、μ=1、d=15。

2.2 提出的水平集能量方程

为了应对噪声，传统活动轮廓模型采用的图像平滑方法一般为高斯滤波^[16]卷积：

$ \hat I(x, y) = I(x, y) * {G_\sigma }(x, y) $

(15)

但高斯滤波在平滑图像的同时，也会对边界造成一定的模糊效果，与其他边缘检测滤波相比，贝塞尔滤波有更好的整体性能，且具有尺度不变性，因此本文采用式(14)提出的贝塞尔滤波窗口函数对图像进行平滑处理：

$ \hat I(x, y) = I(x, y) * B $

(16)

其中：I(x, y)为原始图像，$ \hat I(x, y)$为平滑后的图像。由于贝塞尔与高斯滤波有相似的波形，达到了相似的平滑效果，但在实际操作中有着更高的效率，且更好地保持了目标边界。然后分别对图像从x、y方向进行卷积：

$ \left\{ \begin{array}{l} {I_x} = \hat I(x, y) * B(x)\\ {I_y} = \hat I(x, y) * B(y) \end{array} \right. $

(17)

新的边缘停止函数为：

$ b(x, y) = \frac{1}{{1 + I_{^x}^2 + I_{^y}^2}} $

(18)

于是构建出新的能量函数为：

$ E(\mathit{\boldsymbol{C}})=\oint_{C}{g(\mathit{\boldsymbol{C}})\rm{d}\mathit{s}+\iint_{\mathit{inside}(\mathit{\boldsymbol{C}})}{\mathit{b}(\mathit{x}, \mathit{y})\rm{d}\mathit{x}\rm{d}\mathit{y}}} $

(19)

其中:等号右边第一项g(C)表示曲线C的梯度信息，对g(C)沿C的环路积分表示曲线自身能量，通过能量最小化可以控制曲线运动；等号右边第二项对边缘停止函数在曲线内部取积分，使得能量函数在边缘处达到最小值，从而在能量最小化过程中实现曲线运动停止到边缘位置。采用水平集的方法，用隐式表示曲线：

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{C}} = \{ x \in \mathit{\Omega }:{\mathit{\boldsymbol{\phi}}} (x) = 0\} \\ in(\mathit{\boldsymbol{C}}) = \{ x \in \mathit{\Omega }:{\mathit{\boldsymbol{\phi}}} (x) > 0\} \\ out(\mathit{\boldsymbol{C}}) = \{ x \in \mathit{\Omega }:{\mathit{\boldsymbol{\phi}}} (x) < 0\} \end{array} \right. $

(20)

可得能量函数的水平集表示法：

$ \begin{align} &E({\mathit{\boldsymbol{\phi}}} )=\iint_{\mathit{\Omega }}{\mathit{g}|\nabla {{H}_{\varepsilon }}({\mathit{\boldsymbol{\phi}}} )|}\rm{d}\mathit{x}\rm{d}\mathit{y}+ \\ &\ \ \ \ \ \ \ \ \ \ \ c\iint_{\mathit{\Omega }}{[1-{{H}_{\varepsilon }}({\mathit{\boldsymbol{\phi}}} )]}b(x, y)\rm{d}\mathit{x}\rm{d}\mathit{y} \\ \end{align} $

(21)

其中，c为受到边缘停止函数控制的收缩力，一般设置为常数，g为梯度算子。H_ε(ϕ)为阶跃函数，通常用其正则化表示形式：

$ {{H}_{\varepsilon }}({\mathit{\boldsymbol{\phi}}} )\rm{=}[1+2\arctan ({\mathit{\boldsymbol{\phi}}} /\varepsilon )/\rm{ \mathsf{ π} } ]/\rm{2} $

(22)

其中ε＞0，为较小的参数。对提出的能量函数式(21)应用梯度下降法，得到最小化能量方程：

$ \frac{\partial {\mathit{\boldsymbol{\phi}}} }{\partial t}={{\delta }_{\varepsilon }}({\mathit{\boldsymbol{\phi}}} ){\rm{div}}(g\nabla {\mathit{\boldsymbol{\phi}}} /|\nabla {\mathit{\boldsymbol{\phi}}} |)+{{\delta }_{\varepsilon }}({\mathit{\boldsymbol{\phi}}} )cb $

(23)

其中:div(·)表示散度，δ(ϕ)为狄拉克函数。δ(ϕ)正则化形式表示为：

$ {{\delta }_{\varepsilon }}({\mathit{\boldsymbol{\phi}}} )=H_{\varepsilon }^{'}({\mathit{\boldsymbol{\phi}}} )=\frac{1}{\rm{ \mathsf{ π}} }\centerdot \frac{\varepsilon }{{{\varepsilon }^{2}}+{{{\mathit{\boldsymbol{\phi}}} }^{2}}} $

(24)

对于水平集正则化处理方式，本文加入了式(10)中的扩散项，得到的正则化水平集方程为：

$ \frac{{\partial {\mathit{\boldsymbol{\phi}}} }}{{\partial t}} = {\lambda _1}{\delta _\varepsilon }({\mathit{\boldsymbol{\phi}}} )\left( {{\rm{div}}(g\nabla {\mathit{\boldsymbol{\phi}}} /|\nabla {\mathit{\boldsymbol{\phi}}} |) + cb} \right) + {\lambda _2}\Delta {\mathit{\boldsymbol{\phi}}} $

(25)

其中:等号右边第一项为反应项，控制曲线的运动与收敛；等号右边第二项为扩散项，用于保持水平集函数的稳定性。λ₁、λ₂分别为大于零的权重系数，该方程结合了GAC模型与RD方程的优点，同时加入了新提出的边缘停止函数，在实验中显示了较好的效果。

3 实验结果与分析

本文实验是在Matlab R2014a软件环境下完成的，电脑处理器为AMD Athlon Ⅱ X4 3.10 GHz，内存为4 GB，运行的操作系统为Windows10专业版。通过多幅图片从分割效果和分割效率上分别对比了GAC模型、RD模型与本文提出的模型。

3.1 边缘停止函数对比

图 1为边缘检测函数针对不同图像的效果对比，即式(2)与式(18)对比，高斯滤波参数为：σ=1.5，窗口大小为15;贝塞尔滤波参数为：ε=2，μ=1，窗口大小为15。图 1中，灰度不均匀的合成图像，图像像素大小为128×128；T形合成图像，灰度分布不均匀，且在部分区域存在弱边界，图像像素大小为127×96；弱边界的小汽车真实图像，图像像素大小为118×93；玩具真实图像，图像像素大小为160×120。

由图 1可以直观地看出，传统的边缘停止函数对于灰度不均匀图像虽然有效地检测出了对象边缘，但在部分区域受背景干扰较为严重，边界区分不明显，如T形图、玩具图。而对于弱边界图像，由于图像梯度不明显，传统的边缘停止函数并不能很好地区分出目标边界，为下一步的分割操作带来困难，如小汽车图像。另外由于高斯滤波的平滑效果，边缘区域较宽，信息有损失。而本文提出的方法可以很好地区分出边界区域，尤其是针对弱边界对象也有很好的效果，且背景干扰信息较少、区分度较高，从而可以从根本上提高分割结果的精确度。为了进一步显示提出的边缘停止函数的效果，将提出的边缘停止函数分别加入到不同的模型中进行仿真对比，结果如图 2所示。

图 2为不同的活动轮廓模型应用不同的边缘停止函数后的分割过程及结果的对比，对比图像为合成的噪声图像，且边界较为模糊，图像像素大小84×84。参数设置为：λ₁=0.1，λ₂=0.001，收缩力c=0.5，由于DRLSE模型容易穿过弱边界，需要设置较小的收缩力c=0.2，其他对比参数均采用文献[13]默认设置，且同一模型均采用相同的参数。每幅图片中最外层曲线为初始曲线，最内侧为最终分割曲线，中间为曲线演变过程。

由图 2中轮廓线迭代次数(每迭代200次显示一条闭合曲线)可以看出，在相同的参数下，本文所提方法作用于不同的模型后，轮廓线均能够更快地到达目标边界，有效地减少了迭代次数，提高了分割的时间效率。其中原始GAC模型由于自身局限，曲线提前收敛，未能正确地进行分割，改进后可以更好地拟合目标边界。原始的DRLSE模型对于弱边界敏感，其轮廓线在部分区域穿过了目标边界，而在部分区域未能正确地收敛，改进后在弱边界和凹陷区域均能够正确收敛。原始RD模型具有相对较好的分割效果，但在仿真实验中其速度较慢，改进后迭代次数明显减少，提高了时间效率。

3.2 提出的模型对比

为了验证本文提出模型的效果，实验采用图 1(a)中的图像将本文模型与多种不同的模型进行仿真分割对比，结果如图 3~4所示。

图 3为合成图像仿真对比，灰度不均匀图像部分区域存在梯度不明显，参数设置为：λ₁=1，λ₂=0.01，c=0.05；T形图像为弱边界图像，边界较为模糊，目标前景背景分界不明显，参数设置为：λ₁=0.1，λ₂=0.001，c=0.2。图 3(a)为统一设置的初始轮廓，图 3(b)~(e)为不同模型的分割结果，依次为GAC模型、DRLSE模型、RD模型与本文提出的模型。从图 3中可以看出，在相同条件下：对于灰度不均匀图像中间梯度不明显的区域，对比的模型均存在边界泄漏情况，而本文模型取得了较好的效果；对于T形图像，GAC模型未能正确分割，DRLSE模型分割结果穿过了T形低部弱边界区域，RD模型与本文模型取得了较好的分割效果。而从表 1数据可以看出，本文模型对于两幅图像均能够快速收敛到目标位置，减少了迭代次数，具有较高的时间效率。

图 4为真实图像分割对比，由图 4中可以看出:GAC模型对于两幅图像的分割结果均存在偏差，轮廓线没有收敛到正确的位置; DRLSE模型对于弱边界区域的效果不好，如小汽车尾部等区域边界梯度不明显，轮廓线穿过了该区域; RD模型对弱边界图像取得了较好的分割效果，但在玩具图像中不够精确，一些凹陷区域没有收敛到位; 本文提出的模型与其他模型相比，取得了更好的分割结果，在保证分割精确度的同时，轮廓线更快地到达了真正的目标边界，取得了较高的时间效率。

3.3 量化分析

为了更精确地对提出的模型进行度量，本文选择了Weizmann分割评估数据库进行量化分析，从中选取了几幅代表性的图片进行了对比, 结果如图 5所示，同时采用了F-Score算法作为评价标准进行客观分析，该算法公式表示如下：

$ \left\{ \begin{array}{l} P = TP/\left( {TP + FP} \right)\\ R = TP/\left( {TP + FN} \right)\\ F = 2 * PR/\left( {P + R} \right) \end{array} \right. $

(26)

其中：TP为正确的目标前景分割样本；FP为将目标背景误分割为前景的样本；FN为将目标前景误分割为背景的样本；P为精确率，表示正确的分割前景与总分割前景的比重，其值与是否过分割相关联；R为召回率，表示正确的分割前景与标准的目标前景样本的比重，其值与是否漏分割相关联；F为准确率，是对分割结果是否准确的整体评价指标。

图 5为DRLSE模型、RD模型、本文提出的模型与标准分割数据集的对比。从图 5中可以看出，本文提出的模型与另外两个模型相比，更接近标准的分割结果，出现错误分割或漏分割的情况较少。表 2为不同模型分割结果的准确率的对比，从数据中可以看出，本文模型的精确率P与召回率R均达到了90%以上，没有出现较大的漏分割或过分割，而整体的准确率F均达到了95%以上，与另外两个模型相比有更好的分割效果。

本文模型仍存在一些缺陷，如由于所基于的GAC模型的非凸性，水平集函数对初始化较为敏感，不同的轮廓初始位置会影响到最终的分割结果。如图 5中一些图像的分割，需要将初始位置设置在目标前景区域，并适当调整相应的参数才能达到最佳的分割效果，需要进一步的改进。

4 结语

本文在GAC模型与RD模型的基础上，首先利用贝塞尔滤波对图像进行平滑处理，然后利用贝塞尔边缘检测函数构建边缘停止项，同时结合反应扩散方程提出一种新的测地活动轮廓模型。该模型不仅利用了贝塞尔滤波在边缘检测上快速精确的优势，同时也高效融合了反应扩散方程的特点，避免水平集重新初始化过程，从分割效果和时间复杂度上都得到了优化。今后的主要工作是将新的边缘停止函数与更多活动轮廓模型结合，以适用于更多不同的图片分割。