基于增强补种的混合水平集CT图像分割

姚倩　侯筱婷 张娟　王丽

摘 要：针对医学图像CT图像像素不均匀对图像局部分割算法影响较大的问题，提出一种基于Lagrangian粒子增强补种算法的混合水平集医学CT图像分割算法。首先，针对局部图像的非均匀性，通过在计算水平集公式前先计算Lagrangian标记粒子来重建内嵌交接界面，从而提高水平集算法的质量守恒特性；其次，针对传统粒子方法在处理界面奇异性和复杂几何相关问题上的不确定性，通过增加速度矢量和单位法向量来促进奇异点和拓扑变化点速度场的收敛；最后，通过在合成数据测试集和真实CT图像上的仿真测试表明，所提算法在边缘分割收敛精度及运算速度上均要优于对比算法。

关键词：补种算法；奇异性；水平集； CT图像；分割

中图分类号：TP391 文献标识码：A

1 引 言

图像处理方向学者Osher与Sethian等[1]联合研究并提出了水平集（Level Set， LS）分割算法，主要思路是将低维闭合二维曲线映射到高维的水平集上进行处理，从而实现了图像分割算法稳定性的增强。尽管如此，图像界面奇异性和复杂几何相关性对水平集算法影响仍然很大，如何处理好上述问题，对于提高水平集分割算法在CT、X-ray等医学成像中的应用效果，至关重要，有利于增强算法在处理对比度差、像素不均匀、边界模糊等问题并存的医学影像成像分割问题上的表现[2-3]。

针对锋利边缘和极端变形等情况，不同学者提出各种改进方式，如文献[4]基于欧拉-拉格朗日粒子的方法对水平集算法进行改进，该方法在保持分割界面质量上表现不错。此外，还有局部水平集方法[5-6]、间断Galerkin方法[7-8]、梯度增强水平集方法[9-10]等。上述方法在处理锋利边缘和极端变形等情况有一定效果，但是当像素不均匀或交接界面存在奇异点或拓扑变化时，上述算法效果不佳。

对此，本文在文献[4]基础上，采用Lagrangian粒子增强补种算法与水平集算法相结合，通过增加速度矢量和单位法向量来促进奇异点速度场的收敛，提出了一种基于Lagrangian粒子增强补种算法的混合水平集医学CT图像分割算法（LPRLS）。

2 基于粒子的水平集算法

2.1 水平集算法

2.2 混合粒子水平集方法

该方法融合了欧拉水平集算法[16]和拉格朗日粒子算法[17]，采用拉格朗日无质量粒子来纠正图像界面待求解区域的水平集函数。两组粒子被随机放置在窄带的界面附近，并吸引界面到正确的一边，其中阳性粒子趋向>0的一边，而阴性粒子趋向≤0的一边，格朗日粒子将在给定速度下产生对流：

式中，xp为粒子位置；uxp为粒子的速度向量。由于演化方程的耗散缺失，因此流动的特征信息被完全保存。然后，利用三阶TVD Runge Kutta方法来求解时间导数的演化方程。

每个粒子都有一个位置和范围半径，用来弥补因水平集函数的局限性所引起的的位置误差。每个粒子的半径由基于网格的最大值和最小值进行界定。根据公式（6）粒子半径的最大和最小界定值，可由下式进行计算：

在界面清晰像素均匀的较容易解决图像表面地区，水平集函数的解是足够准确的，并且粒子不漂移距离交接界面区域过远。然而在像素不均匀、界面模糊的不容易解决区域，水平集函数计算会产生质量损失，粒子会漂移距离交接界面过远。界面上逃脱的正（或负）粒子，基于网格点、粒子半径和粒子符号可定义局部水平集函数。水平集函数集的纠正可通过对比网格水平集函数值和逃脱粒子的局部水平集函数值实现。基于网格水平集函数的符号，则逃脱阴阳离子的局部水平集函数可分别定义如下：

通过将>0和≤0的区域分别与网格水平集函数值+和-进行比较，然后基于误差校正方法对水平集函数进行重建。上述水平集函数值+和-分别位于阳性和阴性区域，然后通过逃脱粒子计算定义的局部水平集函数值。x>0和x≤0区域的每个角落，可通过下式计算：

