图像分割是图像分析和理解的重要步骤之一。由于图像中存在复杂噪声、灰度不均匀以及低对比度的问题，使得图像分割面临诸多难题。而活动轮廓模型^[1]以其自适应性、亚像素精度^[2]等优点成为研究热点。

活动轮廓模型以光滑封闭的曲线演化恢复目标边界，根据使用约束条件的不同可以分为两大类：基于边界^[3-5]和基于区域^[6-8]的活动轮廓模型。基于边界的活动轮廓模型主要以图像的边缘梯度信息为约束条件，驱动活动曲线靠近并停止于目标边界。基于区域的活动轮廓模型则根据区域内图像统计信息，如灰度、纹理、色彩等信息为约束，引导曲线移动从而分割图像。一般来说，基于区域的活动轮廓模型考虑图像全局信息，能识别弱边界或无边界的目标，对噪声和初始位置具有鲁棒性。

C-V(Chan-Vese)^[3]模型是经典的基于区域的活动轮廓模型，它利用图像的全局信息，假定每个区域内的灰度是同质的，用二值分段常数函数来逼近图像完成分割。该模型能分割灰度均匀图像，对初始轮廓有强鲁棒性，但不能有效分割灰度非同质图像。受上述模型启发，研究人员首先提出了基于全局灰度的分段光滑(PS)活动轮廓模型^[6-7]，但该模型计算复杂度高，且由于引入非凸的能量泛函，容易陷入局部极小解。为解决上述问题，Li等人^[2]提出基于局部灰度信息的LBF(local binary fitting)模型，其核心是利用局部K-means方法对图像进行聚类分割，能够对灰度不均匀图像有较好的分割结果。受其启发，Zhang等人^[8]提出LIF(local image fitting)模型和Wang等人^[9]提出的LGF(local Gaussian fitting)模型，提高了LBF模型的计算效率。但LIF模型、LGF模型以及LBF模型只考虑了局部K-mean聚类，都对初始轮廓敏感，容易陷入局部极小值。为了克服这些缺点，研究人员提出一些混合模型^[10-14]，这些模型采用局部和全局灰度信息的加权组合，重新定义能量泛函，在一定程度上解决了LBF模型和LIF模型陷入局部极小值问题。但上述模型对含有非高斯噪声图像的分割效果并不理想。

在实际的医学图像分割中，观测图像经常受到强度不均匀现象和各种复杂噪声的影响，所以精确分割仍然是一个非常困难的任务。为此，Wang等人^[14]利用高斯函数提取图像中的局部信息，结合相关熵准则度量像素点到聚类中心的距离，提出了LCK (local correntropy-based K-means)^[15-16]模型，可以抵抗未知复杂噪声的影响，对非同质图像具有较好的分割效果。但该模型将聚类算法应用于传统的区域活动轮廓模型，采用梯度下降法求解，具有计算复杂度高、对初始轮廓敏感，收敛速度慢的缺点。

为解决上述LCK模型的缺陷，提出一种新的以LCK模型为基础的两步快速分割算法，该模型分两步进行：首先利用下采样减少数据量并使用LCK模型对采样后的粗尺度图像进行分割，快速获得粗分割结果及其相应的粗水平集函数。然后在水平集函数光滑性约束下，将粗分割结果及其对应的粗水平集函数上采样到原始图像的尺度，将上采样后的粗水平集函数作为细分割的初始值，利用LCK模型对原始图像进行精细分割。实验结果表明，本文算法能够在极少的迭代次数内准确分割非同质图像和噪声图像。

LCK模型是在基于区域的活动轮廓模型框架下加入聚类算法，利用交互熵度量像素点到聚点中心的距离。

令$\mathit{\pmb{\Omega}} \subset {{\mathbb{R}}^{2}}$表示图像区域，$\mathit{\pmb{\Omega}} \subset \mathbb{R}$表示输入图像。$\phi :\mathit{\pmb{\Omega}} \subset \mathbb{R}$表示水平集函数，分割轮廓、前景和背景区域分别表示为$\boldsymbol{C}=\left\{ x\left| \phi \left( x \right)=0 \right. \right\}$，${{\mathit{\pmb{\Omega}} }_{f}}=\left\{ x\left| \phi \left( x \right)>0 \right. \right\}$和${{\mathit{\pmb{\Omega}} }_{b}}=\left\{ x\left| \phi \left( x \right)<0 \right. \right\}$，关于$\phi $的Heaviside函数$H\left( \phi \right)$表示演化曲线的内部，导数表示为Dirac函数$\delta \left( \phi \right)$。为了计算的稳定性，通常用光滑函数^[3]近似替代Heaviside函数和Dirac函数。

