基于改进水平集方法的Micro-CT鼠脑三维自动化分割

【作 者】陈诗烨，王成＊，郑秀娟2,

1 上海交通大学医学院生物医学工程系，上海，200025

2 上海交通大学医学院附属仁济医院核医学科，上海，200127

0 引言

要开展鼠脑Micro-PET/CT图像自动定量分析，首要步骤是从头部图像中实现对全脑组织的分割和提取。这一过程通常是利用CT图像提取大脑结构，在此基础上分析PET功能图像。尽管在人类的CT或MRI图像中分割并提取大脑组织的方法上已经相当成熟，然而与人类大脑的提取相较，在Micro-PET/CT图像中的分割提取鼠脑的难度更大。其主要原因在于：①需要提取的鼠脑组织非常小，且结构复杂，尤其是前脑部分，即使用人眼也很难分辨；② 虽然Micro-CT拥有极高的分辨率，但图像边缘依然比较模糊，且不连续，存在很多断点，给分割带来很大难度。因此，为了准确分割鼠脑图像，就不能简单套用既有的相对成熟的人脑分割方法。本文提出一种基于水平集的方法，能有效地实现对鼠脑组织的三维自动化分割。

水平集方法最早是由Osher和Sethian[1]在1988年提出，现已在图像处理和计算机视觉等领域得到了广泛的应用[2-3]。2008年，Zhang等人[4]提出了一种基于水平集的混合分割方法，其特点是同时结合了边缘和区域信息，以达到准确分割的效果。其中，边缘信息可以帮助检测目标对象的准确位置，而利用区域信息则可以防止边缘泄漏问题。2009年，Uberti MG[5]提出了一种基于水平集方法并且针对小鼠MRI图像提取大脑组织的半自动方法，其分割效果也较为显著。但是，使用这些近年来所发展的水平集方法来分割鼠脑Micro-CT图像仍然有较多不足：① 水平集初始表面轮廓仍需要手动设置；② 由于图像边缘模糊，存在边缘泄露问题；③ 由于鼠脑组织结构复杂，水平集在前脑等窄细结构处的演化能力不够；④ 停止条件（每次迭代的持续时间和迭代次数）仍然需要手动设置，存在诸多不便。

本文在水平集方法的基础上做了三个方面的改进，提出了一种针对大/小鼠Micro-CT图像的脑组织分割方法，实现了图像的自动化分割处理。

1 方法

1.1 改进水平集方法的框架

本文基于改进水平集方法，提出了一种针对大/小鼠Micro-CT图像的脑组织分割框架，算法流程如图1所示。首先使用模糊C均值聚类方法(Fuzzy-C-Means,FCM)，实现初始水平集表面轮廓的自动化设置；接着使用梯度矢量流(Gradient Vector Flow, GVF)增加外向演化能力，使得演化表面轮廓能更准确地进入窄细处等复杂结构；最后提出一种判断演化表面轮廓附近平均带宽能量是否达到最大化(Average Bandwidth Energy Maximization, ABEM)的方法，作为停止条件，防止泄漏且能同时实现自动化停止。

图1 算法总流程Fig.1 Algorithms of the proposed method

1.1.1 模糊C均值聚类方法的预分割

FCM方法是由Bezdek等[6]提出，该算法具有良好的局部收敛性和分割效果，被成功地用于图像分割中[7-8]。因为脑组织和周围组织存在差别，所以使用FCM方法，并结合阈值和形态学方法，可以获取大致的脑部形态结构，作为自动预分割。先通过设置一个参考阈值（该阈值由已经发展成熟的人脑不同组织的CT值等先验知识可以确定[9]），分割出部分头骨组织；再通过膨胀、腐蚀和填充孔洞等形态学方法，提取包含整个脑组织的图像区域，去除部分周围组织的干扰；接着利用FCM方法对阈值和形态学方法分割得到的一次结果进行再分割，获得更为精确的脑组织区域。但是，经过FCM分割的图像边缘不够连续，周边有部分噪声，图像上有很多孔洞。为了有效解决上述问题，再次利用形态学方法对图像进行处理，其结果作为下一步水平集演化的初始图像。

1.1.2 梯度矢量流的水平集演化

经典水平集演化是以边界上点的演化速度作为动力的。在图像处理上，常以灰度梯度作为计算速度的要素，据此水平集方法的演化方程可以简化为：

