基于变分模型的图像分割方法^[1-2]能够提供光滑封闭的曲线，并能够结合先验信息获得非同质的目标^[3]区域，因而被广泛地应用于目标检测、视觉跟踪、场景理解、图像检索以及医学图像分析^[4]等领域。

基于变分模型的图像分割方法一般分为两大类：基于边的分割方法以及基于区域的分割方法。基于边的分割方法主要利用图像的局部梯度吸引活动轮廓朝着图像的边缘方向进化^[5-6]。虽然该方法在特定应用环境下能够获得较好的分割结果，但对初始标记点以及初始轮廓的位置比较敏感，且在缺乏全局信息的约束下，容易受噪声的干扰。基于区域的分割方法，它假设每个区域内部的统计密度是同质的，依靠全局信息来引导活动轮廓的演化。相对基于边的分割方法^[7-9]而言，它对初始轮廓的位置不敏感，因此具有较强的抗噪声能力。但它丢失了图像的全局刻画，导致分割的结果容易出现虚假的目标区域，且捕获的目标边缘不够光滑。

文献[10]提出将边缘信息(梯度)与区域信息相结合的方法，但在进行边缘结合时，并未引入边缘检测能力更强的测地线项。且自然图像通常包含多个非同质目标区域^{[7, 11]}，各个区域内部的概率密度变化具有非线性及连续性，若采用恒定的聚类中心描述，则难以实现可靠、有效的多类自然图像分割。

本文通过自由度调和的多变量混合学生-t概率密度分布来提高多类非同质目标区域的非线性以及连续性刻画^[12]。为了提高抗噪能力及光滑边缘的检测能力，本文引入了测地线边缘外力项。由于本文提出的多类变分活动轮廓模型的最小化问题是NP难问题，如果采用多类水平集的方法^[13]进行求解，其收敛较慢，容易陷于局部最小。为了打破此瓶颈，本文提出利用微分几何与积分几何的相关理论，建立最小割与测地线之间的关系，将多类变分活动轮廓模型的最小化问题转化为多层图割模型^[14-16]的最大流/最小割问题，这不仅能够提高自然图像的分割速度，而且能够求得全局近似最优解。

假设I₀是一幅自然图像，其对应的图像域$\Omega :{R^2} \to R$。$\Omega $对应的分割区域边界子集为C，它将自然图像I₀分割为若干互不联通的子区域${\Omega _l}$，满足：$\Omega = \mathop \cup \limits_{l = 1}^{N{\kern 1pt} (\Omega )} {\Omega _l} \cup C$，且$\mathop {{\Omega _l} \cap {\Omega _k}}\limits_{l \ne k} = \emptyset $，其中$N(\Omega )$表示所有不连通的分割区域数。图像I₀常被假设为m个类别区域$\{ \Omega {{\kern 1pt} _l}\left| {l = 1,2, \cdots ,m} \right.\} $，如果每个区域$\Omega {{\kern 1pt} _l}$内部用$n{{\kern 1pt} _l}$个恒定聚类中心描述，则它丢失了区域$\Omega {{\kern 1pt} _l}$内部存在的非线性与连续性变化。为了更好地描述特征空间的变化，本文引入了多变量混合学生-t分布进行描述，它通过自由度调节参数v来调节概率密度分布的形状变化。在$\nu \to + \infty $时，多变量混合学生-t分布退化为多变量高斯分布。则本文提出的多类变分模型可表达为：

