Print

发布时间: 2016-07-25
摘要点击次数: 288
全文下载次数: 39
DOI: 10.11834/jig.20160710
2016 | Volumn 21 | Number 7




    医学图像处理    




  <<上一篇 




  下一篇>> 





形状约束下活动轮廓模型冠脉血管图像多尺度分割
expand article info 郭笑妍, 梅雪, 李振华, 曹佳松, 周宇
南京工业大学电气工程与控制科学学院, 南京 211816

摘要

目的 由于计算机断层血管造影(CTA)图像的复杂性,临床诊断冠脉疾病往往需要经验丰富的医师对冠状动脉进行手动分割,快速、准确自动分割出冠状动脉对提高冠脉疾病诊断效率具有重要意义。针对双源CT图像特点以及传统单一基于区域或边界的活动轮廓模型的不足,研究了心脏冠脉3维分割算法,提出一种基于血管形状约束的活动轮廓模型分割方法。 方法 首先,利用改进的FCM(fuzzy C-means)对心脏CT图像感兴趣区域初分割,其结果用于初始化C-V模型水平集演化曲线及控制参数,提取感兴趣区域轮廓。接着,由3维心脏图像数据获取多尺度梯度矢量信息构造边界型能量泛函,然后利用基于Hessian矩阵的多尺度血管函数对心脏感兴趣区域3维体数据增强滤波,获取血管先验形状信息用于约束能量泛函。最后融合边界、区域能量泛函并利用变分原理及水平集方法得到适合冠脉血管分割的水平集演化方程。 结果 由于血管图像的灰度不均匀,血管末端区域更为细小,所以上述算法的实施是面向被划分多个子区域的血管,在缩小的范围内进行轮廓的演化。相比于传统的血管分割方法,该方法充分融合血管图像的先验信息及梯度场信息,能够从灰度及造影剂分布不均匀的冠脉血管图像中准确分割出冠状动脉,对于细小的血管结构亦能获得较好的分割效果。实验结果表明,该方法只需在给定初始轮廓前提下,有效提取3维冠脉血管。 结论 对多组心脏CT图像进行分割,本文基于血管先验形状约束的活动轮廓模型可以准确分割出冠脉结构完整轮廓,并且人工交互简单。该方法在双源CT冠脉图像自动分割方面具有较好的正确率与优越性。

关键词

冠脉分割, 心脏CT图像, 活动轮廓模型, 形状约束, 水平集方法

Multi-scale segmentation of coronary artery image based on active contour model with shape constraints
expand article info Guo Xiaoyan, Mei Xue, Li Zhenhua, Cao Jiasong, Zhou Yu
College of Electrical Engineering & Control Science, Nanjing Tech University, Nanjing 211816, China

Abstract

Objective Owing to the complexity and randomness of CTA images in the clinic, clinicians are required to manually segment the coronary artery to diagnose a patient suffering from coronary artery diseases. How to segment the coronary artery quickly, accurately, and automatically is important for diagnostic efficiency. According to the features of dual-source CT image, a three-dimensional segmentation algorithm for cardiac coronary artery is investigated. To solve the difficulty of single region or boundary-based active contour model, the segmentation method of active contour model-based vessel prior shape is proposed. Methods First, the improved fuzzy C-means algorithm is used for the initial coarse segmentation of the region of interest (ROI) in cardiac CT images to extract the ROI. The coarse segmentation result is applied to initialize the C-V model level set automatically and estimate the controlling parameters for level set evolution. Next, the 3D heart volume data are used to obtain multiscale gradient vector information and determine the structure of the boundary-based energy function. Then, the vascular structures within the 3D ROI volume are enhanced using the Hessian matrix-based multiscale vascular function. After vascular filtering, the vessel's prior shape information can be obtained to construct the region-based energy function. Finally, the boundary- and region-based energy functions are fused to construct the hybrid energy function, and the variation principle and level set method are used to derive the level set evolution equation for coronary vessel segmentation. Results The gray scales of vessel images are uneven, and the area of the terminal vessel is smaller than the others; thus, the aforementioned algorithm will be divided into several small-scale subareas of the vessels for the evolution of the contour. The vessel prior information and gradient information are fused to segment accurately the coronary artery from the image with uneven distribution of the gray and contrast agents, better than the conventional vessel segmentation method. Furthermore, good segmentation effect can further be obtained even for a small vessel structure. Experimental results show that the proposed method is only dependent on the initial contour to extract 3D vascular images effectively. Conclusion Several groups of heart volumes are tested to evaluate its performance in this study. The results show that active contour model-based vessel prior shape can accurately segment the complete coronary structures with a simple human interaction. The method has good accuracy and superiority in terms of the automatic segmentation of a dual-source CT coronary artery image.

Key words

coronary artery segmentation, cardiac CT image, active contour model, shape constraint, level set method

0 引 言

冠状动脉的分割对冠脉钙化、动脉瘤及冠脉狭窄等冠脉疾病的诊断、病灶检测、疾病的评估以及外科手术方案的制定提供了良好的数据基础。临床诊断冠脉疾病往往需要经验丰富的医师对冠状动脉进行手动分割,因此快速、准确自动分割出冠状动脉对提高冠脉疾病诊断效率具有重要意义[1-2]