水平集函数值+和-可通过设置+和-的等价变量来合成为一个单一的水平集函数，该水平集都具有在每个网格点的最小值。变量可定义如下：

3 粒子补种方法及改进

3.1 粒子补种方法

在存在拉伸和撕裂的流动界面中，特别是不均匀区域，这样的区域在执行上述混合粒子水平集算法时会缺乏足够数量的粒子产生。为了保持界面的准确性，在极度变形的交接界面附近进行粒子补种以保持适当的粒子分布是很重要的步骤。对此提出一种粒子补种算法，阴阳粒子补种范围，如图2所示：

如果粒子在距离交接界面3个网格点以内的总数量大于最大值1.5Np，或者小于最小值0.5Np，式中Np取值可定义为：

Np=16， for 2D screen64， for 3D screen（15）

粒子被随机删除（或添加）直到颗粒的数目在规定的最大值和最小值之间，如图2所示。补种过程不仅加强了区域特征信息实现对运动区域的完全描述，而且消除了不良粒子诱导。

在处理相关性问题中，在交接界面附近补种适量的粒子数量是很重要的。当基于水平集函数的几何信息进行补种时，表面的数量（面积或体积）可由下式给出：

式中，Ω为操作区域；H为Heaviside辅助函数，可定义如下：

这种基于计算数或几何信息补种程序方法最为简单，但这种补种标准并不明确。如果特征被冲击波或可压缩流动等外界干扰破坏，则原有的粒子补种机制无法删除含有错误信息的污染颗粒，从而影响到水平集切面的准确性。频繁的粒子补种操作并不是解决问题的根本方案，补种操作并不能找出导致界面错误的非正常粒子，反而会删除含有正确特征的粒子。

3.2 改进粒子补种算法

原始粒子补种算法对于处理界面奇异问题和复杂的几何相关性问题时效果不佳，比如单独或联合区域划分问题，因此需要一种更为严格和可靠的粒子补种方法。endprint

速度场和几何流会引起界面的奇异性，比如角落或边缘。通过分析速度场和界面之间的相关性可知包含错误信息的特定区域。当速度场（或特征线）收敛于界面上的一点，由于拉格朗日粒子特性保持的完全性，那么逃逸粒子将会沿该奇异点堆积。这些逃逸的受污染粒子将会通过补种算法清除污染粒子前的校正阶段影响交接界面划分的准确性。由于粒子方法无法区分拓扑事件，例如，合并或分割，因此拓扑结构的变化是影响水平集交接界面准确性的另一个因素。

图3 拓扑变化示意图

本文求解思路是在处理奇异性和拓扑合并或破裂问题时（如图3所示），采用求解双曲型方程的柔性解决方案。采取的策略是：当检测到区域中含有奇异点或拓扑结构变化时，粒子的校正阶段停止执行，水平集函数值保持原始柔性计算值。奇异点是由速度场收敛而出现，并根据相应的交接界面和检测的速度向量和单位法向量的角度值进行评测，如下：

式中，θ表示交接界面附近节点特征聚集程度，如果所有的θ都接近于0，则说明此处存在奇异点，在奇异点处的特征被融合并且传统的粒子算法不能正确计算。在这里假定θ∈-10o，10o时，该处为奇异点。

在图3中，拓扑合并和分裂时发生在两界面彼此非常接近时，这种拓扑结构的变化可通过检查周围环境和水平集值进行检测。如果像素点的水平集值在0，0.1Δx范围内，且水平集值为负值，则该位置附近将发生合并或者分裂的拓扑结构变化。例如，在图3中，节点Ei，j为发生拓扑结构变化的可疑点。

综上所述，增强粒子补种算法的计算步骤如下：

Step1：（初始化）读取水平集待分割图像，并初始化算法参数；

Step2：（奇异点检测）根据3.3节检测当前区域是否属于奇异点或拓扑结构变化点，若是转Step4，否则转Step3；

Step3：（粒子补种）根据3.1节，执行粒子补种算法；

Step4：（粒子状态检测）检查环绕点Ai-1，j-1、Bi+1，j-1、Ci+1，j和Di-1，j的粒子符号状态；

