一种新的基于二值水平集和形态学的局部分割方法

郑 强 董恩清

(山东大学(威海)机电与信息工程学院 威海 264209)

1 引言

图像分割是图像处理领域中的热点问题之一。主动轮廓模型(Active Contour Model, ACM)[1]以其强大的数学理论基础和基于水平集(level set)的高效数值方案吸引了众多学者。主动轮廓模型分为两类：基于边缘信息的模型[2,3]和基于区域信息的模型[4,5]，这两类模型的典型代表分别是测地主动轮廓(Geodesic Active Contours, GAC)模型[2]和 CV(Chan-Vese)模型[4]。

在实际应用中，经常只需分割感兴趣的特定区域，即局部分割，如医学实践中常常需要对具有解剖意义的壳核、小脑等特定组织进行分割。基于边缘信息的模型虽能实现局部分割，但对噪声敏感，对边缘较弱图像的分割能力较差；基于区域信息的模型虽对噪声、弱边缘等问题具有较强的鲁棒性，但其模型是一种全局分割模型，会将图像中具有相近亮度的区域都分割出来，无法实现局部分割。

特别强调的是，本文中局部分割的概念与基于局部统计信息的分割模型是不同的。局部分割(local segmentation)指的是从复杂图像中分割局部区域，而基于局部统计信息的分割模型(local region-based model)是一种全局分割(global segmentation)模型。两者之间的相同点是，本文的局部分割模型和基于局部统计信息的分割模型都会用到局部统计信息。

局部分割是一个具有挑战性的课题，局部分割精度(Local Segmentation Precision, LSP)不足和窄带控制不稳定是当前局部分割方法中存在的两个主要问题。近年来，局部分割逐渐引起众多学者的重视[6-11]，其中 Lankton等人[6]的研究最具代表性。他们提出一种局部分割框架 LRBAC (Localizing Region Based Active Contours)，该框架允许任意基于区域信息的模型以局部的形式表达，将曲线进化过程控制在零水平集附近的窄带范围内。实验结果表明，该框架是一种非常有效的局部分割方法，对灰度不均匀图像也能取得较好的分割结果。然而，该方法存在以下不足：(1)水平集函数需要不断重新初始化为符号距离函数作为保证，并且为了保证算法的稳定性，进化步长必须满足 CFL(Courant-Friedrichs-Lewy)条件，这些将导致该方法的计算量较大；(2)水平集函数重新初始化过程中存在误差导致窄带控制方案不稳定。

基于LRBAC方法，Mille[7]提出一种窄带主动轮廓模型，该模型能够更好地控制窄带，并且构建了基于曲率的模型；胡玉晖[8]提出一种基于局部子区域的活动轮廓模型，较好地解决了复杂多纹理图像中的目标分割问题。然而，两者都没有解决LRBAC方法[6]中符号距离函数所带来的困扰。Zhang等人[9]提出一种可以选择局部分割或是全局分割的图像分割方法 SBGFRLS(Selective Binary and Gaussian Filtering Regularized Level Set)方法，该方法用可选择的二值水平集函数实现局部分割，并且用高斯滤波实现曲线平滑。然而，该方法对灰度不均匀图像分割能力不足，且局部分割精度受高斯滤波影响较大。

本文在LRBAC方法的基础上，提出一种新的局部分割方法 BSMORLS(Binary and Selective Morphological Operation Regularized Level Set)。主要思路如下：(1)用二值水平集取代传统的符号距离函数以降低计算量，并通过严格保持水平集函数的二值性以确保窄带控制的稳定性和局部分割的精度；(2)引入局部统计信息，增强对灰度不均匀图像的分割能力；(3)用可选择的形态学算子改进曲线平滑项，增加曲线平滑的灵活性；(4)本文方法用稀疏场算法实现以提高计算效率。在合成图像和医学图像上的实验结果表明，提出的方法能更好地实现图像局部分割。

2 局部分割模型(LRBAC)

文献[6]提出的局部分割框架LRBAC是图像局部分割模型的典型代表。设图像I的定义域为Ωx，变量x表示定义域Ωx上的任意一点，即x∈Ωx。C是用符号距离函数的零水平集表示的闭合曲线，即C= {x|φ(x)=0 }。用海维塞德函数(Heaviside Function)指定闭合曲线C的内部部分：