基于活动轮廓模型(ACM)[3]的图像分割方法由于其能够自动处理图像中拓扑结构的变化,具备高维曲线描述能力、能量函数模型易于建立、模型易于引入目标先验信息(目标形状、亮度等)等优势,在血管图像分割领域的应用日趋广泛[4-6]

依据活动轮廓线的表达方式和实现方法的不同,活动轮廓模型主要分为参数活动轮廓模型和几何活动轮廓模型两种。其中参数活动轮廓模型直接以曲线或曲面的参数化形式来表达曲线或曲面的变形,这种表达形式允许和模型的直接交互,并且模型表达紧凑,利于模型的快速实时实现,但是难以处理模型的拓扑结构的变化,比如曲线的合并或分裂等。几何活动轮廓模型则是基于波前进化理论和水平集方法,将平面闭合曲线间接地表达为高1维水平集函数的零水平集的形式,再通过曲面的演化来隐含地求解曲线的演化。该模型能够自然地处理拓扑结构的变化,曲线的参数化仅仅是在模型变形后用于显示。

Snake模型[7]和测地线活动轮廓(GAC)模型[3]作为基于边界的活动轮廓模型的典型代表,其显著特点是利用由图像梯度矢量场构造出的边缘检测函数来引导曲线的演化。虽然对灰度不均图像具有较好的边缘识别、分割能力,但是对噪声较敏感。典型的基于区域的活动轮廓模型,例如C-V模型[8],考虑了图像灰度的全局特征,因此具有较好的抗噪声性能并且分割结果不受初始轮廓曲线位置的影响。此外,不少学者相继研究出了基于图像局部特征的区域型活动轮廓模型,Zhang等人[9]通过引入局部图像拟合(LIF)能量项获取演化曲线附近图像局部信息;Liu等人[10]在C-V能量函数上附加局部能量项,构造出一种局部C-V模型(LCV)。上述改进模型在较大程度上改善了传统C-V模型对强度非均匀图像的分割性能,但是在分割弱边界图像时灰度信息会对演化曲线的运动做出错误的引导,出现过分割或漏分割现象,致使演化曲线不能收敛到真实目标轮廓。

Klein和Amini等人[11]采用了Snake算法成功实现X射线血管造影图像中血管结构的分割提取。该方法首先通过人机交互的方式指定出待提取血管结构的初始轮廓,然后基于样条曲线的连续性及平滑性构造能量函数,通过极小化能量函数最终使得轮廓曲线逼近待提取血管的边缘。Lorigo等人[4, 12]提出了一种CURVES模型,利用1维曲线的曲率替代传统的平均曲率并正则化3维空间演化曲面,在3维管状结构血管提取上取得了较好的效果。Yan等人[5]提出了一种基于毛细管原理的活动轮廓模型。该模型在传统测地线活动轮廓(GAC)模型基础上增加湿壁能量项与体积能量项,提升算法对磁共振血管造影(MRA)图像细小血管分割的准确性。而Shang等人[6, 13]则结合血管灰度分布先验信息、多尺度血管矢量场和血管几何曲率特征,提出了基于血管几何特征的3维血管自动提取模型。

CT造影图像中往往存在噪声且背景极为复杂,其图像目标处像素具有非同质性,不同区域血管结构及灰度分布亦有差异性,仅考虑图像的灰度信息或梯度信息难以捕捉真实的目标轮廓。针对以上问题,首先利用改进的FCM[14]和C-V模型水平集提取出心脏CT图像感兴趣区域轮廓。接下来,由3维心脏图像数据获取多尺度梯度矢量信息构造边界型能量泛函,然后利用基于Hessian矩阵的多尺度血管函数对心脏感兴趣区域3维体数据增强滤波,获取血管先验形状信息用于约束能量泛函,最后融合边界、区域能量泛函并利用变分原理及水平集方法得到适合冠脉血管分割的水平集演化方程。由于血管图像的灰度不均匀,血管末端区域更为细小,所以上述算法的实施是面向划分为多个子区域的血管,在缩小的范围内进行轮廓的演化。最后通过对多组CT冠脉图像进行分割,验证了本文算法的有效性。

1 几何活动轮廓模型

由Caselles等人[3]提出的GAC模型作为基于边界的几何活动轮廓模型的典型代表,利用图像的梯度信息构建演化曲线的停止函数来引导轮廓曲线逼近真实目标边界。GAC模型通过极小化如下能量泛函,使得在分割目标处以g(·)为加权系数的闭合曲线最短,能量泛函为

${{E}^{GAC}}=\int_{0}^{L\left( C \right)}{g\left( \left| \nabla I\left[ C\left( s \right) \right] \right| \right)dS}$ (1)

式中,L(C)表示闭合轮廓曲线C(s)的弧长,Δ I为图像I的梯度,g(·)为边缘停止函数。

最小化式(1)对应的水平集演化方程为

$\frac{\partial \phi }{\partial t}=g\left( \left| \nabla I \right| \right)\left( div\left( \frac{\nabla \phi }{\left| \nabla \phi \right|} \right)+{{v}_{0}} \right)+\nabla g\cdot \nabla \phi $ (2)