$ \begin{align} & {{E}_{\text{LCK}}}\left( \phi \right)=-{{\lambda }_{1}}\iint_{\Omega }{{{K}_{\varepsilon }}\left( y-x \right){{\sigma }^{2}}{{K}_{\sigma }}\left( I\left( y \right)- \right.} \\ & \ \ \ \ \ \ \ \ \ \ \ \ \left. {{\mu }_{1}}\left( x \right) \right)H\left( \phi \left( y \right) \right)\text{d}y\text{d}x- \\ & \ \ \ \ \ \ \ \ \ \ \ \ {{\lambda }_{2}}\iint_{\Omega }{{{K}_{\varepsilon }}\left( y-x \right){{\sigma }^{2}}{{K}_{\sigma }}\left( I\left( y \right)- \right.} \\ & \ \ \ \ \ \ \ \ \ \ \ \ \left. {{\mu }_{2}}\left( x \right) \right)\left( 1-H\left( \phi \left( y \right) \right) \right)\text{d}y\text{d}x+ \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \upsilon \int_{\Omega }{\left| \nabla H\left( \phi \left( y \right) \right) \right|}\text{d}x+ \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \kappa {{\int_{\Omega }{\left( \left| \nabla \phi \left( x \right) \right|-1 \right)}}^{2}}/2\text{d}x \\ \end{align} $

(1)

式中，第1项和第2项称为外部能量项；第3项为长度约束项；第4项是正则项。

$ {{K}_{\sigma }}=\exp \left( \frac{-{{x}^{2}}}{2{{\sigma }^{2}}} \right) $

为Gauss核函数，参数$\sigma $控制着核函数的尺度大小，即

$ {{K}_{\varepsilon }}=\left( \frac{1}{2\text{ }\!\!\pi\!\!\text{ }\sigma _{d}^{2}} \right)\exp \left\{ \frac{{{\left\| y-x \right\|}^{2}}}{2\text{ }\!\!\pi\!\!\text{ }\sigma _{d}^{2}} \right\} $

是用来提取局部信息的高斯核函数，参数${{\sigma }_{d}}$控制着核函数的尺度大小。${{\mu }_{1}}$，${{\mu }_{2}}$是局部区域像素均值。${{\lambda }_{1}}$，${{\lambda }_{2}}$，$\upsilon $，$\kappa $是非负参数。

式(1)中含有高斯核函数，是非线性的，在求解最小值时，采用加权迭代(IR)算法。在IR迭代算法中，第t次迭代的相关熵距离采用均方差来衡量，即

