基于snake模型的IVUS图像序列三维分割方法

孙 正， 杨 宇

（华北电力大学电子与通信工程系，河北 保定 071003）

近年来，血管内超声(Intravascular Ultrasound,IVUS)成像技术已经越来越广泛地被应用到血管疾病，特别是冠心病的临床诊断和介入治疗中。它不仅能够反映血管内腔的变化，同时也能反映含斑块在内的管腔横断面结构、血管壁厚度、形态以及斑块成分等[1]。血管壁在IVUS图像上呈现回声不同的三层环形结构：内层为一强回声亮环，在组织学上由内膜及内弹性膜组成；中层相当于中膜，呈现为低回声暗区；外层为较亮的强回声带，由外弹性膜和外膜构成。图像序列中的前后帧之间非常相似，具有很强的相关性。

对IVUS图像进行分割，提取出血管内腔边界和中-外膜边界（包括可能存在的斑块边缘），是定量分析（包括测量血管腔径、容积和斑块体积、估计狭窄百分比、测量血流动力学参数等）[2]和血管三维重建[3]的重要步骤，分割质量直接决定定量测量和三维重建的精度。目前临床常采用手动分割，由富有经验的医生用鼠标在屏幕上逐帧勾画。由于IVUS图像序列常包括上千帧，因而该方法不但耗时，工作量大，而且可重复性差，分割结果在很大程度上取决于操作者的经验和主观因素。若采用数字图像处理的方法进行自动提取，则由于IVUS图像一般受噪声污染严重（包括斑点噪声、图像伪影和部分血管壁钙化阴影等的影响），常规的灰度图像轮廓提取方法，如微分法等，难以获得满意结果。已经提出的半自动提取方法主要包括图搜索[4]、活动轮廓(snake)模型[5]及其相关改进算法[6]等。这些方法的分割效果较常规方法有很大提高，但仍需操作者一定程度的手动参与，从而延长处理时间或引入误差。全自动的分割方法目前已成为该领域的研究热点，例如，Brusseau等[7]的基于统计活动轮廓模型的方法、Filho等[8]的基于模糊聚类的方法、Giannoglou等[9]的自动获得初始轮廓的 snake方法、Ruiz等[10]的概率统计分割内腔方法、Papadogiorgaki 等[11]的基于纹理分析的方法、Bucher等[12]的形状驱动方法、Bovenkamp等[13]的多智能体方法、Zhang等[14]的基于 contourlet变换的snake方法等。上述方法都是在二维空间中进行串行处理，即对IVUS序列中的各帧图像进行逐帧分割。近年来，三维分割技术[15]成为该领域新的研究方向。三维分割与二维分割相比具有以下优点：第一，可以利用整个序列的信息；第二，可同时完成对各帧图像的分割，实现对整个序列的并行处理，从而大大提高处理效率，缩短处理时间。

本文提出一种基于snake模型的IVUS图像序列三维分割方法。在对原始图像进行滤波去噪和去除环晕伪像的预处理的前提下，首先获取IVUS序列的四个纵向视图，并从中提取出内腔和中-外膜边界；然后，将各纵向视图中的边界曲线映射到各IVUS帧中，得到各横向视图中的初始轮廓；将此初始轮廓作为snake模型的初始形状，通过使预先设定的能量函数最小，模型在内外力的共同作用下不断变形，最终得到各IVUS帧中的内腔边界和中-外膜边界，从而实现对各帧图像的并行分割。

1 方 法

如图1所示，本文的分割方法包括以下步骤：① 原始图像的预处理，包括滤波去噪和去除环晕伪像；② 获取 IVUS图像序列的四个纵向视图；③ 提取纵向视图中的内腔边界和中-外膜边界；④ 将步骤③中获得的边界曲线映射到各帧IVUS切片图像中，得到各 IVUS帧的初始内腔边界和中-外膜边界；⑤ 在各IVUS帧中，初始轮廓演化变形，得到最终的内腔和中-外膜轮廓。

图1 IVUS图像序列三维分割方法流程图

1.1 图像预处理

首先，采用中值滤波和高斯平滑两种通用预处理方法，减少IVUS图像中的椒盐噪声和随机噪声。

在IVUS成像过程中，超声导管在图像中心区域产生一个无回声暗区，围绕导管黑色影像的厚薄不一的白色环状影像称为环晕伪像(ring-down artifact)[16]，即图2中箭头所指。它与传感器声波振荡引起的高振幅信号有关，不是导管本身的影像，应该在预处理阶段加以去除，否则会影响图像分割的准确度。根据环晕伪像在IVUS极坐标视图中的分布特性，本文采用极坐标变换法对其加以去除。

图2 环晕伪像在两种视图上表示

