第22章：颈动脉超声图像分割（Segmentation of Carotid Ultrasound Images）

1 Introduction

世卫组织数据显示，中高收入国家中脑血管疾病是主要死因之一。预防是降低全球死亡人数的关键，影像技术在其中扮演重要角色。然而相对心肌梗死而言，卒中医学进展并不显著。心肌梗死通常由冠状动脉内血栓斑块破裂引起，而卒中的成因则更复杂，可由大血管栓塞或血栓性闭塞（缺血性卒中）、血管破裂（原发性脑出血）或小血管病变（腔隙性缺血或原发性脑出血）造成。在症状前人群中精确区分这些病理过程不仅可揭示卒中的病理生理学机制，也为预防性治疗打开窗口。大血管疾病与动脉壁退化和脂质等血液成分在血管壁内沉积的动脉粥样硬化过程相关，几乎涉及所有血管床。

动脉粥样硬化的诊断是预防心脑血管事件（如心肌梗死和卒中）最重要的医学检查之一。颈动脉内中膜厚度（IMT）是衡量大血管病变的指标，与脑认知功能相关；IMT 提供动脉粥样硬化的指数，用于心血管风险评估。IMT 可通过超声成像测量，超声使用安全的非电离辐射。超声图像的特点是信噪比低且存在多种成像伪影，使其解读和分割非常困难。图像受三类噪声影响：（1）斑点（speckle），一种赋予图像颗粒感的乘性噪声；（2）回声混响，引入虚假边界；（3）声影或回声衰减，可能掩盖边界或边界段。另一方面，超声 IMT 测量具有高度可重复性、无创性，可用于监测病理和药物疗效。如图 22.1 所示，B 型超声显示血管壁为与解剖层相关的规则模式：壁层表现为两条几乎平行的回声线（内膜层和外膜层），中间被低回声的中膜层分隔，形成所谓"双线"模式。这一模式的内膜-中膜（IM）部分代表从动脉腔-内膜（LI）边缘开始到中膜与外膜之间（MA）边界结束的组织区域。IMT 测量通常在颈总动脉（CCA）远壁（FW）进行，因为该处"双线"模式更易观察；但对于某些应用（如三维动脉表面重建），近壁（NW）分割也可能需要。

2 Segmentation of Carotid Ultrasound Images. Brief Survey

颈动脉图像分割的挑战在于获得准确的 IMT 测量，这要求同时检测 MA 和 LI 边界。过去二十年里发表了多种颈动脉壁分割和 IMT 测量方法，包括：边缘检测 [13–19]；灰度密度分析 [20]；动态规划 [21–24]；snakes [25–36]；离散动态轮廓与多尺度分析 [37]；Hough 变换 [36, 38, 39]；模板匹配 [40, 41]；Nakagami 建模 [42]；特征提取、拟合与分类 [43–45]；特征提取、三次样条、几何轮廓与动态规划 [46, 47]；分水岭变换 [48]。还有使用 inter-greedy 技术融合不同分割方法的研究 [49–51]，但需要手工描迹作为真值来估计每次迭代的误差。多数方法的综述可在 [52] 中找到。

2.1 Edge Detection

1998 年 Pignoli 和 Longo 发表了第一个计算机辅助的二维 B 型颈动脉图像分割方法 [13]，该方法通过从管腔向远壁方向移动检测血管结构边缘，LI 和 MA 界面对应的边缘用于测量 IMT。Touboul 等 [14] 采用类似方法，证明了半自动计算机辅助分割可降低观察者内和观察者间变异性。

Selzer 等 [15, 16] 用基于边缘检测的方法测量 IMT 并评估其随时间的变异性。操作者用鼠标在第一帧图像的 LI 和 MA 界面上标识若干点，系统通过在这些点附近的平滑曲线邻域搜索边缘来检测界面。后续帧重复同一过程，以前一帧的 LI 和 MA 界面作为引导，该过程可能需要对检测错误进行手工校正。

Liguori 等 [18] 提出另一种基于边缘检测的方法，采用统计阈值化降低噪声后再计算强度梯度，但需要手工选择感兴趣区域（ROI）。

Stein 等 [17] 使用边缘检测方法证明计算机辅助的 IMT 测量比手工测量更快、可重复性更高、更准确。ROI 必须包含一段管腔及对应的内膜、中膜、外膜层；操作者定义 ROI 长度并在 ROI 左限的管腔区域内手工标识一个点，然后自动检测 LI 和 MA 界面，对错误检测的边界进行手工校正。

Faita 等 [19] 提出一种最准确的边缘检测方法之一，通过使用一阶绝对矩边缘算子（FOAM）和模式识别方法改善噪声存在时的边缘检测。该方法的优点是实时，但有两个重要局限：首先，难以处理弯曲或非水平的血管；其次，需要手工选择 ROI。

2.2 Gray Level Density Analysis

Gariepy 等 [20] 使用基于灰度强度和组织识别的计算机化方法，证明超声和组织学 IMT 测量之间存在强相关性，并证明高血压与大动脉 IMT 异常增加相关。为避免双线动脉模式受斑块干扰，斑块患者被排除在研究之外。ROI 手工选择，由包含待分析远壁段的矩形组成。该方法检测两壁的 LI 界面和远壁的 MA 界面，用以计算平均管腔直径和远壁平均 IMT。该计算机化方法未与任何手工描迹评估，主要用于降低观察者内变异性。

2.3 Dynamic Programming

Gustavsson 等 [21] 和 Wendelhag 等 [22] 提出相关的动态规划过程，用于 B 型图像中颈动脉边界的自动检测。其方法将局部强度、强度梯度和边界连续性三种测量线性组合为代价函数，代价函数低值对应目标边界像素，高值对应其他像素。代价函数三项（强度、梯度和连续性）的权重在训练阶段以手工描迹为真值进行估计。过程完全自动，但允许操作者交互以修正错误检测。为降低计算负担，Liang 等 [23] 引入 [22] 中动态规划方案的多尺度版本，在粗尺度计算动脉的全局位置，再逐尺度向下细化。这些方法的主要缺点是需要训练阶段，且在更换扫描仪或其设置时系统可能需要重新训练。

