改进型骨骼细化算法提取冠状动脉中心线

张子恒，祝 磊，马 骏，徐 平，刘亦安，薛凌云

(杭州电子科技大学自动化学院，浙江 杭州 310018)

冠心病(coronary artery disease, CAD)在中国心血管疾病死因中居第2位[1]。诊断CAD主要依靠CT血管造影术(computed tomography angiography, CTA)和冠状动脉造影(coronary artery angiography, CAG)。CTA易于三维可视化，扫描时间短、空间分辨率较高，已越来越多用于诊断CAD。CTA诊断CAD时，需要对图像的血管结构进行包括分割血管区域，提取中心线、分叉点以及测量血管直径等精确处理，其中提取血管中心线最为重要[2]。

现有中心线提取方法分为全自动、半自动和人工提取3类。人工提取中心线简单、可靠，但耗时、费力[3]。全自动方法一般指提取中心线过程中不施加任何人工干预，较经典者如利用圆柱体模型拟合提取冠状动脉血管后进一步提取中心线[4]；半自动方法指在血管提取过程中需人为指定1个或数个点作为参照点，如分割主动脉和冠状动脉并行三维重建后利用血管局部灰度值及走向选取迭代起始点[5]。此外，还有利用DBSCAN聚类算法进行分割后利用冠状动脉主曲线提取中心线[6]，利用改进型Frangi滤波器消除心腔边界处无关阶跃响应后自动提取中心线[7]，利用快速行进算法提取图像中初始点之间的最短路径，据此在原始图像上计算冠状动脉中心线[8]，以及同时使用基于边界的距离变换和基于源点(种子点)的距离变换[9]等。随着人工智能的兴起，可利用神经网络提取中心线，如手动设置种子点，利用卷积神经网络(convolutional neural network, CNN)结合冠状动脉半径沿血管不同方向预测血管中心线走向而提取中心线[10]。针对传统骨骼细化方法在血管连接处易产生细小分支的问题，本研究结合Dijkstra算法提出一种基于骨骼细化的全自动中心线提取算法，旨在使所提取的骨架线更加准确、平滑。

1 提取冠状动脉中心线

1.1 图像预处理 原图中分割出冠状动脉进行三维重构，之后去除噪声和增强血管区域。

各向异性扩散滤波器(anisotropic diffusion filter)[11]系空域滤波器，平滑程度取决于扩散函数，在图像各个方向上有所不同。各向异性扩散滤波见公式1：

(1)

(2)

其中，k为与噪声级和边界强度相关的常数，这里选取k=max|u|。

若三维图像I中存在一点A，想要判断点A是否为血管上一点就需要分析图像I在A点的局部特征，利用Hessian矩阵就可以判断该点的位置。点A的Hessian矩阵H(A)见公式3：

(3)

式中Ixx(A)，Ixy(A)等分别代表图像I在点A处的二阶偏导，并且H(A)三个特征值符合|λ1|≤|λ2|≤|λ3|的关系。若点A属于血管上一点利用基于Hessian矩阵的Frangi血管增强函数[12]就可将血管区域增强，该血管增强函数见公式4：

(4)

1.2 提取中心线 三维细化算法属三维形态学操作，可用“在草地上放火”加以解释。将图像视作一块草地，在其两端各位置同时点火，草地上会形成两面火墙，逐渐向草地中央逼近，来自草地两端火墙相遇位置的连线就是该草地的骨架线，各点火点即为算法的起始点，被烧掉区域即为算法中可删除的点。

三维空间中存在一点p，根据该点到其相邻一点的欧氏距离的长度，可将其关系分为6-邻接、18-邻接和26-邻接，见图1。用Nj(p)(j=6,18,26)表示点p和j邻接点的集合即邻域(neighbor)。图1中用U、D、W、E、S、N和点p表示N6(p)，6-邻接中的6个方向即为三维立体网格中的6个主要方向；用N6(p)和12个圆圈标识表示N18(p)；用N18(p)和6个五角星标识表示N26(p)。三维细化算法通常在二值化图片上进行操作，而二值化图像只有0，1两种值，将图中值为1的点定义为目标点(black point)，值为0的点定义为背景点(white point)。如果集合N26(p)/{p}恰好包含一个目标点，则目标点p是图片中的曲线端点(curve end-point)；若N6(p)包含至少1对相对背景点，则目标点p是图片中的表面端点(surface end-point)。如果图像中任意一目标点至少与1个背景点6-邻接，则称为边界点(border point)。简化点(simple point)则是可删除(将该点的值由1置0)而不改变图形拓扑结构的目标点，应只由某些类型边界点构成；曲线或曲面的端点与图形的拓扑结构相关，因而不能被删除。任何细化算法中，可删除的点只有简化点，而传统三维细化算法大多是从体素的6个主要方向并行操作，可能会将一些曲线或曲面的端点当作简化点删除，导致原图像的拓扑结构发生改变。为保证细化算法不会误删某些端点，以保证原图形的拓扑结构，本研究改进传统算法，使每个主迭代步骤由多个并行子迭代过程组成，其中每个子迭代过程只能删除某种特定的边界点。