式中,v0为常数,作用是加快曲线在图像平坦区域向内部收缩的速率,div为散度算子,g(·)是关于图像梯度模值 I的递减函数, g·Φ为增加的“吸引因子”项,用来对轮廓曲线的运动进行校正,提高模型对弱边缘目标边界分割能力,避免发生泄漏。通常边缘函数g(·)可定义为

$g\left( \nabla I \right)=\frac{1}{1+{{\left| \nabla I \right|}^{p}}},p=1或2$ (3)

式中,I=I*Gσ,Gσ为方差为σ的高斯核函数的。

从式(3)可以看出,边缘函数g(I)是来自单一尺度σ下的近似图像梯度。但选择的尺度过大,会使得两个空间相距较近的目标难以分离;反之,若尺度过小,对含有噪声的图像将存在大片图像区域梯度为零,使得演化曲线校正项为零,从而降低模型对若边界及含噪声图像的分割能力。

2 融合血管形状约束的活动轮廓模型

血管造影图像(如X射线、CT、磁共振成像等)受成像技术、外界噪声、造影剂分布不均以及组织病变等因素的影响,难以保证血管图像边界清晰、灰度分布均匀。因此仅依靠图像梯度信息实现血管结构的分割具有很大的局限性与缺陷,而且利用传统的活动轮廓模型分割形状复杂的血管方法,也难以获得满意的分割结果。所以在分割模型中引入血管形状信息具有重要意义。

2.1 模型介绍

对已提取3D心脏感兴趣区域,利用多尺度梯度场的Hessian矩阵获取血管先验信息函数Pvessel(x,y,z),进而获得基于区域的能量函数Evessel,构建融合边界和区域信息的能量函数。通过极小化该能量函数,利用水平集方法演化曲线,分割出血管真实轮廓。由于不同区域血管结构及灰度分布的差异性,以及血管尾端或枝节部分面积比例很小,对整个3D心脏区域数据采用固定的阈值参数难以保证较好的分割效果。所以先将血管尾端和枝节等分为几个子区域,在各子区域内进行独立的曲线演化后再融合。完整的3D心脏CT冠脉分割策略如图 1所示。

图 1 3D心脏CT冠脉分割示意图
Fig. 1 Segmentation of coronary artery in the 3D heart CT data

2.1.1 融合血管先验信息的能量函数

CT血管图像中由于造影剂的引入,冠脉血管呈现高亮度特性,背景为暗色。在血管中心线区域血管相似性函数取值接近于1,而在血管边缘及背景区域血管相似性函数取值为0,利用多尺度Frangi血管相似性函数可有效确定空间体素隶属于血管结构的可能性。定义血管先验信息函数为

$\begin{array}{l} {P_{vessel}}\left( {x,y,z} \right) = \\ \left\{ {\begin{array}{*{20}{c}} 1&{V\left( {x,y,z} \right) \ge a}\\ {V\left( {x,y,z} \right)}&{b \le V\left( {x,y,z} \right) < a}\\ { - \left( {1 - V\left( {x,y,z} \right)} \right)}&{V\left( {x,y,z} \right) < b} \end{array}} \right. \end{array}$ (4)

式中,a、b为控制血管先验信息函数函数对血管相似函数V(x,y,z)敏感性的阈值,a,b∈[0,1];(x,y,z)为3维体素坐标。Pvessel(x,y,z)∈[-1,1]为分段线性函数,当Pvessel取值接近于1时,说明该体素位于血管区域内;当Pvessel取值为接近-1时,说明体素点位于血管区域外;当Pvessel取值为远小于1的正数时,该体素点既可能在血管区域内,也可能在血管区域外。因此,根据函数Pvessel(x,y,z)的取值情况可有效判断体素点位于血管区域内的可能性,Pvessel(x,y,z)取值的正负可以分别控制演化函数的扩张与收缩。

由3D空间内血管区域信息构造血管能量函数,即

${{E}_{vessel}}=-\iiint_{R}{{{P}_{vessel}}\left( x,y,z \right)dxdydz}$ (5)

从式(5)可以看出,能量函数Evessel是基于3维体素区域项的,通过统计体数据内血管函数取值正负情况,寻找最优边界R使得能量函数Evessel最小。由于在血管区域内与区域外Pvessel(x,y,z)分别取正、负值,当演化曲线包围并逼近血管区域真实边界时,能量函数达到极小值。

利用变分法与梯度下降流对能量函数Evessel极小化,得到曲线演化方程为

$\frac{\partial C}{\partial t}={{P}_{vessel}}\left( x,y,z \right)\cdot N$ (6)

式中,N为演化曲线(或者曲面)的外法向量。当演化曲线(面)位于血管内部时,血管区域函数Pvessel(x,y,z)取值接近于1,亦即速度函数取值较大,这时会产生较大的扩张力,加速曲线的膨胀;当演化曲线(面)位于血管外部时,血管区域函数Pvessel(x,y,z)取值接近于-1,这时会加速曲线的演化的收缩。由基于区域项的能量函数Evessel引导的演化曲线在血管内外部膨胀、收缩运动如图 2所示。

图 2 曲线演化示意图
Fig. 2 Curve propagation

根据水平集演化原理可知,对于速度函数式Pvessel,水平集嵌入函数Φ按照演化公式为

$\frac{\partial \phi }{\partial t}={{P}_{vessel}}\left( x,y,z \right)\cdot \left| \nabla \phi \right|$ (7)

利用上述血管形状能量函数构建的区域型活动轮廓模型对X射线血管图像进行分割。如图 3所示,原始图像经过多尺度血管增强提取出的血管结构如图 3(b)。当分割血管设置的初始水平集演化曲线如图 3(c)中所示,血管信息函数Pvessel上下限阈值ab分别为0.95、0.05时,其对应的最终分割结果如图 3(d)所示;当初始水平集曲线位置如图 3(e),阈值ab取值分别为0.90、0.15时,其对应的最终分割结果如图 3(f)。可以看出,尽管在图像灰度分布不均,目标与背景灰度差异性较小且存在噪声的情况下,本文给出的基于血管先验形状的区域型活动轮依然能够逼近血管边缘。更重要的是该算法对初始轮廓位置不敏感,无论初始轮廓位于目标处或严重偏离目标均可实现对血管结构的分割。然而,该算法也存在一定的不足,主要包括以下两个方面:1)由于能量函数仅基于血管信息函数,而血管信息函数的取值受阈值ab的影响,导致模型对血管图像分割结果会因阈值选择的不同而存在差异;2)能量函数实质是通过统计水平集演化曲线内部、外部血管信息函数的取值的正负进而实现对血管结构的分割,因此轮廓曲线不具备较好的平滑性,曲线存在“毛疵”。