最近，[24] 中引入了一种不同的基于动态规划的方法。瞬时变异系数（ICOV）用于改善超声图像中的边缘检测，因为超声图像的噪声具有乘性特征。基于 ICOV 的边缘检测、强度梯度、稳健统计和关于动脉边界空间分布的先验知识被用于计算 LI 和 MA 界面的模糊评分图，然后这些评分图被传递给动态规划过程，计算具有最大累积得分的轮廓。该方法能够对健康和动脉粥样硬化动脉的近壁和远壁 LI 和 MA 界面进行检测，适用于多种类型和大小的斑块，但对远壁的检测结果显著更好，其自动检测的准确性与手工检测相当。主要缺点是需要手工选择 ROI 和训练阶段。

2.4 Snakes

传统 snakes 常因噪声影响被错误边界吸引。为克服这一局限并防止 snake 在 LI 和 MA 界面之间被捕获，Cheng 等 [25, 27] 使用更稳健的 snake 外能公式，但需要手工初始化 snake。在另一项研究 [26] 中，同一作者表明使用基于 snake 方法估计的轮廓代替手工描迹可显著降低观察者内和观察者间变异性。

Loizou 等 [28–30] 表明如果在 snake 之前进行强度归一化和去斑，可改善基于 snake 方法的结果。强度归一化将管腔区域的强度中位数设置为 0 到 5 之间，将外膜层的强度中位数设置为 180 到 190 之间。最佳去斑方法包括对基于 7×7 像素窗口统计量的线性缩放滤波器迭代五次。强度归一化和 snake 初始化都需要用户交互。作者呈现了与手工描迹对比的大量统计验证。

Delsanto 等 [32, 34] 提出一种完全独立于用户的方案，结合局部统计与 snake。局部统计首先用于定位颈动脉管腔，然后基于 snake 检测 LI 和 MA 界面。其假设是具有低强度均值和标准差的邻域像素通常属于管腔，因此通过聚类将图像分为两类，将具有这两个性质的像素与其他像素分开。均值和标准差的阈值均通过经验确定，使用在 10×10 像素窗口上计算的表示强度标准差随强度均值变化的二维直方图。一旦定位管腔区域，通过在每列图像上沿管腔下方搜索高于给定阈值的强度峰值来估计外膜层；将纵向连接这些强度最大值的轮廓作为初始化 snake 以细化 MA 界面估计。LI 界面检测采用类似过程，但使用管腔与 MA 界面之间最大的强度峰值来初始化 snake。该方法与手工描迹对比验证，结果良好，但存在若干参数经验选择、计算负担过重不适用于实时应用，以及 10% 病例因噪声导致错误检测等局限。同作者在 [31, 33] 中提出改进方案，处理近壁和病变血管。考虑三个聚类：（i）管腔；（ii）内膜和中膜层；（iii）外膜层。聚类使用模糊 K-means 分类器实现。类（i）和（ii）之间的边界作为 LI 界面的初始猜测，类（ii）和（iii）之间的边界作为 MA 界面的初始猜测。这些边界用于初始化 snake 以细化 LI 和 MA 界面估计。该章报告了 IMT 测量误差的改善以及对斑块分割的稳健性，但对低回声斑块的结果较差。该方法在 8% 的图像中失败，原因包括低信噪比或 ROI 选择失败。同方法更近期的版本 [35] 在用三种不同超声扫描仪采集的 200 张图像集合上测试，证明分割性能不依赖于所用扫描仪，并报告了 IMT 测量误差相对初始版本的改善。

2.5 Discrete Dynamic Contour and Multi-scale

Gutierrez 等 [37] 在多尺度方案中使用离散动态轮廓自动检测 B 型图像中的颈动脉边界。轮廓演化被实现为三种加权力的线性组合函数：（1）内力，与轮廓曲率成正比；（2）外力，等于强度梯度的局部幅值；（3）阻尼力，与轮廓每个顶点的速度成正比。结果显示 IMT 测量准确度较低。

2.6 Hough Transform

Golemati 等 [38, 39] 提出基于 Hough 变换的方案。假设 LI 和 MA 边界在纵向截面中为水平直线，在横截面中为圆形。该方法完全自动、实时，可分割颈动脉的纵向和横截面。然而需要 LI 和 MA 界面为水平直线或圆形是一个严重局限，因为这些边界的外观常不规则。另一缺点是大量参数经验设置。需要对超声图像进行归一化以最小化不同操作者、扫描仪和增益设置引入的变异性。最近，同一作者 [36] 引入对颈动脉横截面边界的基于 snake 的细化，其中 snake 由 Hough 变换方法计算出的圆初始化。结果通过自动轮廓与手工描迹之间的 ROC 分析度量进行验证。

2.7 Template Matching

Rossi 等 [40, 41] 引入一种有趣的参数化模板匹配方法，用于颈总动脉纵向截面中的自动管腔检测。该方法不计算 LI 和 MA 界面，因为设计目的仅为估计管腔位置。但该技术可作为 LI 和 MA 检测完全自动化的第一步，因为它避免了手工定位管腔的需要。为降低计算负担，对图像列进行抽取。所提出的模板匹配方案基于对预期直径范围和动脉壁-管腔复合体典型强度模式的先验知识。在帧序列上执行空间和时间聚类以拒绝错误估计。该方法实时但需要训练阶段。方法的性能通常较高，但当颈静脉存在或管腔区域信噪比较低时可能显著下降。算法通过比较自动估计的管腔中心位置与操作者描迹的参考位置进行验证，若自动识别与参考位置偏差不超过 2 mm 则视为正确。

2.8 Nakagami Modeling

Destrempes 等 [42] 使用基于随机优化的分割方案，假设在包含 CCA 壁的小 ROI 内强度分布服从 Nakagami 模型。管腔和外膜分别与均值最低和最高的分布相关，内膜-中膜复合体被视为混合。虽然报告了良好性能，但该方法对含有病理结构（如斑块）的图像可能用处有限，因为壁层的统计分布假设可能不成立。例如在钙化斑块存在的情况下，外膜层可能不是均值最高的一层。另一缺点是需要手工初始化。

2.9 Feature Extraction, Fitting, and Classification

