基于颈部超声图像的主颈动脉血管提取和厚度测量方法

摘要

本发明公开了一种基于颈部超声图像的主颈动脉血管提取和厚度测量方法，具体为：读取颈部超声三维体数据，基于主颈动脉分叉点标记主颈动脉中轴；对颈部超声三维体数据按三视图方向，依次投影切分得到二维横断面、冠状面和矢状面序列图像；分别对二维横断面、冠状面和矢状面序列图像做预处理、分割、重建或内中膜厚度统计，得到最终的颈部超声主颈动脉血管壁内外轮廓和血管壁厚度等相关信息。本发明克服现有计算机辅助诊断中血管分割方法计算复杂度大，不能准确测量血管管壁厚度，主观因素易造成误差等缺点，能够完整、快速、准确地得到颈部超声主颈动脉血管壁内外轮廓和血管壁厚度。本发明与手动分割方法相比，操作快捷，可用于颈部粥样硬化以及心血管疾病的辅助诊断和防治。

说明书

技术领域

本发明属于生物医学工程和图像处理交叉领域，具体涉及一种基于颈 部超声图像的主颈动脉血管提取和厚度测量方法。

背景技术

世界卫生组织最新统计数据显示，心血管疾病(Cardiovascular Diseases， CVDs)是世界上三大致死疾病之一。动脉血管壁增厚、形成斑块，进而导致 血管狭窄，是动脉粥样硬化的典型病征之一。颈动脉粥样硬化是一种能引 发心脏病、中风的心脑血管疾病，严重危害了人类身体健康。因此，对其 早期预防、诊断、治疗和监控有着重要意义。

2010年1月20日，美国心脏学会（American Heart Association，AHA） 战略规划工作委员会发布的十年健康战略目标，即AHA 2020健康战略—— 《定义和制定促进心血管健康和减少疾病的国家目标》，首次提出“以改善 健康水平为主要目标”，标志着美国心脏学会AHA将心血管疾病CVDs的 预防战线进一步前移，不仅针对已具有危险因素的高危人群和患病人群， 而且要改善普通人群的健康水平。研究表明，预防和控制心血管疾病CVDs 的关键是早检测、早治疗。因此，对动脉粥样硬化的早期检测与防治，对 降低心血管疾病CVDs死亡率有着极其重要的临床意义。

主颈动脉(Common Carotid Artery，CCA)的血管壁增厚程度可作为衡 量病变的重要指标。超声成像技术所特有的“实时、经济、可靠、安全” 优点，使得基于该技术的内中膜厚度(Intima-media Thickness，IMT)成为评 估颈动脉粥样硬化程度的常用指标之一。超声图像血管提取和厚度测量成 为近年的研究热点。

首先，就血管提取而言，中国专利申请号为200910106119.6和 201010297322.9的两个专利提出了超声图像血管提取的方法，前者针对非 序列单张二维超声血管灰度图像，后者应用于血管内超声序列图像，两者 均缺乏最终详细的血管信息。如：未针对序列图像有针对性地完成分割和 方法评价等工作、无法根据主颈动脉CCA内部解剖结构获得最接近真实的 血管厚度值、在国产超声机软件升级换代中难以产业化应用。

其次，就厚度测量而言，CULEX(Completely User-independent Layer  Extraction)和CALEX(Completely Automatic Layer Extraction)为目前最新颖、 最智能的厚度测量方法。两者均可以达到全自动厚度测量，但实现难度大、 计算复杂度高。CULEX和CALEX利用不一样的图像特征，采用完全不 同的思想——前者利用图像像素的局部统计值来分辨血管腔内像素和组织 像素，分割方法则是基于梯度和活动模型的结合；而后者集特征提取、线 性拟合和分类于一体。对比两种方法的分割效果可以看出，CALEX对血管 腔-内膜LI轮廓的分割不完整，但对中膜-外膜MA轮廓的分割优于CULEX； 同时CULEX受图像噪声和图像伪影影响较大，而且CALEX的执行效率远 高于CULEX。

发明内容

本发明的目的在于提供一种能全方位、多角度、快速、准确且易于实 现、操作方便的主颈动脉CCA血管提取和厚度测量方法。

基于颈部超声图像的主颈动脉血管提取和厚度测量方法，包括以下步 骤：

读取颈部超声三维体数据，基于主颈动脉分叉点标记主颈动脉血管中 心轴；

依据中心轴，按三视图方向投影，切分颈部超声三维体数据，得到二 维横断面、矢状面和冠状面序列图像；

在二维横断面序列图像的每一张图像中，分别选取各主颈动脉血管感 兴趣区域，并对各主颈动脉血管感兴趣区域预处理；在预处理后的各主颈 动脉血管感兴趣区域中，分别分割得到各主颈动脉血管的内、外轮廓；依 据各二维横断面图像的主颈动脉血管感兴趣区域的内、外轮廓，按其空间 位置关系三维重建，得到三维主颈动脉血管轮廓；

在二维矢状面序列图像的每一张图像中，分别选取各主颈动脉血管感 兴趣区域，并对各主颈动脉血管感兴趣区域预处理；在预处理后的各主颈 动脉血管感兴趣区域中，分别分割得到各主颈动脉血管的内、外轮廓；分 别计算各主颈动脉血管的内、外轮廓之间的厚度即内中膜厚度；统计各二 维矢状面图像的内中膜厚度均值；

在二维冠状面序列图像的每一张图像中，分别选取各主颈动脉血管感 兴趣区域，并对各主颈动脉血管感兴趣区域预处理；在预处理后的各主颈 动脉血管感兴趣区域中，分别分割得到各主颈动脉血管的内、外轮廓；分 别计算各主颈动脉血管的内、外轮廓之间的厚度即内中膜厚度；统计计算 各二维冠状面图像的内中膜厚度均值；

进一步地，在所述二维横断面序列图像的每一张图像中的主颈动脉血 管感兴趣区域中，按照如下方式分割得到各主颈动脉血管的内、外轮廓：

A1、在二维横断面图像预处理后的主颈动脉血管感兴趣区域中，提取 主颈动脉血管的内轮廓，并做平滑处理；

A2、依据主颈动脉血管的内轮廓进行椭圆拟合，将拟合得到的椭圆放 大后，作为主颈动脉血管的初始外轮廓；