其中，I是需要分割的图像，μ表示灰度级下限的预定义参数，φ代表活动轮廓曲线C={x│φ(x)=0}的初始水平集，H(φ)是Heaviside 方程，其定义为是图像梯度图，Ω表示图像域，α和β是用来调节方程两个因子比重所预定义的参数。

1.1.3 平均带宽能量最大化 (ABEM)

为了能自动获得最佳演化结果，本文提出了一种判断演化表面轮廓周围一定带宽区域的平均能量是否达到最大化，作为停止水平集继续演化的条件。平均带宽能量(Average Bandwidth Energy，ABE) 被定义为

其中，I(i)代表第i个像素的灰度级，n代表带宽中的像素数。带的宽度可根据图像来设置，以期获得最佳效果。通常是带宽越宽，带宽中所包含的图像信息也越大，随之计算量也会更大，但是会更加准确。

以一幅二维采样图像为例（见图2所示），需要分割的边缘为像素高亮的脑骨骼组织，标示为白色的演化曲线不断向边缘靠近，它周围的黑色区域为该演化曲线的窄带。因为边缘的像素往往具有更高的灰度等级，且与周围组织对比明显，所以当演化曲线不断向图像边缘靠近时，窄带中所包围的平均能量也越大。当其达到最大化时，窄带中所包含的平均能量达到最大(图3为该层图像的平均带宽能量曲线，最大值如五角星所指处)，即达到最接近图像边缘的程度。

图2 平均带宽能量最大化的原理Fig.2 Principle of ABEM

1.2 算法分割性能评价

本文采用基于区域法的评价方法，即计算两个对应区域的重叠比例，这也是国际上公认的评价方法之一（Uberti MG[5]）。重叠指数定义为

图3 平均带宽能量曲线Fig.3 Energy curve of ABEM

其中 V1和 V2为两幅需要对照的分割结果图像，Vo是两幅图像的重叠区域。 OM 的范围在0到1之间，0代表两幅图像之间没有重叠部分，而1则代表两幅图像完全重叠。

为了评价本文所提出方法的分割效果，将手动分割的结果作为评价本文所提出方法是否分割准确的标准，公式(5)中的V2由手动分割获得，采用了Analyze®软件[12]中勾画感兴趣区域的功能进行一层层的分割并整合而得。V1则分别采用本文所提出的方法和Zhang等人在其论文中所提到的混合分割方法(Hybrid Level-Set，HLS)而得。2D/3D image segmentation toolbox®是Zhang等人在其论文中所提到的基于Matlab所开发的工具箱，用该方法获得的图像作为本文所提出的图像处理方法的对照。为了使用的方便性，其中初始水平集表面轮廓的获取也采用了利用FCM方法，其他步骤和方法基本参照原方法，没有变动。方法中所设定的参数也通过反复试验调节所得。

2 实验结果

2.1 初始水平集表面轮廓的获取

从图4可以看出，通过使用基于FCM方法自动获取的初始水平集表面轮廓（图4中以深色曲线表示），已经非常接近脑组织边界。但是，前脑部分以及脑部和颈部连接处的分割效果还是不够。以大鼠一为例，经过FCM预分割的初始水平集表面轮廓与手工分割的结果作对照，其重叠指数为0.628。此与重叠指数为0.881的最终分割结果相比，还有一定距离，所以需要利用水平集方法作进一步的演化处理。

2.2 水平集演化的最终结果

基于FCM预分割得到的初始水平集表面轮廓后，使用本文提出的GVF水平集演化，并通过判断平均带宽能量最大化作为停止的条件，分别获得如图5和图6的三维和二维图像。

图4 初始水平集分割结果的二维采样例图Fig.4 Sampled axial slices of the obtained initial level-set surface

图5 最终分割结果的三维例图Fig.5 Segmentation result in 3D

图6 最终分割结果的二维例图Fig.6 Sampled axial slices of segmentation result

2.3 分割结果评价

重叠指数的比较结果如图7所示，OM_3D为对三维结果计算得到的重叠指数，OM_2D为分别对各层横截面计算的平均重叠指数。从三维对比的结果看，使用本文提出的分割方法与手动分割的法计算所得的大鼠和小鼠的重叠指数平均值分别为0.883和0.874，标准偏差分别为0.003和0.011。与之相对应的，使用Zhang等所提出的方法分割的结果，重叠指数平均值分别为0.663和0.819，标准偏差分别为0.008和0.009。比较而言，大鼠的分割改善效果小鼠较好。相比Zhang的方法，大鼠的分割准确率提高了33%，小鼠提高了6.7%。可以看到，本文所提出的方法在准确性上有显著改善。