式中，$f({I_0}(x)\left| {{\Theta _k}} \right.)$为多变量学生-t分布，它可对区域${\Omega _k}$内部的特征进行非线性与连续性概率密度分布描述；m为总区域数，即图像中非同质目标区域的个数；${\omega _k}$为对应的区域权重；${\Theta _k} = \{ {\nu _k},{\mu _k},{\Sigma _k}\} $为区域$\Omega {{\kern 1pt} _k}$内部的相关统计参数；${\mu _k}$是区域${\Omega _k}$的均值向量；${\Sigma _k}$是对应的协方差矩阵；${\nu _k}$是自由度调节参数；$\Gamma (x) = \int_0^{ + \infty } {{t^{x - 1}}{{\rm{e}}^{ - t}}{\rm{d}}{\kern 1pt} t} $为Gamma函数。对于区域数m及相关的统计参数$\Theta {\kern 1pt} = \left\{ {{\omega _k},{\Theta _k}\left| {k = 1,2, \cdots ,m{\kern 1pt} } \right.} \right\}$，可利用CEM³ST算法^[9]进行初始化。

式(1)中，$\gamma \oint {_C{\rm{d}}s} $虽对边缘项进行了长度约束，但它对噪声干扰以及区域间边界的捕获能力有限，因此极易陷入局部最小。而文献[5]提出了geodesic active contour(GAC)测地线活动轮廓模型，它具有较少的量化参数，且能够吸引活动轮廓沿边界法线的方向进化，因而被广泛地应用于灰度图像分割，其对应的能量函数为：

式中，$g(C(s))$为单调递减的边缘检测函数，它在区域边界和脊岭处的边界约束能力较强，通常$g(x) = {(1 + \nabla \left| {{I_0}(x)} \right|)^{ - 1}}$。因此，可将$g(C(s))$区域外力扩展应用于改进的多类变分模型式(1)中，即：

通过区域项与边缘项的结合，式(4)中所对应的能量函数不仅能够有效约束活动轮廓的边界长度，也可以进行区域边界的有效检测，共同提高非同质目标区域的整体描述能力。但式(4)所提出的多类变分活动轮廓能量函数是变分形式，计算复杂且最小化是NP难问题，所以很难直接对其进行优化。文献^[7]提出了多类水平集优化方式，但它需将目标能量函数嵌入到高维空间，且要求被处理的区域数必须是偶数。因此，很难将其应用于优化式(4)所对应的能量函数。而本文提出了利用积分几何与微分几何之间的关系对所提出的多类变分活动轮廓模型进行离散化推导与表达，进而建立测地线与最大流/最小割间的联系，通过构建多层图割模型来进行快速优化求解。

式(4)中对应的多类活动轮廓能量函数为变分形式，由于变分能量函数复杂，难以直接将其运用于二维离散图像，因此首先需要进行离散化表达。

假设图像I₀可表示为二维离散格图$S = \left\{ {s(i,j)\left| {i \in \left\{ {1,2, \cdots ,W} \right\},j \in \left\{ {1,2, \cdots ,H} \right\}} \right.} \right\}$，W与H分别对应于图像的宽和高。为了便于离散化描述，首先引入一个辅助函数$\delta (x)$，如果$x \ne 0$，$\delta (x) = 0$；否则$\delta (x) = 1$。在二维格图S上，位于s位置的类别标签为$\Upsilon (s)$，则S对应的标签图可定义为$\left\{ {\Upsilon (s){\kern 1pt} \left| {\Upsilon (s){\kern 1pt} \in \left\{ {1,2, \cdots ,m} \right\} \cap s \in S} \right.} \right\}$，则式(4)中的区域项可离散化为：

式中，函数$f( \cdot )$描述了s位置的像素隶属于第k类区域的概率相似度。对于边缘项，文献[17]提出了切分准则与欧几里得长度之间的关系，本文可将其作进一步的扩展，用于求带测地线边缘项的离散化近似。在二维格图上边界曲线C的长度，可利用Cauchy-Crofton公式进行近似，即${\kern 1pt} \oint {_C} {\rm{d}}s = 2{\left| C \right|_E}$。对应于不同的邻域系统Q，可得到边界曲线C的不同精度逼近，图 1为4领域系统中边界曲线C的长度近似。其具体的近似过程如下：假设在二维平面上所有直线可表示为${\kern 1pt} \left\{ {(\rho ,\theta )\left| {x{\rm{cos}}} \right.\theta + y{\rm{sin}}\theta = \rho } \right\}$(采用极坐标方式)，即边界曲线C的欧几里得长度可计算为：