首先，对IVUS图像进行极坐标变换，得到其极坐标视图，如图2(b)所示。可见，环晕伪像固定的位于极坐标视图的上部。采用阈值法，按照下式去除极坐标视图中的环晕伪像

其中 I( r,θ)和I'(r,θ)分别为原极坐标视图和去除环晕伪像后的极坐标视图中像素点(r,θ)处的灰度值；Icatheter是原 IVUS图像中导管区域的像素灰度值；r为像素点的极径；t=ε×ImageHeight 2；ε为权重参数，经实验确定其取值范围为[0.1,0.35]。在极坐标视图中消除环晕伪像后，再经过极坐标逆变换，即可得到直角坐标系下去除环晕伪像后的IVUS图像。

1.2 获得IVUS序列纵向视图

采集IVUS图像时，随着超声导管的缓慢回撤，获得的是一系列血管横截面切片图像。IVUS纵向视图是指IVUS图像序列沿平行于导管方向的纵向截面视图。如图3所示，本文沿血管长轴（纵向）方向取四个纵向视图，即垂直切面A、水平切面B、左对角线切面C和右对角线切面D。

图3 获取IVUS纵向视图

1.3 纵向视图中血管壁边缘的提取

如图4所示，IVUS纵向视图中的每一行像素即为一帧切片图像与该纵向视图的交线。本文根据纵向视图中的内腔边界和中-外膜边界的灰度同其周围像素灰度相差较大的特征，采用灰度阈值法进行提取。算法的基本原理是：遍历纵向视图中的各像素，判断识别目标边界点。由于纵向视图中的目标边界分布于图像中部，所以从图像中轴线开始，分别向左和向右进行逐行遍历，如图 5所示。对于当前像素(i, j)，其灰度值为I( i, j)，若 I ( i, j + 1 )- I ( i, j) ≥η ，η为阈值，则为目标边界点，否则不是。将每行左右两部分的像素中，第一个符合上述条件的像素记为内腔边界点，第二个记为中-外膜边界点。经实验确定η的取值区间为[10, 20]。

图4 纵向与横向视图的对应关系示意图

图5 逐行遍历纵向视图像素

1.4 横向视图中初始轮廓的获取

纵向视图的轮廓线对应到横向视图中为轮廓点，四条纵向视图轮廓线映射到横向视图中为四个轮廓点，如图6所示。在得到各帧横向视图中的内腔和中-外膜初始边界点以后，采用DDA(Digital Differential Analyzer)算法[17]依次连接各初始边界点，如图7所示，同时获得各帧图像中的初始内腔边界和中-外膜边界。

1.5 基于snake的初始轮廓演化

snake模型是Kass等[5]在1989年提出的一个能量最小化样条模型，它通过将由图像数据直接获得的约束力和关于目标的先验知识结合起来，完成对图像中目标的分割、匹配和跟踪[18]。本文将1.4节中获得的各帧IVUS图像中血管壁的初始轮廓作为snake模型的初始形状，通过使预先设定的能量函数最小，snake不断变形，当能量函数取得全局最小值时，snake模型即停留在目标轮廓处，从而完成对各帧图像中内腔和中-外膜边界的并行提取。

图6 纵向与横向视图轮廓对应关系示意图

图7 由初始点获取初始轮廓

将初始轮廓离散成由 N个点组成的有序点集，则snake能量函数的离散表达式为

其中 内部能量Eint使 snake在变形过程中保持连续和平滑，其归一化表达式为

其中 ci( xi, yi)(i=1,2,…, N-2)是第 i个 snake点；和maxd分别是相邻snake点之间的平均距离和最大距离，在每次迭代结束时，它们的值都被更新；α和β是权重参数。式(3)中的第一项保证变形过程中snake的连续性，使snake点均匀分布，不致产生收缩的现象；第二项是曲线二阶导数的离散形式，保证变形过程中snake的光滑性。由于进行了归一化，故式(3)中两项的取值范围都在[0, 1]区间内，权重α和β的取值区间也是[0, 1]。

外部能量Eext是保证snake收敛的外部力，它决定snake的移动方向。本文将其定义为

其中 I( xi, yi)和 ∇ I ( xi, yi)分别是像素(xi, yi)的灰度和灰度梯度值。由于8位灰度图像的灰度和灰度梯度的取值范围分别为[0, 255]和别进行了归一化，使其取值范围为[0, 1]区间。γ, λ ∈ [ 0,1]是权重参数。

通过使能量函数E最小，snake在内外力的共同作用下从初始形状开始不断变形，最终停留在能量函数取得全局最小值的最优位置，即为目标轮廓。能量函数的全局最优采用Williams快速算法[19]完成，不仅计算速度快，所需存贮空间小，而且非常适合于编程实现。

2 实验结果及分析

2.1 实验结果