Molinari 等 [43, 44] 提出由特征提取、拟合和分类组成的集成方法。技术分为两个主要步骤：（1）CCA 的自动定位；（2）远壁 LI 和 MA 界面的检测。第一步从种子点定位开始，种子点定义为每列图像中高于给定阈值的局部强度最大值，将线段拟合到这些种子点；对齐且接近的线段被融合，其他线段被消除，产生外膜层位置的估计。第二步对每列图像的强度剖面应用模糊 K-means 分类器，将像素分配到管腔、内膜-中膜结构或外膜层，分隔这三类的边界被视为 LI 和 MA 界面的估计。还有一个预处理模块使用低通滤波器减少斑点。与手工描迹的大量验证显示 MA 界面检测非常准确，但 LI 界面检测较差。另一缺点是第一步中需要训练阶段。

最近，同一作者 [45] 引入该方法的改进版，其中第二步中 LI 和 MA 界面的检测基于 Faita 等 [19] 之前用过的 FOAM 边缘算子。该方法在来自两个不同机构的 300 张颈动脉图像数据库上与手工描迹对比验证，对远壁 LI 和 MA 界面均产生非常准确的检测。作者报告平均计算时间为每张图像 2.3 秒，失败率为 4%。

2.10 Watershed Transform

Molinari 等 [48] 采用与 [43, 44] 中类似的方法，但在 CCA 定位步骤使用基于分水岭变换的算法。在用 12 像素盘形结构元素腐蚀后，图像被针对原图进行重建。通过阈值化计算黑白图像，白色区域用于执行基于标记的分水岭分割。然后将分割图像中自定义评分函数值最大的区域选为 CCA 位置。CCA 区域包括管腔和动脉壁，但管腔远暗于且大于壁。颈静脉（若存在）具有类似性质但出现在 CCA 上方，因此评分函数与区域中暗像素百分比成正比，与区域平均强度和区域平均高度成反比。该方法与手工描迹对比测试。报告的失败率（9.5%）仅略低于 [43, 44] 中方法。在 165 张图像的数据库上的评估显示相对 [43, 44] 的准确度有改善，但低于 [45]。平均分割时间为每张图像 18 秒。

2.11 Limitations of the Current Solutions

上述方法存在若干重要局限。半自动方法 [13–22, 24–30, 37, 42] 不适合处理大型数据库，结果受用户交互影响。若干自动方法计算负担较重 [31–35, 43, 44, 48]，使其在临床实践中缺乏吸引力。某些研究 [13–16, 20, 21, 25, 27, 37] 未考虑斑块的分割。少数研究 [18, 24, 31, 33, 46, 47] 测量近壁 IMT。

已有一些 IMT 测量的商业系统，包括通用电气公司的 Vivid 7、SonoSite 公司的 SonoCalc IMT 和 Intelligence in Medical Technologies 公司的 M'Ath。这些系统计算指定段的 IMT 测量（如均值和最大值），但也存在局限：它们都需要用户交互；Vivid 7 和 SonoCalc 需要手工选择包含动脉壁段的矩形 ROI；M'Ath 需要用户在管腔内绘制纵向线以定位管腔并定义待分割壁段。它们仅在待检测边界在指定壁段中良好可见时正常工作；至少在 M'Ath 的情况下，该段不能包括斑块。

3 Cubic Splines and Dynamic Programming for Carotid Image Segmentation

所提出的颈总动脉 B 型超声图像近壁和远壁内膜-中膜区域分割方法从边缘检测步骤开始。最终边缘图由与待检测颈动脉壁边界兼容的边缘组成。MA 界面通过 RANSAC 搜索将三次样条最佳拟合到具有与 MA 边界兼容特征的边缘来检测；使用 MA 界面的若干判别特征来减少自动轮廓被其他边缘吸引。通过动态规划、平滑强度阈值化表面和几何 snake 获得 LI 界面的平滑估计。

3.1 Edge Estimation

边缘估计包括确定表示像素集合的边缘图，这些像素的属性（如梯度方向）与待检测的颈动脉壁边界兼容。边缘按四步估计：边缘检测（第 3.1.1 节）；边缘处主导梯度方向的估计（第 3.1.2 节）；最终边缘图的确定（第 3.1.3 节）；谷值边缘图的确定（第 3.1.4 节）。

3.1.1 Edge Detection

[46, 47] 中提出一种用于 B 型图像边缘检测的新型平滑滤波器，可在平滑噪声的同时产生良好定位的边缘并保留重要的弱边界。该滤波器利用瞬时变异系数（ICOV）边缘检测器 [53, 54]，后者非常适合超声图像。像素 (x,y) 处的 ICOV 值由 (22.1) 式计算，使用其 8 邻域内像素的强度。在 (22.1) 中，I(x,y) 表示像素 (x,y) 处的图像强度，∇I(x,y) 是像素 (x,y) 处的强度梯度，∥(u,v)∥ = √(u²+v²) 是向量 (u,v) 的范数。如 [53] 所述，ICOV 是斑点图像的边缘检测器，结合了归一化梯度幅值算子和归一化拉普拉斯算子。在边缘像素处，拉普拉斯项为零，梯度项占主导，ICOV(x,y) → |∇I(x,y)|/I(x,y)。这种梯度幅值的归一化使 ICOV 能够在具有乘性噪声的图像中检测亮区域和暗区域中的边缘。如图 22.2 所示，ICOV 在强度变化最剧烈的边缘像素点处产生局部最大值。

如 [55] 中那样，稳健统计用于决定扩散应该发生的位置以及应该被抑制的位置。像素 (x,y) 和时间 t 处的扩散系数是由 Tukey 函数 [56] 给出的（式 22.2），其中 σₛ = √5σₑ，σₑ 是图像边缘尺度 [55]，由 (22.3) 式计算。MAD 表示中位数绝对偏差，med(r) 是 r 在图像域上的中位数，C = 1.4826 是常数。

所提出的滤波器的优点是即使在重要解剖边界具有低 ICOV 时也能保留它们。它基于总变分理论 [57, 58] 的概念并嵌入曲率信息，如 (22.4) 式所述，其中 c(x,y;t) 是 Tukey 函数（式 22.2），∇I 是强度梯度，I₀ 是初始图像（在 t=0 时），∂Ω 是图像边界，n 是图像边界处的外法线。κ(x,y) 是平均曲率，在每个时间步更新，由 (22.5) 式给出。