图1

图1 4领域系统中曲线C的近似长度

式中，${\kern 1pt} {n_c}(\rho ,\theta )$表示二维平面上直线$x{\kern 1pt} {\rm{cos}}\theta + y{\kern 1pt} {\rm{sin}}\theta = \rho $与边界曲线C相交的总次数；${\kern 1pt} (\rho ,\theta )$对应于一个直线族(族内部的直线间彼此相互平行)。对于$Q = {n_G}$邻域，图 1中二维格图对应的空间平面邻域系统可表示为${N_Q} = \left\{ {{\kern 1pt} {e_k}{\kern 1pt} {\kern 1pt} \left| {{\kern 1pt} 0 \le {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} k{\kern 1pt} {\kern 1pt} {\kern 1pt} \le {n_G}{\kern 1pt} } \right.} \right\}$，最近邻两个向量之间的夹角可设为$\Delta {\theta _k} = 2\pi /{n_G}$，如4邻域系统的$\Delta {\theta _k} = \pi /4$。利用Cauchy-Crofton公式可近似计算曲线长度$\left| {{\kern 1pt} {\kern 1pt} C{\kern 1pt} {\kern 1pt} } \right|{\kern 1pt} {{\kern 1pt} _E}$，则式(6)中的Q邻域欧几里得长度可离散化为${\left| C \right|_E} = \sum\limits_{k = 1}^{{n_G}} {{n_c}(k)} {\kern 1pt} {\delta ^2}\Delta {\theta _k}/(2{\kern 1pt} \left| {{e_k}} \right|)$。这里，$\left| {{\kern 1pt} {e_k}} \right|$为二维格图上方向向量的长度。对于方向相同的直线，其权重${\omega _k} = {\delta ^2}\Delta {\theta _k}/\left( {2\left| {{e_k}} \right|} \right)$，即曲线长度${\left| C \right|_E}$可简化为$\left| {{\kern 1pt} {\kern 1pt} C{\kern 1pt} {\kern 1pt} } \right|{\kern 1pt} {{\kern 1pt} _E}{\kern 1pt} = \sum\limits_{k = 1}^{{n_G}} {{\kern 1pt} {\kern 1pt} {\kern 1pt} {n_c}(k)} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\omega _k}$。

考虑二维网络流图$G = (V,E,W)$，其中V是顶点的集合，E是边的集合，W是对应于边集E的非负权重集合。在V中包含两个特殊的端点S和T(S为源点，T为汇点)。根据切分准则，边界曲线C可将S与T切割开，则对应的切割代价可表示为：

式中，${(s,{\kern 1pt} r)}$表示图G被边界曲线C切割开的边，其对应的边权重为${{w_{s,r}}}$；而${\left| C \right|_G}$表示边界曲线C切割开边的总代价。由于${\left| C \right|_E}$与${\left| C \right|_G}$都表示边界曲线C的总切割代价，则在平面离散化的情况下，它们近似相等，即：

在多类的情况下，可根据图G对应的多个分割边界曲线C，以及对应的Q邻域标签区域来共同离散化$\oint {_C} {\kern 1pt} {\rm{d}}s$，即：

式中，${N_Q}(s)$表示s点的Q邻域点集，利用式(7)可离散化式(4)中的边缘项能量函数，即：

边缘项${E_{{\rm{edge}}}}$中增加了测地线边缘检测函数${g(\left| {\nabla C(s)} \right|)}$，主要用于相邻两个像素被分配不同类别标签时利用边缘梯度信息来分配惩罚强度。如在平坦区域，当相邻两个像素被分配不同的类别标签时，其惩罚强度最大。关于${g(\left| {\nabla C(s)} \right|)}$函数的离散化，可通过二维格图上Q邻域点对之间的梯度计算(本文选取Q为4)，则格点s(i,j) 处的梯度向量为：