图 3 本文血管先验函数活动轮廓模型X射线血管图像分割
Fig. 3 X-ray vascular image segmented by the vessel prior energy function based GAC model ((a)X-ray vascular image; (b)Frangi vascular function enhanced result; (c)initial contour Ⅰ; (d)segmentation result Ⅰ; (e)initial contour Ⅱ; (f)segmentation result Ⅱ)

2.1.2 融合后模型的设计

定义总能量函数为

$E=\alpha {{E}_{vessel}}+\beta {{E}^{GAC}}$ (8)

通过式(1)(5)进一步定义能量函数为

$E=-\alpha \iiint_{R}{{{P}_{vessel}}}\left( x,y,z \right)dzdydz+\beta \int_{0}^{L\left( C \right)}{g}\left( \left| \nabla I\left[ C\left( s \right) \right] \right|ds \right)$ (9)

式中,αβ分别为分割过程中血管先验形状项和边界项的权值系数,均为正常数。当图像的灰度严重不均或存在较多噪声时,图像的形状约束信息应占主导,α取较大值;当图像目标内部梯度比较明显时,β取较大值。

由于水平集函数构造越平滑越好,使用符号距离函数作为水平集函数。在曲线演化过程中,为保证水平集函数为符号距离函数,需满足|Φ|=1,因此,在每一次水平集求解后需要重新初始化确保水平集函数满足以上条件,这就加大了计算量和运行时间。为解决水平集函数在每一次演化过程中都需要对水平集函数进行重新初始化的问题,在能量函数中引入了一个惩罚项[15]来驱使水平集函数接近于符号距离函数。总能量泛函为

$\begin{align} & E=-\alpha \iiint_{R}{{{P}_{vessel}}\left( x,y,z \right)dxdydz}+\beta \int_{0}^{L\left( C \right)}{g}\left( \left| \nabla I\left[ C\left( s \right) \right] \right|ds \right)+ \\ & \ \ \ \ \ \ \ \ \gamma \int_{\Omega }{0.5}{{\left( \left| \nabla \phi \right|-1 \right)}^{2}}d\Omega \\ \end{align}$ (10)

利用变分原理及梯度下降流极小化总能量泛函,推导水平集演化方程为

$\begin{array}{l} \frac{{\partial \phi }}{{\partial t}} = - \alpha {P_{vessel}}\left( {x,y,z} \right) \cdot \left| {\nabla \phi } \right| + \\ \beta \left[ {g\left( {\left| V \right|} \right)\left( {div\left( {\frac{{\nabla \phi }}{{\left| {\nabla \phi } \right|}}} \right) + {v_0}} \right) + \nabla g \cdot \nabla \phi } \right] + \\ \gamma \left[ {\Delta \phi - div\left( {\frac{{\nabla \phi }}{{\left| {\nabla \phi } \right|}}} \right)} \right] \end{array}$ (11)

2.2 形状约束下的血管区域增强

2.2.1 图像的Hessian矩阵

血管的3D影像呈现一种局部管状或线状特性,判断空间体素点是否属于血管,需要考察该点的局部特征。在数学上,可以用泰勒展开式来分析影像局部特征,对于二阶(Hessian 矩阵)项来讲,当梯度很小时仍然可以表现出局部结构信息。图像在像素点x0周围邻域内的泰勒展开式可描述为

$I\left( {{x_0} + \Delta x,\sigma } \right) \approx I\left( {{x_0},\sigma } \right) + \Delta {x^T}\nabla I\left( {{x_0},\sigma } \right) + \Delta {x^T}H\left( {{x_0},\sigma } \right)\Delta x$ (12)