使用该滤波器，扩散在 ICOV 高和曲率小的位置都被抑制。这样，由于 LI 和 MA 边界通常曲率较低，即使它们具有低 ICOV 也可以被保留。噪声通常具有高曲率和低 ICOV，因此被强烈平滑。图 22.3 显示了该新滤波器相对 [55] 中先前提出的相关滤波器的改进。图 22.3c 显示 [55] 中滤波器产生的弱边缘侵蚀，图 22.3e 显示新滤波器更好地保留了颈动脉壁边界。

边缘图从新滤波器产生的图像平滑版本计算，使用 ICOV 作为边缘强度度量、非极大值抑制和滞后 [59]。滞后低阈值设为 T₁ = σₑ，高阈值设为 T₂ = 0.4T₁。形态学细化 [60] 用于确保边缘为单像素宽度。

3.1.2 Dominant Gradient Direction

[46, 47] 中通过使用局部主导梯度方向降低梯度方向误差，计算方法如下：设 ∇Iₙ(x,y) 为像素 (x,y) 在迭代 n 处的强度梯度，∇Iₖⁿ⁻¹ 为像素 (x,y) 的 8 邻域中第 k 个像素在迭代 n−1 处的梯度。那么 ∇Iₙ(x,y) 通过对 ∇Iₖⁿ⁻¹（k = 1, 2, …, 9）取平均来计算，仅考虑与中心像素梯度方向夹角小于 45° 的邻域，以避免方向差异较大的相邻轮廓的干扰。

迭代过程的停止准则基于梯度方向的稳定性。设 α 为连续迭代之间每个边缘像素处梯度方向的角度变化；离群点可定义为 α 值未稳定的边缘像素。不再有 α 内点预期的阈值 σα 可估计为 σα⁻ = √5σα [56]，其中 σα = C·MAD(α) + med(α) 是离群点开始出现的阈值的稳健估计 [55]。因此当 σα⁻ < ε 时迭代停止，设 ε = 0.1° 以确保所有内点的良好稳定性。

3.1.3 Final Edge Map

最终边缘图通过使用基于其梯度方向、其到管腔轴的距离以及其到管腔边界的符号距离（SDL）值的判据，仅选择与待检测颈动脉壁边界兼容的边缘来确定。不兼容边缘的消除降低了计算成本和自动轮廓被其他边缘吸引的可能性。

由于在颈动脉壁边界处梯度应指向外部，因此指向动脉内部的边缘应被移除。如果 θ(x,y) 是像素 (x,y) 处强度梯度与到估计的中轴的距离图梯度之间的夹角，则所有满足 θ(x,y) < θ_max 的边缘像素从边缘图中移除，其中 θ_max 是找到颈动脉 MA 界面边缘概率几乎为零以上的阈值。

ROI 中所有到管腔轴距离大于某阈值 d_max 或满足 SDL_min < SDL < SDL_max 的边缘也从边缘图中移除。

3.1.4 Valley Edge Map

MA 界面常与谷形强度剖面相关（图 22.4），称为"双线"模式 [10]。

如 [46, 47] 所述，谷值边缘图是第 3.1.3 节所述最终边缘图的子集。谷值边缘图的确定从沿每个边缘点 e 及其对应强度梯度 ∇I(e) 定义的线在两个方向上搜索一定距离 L 内的第一个局部强度最大值（图 22.4）开始。谷形边缘的强度剖面在位置 a 和 b 处有两个强度峰值，其中之一通常较低。只有具有强较低峰的剖面才应被分类为谷形边缘。如第 3.1.1 节使用滞后，并用 A 表示较低峰的振幅，则高阈值可设为 T_A = C·MAD(A) + med(A)，低阈值设为 0.4T_A。然而实验表明，若当 A > 0.4T_A 时将边缘像素分类为谷形边缘，则谷形边缘检测更佳。

3.2 Estimation of the Media–Adventitia Interface Using RANSAC

[46, 47] 中提出的 MA 边界分割方法基于根据指定代价函数对三次样条 [61] 模型进行最佳拟合的 RANSAC 搜索。RANSAC 算法允许从含有大量离群点的数据集中估计模型参数，其工作原理是反复提取包含确定模型参数所需最少数数据点的随机样本，然后评估模型对其余数据的共识，选择具有最佳共识的模型。当对至少抽到一个好样本具有高置信度时终止过程。

选择三次样条 [61] 模型用于 MA 边界，因为其产生平滑曲线、实现相对容易、行为稳定，且结果表明它能够充分跟踪纵向截面中的 MA 界面。设计了专门的增益函数来评估样条共识。

3.2.1 Sample Generation and the Adventitia Model

在 CCA 的纵向截面中，MA 边缘的样本必须具有不同的横坐标。因此随机抽取一组 n 个不同横坐标，分别用于在管腔轴上方和下方确定 n 条垂直线。好的横坐标是对应垂直线包含 MA 边界好的边缘点（即 MA 边界的一个边缘点）的那些。

为检测近壁 MA，从每个由 n 个不同横坐标的边缘点样本构建三次样条，位于管腔轴上方。由于每个横坐标通常有几个边缘点，算法对每个由 n 个横坐标的样本拟合的所有样条进行评估，然后根据预定义准则选择最佳样条。检测远壁 MA 界面使用类似过程。设 n = 5 在不过度损失对噪声的稳健性的前提下为样条提供一定灵活性。

3.2.2 Model Consensus

拟合样条的共识通过集成对颈动脉边界若干判别特征响应的增益函数来度量。

一个特征与边缘类型相关。由于谷形边缘并非总是存在，因此也考虑阶梯形边缘，对谷形边缘和更接近已检测管腔边界的边缘赋予更大权重。另一特征是 SDL。对 SDL 较大绝对值施加更强惩罚以减少其他解剖边界的影响。然而斑块的可能存在要求该惩罚对 SDL 正值（管腔外）增长更慢。第三个特征是梯度方向，必须与颈动脉 MA 界面模型的法线方向一致。总之增益函数必须集成以下特征：

拟合样条到任何边缘点的距离 dₑ
拟合样条到谷形边缘点的距离 d_ve
拟合样条法线方向与强度梯度之间的夹角 δ
拟合样条到管腔边界的符号距离 SDL