即对应于$\nabla C(s)$的离散测地线函数为：

为方便计算，式(5)～式(8)中s与r都未写成格点形式(实际执行时需采用格点形式计算)。为了捕获光滑的边界，提高抗噪能力，需要将离散的区域内力与边缘外力相结合，即离散化的多类变分活动轮廓为：

式中，${w_{s,r}}$用于描述邻域像素对之间的差异性与空间约束，采用欧氏距离度量，即：

式(7)中对应的多类离散能量函数具有凸函数形式，且边缘项为多标签的Potts模型，因此很难直接对其最小化求解。而文献[18]提出的craph cuts优化方法，能够快速求解凸函数的最小化问题，因此可将其扩展于求解式(9)对应的多类能量函数。对于能量函数的最小化问题可将其转化为对应多层图割模型的最大流/最小割问题。首先，构建一个多层图$G = (V,E,W)$(多层图的网络流从源点(source)流向汇点(sink))，它的每一层对应于一个二维格图S，其上每个点s对应于自然图像I₀中的一个像素，层与层之间的边表示像素隶属于某一类的相似度。对于s点，它可能被分配到m类中的任意一类，因此多层图G就对应于一个三维格图，它可定义为$\left\{ {(i{\kern 1pt} ,{\kern 1pt} s{\kern 1pt} ) \in {R^2} \times R{\kern 1pt} {\kern 1pt} {\kern 1pt} \left| {{\kern 1pt} {\kern 1pt} i \in \left\{ {1,2, \cdots ,m - 1} \right\},s \in S} \right.} \right\}$。对于多层图上任意一点(i,s)可表示为v_i,s，相邻图层间的边$e({v_{i{\kern 1pt} ,{\kern 1pt} {\kern 1pt} s{\kern 1pt} }},{v_{i + 1{\kern 1pt} ,{\kern 1pt} {\kern 1pt} s{\kern 1pt} }})$表示t-link，同层格图上的边$e({v_{i{\kern 1pt} ,{\kern 1pt} {\kern 1pt} s{\kern 1pt} }},{v_{{\kern 1pt} i,{\kern 1pt} {\kern 1pt} r{\kern 1pt} }})$表示n-link。对应于E的顶点集V可表示为$\left\{ {{\rm{source}}} \right\} \cup _{i = 1}^m\left\{ {{v_{i,{\kern 1pt} {\kern 1pt} s}}\left| {s \in S} \right.} \right\} \cup \left\{ {{\rm{sink}}} \right\}$。因此，式(7)中的区域项能量E_reg与边缘项能量E_edeg，可利用顶点集V中的点分别表示为t-link边集Edge__T与n-link边集(边缘项边集)Edge__N。则区域项边集E__T可表示为：

${\rm{Edg}}{{\rm{e}}_{\_T}} = \cup \left\{ {{E_s}\left| {s \in S} \right.} \right\}$

${E_s} = e({\rm{source}},{v_{1,{\kern 1pt} {\kern 1pt} s}}) \cup _{i = 1}^{m - 1}e({v_{i,{\kern 1pt} {\kern 1pt} s{\kern 1pt} }},{v_{i + 1,{\kern 1pt} {\kern 1pt} {\kern 1pt} s{\kern 1pt} }}) \cup e({v_{m - 1,{\kern 1pt} {\kern 1pt} s{\kern 1pt} {\kern 1pt} }},{\rm{sink}})$

式中，数据项边$e({v_{{\kern 1pt} i,{\kern 1pt} {\kern 1pt} s{\kern 1pt} }},{v_{i + 1{\kern 1pt} ,{\kern 1pt} {\kern 1pt} s{\kern 1pt} }})$对应的权重为$w(e({v_{{\kern 1pt} i,{\kern 1pt} {\kern 1pt} s{\kern 1pt} }},{v_{i + 1{\kern 1pt} ,{\kern 1pt} {\kern 1pt} s{\kern 1pt} }})) = - {\omega _{i + 1}}\log (f({I_0}(s)\left| {{\Theta _{i + 1}}} \right.))$，它刻画了特征s到第i+1类的区域相似性，起到相似性特征聚集的作用。而边缘项边集可表示为：