图1 6-邻接、18-邻接和26-邻接示意图

采用12方向对三维重构后的血管进行细化处理，每个迭代步骤根据图2所示12个方向进行并行细化。选择删除方向的顺序为(US，NE，WD)、(ES，UW，ND)、(SW，UN，ED)、(NW，UE，SD)，据此又可将这12个方向分为4组，每组3个方向，每个小组包含6个主要方向，即在第1个子迭代过程US方向上可删除某些U和S方向上的边界点，其余各迭代过程以此类推。需要被细化的冠状动脉会在每个方向上均匀收缩，即从外层开始，向内被层层腐蚀。当12个子迭代过程中不再有可被删除的目标点时，主迭代过程结束，即中心线提取完成。

图2 12方向示意图

迭代过程中，利用匹配模板可判断某点是否为可被删除的简单点。若图像上一点p的3×3×3邻域与一组匹配模板中的1个模板相符，则认为点p是可在细化过程中删除的简单点。以图2中US方向为例，展示该方向上的14个基本模板(图3)，可用这些模板经过一系列旋转或反转操作获得其他方向的模板。

1.3 剔除骨架线分支 血管为不规则管状结构，进行细化算法时会产生细小分支，影响后续测量冠状动脉直径及判定是否产生狭窄的准确性。为消除干扰，提出基于Dijkstra算法的分支消除算法。

Dijkstra算法是一种解决有向加权图中单源最短路问题的算法，特点是迭代从起点向外层层扩张，直至到达终点[13]。首先将上一步提取出的三维骨架线转换为由节点和有向边组成的有向图，图中每个节点代表被提取出的体素(voxel)，不同体素间欧氏的距离由有向边表示。根据Dijkstra算法，将所有体素设为全集U；用点集G表示已找到最短路径的点，初始状态下仅有起点S∈G；设点集O表示还未找到最短路径的点的集合，则可以推出O=U-G；设Lk为当前情况下，起点经过G中若干点到点k(k∈U)的最短距离，初始状态下LS=0，其他均为+∞。当算法开始时，从起点S开始，沿某条有向边(设权值为arcs)找到起点的一个邻居n；此时令Ln=min{Ln,LS+arcs}；按此方式更新起点所有邻居；在集合O中找到Lk最小的点v，则Lv即起点到v的最短路径长度；将点v从O中取出加入G，对点v重复上述所有操作；如此重复，直到G=U，即O=φ时，算法结束，Lk即为从起点到各点的最短路径长度。冠状动脉中心线的有向图模型见图4。最后再将有向图转换为三维骨架线，即可得到去除细小分支的冠状动脉中心线。

2 结果与评估

2.1 运行环境与实验数据 算法在Matlab R2019a平台下编译、运行。电脑配置为64位Windows 10操作系统，CPU主频2.20 GHz，内存16 GB。所用数据均来自鹿特丹冠状动脉算法评价框架(Rotterdam Coronary Artery Algorithm Evaluation Framework, RCAAEF)提供的公共数据集，包含32例，于2005年6月—2006年6月采自荷兰鹿特丹伊拉斯谟医学中心(Erasmus MC)，其中20例以Siemens 64排螺旋CT、12例以Siemens双源CT采集。

2.2 评价标准 用2个指标评价中心线提取的准确率[14]。第1个指标为重叠率(rate of overlap, OV)，计算见公式5：

图3 US方向上的14个基本模板 点p为图中各模板中心点，图中符号“○”表示背景点，“●”表示目标点，“·”表示该位置上可为目标点也可为背景点；字母X表示在这些位置上至少有1个为目标点，V和W表示这些位置上至少有1个为背景点；如2个不同位置标注字母Z，则表示这2个位置上必须为不同的点(如1个为目标点，另1个为背景点)，但其位置不固定