这些特征以反映每个样条点属于颈动脉边界概率的方式集成到 (22.6) 式所示增益函数中。如果样条点 Pₖ 接近谷形边缘或阶梯形边缘，并且其强度梯度方向与样条法线方向相似，且落在与管腔轴的预期距离限制内，则其属于颈动脉边界的概率较高。每个样条点成为好候选的概率应随每个特征的增强而增加。因此最佳样条拟合选择为对以下增益函数得分最高者：

G = (1/(2m)) Σₖ₌₁ᵐ [g₁(Pₖ) + g₂(Pₖ)]·g₃(Pₖ)·g₄(Pₖ) (22.6)

其中 m 是数字样条点数 Pₖ，gⱼ(Pₖ)（1 ≤ j ≤ 4）是表示特征 j 在 Pₖ 处贡献的模糊函数。在 (22.6) 中，g₁(Pₖ) = f(dₑ(Pₖ))，g₂(Pₖ) = f(d_ve(Pₖ))，g₃(Pₖ) = f(δ(Pₖ))，其中 f(z) 是 Tukey 函数（式 22.7），特征 dₑ 和 d_ve 的尺度为 σ = σ_d，特征 δ 的尺度为 σ = σ_δ。

每个模糊函数的尺度 σ 表示相应特征的阈值，超过该阈值则不再预期找到 MA 边界像素。函数 g₄ 给偏好更接近管腔边界的曲线，并对管腔外的曲线给予更大容差，由 (22.8) 式表示，其中 f⁻(z) 和 f⁺(z) 由 (22.7) 式给出，尺度分别为 σ = σ⁻ 和 σ = σ⁺。

增益函数产生 [0, 1] 范围内的值，其中 1 表示完美拟合。其得分反映样条路径上好点的百分比，可视为抽到好横坐标的概率估计。

为降低计算成本，使用两个 bailout 测试和数字样条。数字样条避免插值。第一个 bailout 测试拒绝沿图像列未良好展开的 n 个横坐标的样本，以保证样条模型的良好支撑。通过仅选择连续横坐标之间距离不小于规定距离 Δ = (m−1)/(2(n+1)) 的样本来获得每个样本中横坐标的良好展开，其中 n 是样本大小，m 是图像列数。对于 n = 5（4 个多项式），Δ = (m−1)/12。第二个 bailout 测试在样本中任何点的强度梯度和样条法线之间的夹角 δ 大于阈值 σ_δ 时拒绝任何 n 点样本。

3.2.3 Stopping Criterion

应检查的样本最小数 k 可通过将若干标准差 σ 加到获得好样本所需的预期样本数 ν 上来确定（[62]），即：

k > ν + Nσ = ω⁻ⁿ · [(1−ω)ⁿ/(ωⁿ)] + N·σ_ωⁿ (22.9)

其中 N 是加到均值的标准差数，n 是样本大小，ω 是数据集中内点的比例。

处理每个 n 个横坐标样本后，ω 更新为截至该时刻为增益函数找到的最高值。当横坐标样本数 k 超过 (22.9) 右端给出的数时，过程终止。图 22.5 显示了 MA 界面检测的两个示例。

3.3 Estimation of the Lumen–Intima Interface Using DP

所提出的 LI 边界检测方法使用动态规划，并假设颈动脉 MA 界面先前已被检测（使用第 3.2 节描述的方法）。该方法概述如下步骤：

若 (x,y) 是 ∇I(x,y) 方向上 ICOV 的局部最大值，则设 E_all(x,y) = 1；否则 E_all(x,y) = 0。E_all 是包含所有边缘的边缘图（图 22.6b），如第 3.1 节所述但不对边缘强度施加任何约束（由 ICOV 度量）。该边缘图很重要，因为管腔边界的一些边缘可能相当弱。
若 ICOV(x,y) > T，则设 E_strong(x,y) = E_all(x,y)，其中 T 是用稳健统计 [12] 自动估计的阈值；否则设 E_strong(x,y) = 0。E_strong 是仅包含强边缘的边缘图（图 22.6c）。
若 (x,y) ∈ Ω，其中 Ω 是由检测到的 MA 轮廓定界的内部区域，则设 E(x,y) = E_all(x,y)；否则 E(x,y) = 0。
若 θ(x,y) ≤ 90°，则设 E(x,y) = 0 且 E_strong(x,y) = 0，其中 θ 是 ∇I(x,y) 与 ∇D(x,y) 之间的夹角，D 是到管腔中轴的距离图。该步骤移除 E 和 E_strong 中具有指向动脉内部的梯度的边缘，即与 LI 界面不兼容的边缘。图 22.6d 说明了该步骤对 E_strong 的结果。
若 E(x,y) = 1，计算 ICOV(x,y) = ICOV(x,y)/max_y(ICOV(x,y))，即像素 (x,y) 在垂直方向的归一化 ICOV。ICOV 在每个横坐标的管腔轴上方和下方分别计算。该过程使通常较弱的管腔边界边缘有机会与 MA 界面边缘竞争。否则 DP 轮廓倾向于被吸引到 MA 界面（ICOV 强得多）。
若 E_strong(x,y) = 1 且在 (x,y) 与管腔轴之间的垂直方向 E_strong 中没有其他边缘，则设 E(x,y) = 1（图 22.6e, f）。这是因为当内膜-中膜区域不可见时，颈动脉 MA 界面是管腔边界的最佳估计。对这些边缘设 ICOV*(x,y) = 1。
使用 DP，在 E 中从 ROI 第一列到最后一列寻找使代价函数

C_t = Σⱼ₌₁ᴺ c(xⱼ, yⱼ) (22.10)

最小化的路径，其中若 E(xⱼ,yⱼ) = 1，则 c(xⱼ,yⱼ) = 1 − ICOV⁻*(xⱼ,yⱼ)；否则 c(xⱼ,yⱼ) = 1。N 是 ROI 的列数。DP 算法直接应用于 E，在管腔轴上方和下方分别独立搜索，以检测近壁和远壁 LI 界面。

3.4 Smoothing and Improving the Detected Lumen Boundaries

