一种普适的基于多尺度滤波和统计学混合模型的血管分割方法

陆 培 王 磊 李志成 周寿军

(中国科学院深圳先进技术研究院生物医学与健康工程研究所，广东 深圳 518055)

一种普适的基于多尺度滤波和统计学混合模型的血管分割方法

陆 培 王 磊 李志成 周寿军*

(中国科学院深圳先进技术研究院生物医学与健康工程研究所，广东 深圳 518055)

血管的精确提取和定位，是实现心脑血管介入手术的关键。多尺度滤波算法可以增强血管目标，同时抑制背景噪声，但并没有把血管从图像背景中区分出来。基于统计学的分割算法，通过对直方图进行拟合实现血管的分类，但需要调整混合模型去拟合特定的图像直方图。为了克服上述问题，提出一种具有固定模型的普适的血管分割方法。首先，利用多尺度滤波算法进行图像预处理。其次，针对滤波增强后数据的直方图曲线，用由3个概率分布函数(1个高斯和2个指数)组成的混合模型进行拟合。期望最大化算法用于混合模型参数的估计。最后，通过最大后验概率分类算法将血管分离出来。为了验证上述方法的有效性，分别在仿真(phantom)数据、磁共振血管造影(MAR)数据和计算机断层血管造影(CTA)数据上进行实验测试。结果表明，所提出的方法在多套仿真数据上的分割误差低于0.3%，同时对于不同模态的血管图像具有很好的分割效果及较强的鲁棒性。

血管分割；多尺度滤波；混合模型；多模态图像

引言

目前，心脑血管疾病已经严重威胁着人类的健康。随着CTA(CT血管造影)、MRA(磁共振血管造影)等成像技术的发展，医学图像的后处理技术变得越来越重要，精确的血管分割是血管类疾病诊断与治疗的关键[1]。

在血管分割中，多尺度滤波、形变模型、统计模型及混合法得到了广泛的研究和应用[2]。多尺度滤波通过抑制背景噪声增强血管，但是血管并没有与背景区分开[3]。形变模型有效地整合了自底向上的信息和自上而下的先验知识，但分割质量主要依赖于模型参数[4]。统计学模型主要是根据像素点的强度信息进行分类，但其模型选择受限于图像的统计学特性[5]。混合方法是将多种分割算法相结合，以降低血管误分割的概率[6]。

到目前为止，基于统计学的分割算法已经引起了越来越多的关注，而模型选择则是基于统计学的分割算法的关键。针对脑部MRA图像，Wilson和Noble于1999年提出了一种基于统计学模型的血管分割方法[7]。根据灰度值特性，他们用一个均匀分布拟合血管类，用两个高斯分布拟合背景类。2006年，Hassouna等对模型进行了改进，他的方法中用一个瑞利函数和两个高斯函数拟合背景噪声，同时用一个高斯函数拟合血管，并取得了更好的血管分割效果[8]。2013年，Zhou等将统计学模型和多模式邻域系统相结合，提取出了更精细的血管目标[9]。然而，不同模态的血管图像或者不同器官的血管成像具有不同的统计学特性，很难用同样的混合模型去精确拟合。

针对上述问题，笔者将多尺度滤波和统计学方法相结合，提出了一种普适的血管分割方法。首先，将不同模态的图像数据通过多尺度滤波进行血管增强，增强后的数据具有一个新的相近的统计学特性；其次，针对滤波后的数据，用一个高斯函数拟合血管类，同时根据直方图特性，选用两个指数函数拟合背景噪声类，模型参数选用最大期望(EM)算法进行估算；最后，利用最大后验概率分类，将血管和背景标记出来。

1 方法

图1为笔者所提出的血管分割算法的流程。首先，用多尺度滤波进行图像预处理；其次，通过直方图分析，选用一个混合模型对滤波增强后的数据进行拟合；最后，通过最大后验概率将血管类标记出来。

图1 本算法流程Fig.1 The framework of the proposed method.

1.1 血管增强

由直方图曲线可以看出，不同模态的图像数据具有不同的统计学特性，因此很难找到一个固定的混合模型去拟合所有的直方图曲线。本研究利用Frangi的多尺度变换方法，通过Hessian矩阵的特征值计算局部血管出现的概率，从而得到一个具有新的统计学特性的血管增强图像[10]。

图像I在点x0处的泰勒展开可表示为

(1)

这里，图像的微分定义为与高斯函数的卷积，有

(2)

式中，参数γ为Lindeberg引进来的一个规则化的倒数族[11]。

D-维高斯函数定义如下：

(3)

(4)

由特征值可以获得3个不受尺度因素影响的正交方向，根据管状结构特点，由下式可以鉴别出血管结构：

(5)

在多尺度滤波中，引入了3个测度函数如下：