本文采用临床采集的冠状动脉IVUS图像序列对所提出的分割方法进行了实验验证，在这一节中给出部分实验结果。实验序列是采用Jomed Endosonic超声成像仪在临床获取的，超声探头为2.9F Jomed单轨机械探头，频率为30MHz。在回撤导管的过程中，探头导管以1800转/分作360°旋转时可连续获得30帧/秒的血管横轴实时切面图像。采用Visual C++6.0编程实现算法，所用计算机配置为2.6GHz Intel Pentium III Xeon Tanner processor、2 GB内存。

图8和图9是对IVUS图像去除环晕伪像的实验结果。图10是从一个IVUS序列中选出的7帧，均已完成了对环晕伪像的去除。图 11是从该序列获取的四个纵向视图，其中提取出了血管壁的内外膜边缘；将其映射到各帧切片图像上，即可得到血管壁内外膜的初始轮廓点，如图12(a)～(d)所示，然后采用 DDA算法用直线段依次连接各轮廓点即可得到初始轮廓，如图12(e)～(h)所示。之后，通过使 snake能量函数最小，初始轮廓不断变形，最终停留在使能量函数取得全局最小值的地方，即内腔和中-外膜边界，如图13和图14所示，其中式(3)和式(4)中的权重参数设定为 α=β=γ=λ=1，即各能量项的权重相等。

图8 一帧IVUS图像的环晕伪像去除结果(ε = 0.21)

图9 四帧IVUS图像的环晕伪像去除结果(ε=0.21)

图10 去除环晕伪像的7帧IVUS图像

图11 完成了管壁轮廓提取的IVUS序列纵向视图

图12 IVUS横向视图中血管壁的初始轮廓 (a)～(d)初始轮廓点；(e)～(h)初始轮廓

图13 IVUS横向视图中初始轮廓的变形结果1

图14 IVUS横向视图中初始轮廓的变形结果2

2.2 分割精度和速度

目前还没有一种公认的客观、定量评价IVUS图像分割方法精度的体系或标准，一般采用将自动分割结果与医学专家手动分割结果进行比较的方法来衡量自动方法的分割质量[20]。作者从临床获取的冠状动脉内超声图像序列中随机选取了一个序列，由200帧组成，由医生手动描绘出各帧图像中的内腔边界和中-外膜边界，并将其作为“金标准”。然后，分别采用本文方法和二维方法对图像进行自动分割，并记录处理时间。对各帧图像，计算自动分割出的轮廓与“金标准”之间的欧氏距离作为误差[21]，定量估计分割算法的精度。其中二维方法采用 Giannoglou等[9]提出的自动获取初始轮廓的 snake模型法。采用两种方法进行分割之前，均对原始图像进行1.1节中所述的预处理，即滤波去噪和去除环晕伪像（参数设置为 ε=0.21）。处理时间均是在2.6GHz Intel Pentium III Xeon Tanner processor、2GB内存的计算机参数下测得的，且均不包含预处理（去噪和抑制伪影）的时间。

对200帧图像的分割误差以及处理时间的统计结果如表1所列。由表可知，无论对内腔边界还是中-外膜边界，本文方法的平均误差与二维方法均十分接近。而且，对内腔边界的分割结果均明显优于中-外膜边界。这是由于IVUS图像中血管内腔与内膜的灰度值差异通常比较明显，因而在去除环晕伪像的情况下，便可较准确地提取出内腔边界。而中-外膜边界的精确提取依赖于图像中中膜的形态及初始轮廓的位置，若中膜在图像中表现为较厚的低灰度带，且初始轮廓较接近真实轮廓，则最终分割结果较好，否则易被伪边缘吸引[22]。对于处理时间，二维方法是进行逐帧处理的串行分割，对每帧图像的处理时间约为0.26s，则处理200帧共需约52s。而采用本文方法，由于对各帧的分割是同时进行的，因而对200帧图像的处理时间仅需 16.8s，从而大大提高了处理效率。大大提高处理效率，缩短处理时间。

表1 对200帧图像的分割误差及处理时间统计结果

3 结 论

对IVUS图像进行分割，提取出管腔横截面图像中的内腔边界和中-外膜边界，是对临床重要的血管形态参数进行定量测量和血管三维重建等的重要步骤，对血管疾病的临床诊治具有重要意义。本文提出了一种基于snake模型的IVUS图像序列的三维分割方法，不仅可充分利用整个图像序列的信息，而且可实现对图像序列的自动、并行分割。该方法在对原始图像进行去噪和抑制环晕伪像等预处理的基础上，将从四个纵向视图中提取出的血管壁边界映射到各帧IVUS横向视图中，获得各帧图像中的初始轮廓；然后通过使预先设定的snake能量函数最小，初始轮廓不断演化得到最终轮廓。实验结果表明，在获得与现有逐帧处理的二维串行分割方法相类似的分割精度的情况下，采用本文的三维分割方法可