DP 算法的输出需要平滑，因为噪声常导致其不规则。该平滑未通过在 (22.10) 中添加几何项来实现，因为通常需要较强的几何项，这会导致轮廓在 LI 边界高曲率部分绕过组织区域。平滑通过几何 snake 引入，也可改善存在深凹陷或明显凸起时的检测，如图 22.6b 上管腔边界左部分所示。为此目的，[46, 47] 中引入了 Chan-Vese 分段常数分割模型 [63, 64] 的混合版本。

Chan-Vese 几何 snake 是成功的基于区域的水平集框架 [65, 66] 主动轮廓，不依赖于梯度，使其对主动轮廓的初始位置和边界中的小间隙具有稳健性。它还具有自动检测内部轮廓的能力，其拓扑变化以自然方式处理，其在空间中演化的方向由图像数据自动确定。

Chan-Vese 主动轮廓假设两阶段分段常数图像，由函数 F(c₁, c₂, C) 描述（式 22.11），其中 μ、λ₁ 和 λ₂ 是正常数，u₀: Ω → ℝ 表示输入图像，c₁ 和 c₂ 分别是 u₀ 在由 C 表示的区域边界内外的平均值。在 [63] 中为 (22.11) 引入了水平集公式，其中主动轮廓 C 被表示为隐函数的零水平集。常数 c₁ 和 c₂ 可被视为随主动轮廓在空间中演化的强度阈值的函数。

将该主动轮廓应用于颈动脉 B 型图像的主要困难在于其非分段常数特性。[46, 47] 中通过引入 Chan-Vese 模型的一些修改来解决该问题，使其也快于原始模型。Chan-Vese 模型的强度阈值估计被阈值化表面替代，保留其他有吸引力的特性。混合模型的实现可分为两个步骤：

在主动轮廓初始化之前，通过沿 DP 算法产生的轮廓上的边缘处的强度插值得到最优平滑阈值化表面 T_opt(x,y)。
然后图像由 Chan-Vese 两阶段分段常数主动轮廓的修改版本处理，其中 c₁ 和 c₂ 被确定为 T_opt(x,y) 的函数，使得每个像素的强度阈值由阈值化表面给出。

在混合模型中，c₁(x,y) = 2·T_opt(x,y) 且 c₂(x,y) = 0，以保持强度阈值作为 c₁ 和 c₂ 的平均值。c₁ 和 c₂ 的值在轮廓演化期间不更新。主动轮廓仅用于平滑并改善估计边界的定位精度。

[46, 47] 中引入一种新方法来估计阈值化表面 T_opt(x,y)，通过插值边缘像素处的图像强度，因为这些位置预期能找到好的局部阈值。如 [67] 中，插值表面通过求解 Laplace 方程 ∂²u/∂x² + ∂²u/∂y² = 0 (22.12) 得到插值表面 u。

(22.12) 式分两步计算。首先，使用快速算法获得近似解，将插值像素处的强度值传播到图像其余部分。然后使用以下数值格式求解 Laplace 方程：

uⁿ⁺¹(x,y) = 0.25·[uⁿ(x+1,y) + uⁿ(x−1,y) + uⁿ(x,y+1) + uⁿ(x,y−1)] uⁿ⁺¹(x,y) = β·uⁿ⁺¹(x,y) + (1−β)·uⁿ(x,y)；n = 0, 1, 2, … (22.13)

其中 uⁿ(x,y) 是像素 (x,y) 和迭代 n 处的阈值化表面值，β = 1.5，初始解 u⁰ 是第一步中的表面。

(22.13) 中的第二个等式引入超松弛 [61] 以加速收敛。当每个像素在两次连续迭代之间的相对强度差 < 1% 时迭代停止。

图 22.7 显示了最终平滑轮廓的两个示例。

3.5 Parameter Settings

参数要么自动估计（σₑ, T_A 等），要么通过实验测试确定。例如，ROI 大小尽可能小但足够大以包含颈动脉区域。

非线性平滑滤波器 (22.4) 使用 [68] 中提出的数值格式离散化。该格式的两个参数是图像网格大小 h（归一化为 1）和时间步 Δt（设为 0.25 以获得良好收敛速度且不失数值稳定性）。边缘尺度 σₑ（式 22.3）的值以及相应曲线的斜率随迭代次数增加而减小。当 σₑ 的变化率降至 5×10⁻⁴ 以下时迭代停止，之后平滑增加非常缓慢。

在 RANSAC 相关文献中（例如 [62, 69]），选择参数 N（式 22.9）使 1/N² ≤ 0.05 是常见的，这给出 N ≤ 4.5；因此选择 N = 5。

在所用数据集中，颈动脉 MA 界面像素到管腔中轴的距离小于 70 像素，因此使用 d_max = 90 像素以保留颈动脉边界外的一些图像数据。

θ 的值在每列图像中沿手工描迹的颈动脉边界对最接近这些边界的边缘像素计算。观察到在至少 99% 的情况下 θ < 30°，这意味着 θ_max = 30° 是合适的值。

参数 SDL_min 表示 SDL 的阈值，低于该阈值不预期能找到颈动脉 MA 界面的边缘像素。观察到数据集中每张图像 SDL > −6.4，因此选择 SDL_min = −7 和 σ⁻ = −SDL_min = 7。参数 SDL_max 和 σ⁺ 表示 SDL 的阈值，高于该阈值不预期能找到 MA 界面的边缘像素。该阈值的自然值是 d_max，因此 SDL_max = σ⁺ = d_max。

数据集中谷形边缘的宽度测量显示 L = 10 足以满足属于 MA 界面的谷形边缘。

为估计 σ_d 和 σ_δ 的值，对一位专家手工描迹的 MA 界面（MA1）的每点计算 dₑ 和 δ 的值。参数 δ 由使用局部主导梯度方向滤波器从图像获得的梯度图计算。参数 dₑ 由第 3.1 节描述的边缘图计算。参数 σ_d 和 σ_δ 是 Tukey 函数的尺度，因此可分别计算为 σ_d = √5·[C·MAD(dₑ) + med(dₑ)] 和 σ_δ = √5·[C·MAD(δ) + med(δ)]，给出 σ_δ ≈ 11° 和 σ_d ≈ 4。

阈值 σ_δ 还用于其中一个 bailout 测试中，作为强度梯度与拟合样条法线之间夹角 δ 的限制。