$ \left\{ \begin{align} & {{\sigma }^{2}}{{K}_{\sigma }}\left( I\left( y \right)-{{\mu }_{1}}\left( x \right) \right)= \\ & -\omega _{1}^{t}\left( x \right){{\left\| I\left( y \right)-\mu _{1}^{t}\left( x \right) \right\|}^{2}} \\ & {{\sigma }^{2}}{{K}_{\sigma }}\left( I\left( y \right)-{{\mu }_{2}}\left( x \right) \right)= \\ & -\omega _{2}^{t}\left( x \right){{\left\| I\left( y \right)-\mu _{2}^{t}\left( x \right) \right\|}^{2}} \\ \end{align} \right. $

(2)

$ \left\{ \begin{align} & \omega _{1}^{t}\left( x \right)={{K}_{\sigma }}\left( I\left( y \right)-\mu _{1}^{t-1}\left( x \right) \right) \\ & \omega _{2}^{t}\left( x \right)={{K}_{\sigma }}\left( I\left( y \right)-\mu _{2}^{t-1}\left( x \right) \right) \\ \end{align} \right. $

(3)

y点像素权重是其周围所有点局部权重的加和，在阈值$\tau =0.1$时，定义y点像素权重为

$ \begin{align} & {{\omega }^{t}}=\int_{\Omega }{{{K}_{\sigma }}\left( I\left( y \right)-\mu _{1}^{t-1}\left( x \right) \right)H}\left( {{\phi }^{t-1}}\left( y \right) \right)\text{d}y+ \\ & \int_{\Omega }{{{K}_{\sigma }}\left( I\left( y \right)-\mu _{2}^{t-1}\left( x \right) \right)\left( 1-H\left( {{\phi }^{t-1}}\left( y \right) \right) \right)}\text{d}y \\ \end{align} $

(4)

$ 且{{w}^{t}}\left( \phi \right)=\left\{ \begin{align} & {{w}^{t}}\left( \phi \right)\ \ \ \ \ {{w}^{t}}\left( \phi \right)\ge \tau \\ & 0\ \ \ \ \ \ \ \ \ \ \ {{w}^{t}}\left( \phi \right)<\tau \\ \end{align} \right. $

(5)

使用梯度下降法最小化能量泛函${{E}_{\text{LCK}}}\left( \phi \right)$，求解过程分以下两步进行：

1)极小化能量泛函${{E}_{\text{LCK}}}\left( \phi \right)$可得到$\mu _{1}^{t}\left( \phi, x \right)$、$\mu _{2}^{t}\left( \phi, x \right)$的更新方程为

$ \mu _{1}^{t}\left( \phi ,x \right)=\frac{{{K}_{\sigma }}*\left[ {{\omega }^{t}}H\left( {{\phi }^{t-1}} \right)I\left( x \right) \right]}{{{K}_{\sigma }}*\left( {{\omega }^{t}}H\left( {{\phi }^{t-1}} \right) \right)} $

(6)

$ \mu _{2}^{t}\left( \phi ,x \right)=\frac{{{K}_{\sigma }}*\left[ {{\omega }^{t}}\left( 1-H\left( {{\phi }^{t-1}} \right) \right)I\left( x \right) \right]}{{{K}_{\sigma }}*\left( {{\omega }^{t}}\left( 1-H\left( {{\phi }^{t-1}} \right) \right) \right)} $

(7)

2)固定$\mu _{1}^{t}\left( \phi, x \right)$、$\mu _{2}^{t}\left( \phi, x \right)$，更新水平集函数通过梯度流方程，即

$ \begin{align} & \frac{\partial \phi }{\partial t}=\delta \left( \phi \right)H\left( \phi \right)\times \\ & \int_{\Omega }{{{K}_{\varepsilon }}\left( y-x \right)\left[ -{{\lambda }_{1}}{{\left( I\left( y \right)-{{\mu }_{1}}\left( x \right) \right)}^{2}} \right.}+ \\ & \ \ \ \ \ \ \ \ \ \left. {{\lambda }_{2}}{{\left( I\left( y \right)-{{\mu }_{2}}\left( x \right) \right)}^{2}} \right]\text{d}y+ \\ & \ \ \ \ \ \ \ \ \ \ \nu \delta \left( \phi \right)\text{div}\left( \nabla \phi /\left| \phi \right| \right)+ \\ & \ \ \ \ \ \ \ \ k\left( {{\nabla }^{2}}\phi -\text{div}\left( \nabla \phi /\left| \phi \right| \right) \right) \\ \end{align} $

(8)

LCK模型是能量泛函为非凸函数的局部区域活动轮廓模型，在求解最小值时对初始值敏感，即该分割模型对初始轮廓敏感。受上述启发，结合下采样方法降低数据，快速定位粗分割轮廓。提出了一种新的基于LCK模型的两步快速分割算法，既克服了LCK模型收敛速度慢，计算复杂度高的缺点，又提高了算法分割的精确度。该算法包括粗分割和细分割两个步骤。

对图像的灰度信息进行下采样可以减少数据量，但直接下采样容易导致丢失过多图像灰度信息进而使得采样的结果不具有代表性。首先将原始图像的灰度值与高斯核函数做卷积，用该像素邻域的权重均值来代替像素值得到新的图像的灰度信息，对新的图像进行下采样，既减少样本量，又避免丢失过多灰度信息，得到粗尺度图像用${{I}_{h}}:\Omega \subset \mathbb{R}$表示并代替式(1)中的I，对I_h采用LCK模型进行分割，得到粗水平集函数用${{\phi }_{h}}$表示。图 1(a)是原始图像，灰度值数据量为90×100；图 1(b)是1/2下采样粗尺度图像及迭代20次对应的粗分割结果。图中粗尺度图像的灰度值数据量为45×50。数据量减少到原始图像数据量的1/4。实验结果表明，粗分割数据量小，速度快；结果不精确，但接近真实结果。

粗分割是对下采样后粗尺度图像进行的分割，在水平集函数光滑性约束下，采用双三次样条插值将粗水平集函数${{\phi }_{h}}$及粗分割结果上采样到原始图像的尺度。${{\phi }_{H}}$用来表示上采样后的粗水平集函数。图 1(c)是粗水平集函数${{\phi }_{H}}$对应的分割结果与真实目标边界很接近，但仍有一定差距。为解决上述问题，将水平集函数${{\phi }_{H}}$作为初始水平集函数，采用LCK模型进行细分割。

细分割是对原始图像进行的分割，有较好的初始值，改进LCK模型为

$ \begin{align} & E\left( \phi \right)= \\ & -{{\lambda }_{1}}\iint_{\Omega }{{{K}_{\varepsilon }}}\left( y-x \right){{\sigma }^{2}}{{K}_{\sigma }}\left( I\left( y \right)-{{\mu }_{1}}\left( x \right) \right)H\left( \phi \left( y \right) \right)\text{d}y\text{d}x- \\ & {{\lambda }_{2}}\iint_{\Omega }{{{K}_{\varepsilon }}}\left( y-x \right){{\sigma }^{2}}{{K}_{\sigma }}\left( I\left( y \right)-{{\mu }_{2}}\left( x \right) \right)\left( 1-H\left( \phi \left( y \right) \right) \right)\text{d}y\text{d}x+ \\ & \nu \int_{\Omega }{\left| \nabla H\left( \phi \left( x \right) \right) \right|}\text{d}x+\kappa \int_{\Omega }{{{\left( \left| \nabla \phi \left( x \right) \right|-1 \right)}^{2}}}/2\text{d}x+ \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \gamma /2\int_{\Omega }{{{\left| \phi \left( x \right)-{{\phi }_{H}}\left( x \right) \right|}^{2}}\text{d}x} \\ \end{align} $

(9)

式中，$\gamma $是用来控制精确度的一个正参数。最小化能量泛函$E\left( \phi \right)$，可以得到水平集函数的演化方程，即

$ \begin{align} & \frac{\partial \phi }{\partial t}=\delta \left( \phi \right)H\left( \phi \right)\int_{\Omega }{{{K}_{\varepsilon }}\left( y-x \right)\left[ -{{\lambda }_{1}}\left( I\left( y \right)- \right. \right.} \\ & \ \ \ \ \left. {{\left. {{\mu }_{1}}\left( x \right) \right)}^{2}}+{{\lambda }_{2}}\left( I\left( y \right)-{{\mu }_{2}}\left( x \right) \right)2 \right]\text{d}y+ \\ & \ \ \ \ \ \ \ \ \ \ \ \ \nu \delta \left( \phi \right)\text{div}\left( \nabla \phi /\left| \phi \right| \right)+ \\ & \ \ \ \ \ \ \ \ \ k\left( {{\nabla }^{2}}\phi -\text{div}\left( \nabla \phi /\left| \phi \right| \right) \right)+ \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \gamma \left| \phi -{{\phi }_{H}} \right| \\ \end{align} $

(10)

本文方法的具体实现步骤可归纳如下：

1)输入原始图像的灰度信息I。

