0 引言

图像分割是图像分析到图像理解的关键步骤^[1]。活动轮廓模型是获取目标边界的重要工具之一，具有良好的封闭性与平滑性，被广泛应用于图像分割领域中。参数活动轮廓模型又称为Snake模型，自Kass等^[2]提出以来，已广泛应用于数字图像分析和计算机视觉领域^[3]。其基本原理是将图像分割视为轮廓曲线演化问题，通过最小化曲线能量泛函得到最终的目标边缘。

在Snake模型中，外部力在整个曲线收敛过程中具有重要作用。传统Snake模型的外部力为高斯势能力，存在两大难点：一是初始轮廓线必须靠近真实边界，二是难以收敛到凹部边界^[4-5]。众多学者对此进行了一系列的研究。Xu等^[4]针对这两个难点提出的梯度向量流(Gradient Vector Flow, GVF)Snake模型和广义梯度向量流(Generalized GVF, GGVF)Snake模型^[6]，均采用扩散的方式。这两个模型都将梯度信息扩散至平滑区域，具有捕获范围较大、能收敛到一般凹部边界的优势，但都不能收敛到细小凹部边界^[7]。在GVF模型的扩散过程中，法向外力分量在曲线收敛到凹部边界过程中起着更加重要的作用。Ning等^[8]提出的法向梯度向量流(Normalized GVF, NGVF)Snake模型，仅保留沿边缘法向的外力分量，缺失了切向的外力分量，易消除弱边缘且噪声鲁棒性差。Bhan等^[9]延用GGVF模型思想，提出了改进的GVF模型，能够促使轮廓曲线运动到目标凹部，但计算复杂度高。

传统Snake模型及其改进模型获取外力场的能量泛函都仅仅利用梯度大小，而忽略了外力场的方向信息。为此学者们做了大量的尝试。例如，Li等^[10]通过添加边缘保真项，使外力场方向与边缘方向一致，提出的边缘保护梯度向量流(Edge Preserving GVF, EPGVF)Snake模型有效防止了曲线穿过弱边缘，但这种模型的角点定位仍不够精确。胡学刚等^[11]在EPGVF模型的基础上通过构造新的边缘保真项提出了一种新的分割模型，提高了角点定位的精确率，但噪声鲁棒性较差。分析发现，以上模型仅改进保真项，它们的平滑项中拉普拉斯算子对外力场都产生各向同性光滑作用，不利于保护弱边缘。针对这些问题，研究者们又提出了一系列改进算法。例如, Wu等^[12]提出的Snake分割模型具有保护弱边缘、抗噪性能好等优势，但对角点的定位仍需改进。Amin等^[13]通过对医学图像的分析，在文献[12]的基础上引入新的权重系数提出了改进的自适应扩散流模型，该模型对医学图像具有计算复杂度较低，噪声鲁棒性好的优势，但易陷入局部极小。

如何在保留现有模型优势的同时，提高分割效果和分割效率，又具有保护弱边缘、抗噪性能好等优势是目前图像分割亟待解决的问题。寻找恰当的平滑项，构造新的Snake模型是解决该问题的关键。针对这些问题，本文提出了一种基于Snake模型的图像分割新算法。该算法构造新的平滑项替代传统拉普拉斯算子，从而强化沿边缘法线方向的外力；再考虑外力场方向信息，借助边缘保护项使外力场方向与边缘方向一致。实验结果表明，所提出的模型不仅能防止弱边缘泄露，还能促使轮廓线收敛到细小深凹边界，明显提高了抗噪性能和角点定位精度，具有更好的分割效果。

1 Snake模型