(6)

式中，RA用来区分片状和线状结构，RB用来区分点状结构和线状结构，S用于区分背景像素。

因此，在多尺度滤波中的血管增强函数如下：

(7)

这里，α、β和c作为阈值，用于控制上述血管增强算法对RA、RB和S的敏感性。

1.2 混合模型及参数估计

针对滤波增强后的数据，通过直方图分析，用一个固定的混合模型来进行拟合。经过多次测试，高灰度值的血管部分用一个高斯分布来拟合，低灰度值部分用两个指数分布来拟合。因此，本研究所提出的混合模型可以描述为

(8)

式中，f(x)为增强后图像整体的混合概率密度函数，wEl(l=1,2)和wG为三类分布所占的比例，fEl(x)(l=1,2)和fG(x)分别表示指数分布和高斯分布，有

(9)

在进行模型参数估计之前，先要将滤波增强的数据范围变换为与原始图像一致，随后，通过以下几步完成参数的估计：

第一步，利用k均值聚类进行初步分类。如果原始图像的直方图曲线中峰值点超过3个，则用前3个峰值点进行k均值初始化；未超过3个，则用如下规则进行k均值算法的初始化，即

(10)

式中，Imax是原始图像中像素最大值，μk1，init、μk2，init和μk3，init表示k均值算法的初始聚类中心，μkl、σkl、wkl(l=1,2,3)是由k均值聚类算法估计得到的3个分布的模型参数。

第二步，用EM算法进行混合模型参数优化。本研究用k均值的初步分类结果进行EM迭代算法的初始化，有

(11)

式中，λE1，init、wE1，init、λE2，init、wE2，init、μG，init和wG，init表示EM算法迭代的初始值。

用EM算法进行各个模型参数的迭代更新如下：

(12)

(13)

式中，i表示迭代次数，xj是第j个体素的强度，N是体素的总数，El表示第l(l=1,2)个指数模型，G表示高斯模型。

后验概率fi(El|xj)和fi(G|xj)可以由贝叶斯公式计算如下：

(14)

1.3 最大后验概率

根据最大后验概率估计法则，由下式可以把血管从背景中标记出来，有

(15)

2 结果

为了验证所提出算法的有效性，本研究分别进行了一系列Phantom、MRA、CTA的血管分割实验及与其他分割算法的对比实验。

2.1 仿真数据的血管分割

为了定量地验证本研究所提出算法的有效性，构建了多套血管仿真数据进行分割测试，并通过对比加有多种高斯噪音的Phantoms的分割结果与仿真数据中提前设定的血管真实值，以得到血管分割的误差。图2给出了两套仿真数据的3D效果和切片显示。其中，PhantomⅠ只包含一条血管结构，PhantomⅡ包含多条血管结构，且两套仿真数据中的血管直径像真实血管一样由粗到细逐渐变化。两套仿真数据的大小分别为128像素×128像素×128像素和240像素×240像素×240像素。

图2 两套仿真数据。(a)真实值；(b)真实值切片；(c)加了噪音后的切片Fig.2 Two phantoms.(a)Ground truth；(b)Slices of ground truth；(c)Slices of the phantoms after the noises are added

2.2 不同血管造影图像的分割