同理，闭合曲线C的外部部分用(1 -H(φ(x)))表示,并用海维塞德函数的导数δ(φ(x))来指定围绕在进化曲线附近的窄带区域：

为了引入变量y，定义如下特征函数

当y在以点x为圆心，r为半径的圆的内部时，其取值为1，否则为0。Ωy为变量y的定义域，即y∈Ωy。函数B(x,y)用来标记局部邻域，且被闭合曲线C分为局部内部区域和局部外部区域。

如图1所示，取δ(φ(x))控制的窄带区域上任意一点为圆心，r为半径做圆，则圆被闭合曲线C分为局部内部区域(图1闭合圆中的“1”标记部分)和局部外部区域(图1闭合圆的“2”标记部分)。

因此，LRBAC方法的局部分割模型如下：

函数F(I(y),φ(y))为任意基于区域信息的能量模型，文献[6]分析了 3种能量模型：UM(Uniform Modeling)能量模型，MS(Mean Separation)能量模型和HS(Histogram Separation)能量模型。本文以UM能量模型中的典型代表CV模型为例，式(4)的梯度下降流为

图1 局部化模型示意图

3 本文方法及相关问题分析

本文针对局部分割中存在的两个问题，引入二值水平集函数取代传统的符号距离函数，并用二值水平集函数的取代式(5)中δ(φ(x))以精确控制围绕在进化曲线附近的窄带区域，并且在每次迭代时都严格保证水平集函数的二值性。与LRBAC[6]方法相同，本文引入局部统计信息，增加对灰度不均匀图像的分割能力。同时，本文提出用可选择的形态学算子(开运算、闭运算和中值集算子)作为新的曲线平滑方案，增加曲线平滑的灵活性。为了提高计算效率，本文方法采用稀疏场算法实现。特别说明的是，二值水平集的用中心差分实现，后文中若无特别说明均以此为准。

另外，Zhao等人[12]也曾将曲线进化公式中的δ(φ(x))用取代，取代的目的是为了扩大曲线进化范围。因为当进化曲线离图像边缘较远时，δ(φ(x))函数会抑制离进化曲线较远的图像边缘的检测，而符号距离函数满足，可以消除δ(φ(x))函数对非零水平集的抑制。这与本文方案是有本质区别的。在本文方案中，二值水平集函数的并没有改变δ(φ(x))函数窄带控制的原始思想，采用二值水平集函数的取代δ(φ(x))函数是为了使得窄带控制更加稳定和精确，从而更加有利于实现图像局部分割。

3.1 本文方法(BSMORLS)

综上，本文提出的图像局部分割方法如下：

(1)初始化水平集函数φ(x)

其中Ω0是图像定义域Ω的一个子集。

其中ux和vx用式(7)和式(8)更新，用中心差分实现。值得注意的是，式(10)除了将式(5)中的δ(φ(x))函数用函数取代以更加准确而稳定的选取窄带内的点以外，还将式(5)中的δCV(φ(y))函数去掉。

Chan等人[13]指出，由于CV模型用非紧支撑且光滑的函数逼近海维塞德函数HCV(φ(y))，因此，CV模型的梯度下降流

与式(12)具有相同的稳态解

因此，将式(5)中的δCV(φ(y))函数去掉是合理的。

(3)对步骤(2)中更新的窄带内的点取如下运算

(4)用可选择的形态学运算平滑曲线，主要包括开运算、闭运算和中值集运算。

用可选择的形态学运算平滑曲线不但更加有利于保持水平集函数的二值性，并且能够增加曲线平滑的灵活性，使其更加有利于实现图像局部分割。

(5)判断进化曲线是否收敛，若没有收敛，则转到步骤(2)，若收敛，则停止曲线进化。

为提高运算速度，本文方法采用稀疏场算法实现。稀疏场算法一般用数据链表的形式实现，本文提出用二值水平集函数的作为新的数据链表。如图 2,Lout和Lin是由控制的窄带的非零数值区域，数据宽度是2ρ。以闭合曲线为分界线，Lout和Lin分别为数据链表的外部部分和内部部分。一般假定Lout和Lin的宽度相同，且都为ρ。本文算法在实现过程中只提取控制的窄带之内的点(Lin和Lout)进行更新，此举提高了计算效率。