假设轮廓曲线为C(s)=(x(s, y(s), s∈[0, 1]，s表示弧长参数，该曲线在式(1)能量泛函的驱使下得到分割结果。

$ \begin{align} & {{E}_{\rm{Snakes}}}=\int_{0}^{1}{\left\{ \frac{1}{2} \right.}\left[ \alpha {{\left| {{C}^{'}}(s) \right|}^{2}}+\beta {{\left| {{C}^{''}}(s) \right|}^{2}} \right]+ \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left. {{\mathit{E}}_{\rm{ext}}}(\mathit{C}(\mathit{s})) \right\}\rm{d}\mathit{s} \\ \end{align} $

(1)

其中：α和β分别是控制Snake模型的弹性能量与刚性能量的权重系数。积分的第一项称为内部能量，只与曲线自身相关，使曲线在形变过程中保持连续性与光滑性；第二项称为外部能量，仅来源于图像，能够驱使曲线不断向目标边界运动。

由变分原理可知，式(1)最小化必要条件是要求C(s)必须满足下列欧拉方程：

$ \alpha {{C}^{''}}(s)-\beta {{C}^{''''}}(s)-\nabla {{E}_{\rm{ext}}}=0 $

(2)

引入时间t得到如下动态方程：

$ {{C}_{t}}(s, t)=\alpha {{C}^{''}}(s, t)-\beta {{C}^{''''}}(s, t)-\nabla {{E}_{\rm{ext}}} $

(3)

当式(3)求得稳态解，即为式(2)的解。

GVF模型引入梯度向量流V(x, y)=(u(x, y), v(x, y))代替传统Snake模型的外部力-▽E_ext，通过最小化下列能量泛函获得外力场：

$ \begin{align} & \varepsilon =\iint{\mu (u_{x}^{2}+u_{y}^{2}+v_{x}^{2}+v_{y}^{2})}+ \rm{ }{{\left| \nabla \mathit{f} \right|}^{2}}{{\left| \mathit{V}-\nabla \mathit{f} \right|}^{2}}\rm{d}\mathit{x}\rm{d}\mathit{y} \\ \end{align} $

(4)

其中：μ是控制外力场平滑程度的权重系数；f(x, y)是图像的边缘映射。当曲线远离目标轮廓时，第一项(平滑项)起主导作用; 反之，第二项(数据项)起主导作用，从而扩大捕获范围。

根据变分原理，通过解下列欧拉方程得到GVF场：

$ \left\{ {\begin{array}{*{20}{c}} {\mu {\nabla ^2}u - (u - {f_x})(f_x^2 + f_y^2) = 0}\\ {\mu {\nabla ^2}v - (v - {f_x})(f_x^2 + f_y^2) = 0} \end{array}} \right. $

(5)

式中▽²表示拉普拉斯算子。

由文献[8]可知，式(5)中的拉普拉斯算子对外力场产生各向同性光滑作用，不能保护边缘。根据图像局部结构，将拉普拉斯算子分解成法向和切向两个扩散分量。其中法向扩散力在曲线收敛到凹部边界过程中起主要作用，由此得到如下NGVF模型相应的欧拉方程：

$ \left\{ {\begin{array}{*{20}{c}} {\mu {u_{nn}} - (u - {f_x})(f_x^2 + f_y^2) = 0}\\ {\mu {v_{nn}} - (v - {f_x})(f_x^2 + f_y^2) = 0} \end{array}} \right. $

(6)

式中u_nn和v_nn是法向扩散力。

针对GVF模型与NGVF模型难以检测到弱边缘的问题，EPGVF模型在数据项中增加了μ|J_VP|²项，这里J_V是外力场V的雅可比矩阵，并且保留了水平方向扩散力。通过最小化式(7)获得外力场：

$ \begin{align} & \varepsilon =\iint{\mathit{g}(x, y){{\left| \nabla V \right|}^{2}}+h(x, y)(\mu {{\left| {{J}_{V}}{\mathit{\boldsymbol{P}}} \right|}^{2}}} +\\ & \rm{ }{{\left| \mathit{V}-\nabla \mathit{f} \right|}^{2}})\rm{d}\mathit{x}\rm{d}\mathit{y} \\ \end{align} $

(7)

其中:μ、g和h是权重函数; $ {\mathit{\boldsymbol{P}}}={{\left[-{{{I}_{y}}}/{\sqrt{I_{x}^{2}+I_{y}^{2}}}\;, {{{I}_{x}}}/{\sqrt{I_{x}^{2}+I_{y}^{2}}}\; \right]}^{\rm{T}}}$表示边缘方向。

式(7)对应的欧拉方程如下：

$ \left\{ {\begin{array}{*{20}{c}} \begin{array}{l} g(x, y){\nabla ^2}u - h(x, y)(u - {f_x} - \mu P_1^2{u_{xx}} - \\ \;\;\;\;\mu P_2^2{u_{yy}} - 2\mu {P_1}{P_2}{u_{xy}}) = 0 \end{array}\\ \begin{array}{l} g(x, y){\nabla ^2}v - h(x, y)(v - {f_x} - \mu P_1^2{v_{xx}} - \\ \;\;\;\;\mu P_2^2{v_{yy}} - 2\mu {P_1}{P_2}{v_{xy}}) = 0 \end{array} \end{array}} \right. $

(8)

解式(8)可得到外力场。

2 本文模型

EPGVF模型的扩散过程主要依赖于拉普拉斯算子，而拉普拉斯算子对外力场具有很强的各向同性光滑作用，不利于保护弱边缘。针对这个问题，本文选取(1/(q(|▽f|)))·(1+|G_σ⊗▽V|²)^q(|▽f|)/2代替EPGVF模型平滑项中的|▽V|²，其中q(x)为单调递减函数，取q(|▽f|)=1+1/(1+(|▽G_σ*f(x)|))。在图像平滑区域中，由于|▽G_σ*f(x)|→0，因此q(|▽f|)→2，外力场具有各向同性扩散作用；在图像边缘处，由于|▽G_σ*f(x)|→∞，因此q(|▽f|)→1，$ \sqrt {1 + {{\left| {{G_\sigma } * \nabla V} \right|}^2}} \to \left| {\nabla V} \right|$，外力场只沿着与边界一致的方向进行扩散，具有防止弱边缘泄露并提高抗噪性能的作用。

EPGVF模型的边缘保护项μ|J_VP|²可以使外力沿边缘方向平滑，却易陷入局部极小并且角点定位不精确。对此，本文首先构造新的泛函促使曲线收敛到细小深凹边界。由NGVF模型可知，拉普拉斯算子可分解成沿法线方向和切线方向两个分量，并且法向扩散分量具有促使轮廓线收敛到深凹部位的作用。为此，本文利用如下p-拉普拉斯泛函鼓励法向外力场分量：

$ {E_p}(V) = \frac{1}{p}{\int_\mathit{\Omega } {\left| {{G_\sigma } \otimes \nabla V} \right|} ^p}{\rm{d}}\mathit{\Omega } $

(9)

最小化式(9)必须满足下列欧拉方程：

$ {\rm{div}}({\left| {\nabla V} \right|^{p - 2}}\nabla V) = 0 $

(10)

受文献[11]的启发，本文给出如下雅可比矩阵来改进边缘保护项，以提高角点定位精度:

$ \mathit{\boldsymbol{P}}=\left[\begin{matrix} {-{{I}_{xy}}}/{\sqrt{I_{xx}^{2}+I_{xy}^{2}}}\; & {{{I}_{xx}}}/{\sqrt{I_{xx}^{2}+I_{xy}^{2}}}\; \\ {-{{I}_{yy}}}/{\sqrt{I_{yx}^{2}+I_{yy}^{2}}}\; & {{{I}_{yx}}}/{\sqrt{I_{yx}^{2}+I_{yy}^{2}}}\; \\ \end{matrix} \right] $

(11)

由于EPGVF模型平滑项的权重系数均通过人工选取，具有较大的主观性与局限性，影响实际应用。本文引入与图像梯度相结合的权重系数m和1-m，对于平滑项参数实现了自动选取，拓宽模型的适用范围并提高分割精度。

$ m=\left\{ \begin{align} & [1-{{{f}^{2}}}/{6{{K}^{2}}{{]}^{2}}, {\ {{f}^{2}}}/{6}\;\le {{K}^{2}}}\; \\ & 0, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \rm{其他} \\ \end{align} \right. $

(12)

式中K=1.482 6E(||▽f|-E(|▽f|)|), E(·)为均值。

根据以上分析，得到本文的能量泛函如下：

$ \begin{align} & \varepsilon =\iint{g(\left| \nabla f \right|)(-m\left| {{G}_{\sigma }}*\nabla V \right|_{{{L}^{\infty }}}^{2}+(1-m)} \\ & \;\;\;\;{{({{{(1}/{(q(\left| \nabla f \right|)))(1+\left| {{G}_{\sigma }}\otimes \nabla V \right|}\;}^{2}})}^{{q(\left| \nabla f \right|)}/{2}\;}})) +\\ &\;\;\;\; \rm{ }\mathit{h}(\left| \nabla \mathit{f} \right|)(\mathit{\mu }{{\left| {{\mathit{J}}_{\mathit{V}}}\mathit{\boldsymbol{P}} \right|}^{2}}+{{\left| \mathit{V}-\nabla \mathit{f} \right|}^{2}})\rm{d}\mathit{x}\rm{d}\mathit{y} \\ \end{align} $

(13)

式中：g、h、m和μ均为权重系数，g和h的选取方式与GGVF模型[6]一致；$ \iint{{}}$|Gσ*▽V|L∞2表示p→∞时的式(9)拉普拉斯泛函。若f2/6＞K2, 则m=0，此时，(1/(q(|▽f|)))(1+|Gσ⊗▽V|2)q(|▽f|)/2项控制外力场的扩散。当m=1时，本文模型平滑项近似于NGVF模型平滑项；若f2/6≤K2, 则m→1，|Gσ*▽V|L∞2项在外力场扩散过程中起主要作用，强化了深凹边界处的法向分量，此时力场向量的方向不是相互对立而是向下的，从而促使轮廓曲线更快的收敛到细小深凹边界。

由变分原理可得式(13)相应的欧拉方程组：

$ \left\{ \begin{align} & g(\left| \nabla f \right|)[(1-m)({{(1+{{\left| \nabla u \right|}^{2}})}^{{(q-2)}/{2}\;}}{{u}_{tt}}+((q-2) \\ & \ \ {{\left| \nabla u \right|}^{2}}{{(1+{{\left| \nabla u \right|}^{2}})}^{{(q-4)}/{2}\;}}+{{(1+{{\left| \nabla u \right|}^{2}})}^{{(q-2)}/{2}\;}}){{u}_{nn}})+ \\ & \ \ m{{u}_{nn}}]-h(\left| \nabla f \right|)(u-{{f}_{x}}-\mu P_{11}^{2}{{u}_{xx}}-\mu P_{12}^{2}{{u}_{yy}}- \\ & \ \ 2\mu {{P}_{11}}{{P}_{12}}{{u}_{xy}})=0 \\ & g(\left| \nabla f \right|)[(1-m)({{(1+{{\left| \nabla v \right|}^{2}})}^{{(q-2)}/{2}\;}}{{v}_{tt}}+((q-2) \\ & \ \ {{\left| \nabla v \right|}^{2}}{{(1+{{\left| \nabla v \right|}^{2}})}^{{(q-4)}/{2}\;}}+{{(1+{{\left| \nabla v \right|}^{2}})}^{{(q-2)}/{2}\;}}){{v}_{nn}})+ \\ & \ \ m{{v}_{nn}}]-h(\left| \nabla f \right|)(v-{{f}_{y}}-\mu P_{21}^{2}{{v}_{xx}}-\mu P_{22}^{2}{{v}_{yy}}- \\ & \ \ 2\mu {{P}_{21}}{{P}_{22}}{{v}_{xy}})=0 \\ \end{align} \right. $

(14)

类似于GVF模型获得数值解的方法，容易得到如下的本文模型的数值解的迭代公式：

$ \left\{ \begin{align} & u_{i, j}^{n+1}=(1-\Delta th(\left| \nabla f \right|))u_{i, j}^{n}+\Delta tg(\left| \nabla f \right|)[(1-m) \\ & \ \ ({{(1+{{\left| \nabla u \right|}^{2}})}^{{(q-2)}/{2}\;}}{{u}_{tt}}+((q-2){{\left| \nabla u \right|}^{2}} \\ & \ \ {{(1+{{\left| \nabla u \right|}^{2}})}^{{(q-4)}/{2}\;}}+{{(1+{{\left| \nabla u \right|}^{2}})}^{{(q-2)}/{2}\;}}){{u}_{nn}})+ \\ & \ \ m{{u}_{nn}}]+\Delta t{{c}_{1}}+ \\ & \ \ \ \Delta t\mu h(\left| \nabla f \right|)(\mu P_{11}^{2}{{u}_{xx}}+\mu P_{12}^{2}{{u}_{yy}}+2\mu {{P}_{11}}{{P}_{12}}{{u}_{xy}}) \\ & v_{i, j}^{n+1}=(1-\Delta th(\left| \nabla f \right|))v_{i, j}^{n}+\Delta tg(\left| \nabla f \right|)[(1-m) \\ & \ \ \ ({{(1+{{\left| \nabla v \right|}^{2}})}^{{(q-2)}/{2}\;}}{{v}_{tt}}+((q-2){{\left| \nabla v \right|}^{2}}(1+ \\ & \ \ \ {{\left| \nabla v \right|}^{2}}{{)}^{{(q-4)}/{2}\;}}+{{(1+{{\left| \nabla v \right|}^{2}})}^{{(q-2)}/{2}\;}}){{v}_{nn}})+m{{v}_{nn}}]+ \\ & \ \ \Delta t{{c}_{2}}+\Delta t\mu h(\left| \nabla f \right|)(\mu P_{21}^{2}{{v}_{xx}}+\mu P_{22}^{2}{{v}_{yy}}+ \\ & \ \ \ 2\mu {{P}_{21}}{{P}_{22}}{{v}_{xy}}) \\ \end{align} \right. $

(15)

式中：c₁=h(|▽f|)f_x(x, y)；c₂=h(|▽f|)f_y(x, y)；Δt为迭代时间步长；(u_{i, j}ⁿ, v_{i, j}ⁿ)表示迭代n次时坐标为(i, j)的外力场。

3 实验结果及分析

以人工合成图像、心脏CT图像、脑部CT图像、鲜花和鸟自然图像作为实验对象，分别用文献[5]模型、文献[14]模型、文献[11]模型与本文模型进行对比实验。实验平台为Windows 7操作系统，IntelI Celeron E5700处理器，主频2.60 GHz，3.00 GB内存，Matlab R2012a。各模型参数的选取方法见下列对应实验，其中对比模型的参数按相关文献的要求选取。

首先，选取64×64含有弱边缘的合成图像为实验对象，如图 1(a)所示，其右边边缘为弱边缘，不远处还有一条黑色的强边缘线。参与比较的四种模型均选取参数α=0.04，β=0.4，τ=0.5；本文模型其余参数δ=0.01、μ=0.1；文献[5]模型μ=0.1、k=0.36；文献[14]模型μ=0.2；文献[11]模型μ=0.1、k=0.08。图 1(a)的圆形曲线为初始轮廓线，图 1(b)为分割的金标准，图 1(c)~(f)分别为文献[5]模型、文献[14]模型、文献[11]模型和本文模型的分割结果。从图 1可以看出：文献[5]和[11]模型在拐角处都未能准确分割；文献[14]模型受到了强边缘的干扰导致分割错误；而本文模型能够较为准确地分割出拐角边界并且很好地检测出了弱边缘。

其次，选取128×128合成图像，并加入了均值为0、方差为30的高斯噪声，如图 2(a)所示, 以此验证模型的抗噪性能。比较的四种模型均选取参数α=0.5，β=0.5，τ=0.5；本文模型其余参数δ=0.6，μ=0.1；文献[5]模型其余参数μ=0.18，k=0.36；文献[14]模型μ=0.01；文献[11]模型μ=0.1，k=0.08。图 2(a)的圆形曲线为初始轮廓线，图 2(b)为分割的金标准，图 2(c)~(f)分别为文献[5]模型、文献[14]模型、文献[11]模型和本文模型的分割结果。从图 2可以看出：文献[5]模型在右上角处受噪声影响出现误分割现象；文献[14]模型在上方、下方、左方、右方均存在欠分割并且在斜方向的拐弯连接处出现过分割现象；文献[11]模型在下方和右上角方向存在欠分割，在右上角出现过分割现象；而本文模型的分割结果没有出现这些缺陷，分割效果最佳，很好地保护了图像细节，并具有良好的抗噪性能。

图 3(a)中的128×128人体左心室的CT图像带有噪声、图像灰度不均匀且存在弱边缘，选取该图为实验对象，以验证模型的有效性和实用性。按照图 2对应实验选取四种模型参数α、β和τ以及文献[5]、[11]、本文模型的其余参数，另外文献[14]模型其余参数μ=0.1。图 3(a)的圆形色曲线为初始轮廓线，图 3(b)为分割的金标准，图 3(c)~(f)分别为文献[5]模型、文献[14]模型、文献[11]模型和本文模型的分割结果。从图 3可以看出：文献[5]模型在左上角处受强边缘影响而导致错误分割；文献[14]模型收敛的最终轮廓不够流畅并且在左方存在欠分割现象；文献[11]模型在右下角和左下角都受到噪声影响，导致未能准确分割；而本文模型能收敛到拐角处，分割效果好，既不受噪声影响，也能够检测出弱边缘。

再选取443×443人体脑部CT图像作为实验对象。除本文模型中参数δ取2以及文献[14]参数μ取0.17之外，四种模型的其余参数同图 2对应实验一致。图 4(a)的圆形曲线为初始轮廓线，图 4(b)为分割的金标准，图 4(c)~(f)分别为文献[5]模型、文献[14]模型、文献[11]模型和本文模型的分割结果。从图 4可以发现：文献[5]模型在正下方受到中间无关边缘的影响，出现过分割现象；文献[14]模型在左下方和右下方存在欠分割，并受到强边界影响，左下方出现部分错误分割；文献[11]模型在拐角的定位上较为精准，但也受到下方强边界干扰，出现了过分割现象；而本文模型能够正确分割四个长窄凹面，不受过渡区域周围强边缘和无关边缘的影响，准确分割出边缘，并且在拐角处定位较为精准，得到了很好的分割效果。

图 5(a)是鲜花自然图像，本次选取的图像具有较为复杂的背景，且背景虚化出现了很多的弱边缘。在实验中，四种模型均选取参数α=0.5，β=1.5，τ=0.5；本文模型其余参数δ=0.35，μ=0.1；文献[5]和[11]模型参数μ、k均与图 2对应实验一致；文献[14]模型参数μ=0.25。图 5(a)的圆形曲线为初始轮廓线，图 5(b)为分割的金标准，图 5(c)~(f)分别为文献[5]模型、文献[14]模型、文献[11]模型和本文模型的分割结果。显然，本文模型分割效果最好，能够很好地保护图像细节，正确分割花瓣形状，不仅不受弱边缘的影响，而且能够深入花瓣之间的凹缝。

此外，本文方法也适合对彩色自然图像的分割。按图 2实验的方法选取参数α、β、μ，τ和δ分别取0.7和0.5，用本文模型对其进行分割实验，结果如图 6所示。从图 6可以看出，本文模型对于颜色丰富、背景复杂的彩色图像也可以得到理想的分割结果。

现在，利用边界位移误差(Boundary Displace Error, BDE)客观指标评价分割精度^[15]。BDE的计算公式为：

$ {BDE \rm{= }}\frac{1}{N}\sum\limits_{C \in {B_2}} {\mathop {\min }\limits_{e \in {B_1}} } \{ d(c, e)\} $

(16)

其中：B₁为真实的边界曲线像素点集合；B₂为某种模型分割得到的边界像素点集合；N为B₂中像素点的个数；d(i, j)表示两个像素点i和j的欧氏距离。金标准分割BDE值为0，BDE值越小分割效果越好。表 1给出了图 1~5对应实验的BDE值，从中可以看出，本文模型的BDE值明显小于其他三种模型，与上述实验的视觉比较结果一致。

关于模型的计算复杂度，表 2和3分别给出了在相同软硬件实验平台上的关于图 1~5对应实验的模型迭代次数和收敛时间。容易发现，本文模型收敛速度快，其迭代次数明显少于其他三种模型，平均耗时低于文献[11]模型，比另两种模型减少了26.7%。

4 结语

本文提出了一种基于Snake模型的图像分割新算法，首先，该算法新构造的平滑项强化了沿边缘法线方向的外力；其次，利用边缘保护项保留了外力场方向信息；最后，自动选取平滑项参数，克服了传统人工选取参数的缺陷。实验结果表明，所提出的模型不仅保持了现有模型捕获范围大、初始位置不敏感等优势，还有效防止弱边缘泄露且可以收敛到细小凹部边界，明显提高了抗噪性能和角点定位精度，具有更好的分割效果。