2)粗分割：

(1)将I与高斯核函数卷积下采样，得到粗尺度图像的灰度信息I_h。

(2)初始化水平集函数

$ {\varphi _0}\left\{ \begin{array}{l} - \rho \;\;\;\;\;i \in {\Omega _f}\\ 0\;\;\;\;\;\;\;\;i \in C\\ \rho \;\;\;\;\;\;\;i \in {\Omega _b} \end{array} \right. $

(11)

(3)由式(4)(5)计算权重${{\omega }^{t}}\left( \phi \right)$；

(4)通过式(6)(7)计算局部聚点$\mu _{1}^{t}\left( \phi, x \right)$，$\mu _{2}^{t}\left( \phi, x \right)$，即

$ \begin{align} & \mu _{1}^{t}=\min \left\{ \mu _{1}^{t}\left( \phi, x \right), \mu _{2}^{t}\left( \phi, x \right) \right\} \\ & \mu _{2}^{t}=\max \left\{ \mu _{1}^{t}\left( \phi, x \right), \mu _{2}^{t}\left( \phi, x \right) \right\} \\ \end{align} $

(5)根据式(8)更新水平集函数

$ {\phi ^t} = {\phi ^{t - 1}} + \eta \frac{{\partial {\phi ^{t - 1}}}}{{\partial t}} $

(12)

式中，$\eta $是时间步长；

(6)检查演化方程是否稳定，输出粗水平集函数${{\phi }_{h}}$，否则返回步骤(3)。

3)细分割：