图2 数据链表及局部分割精度示意图

3.2 相关问题分析

本文主要对局部分割精度、窄带控制的稳定性、曲线平滑方案等问题进行分析。

(1)局部分割精度及窄带控制稳定性分析

定义局部分割精度(Local Segmentation Precision, LSP)ρ是曲线法线方向上参与水平集函数进化，且对曲线进化有贡献的数据宽度。如图 2所示，参与水平集函数进化的数据宽度是2ρ，但是由于曲线上的点只能沿一个方向进化，即向内或者向外，所以实际对曲线进化有贡献的数据宽度是ρ。数据宽度越宽，LSP越低；数据宽度越窄，LSP越高。

LSP对局部分割有重要意义。如图3所示为两目标由远及近的情况下的局部分割结果，其中方形闭合曲线为初始轮廓。图3(a)和图3(b)都采用较低的LSP；图3(c)和图3(d)都采用较高的LSP。实验结果表明，当两目标距离较远时，LSP对局部分割影响较小；当两目标距离较近时，若LSP较低，则两目标之间会互相干扰，导致无法实现局部分割，只有当LSP较高时，才能较好的实现局部分割。

图3 局部分割精度分析

结合以上定义，分析LRBAC方法中的δ(φ(x))函数，SBGFRLS方法中的函数与本文BSMORLS方法的函数发现，三者控制的都是 LSP。δ(φ(x))函数建立在符号距离函数的基础上，函数建立在二值水平集函数的基础上。为方便图示，假定控制 LSP的函数为L(φ(x))，则L(φ(x))可 能 为δ(φ(x)),或。分析图4发现，δ(φ(x))函数是一平滑钟形函数(ε=1.5)，而(高斯滤波器大小为3×3)和函数类似于理想的低通滤波器。显然，只要能够保持水平集函数的二值性，比δ(φ(x))能更精确更稳定地控制窄带区域，且比函数的 LSP更高。根据 LSP的定义，若用中心差分实现，则本文方案的LSP可以稳定地控制在一个像素的精度。一个像素宽度是较为理想的LSP，如果LSP小于一个像素宽度则会影响进化速度，而如果LSP大于一个像素宽度，对于图像中存在相距较近目标的情况下，图像局部分割结果则有可能出现目标之间互相干扰的现象。

图4 局部分割精度控制函数

在SBGFRLS方法中，在使用3×3大小的高斯滤波器的情况下，LSP不能达到一个像素的精度。下面特别说明，当使用更小尺寸如2×2大小的高斯滤波器的情况下，LSP也不能稳定地达到一个像素的精度，从而证明不论高斯滤波器的大小如何，SBGFRLS方法中LSP都不能稳定地达到一个像素的精度。设定φ(x)为如下二值水平集，且考虑零水平集处于竖直状态下的情况，G(x)是大小为2×2方差为1.0的高斯滤波器。

二者卷积后得

若用中心差分由φnew(x)求得则LSP必大于 1。因此，不论高斯滤波器的大小如何，SBGFRLS方法中LSP都不能稳定地达到一个像素的精度。

(2)曲线平滑方案分析 本文提出用可选择的形态学算子(开运算、闭运算和中值集运算)作为曲线平滑方案基于以下原因：

(a)形态学运算不但具有较好的平滑作用，而且更容易保持水平集函数的二值性。LRBAC[6]采用曲率作为平滑方案，SBGFRLS[9]采用高斯滤波作为平滑方案，都与形态学中值集运算的平滑效果基本相同，但是形态学算子在运算前后更加容易保持水平集函数的二值性。在本文方法中严格保持水平集函数的二值性是一个核心问题。

(b)用可选择形态学算子作为曲线平滑方案是基于图像特征的考虑。开运算只磨光凸向图像外部的边角，而能够保留凸向图像内部的边角。闭运算只磨光凸向图像内部的边角，而能够保留凸向图像外部的边角。对于没有特别明显或者特别尖锐边角的图像，则可以用中值集运算平滑曲线。这极大增加了曲线平滑方案的灵活性。