式中,I(x0,σ)和H(x0,σ)分别表示图像在像素点x0处选取的尺度为σ时计算出来的梯度向量和Hessian矩阵;图像数据点x0=(x,y) (2D图像)或x0=(x,y,z) (3D图像)。在给定尺度σ条件下,有

$\begin{array}{l} G\left( {x,y,z,\sigma } \right) = \frac{1}{{\sqrt {2\pi \sigma } }} \times \\ \exp \left( { - \frac{{\left\| {{x^2} + {y^2} + {z^2}} \right\|}}{{2{\sigma ^2}}}} \right) \end{array}$ (13)

${I_{ab}}\left( {x,y,z;\sigma } \right) = \frac{{{\partial ^2}}}{{\partial a\partial b}}G\left( {x,y,z;\sigma } \right)*I\left( {x,y,z} \right)$ (14)

式中,I(x,y,z)表示图像在体素点(x,y,z)处的灰度值,G(x,y,z,σ)为高斯函数。空间某体素点x0处的Hessian矩阵可以由图像I在该点的二阶偏导数求得,对于3维体素点x0

$H\left( {{x_0}} \right) = \left[ {\begin{array}{*{20}{c}} {{I_{xx}}}&{{I_{xy}}}&{{I_{xz}}}\\ {{I_{yx}}}&{{I_{yy}}}&{{I_{yz}}}\\ {{I_{zx}}}&{{I_{zy}}}&{{I_{zz}}} \end{array}} \right]$ (15)

式中,IxxIxy等分别表示体素点处的二阶偏导数。设λk为Hessian矩阵的第k个特征值,vk为对应的特征向量,按照特征值和特征向量的定义得出

$H\left( {{x_0}} \right){v_k} = {\lambda _k}{v_k} \Rightarrow v_k^TH\left( {{x_0}} \right){v_k} = {\lambda _k}$ (16)

由式(16)可知,图像点x0处Hessian 矩阵的特征向量互相正交并且在空间内相互垂直。因此可以用Hessian 矩阵的特征值和特征向量来描述图像点二阶导数的大小和方向。

2.2.2 血管增强

Frangi在文献[16]中首次提出的基于Hessian矩阵特征值的血管结构判断方法,通过分析矩阵特征值λσ,k(k=1,2,3;σ为对应的尺度)及与其对应的特征向量vσ,k来判断体素是否为血管区域并确定血管的方向,达到增强管状结构、抑制其他结构的目的。对于3D血管影像,λσ,k中最大的特征值对应的特征向量,表示此点为曲率变化最大的方向,与血管走向垂直;最小的特征值对应的特征向量表示曲率变化最小的方向,为血管真实走向。设最大特征值λσ,1对应的特征向量为vσ,1,则各特征值与血管结构关系如图 4所示。

图 4 Hessian矩阵特征值与血管结构关系示意图
Fig. 4 Relationship between Hessian matrix eigenvalue and vascular structures

在理想管状结构处有,|λσ,1|≈0,|λσ,1|<<|λσ,2|,λσ,2λσ,3,血管模式与Hessian矩阵特征值之间的关系如表 1所示。据此,Frangi定义多尺度Hessian矩阵3D血管增强函数为

表 1 2D/3D图像几何结构与Hessian矩阵特征值(大小及符号)对应关系
Table 1 Relationship between 2D/3D Image geometryand Hessian matrix eigenvalues

下载CSV
2D 3D 几何结构
λ1 λ2 λ1 λ2 λ3
N N N N N N背景/噪声
L L H- 盘状结构(bright)
L L H+ 盘状结构(dark)
L H- L H- H- 管状结构(bright)
L H+ L H+ H+ 管状结构(dark)
H- H- H- H- H- 团状结构(bright)
H+ H+ H+ H+ H+ 团状结构(dark)
注:H=high, L=low, N=noise, +/-代表特征值的符号。