(1)将粗水平集函数${{\phi }_{h}}$采用双三次样条插值上采样到原图像比例作为此次分割的初始水平集${{\phi }_{H}}$；

(2)重复上述步骤2)中的步骤(3)(4)；

(3)根据式(10)(12)更新水平集函数；

(4)检查方程的稳定性，否则返回步骤(2)。

本节将提出的模型应用到人工合成图像和医学图像中，在以下所有实验中，对比在高斯噪声、泊松噪声和椒盐噪声下，新模型相比较与原始的LCK模型分割效果较好及迭代次数更少。每次试验在相同的迭代步长$\eta =0.1$下进行。比较时采用F-score评价指标

$ {F_{{\rm{score}}}} = \frac{{2TP}}{{2TP + FN + FP}} $

(13)

式中，TP表示算法分割后图像目标区域判定为原始图像目标区域的像素点个数，FP表示算法分割后图像目标区域未判定为原始图像目标区域的像素点个数，FN表示算法分割后非图像目标区域判定为原始图像目标区域的像素点个数，${{F}_{\text{score}}}$值越大，则分割的结果越好。实验中用r表示迭代次数，迭代次数越少，模型的收敛速度越快。

实验1分割合成图像，将提出的模型与LCK模型进行比较。表 1中分别列出了LCK模型和本文方法分割结果的${{F}_{\text{score}}}$值和迭代次数；图 2是LCK模型和本文方法对原始图像及泊松噪声图像的分割结果。对比发现，本文方法的分割效果明显优于LCK模型，尤其是对泊松噪声污染的图像。而且，本文方法迭代次数远低于LCK模型的迭代次数，明显提高了模型的收敛速度。

实验2利用两步分割方法分割一幅典型的灰度非同质T-形图像。图 3(a)给出了原始的T-形图像，图 3 (b)给出了第1步分割的结果双三次样条插值到原图形比例后应用于原始图像，图 3(c)(d)分别给出了第2步分割第一次迭代后的分割结果和第五次迭代后的最终分割结果。${{S}_{b}}/{{S}_{d}}=0.9779$，${{S}_{c}}/{{S}_{d}}=0.9864$，S_b、S_c、S_d分别表示图 3(b)(c)(d)图像分割结果目标域的面积；综上可知，第2步分割初始值和最终值差别不大且五次迭代后就达到较好的分割结果。

实验3分割灰度非同质T-形图像的高斯噪声、泊松噪声和椒盐噪声图像。实验中，令尺度参数$\sigma =5$。图 4(a)是原始图像及其在高斯噪声水平为0.02、泊松噪声和椒盐噪声水平为0.2严重污染下的图像；图 4(b)(c)表明本文方法的分割效果明显优于LCK模型。且从表 2中可以发现迭代次数在50次以内，提高了分割的效率。

实验4分割医学图像和相应的泊松噪声与椒盐图像。图 5(a)分别为原始灰度非同质图像、高斯噪声水平为0.03和泊松噪声严重污染下的图像，图 5(b)(c)分别列出了LCK模型和本文模型的分割结果。结果表明，本文模型对灰度非同质图像和噪声污染的图像均能有较好地分割效果。

根据下采样以减少数据量的原理，将下采样的思想引入到LCK模型中，得到基于LCK模型的一种快速分割算法。为了克服直接采样得到的样本，具有代表性，结合像素信息与高斯核做卷积后得到局部权重的特点，提出了基于LCK模型的两步分割算法：第1阶段：粗分割，将全部的像素信息与高斯核做卷积后下采样得到新的图像，然后用LCK模型对该图像进行分割，得到粗分割结果及其对应的粗水平集函数；第2阶段：细分割，在水平集函数光滑性约束下，将粗分割结果及其对应的粗水平集函数上采样到原始图像的尺度，然后将上采样后的粗水平集函数作为细分割的初始值，利用LCK模型对原始图像进行精细分割。这种算法保留了LCK模型对灰度非同质图像的分割效果，给模型提供了较好的初始轮廓，从而提高了分割速度。实验结果表明，本文方法具有较高的分割精度和更快的收敛速度。但是，整个算法依然是以LCK模型为主，依然无法克服其对初始轮廓敏感和尺度参数依赖明显的特点，在后续的研究中，可以将图像的全局信息结合到算法，进一步提高算法的准确率。