利用本研究所提出的算法，分别在网上下载的10套脑部血管造影数据(http://www.insight-journal.org/midas /community/view/21)、深圳市第二人民医院提供的3套胎盘数据和5套腿部数据上进行了分割实验。实验中相同部位的数据来自于相同的成像设备，因而具有较为一致的统计学特性，因此该实验部分仅给出了同一部位的1套数据的分割结果。图4为进行血管分割的三维图像数据的切片。

图3 两套仿真数据的分割结果。(a)某一切片；(b)三维效果Fig.3 Segmentation results of the two phantoms.(a)A slice；(b)3D view

图4 血管图像的切片。(a)脑部；(b)胎盘；(c)腿部Fig.4 Slices of different vascular images.(a)brain；(b)placenta；(c)legs

原始的三维图像经过多尺度滤波后，背景噪声被抑制，同时血管目标得到加强。图5分别为图4中3个图像数据的原始直方图曲线，以及对应的滤波增强后的直方图曲线。可以看出，具有不同统计学特性的原始图像数据经过滤波增强后，可以得到一个新的具有近似统计学特性的数据。

图5 直方图曲线 (每行左为原始直方图曲线，右为增强后直方图曲线)。(a)脑部；(b)胎盘；(c)腿部Fig.5 Histogram curves of the original and the corresponding enhanced data(In each the left is the histogram curves of the original image, the right is the histogram curves of the corresponding enhanced data).(a)Brain；(b)Placenta；(c)Legs

针对滤波增强后的数据，经过多次测试，选用由两个指数函数和一个高斯函数组合的混合模型进行拟合。在图6中，黑色实线表示增强后数据的直方图曲线，加点黑色实线表示混合模型拟合后的直方图曲线。可以看出，本研究所提出的混合模型对于不同的血管造影图像均有很好的拟合效果。

图6 直方图拟合。(a)脑部；(b)胎盘；(c)腿部Fig.6 Three estimated distributions together with the total fitting curve.(a)Brain；(b)Placenta；(c)Legs

图7 血管分割结果(每行至上而下为图像切片、Z轴投影和三维视图)。(a)脑部；(b)胎盘；(c)腿部Fig.7 The results of the vessel segmentation(The three lines from the first one to the third one are the image slices, the Z-axis projections, and the 3D views).(a)Brain；(b)Placenta；(c)Legs

在图7中，第一行为与图4对应的切片的血管分割结果，第二行为血管分割后的三维数据的Z轴投影图像，第三行为分割结果的三维视觉效果。可以看出，本研究所提出的方法对于不同的血管造影图像具有较好的分割效果。

2.3 同一图像不同分割算法的对比

最后，为了进一步测试本研究提出的分割算法，针对广州军区总医院提供的同一套脑部MRA数据和胸部CTA数据，分别用文献[9]的分割算法及本研究所提出的方法进行血管分割。脑部MRA数据大小为561像素×561像素×361像素，数据间距为0.50像素×0.50像素×0.50像素；胸部CTA的数据大小为512像素×512像素×338像素，数据间距为0.36像素×0.36像素×0.4像素。图8、9给出了两套数据在两种方法下的直方图拟合结果，以及血管分割效果的Z轴投影。

图8 脑部MRA数据模型拟合及分割结果。(a)文献[9]方法；(b)本方法Fig.8 Histogram fitting from brain MRA and the corresponding segmentation results.(a)Method in reference[9]; (b)The proposed method

图9 胸部CTA数据模型拟合及分割结果。(a)文献[9]方法；(b)本方法Fig.9 Histogram fitting from thorax CTA and the corresponding segmentation results.(a)Method in reference[9]; (b)The proposed method

由两套数据的拟合效果可以看出，用文献[9]所提出的方法进行的模型拟合受限于原始图像的统计学特性，而本研究中所提出的混合模型能够较好地拟合滤波增强后的数据。新提出的分割算法在两套数据上均优于传统的分割算法，验证了该算法的有效性，以及在不同模态的血管造影图像中具有较强的鲁棒性。

3 讨论

基于统计学的血管分割方法，根据原始的血管造影图像的直方图特性，构建相应的混合模型，对其直方图曲线进行拟合，最终完成血管的分类。然而，在利用该方法进行血管分割时，需要根据数据的直方图特性选择相应的统计学混合模型。为了克服这个限制，笔者针对血管造影图像，首先用多尺度滤波算法进行预处理，以实现血管的增强和背景的抑制，再利用固定的统计学混合模型，对滤波增强后的数据进行血管的提取。笔者所提出的普适的血管分割方法特点在于：一是具有不同统计学特性的血管造影图像通过多尺度滤波进行图像预处理，能够得到具有相对一致的统计学特性的血管增强数据；二是针对增强后的数据，可以用由一个高斯分布和两个指数分布组成的固定混合模型来进行直方图的拟合；三是本方法适用于不同部位及不同模态的血管造影图像。

实验1通过将所提出的分割方法在加有不同噪声的仿真数据上进行测试，给出了本研究所提出方法的定量评估结果。实验表明，本方法在仿真数据上的血管分割结果与其预设的真实值相比误差低于0.3%，说明具有较高的分割精度。在实验2中，对不同部位的血管造影数据进行了分割测试。由图5可以看出，不同部位的血管造影图像具有不同的统计学特性，而经过多尺度滤波算法增强后数据的直方图曲线具有较为一致的统计学特性。由图6直方图拟合结果及图7血管分割结果可以看出，本方法对于不同部位的血管造影图像均能给出较好的拟合和分割结果。实验3给出了本方法及经典的血管分割方法在两套不同模态的血管造影图像中的血管分割结果。文献[9]中所提的方法用由一个瑞利分布和3个高斯分布组成的混合模型对原始血管图像直方图进行拟合。由于实验3中的两套血管造影数据具有不同于文献[9]中图像数据的统计学特性，因此文献[9]中所提出的混合模型在这两套数据上表现出较差的分割性能。而本研究中采用一个高斯分布和两个指数分布组成的混合模型，对增强后的血管数据进行拟合，因此具有更广的应用范围。实验3的分割结果表明，本方法对于不同模态的血管数据具有较好的分割结果和更强的鲁棒性。

本研究将多尺度滤波与统计学混合模型相结合，很好地解决了统计学模型在进行血管分割过程中所存在的局限性问题。在后续的研究中，将会在更多的血管造影图像中测试本算法，并进行不断的改进和优化，在保证血管分割质量的前提下，扩大其应用范围。

4 结论

本研究针对不同部位、不同模态的血管造影图像，提出了一种普适的血管分割方法，可以广泛应用于血管类疾病的诊断和治疗，在手术路径规划、血管介入手术导航方面具有重要的应用价值。笔者的日后工作将集中于通过模型参数进行优化估计，从而进一步提高血管分割的精度。

[1] 金可, 刘苏. 血管影像分割技术综述[J]. 中国制造业信息化, 2012, 41(11): 61-64.

[2] Mesejo P, Valsecchi A, Marrakchi-Kacem L, et al. Biomedical image segmentation using geometric deformable models and metaheuristics[J]. Computerized Medical Imaging & Graphics the Official Journal of the Computerized Medical Imaging Society, 2015, 43: 167-178.

[3] Ye DH, Kwon DJ, Yun ID, et al. Fast multiscale vessel enhancement filtering[J]. Processings of the SPIE, Medical Imaging, 2008, 6914: 691423-691428.

[4] Gao Xin, Uchiyama Y, Zhou Xiangrong, et al. A fast and fully automatic method for cerebrovascular segmentation on time-of-flight (TOF) MRA image[J]. Digit Imaging, 2011, 24 (4): 609-625.

[5] 曹容菲, 张美霞, 王醒策, 等. 基于高斯-马尔科夫随机场模型的脑血管分割算法研究[J]. 电子与信息学报, 2014, 36(9): 2053-2060.

[6] 王醒策, 文蕾, 武仲科, 等. 面向时飞磁共振血管造影术的脑血管统计分割混合模型[J]. 光学精密工程, 2014, 22(2): 497-507.

[7] Wilson DL, Noble JA. An adaptive segmentation algorithm for time-offlight MRA data[J]. IEEE Trans Med Imaging, 1999, 18 (10): 938-945.

[8] Hassouna MS, Farag AA, Hushek S, et al. Cerebrovascular segmentation from TOF using stochastic models[J]. Med Image Anal, 2006, 10: 2-18.

[9] Zhou Shoujun, Chen Wufan, Jia Fucang, et al. Segmentation of brain magnetic resonance angiography images based on map-mrf with multi-pattern neighborhood system and approximation of regularization coefficient[J]. Medical Image Analysis, 2013, 17(8): 1220-1235.

[10] Frangi AF, Niessen WJ, Vincken KL, et al. Multiscale vessel enhancement filtering[C]//Proceedings of the International Conference on Medical Image Computing Computer Assisted Intervention. Lect: Notes Comp Sci, 1998, 1496: 130-137.

[11] Lindeberg T. Edge detection and ridge detection with automatic scale selection[J]. International Journal of Computer Vision, 1998, 30(2): 117-156.

A Universal Vessel Segmentation Method Based on Multi-Scale Filtering and Statistical Mixture Model

Lu Pei Wang Lei Li Zhicheng Zhou Shoujun*

(InstituteofBiomedicalandHealthEngineering,ShenzhenInstitutesofAdvancedTechnology,ChineseAcademyofSciences,Shenzhen518055,Guangdong,China)

Accurate extraction and localization of blood vessels are the keys to the intervention operation of cardiac and cerebral vessels. Multi-scale filtering strengthens the vessels while weaken the background voxels, but the vessels are still not marked out. Statistical based segmentation method classifies the vessels through model fitting for the histogram curve, but it needs to adjust its model to fit a certain image histogram. To overcome these problems, a universal vessel segmentation method with a fixed model has been proposed in this paper. Firstly, the original image was preprocessed with multi-scale vessel enhancement algorithm. Secondly, a mixture model formed by three probabilistic distributions (one normal distribution and two exponentials) was built to fit the enhanced data. Expectation maximization algorithm has been used for parameters estimation. Finally, the vessels were segmented by maximum a posteriori classification. To test the effectiveness of the proposed method, experiments have been done on a series of phantoms, magnetic resonance angiography (MRA) data and computed tomography angiography (CTA) data. As a result, the segmentation errors of the phantoms are less than 0.3%. Meanwhile, the proposed method performed well on multi-modality images with strong robustness.

vessel segmentation; multi-scale filtering; mixture model; multi-modality images

10.3969/j.issn.0258-8021. 2016. 05.002

2016-01-12， 录用日期:2016-06-28

国家自然科学基金国家高技术研究发展计划(863计划)(2015AA043203)；广东省创新研究团队项目(2011S013)

R318

A

0258-8021(2016) 05-0519-07

*通信作者(Corresponding author)， E-mail: shoujz@163.com