图7 两种分割方法重叠指数比较Fig.7 Comparison of the proposed method versus HLS by using OM

2.4 算法的处理速度

本文提出的算法是基于Matlab平台开发的，所用的计算机配置为Intel Core 2 Quad CPU: 2.83GHz，RAM: 2 GB，Matlab版本为R2010a。图像大小不同，其处理速度也不同，比如处理大小约为160×200×250像素的大鼠图像，水平集演化过程需要8 min左右的时间。根据处理时间与图像大小的比值得到大/小鼠的处理效率，详细情况如表1所示。

表1 本文算法的处理效率Tab.1 The processing ef fi ciency of the algorithm

3 讨论

随着成像技术的发展，利用Micro-PET/CT技术来对动物模型做研究正在被越来越多研究者所采用。然而，就我们所知，还没有专门针对大/小鼠Micro-CT图像的脑组织自动化提取方法。因此，本文提出了基于改进水平集方法的算法框架，可以实现对鼠脑组织的Micro-CT图像自动化分割，而且在准确性上有较大改善，可以作为后续进一步量化分析的基础。

3.1 关于框架流程的考虑

本文所提出的框架是在水平集演化之前加入了基于FCM方法的预分割，虽然这个部分需要耗费额外的处理时间，但是对于整个结果的准确性和便捷性都是很有用的。使用FCM方法并结合阈值和形态学的方法来获取初始水平集表面轮廓，不仅避免了手动设置，而且获取的初始水平集表面轮廓已经非常接近目标边缘，因此为后续演化过程节约了大量时间，并且实现了全自动化获取，不需要手动勾画。

3.2 使用梯度矢量流的改善效果

引入梯度矢量流，增强了水平集演化的外部张力，使在图像的窄细等复杂结构处的演化能力大大增强。通过测试比较发现，在同样使用FCM预分割得到的结果作为初始水平集表面，并且其他参数均设置相同的条件下，以水平集表面轮廓成功演化到前脑部分以及脑部和颈部等窄细结构处为参照标准，使用GVF作为梯度图通常比使用一般梯度方程要快大约10个迭代过程。从分割的整体效果看，使用一般梯度方程在演化表面轮廓成功到达前脑等区域的同时，其横截面的其他大部分层面均出现了较为严重的泄漏，但是观察使用GVF的结果图，泄漏情况大为改善，且各个层面的分割进度较为统一。所以实验结果证明了GVF扩张梯度场张力的作用，可以加快水平集表面轮廓更快地演化到远端，同时加强了在窄细结构处的演化能力。

3.3 平均带宽能量最大化 的改善效果

通过使用平均带宽能量最大化作为停止条件，省去了人工设置的麻烦，较好地解决了部分层面泄漏的问题。这个过程是在二维层面完成的，因为在水平集演化过程中，不同层面的演化进度是不一样的，虽然通过使用GVF扩张了梯度场张力，使得水平集表面轮廓可以更快地演化到远端，例如前脑等窄细结构处，但是还是有部分中间层面出现了泄漏的情况。为了使水平集表面轮廓在成功演化到远端的同时，控制中间层面的演化进度，就需要计算各个横截面的平均带宽能量，并分别判断达到最大化的时间和结果，最终得到优化的结果。虽然在二维层面处理增加了处理的时间，但是与水平集演化所需要的时间相比，这部分时间可以忽略。

3.4 参数的优化选择

在水平集演化方程中，主要涉及参数有μ、α和β。如上文所介绍的，μ为全局阈值，用来表示目标对象边界的低灰度级。演化方程(1)的第一部利用图像区域信息，推动曲面轮廓向大于 的方向演化。μ的选择不是一件容易的事，如果选择太小，则可能会导致演化曲面轮廓越过边界从而过度分割；如果选择太大，则可能会导致演化曲面轮廓始终到达不了目标边缘。此外，由于图像不均匀，很难确定一个特定的值赋值给。本文采用LBF(local binary fi tting energy)[13-14]，—— 一种局部二值拟合能量泛函，根据局部图像信息设定该阈值。与全局阈值法相比，使用LBF能够提取更准确的图像信息。