图4 冠状动脉有向图模型 虚线部分为提取的骨架线模型，点S为中心线的起点，点E为终点，A为骨架细化后产生的细小分支的端点。虚线SE表示实际的中心线走向，和表示三点的有向边其值等于两点之间的欧氏距离，所以图4中的最短路径即为间的体素点得以保留，A点则被删除

(5)

式中‖·‖表示为这一类点的集合，TPMov表示本文算法所求出的中心线上的一点，该点到实际中心线的位置小于参考标准中给出的误差半径，与之对应的实际中心线上的点记为TPRov；若求出的点位于误差半径之外，则将其记为FPov，与之对应的实际中心线上的点记为FNov，这4种点的位置关系见图5。

第2个指标为平均距离(average distance， AD)，即所有由本文算法提取出的中心线的位置与实际中心线之间的距离，用distance(i)表示第i个点到该点的实际位置的距离，则AD见公式6：

(6)

图5 评价中所用术语的位置说明 血管分为A、B、C、D四段，代表中心线实际位置与算法提取的中心线之间的不同结果，位于图片顶端的字符表示该段中算法提取的中心线类型，中心线下方的字符代表实际中心线类型。TPMov、TPRov、FNov均只存在一种情形，而FPov有可能存在于未出现冠状动脉的区域(D段)

图6 图片预处理过程 A.原始CTA图像； B.各向异性扩散滤波处理； C.血管区域； D.三维重构血管

冠状动脉为直径逐渐缩小的管状结构，近心段冠状动脉直径较大，而其分支血管远端直径较小，临床意义较小，故可以血管中首次出错前的重叠率(overlap until first error, OF)来评估程序算法的稳定性。OF见公式7：

(7)

如图5所示，TPRof存在于A段中，C段之前因为已经发生错误，故该段无TPRof。

2.3 实验结果 本研究冠状动脉CTA图像为512×512，每例包含近300幅图像，根据不同长度，每条血管会在40～130幅图像中连续出现，图6A为原始图像。利用各向异性扩散滤波器对原始图像进行去噪处理(图6B)，可明显抑制噪声；根据数据集中提供的每条冠状动脉的中心点位置，将无关区域删除(像素值置0，不改变图像大小)，保留原图信息的区域包含该血管在此层切面的完整图像；再利用Frangi血管增强函数，增强被保留图像中的血管区域，使其与背景形成较明显的差异，为区域生长法提供更易区分的灰度值；最后在该区域中选取血管区域的像素作为种子点进行区域生长，获得完整血管区域(图6C)。对每例所有存在该条血管的层切面进行上述处理，可在二维图像中分割出整条血管；再根据原始图像中像素与世界系坐标的对应关系进行三维重建，即得到该血管的三维数据(图6D)。

完成图片预处理后，利用改良骨骼细化算法提取血管中心线，删除细小分支(图7)。

为保证算法具有较强的鲁棒性，分别选取6个不同病例的血管进行评价(表1、2)。

表1 未改进算法提取中心线部分结果

图7 提取中心线结果 A.右冠状动脉； B.左前降支浅 (灰色部分为重构的冠状动脉，中间蓝色曲线为提取中心线)

表2 改进算法提取冠状动脉中心线部分结果

相比未移除细小分支的骨骼细化算法，本算法将OV提高了2%，AD减少了38.2%，OF提高了16%；中心线提取的平均时间为0.48 s。计算数据集中所有血管的OV、AD、OF，平均重叠率、平均距离和平均第1次出错前重叠率分别为0.985、1.974和0.832。

3 讨论

实验结果表明本算法用于提取冠状动脉中心线有效，相比文献[4-5]提出的重叠率分别为0.847和0.670的2种方法，本算法精度较高。提取冠状动脉中心线在计算机辅助诊断CAD中具有重要作用，提取速度直接影响运行效率。本算法提取中心线用时较短，可满足临床需求。

本算法仍存在需改进之处：未对提取完整冠状动脉血管树中心线性能进行测试，后续将扩大数据样本进一步完善；以区域生长法分割血管区域，要求较高，且较依赖图像阈值，即图像中因加入对比剂而产生的高亮区域，若远端血管无强化，则分割将会产生较大误差，从而影响提取冠状动脉中心线。后续研究中拟与人工智能相结合，利用训练完成的神经网络自动分割血管区域。