$\begin{array}{c} {\rm{Edg}}{{\rm{e}}_{\_N}} = \{ e({v_{i,s}},{v_{i,r}})|i \in \left\{ {1,2, \cdots ,m - 1} \right\},{\kern 1pt} {\kern 1pt} \\ {\kern 1pt} s,r \in M \cap {\kern 1pt} {\kern 1pt} r \in {N_Q}(s)\} \end{array}$

式中，边缘项边$e({v_{{\kern 1pt} i,{\kern 1pt} {\kern 1pt} s{\kern 1pt} }},{v_{i + 1{\kern 1pt} ,{\kern 1pt} {\kern 1pt} s{\kern 1pt} }})$的权重$w(e({v_{i,s}},{v_{i{\kern 1pt} ,r}})) = 2\gamma {w_{s,r}}(1 - \delta (\Upsilon (s) - \Upsilon (r)))g(\left| {\nabla C(s)} \right|)$，它引入了邻域空间约束关系及边界的检测，在分割过程中它可吸引活动轮廓尽量向目标的真实边界移动。图 2a提供了4类非同质目标区域的自然图像，通过前面t-link边集Edge__T与n-link边集Edge__N的描述，构建了3层带权有向图G，如图 2b所示。其中，层与层之间较粗的t-link边表示较大的相似度，反之亦然。位于同一图层中的n-link边表示空间领域内像素间的空间约束关系，当局部相邻特征点之间具有较大的差异时，则被分配到相同类别的可能性较小，于是n-link边较细。式(7)中的多类能量函数可按上述思想，构建得到对应的多层图，利用最大流/最小割^[15]方法可求得全局近似最优解。

图2

图2 3×3的自然图像构建的3层图

在多类变分活动轮廓演化的过程中，图像区域被活动轮廓分割为多个标签区域，为了及时将分割后的标签区域作用于区域${\Omega _k}$内部，需对多变量学生-t分布及时进行相关统计参数的更新。由于样本之间相互独立，t次迭代时整个图像的m类最大期望/最大似然能量为：

$\begin{array}{c} E_{{\rm{ML}}}^{(t)} = \mathop {{\rm{Max}}}\limits_{{\Theta ^{(t)}}} \left\{ {\sum\limits_{s \in S} {( - \delta (\Upsilon (s) - k){\kern 1pt} } } \right. \times \\ {\kern 1pt} \log (\sum\limits_{k = 1}^m {\omega _k^{(t)}f({I_0}(s)\left| {\Theta _k^{(t)}} \right.)} )) + \left. {\lambda (\sum\limits_{k = 1}^K {\omega _k^{(t)}} - 1)} \right\} \end{array}$

式中，第k类标签区域对应的相关统计参数$\omega {\kern 1pt} _{_k}^{(t)}$和$\Theta {\kern 1pt} _{_k}^{(t)} = \left\{ {\nu _{_k}^{(t)},\mu _{_k}^{(t)},\Sigma {\kern 1pt} _{_k}^{(t)}} \right\}$，可通过偏导数${\partial E_{\text{ML}}^{(t)}}/{\partial }\;\omega _{_{k}}^{(t)}$、${\partial E_{\text{ML}}^{(t)}}/{\partial }\;\nu _{_{k}}^{(t)}$、${\partial E_{\text{ML}}^{(t)}}/{\partial }\;\mu _{_{k}}^{(t)}$、${\partial E_{\text{ML}}^{(t)}}/{\partial }\;\Sigma _{_{k}}^{(t)}$求得。首先，引入中间变量$p{\kern 1pt} _{_{s,k}}^{(t)}$及$u{\kern 1pt} _{_{s,k}}^{(t)}$：