A3、依据主颈动脉血管的初始外轮廓，演化得到主颈动脉血管的最终 外轮廓。

所述步骤A1具体为：在二维横断面图像预处理后的主颈动脉血管感兴 趣区域中，提取主颈动脉血管的圆心位置和半径，依据圆心位置和半径提 取内轮廓。

所述步骤A2中椭圆放大比例系数为1.02~1.08。

所述步骤A3采用的演化方法为Active Contour Model、GVF-Snake、Fast Marching Method、Level Set、Sparse Field Algorithm任意一种。

进一步地，在所述二维矢状面序列图像的每一张图像中的主颈动脉血 管感兴趣区域中，按照如下方式分割得到各主颈动脉血管的内、外轮廓：

B1、对二维矢状面图像预处理后的主颈动脉血管感兴趣区域内的每一 列的点，从上到下依次进行序号标记作为每一列点的索引值；

B2、感兴趣区域的粗分割：

B21、获取下边界轮廓DB：从左到右，顺序扫描二维矢状面图像的主 颈动脉血管感兴趣区域的每一列，获得每一列的最大灰度值、最小灰度值 以及满足阈值条件的下边界候选点；在下边界候选点中，标记索引值最大 的点作为唯一下边界轮廓点；顺次连接所有下边界轮廓点构成下边界轮廓 DB；下边界候选点灰度值GV_candi应满足的灰度阈值条件为： GV_candidate≥a×ΔGV+GV_min＝a×(GV_max-GV_min)+GV_min，权重系数a的取值范围为 0.87~0.93；

B22、获取上边界轮廓UB：将下边界轮廓DB在二维矢状面图像的主 颈动脉血管感兴趣区域中向上平行移动20~30个像素作为上边界临时轮廓 UB_t；选取大小为X×X的模板，X取8~15像素，利用模板从左到右，从上 到下对上边界临时轮廓UB_t和下边界轮廓DB之间的区域进行遍历，计算 模板覆盖区域的像素均值EX和方差DX；若满足均值方差阈值条件 EX≤b且DX≤c，b取值范围为0.05~0.09，c取值范围为0.11~0.15，则模板覆 盖区域的中心为上边界候选点；若不满足均值方差阈值条件，则对应UB_t 轮廓点为上边界候选点；在每一列的上边界候选点中，标记索引值最大者 作为唯一上边界轮廓点；顺次连接所有上边界轮廓点构成上边界轮廓UB；

B3、感兴趣区域的细分割：在下边界轮廓DB和上边界轮廓UB之间 的区域细分割出血管腔、内膜和中膜、外膜和外膜以外组织；

B4、提取血管腔与内膜的交界面即为内轮廓，提取中膜与外膜的交界 面为外轮廓。

所述感兴趣区域的细分割步骤B3中采用C均值、模糊C均值、支撑 向量机SVM、AdaBoost算法中的任意一种。

进一步地，在所述二维冠状面序列图像的每一张图像中的主颈动脉血 管感兴趣区域中，按照如下方式分割得到各主颈动脉血管的内、外轮廓：

C1、采用数学形态学方法在二维冠状面图像预处理后的主颈动脉血管 感兴趣区域中初始化内轮廓，依据初始化内轮廓演化生成最终的内轮廓；

C2、通过初始内轮廓下移得到初始外轮廓，依据初始外轮廓演化生成 最终的外轮廓；

所述步骤C1和C2采用的演化方法为经典蛇形算法、GVF-Snake算法、 水平集、MS模型、CV模型演化方法中的任意一种。

本发明的有益效果：

与传统的血管提取和厚度测量方法相比，本发明提供的基于颈部超声 图像的主颈动脉CCA血管提取和厚度测量方法，有四点区别：

（1）不同于传统方法在某一个面上进行血管提取和厚度测量，本发明 将颈部超声三维体数据，按三视图方向投影切分，分别得到二维横断面、 矢状面和冠状面的序列图像，并在各投影面上分别进行处理；不同于传统 方法在某一个面的单帧图像上进行血管提取和厚度测量，本发明有针对性 地对三个面的序列图像分别进行处理，经血管分割、提取、测量和计算后， 可充分利用原有的完整三维数据信息，获得全方位的血管参数，如：血管 厚度、面积、体积和斑块大小、数量等；

（2）在二维横断面序列图像的每一张图像中的主颈动脉血管感兴趣区 域中，分割得到各主颈动脉血管的内轮廓时，引入Hough变换圆检测，获 得血管圆心位置和半径参数，从而可指导后续数学形态学中的结构基元大 小改变，优化内轮廓分割结果；在二维横断面序列图像的每一张图像中的 主颈动脉血管感兴趣区域中，演化得到各主颈动脉血管的外轮廓前，引入 椭圆拟合策略作为先验知识，并将拟合的椭圆作为初始外轮廓，更符合血 管的生理形状；在二维横断面序列图像的每一张图像中的主颈动脉血管感 兴趣区域中，演化得到各主颈动脉血管的外轮廓时，引入主动轮廓模型 ACM方法，以解决血管外轮廓的弱边界、难以区分的问题；

（3）在二维矢状面序列图像的每一张图像中的主颈动脉血管感兴趣区 域中，分割得到各主颈动脉血管的内、外轮廓时，采用“粗-细”两层分割 结构；在“粗分割”时，利用灰度阈值条件和均值方差阈值条件加以约束； 在“细分割”时，首先引入模糊C均值FCM聚类算法；其次，将血管组织 划分为多类；最后再分别合并为三类——“血管腔”、“内中膜”和“外 膜及外组织”，获得各交界面轮廓，以提高矢状面内中膜厚度测量的准确 度；

（4）在二维冠状面序列图像的每一张图像中的主颈动脉血管感兴趣区 域中，分割得到各主颈动脉血管的内、外轮廓时，利用数学形态学方法获 得其初始内、外轮廓；再通过演化得到最终内、外轮廓；在二维冠状面图 像的内中膜分割中，在Snake算法的基础上，再次引入GVF-Snake算法， 解决深度凹陷问题，以提高冠状面内中膜厚度测量的精确度。