图5所示是曲线平滑方案对比分析图，其中方形闭合曲线为初始轮廓。图5(a)和图5(b)分别用形态学中值集算子和开运算作为平滑方案的分割结果；图5(c)和图5(d)分别用形态学中值集算子和闭运算作为平滑方案的分割结果。实验结果表明，开运算能保留凸向图像内部的边角，闭运算能保留凸向图像外部的边角，增加曲线平滑方案的灵活性将更加有利于实现图像局部分割。

4 实验结果及对比分析

由于 SBGFRL[9]方法基于全局统计信息，对灰度不均匀图像的分割能力不足，为了增加可比性，本文将该方法的数据驱动力改为与本文相同，在相同数据驱动力的前提下更能证明本文方案在局部分割中的优势。此处特别标记，原文献[9]中的方法为SBGFRL(G)，修改数据驱动力后的方法为SBGFRL(L), L和G分别指基于局部的统计信息和基于全局的统计信息。

为了分析和验证本文提出的局部分割算法BSMORLS的实际效果，试验中同时采用了LRBAC[6]方法、SBGFRLS(G)[9]方法和 SBGFRLS(L)方法一同对灰度不均匀图像进行局部分割，其中方形闭合曲线为初始轮廓。

图6和图7是对灰度不均匀图像的分割结果。图6(a)为LRBAC 方法的分割结果，局部邻域半径为15；图6(b)为SBGFRLS(G)方法的分割结果，高斯滤波器大小为5×5，方差为 1.0；图 6(c)为SBGFRLS(L)方法的分割结果，局部邻域半径为15，高斯滤波器大小为5×5，方差为 1.0；图 6(d)为本文BSMORLS方法的分割结果，局部邻域半径为15，曲线平滑方案选为形态学闭运算。实验结果表明，LRBAC 方法受窄带控制不稳定影响，分割过程中出现窄带范围之内搜索不到点参与曲线进化，从而迫使曲线进化停止的现象；SBGFRLS(G)和SBGFRLS(L)方法由于LSP不够，相邻目标之间会出现互相干扰的现象，此外，对比图6(b)和图6(c)发现，SBGFRLS(G)对灰度不均匀图像的分割能力不足；本文提出的BSMORLS方法由于LSP高，对灰度不均匀图像的分割能力较强，曲线平滑方案选取灵活，因此能够达到较好的分割效果。

图5 曲线平滑方案对比分析

图6 灰度不均匀图像分割结果

图7 含有壳核区域的整个脑部切片图及解剖结构

图7(a)为取自横轴位质子和 T2加权快速自旋双回波序列(TSE-PD/T2)中的含有壳核区域的一张脑部切片图，磁共振图像来自于威海市市立医院；图7(b)为相应的解剖结构。在MR脑图像中，壳核的分割受壳核外侧的屏状核及脑岛影响较大，因为壳核外侧与屏状核之间仅有一层很薄的外囊相隔，距离非常近，在MR脑图像中一般只有一到两个像素的宽度，而脑岛和屏状核之间距离也非常近，所以在分割壳核区域时分割结果常会受到屏状核甚至脑岛区域的干扰，造成误分割。另外，尾状核对壳核的分割也有一定的干扰。图8给出相应实验结果。

图8为对几种算法对壳核区域的分割结果对比图。图8(a)为LRBAC 方法的分割结果，局部邻域半径为11；图8(b)为SBGFRLS(G)方法的分割结果，高斯滤波器大小3×3，方差 1.0；图 8(c)为SBGFRLS(L)方法分割结果，邻域半径 11，高斯滤波器大小3×3，方差 1.0；图 8(d)为本文BSMORLS方法的分割结果，邻域半径11，曲线平滑方案用形态学闭运算。实验结果表明，LRBAC 方法受窄带不稳定影响，分割过程中出现窄带范围之内搜索不到点参与曲线进化，从而迫使曲线进化停止的现象；SBGFRLS(G)和SBGFRLS(L)方法由于LSP不够，尾状核、屏状核和脑岛对壳核的分割都存在干扰现象；本文提出的 BSMORLS方法由于LSP高，对灰度不均匀图像的分割能力较强，曲线平滑方案选取灵活，因此能够达到较好的分割效果。