混合 Chan-Vese 模型中有若干参数。除确定轮廓弹性的参数 μ 外，其他所有参数都按 [63, 64] 中的建议设置。通过选择 λ₁ = λ₂ = 1 给图像中两个阶段同等重要性。使用 [63, 64] 中提出的半隐式数值格式求解 (22.11)，设 Δx = Δy = 1（其中 (Δx, Δy) 表示图像网格大小），时间步 Δt = 0.1·Δx。该数值格式迭代直到主动轮廓在两次连续迭代之间覆盖的最大距离小于 0.01·Δx。

由于参数 μ 可被视为尺度参数，因此将其定义为分割管腔尺度的函数，由 DP 算法产生的 LI 界面长度 L 度量。在纵向截面中，L 是由两个估计轮廓（一个在中轴上方，一个在中轴下方）定界区域的周长。尺度参数设为 μ = ν·L²，其中 ν = 10⁻⁴ 经验确定为满意平滑水平与数据保真度之间的折衷。[63, 64] 中也使用的因子 255² 是必要的，以保持水平集方程中的一致性，因为图像强度表示在 {0, …, 255} 范围内。

4 Image Database

图像数据库由 47 张颈总动脉纵向 B 型图像组成，使用 Philips HDI 5000 超声系统采集，以 256 灰度级记录。图像集来自 24 名不同的症状性患者，包含多种斑块类别（II 至 IV 类）、大小和形状。像素大小归一化为 0.09 mm，这是临床实践中的常用值。

所有图像由两位医学专家 A 和 B 以及其中一位专家（A）在相隔 1 年的两个时刻各手工分割一次。这些分割由 LI 和 MA 界面的手工描迹组成，被视为真值。此后专家 A 和专家 B 的手工分割以及自动分割将分别称为 MA1、MA2、MB1 和 A。

5 Evaluation Methodology

[46, 47] 中评估时使用了若干统计量和统计分析。它们基于每条手工和自动轮廓共有横坐标处轮廓之间的垂直距离。

一个统计量是 D_max，即两次不同检测获得的对应轮廓之间的最大垂直距离。另一统计量是 D_mean，即两次不同检测获得的同源轮廓之间的平均垂直距离。其他使用的统计量是 IMT_min、IMT_mean 和 IMT_max，分别是每个分割的内膜-中膜区域的最小、平均和最大 IMT。

由 IMT_min、IMT_mean 和 IMT_max 测量的差异的合并均值 N_x 和标准差 sd 计算得到方法间（手工 vs 自动）、观察者内和观察者间误差，定义为 se = sd/√2 [21]。IMT_min、IMT_mean 和 IMT_max 的变异系数 CV 也使用以下公式估计：

CV = 100·se/N_x % (22.14)

方法间和观察者间的一致性通过箱线图和 Bland-Altman 图 [70] 评估。

6 Results

图 22.8 显示了内膜-中膜区域良好分割的示例，图 22.9 显示了缺陷分割的示例。为比较起见，还显示了相应的 MA1 分割。

在测试的 47 张 CCA B 型图像中，86.2% 的 MA 边界（94 个中的 81 个）在自动轮廓与任一相应手工版本之间呈现 D_max < 1 mm。图 22.10 显示了 D_max < 1 mm 的 MA 界面检测的 D_max 统计分布。

LI 界面的统计分析仅针对 MA 界面 D_max < 1 mm 的 81 个内膜-中膜复合体子集计算。对于该 LI 边界子集，81.5%（上述 81 个中的 66 个）呈现 D_max < 1 mm。D_max < 1 mm 的 LI 界面子集的 D_max 和 D_mean 分布可在图 22.11 和 22.12 中找到。

以下统计分析（评估观察者内、观察者间和方法间误差以及方法间一致度）仅针对 D_max < 1 mm 的 66 个 LI 边界计算。

表 22.1 显示了方法间（手工 vs 自动）、观察者内和观察者间误差，以及 IMT_min、IMT_mean 和 IMT_max 的变异系数 CV。

图 22.13 显示了两种给定分割的平均 IMT 差异相对其平均值的 Bland-Altman 图。标准差用 SD 表示，一致性限值（形式为 mean ± 2SD mm）在 MB1 和 A 之间为 0.03 ± 0.22，在 MA2 和 A 之间为 −0.07 ± 0.24，在 MA2 和 MB1 之间为 −0.10 ± 0.20，在 MA1 和 A 之间为 0.11 ± 0.24，在 MA1 和 MB1 之间为 0.08 ± 0.24，在 MA1 和 MA2 之间为 0.18 ± 0.24。

使用 Matlab 在配备 2.13 GHz Intel Core 2 Duo 处理器的计算机上，估计的中位计算时间为：DP 算法 0.0470 秒，Chan-Vese 算法 18.5 秒，整个 LI 界面分割过程 18.9 秒。分割 MA 界面所需的中位计算时间为 28.5 秒。

7 Discussion

如上一节所示，[46, 47] 中提出的方法能够分割 CCA 的 B 型图像中的颈动脉壁，即使存在不同大小、形状和类别的斑块。它可以分割远壁和近壁，仅需要最少的用户交互。自动检测的界面未进行任何手工校正。自动测量的误差分散略大于手工测量，特别是在低质量图像和近壁的情况下。然而统计分析显示了有希望的结果。

对于 MA 界面，D_max 的分布（图 22.10）引出三个主要观察。首先，近壁边界的 D_max 值及其变异性大于远壁边界，这与这些边界的可见性较低以及医学专家和计算机检测的较大难度一致。其次，当其中一个对比轮廓是自动轮廓时 D_max 倾向于更高，虽然差异通常很小，对远壁边界尤其小。在至少 75% 的近壁案例中，对比两个手工 MA 轮廓时 D_max < 0.40 mm，而当其中一个对比轮廓是自动的时 D_max < 0.51 mm。对远壁 MA 边界考虑相同百分位，手工轮廓对的 D_max < 0.40 mm，其他对的 D_max < 0.43 mm。第三，某些自动检测距离手工版本过远。在这些情况中的某些情况下，不同手工轮廓之间的距离也很高，可能由于图像质量较差。然而 13.8% 的 MA 边界自动检测 D_max > 1 mm。