$p_{s,k}^{(t)} = \frac{{\omega _k^{(t - 1)}f({I_0}(s)\left| {\Theta _k^{(t - 1)}} \right.)}}{{\sum\limits_{k = 1}^K {\omega _k^{(t - 1)}f} ({I_0}(s)\left| {\Theta _k^{(t - 1)}} \right.)}}$

$u_{s,k}^{(t)} = \frac{{\nu _k^{(t - 1)} + D}}{{\nu _k^{(t - 1)} + {{({I_0}(s) - \mu _k^{(t - 1)})}^{\rm{T}}}\Sigma {{_k^{(t - 1)}}^{^{ - 1}}}({I_0}(s) - \mu _k^{(t - 1)})}}$

由于上面相关参数的偏导数等于0，可得：

$\omega _k^{(t)} = \frac{{\sum\limits_{s \in S} {\delta (\Upsilon (s) - k)p_{s,k}^{(t)}} }}{{\sum\limits_{s \in S} {\sum\limits_{g = 1}^m {\delta (\Upsilon (s) - g)p_{s,g}^{(t)}} } }}$

$\mu _k^{(t)} = \frac{{\sum\limits_{s \in S} {\delta (\Upsilon (s) - k)p_{s,k}^{(t)}u_{s,k}^{(t)}{I_0}(s)} }}{{\sum\limits_{s \in S} {\delta (\Upsilon (s) - k)p_{s,k}^{(t)}} u_{s,k}^{(t)}}}$

$\begin{array}{l} \Sigma _k^{(t)} = \\ \frac{{\sum\limits_{s \in S} {\delta (\Upsilon (s) - k)p_{s,k}^{(t)}u_{s,k}^{(t)}({I_0}(s) - \mu _k^{(t)}){{({I_0}(s) - \mu _k^{(t)})}^{\rm{T}}}} }}{{\sum\limits_{s \in S} {\delta (\Upsilon (s) - k)p_{s,k}^{(t)}u_{s,k}^{(t)}} }} \end{array}$

且自由度$\nu {\kern 1pt} _{_k}^{(t)}$满足等式：

$\begin{array}{c} \sum\limits_{s{\kern 1pt} \in S} {\delta (\Upsilon (s) - k)p_{s,k}^{(t)}\left[ {\Psi \left( {\frac{{\nu _k^{(t - 1)} + {\rm{D}}}}{2}} \right) - \Psi \left( {\frac{{\nu _k^{(t - 1)}}}{2}} \right) + } \right.} {\kern 1pt} \\ {\kern 1pt} \log \left( {\frac{2}{{\nu _k^{(t - 1)} + ({I_0}(s) - \mu _k^{(t - 1)})\Sigma _k^{{{(t - 1)}^{ - 1}}}({I_0}(s) - \mu _k^{(t - 1)})}}} \right) + \\ \left. {\log \left( {\frac{{\nu _k^{(t)}}}{2}} \right) - u_{s,k}^{(t)} - \Psi \left( {\frac{{\nu _k^{(t)}}}{2}} \right){\rm{ + }}1} \right]{\rm{ = }}0 \end{array}$

式中，\$\Psi ( \cdot )$为双伽玛函数，即$\Psi (x) = \frac{{\partial \ln (\Gamma (x))}}{{\partial x}}$。

自由度调控参数$\nu {\kern 1pt} _{_k}^{(t)}$可通过牛顿-拉夫逊迭代方式计算得到。为了方便、快捷地执行本文提出的多类变分图像分割方法，具体的算法执行过程如下：

1) t=0，利用CEM³ST算法^[10]初始化活动轮廓区域${\Omega _k}$的总数m，以及各个区域内部对应的初始统计参数$\omega {\kern 1pt} _{_k}^{(t)}$和$\Theta {\kern 1pt} _{_k}^{(t)} = \left\{ {\nu _{_k}^{(t)},\mu _{_k}^{(t)},\Sigma {\kern 1pt} _{_k}^{(t)}} \right\}$。