综上，本发明提供的基于颈部超声图像的主颈动脉CCA血管提取和厚 度测量方法，首先能有效地应对超声图像的多样性，达到准确、快速分割 血管和测量厚度的目的；其次，经定量分析与比较，本方法与手动分割、 测量方法误差相当，其应用可直接减轻医务工作者海量的图像手动标记工 作量；最后，基于本发明方法分析所得的临床参数（如：血管厚度、面积、 体积和斑块大小、数量等），不仅可以直观定量地反映血管病变，而且能 为颈动脉粥样硬化的早期防治提供更多的指导意义。

附图说明

图1为本发明中的具体步骤流程图；

图2为本发明中二维序列图像处理的流程图；

图3为本发明具体步骤的流程图；

图4为标记三维体数据3D_Data分叉点和中轴方向示意图；

图5为主颈动脉根据中轴三视图投影切分示意图；

图6为二维横断面序列图像（2D_Data_Z）；

图7为二维横断面序列图像中的第13张原始图像（2D_Data_Z_k，k=13）；

图8为二维矢状面序列图像（2D_Data_Y）；

图9为二维矢状面序列图像中的第2张原始图像（2D_Data_Y，j＝2）；

图10为二维冠状面序列图像（2D_Data_X）；

图11为二维冠状面序列图像（2D_Data_X）中的第1张原始图像 （2D_Data_X_i，i=1）；

图12为二维横断面序列图像分割流程图；

图13为该二维横断面图像的感兴趣区域ROI图；

图14为该感兴趣区域ROI预处理逐步结果图：图14-（a）为分段线性 拉伸结果图；图14-（b）为SRAD非线性滤波结果图；

图15为对预处理感兴趣区域Hough变换圆检测的结果图；

图16为对预处理感兴趣区域Canny算子边缘检测的结果图；

图17为数学形态学处理结果图；

图18为提取的单像素点内轮廓结果图；

图19为内轮廓展开点集图；

图20为内轮廓展开点集的拟合曲线图；

图21为内轮廓复原点集图；

图22为该二维横断面内轮廓分割结果图：其中实心圆形点表示手动分 割轮廓；实线表示本发明分割轮廓；

图23为内轮廓椭圆拟合结果图；

图24为该二维横断面外轮廓分割结果图：其中实心菱形点表示手动分 割轮廓；实心方形点表示本发明分割轮廓；

图25为二维矢状面预处理结果图：图25-（a）为归一化结果图；图25- （b）为滤波结果图；

图26为二维矢状面感兴趣区ROI预处理结果图；

图27为二维矢状面感兴趣区ROI粗分割结果图：图27-（a）为原始大 小；图27-（b）为放大一倍显示图；图27-（c）为放大两倍显示图；

图28为二维矢状面感兴趣区ROI细分割结果图：图28-（a）为原始大 小；图28-（b）为放大一倍显示图；图28-（c）为放大两倍显示图；

图29为该二维冠状面提取感兴趣区域ROI结果图：图29-（a）为原始 图；图29-（b）为对应ROI图；

图30为该二维冠状面感兴趣区域ROI预处理结果图；

图31为该感兴趣区域阈值化结果图；

图32为该感兴趣区域填补空洞结果图；

图33为该感兴趣区域内膜表面去毛刺结果图；

图34为该感兴趣区域内膜修正结果图；

图35为内膜的初始轮廓（虚线）和最终轮廓（实线）结果图：图35- （a）为原始大小；图35-（b）为放大一倍显示图；图35-（c）为放大两倍 显示图；

图36为外膜的初始轮廓（虚线）和最终轮廓（实线）结果图：图36- （a）为原始大小；图36-（b）为放大一倍显示图；图36-（c）为放大两倍 显示图；

图37为内膜（上）和外膜（下）的初始轮廓（虚线）和最终轮廓（实 线）结果图：图37-（a）为原始大小；图37-（b）为放大一倍显示图；图 37-（c）为放大两倍显示图；

图38为该二维冠状面手动测量内中膜厚度金标准结果图。

具体实施方式

下面结合具体的实施示例及附图说明，对本发明做进一步的详细描述， 并给出了本发明的方案验证。

一种基于颈部超声图像的主颈动脉CCA血管提取和厚度测量方法，包 括以下五个步骤和方案验证，如图1所示；具体的，对每个二维序列图像 处理的流程图，如图2所示；总流程图，如3所示。

步骤（1）读取颈部超声三维体数据3D_Data，标记主颈动脉分叉点“BF” 和中心轴，如图4所示。

（1.1）定位分叉点“BF”

在三维超声体数据3D_Data中，识别主颈动脉CCA、内颈动脉ICA与 外颈动脉ECA，并确定颈动脉分叉点“BF”的位置（参见图4）。

（1.2）定位主颈动脉CCA的中心轴

参考分叉点“BF”的位置，根据主颈动脉CCA生理解剖结构，初步估 计中心轴位置，在估计的中心轴上任意选择两点——“点1”和“点2”（参 见图4），“点1”和“点2”的连线作为主颈动脉CCA的中心轴。

步骤（2）依据中心轴，按三视图方向投影，如图5所示，切分颈部超 声三维体数据3D_Data，分别得到二维横断面序列图像（2D_Data_Z）、二 维矢状面序列图像（2D_Data_Y）和二维冠状面序列图像（2D_Data_X）；

（2.1）二维横断面序列图像(2D_Data_Z)

按垂直于颈动脉中心轴方向进行切分，得到二维横断面序列图像。一 般的，二维横断面序列图像至多有13~15张，如图6所示；本发明中，取 第13张作为示例说明，如图7所示。

（2.2）二维矢状面序列图像(2D_Data_Y)

按平行于颈动脉中心轴方向、从侧面投影进行切分，得到二维矢状面 序列图像。一般的，二维矢状面序列图像至多有4~6张，如图8所示；本 发明中，取第2张作为示例说明，如图9所示。

（2.3）二维冠状面序列图像(2D_Data_X)

按平行于颈动脉中心轴方向、从正面投影进行切分，得到二维冠状面 序列图像。一般的，二维冠状面的序列图像至多有3张，如图10所示；本 发明中，取第1张作为示例说明，如图11所示。