$\begin{array}{l} V\left( {x,y,z,\sigma } \right) = \\ \left\{ \begin{array}{l} 0\;\;\;\;\;\;\;{\lambda _2} \ge 0,{\lambda _3} \ge 0\\ \left[ \begin{array}{l} \left[ {1 - \exp \left( { - \frac{{{R_a}}}{{2{\alpha ^2}}}} \right)} \right] \times \\ \exp \left( { - \frac{{{R_b}}}{{2{\beta ^2}}}} \right) \times \\ \left[ {1 - \exp \left( { - \frac{{{S^2}}}{{2{\gamma ^2}}}} \right)} \right] \end{array} \right]\;\;\;\;其他 \end{array} \right. \end{array}$ (17)

式中,Ra、RbS分别用来抑制片状、点状及背景噪声,$ {R_a} = \frac{{\left| {{\lambda _2}} \right|}}{{\left| {{\lambda _3}} \right|}}$${R_b} = \frac{{\left| {{\lambda _1}} \right|}}{{\sqrt {\left| {{\lambda _2}{\lambda _3}} \right|} }}$$S = \sqrt {\lambda _1^2 + \lambda _2^2 + \lambda _3^2} $α、βγ 分别表征对测度其敏感性经验阈值,V(x,y,z,σ)∈[0,1],在血管结构处其取值接近于1,而在非血管结构处取值趋近于0。

分别以2维X射线及CT血管造影图像为例,对原始血管图像5(a)(c)利用多尺度Hessian矩阵增强后的图像分别如图 5(b)(d)所示。在血管内部血管函数取值接近于1,增强后的图像亮度较高;在血管外部血管函数值接近于0,增强后的亮度较低。

图 5 Frangi多尺度血管增强
Fig. 5 Frangi multi-scale vascular enhance ((a) X-ray vascular image; (b)vascular enhanced result of fig.5(a); (c)CTA vascular image; (d)vascular enhanced result of fig.5(c))

2.2.3 多尺度矢量场边界活动轮廓模型

借鉴文献[17]中提出的基于多尺度梯度矢量场的边界活动轮廓模型,用多尺度梯度矢量V取代GAC模型中单一尺度下平滑图像的梯度矢量I,改善模型对模糊边界分割能力的同时有效分离两相邻目标。

基于多尺度矢量场的边界活动轮廓演化方程为

$\frac{{\partial \varphi }}{{\partial t}} = g\left( {\left| V \right|} \right)\left( {div\left( {\frac{{\nabla \varphi }}{{\left| {\nabla \varphi } \right|}}} \right) + {v_0}} \right) + \nabla g \cdot \nabla \varphi $ (18)

式中,V为多尺度矢量场,当曲线按照该模型演化时,受到多尺度梯度场V所提供的标量场g(V)g(V)的共同作用,使曲线C朝着图像中目标的边缘靠近,最终稳定在边缘处。

图 6是采用多尺度矢量场GAC模型对心脏CT切片图像右冠状动脉提取的结果,其中V是尺度σ=0,1,2,3下依次递进地进行高斯平滑后的梯度矢量场的综合。在相同的初始轮廓下,传统GAC模型迭代200次后已发生边界泄漏问题;而多尺度下的GAC模型对模糊边界具有较好的鲁棒性,迭代收敛至右冠脉区域真实边界。实验结果表明,利用多尺度矢量场替代单一的矢量场能够改善传统GAC模型对医学图像边界提取能力。

图 6 心脏CT切片右冠脉分割结果
Fig. 6 Segmentation results of right coronary artery in the heart CT slice image ((a)image of heart region; (b)initial curve; (c)segmentation result of the traditional GAC model; (d)segmentation result of the multiscale GAC model)

3 仿真实验及结果分析

将本文算法用于双源CT心脏冠脉图像分割实验。实验中8组心脏CT图像序列均由南京某医院提供,采样仪器为西门子SOMATOM Definition 双源CT心脏扫描机,图像在水平面上的像素大小为0.50 0.92 mm,层厚0.6251.00 mm,空间分辨率为0.33×0.33 mm,扫描层数220522层。

3.1 参数设置

实验中以下参数可以设置为相同的数值:迭代时间步长τ=0.1,演化常数v0=1,能量惩罚项系数γ=1。血管区域信息函数Pvessel上下限阈值αβ∈(0,1)的取值根据血管图像模糊程度、噪声及血管大小进行设定。一般情况下当图像背景清晰、血管结构明显时,多尺度血管增强获得的血管相似函数V的值更均匀(血管内部像素点对应的V的取值几乎都趋近于1,而血管外部V的取值接近于0),此时α可取较大的小数(如0.98),β取较小的值(如0.05);反之,α取值应适当减小,β适当增大。边界项和血管先验形状项的权值系数αβ,依据待分割图像特征、质量等因素来确定。因为血管先验项引导的曲线演化往往比由梯度信息控制的曲线演化效率更高,尤其是在图像存在噪声等干扰时。在大多数情况下,α≥β,且α取值应大于0.5。更重要的是较大的α可以使得活动轮廓模型对初始化曲线的位置不敏感,其原因在于血管先验项能量函数利用血管函数获取血管位置,能够产生一种“约束力”,使得演化曲线固定在血管周围,最终在与梯度项能量函数共同作用下逼近血管结构真实轮廓。

3.2 3D冠脉分割实验

实验1 实验中的3维CT心脏图像大小为1 028×878×276像素,像素大小为0.50 mm,层厚1.00 mm。在心脏区域3D体数据中分出6个子区域,选定各子区域中左、右冠脉水平集初始轮廓,利用本文模型进行分割实验。其中图 7(a)(b)分别为右、左冠脉区域水平集迭代50次分割结果,当水平集迭代收敛后分割出完整的右冠状动脉及左管状动脉分别如图 7(c)(d)

图 7 本文算法冠脉分割结果
Fig. 7 Coronary artery segmentation results of the proposed method ((a)right coronary segmentation result with 50 iterations; (b)left coronary segmentation result with 50 iterations; (c) segmentation result of right coronary with iteration convergence;(d)segmentation result of right coronary with iteration convergence)

图 8为右冠脉血管末端子区域分割切片的部分结果,可以看出,血管在这一区域已经变的很细小,若对整个3D血管采用同一曲线进行演化,这一段将无法分割出来。

图 8 右冠脉分割切片显示结果
Fig. 8 Segmentation results of right coronary artery ((a)the 40th slice;(b)the 50th slice; (c)the 60th slice; (d)the 70th slice;(e)the 80th slice;(f)the 90th slice)

图 9为对比实验结果,其中图 9(a)为本文算法提取出的完整冠状动脉结构;图 9(b)为采用传统区域生长算法分割出的冠脉结构,受算法的局限性及阈值设置不当的问题,已出现严重的过分割和欠分割的情况;图 9(c)为采用文献[18]区域边界混合模型分割的结果,该算法的能量函数虽然融合了区域项与边界项,但是由于其区域项仅依靠灰度强度阈值作为引导,并没有引入冠脉先验知识作为水平集演化约束,因此分割结果中也包含了像素较高的升主动脉、心房、心室等结构,且无法分割出灰度不均的冠脉末端。由实验结果可知,本文算法未发生分割泄漏问题,且人工交互简单,只需手动选定左右冠脉分割初始轮廓曲面位置即可,这表明了本文算法在双源CT冠脉图像分割方面具备较好的正确率与优越性。

图 9 冠脉分割对比结果
Fig. 9 Contrast results of coronary segmentation ((a) segmentation result of proposed method; (b) segmentation result of region growing method; (c) segmentation result of the hybrid boundary and region model)

实验2 为定量评价本文冠脉分割算法的性能,从所有冠脉CT数据中随机选取5组,利用Mimics软件,手动选择种子点和阈值,分割出全部3维冠脉血管,将其分割结果作为标准参考图像(Ground Truth)与本文提出的冠脉分割算法获得的结果相比较。目前,国内外对血管分割效果较好的算法和软件是Mimics软件,但利用该软件进行分割时,需要手动选择的种子点数目较多,(心脏右冠脉为3040个),阈值的选择也需要有经验的医师进行。图 10(b)为Mimics软件手动分割的效果,图 10(a)为利用本文算法分割出的冠状动脉的示例。可以看出本文方法丢失了血管部分小的枝节,但未发生分割泄漏问题,且人工交互简单,只需手动选定左右冠脉分割初始轮廓曲面位置即可,表明本文算法在双源CT冠脉图像自动分割方面具有优越性,同时也保持了较高的正分割确率。

图 10 手动分割与本文方法分割结果示例
Fig. 10 Examples of both manual segmentation and the proposed method ((a)segmentation result of proposed method; (b)manual segmentation result of manual)

实验中,3维体素最终被分割为血管体素或非血管体素,可以用4个事件来表述:

1)真阳性(TP)为血管点被正确分割为血管点的体素数目;