α和β是用来调节方程两部分比重所预定义的参数。α为演进权重指数，α越大，则水平集演化所利用的区域信息比重越大，且受μ值的影响也大。演化方程（1）的第二部分利用梯度信息，使演化曲面轮廓向梯度值大的图像边缘靠近，β为该部分的权重指数。关于α和β的选取，以大鼠为例，通过穷尽法选取正确的数量级区间，观测3例图像分割结果的平均重叠指数效果变化，比较好的结果出现在β=0.001，α在0.001到0.1之间，针对这个结果，再进行一定范围的测试，观察重叠的效果变化，得到最优近似值为α=0.009，β=0.001 。利用同样的方法对小鼠的和 也做了选取，得到最优近似值为α=0.009，β=0.001。

4 结论

本文对大/小鼠脑组织提取框架是基于水平集方法的改进，并被验证具有较好的准确性和便捷性。通过与已有经典方法的对比，本文所提出的方法被证明具有更好的效果。利用该方法，研究人员可以很容易地对大/小鼠Micro-CT图像进行自动化处理，为进一步与图谱（例如Digimouse[15]）配准，获取丰富的脑区信息或不同个体脑组织之间的比对，打下了基础。

致谢

感谢上海交通大学医学院附属瑞金医院小动物PET中心和上海交通大学MED- X研究院在本研究中所提供的帮助。

[1] Osher S, Sethian JA. Fronts propagating with curvature dependent speed: algorithms based on the Hamilton-Jacobi formulation [J].Journal of Computational Physics, 1988, 79(1): 12-49.

[2] Caselles V, Catte F, Coll T, Dibos F. A geometric model for active contours in image processing [J]. Numeric Math, 1993, 66(1): 1-31.

[3] Malladi R, Sethian JA, Vemuri BC. Shape modeling with front propagation: a level set approach [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1995, 17(2): 158-175.

[4] Zhang Y, Matuszewski BJ, Shark L, et al. Medical image segmentation using new hybrid level-set method [A]. In:Proceedings of the 2008 Fifth International Conference BioMedical Visualization: Information Visualization in Medical and Biomedical Informatics[C], 2008, 71-76.

[5] Uberti MG, Boska MD, Liu Y. A semi-automatic image segmentation method for extraction of brain volume from in vivo mouse head magnetic resonance imaging using Constraint Level Sets [J]. J Neurosci Methods, 2009, 179(2): 338-344.

[6] Bezdek JC. Pattern recognition with fuzzy objective function algorithms [M]. New York: Plenum Press, 1981, 95-107.

[7] Pham D, Prince J. An adaptative fuzzy C-means algorithm for image segmentation in the presence of intensity inhomogeneities [J].Pattern Recognition Letter, 1999, 20(1): 57-68.

[8] W.J. Chen, M.L. Giger, U. Bick. A fuzzy c-means (FCM)-based approach for computerized segmentation of breast lesions in dynamic contrast enhanced MRI images [J]. Acad. Radiol, 2006,13(1): 63-72.

[9] 李传富, 周康源, 黄丹, 等. 基于先验知识的颅脑CT图像自动化分割[J].中国医疗器械杂志, 2004, 28(3):168-171.

[10] Chenyang Xu and Jerry L. Prince. Gradient Vector Flow: A New External Force for Snakes [J]. IEEE Transactions on Image Processing, 1998, 7(3): 359-369.

[11] Chenyang Xu, Prince, JL. Snakes, shapes, and gradient vector fl ow[J]. IEEE Transactions on Image Processing, 1998, 7(3):359-369.

[12] Robb RA. The biomedical imaging resource at mayo clinic [J].IEEE Transactions on Medical Imaging, 2001, 20(9):854-867.

[13] C M Li, C Kao, J Gore, et al. Implicit active contours driven by local binary fi tting energy [A]. In: IEEE Conference on Computer Vision and Pattern Recognition [C]. Minneapolis, 2007, 1-7.

[14] Qingqi Hong, Qingde Li, Jie Tian. Local hybrid level-set method for mra image segmentation[A]. In: 2010 10th IEEE International Conference on Computer and Information Technology (CIT 2010).Bradford [C], 2010, 1397- 1402.

[15] Dogdas B, Stout D, Chatziioannou AF, et al. Digimouse: a 3D whole body mouse atlas from CT and cryosection data [J]. Physics in Medicine and Biology, 2007, 52(3), 577-587.