Step5：（粒子删除）判断是否满足删除条件，环绕点符号是否为“-”，若是则删除节点Ei，j附近的粒子。否则转Step6；

Step6：（水平集分割）根据2.1～2.2节执行混合粒子水平集分割算法。判断是否满足终止条件，是则终止算法，并输出结果，否则转Step2，继续执行粒子的补种和删除操作。

4 仿真实验与分析

仿真设备参数：处理器为i7 2.4GHz，内存为 4G ddr1333。仿真平台选取：matlab2012a。实验对象：选取常用医学用CT图像分别为血管CT、肿瘤CT、乳腺CT、脑CT以及胸腔CT五幅图像。评价指标选取图像分割算法的识别率、运算时间及方差。

仿真参数设置：角度θ∈-10o，10o、公式（7）参数a=0.1，b=0.5，粒子补种参数Np=16。所选取的仿真对比算法为标准LS算法以及本文所提出的LPRLS算法。仿真结果如图4（a）～（e）所示，通过对matlab水平集工具箱函数进行研究，并且结合本文所提算法理论，经过编写接口程序实现LPRLS算法。除上述参数设置外，其他水平集参数设置同工具箱内水平集算法参数设置。水平集算法的分割精度、时间消耗以及方差等评测指标仿真结果如表1所示，表1中的仿真结果为上述两种算法分别独立运行20次所求取的评价指标均值。

图4（a）～（e）显示的分别为标准LS算法和LPRLS算法在血管CT、肿瘤CT、乳腺CT、脑CT以及胸腔CT五组图像上的分割效果。（a）图血管CT图像分割中，LPRLS算法比LS算法在主干血管的识别率上要明显更高，但是存在的问题是，由于灵敏度过高，导致存在一定的识别噪声点，如图中散落的小圆圈；（b）图肿瘤CT图像分割中，LPRLS算法比LS算法更全面的对肿瘤轮廓进行分割，并且肿瘤内部较精细的轮廓也能识别出来；（c）图乳腺CT图像分割中，LPRLS算法对整个乳腺的识别全面性要明显高于LS算法，并且乳腺内部也可较为精确的进行分割识别；（d）图脑CT图像分割中，脑CT图像的纹理结构更加复杂，LS算法的分割效果不是很理想，从分割的精细程度上要明显不如LPRLS算法；（e）图胸腔CT图像分割中，LS算法在左胸中部轮廓识别中存在杂乱现象，整体轮廓也较杂乱，而LPRLS算法对于胸腔轮廓识别要更圆润，识别效果更好。图4从视觉直观性上给出了标准LS算法和LPRLS算法的图像分割效果对比，直观上LPRLS算法要明显优于标准LS算法。表1分别给出标准LS算法和LPRLS算法在血管CT、肿瘤CT、乳腺CT、脑CT以及胸腔CT五组图像上的分割数值对比结果。仿真对比指标

（a）血管CT图像

（b）肿瘤CT图像

（c）乳腺CT图像

（d）脑CT图像

（e）胸腔CT图像

选取识别率、运算时间和识别率方差三种，算法各允许20次求均值。在算法分割识别率指标上，LPRLS算法要明显高于标准LS算法，并且这种优势非常明显，在该指标上LPRLS算法要平均高出20%～30%，如在血管CT、肿瘤CT、乳腺CT三组图像中，标准LS算法的识别率近50%～60%，而LPRLS算法在这三幅图像上的分割识别率均达到85%以上；在脑CT和胸腔CT两组图像分割中，标准LS算法略高达到70%以上，但是LPRLS算法在这两组分割图像的分割识别率均达到90%以上。在计算时间指标对比中，标准LS算法在这五组图像上的运算时间在20s～35s之间，而LPRLS算法在这五组图像上的分割时间则分布在16s～20s之间，可看出在计算时间方面LPRLS算法也要优于标准LS算法。在代表算法稳定性方面的算法方差指标中，两者相差不大，LPRLS算法略优于标准LS算法。

5 结束语endprint