对于 LI 界面，检测结果（图 22.11 和 22.12）在近壁比远壁更差。与 MA 界面相比 LI 界面的 D_max 较大值是 LI 边界可见性较差的结果。这也解释了为何 LI 界面自动检测中 D_max > 1 mm 的百分比（18.5%）高于 MA1。

IMT_min 的变异性（表 22.1 中 se 和 CV 值）大于 IMT_mean 和 IMT_max 可由一些非常薄的内膜-中膜区域中 LI 边界的较差定义来解释。当 LI 边界靠近 MA 界面处存在间隙时，专家可以推断 LI 界面的正确位置，而自动过程仅寻找最近的兼容边缘，因此相应的方法间变异性较高。与 IMT_min 相比，IMT_mean 和 IMT_max 的方法间变异性对动脉粥样硬化的诊断更重要，仅略大于手工分割的同源变异性。

手工分割的 IMT 均值的 Bland-Altman 图（图 22.13）显示出高度的观察者内和观察者间一致性，平均差异接近零（−0.10 mm ≤ mean ≤ 0.18 mm），标准差值较小（0.10 mm ≤ SD ≤ 0.12 mm），几乎所有差异都在一致性限值内（[mean − 2SD; mean + 2SD]）。这些图显示自动分割和手工分割之间良好的一致性，平均差异值（−0.07 mm ≤ mean ≤ 0.11 mm）和标准差值（0.11 mm ≤ SD ≤ 0.12 mm）相似。

8 Summary

本章综述了颈动脉超声图像分割的方法学。还描述了一种在 B 型超声图像中分割 CCA 近壁和远壁内膜-中膜区域的方法。MA 界面通过将三次样条最佳拟合到具有与 MA 边界兼容特征的边缘来检测。使用 MA 界面的若干判别特征来减少自动轮廓被其他边缘吸引。通过动态规划、平滑强度阈值化表面和几何 snake 获得 LI 界面的平滑估计。

该方法使用来自 24 名不同症状性患者的 47 张图像（包括 II–IV 类、大小和形状的多种斑块）的图像集进行统计评估。结果表明所提出的方法能够产生准确度可与医学专家手工描迹相比的分割。

尽管结果有希望，但在质量非常差的图像中 LI 界面的检测可能失败。在许多情况下，斑块只能用补充信息（如功率多普勒成像或其他角度的 B 扫描）检测。因此，可能需要集成此类补充数据，以将自动检测误差降低到与专家手工描迹相当的水平。

本章个人批注

本章是这本书中第一个明确把"算法 + 实现细节 + 实验验证"完整呈现的技术性章节，与前面几章偏向临床或影像学的风格形成对照。从写作结构上看，作者基本上是按照"经典方法综述 → 提出自己的方法（cubic spline + RANSAC + DP + Chan-Vese）→ 数据库描述 → 评估方法论 → 结果 → 讨论 → 总结"的 IEEE 期刊论文模板展开，对一个想复现该工作的研究者非常友好。

值得注意的是 §2 综述部分基本把 1998 到 2011 年间颈动脉 IMT 分割领域的主流方法按类别分块呈现（edge detection、gray level density、dynamic programming、snakes、Hough、template matching、Nakagami、feature extraction、watershed），每类给出代表性文献和关键缺陷。这其实是一份不错的领域入门地图。缺陷的总结也很有规律性：半自动方法需要用户交互、计算负担重的自动方法不适合临床、不考虑斑块的方法缺乏实用性、近壁 IMT 测量研究的稀少。这些限制基本都对应作者自己方法的卖点。

§3 是本章的核心。从工程实现角度，方法包含四个串行模块：edge estimation（含 ICOV 滤波器、主导梯度方向、最终边缘图、谷形边缘图）→ RANSAC cubic spline 拟合 MA → DP 检测 LI → Chan-Vese 混合模型平滑 LI。每一节都给出关键公式（如 ICOV 的定义 (22.1)、RANSAC 停止准则 (22.9)、DP 代价函数 (22.10)、Laplace 方程数值解 (22.13)）。作为读者，如果想实现，§3.5 几乎把全部参数都给到了具体数值（包括 ROI 大小策略、Δt、β、ν、SDL_min、σ_d、σ_δ 等），这在实际工程文献里很少见。

Bland-Altman 的数值（mean ± 2SD）说明自动分割与手工描迹的一致性基本和两位专家之间的一致性相当（−0.07 到 0.11 mm vs −0.10 到 0.18 mm），这是评估 IMT 自动算法最关键的指标之一。远壁 IMTmean 的 CV 在 8% 左右，与手工测量的变异系数在同一量级。这个数值水平对临床应用是可接受的。

但 §7 的讨论也提到 13.8% 的 MA 检测和 18.5% 的 LI 检测 D_max > 1 mm，说明在困难样本上仍有失败，且对近壁和低质量图像的鲁棒性不足。本章没有给出针对失败模式的细致分类或针对性改进，只提到"未来可能需要集成多普勒或多角度 B 扫描"——这其实是一个开放的工程问题。

本章与上一章（Ch21 关于 level set 测量 IMT）形成有趣对照：Ch21 主要解决 NMI（normalised mean intensity）类特征驱动的分割，Ch22 是经典的 edge-based + region-based 混合方案。两者都依赖手工 ROI 初始化，都报告远壁比近壁准确度高很多。这两个反复出现的现象（ROI 初始化需求、近壁弱）实际上贯穿了本书中所有 IMT 自动测量方法，是这个领域的结构性局限。

与上下章的衔接（一段话）

本章紧接 Ch21 同样以 IMT 为目标的算法章节之后，从 level set/NMI 路线切换到 edge-based + region-based 的混合路线；同时本章是"经典算法综述 + 自主方法"的典型模式，作为后续超声自动分割章节（Ch17–Ch19）的算法基础——前几章侧重具体实现（CALEX 测量、CARS 识别、plaque 分割），本章则提供了从 1998 到 2011 的完整方法地图。本章也是 Part IV "Ultrasound" 部分的最后一章，下一章（Ch23）转入主动脉瘤的影像学，主题从 IMT 测量转向血管几何形态评估，标志着 Part IV 从"颈动脉粥样硬化"焦点转向其他大血管病变。