2) 将式(7)对应的多类变分活动轮廓能量转化为对应的多层图G。

3) 利用最大流/最小割算法^[19]进行多层图割优化，令t=t+1，并根据分割后标签区域${\Omega _k}(k = 1,2, \cdots ,m)$更新各个标签区域对应的相关统计参数$\omega {\kern 1pt} _{_k}^{(t)}$和$\Theta {\kern 1pt} _{_k}^{(t)} = \left\{ {\nu _{_k}^{(t)},\mu _{_k}^{(t)},\Sigma {\kern 1pt} _{_k}^{(t)}} \right\}$，并重新计算相似度$p_{s,k}^{(t)}$。

4) 根据活动轮廓区域${\Omega _k}(k = 1,2, \cdots ,m)$的整体区域变化，确定分割过程是否稳定，如果稳定，则退出；否则，转到步骤2)。

为了合理地进行实验评估与量化分析，将采用改进的多段恒定变分能量模型(MMPC-ACM)^[6]方法以及多类彩色纹理图像分割方法(CEM³ST)^[10]与本文的方法进行实验对比分析。在MMPC-ACM方法中，作者假定每个区域$\Omega {\kern 1pt} {{\kern 1pt} _l}$内部采用相同数量的恒定密度中心描述，即$n{{\kern 1pt} _l} = 5$。在CEM³ST方法中，平滑因子$\lambda $与除噪常数$\tau $分别为10和5。在本文方法中，边缘外力项调节因子$\gamma $设置为5。为了对本文提出的方法进行有效地验证，将采用具有非线性密度变化的自然图像进行实验对比与量化分析。

图 3a提供了自然图像，它包含多个非同质目标区域，且不同区域间具有明显的边界差异性。如天空、草地、岩石、山林等。图 3中第1列、第2列、第3列分别表示分割边界、标签均值图、分割结果图。由图 3c提供的分割区域与分割边界可见，MMPC- ACM方法虽能将多个非同质目标区域分割出来，但在目标区域的边界处容易出现边界模糊及误分割现象，此外，部分边界出现重叠。这种分割结果出现的原因在于MMPC-ACM方法采用恒定聚类中心来描述每个区域内部的概率密度，容易出现过拟合或欠拟合现象。当区域间的边界具有多样变化的密度分布时，MMPC-ACM方法容易分割出很多无意义的边界，影响最终分割的整体效果。而CEM³ST方法与本文提出的方法能够准确的将多个非同质目标区域分割开来。相比而言，本文的方法比CEM³ST方法捕获的边缘更加准确，且最终分割的边界更加光顺。

图3

图3 自然场景图像的分割边界、标签均值图以及分割结果对比

为了进一步验证本文方法对概率密度的非线性及连续性描述，图 4提供了猕猴图像，它包含的特征信息在特征空间具有某种非线性变化，图中第1列、第2列、第3列分别表示分割边界、标签均值图、分割结果图。图 4b与图 4d提供了利用自由度调控的多变量混合学生-t分布进行密度描述的分割结果，其能够更加准确的将同质目标区域分割出来，虽然图 4d具有较好的目标区域整体性，但是对于阴影区域与树叶区域它难以进行有效的区分，而本文方法能够很好地将猕猴区域、树叶区域以及阴影区域分割开来。且分割效果更加接近于地面真实分割。而MMPC-ACM方法采用5个恒定聚类中心来描述每个区域$\Omega {{\kern 1pt} _l}$，其分割结果如图 4c所示，出现了很多离散的、零碎的小区域，且部分目标同质区域被细分，因此它丢失了同质目标区域的非线性刻画。此外，由分割结果可见，利用多段恒定聚类中心难以有效描述具有复杂性、非线性的猕猴区域与植物区域，它容易将植物区域分割为多个零碎区域，且对猕猴嘴部的小阴影区域比较敏感。而采用本文提出的自由度调控多变量混合学生-t分布描述，它不仅具有较好的非线性描述能力，而且能够有效地刻画目标区域间的非线性密度变化。因此，对于大多数具有非线性密度分布的图像而言，利用自由度调控的多变量混合学生-t分布描述更加适合，它具有更强的非线性描述能力。