[1]杨青苗, 吕树铮. 血管内超声在冠心病诊断中的临床应用进展[J]. 心肺血管病杂志, 2010, (01):73-76.

[2]Baldewsing RA, Danilouchkine MG, Mastik F, et al.An inverse method for imaging the local elasticity of atherosclerotic coronary plaques [J]. IEEE Transactions on Information Technology in Biomedicine, 2008, 12(3): 277-289.

[3]孙 正. 应用血管内超声与X射线造影图像融合的血管三维重建[J]. 工程图学学报, 2010, 31(1):116-123.

[4]Olszewski ME, Wahle A, Vigmostad SC, et al.Multidimensional segmentation of coronary intravascular ultrasound images using knowledgebased methods [C]//Proceedings of SPIE Conference on Medical Imaging and Image Processing, 2005,5747: 496-504.

[5]Kass M, Witkin A, Terzopoulos D. Snakes: active contour models [J]. International Journal of Computer Vision, 1987, 1(4): 321-331.

[6]Klingensmith J D, Nair A, Kuban B D, et al.Segmentation of three-dimensional intravascular ultrasound images using spectral analysis and a dual active surface model [C]//Proceedings of 2004 IEEE Ultrasonics Symposium, 2004, 3: 1765-1768.

[7]Brusseau E, de Korte C L, Mastik F, et al. Fully automatic luminal contour segmentation in intracoronary ultrasound imaging-a statistical approach [J]. IEEE Transactions on Medical Imaging,2004, 23: 554-566.

[8]Filho ES, Yoshizawa M, Tanaka A, et al. Detection of luminal contour using fuzzy clustering and mathematical morphology in intravascular ultrasound images [C]//Proceedings of 27th IEEE-EMBS Annual International Conference on Engineering in Medicine and Biology, 2005: 3471-3474.

[9]Giannoglou G D, Chatzizisis YS, Koutkias V, et al. A novel active contour model for fully automated segmentation of intravascular ultrasound images: in vivo validation in human coronary arteries [J].Computers in Biology and Medicine, 2007, 37(9):1292-1302.

[10]Mendizabal-Ruiz G, Rivera M, Kakadiaris IA. A probabilistic segmentation method for the identification of luminal borders in intravascular ultrasound images [C]//Proceedings of IEEE Conference on Computer Vision and Pattern Recognition, 2008: 1-8.

[11]Papadogiorgaki M, Mezaris V, Chatzizisis Y S, et al.Image analysis techniques for automated IVUS contour detection [J]. Ultrasound in Medicine and Biology, 2008, 34(9): 1482-1498.

[12]Bucher S, Carlier S, Slabaugh G, et al. Shape-driven segmentation of the arterial wall in intravascular ultrasound images [J]. IEEE Transactions on Information Technology in Biomedicine, 2008, 12(5):335-347.

[13]Bovenkamp E G P, Dijkstra J, Bosch J G, et al.User-agent cooperation in multiagent IVUS image segmentation [J]. IEEE Transactions on Medical Imaging, 2009, 28(1): 94-105.

[14]Zhang Q, Wang Y, Wang W, et al. Automatic segmentation of calcifications in intravascular ultrasound images using snakes and the contourlet transform [J]. Ultrasound in Medicine and Biology,2010, 36(1): 111-129.

[15]Cardinal M H R, Meunier J, Soulez G, et al.Intravascular ultrasound image segmentation: a three-dimensional fast-marching method based on gray level distributions [J]. IEEE Transactions on Medical Imaging, 2006, 25(5): 590-602.

[16][美]Paul Schoenhagen, et al. 轻松掌握血管内超声[M].刘茜蒨, 刘 健, 陈 芸译. 北京: 人民军医出版社, 2009: 20-30.

[17]范幸义. 计算机图形学[M]. 重庆: 重庆大学出版社, 2008: 110-120.

[18]孙 正. 变形模型技术在冠状动脉造影图像序列后处理中的应用[J]. 工程图学学报, 2007, 28(4):90-95.

[19]Williams D J, Shah M. A fast algorithm for active contours and curvature estimation [J]. CVGIP: Image Understanding, 1992, 55(1): 14-26.

[20]Papadogiorgaki M, Mezaris V, Chatzizisis Y S, et al.Image analysis techniques for automated IVUS contour detection [J]. Ultrasound in Medicine and Biology Journal, 2008, 34(9): 1482-1498.

[21]Gil D, Hernandez A, Rodriguez O, et al. Statistical strategy for anisotropic adventitia modelling in IVUS [J].IEEE Transactions on Medical Imaging, 2006, 25(6):768-778.

[22]张 麒, 汪源源, 王威琪, 等. 活动轮廓模型和Contourlet多分辨率分析分割血管内超声图像[J].光学精密工程, 2008, 16 (11): 2303-2311.