2)假阳性(FP)为血管点被错分割为非血管点的体素数目;

3)真阴性(TN)为非血管点被正确分割为非血管点的体素数目;

4)假阴性(FN)为非血管点被错分割为血管点的体素数目。

采用基于体素统计的测度来评价本文分割算法的性能,分别为正确率、灵敏性和特异性,计算公式为

$\begin{array}{l} acc = \frac{{TP + TN}}{{TP + FP + TN + FN}}\\ sen = \frac{{TP}}{{TP + FN}}\\ spe = \frac{{TN}}{{TN + FP}} \end{array}$ (19)

表 2为本文算法冠脉分割结果与医师利用软件手动分割结果的对比。从表 2中可以看出,本文算法的特异性明显高于正确率及敏感性,表明将非血管区域错误分类到血管区域(即发生分割泄露)的情况相对较少;第4组与第5组数据分割结果各项指标均高于前3组,这也间接说明分割算法的性能与待分割心脏CT图像大小及质量等因素有关,图像分辨率越高、造影越明显分割效果越好。从整体上看,本文冠脉血管分割算法各项指标都较高,基本上满足冠脉分割的要求。

表 2 本文算法与手动分割结果对比
Table 2 Comparison results of the proposed method and manual segmentation

下载CSV
数据尺寸/像素正确率敏感性特异性
第1组512×512×2750.874 10.875 20.942 5
第2组512×512×2860. 866 50.853 20.885 5
第3组512×512×2690.797 30.840 50.931 5
第4组1 028×878×2660.914 50.889 50.943 7
第5组1 028×878×2790.920 10.903 60.941 9
平均 0.874 50.872 40.927 0

4 结 论

针对心脏CT图像血管结构的复杂性,提出了基于形状约束的活动轮廓模型血管分割方法,利用Frangi多尺度血管函数对心脏感兴趣区域增强滤波获取血管先验形状信息构造区域项能量函数,由原始图像的梯度矢量场构造边界项能量函数。该方法充分融合血管图像的灰度信息及梯度场信息,克服了传统方法在血管分割上的不足。针对冠脉血管结构复杂,末端存在诸多细小枝节的情况,对整体血管进行了分区,曲线演化是针对各个子区域独立完成的。实验结果表明本文方法能够有效提取3D血管各部分,且仅需给出初始的轮廓线,人工交互简单。对于血管细枝末节的提取,与手动分割效果相比,不是很完美。进一步,考虑可以根据轮廓演化得到的细小血管初始区域,自动获取种子点,结合区域生长方法达到对其精准分割的目的。