此外，在MR脑图像中，壳核上半部分成像比较清晰，下半部分靠近苍白球的一侧成像也比较清晰，而下半部分靠近屏状核的区域比较模糊，在分割过程中闭合曲线容易出现凸向外部的边角，因此分割壳核区域时曲线平滑方案采用形态学闭运算将更加有利于局部分割。图9给出相应的对比分析。

图9为不同曲线平滑方案对分割影响对比图。图 9(a)-9(d)都用本文方法进行分割，局部邻域半径为11，但每个图所对应的曲线平滑方案不同。图9(a)为无任何曲线平滑方案时分割结果；图 9(b)为用形态学中值集运算作为曲线平滑方案时分割结果；图 9(c)为用形态学开运算作为曲线平滑方案时分割结果；图9(d)为用形态学闭运算作为曲线平滑方案时分割结果。实验结果表明，针对壳核区域的分割，选用形态学闭运算更加有利于分割。由此可见，灵活选取曲线平滑方案将更加有利于局部分割。

5 结束语

图像分割是图像处理中的研究热点，当前的分割方法大多热衷于全局分割，而在实际应用尤其是医学图像中却更多地要求局部分割。局部分割是一个具有挑战性的课题，窄带控制不稳定和局部分割精度不足是当前局部分割方法中存在的两个主要问题。本文提出一种新的局部分割方法，能够较好地解决以上问题，并且进一步增加了曲线平滑的灵活性，稀疏场运算的引入使得本文方法实现更加简单。应用合成图像和医学图像对几种方法进行对比分析，本文方法能更好地实现图像局部分割。

图8 几种算法对壳核区域分割对比图

图9 不同曲线平滑方案对分割影响对比图

[1]Jayadevappa D, Kumar S S, and Murty D S. Medical image segmentation algorithms using deformable models: a review[J].IETE Technical Review, 2011, 28(3): 248-255.

[2]Caselles V, Kimmel R, and Sapiro G. Geodesic active contours[J].International Journal of Computer Vision, 1997,22(1): 61-79.

[3]Li C M, Xu C Y, Gui C F,et al.. Distance regularized level set evolution and its application on image segmentation[J].IEEE Transactions on Image Processing, 2010, 19(12):3243-3254.

[4]Chan T and Vese L. Active contour without edges[J].IEEE Transactions on Image Processing, 2001, 10(2): 266-277.

[5]Wang X F, Huang D S, and Xu H. An efficient local Chan-Vese model for image segmentation[J].Pattern Recognition, 2010, 43(3): 603-618.

[6]Lankton S and Tannenbaum A. Localizing region-based active contours[J].IEEE Transactions on Image Processing,2008, 17(11): 2029-2039.

[7]Mille J. Narrow band region-based active contours and surfaces for 2D and 3D segmentation[J].Computer Vision and Image Understanding, 2009, 113(9): 946-965.

[8]胡玉晖. 基于局部子区域的活动轮廓图像分割方法[J]. 通信技术, 2010, 43(2): 68-70.Hu Yu-hui. Sub-localizing region-based active contours[J].Communications Technology, 2010, 43(2): 68-70.

[9]Zhang K H, Zhang L, Song H H,et al.. Active contours with selective local or global segmentation: a new formulation and level set method[J].Image and Vision Computing, 2010,28(4): 668-676.

[10]Shi Y H, Dong E Q, Li Z Z,et al.. Research on the segmentation of tiny multi-target in brain tissues based on support vector machines[C]. IEEE International Conference on Complex Medical Engineering, Harbin, China, May 22-25,2011: 478-482.

[11]刘晓军, 董恩清, 吕成林, 等. 磁共振脑部图像壳核区域有监督蚁群分割算法研究[J]. 中国生物医学工程学报, 2010, 29(5):683-687.Liu Xiao-jun, Dong En-qing, Lv Cheng-lin,et al.. Supervised segmentation algorithm based on ant colony for putamen region in brain MRI[J].Chinese Journal of Biomedical Engineering, 2010, 29(5): 683-687.

[12]Zhao H K, Chan T, Merriman B,et al.. A variational level set approach to multiphase motion [J].Journal of Computational Physics, 1996, 127(1): 179-195.

[13]Chan T, Esedoglu S, and Nikolova M. Algorithms for finding global minimizers of image segmentation and denoising models [J].Journal on Applied Mathematics, 2006, 66(5):1632-1648.