步骤（3）依据二维横断面序列图像，提取三维主颈动脉血管轮廓。本 步骤分割各二维横断面序列图像的内、外血管轮廓，具体分割步骤流程图， 如图12所示。

在本步骤中，选择二维横断面序列图像（2D_Data_Z）中的一张为例说 明本步骤（参见图7），其它图像均采用同样的方法加以处理，最终获得所 有二维横断面序列图像的血管内、外轮廓，从而进行面积、体积等参数计 算，提供临床诊疗指标依据，定性地衡量临床病情。鉴于超声序列图像质 量整体较差，需经由预处理加以改善；血管提取的关键就是内、外轮廓的 分割，本发明中结合血管形态，引入圆、椭圆作为形状先验信息，更加有 效的分割血管；为了更好的呈现分割结果，最后再有针对性对内轮廓进行 后处理。

（3.1）在二维横断面序列图像（2D_Data_Z）中，分别选取各主颈动 脉血管感兴趣区域（2D_Data_Z_k_ROI），并对各主颈动脉血管感兴趣区 域预处理。预处理的技术思路为：先灰度变化，再滤波降噪。

原始的二维横断面图像灰度偏暗、对比度低、噪声较大，因此需要对 其做预处理。本发明是为了提取血管，因此只需将二维横断面图像中包含 血管的子区域作为感兴趣区域ROI进行预处理，图7的感兴趣区域ROI如 图13所示，每个预处理步骤结果如图14所示，预处理过程具体如下：

（3.1.1）灰度变换（图14-a）

灰度变换是为了提高图像处理时灰度级的动态范围，以提高图像的亮 度和对比度，可采用线性拉伸、非线性拉伸、图像增强等方法。本实例选 用最简单的分段线性变换函数进行说明，该方法由两个基本操作组成：

（3.1.1.1）确定对图像进行灰度拉伸的两个拐点；

对二维横断面图像作直方图统计，得到L灰度级的图像。把灰度变换 处理前后的灰度值分别用r和s分别定义，L灰度级的图像的灰度变换函数 表示为s＝T(r)。假设P1，P2是分段线性变换函数T(r)的两个拐点，其变换前 后的灰度值分别为（r₁，s₁）和（r₂，s₂）。

（3.1.1.2）灰度变换

本实例统计了10名匿名病人三维体数据，共获得二维横断面图像300 张，经统计分析后，设置所关注的原图的灰度级范围为[r_min,r_max]，则P1， P2两个拐点把变换函数T（r）分成3段：r_min≤r＜r₁，r₁≤r≤r₂，r₂＜r≤r_max。T(r) 的表达式如下所示：

$T\left(r\right)=\left(\begin{array}{c}(r-r_{\min })\times \frac{s_1}{r_1},r_{\min }\le r\le r_1\\ s_1+(r-r_1)\times \frac{s_2-s_1}{r_2-r_1},r_1\le r\le r_2\\ s_2+(r-r_2)\times \frac{(L-1)-s_2}{r_{\max }-r_2},r_2<r\le r_{\max }\end{array}\right)$

如此一来，就把灰度级由原来所关注的范围线性拉伸至饱和范围 [0，L-1]，L为变换后的灰度级。图13的灰度变换结果如图14-a所示。

（3.1.2）滤波降噪（图14-b）

滤波降噪可选用非局部均值、高斯滤波、斑点噪声各向异性扩散模型 SRAD等方法。其中，SRAD非线性滤波器由于在不同的扩散方向上采用不 同的扩散系数，故具有增强对比度、保留细节等优点。SRAD滤波器是基 于局部统计滤波器以及各向异性平滑提出的，相比于其他传统的基于局部 统计滤波器以及各向异性平滑滤波器，具有更好的同质区域平滑性能以及 更好地保留细节和边界增强的优点。图14-a的滤波降噪结果如图14-b所示。

（3.2）在预处理后的各主颈动脉血管感兴趣区域中，分割得到各主颈 动脉血管的内、外轮廓。

主颈动脉血管的内、外轮廓提取的技术思路为：首先提取内轮廓；其 次，拟合得到初始外轮廓；最后，演化得到主颈动脉血管的最终外轮廓。

（3.2.1）在预处理后的各主颈动脉血管感兴趣区域中，提取各主颈动 脉血管的内轮廓，并平滑处理。

提取主颈动脉血管的内轮廓的技术思路为：首先提取血管参数信息， 其次，依据血管参数信息，提取血管单像素点内轮廓；最后，对内轮廓作 平滑处理。

（3.2.1.1）在预处理后的各主颈动脉血管感兴趣区域中，提取血管参数 信息。

本步骤是为了提取血管的参数信息即圆心位置和半径，本实例采用霍 夫Hough变换圆检测方法示例说明（不局限该方法）。Hough变换圆检测 中可获得O(centre,R_max,R_min)，分别为圆心坐标，最大半径、最小半径，从而为 后续步骤的基元尺寸大小提供参考依据。Hough变换的实质是将图像空间 内具有一定关系的象元进行聚类，建立能把这些象元用某一解析形式联系 起来的参数空间，寻找累计对应点。在Hough变换检测圆中，由Canny算 子得到的边界点即为象元，圆心坐标和圆半径这三个参数即为其相应解析 形式的参数。另在本例中，为了减少程序运行时间，人为地选定一定的圆 半径范围和圆心位置。图15给出了对图14-b进行变换圆检测的结果。

（3.2.1.2）依据血管参数信息，提取血管单像素点内轮廓

基于数学形态学方法、主动轮廓模型、水平集方法或Fast Marching， 均可分割得到二维横断面图像中预处理后感兴趣区域的主颈动脉血管内轮 廓。本实例采用数学形态学方法示例说明。

（3.2.1.2.1）提取感兴趣区域内的边缘，本实例采用Canny算子边缘， 图13-b的边缘提取结果，如图16所示；

（3.2.1.2.2）设置结构基元大小及形状

其中结构基元大小SE_si中， R_maxR_m分别为（3.2.1.1）中计算所得圆的最大半径、最小半径（参见图15）； 且结构基元大小SE_size、形状SE shape可调。如：常用的大小为3、5、7等； 常用的形状为圆形、十字形、方形等。本发明中，结构基元形状SE_shape采 用八角圆盘，大小为3的整数倍数。

