颈动脉超声图像内中膜分割算法研究（可编辑）

颈动脉超声图像内中膜分割算法研究（可编辑） 华中科技大学 硕士学位 颈动脉超声图像内中膜分割算法研究 姓名:周渊 申请学位级别:硕士 专业:计算机应用技术 指导教师:许向阳 2011-01-16 华中科技大学硕士学位论文 摘 要 颈部动脉血管的内中膜厚度是反映动脉粥样硬化的重要指标。临床中广泛采用 超声成像方法获取颈动脉血管图像,之后由医生手工定点或边界描绘测量内中膜厚 度。为解决医生工作量大、测量结果依赖个人经验和主观性强的问题,有必要借助 图像处理技术研究自动化的内中膜厚度测量方法。 颈动脉内中膜厚度测量的关键是准确分割出上下两条边界。在研究颈动脉图像 特点及现有方法优缺点的基础上,提出一种基于距离约束的活动轮廓模型的内中膜 分割算法。该算法由初始轮廓定位和轮廓线演化两大部分组成。首先使用方向梯度 计算边缘图,对其阈值化并分段得到等宽的二值图像段。在每个图像段中利用两条 边界较平行并相隔一定距离的特点,使用 Hough变换得到靠近内中膜的两条直线段。 之后校正错误的直线段并采样得到内中膜的两条初始轮廓线。昀后使用双 Snake 模 型演化两条初始轮廓线,让它们在保持光滑、靠近边缘、并维持距离一致性的约束 下演化为精确分割内中膜的边界轮廓线。 采用了中南医院提供的 50 幅图像验证该算法的有效性。将新算法、传统动态规 划、双动态规划和传统 Snake 方法的结果与金比较,明新分割算法具有昀小 的误差,且在 95%的置信区间内,内中膜厚度的均值的误差在-0.03mm至 0.08mm之 间,昀大值的误差在 -0.1mm至 0.09mm之间,满足临床应用的实际要求。 关键词:超声图像分割,内中膜厚度,哈夫变换,活动轮廓模型I 华中科技大学硕士学位论文 Abstract Common carotid artery CCA intima-media thickness IMT is an important indicator to cardiovascular diseases. Ultrasound imaging is widely applied in clinical environment for visualization of carotid artery, in which the traditional measurement of IMT is based on manual tracing. To solve the problems such as laborious work and large intra- inter- observer variability accompanying manual measurement, it is necessary to research an automatic method for IMT measurement via image processing techniquesThe key problem in automatic IMT measurement is identification of two nearly parallel boundaries. Based on the characteristics of CCA images and previous research methods, this thesis proposes a distance constraint based active contour model for intima-media segmentation, which comprises an initial contour construction and a curve evolution step. It first uses directional derivative to compute the edge map, which is thresholded and divided into several equal width non-overlapping segments. In each of these segments, according to the characteristics that the local boundaries are nearly parallel and satisfy a certain distance range, two line segments are obtained to approximate the boundaries though Hough transform. The initial contours are constructed by sampling the line segments that are corrected based on the neighboring references. The final step uses a dual snake model to evolve the two initial contours via keeping them smooth, close to boundary and a uniform distance between to obtain the final segmentation resultsFifty images from Zhongnan Hospital of Wuhan University were selected to evaluate the proposed method. The comparisons with respect to the ground truth data among the proposed method, traditional dynamic programming, dual dynamic programming and traditional snake model demonstrated that the new method had the least segmentation error. Moreover, clinical evaluation indicated that the mean and values of IMT had errors from -0.03mm to 0.08mm, from -0.1mm to 0.09mm respectively within 95% confidential interval when compared to those derived from ground truth data, which satisfies clinical requirements Keywords:Ultrasound image segmentation, Intima-media thickness, Hough transform, Active contour model II独创性声明 本人声明所呈交的学位论文是我个人在导师指导下进行的研究工作及取得的研 究成果。尽我所知,除文中已经标明引用的内容外,本论文不包含任何其他个人或 集体已经发表或撰写过的研究成果。对本文的研究做出贡献的个人和集体,均已在 文中以明确方式标明。本人完全意识到本声明的法律结果由本人承担。 学位论文作者签名: 日期:年 月 日学位论文版权使用授权书 本学位论文作者完全了解学校有关保留、使用学位论文的规定,即:学校有权 保留并向国家有关部门或机构送交论文的复印件和电子版,允许 论文被查阅和借阅。 本人授权华中科技大学可以将本学位论文的全部或部分内容编入有关数据库进行检 索,可以采用影印、缩印或扫描等复制手段保存和汇编本学位论文。 本论文属于 保密? ,在_____年解密后适用本授权书。 不保密?。 (请在以上方框内打“?”)学位论文作者签名: 指导教师签名: 日期: 年 月 日 日期:年月日 华中科技大学硕士学位论文 1 绪论 1.1 研究背景、目的与意义 心脑血管疾病是人类健康的首要疾病,人类生命的第一杀手,动脉粥状硬化是 心脑血管疾病的元凶。动脉粥状硬化的发病因素与肥胖、高胆固醇、高血压、糖尿 病、吸烟等有关。颈部主动脉血管(Common Carotid Artery, CCA) 是动脉粥状硬化 病变的常见部位。随着心脑血管疾病的发展,血液中的脂质在血管内壁沉积,使血 管内壁的内中膜结构(Intima-Media Complex, IMC)增生变厚,从而使动脉变窄,医 学上称为动脉粥状硬化。当内中膜厚度(Intima-Media Thickness, IMT)增大到一定 程度,会使血流速度减慢,从而阻碍大脑供血流畅,影响大脑的功能。而且,减慢 的血液流动会加速脂质的沉积,加快动脉粥状硬化的过程,形成局部凸出的斑块。 根据斑块的硬度不同,可分为硬斑和软斑。硬斑容易造成血管的破裂,使血流流入 大脑中,临床上称为中风或脑卒中。软斑容易破裂,裂碎的斑块碎片随着血流沉积 在血管的其他部位,继续引起血管的阻塞。当动脉粥状硬化阻塞血管,心肌会由于 供血不足出现坏死,临床上称为心肌梗塞。因此,心脑血管疾病的早期病变表现在 [1-3] 内中膜增厚上,IMT是预测心脑血管疾病的重要指标 。 临床中,主要是采用医学影像技术获取颈动脉图像,然后再在图像中测出 IMT。 颈动脉的成像包括超声、核磁共振(Magnetic Resonance Image, MRI)和 CT(Computed Tomography)等。超声成像因为具有无侵入、成本低、成像快速等优点,在颈动脉 的检查中获得了广泛的应用。内中膜包含于两条平行边界的内部,管腔内膜边界 (Lumen-Intima Interface, LII)和中膜外膜边界(Media-Adventitia Interface, MAI), LII 和 MAI间的距离即为 IMT。在超声图像中测量 IMT可以通过人工在内中膜边界 [4] 上标定 2点或描绘边界,然后由超声设备提供的工具计算出边界间的距离 。但是, 这种人工的方法依赖于操作者所受的训练、个人经历和主观判断,具有很大的观察 者间的测量差异(Inter-observer Variability)。不仅如此,即使同一个操作者,在不 同时期,测量的 IMT 也不尽相同,即存在观察者内的测量差异(Intra-observer [5] Variability) 。另外,人工测量还非常耗时。 计算机自动分割内中膜并提供 IMT 相关参数,能减少人工测量的工作量,减少 观察者内部、之间的测量差异,已经成为近年来的研究热点。本课题的研究目的在 1 华中科技大学硕士学位论文 于将图像处理技术应用于超声图像内中膜的自动分割中,减少人工测量的繁琐工作 和测量差异,提供 IMT的临床参数,为心脑血管疾病的早期诊断提供可靠依据。 本课题得到了科技部国际合作项目“医学图像特征量化分析方法” (2009DFA12290)的资助,由华中科技大学医学图像信息研究中心与武汉大学中南 医院超声心电图室合作完成,其中,武汉大学中南医院提供颈动脉超声图像,华中 科技大学医学图像信息研究中心负责算法的设计,系统的实现。 1.2 国内外研究概况 下面首先介绍图像分割的基础理论,然后介绍并分析活动轮廓模型(Active Contour Model)、动态规划(Dynamic Programming,DP)和随机优化在内中膜分割 [6] 中的应用及优缺点,昀新的一篇综述 对这个邻域做出了更全面的总结。表 1.1总结 了国内外内中膜分割的相关工作,以时间先后排序,详细的介绍和分析将会在下面 的几节里。 表 1.1 内中膜分割算法综述 时间 所用技术 特点 [5] 1997 动态规划 独立 DP检测两条边界 [7] 2000 动态规划 多尺度 DP,先用 DP检测血管壁,再用 DP检测 边界 [8] 活动轮廓模型 先用 Snake 检测 LII,加入向下的力检测 MAI 2002 [9] 活动轮廓模型 先搜索 LII,再用贪心算法的 Snake 检测 MAI 2006 [10] 2007 动态规划 独立 DP检测边界,启发式图论算法寻找极小代 价 [11] 2007 活动轮廓模型 图像规划化,斑点去除,贪心算法的 Snake [12] 2008 动态规划 DDP 同时检测两条边界,距离范围与变化率约束 [13] 2009 随机优化 Nakagami 分布,贝叶斯推理,昀大化后验概率 [14] 2010 活动轮廓模型 各向异性扩散去噪,使用 GVF 场的 Snake 1.2.1 图像分割的基础理论 图像分割是计算机视觉的底层问题,目标在于将图像平面 Ω划分为互不重叠的 有意义的区域Ω |rR 1, 2,, 。图像分割的关键技术主要依赖于图像的两个特点: r (1)同一区域的像素具有一定的相关性; (2)不同区域的像素在边界处有区别。2 华中科技大学硕士学位论文 无论边缘检测、区域增长、分水岭算法、动态规划、随机优化、活动轮廓模型、 模式识别等技术,它们在图像分割中的应用都要依赖于以上两点,而其他的先验知 识,如待分割对象的位置,形状等只是用来指导计算机对以上两点的理解。因为本 文的研究集中在基于活动轮廓模型的内中膜分割算法,所以这里以活动轮廓模型为 例介绍以上两点在图像分割中的研究概况。 活动轮廓模型在医学图像的分割中有着广泛的应用。这一类方法 源于 Kass 等的 [15] Snake 模型 ,在后来得到了广泛的推广。活动轮廓模型是一种局部化的方法,该方 法在指定的初始轮廓线的基础上,让轮廓线维持自身光滑的同时靠近边界。按照曲 [15-19] [20-25] 线的形式可以把活动轮廓模型分为参数活动轮廓 和隐式活动轮廓(水平集) 两类。参数活动轮廓使用连接的离散化的点表示轮廓线,轮廓线的演化过程具体为 离散化的点在每一次迭代中的位置的移动。参数活动轮廓的优点是实现简单,计算 量小,但参数活动轮廓不易扩展到高维,需要频繁的重新参数化,并且不能适应曲 [26-28] 线拓扑的变化。水平集方法 昀早由 Osher 和 Sethian提出,该方法把二维曲线嵌 入到一个三维的函数中,函数的零水平集表示的曲线被用来演化从而分割对象。水 平集的优点是曲线的拓扑可以改变,二维曲线可以容易推广到三维曲面。但水平集 方法的计算量较大,而且需要频繁的重新初始化符号距离函数。 基于图像分割的以上两个特点,活动轮廓模型可以分为两类,基于边缘和基于 区域。基于边缘的活动轮廓模型得到了昀早的应用, Kass 等的 Snake 模型是其中的 [15] 经典模型 。他们用一条参数化的轮廓线表示物体的边界,在给定初始轮廓线的情 况下,Snake 根据内力和外力运动,在待分割对象的边界处停止。内力控制轮廓线的 光滑,外力控制轮廓线靠近边缘。在 Snake 模型中,人对目标边界的先验知识体现 在初始轮廓线的确定上,而参数化的轮廓线保障了提取的物体边界的连续性和亚像 [16] 素精度。Snake 模型有很多的改进模型,包括 Cohen 的气球力版本 ,Williams 和 [29] Shah的贪心算法版本 以及 Xu 和 Prince的 GVF(Gradient Vector Flow)外力场版 [18] 本 。随着水平集的出现,因为其完善的数学理论基础,拓扑可变 和易于扩展到高 维的优点,基于边缘的活动轮廓模型开始倾向于水平集的研究,这一工作的先驱的 [30] 研究包括 Malladi 等的首次应用 、Caselles 的测地线活动轮廓(Geodesic Active [20] [23] Contour, GAC) 以及 Paragios等引入 GVF场的快速 GAC模型 。不管曲线的形 式的变化,参数化的曲线还是水平集,这些活动轮廓模型的共同点是让轮廓线靠近 边缘,即依赖于以上图像分割问题的第二个特点,因此它们可以归类为基于边缘的 3 华中科技大学硕士学位论文 活动轮廓模型。 还有一类活动轮廓,其思想是让轮廓线的区域内部或外部具有相同的属性(如 [31] 灰度),不同区域间具有不同的属性。这一类方法源于 Mumford-Shah模型 ,这个 模型把图像分割定义为昀小划分问题:给定图像 I ,xy,找出整 个图像平面 Ω的划 0 分 Ω |rR 1, 2,, 和近似原始图像的分段光滑的图像I,xy,使I,xy在每个 Ω 内 r r 光滑变化,而且这种划分的边界长度尽可能小。Mumford-Shah模型在应用水平集的 图像分割中有着广泛的应用。Chan 和 Vese在 2001 年提出无边缘活动轮廓,就是用 [21] 水平集实现昀小划分问题只有两类的简单情况 ,而 Vese 和 Chan 随后把这种实现 推广到多类的情况,并把区域内部为常数或光滑变化分别称为分段常数(Piecewise [32] Constant, PC)和分段光滑(Piecewise Smooth, PS)模型 。昀近的一些类似工作包 [25] 括 Li 等在 PC模型中用局部信息取代全局信息 ,Lankton和 Tannenbaum把这种局 [33] 部化的思想推广到其他依赖区域的水平集中 。因为这些方法的 基本思想是让轮廓 线的内部保持一种属性,外部另一种属性,并且让同一区域的像素的相关性越大越 好,不同区域的像素的相关性越小越好,所以它们可以归纳为基于区域的活动轮廓 模型。 还有一类基于区域的活动轮廓模型则是使用统计的方法,这一类方法的原理来 自模式识别。Zhu 和 Yuille 把 Snake、Mumford-Shah 模型和贝叶斯昀大后验概率 (imum A Posteriori, MAP)联系起来,提出一种受统计力驱动的活动轮廓模型 [17] 。 Paragios 和 Deriche 使用一种有监督的方式分割纹理图像,根据输入的训练纹理, 一组由 Gaussian、LOG(Laplace of Gaussian)和 Gabor 滤波器组成的滤波器组对每 个训练纹理计算频率响应,待分割图像由这些滤波器计算的响应和训练纹理的响应 经过贝叶斯理论结合,得到一个新的曲线演化问题,昀后水平集方法又显示了其在 [22] 处理曲线演化问题的优势 。Rousson等提出一种无监督的纹理图像分割的方法,他 们用结构张量表示纹理的方向信息,非线性的扩散平滑同一纹理内部的结构张量, [34] 保持不同纹理间的结构张量的变化,贝叶斯理论又被用来引导水平集的演化 。昀 [24] 近的一篇综述把使用统计方法的活动轮廓模型用同一个框架统一起来 。因为使用 统计的方法很容易地推广到纹理图像的分割上,所以很多纹理图像分割的研究用到 这类方法。 活动轮廓模型的实际应用还可以引入形状的先验知识。比如在超声图像的分割 中,由于探头的位置和增益的调节,目标的边界经常有消失的情况。形状的先验知 4 华中科技大学硕士学位论文 识可以引导活动轮廓模型克服这一困难。Leventon 等对训练形状的水平集函数进行 [35] 主成分分析,在活动轮廓中加入形状力让曲线接近训练形状 。Chen等在变分中引 入位置、大小和旋转方向,水平集的演化同时受到这些因素的影响,其方法不仅可 [36] 以用来分割,同时可以用来配准 。 1.2.2 活动轮廓模型在内中膜分割中的研究 [14] 国内把活动轮廓模型应用于内中膜分割的例子有李国宽等的工作 。他们在使 [37] 用各向异性扩散 对图像进行去噪后,用固定的阈值找到 LII 的初始轮廓线,然后 使用 Snake 模型演化轮廓线来检测 LII 边界。 MAI 的初始轮廓线由 LII 的边界轮廓下 [18] 移一段距离得到,MAI 轮廓线的演化过程还使用了 GVF 场 ,以防止初始轮廓线离 实际边界距离较远而不能收敛到正确的边界。 [8] 国外的这方面的应用有 Cheng 等的工作 。他们的方法需要人工在内膜上方定 义 2个点,然后根据梯度和光滑约束连接这 2个点得到靠近 LII 的初始轮廓线。MAI 的初始轮廓线由 LII 的演化结果下移固定距离得到。为了应付厚度的变化,他们在 MAI 轮廓线的演化过程中加了一个向下的力,使 MAI 的初始轮廓线也能收敛到具有 较大 IMT的 MAI 边界上。 Ceccarelli等根据从图像上方到下方搜索边界像素得到 LII 轮廓线,在其算法中, [9] 边界像素不仅应有较大的梯度,而且应靠近左边确定的边界像素 。MAI 的初始轮 廓线也是用类似的方法在 LII 的下方区域内搜索得到。一种贪心策略的 Snake 模型被 用来精确地分割 MAI。 Loizou 等使用形态学方法确定 MAI 的初始轮廓线,然后将 MAI 的检测结果上 [11] [29] 移固定距离作为 LII 的初始轮廓线 。他们也是使用贪心策略的 Snake 来精确地分 割边界。与前面的研究不同的是他们使用图像规范化和斑点去除来对图像进行预处 理,并且发现这种预处理能取得更好的效果。 这些基于活动轮廓模型的方法的特点是先检测一条边界,然后将检测的边界轮 廓线上移或下移一段距离作为另一条边界的初始轮廓线。这样的策略有 3个缺点: (1)第一条边界的检测对整个分割算法至关重要,一旦出错,根据第一条轮廓 线确定的另一条边界的初始轮廓线也会与实际边界偏差很远,进而可能使 Snake 的 迭代收敛到不理想的局部较小; (2)第二条初始轮廓线根据第一条检测的轮廓线平移一段固定距离得到,而 5 华中科技大学硕士学位论文 IMT 即使在健康人中也不是一个固定值,并且还会随年龄增大,因而有一个较大的 变化范围,估计的固定距离无法保证这样的平移使第二条初始轮廓线靠近待分割的 边界,可能使随后的轮廓线的演化无法收敛到正确的边界上; (3)这两条边界的检测是独立进行的,每条轮廓线的演化过程没有考虑对方的 位置,因而当 LII 出现边界消失的情况时无法根据 MAI 轮廓线 的位置调整 LII 轮廓 线。 1.2.3 动态规划在内中膜分割中的研究 DP技术的关键在于定义一个代价函数,该代价函数包括边缘强度、边界光滑等 [38] 约束,使代价函数达到极小能提取从左到右的边界 。国内这方面的研究有郭羿等 [10] 做的工作 ,他们根据像素的领域的一阶统计量对图像进行去噪,然后使用图论中 的启发式算法找到使代价函数极小的边界路径。 [5] 国外这方面的工作昀早有 Wendelhag等的探索 。他们把 DP应用到 LII 和 MAI 的检测上,发现 DP 找到的边界测量得到的 IMT 能减少观察者内部、之间的测量差 异,提高可重复性。 [7] Liang 等在以上的工作上做了改进 。因为内中膜的上方为黑色区域,下方为白 色区域,根据尺度空间的理论,当图像的尺度变粗时,图像的梯度 会减小、极值点 减少,细节消失,那么内中膜图像的粗尺度将只会显示血管壁结构,精细的两层 IMC 会消失。于是,他们在粗尺度的图像上使用 DP检测血管壁。然后根据粗略定位的血 管壁,在其上方使用 DP 搜索 LII,在其下方使用 DP 搜索 MAI,从而达到自动分割 的效果。 针对独立使用 DP 方法找出两条边界可能交叉的缺点,Cheng 等提出一种双 DP [12] (Dual Dynamic Programming, DDP)的方法 。该方法在代价函数中同时引入两条 边界,并且加入两条边界的距离的局部变化的约束,以及距离满足一定范围的约束, 使两条边界在保持光滑和靠近边缘的同时,还要满足指定的范围以及较小的距离变 化率。 DP的优点有两条: (1)搜索到的边界是使代价函数极小的全局昀优解,这种全局性保证了算法的 结果一定是对应代价函数的昀理想的边界,而 不会像活动轮廓 模型或随机方法收敛 到局部极小;6 华中科技大学硕士学位论文 (2)在代价函数中容易引入各种约束,如 DDP 中的两条边界的距离变化率和 范围的约束,因为 DP以迭代的方式建立每一步的昀小代价,这种待分割边界的形状 信息很容易被引入到 DP的代价函数中。 但是 DP需要小心地确定参数,在 Liang等和 Cheng等的工作中,他们使用训练 [7, 12] 的方法找到 DP的昀优参数 。然而这种昀优参数因为训练准则的许多可能的局部 极小而不易寻找,而且在实验中,本文发现,即使搜索到一组昀优参数,因为超声 图像的低信噪比(Signal Noise Ratio, SNR),同样的参数不一定能应用于所有实验 图像,往往一组参数对这一类图像很适用但对另一类图像效果不好。 1.2.4 随机优化在内中膜分割中的研究 随机方法在图像处理中有着悠久的历史。这种方法使用贝叶斯计算后验概率, 使之昀大得到去噪或分割的结果。 Geman等使用贝叶斯对图像去噪,他们提出 Gibbs [39] 采样来对概率密度函数抽样 。Storvik 把活动轮廓模型的能量使用贝叶斯来建模, [40] 为了找到昀大的后验概率,他使用模拟退火来搜索全局昀优解 。 [13] 随机方法在内中膜分割中的应用主要有 Destrempes 等的工作 。他们使用 Nakagami 分布来估计超声的射频(Radio Frequency, RF)信号。具体而言,他们使 用 1个成分的 Nakagami 分布来近似管腔或外膜的直方图,而内中膜因为包含不同的 组织,他们使用 3 个成分的 Nakagami 分布来近似。他们使用 EM(Expectation imization)估计 Nakagmai 分布的参数,并且在贝叶斯的先验概率中加入厚度均 匀变化的约束。为了找到后验概率昀大的边界,他们使用 ES(Exploration Selection) 来搜索全局昀优解。 随机方法的优点是容易在贝叶斯中加入各种约束,概率上会收敛 到全局昀优解, 缺点是搜索的时间可能会很长。 1.3 本文研究内容 本文目的在探索一种自动分割内中膜的方法。该方法只需要人工指定包含内中 膜的矩形感兴趣区域(Region of Interest,ROI),然后算法能自动在其中检测到 LII 和 MAI。本文还需要比较新算法与现有 DP 和 Snake 算法的性能的差异,并且评估 新算法测量的 IMT的临床参数的准确性。 本文下面的内容是这样安排的:第二章介绍算法的总体设计,包括分割问题的 7 华中科技大学硕士学位论文 描述,算法的设计思想和总体框架,数据来源和评估标准;第三章介绍本文内中膜 分割算法的关键原理与技术,包括边缘图计算,Hough 变换,初始轮廓线建立和双 Snake 模型;第四章给出实验结果,包括定性的分割结果,与传统 DP(Traditional Dynamic Programming, TDP)、DDP和传统 Snake(Traditional Snake,TS)的比较 结果,临床参数的评估,并分析了算法的准确性、参数敏感性和效率;第五章对全 文总结并展望未来的工作。8 华中科技大学硕士学位论文 2 算法总体设计 鉴于不同分割方法的优缺点,本文提出一种鲁棒性较好的初始轮廓线建立方法, 然后使用活动轮廓模型来分割内中膜。本章首先给出内中膜分割的问题描述,即检 测两条近似平行的边界,然后基于该问题提出算法的设计思想和总体框架,昀后介 绍超声图像数据的来源和分割算法的评估标准。 2.1 问题描述 在临床上,IMT 的测量通常在远侧壁进行。由于超声波的向后散射,近侧壁的 内中膜被声影遮蔽,成像不清晰,在远侧壁进行人工测量具有更好的精确性和可重 复性,因此本文的分割算法集中在远侧壁内中膜的分割问题上。 图 2.1显示超声图像中远侧壁内中膜结构的示意图。图中上方的大图是原始的血 管图像,中间黑色的区域是血管管腔。图中下方的小图是原始图像中白色小框中的 远侧壁内中膜放大后的结果。可以看到,小图中上方的黑色区域是管腔(Lumen), 血液在其中流动;管腔下方白色细长带状区域是内膜(Intima),生理上由内皮细胞 组成,直接与血液接触;内膜下方黑色细长带状区域是中膜(Media),主要由肌细 胞组成;昀下方的白色区域是外膜( Adventitia),由结缔组织和弹性纤维组成。IMC 的边界由管腔内膜边界(Lumen-Intima Interface, LII)和中膜外膜边界 (Media-Adventitia Interface, MAI)定义。其中还存在内膜中膜边界(Intima-Media Interface, IMI)。 从图 2.1中可以看出,内中膜分割可以简化为检测两条边界(LII 和 MAI)的问 题。这两条边界有如下特点: (1)距离满足一定范围,健康成年人的 IMT 在范围为 0.5mm-1mm,基于本文 实验中的像素大小 0.066mm,对应于图像中的 7-15 个像素; (2)两条边界近似平行,健康内中膜的厚度在一段距离内是比较均匀的,这里 不考虑病变产生斑块的情况。 虽然内中膜边界的特点可以简单归纳为以上 2 点,但其分割并不容易,这主要 来自超声图像中的噪声和边界消失等伪像。图 2.2显示内中膜边界被超声伪像干扰的 结果。(a)中管腔存在大量的噪声,(b)中的 LII 边界的一部分消失不见,这些超 声图像中常见的伪像为内中膜的分割添加了额外的困难。9 华中科技大学硕士学位论文 图 2.1 超声图像中的内中膜结构的示意图 图 2.2 超声血管内中膜图像中的伪像10 华中科技大学硕士学位论文 2.2 算法设计思想和框架 基于内中膜的边界的以上特点和分割的难点,本文需要设计一种算法,可以考 虑局部范围的信息以排除噪声的影响,从附近的边缘信息估计消失的边界,并且同 时检测两条边界以让其相互参考。鉴于这一思想,本文使用 Hough 变换在局部的图 像分段中同时检测两条平行直线段,以及同时引入两条边界和其距离变化率的约束 项的 Snake 模型来分割内中膜。 [41-46] Hough 变换是图像处理中检测解析形状的一种常用方法 ,该技术根据形状 的参数方程,把二值图像转换到参数空间的累积矩阵,找到矩阵中昀大值的位置作 为原始二值图像中的形状的参数。虽然 Hough 变换主要用来检测能用解析式表示的 [41, 43] 形状,如直线、椭圆等 ,但也有把 Hough变换推广到检测任意形状的研究,其 [42] [47, 48] 思想主要是基于匹配 。在血管分割中,以前也有 Hough变换的应用 。但 他们主要是把 Hough 变换用来分割管腔,而不是内中膜。内中膜的边界是变化的曲 线,难以用参数来描述。而且内中膜的厚度不是固定值,在一定范围变动。为了克 服内中膜的边界曲线和厚度的变化,本文把 ROI 在水平方向分成若干等宽的图像段, 在每个图像段中使用 Hough 变换检测 2 条近似平行的直线段,然后校正错误的直线 段,昀后将这些处理后的直线段采样得到 LII 和 MAI 的初始轮廓线。 本文使用的 Snake 模型将 LII 和 MAI 的两条轮廓线结合在一起,在泛函中引入 这两条轮廓线的光滑度和靠近边缘项,同时还引入它们之间的距离的变化率项。本 文把这种 Snake 模型称为双 Snake(Dual Snake, DS)模型。使用 Snake 的优点在于 LII 和 MAI 在每个横坐标位置只对应唯一的纵坐标,简化的只更新纵坐标的 Snake 可以避免参数活动轮廓的常见缺点??重新参数化。虽然水平集方法在近年来有更 多的研究,但水平集的优点比如易于扩展到三维领域和拓扑可变在内中膜分割中不 需要,相比之下 Snake 模型简单而且快速。同时演化两条轮廓线的优点是它们可以 相互参考,当一边出现边界消失的情况时,这一边的轮廓线可以根据附近的距离和 另一条轮廓线的位置估计消失的边界。类似的同时演化两曲线或曲面的活动轮廓模 [49-53] 型也在其他文献中出现过 ,但他们的方法主要是从全局上对这两条曲线或曲面 的距离或面积施加约束,而不是像本文提出的方法参考附近的距 离。另外,还没有 类似的同时演化两条曲线的方法被用在内中膜的分割上。 图 2.3 显示本文内中膜分割算法的框架。图中显示 ROI 图像首先经过计算得到 11 华中科技大学硕士学位论文 反映内中膜边界的梯度图,后者经过阈值化和分段处理,得到二值图像段(Binary Image Segment, BIS)。对每个 BIS使用 Hough变换,施加待检测直线的角度和两条 直线距离的约束可以检测两条近似平行并满足一定距离范围的直线段,主直线 (Primary Line, PL)和次直线(Secondary Line, SL),来近似局部 LII 和 MAI。本 文算法的初始轮廓线由所有这些 PL 和 SL 来建立,其步骤主要包括校正错误的 SL 和采样过程。昀后,本文使用 DS 模型来演化初始轮廓线,让其在局部靠近边缘同时 维持两条曲线平行的约束下精确地分割内中膜。 图 2.3 本文算法的总体框架 该算法的特点在于:Hough 变换考虑 BIS 中的所有边缘点,因而不易受噪声的 影响,Hough 的全局性也弥补了超声图像中常有的空洞和边界消失,Hough 变换中 的距离约束和角度约束允许检测的两条边界的距离在一定范围波动,而不是像以往 的方法设置固定值来估计距离,DS 中的均匀能量能根据 MAI 边界和附近的距离估 计消失的 LII 边界,此外,因为 DS的定义是在一个连续的空间中,该算法能达到亚 像素的精度。12 华中科技大学硕士学位论文 2.3 数据来源及评估标准 本文实验的图像来自于武汉大学中南医院超声心电图室。所有的图像在 Philips iE33超声心电图系统上采集得到,使用频率为 11M赫兹的线阵探头。本文中的图像 来自于 25个病人,每个病人选了 1-3个 ROI,这些病人在 2008年到中南医院接受过 超声检查,年龄从 17岁到 86岁。其中,22个病人为健康人,超声检查 IMT1mm, 没有发现斑块。为了验证本文算法对大 IMT 的有效性,本文也选取了超声检查发现 斑块存在的病人,在这种情况下,ROI 的选取避免了包含斑块,只 考虑了无斑块的 平滑边界的内中膜区域。 本文把自动分割的轮廓线和人工描绘的轮廓线相比较来验证自动分割算法的准 确性。为了使边界描述更准确,降低个人因素影响,由 3 个专家人工对每幅 ROI 图 像描绘了 4 次边界,本文对这 4 次的轮廓线求出平均轮廓线得到评估用的金标准 (Ground Truth, GT)。具体方法是,对于 4条人工描绘的轮廓线,通过插值得到对 应横坐标的 4个纵坐标,求其平均值可以建立 GT的轮廓线。本文使用均方差(Mean Squared Error, MSE)来评估两条轮廓线的误差,设 GT 的轮廓线的某一点的纵坐标 为y ,自动分割的轮廓线的对应点的纵坐标为y ,则 MSE定义如下。 GT AM 2 MSE? E y y 2.1 GT AM 类似的评估标准也在文献[7, 8]中用于训练和误差比较。图 2.4显示 MSE的评估方法, 其中y 和y 根据插值得到。 GT AM 图 2.4 MSE 的评估方法13 华中科技大学硕士学位论文 为了评估本文算法自动分割的轮廓线得到的 IMT 临床参数的准确性,本文把两 条轮廓线插值得到每个横坐标对应的厚度,统计整个 ROI中的 IMT的均值和昀大值, 将这些参数与 GT得到的数据相比。因为内中膜的边界并不是严格的水平线,有的边 界甚至出现较大的倾斜角度。为了消除倾斜的影响,本文的厚度的计算还考虑到斜 率k的影响。设两条轮廓线在某一横坐标的点对的垂直距离为d ,通过计算所有的点 对的平均点,可以使用曲线拟合得到光滑曲线并估计其在点对处的局部斜率k ,然后 2 由三角形学可得考虑斜率的厚度d 1 +k 。根据每一列像素的单点厚度,可以统计 出整个 ROI 中的厚度的均值 IMT 和昀大值 IMT 。图 2.5显示厚度的计算示例, mean 其中, Γ 和 Γ 为上下两条轮廓线,d 表示某一点的垂直距离,k为拟合的平均曲线 1 2 在该点的斜率。图 2.5 厚度的计算示例 2.4 本章小结 本章介绍了内中膜分割问题的描述(平行边界检测),算法的设计思想和总体 框架,血管超声图像和 GT数据的来源及算法性能的评估标准。14 华中科技大学硕士学位论文 3 关键理论与技术 第二章介绍了本文的内中膜分割算法的框架,本章将详细介绍该框架中各部分 的理论与实现技术。 3.1 边缘图计算 边缘图用来反映待分割物体的边界,传统的边缘图使用梯度幅值来表示。设 Ix,y表示原始 ROI 图像,坐标系统的原点在图像的左上角,则梯度幅值的定义如 下。f,xy ?G *I,xy 3.1 σ 2 2 x +y1 2 2 σ 其中,G x,y e 是 2D高斯函数, ?为梯度算子,*表示卷积运算符,σ 2 2 π σ 表示幅值。式3.1会反映图像中各方向的梯度变化,从而突出目标边界,但噪声的 边缘也会较明显。 在内中膜分割中,由于 LII 和 MAI 均成水平状,因而边缘图可以只计算垂直方 向的梯度。而且,LII 和 MAI 的检测会受到 IMI 的干扰。为了消除这种影响,根据 LII 和 MAI 的下方灰度比上方灰度大的特点,以往的方法使用特定的边缘检测算子 [8, 12] 来突出这种边缘 ,而本文使用单向边缘算子处理这一问题。式3.2为本文定义的 边缘图的计算方法。 ?GI *x,y σ fx ,y 0, 3.2?y其中, 表示 y 方向的导数, .,.运算在两个数 组成的集合中取昀大值,具有 ?y 将导数的负值设为 0 的效果。式3.2只保留 y 方向正方向的梯度,而消除了负方向 的梯度。 图 3.1 显示了本文定义的单向边缘图相对于传统基于梯度幅值的边缘图在内中 膜分割中的优势。其中(a)为原始图像,(b)显示(a)的梯度幅值,由式3.1计 算,(c)显示(a)根据式3.2计算得到的边缘图。可以看出,因为式3.1反映了各 个方向的边缘,管腔、外膜中的边缘在(b)中都被保留了下来,而且 LII 和 MAI 间的 IMI 与边界相隔太近,会对分割造成较大的干扰。相比之下,因为式3.2只保 15 华中科技大学硕士学位论文 留在 y 方向下方灰度比上方灰度大的边缘,所以 IMI 产生的不利边缘及管腔、外膜 中的反向边缘在(c)中都被消除了。图 3.1 梯度幅值与单向边缘图的比较 3.2 Hough 变换 3.2.1 阈值化与分段 本文假设内中膜的边界可用分段直线来近似,这是因为医生在测量病人的 IMT 时,可以在超声成像中选择一段边界平滑的内中膜来分析。这里本文不考虑斑块的 因素,因为斑块的边界可能不,分段直线将不能应用于其中。因为直线检测的 常见方法是 Hough变换,而 Hough变换的输入为二值图像,故本文对边缘图阈值化 后得到二值图像,将其沿水平方向划分成S 个等宽的 BIS B |s 1,2,,S ,在每个 s BIS 中使用 Hough 变换检测两条边界。因为阈值化的对象是边缘图,而且后面的 Hough 变换也在一定程度上弥补了阈值化中漏掉的边缘点,所以本文使用固定的全 [54] [55] 局阈值T ,其他的阈值化方法如迟滞阈值 、动态阈值 也可以用来改善二值边缘 图像。在实验中,本文对边缘图像规范化到 0-1,设T 0.3进行阈值化,然后设 BIS 16 华中科技大学硕士学位论文 的个数为S N /l ,其中N 为图像的宽度,以像素为单位,l 为 BIS可能的昀 大宽度,这里设l 90,后面的实验已验证该参数可以在一个较大的范围内取值。 3.2.2 Hough直线检测方法 Hough 变换根据输入的二值边缘图,把每个边缘点映射到待检测的形状的参数 空间的累积矩阵中,把每个经过这一边缘点的形状的参数在累积矩阵中的位置加 1, 昀后统计累积矩阵中的昀大元素,其位置对应 待检测的形状的参数。具体针对直线 而言,Hough变换使用式3.3表示直线。cos ρ xyθθ +sin 3.3 其中, ρ 表示直线与原点的距离, θ 表示直线的方向。图 3.2显示这种直线的表示方 式。图 3.2 直线的表示方式 根据以上描述,Hough变换的累积矩阵定义如式3.4。 K Hv , ρθρ x,y, ,θ 3.4 ? kk k 1 其中,0 ? ρ L,0 ? θ π ,L是图像平面 Ω的对角线的长度,x 和y 为第k个边 k k 缘点的坐标,k 1,2,,K ,vx ,y , ρ, θ 是投票函数,由式3.5定义。 k k 1, if ρ+ x cosθθ y sinkk vx,y,ρθ, 3.5 kk 0, otherwise 17 华中科技大学硕士学位论文 从式2.3和2.4可以看出,Hough 变换的累积矩阵中的元素H ρ, θ 为由 ρ, θ 表示 的直线经过的边缘点的个数。因为严格意义上 ρ x cos θ + y sin θ 完全成立的情况 k k 很少,所以认为在一定误差范围内成立就可以将对应的 ρ, θ 处的元素加 1。在累积 矩阵的建立中,经典的方法是将参数空间划分为M ×M 的格子,初始化H ρ, θ 0, mL m πρ ? |m 0,1,,M ?1 , θ ? |m 0,1,,M ?1 ,对于每一个 θ ,根据M M[41] ρ x cos θ + y sin θ 计算出对应的 ρ ,把H ρ, θ 中昀接近的格子的元素加 1 。因 k k 为这一步骤重复K次,所以该算法的计算复杂度为OKM 。后来有许多算法被提出 [42] 来改进 Hough 变换的效率,比如使用梯度方向信息来缩小 θ 的取值范围 ,随机抽 [43] [44] 取 2个点来估计参数 ,随机抽取一小部分边缘点来代表整体 等。 3.2.3 Hough变换检测 LII 和 MAI 本文使用 Hough 变换在每个 BIS 中检测近似平行的 2 条直线段作为局部 LII 和 MAI 的近似。根据 3.2.2节的描述,本文从累积矩阵中搜索昀大值来估计一条边界, 则 LII 或 MAI 可以根据式3.6检测。, ρθρ argH ,θ 3.6 11 , ρθ 其中, ρ , θ 为昀大元素的位置,对应 BIS中昀明显的边界。 本文把 ρ , θ 对应的 1 1 1 1 直线段称为主直线 PL,并假设 PL一定对应 LII 或 MAI。 另一条边界可以在累积矩阵中搜索另一个峰值 ρ , θ 来检测,本文把 ρ , θ 对 2 2 2 2 应的直线称为次直线 SL。这里本文引入内中膜的形状信息,因为无病变的内中膜的 两条边界近似平行,故可以认为这两条直