参考文献

  • [1] Kiri?li H A, Schaap M, Metz C T, et al. Standardized evaluation framework for evaluating coronary artery stenosis detection, stenosis quantification and lumen segmentation algorithms in computed tomography angiography[J]. Medical Image Analysis, 2013 ,17 (8) : 859 –876. [DOI:10.1016/j.media.2013.05.007]
  • [2] Montalescot G, Sechtem U, Achenbach S, et al. 2013 ESC guidelines on the management of stable coronary artery disease[J]. European Heart Journal, 2013 ,34 (38) : 2949 –3003. [DOI:10.1093/eurheartj/eht296]
  • [3] Caselles V, Kimmel R, Sapiro G. Geodesic active contours[J]. International Journal of Computer Vision, 1997 ,22 (1) : 61 –79. [DOI:10.1023/A:1007979827043]
  • [4] Lorigo L M, Faugeras O D, Grimson W E L, et al. Curves: curve evolution for vessel segmentation[J]. Medical Image Analysis, 2001 ,5 (3) : 195 –206. [DOI:10.1016/S1361-8415(01)00040-8]
  • [5] Yan P K, Kassim A A. Segmentation of volumetric MRA images by using capillary active contour[J]. Medical Image Analysis, 2006 ,10 (3) : 317 –329. [DOI:10.1016/j.media.2005.12.002]
  • [6] Shang Y F, Wang H, Wang N, et al. Three dimensional vessel extraction model based on tubular characters and active contour model[J]. Journal of Image and Graphics, 2013 ,18 (3) : 290 –298. [ 尚岩峰, 汪辉, 汪宁, 等. 管状特性和主动轮廓的3维血管自动提取[J]. 中国图象图形学报, 2013 ,18 (3) : 290 –298.] [DOI:10.11834/jig.20130307]
  • [7] Kass M, Witkin A, Terzopoulos D. Snakes: active contour models[J]. International Journal of Computer Vision, 1988 ,1 (4) : 321 –331. [DOI:10.1007/BF00133570]
  • [8] Chan T F, Vese L A. Active contours without edges[J]. IEEE Transactions on Image Processing, 2001 ,10 (2) : 266 –277. [DOI:10.1109/83.902291]
  • [9] Zhang K H, Song H H, Zhang L. Active contours driven by local image fitting energy[J]. Pattern Recognition, 2010 ,43 (4) : 1199 –1206. [DOI:10.1016/j.patcog.2009.10.010]
  • [10] Liu S G, Peng Y L. A local region-based Chan-Vese model for image segmentation[J]. Pattern Recognition, 2012 ,45 (7) : 2769 –2779. [DOI:10.1016/j.patcog.2011.11.019]
  • [11] Klein A K, Lee F, Amini A A. Quantitative coronary angiography with deformable spline models[J]. IEEE Transactions on Medical Imaging, 1997 ,16 (5) : 468 –482. [DOI:10.1109/42.640737]
  • [12] Lorigo L M, Faugeras O, Grimson W E L, et al. Codimension-two geodesic active contours for the segmentation of tubular structures[C]//Proceedings of IEEE Conference on Computer Vision and Pattern Recognition. Hilton Head Island, SC, USA: IEEE, 2000, 1: 444-451. [DOI: 10.1109/CVPR.2000.855853]
  • [13] Shang Y F, Deklerck R, Nyssen E, et al. Vascular active contour for vessel tree segmentation[J]. IEEE Transactions on Biomedical Engineering, 2011 ,58 (4) : 1023 –1032. [DOI:10.1109/TBME.2010.2097596]
  • [14] Li B N, Chui C K, Chang S, et al. Integrating spatial fuzzy clustering with level set methods for automated medical image segmentation[J]. Computers in Biology and Medicine, 2011 ,41 (1) : 1 –10. [DOI:10.1016/j.compbiomed.2010.10.007]
  • [15] Li C M, Xu C Y, Gui C F, et al. Level set evolution without re-initialization: a new variational formulation[C]//Proceedings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition. San Diego, California: IEEE, 2005, 1: 430-436. [DOI: 10.1109/CVPR.2005.213]
  • [16] Frangi A F, Niessen W J, Vincken K L, et al. Multiscale vessel enhancement filtering[C]//Proceedings of 1st International Conference Cambridge. Berlin Heidelberg: Springer, 1998: 130-137. [DOI: 10.1007/BFb0056195]
  • [17] Peng J Y, Hao C Y, Qi H, et al. MR image segmentation based on GAC model with multiscale gradient vector field[J]. Journal of Image and Graphics, 2007 ,12 (7) : 1214 –1217. [ 彭进业, 郝重阳, 齐华, 等. 基于多尺度梯度矢量场GAC模型的MR医学图像分割[J]. 中国图象图形学报, 2007 ,12 (7) : 1214 –1217.] [DOI:10.11834/jig.20070723]
  • [18] Zhang Y, Matuszewski B J, Shark L K, et al. Medical image segmentation using new hybrid level-set method[C]//Proceedings of the 5th International Conference on BioMedical Visualization. London: IEEE, 2008: 71-76. [DOI: 10.1109/MediVis.2008.12]