（3.2.1.2.3）数学形态学闭运算操作

将结构基元作用于图16上进行数学形态学闭运算操作，获得封闭的内 轮廓所在局域。图16经闭运算操作后的结果如图17所示。

（3.2.1.2.4）提取单像素点血管内轮廓

血管轮廓可用其一系列的轮廓点加以描述。因此，形态学分割血管内 轮廓后，对其所在局域提取单像素点集，即可作为该二维横断面血管的内 轮廓。图17的单像素轮廓点如图18所示。

（3.2.1.3）血管内轮廓的平滑处理（图19、图20和图21）

血管轮廓的平滑处理包括对内、外轮廓的平滑处理。本发明中，采用“面 包卷式”展开，经光顺、多项式拟合的后，再复原获得新的单像素轮廓点集， 具体过程如下：

（3.2.1.3.1）计算轮廓中心点坐标

计算轮廓点集中所有点横、纵坐标的平均值，并记此点为中心点(x₀,y₀)；

（3.2.1.3.2）排序轮廓点

任选一轮廓点，并以该点为起始点，按轮廓逆时针方向依次排序；

（3.2.1.3.3）计算相对距离

依次计算中心点(x₀,y₀)与排序轮廓点(x_i,y_i)的相对距离Dis；

${\mathrm{Dis}}_i=\sqrt{{(x_0-x_i)}^2+{(y_0-y_i)}^2}$

（3.2.1.3.4）展开内轮廓点集（图19）

以排序序号为横坐标，各对应相对距离为纵坐标，获得原轮廓展开点 集；

（3.2.1.3.5）多项式拟合（图20）

对上述点集进行曲线多项式拟合，去除毛刺，使得轮廓展开拟合曲线 平滑。其中，多项式阶数n₀可调，一般取5~10，优选地，n₀＝9；

（3.2.1.3.6）复原轮廓点集（图21）

将平滑后的曲线作为新轮廓展开曲线，并根据序号，反求新轮廓点 (x_j,y_j)坐标，即轮廓点集复原。

$\left(\begin{array}{c}{x}_j=x_0+{\mathrm{Dis}}_i\times \cos \left(\mathrm{angle}\right)\\ y_j=y_0+{\mathrm{Dis}}_i\times \sin \left(\mathrm{angle}\right)\end{array}\right)$

其中：angle为原轮廓中心点(x₀,y₀)与原轮廓点(x_i,y_i)连线的倾角。 最终，内轮廓结果如图22所示，其中实心圆形点表示手动分割轮廓；实线 表示本发明分割轮廓。

（3.2.2）依据各主颈动脉血管的内轮廓，对其进行椭圆拟合，将各拟 合得到的椭圆放大后，作为各主颈动脉血管的初始外轮廓；

（3.2.2.1）依据血管内轮廓，对其进行椭圆拟合；

本发明以步骤（3.2.1.2.4）确定的内轮廓单像素点集为给定点集，椭圆 为拟合函数，拟合获取椭圆的参数E(centre,p,q,θ)：分别表示椭圆圆心坐标、 椭圆长半轴长、椭圆短半轴长、椭圆长轴与水平线夹角。

本发明中，对内轮廓椭圆拟合的结果如图23所示。

（3.2.2.2）血管的初始外轮廓

基于血管形状的先验知识，利用血管内外轮廓的一致性和对称性，对 获得的血管内轮廓进行椭圆拟合后，放大作为血管外轮廓的初始轮廓。

放大（3.2.2.1）拟合的椭圆，作为血管初始外轮廓。优选地，椭圆放大 比例系数factor在1.02~1.08之间。本发明放大的含义为椭圆的长半轴、短半 轴乘以放大比例系数，椭圆的圆心和倾角不变。

（3.2.3）依据各主颈动脉血管的初始外轮廓，采用主动轮廓模型ACM 演化，得到各主颈动脉血管的（最终）外轮廓

基于（3.2.2）血管初始外轮廓，采用主动轮廓模型（ACM）演化，获 得最终血管外轮廓，如图24所示，其中实心菱形点表示手动分割轮廓；实 心方形点表示本发明分割轮廓。

由于本发明血管的外轮廓通过主动轮廓模型获取，内轮廓采用数学形 态学获取，因此外轮廓相对内轮廓具有较好的平滑性，所以可按照上述（3. 2.1.3）方式只对内轮廓作平滑处理。

（3.3）依据各感兴趣区域的内、外轮廓，按其空间位置关系三维重建， 得到三维主颈动脉血管轮廓。该血管轮廓可辅助临床应用。

按步骤（3）依次对所有二维横断面序列图像处理后，获得每一张的血 管内、外轮廓。将标记了主颈动脉血管的各感兴趣区域，依照空间位置关 系顺序堆叠，三维重建得到三维主颈动脉血管；三维重建后，即可获得主 颈动脉CCA的面积、体积等参数，从而为医生临床诊断提供指导意义，常 见的如：药物评价；若引入对比研究，则可用于手术效果评价等等。

步骤（4）依据二维矢状面序列图像，计算矢状面的内中膜厚度。

主颈动脉CCA血管壁，由血管腔（Lumen）开始，从内到外依次由三 层膜组成——内膜（Intima）、中膜（media）以及外膜（adventitia）。超 声正是利用“血管腔（L）-内膜（I）-中膜（M）-外膜（A）”所形成的不 同声学阻抗界面，进行成像从而诊断病情的。

本发明中，提取的内膜轮廓（简称内轮廓）位于血管腔与内膜交界处， 记为LIB（Lumen Intima Boundary）；提取的外膜轮廓（简称外轮廓）位于 中膜与外膜的交界处，记为MAB（Media Adventitia Boundary）；测量的内 中膜厚度（Intima-media Thickness，IMT）即为血管腔-内膜交界面LIB与 中膜-外膜交界面MAB之间的距离。

本步骤的主要技术思路为：提取二维矢状面序列图像的内、外膜，测 量内中膜厚度IMT，计算矢状面内中膜厚度IMT的均值后，提供血管定量 分析测量结果。取其中一张二维矢状面图像（参见图9）作为本例说明，其 它图像均采用同样的方法加以处理，最终获得厚度测量结果；