图4

图4 具有非线性的自然图像分割结果对比

虽然对本文方法在自然图像上进行了实验对比与结果分析，但都是基于整体分割性能的描述。而关于本文所提出的区域内力项与区域外力项(数据项与边缘)的作用依然不是很清晰。为了合理、有效地分析本文提出的测地线区域外力的边界检测能力，图 5采用了区域间具有复杂跳变的边界以及区域差异性较大的图像进行测试。由未添加测地线项的图 5c分割结果可见，玫瑰花与摩天大楼被分割为破碎的区域，且出现了很多无意义的目标区域，这将严重影响图像分割结果在高层语意场景理解与视觉分析等方面的应用。此外，不同目标区域间出现了不连续、陡然跳变、尖锐的边界。这种结果出现的原因在于未添加边界检查项时，即本文提出的模型退化为普通的多类CV模型，它缺乏边界的检测以及物理空间的局部约束。与之相反，它不仅可以吸引活动轮廓朝着法线方向移动，而且可以通过区域内力的作用，保证同质目标区域内部的特征进行聚集，共同分辨出不同目标间的区域边界，如图 5b所示，活动轮廓朝着法线方向移动，而且可以通过区域内力的作用，保证同质目标区域内部的特征进行聚集，共同分辨出不同目标间的区域边界。

图5

图5 带边缘检测项与不带边缘检测项的分割结果对比

为了客观评价3种对比方法的有效性，本文采用文献[20提出的概率随机检索PRI(probabilistic rand index)进行准确率量化计算，PRI的取值在0～1之间。较大的量化准确率值反映实验分割的结果更加接近于真实的人工分割结果。图 6提供了由伯克利自然图像库随机选取60张自然图像的量化结果。如图 6a所示，本文方法计算的PRI准确率值高于CEM³ST方法与MMPC-ACM方法，且准确率值以58%的比例集中在0.8～1.0之间分布(图 6b所示)，而CEM³ST方法与MMPC-ACM方法分别为45%和35%。此外，表 1提供了3种对比方法的PRI平均均值与平均运行时间。本文方法的PRI平均值达到0.805，而CEM³ST方法与MMPC-ACM方法分别为0.974和0.789，本文方法的平均准确率值较高。通过平均分割时间消耗可见，本文方法对应的分割速度更快。这些量化指标值可进一步说明，本文方法的整体性能要优于CEM³ST方法与MMPC-ACM方法。

图6

图6 按升序排列的PRI准确率值及对应的区间比例分布

表1

表1 3种方法的PRI均值与运行时间对比

表1 3种方法的PRI均值与运行时间对比 本文方法 MMPC-ACM方法 CEM3ST方法 PRI均值 0.805 0.789 0.794 平均时间/s 116 229 126

本文将多类变分活动轮廓模型与多层图割模型相结合，提出了一种新的自然图像分割方法。为了提高图像特征的非线性与连续性刻画，将自由度调节的多变量学生-t分布引入到多类CV变分模型中，它突破了恒定聚类中心假设。同时，针对活动轮廓模型缺乏区域外力的缺陷，引入了测地线区域外力约束项，它不仅能够有效分割出区域间的光滑边界，而且避免了分割出离散、零碎的噪声区域。对于本文提出的多类变分模型的最小化问题是NP难问题，通过离散化多类变分活动轮廓模型，可将能量函数的最小化问题转化为多层图割模型的最大流/最小割优化问题，并可求得全局近似最优解。此外，对于本文提出的分割方法进行了合理的实验验证，通过对比分析表明，本文提出的分割方法不仅具有较高的分割准确率、光滑的边界，而且最终分割的区域视觉效果较好。