利用Lagrangian粒子增强补种算法对水平集算法进行改进，提高算法应对医学图像CT图像像素不均匀所带来的分割精度降低问题。通过在计算水平集公式前先计算Lagrangian标记粒子来重建内嵌交接界面，从而提高水平集算法的质量守恒特性；通过增加速度矢量和单位法向量来促进奇异点和拓扑变化点速度场的收敛，来处理界面奇异性和复杂几何相关问题。算法仿真结果显示，所提算法比原始算法性能大为提升。从仿真结果可看出，存在的不足是，算法存在误识别率，下一步将主要针对此问题进行进一步分析和研究。

参考文献

[1] OSHER S，SETHIAN J A. Fronts propagating with curvaturedependent speed： algorithms based on HamiltonJacobi formulations[J]. Journal of Computational Physics，1988，79（1）：12-49.

[2] 王斌， 李洁， 高新波. 一种基于边缘与区域信息的先验水平集图像分割方法[J]. 计算机学报， 2012， 35（5）： 1067-1071.

[3] 薛维琴， 周志勇， 张涛. 灰度不均的弱边缘血管影像的水平集分割方法[J]. 软件学报， 2012， 23（9）： 2489-2499.

[4] GIBOU F，MIN C，FEDKIW R. High resolution sharp computational methods for elliptic and parabolic problems in complex geometries[J]. Journal of Science Computing， 2013， 54（2）： 369-413.

[5] BOUMEHED M，ALSHAQAQI B，OUAMRI A. Moving Objects Localization by Local Regions Based Level Set： Application on Urban Traffic[J]. Journal of Mathematical Imaging and Vision， 2013， 46（2）： 258-274.

[6] JEFFREY C L，ZACHARY M. Level sets of the Takagi function： local level sets[J]. Monatshefte für Mathematik， 2012， 166（2）： 201-238.

[7] MARCHANDISE E，REMACLE J F，CHEVAUGEON N. A quadraturefree discontinuous Galerkin method for the level set equation[J]. Journal of Computational Physics， 2006， 212（1）： 338–357.

[8] MU L，WANG J P，WANG Y P. A computational study of the weak Galerkin method for secondorder elliptic equations[J]. Numerical Algorithms， 2013， 63（4）： 753-777.

[9] NAVE J C，ROSALES R R，SEIBOLD B. A gradientaugmented level set method with an optimally local coherent advection scheme[J]. Journal of Computational Physics， 2010， 229（10）： 3802–3827.

[10]LEE C，DOLBOW J，MUCHA P J. A narrowband gradientaugmented level set method for multiphase incompressible flow[J]. Journal of Computational Physics， 2014， 273（5）： 12–37.

[11]杨谊， 喻德旷. 基于改进水平集图像分割方法的乳腺超声病灶提取[J]. 计算机应用与软件， 2014， 31（11）： 217-221.

[12]Li Changyang， Wang Xiuying， Eberl S. A Likelihood and Local Constraint Level Set Model for Liver Tumor Segmentation from CT Volumes[J]. IEEE Transactions on Biomedical Engineering， 2013， 60（10）： 2967-2977.

[13]KACHROO P，RATLIFF L，SASTRY S. Analysis of the GodunovBased Hybrid Model for Ramp Metering and Robust Feedback Control Design[J]. IEEE Transactions on Intelligent Transportation Systems， 2014， 15（5）： 2132-2142.

[14]Zhuang Chijie， Zeng Rong， Zhang Bo. A WENO Scheme for Simulating Streamer Discharge With Photoionizations[J]. IEEE Transactions on Magnetics， 2014， 50（2）： #7007904.

[15]刘世兴， 宋端， 贾林. 辛RungeKutta 方法在求解LagrangeMaxwell方程中的应用研究[J]. 物理学报， 2013， 62（3）： #034501.

[16]OSHERA S，Cheng Litien， KANG M J. Geometric Optics in a PhaseSpaceBased Level Set and Eulerian Framework[J]. Journal of Computational Physics， 2002， 179（2）： 622-648.

[17]王金华， 沈永明， 石峰. 基于拉格朗日粒子追踪的渤海冬季与夏季环流及影响因素[J]. 水利学报， 2011， 42（5）： 544-553.endprint