（4.1）在二维矢状面序列图像（2D_Data_Y）中，分别选取各主颈动 脉血管感兴趣区域（2D_Data_Y_j_ROI），并对各主颈动脉血管感兴趣区域 预处理。预处理的技术思路为：先归一化，再滤波降噪。

选取的二维矢状面序列图像包含主颈动脉CCA、分叉点BF、血管腔、 远端血管壁外膜、远端血管壁外膜外组织，从中提取包含远端血管壁外膜 外组织的子区域作为感兴趣区域ROI，并对其进行归一化和滤波的预处理。 归一化函数可采用线性函数转换、对数函数转换、反余切函数转换等等。 滤波可采用低通滤波、均值滤波、中值滤波、频域滤波、带通滤波等等。 下面以线性函数归一化和低通滤波为例进行说明，分别如图25-（a）和图 25-（b）所示，预处理后的ROI结果如图26所示。

（4.1.1）线性函数归一化（图25-（a））

将二维矢状面图像灰度值进行线性函数转换，使得在[0，1]进行归一化， 表达式如下：

GV_aft=(GV_pre-MinValue)/(MaxValue-MinValue)

其中：GV_pre、GV_aft分别为转换前、后的灰度值，MaxValue、MinValue 分别为图中的最大灰度值和最小灰度值。

（4.1.2）低通滤波（图25-（b））

采用高斯滤波，其具体操作为：用一个均值为0，标准差为10模板， 扫描图像中的每一个像素，计算模板所在邻域内所有像素的加权平均，作 为模板中心像素点的灰度值。一般的，模板大小为5、7、9，本发明中，取 5×5。最后，预处理后的ROI结果如图26所示。

（4.2）在预处理后的各主颈动脉血管感兴趣区域中，分割得到各主颈 动脉血管的内膜和外膜，并计算各感兴趣区域的内中膜厚度；

依据步骤（4）所定义的内中膜厚度，在预处理后的二维矢状面图像的 感兴趣区域中，提取血管的内、外膜，进而测量二维矢状面血管的内中膜 厚度IMT。其血管提取的技术思路为：首先粗分割感兴趣区域ROI，定位 内中膜大致位置；其次，细分割感兴趣区域ROI，分别获得内膜和外膜； 最后计算内中膜厚度。

（4.2.1）感兴趣区域ROI的粗分割（图27）

对ROI区域内的图像上每一列的点，从上到下依次标记0,1,2，……序 号，作为每一列点的索引值。在ROI搜索中，为了减少冗余搜索，有必要 将内、中膜潜在的ROI区域进一步缩小，故对感兴趣区域ROI进行粗分割， 以获得下边界轮廓DB和上边界轮廓UB。感兴趣区域ROI的粗分割结果， 如图27所示。

（4.2.1.1）下边界轮廓DB获取

首先，设定下边界候选点灰度值GV_candidate应满足的阈值条件： GV_candidate≥a×ΔGV+GV_min＝a×(GV_max-GV_min)+GV_min，（一般的，权重系数a的范围为 0.87~0.93，优选地a＝0.9）；然后，从左到右，顺次扫描ROI的每一列，并 获得每一列的最大灰度值、最小灰度值以及满足阈值条件的下边界候选点。 在下边界候选点中，标记索引值最大的点作为唯一下边界轮廓点；最后， 顺次连接ROI中所有下边界轮廓点构成下边界轮廓DB。

（4.2.1.2）上边界轮廓UB获取；

首先，将下边界轮廓DB在ROI中向上平行移动Δ_index个像素距离，Δ_index一般取值为20~30个像素，作为上边界临时轮廓UB_t，优选地Δ_index＝25； 其次，取模板大小X×X，从左到右，从上到下，顺次计算上边界临时轮廓 UB_t和下边界轮廓DB之间的区域中的每个像素所在模板的均值EX、方差 DX。若满足阈值条件：EX≤b且DX≤c，则模板中心为上边界候选点；若不满 足阈值条件，则对应UB_t轮廓点为上边界候选点。（一般的，模板大小X 取8~15；b取0.05~0.09；c取0.11~0.15。优选地，X＝10;b＝0.08;c＝0.14。） 在每一列的上边界候选点中，标记索引值最大的作为唯一上边界轮廓点； 最后，顺次连接ROI中所有上边界轮廓点构成上边界轮廓UB。

（4.2.2）感兴趣区域ROI的细分割（图28）

感兴趣区域ROI的细分割，实质上是在下边界轮廓DB和上边界轮廓 UB之间完成的。粗分割所得的图像区域即为内、中膜所在区域。为了将 DB和UB之间正确的划分出血管腔、内膜、中膜、外膜以及外膜以外组织， 可采用C均值、模糊C均值、支撑向量机SVM、AdaBoost算法等方法。 在本发明中，采用基于模糊C均值聚类方法，加以说明血管壁内外轮廓的 提取。感兴趣区域ROI的细分割结果，如图28所示。

根据远端血管壁特性，对血管部分可直观划分为“黑、灰、白”三类，即 可大致对应“血管腔、内中膜、外膜及以外组织”。在本发明最终的分类结果 中，设定灰度值最低的一类像素点为“血管腔”；灰度值最高的一类像素 点为“外膜及以外组织”；其他灰度值的像素点为“内中膜”。

本发明聚类分割中采用模糊划分，使得每个给定数据点用值在[0，1] 区间的隶属度来确定其属于各个类别的程度。一个数据集的隶属度的和总 等于1。设数据集X＝{x₁,x₂,...,x_n}，每个样本有s个特征属性 x_j＝{x_j1,x_j2,...,x_js}(j＝1,2,...,n)。在本发明中以像素灰度值作为唯一特征属性， 即s＝1。对于ROI中的每一列用FCM聚类来实现分割，这里我们以某一列 为例，设该列共有n个像素样本。具体步骤如下：

（4.2.2.1）参数初始化：给定聚类类别数C（我们取C＝10），设定迭 代停止阈值ε（我们取ε＝1e-5），设定最大迭代次数N（我们取N＝100）；

（4.2.2.2）初始化聚类原型分布矩阵P⁰：随机产生一个大小为n×C的 矩阵，第j行对应该列第j个像素样本的隶属度值的分布，即第j个像素样 本隶属于各类的可能性，且满足

（4.2.2.3）用下列公式计算或更新划分矩阵U^(b)，对于若d_ik≠0则 有：

${\mu }_{\mathrm{ik}}^{\left(b\right)}=1/{\Sigma }_{r=1}^c{[d_{\mathrm{ik}}^{\left(b\right)}/d_{\mathrm{rk}}^{\left(b\right)}]}^{\frac{2}{m-1}}$

若d_ik＝0，则有：且对j≠k，

其中式中d_ik表示第i类中的像素样本x_k与第i类的聚类中心之间的距 离，m为加权指数，b为迭代次数。

（4.2.2.4）用（式1）更新聚类分布矩阵P^(b+1)：

$P_i^{(b+1)}=\frac{{\Sigma }_{k=1}^n{\left({\mu }_{\mathrm{ik}}^{(b+1)}\right)}^m*x_k}{{\Sigma }_{k=1}^n{\left({\mu }_{\mathrm{ik}}^{(b+1)}\right)}^m},i=\mathrm{1,2},...,C$（式1）

（4.2.2.5）若某次迭代前后，聚类分布矩阵满足|P^(b)-P^(b+1)|＜ε或b＝N， 则算法停止并输出划分矩阵U和聚类分布矩阵P，否则令b＝b+1，返回步 骤（4.2.2.3）。最终，获得所有“内中膜”的像素点。

（4.2.3）计算该二维矢状面ROI内中膜厚度

经（4.2.2）聚类分割后，可获得该二维矢状面中ROI每列的内、外膜 对应索引值，计算其差，并计算每一列的平均值，即为该二维矢状面对应 内中膜厚度。

（4.3）依据各感兴趣区域的内中膜厚度，统计计算其均值，即为矢状 面的内中膜厚度

对二维矢状面序列图像的每一张，经上述（4.1）和（4.2）的步骤后， 进行结果统计分析，获得矢状面血管内中膜厚度IMT的均值、方差等参数， 作为血管评价指标。IMT值已成为预测心脑血管疾病的常用有效指标。

步骤（5）依据二维冠状面序列图像，计算冠状面的内中膜厚度：

本步骤中的内中膜厚度定义同步骤（4），即测量的内中膜厚度 （Intima-media Thickness，IMT）即为血管腔-内膜交界面LIB与中膜-外膜 交界面MAB之间的距离。

本步骤的主要技术思路为：提取二维冠状面序列图像的内、外膜，测 量内中膜厚度IMT，计算冠状面内中膜厚度IMT的均值后，提供血管定量 分析测量结果。取其中一张二维冠状面图像（参见图11）作为本例说明， 其它图像均采用同样的方法加以处理，最终获得厚度测量结果；

（5.1）在二维冠状面序列图像（2D_Data_X）中，分别选取各主颈动 脉血管感兴趣区域（2D_Data_X_i_ROI），并对各主颈动脉血管感兴趣区域 预处理。预处理的技术思路为：先增强，再滤波。

选取的二维冠状面序列图像包含主颈动脉CCA、血管腔、远端血管壁 外膜、远端血管壁外膜外组织，从中提取包含远端血管壁外膜外组织的子 区域作为感兴趣区域ROI。图11提取感兴趣区域ROI，结果如图29所示。

对其进行图像增强和滤波预处理。图像增强可采用自适应直方图均衡 化、空间域的规定化、局部统计法等；滤波可采用低通滤波、均值滤波、 中值滤波等等。下面以自适应直方图均衡化和自适应均值滤波为例，进行 说明。获得预处理后的感兴趣区域，如图30所示。

（5.1.1）图像增强

采用自适应直方图均衡化（AHE，Adaptive Histogram Equalization） 进行图像增强。与传统的直方图均衡化相比，AHE更加注重图像的局部特 征。AHE根据像素的局部统计特征来决定对比度增强方法。每个像素的灰 度值都通过一个均衡化变换函数得到，而该变换函数是由以该像素为中心 的一个局部字图像直方图得到的，称其为局部对比度增强法。公式为：

x′_i，j＝m_i，j+k×(x_i，j-m_i,j)

其中，k为自适应参考量，表达式为为窗W内的灰 度方差，为整幅图像的灰度方差，k′为比例系数；x_i，j，x′_i，j分别为变换前后的 灰度值；m_i，j为窗W内像素灰度的平均值。

（5.1.2）自适应均值滤波

均值滤波也称为线性滤波，主要用区域平均法来去除噪声。线性滤波 的基本原理是用均值代替原图像中的各个像素值，即对待处理的当前像素 点的灰度值x_i，j，选择一个模板，该模板由其近邻的若干像素组成，求模板中 所有像素灰度的均值，再把该均值赋予当前像素点，作为处理后图像在该 点上的灰度，即x′_i，j。最后，获得预处理后的感兴趣区域，如图30所示。

（5.2）在预处理后的各主颈动脉血管感兴趣区域中，分割得到各主颈 动脉血管的内膜和外膜，并计算各感兴趣区域的内中膜厚度

依据步骤（4）所定义的内中膜厚度，在预处理后的二维冠状面图像的 感兴趣区域中，提取血管的内、外膜，进而测量二维冠状面血管的内中膜 厚度IMT。其血管提取的技术思路为：首先初始化、并提取血管内膜轮廓； 其次初始化、并提取血管外膜轮廓；最后计算内中膜厚度。

（5.2.1）初始化内轮廓

鉴于主颈动脉CCA的结构特点，二维冠状面序列图像极易受到分叉点 和外颈动脉ECA的伪影影响，导致二维冠状面序列图像至多只有3张左右。 从而增加了二维冠状面分割的难度。为了快速获得接近真实的内膜，本发 明中采用基于数学形态学方法来获取初始化内轮廓。

（5.2.1.1）图像阈值化

阈值化的方法主要有四类：基于点的全局阈值方法、基于区域的全局 阈值方法、局部阈值方法和多阈值方法。本发明实施例采用大津算法OTSU （自适应阈值化），使得白色区域包含内膜、中膜、外膜、外膜外组织以 及部分噪声，黑色区域为血管腔，如图31所示；

（5.2.1.2）填充空洞

为了填充图27所示白色区域中的空洞，可采用模板匹配法、逐点扫描 法等等。鉴于二值图像的特点，本发明采用形态学闭运算，结果如图32所 示。

（5.2.1.3）去除毛刺

为了使得血管内膜表面平滑，符合生理特点，采用数学形态学开运算， 去除白色区域表面附加噪声所产生的毛刺，结果如图33所示；

（5.2.1.3）修正初始内膜轮廓

针对图33所示的内膜表面边缘点，对白色区域的上边缘像素点进行间 隙抽取，以抽取的像素点向下提取与其间距3~5个像素距离的像素点，即 可得到一系列的内膜边缘初始坐标点，将系列内膜边缘初始坐标点连线作 为内膜初始轮廓，如图34所示；

（5.2.2）提取内轮廓（图35）

依据感兴趣区域的梯度信息，采用经典蛇形（Snake）方法计算内膜初 始轮廓的内、外作用力。将内、外作用力进行比较，依据比较结果演化轮 廓直到内、外作用力相等，最终确定的轮廓即为提取的内轮廓，如图35所 示，其中内膜的初始轮廓为虚线；内膜的最终轮廓为实线。除经典蛇形 （Snake）方法以外，还可采用水平集、CV模型、GVF-Snak等演化方法。

（5.2.3）初始化外膜轮廓

将（5.2.1）获得的内膜初始轮廓向下平行移动Δ_index个像素距离，可获得 初始外膜轮廓。一般的Δ_index为15~30，优选地Δ_index取25。

（5.2.4）提取外膜轮廓（图36）

依据感兴趣区域的梯度信息、场信息，采用梯度向量场方法GVF-Snake 计算外膜初始轮廓的内、外作用力。将内、外作用力进行比较，依据比较 结果演化轮廓直到内、外作用力相等，最终确定的轮廓即为外膜轮廓，如 图36所示，其中外膜的初始轮廓为虚线；外膜的最终轮廓为实线。除 GVF-Snake算法外，还可采用经典蛇形（Snake）、水平集、CV模型等演 化方法。

（5.2.5）计算该二维冠状面ROI内中膜厚度（图37）

经（5.2.1）-（5.2.4）分割后，可获得该二维冠状面总ROI每列内、 外膜对应像素点位置，计算内、外膜轮廓对应像素点之间的像素点间距， 并将其转化为实际距离，并计算每一列的平均值，即为该二维冠状面对应 内中膜厚度。如图37所示，上方内膜最终轮廓（实线）和下方外膜最终轮 廓（实线）的间距即为所求内中膜厚度。

（5.3）依据各感兴趣区域的内中膜厚度，统计计算其均值，即为冠状 面的内中膜厚度；

对二维冠状面序列图像的每一张，经上述（5.1）和（5.2）的步骤后， 进行结果统计分析，获得冠状面血管内中膜厚度IMT的均值、方差等参数， 作为血管评价指标。IMT值已成为预测心脑血管疾病的常用有效指标。

方案验证：

将（3）（4）（5）的分割结果与手动分割结果比较，计算三视图分割 结果的DSC(Dice Similarity Coefficient)、MAD(Mean Absolute Distance)和 MAXD(Maximum Absolute Distance)，进行算法验证性评价与分析。如图 22、图24、图38和表1所示：

A、相似系数DSC：

相似系数DSC是一种度量分割算法的标准，用于评价面积的相似度。

$\mathrm{DSC}=2\frac{|R_M\cap R_A|}{\left|R_M\right|+\left|R_A\right|}$

其中，R_M和R_A分别表示手动轮廓围绕的区域和自动轮廓围绕的区域。 相似系数DSC越接近1，说明两方法的接近程度越高。

B、平均绝对距离MAD和最大绝对距离MAXD

平均绝对距离MAD和最大绝对距离MAXD，都是用来度量分割算法 的标准，用于评价距离相似度。

假设手动轮廓点为{m_i：i＝1，廓点为准，对应自动轮廓点为{a_i：i＝1，轮廓点为。 则

$\mathrm{MAD}=\frac{1}{k}\underset{i=1}{\overset{k}{\Sigma }}\left|d(m_i,A)\right|$

$\mathrm{MAXD}=\underset{i\in \{1,k\}}{\max }\left\{\right|d(m_i,A)\left|\right\}$

其中，d(m_i,A)为手动轮廓与自动轮廓的对应点间的相对距离。平均绝对 距离MAD和最大绝对距离MAXD越小，说明两方法的接近程度越高。

表1本发明与手动分割结果比较

  DSC（%）   MAD(mm)   MAXD(mm)   二维横断面内轮廓提取血管   94.1±3.6   0.32±0.21   0.77±0.47   二维横断面外轮廓提取血管   92.8±2.5   0.37±0.19   0.84±0.39   二维矢状面IMT测量厚度   0.73±0.47   1.01±0.75   二维冠状面IMT测量厚度   1.02±0.76   1.89±0.34

上表中，测量厚度均值、标准差，单位为mm。

从表中可以看出，（1）本发明中所采用的分割方法，与专家手动分割 方法相当，验证了其可行性与鲁棒性；（2）在二维横断面的序列图像分割 中，本发明方法与专家手动分割方法有较高的相似度，相似系数DSC均大 于90%，平均绝对距离MAD和最大绝对距离MAXD均在误差范围内；（3） 在二维矢状面和冠状面的IMT测量中，两方法均与手动测量IMT金标准相 当，平均误差均在1毫米左右；（4）二维矢状面的测量结果优于二维冠状 面，这与序列图像的数量、质量、具体分割方法均有一定关系。另外，从 操作时间上比较，本发明在3分钟即可完成一个三维体数据的测量；而经 验丰富的医生则往往需要8分钟完成手动测量。

以上所列实施示例仅表达了本发明的几种实施方式，其描述较为具体 和详细，但并不能因此而理解为对发明专利范围的限制。应当指出并强调 的是，对本领域的相关技术人员，在不脱离本发明的前提下，还能有多种 变形和改进，这些都应属于本发明的保护范围。因此，本发明专利的保护 范围应以所附权利要求为准。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已， 并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等 同替换和改进等，均应包含在本发明的保护范围之内。