一种基于CTA图像的冠状动脉分割方法与流程

本发明属于医学图像处理领域，具体涉及一种非诊断目的的基于cta图像的冠状动脉分割方法，是一种在冠状动脉ct血管造影(computedtomographyangiography，cta)图像上左右冠状动脉支起点识别以及冠状动脉分割的方法。

背景技术：

心脏是一个通过循环系统将血液泵送到全身的器官，它保证了全身各组织可通过血液获得足够的氧气和营养，并清除二氧化碳和其他废物，而冠状动脉则是包裹在心脏外周为心脏自身提供补给的血管，对心脏的正常工作起到关键性的作用。整个冠状动脉树拓扑结构复杂，分叉多且细小，其直径范围一般在2mm到7mm之间。从升主动脉根部出发，冠状动脉可以分为左主冠状动脉(leftmaincoronaryartery，lmca)和右冠状动脉(rightcoronaryartery，rca)两个主要分支，lmca主要负责向左心房和左心室供应血液，其向下延伸还可再细分为左旋支(leftcircumflexartery，lcx)和左前降支(leftanteriordescendingartery，lad)，rca主要负责向右心房、右心室、窦房结以及房室结提供血液，其继续向下延伸还可细分为右后降支(posteriordescendingartery，pda)和右边缘支(rightmarginalartery，ma)。

cta图像是一种通过静脉注射造影剂后，利用ct设备对患者心脏进行断层扫描形成的一个二维图像序列。随着成像分辨率的不断提升以及技术本身无创性和便捷性的优点，cta逐渐成为临床上早期筛查冠状动脉血管疾病的有效方式。基于cta的冠状动脉精准分割是后续冠脉狭窄定量分析和斑块分类等多种病灶识别的重要前置步骤。目前，心脏病专家的手动分割结果是最准确的，但手动分割非常耗时，并且不同专家的人为主观因素也会导致分割结果存在细微差异。近年来，各种计算机辅助诊断技术(computeraideddiagnosis，cad)被用于冠状动脉分割，但是大部分仍属于半自动分割方法，需要手动设定左右冠状动脉各段分支起始种子点以及相应的阈值才能完成后续的冠状动脉分割，且容易遗漏冠状动脉的细小分支。少部分全自动方法很容易产生过分割，将非冠状动脉的血管也误分割为冠状动脉，从而降低了分割的准确度。

技术实现要素：

本发明的目的在于提供一种基于cta图像的全自动冠状动脉分割方法，能够完整、准确地分割并提取冠状动脉树，辅助医生提高心血管疾病诊断效率。

一种基于cta图像的冠状动脉分割方法，其特征在于包括如下步骤：

步骤一、获取原始cta心脏图像；

步骤二、分割冠状动脉感兴趣区域，抑制肺部血管等非心脏组织并去除大部分与冠状动脉灰度值相近的心房与心室等心脏组织，减少非冠状动脉组织对后续冠脉分割产生的影响：

2-1对cta原始图像进行窗宽窗位调整，去除肺部血管；

获取dicom头文件中的窗宽(windowwidth，ww)和窗位(windowlevel，wl)信息，截取ct值在范围内的原始cta图像并归一化到0～255的灰度级，ct值小于的像素点设置为0，ct值大于的像素点设置为255；映射过程如公式(1)所示：

式(1)中y表示归一化后每个像素点的灰度值,x表示原始cta图像中每个像素点的ct值；

2-2对窗宽窗位调整后的cta图像进行基于改进灰狼优化算法的多级阈值处理，用以去除与冠状动脉灰度值相近的心房与心室等组织，获得冠状动脉感兴趣区域；

2-2-1设定最大迭代次数k；随机初始化n组可行解向量每组可行解向量包含dim个取值范围在0到255之间的灰度阈值，可将0到255的图像灰度范围分为dim+1个阈值区间c0,c1…cdim，其中上述每一组可行解向量都为目标函数可行解；

2-2-2将上述n组可行解向量分别代入到目标函数，将目标函数值按照大小排序，取前三大目函数标值所对应的可行解向量分别赋予最优解次优解第三优解

多级阈值的目标函数：

式(2)中hm是某可行解向量内灰度阈值所构建的第m个阈值区间内的熵函数；

式(3)中h0,h1,…hm为每个阈值区间内的kapur熵，t1,t2,…tm为每组可行解向量中的灰度阈值，p0，p1…pi为阈值区间内每一个灰度值在图像中出现的概率，ω0,ω1…ωm为每一个阈值区间的概率，l＝256；

2-2-3根据公式(4-6)更新参数向量和控制系数a：

a＝2(1-t/k)(6)

其中a根据迭代次数t从2线性递减至0，k为总的迭代次数，是[0,1]内的随机向量,n＝1,2,3；

2-2-4根据当前迭代次数中三个等级依次递减的最优解向量通过公式(7-10)更新其余解向量

式(7)中代表每一组解向量和三个最优解向量之间的距离，由公式(5)得到；

式(8)中是每一组可行解向量根据三个最优解解向量计算得到的候选可行解向量，由公式(4)得到；

根据公式(2)分别计算的目标函数值j1、j2、j3并根据公式(9)计算权重系数w1、w2、w3；

最后利用公式(10)根据权重系数w1、w2、w3更新每一组可行解向量得到

2-2-5将更新后的n组可行解向量分别代入到目标函数，将目标函数值按照大小排序，取前三大目函数标值所对应的可行解向量分别赋予最优解次优解第三优解判断当前迭代次数是否满足最大迭代次数，若满足则结束，输出最优解反之则返回2-2-3继续迭代；

2-2-6根据最优解对窗宽窗位调整后的cta图像保留最高一级灰度阈值区间cdim内的像素，最终得到冠状动脉感兴趣区域；

步骤三、基于多尺度融合的血管特征提取，增强不同直径大小的冠状动脉；首先对cta数据内所有体素点计算不同尺度σ下的hessian矩阵h(a)，如式(11)所示；

式(10)中每一个子项iij(a)表示图像i在点a处的二阶偏导，矩阵h(a)存在三个特征值λ1、λ2、λ3且满足|λ1|≤|λ2|≤|λ3|的关系，以及三个特征值各自对应的特征向量为frangi血管增强函数v(σ,a)可以利用多尺度hessian矩阵计算所得的三个特征值λ1、λ2、λ3对血管状区域进行增强，对非血管区域进行抑制，血管相似性函数v(σ,a)如式(12)所示：

式(11)中，ra、rb用于量度管状、盘状和斑状等不同结构之间的特征差异，s是为对比度系数，用于描述目标和周围背景的亮度差异，s的值就越大，亮度差异越大；frangi血管增强函数v(σ,a)的取值范围在0到1之间，v(σ,a)越接近于1表示该点处的组织与管状结构的相似性最大；α、β和γ为三个调参系数，通过影响ra、rb和s在血管相似性函数中所占的比重来调整式(12)的灵敏度；ra、rb和s定义分别如式(13)、式(14)和式(15)表示：

步骤四、升主动脉分割以及左右冠状动脉起始层识别；结合光流法检测相邻两层升主动脉间光流位移距离最大值的变化和心脏医学解剖结构知识判断左右冠状动脉的起始层；

4-1由于冠状动脉的根节点位于升主动脉根部，升主动脉的分割和提取是识别左右冠状动脉起始层的前提；采用霍夫变换检测类圆形的升主动脉区域，以被检测到的类圆形圆心作为种子点，利用区域生长法分割每一层升主动脉，最后得到完整的升主动脉；

所述的霍夫变换检测类圆形的升主动脉区域属于常规技术，故不详解；

4-2使用lucas-kanade光流法检测相邻两层升主动脉间的最大光流及所有像素中的最大位移距离；

在cta图像序列中，大部份的升主动脉区域都呈现类圆形，但是当某一升主动脉层包含冠状动脉分支时，类圆形升主动脉区域会多一小段凸起的冠状动脉区域，该层相对于前一层在冠状动脉根结点处的光流位移距离会存在巨大的变化，通过对比相邻2层间光流位移距离最大值的变化即可初步筛选出冠状动脉的起始层；光流法假设相邻的2层cta图像强度为i(x,y,t)和i(x+δx,y+δy,t+δt)，并且相邻的2层图像强度不变，可以得到公式(16)

i(x+δx,y+δy,t+δt)＝i(x,y,t)(16)

将公式(16)使用泰勒级数展开后可得：

其中为高次项，大小可以忽略，为此结合公式(16)和(17)可得：

其中和分别为图像相对于x轴和y轴的空间导数，为时间导数，公式(18)可改写为

(ix,iy)(u,v)＝-it(19)

其中，u和v就是图像中每个像素分别沿着x轴和y轴的光流位移距离，但是只有公式(19)一个约束方程无法求解2个未知变量，所以光流法假设目标像素的周围像素也存在相同的光流位移距离，即在目标像素周围m*m的窗口内，每个像素点都具有相似的光流位移距离u和v；简单的来说，就是通过增加约束方程来求解光流矢量u和v，如公式(20)所示；

式(19)中(xn，yn)表示m*m的窗口内第n个像素点；

通过公式(20)求解所得的u和v可以计算出每个像素点的位移距离，如公式(21)所示；

根据每一层distances最大的点，筛选出最大的前10层做为候选冠脉起始层，然后结合三个心脏医学解剖结构先验知识进一步筛选出左右冠状动脉的起始层；

所述的三个心脏医学解剖结构先验知识具体如下：

第一，升主动脉根部只存在左右两个冠状动脉开口，在10层冠脉候选层中，位移距离最大的是左冠状动脉起始层sleft；

第二，左右冠脉起始层间相隔10层以上，且左冠状动脉起始层位于右冠状动脉起始层的z轴上方；

第三，左冠状动脉起始点ol位于该层质心的-90°～45°之间，右冠状动脉起始点or位于该层质心的位于45°～135°之间；

最后在满足上述第二、三条件的剩余候选冠脉起始层中选择层号最小的层作为右侧冠状动脉起始层sright；

步骤五、在步骤三增强后的cta图像上，从步骤四获得的左右冠状动脉起始层开始，使用结合端点检测的区域生长法分割每一层cta图像中的血管区域；具体是从步骤三获得的左右冠状动脉起始层sleft和sright开始，沿着z轴向下或向上利用迭代区域生长法进行冠状动脉的初始分割；

所述的迭代区域生长法由步骤5-1种子点自动选择、步骤5-2迭代区域生长法生长准则、步骤5-3端点检测以及步骤5-4端点再生长四部分组成：

5-1种子点自动选择

利用血管上下层之间的空间连通性，即当前层si与相邻层si+1血管区域内的二维坐标存在部分重叠作为每层种子点的选择依据；由于三维的冠状动脉血管分叉众多，所以在二维的单层冠状动脉cta图像中可能存在多个独立的血管连通区域，当si层分割完成后，先检测血管连通区域个数c，然后分别在si层的每个独立的血管区域中以一定的步长选取m个种子点作为相邻层si+1血管区域生长的种子点；m为经验值；

5-2迭代区域生长法生长准则

利用迭代区域生长法分割每一层的冠状动脉血管区域；

5-2-1首先从5-1中选取的每一层c*m个种子点依次开始对该层血管区域进行分割，每一个种子点会将区域边缘像素点8领域内与种子点灰度值相差在1％以内的像素点不断归并到感兴趣区域中，直到不能生长为止；

5-2-2由于5-1中选取的种子点只有部分会落在血管区域内，所以每个种子点生长完成后，通过计算分割得到的血管感兴趣区域内的像素数量num判断该种子点是否在血管区域内：如果像素数num小于nthresh，则表明该种子点位于血管区域内，保留该种子点分割得到的血管感兴趣区域；如果像素数num大于nthresh，则说明该种子点落在血管区域之外，出现了严重的过分割现象，丢弃该种子点分割得到的区域，nthresh为经验值；然后将所有保留下来的血管感兴趣区域之和作为该层最终的冠状动脉血管分割结果；

5-2-3判断每一层最终的冠状动脉血管分割结果像素数是否为0，如果等于0，则表明无法再继续分割血管，结束迭代，得到初始冠状动脉感兴趣区域；如果不等于0，则返回5-1继续执行下一层种子点的选择以及5-2迭代区域生长法；

5-3端点检测

对5-2中提取的初始冠状动脉感兴趣区域利用细化方法提取血管骨架线；具体是利用血管骨架线的26邻域单连通性，将骨架线上的体素分为端点p1、普通连接点p2、分叉点p3；从各个端点p1和分叉点p3开始，继续使用结合5-1的种子点选择机制和5-2的迭代区域生长法尝试z轴向上和向下的双向生长，用来分割冠状动脉分支中拓扑结构复杂，上下起伏的血管部分；

所述端点p1在该体素的26邻域点中，只有1个邻域点位于骨架线上；所述普通连接点p2，在该体素的26邻域点中，只有2个邻域点属于骨架线上；所述分叉点p3，在该体素的26邻域点中，至少有3个邻域点位于骨架线上；

5-4端点再生长

判断结合端点检测的区域生长法继续分割后的冠状动脉区域是否有新的端点p1和分叉点p3，如果有，则返回5-4继续分割；如果没有，则结束分割。

步骤二中采用改进的灰狼优化算法，式(8)中采用离散的权重w1、w2、w3对每一次迭代中各组可行解向量进行更新，相比于传统灰狼优化算法采用的平均权重有助于提高最优解的搜寻效率。

步骤三中采用基于多尺度融合的血管特征提取，增强不同直径大小的冠状动脉。由于整个冠状动脉树从升主动脉根部出发会变得越来越细，其直径范围一般在2mm到7mm之间。单一尺度σ下的hessian矩阵h(a)特征值只能用于增强特定直径的血管，为此本专利计算多个不同尺度σ下的hessian矩阵h(a)的特征值λ1、λ2、λ3用于提取不同直径血管的特征。最后，frangi血管增强函数v(σ,a)可利用不同尺度σ下hessian矩阵的特征值对不同直径的血管状区域进行增强，对非血管区域进行抑制。

本发明的有益效果如下：

本发明提出了一种全自动的冠状动脉分割方法，克服了半自动方法中需要手动初始化冠状动脉起始种子点的技术难点。本发明首先利用窗宽窗位调整和基于灰狼优化的多级阈值处理，有效去除原始cta图像中非冠状动脉组织可能对后续冠状动脉分割产生的干扰；然后利用升主动脉和冠状动脉的空间位置关系，采用光流法识别冠状动脉的起始层；最后从左右冠状动脉的起始层开始，采用结合端点检测的区域生长法对冠状动脉进行分割，相比于传统区域生长法，对于拓扑结构复杂的细小分支有更好的分割能力和准确性。

附图说明

图1a为原始cta图像；

图1b为原始cta图像进行窗宽窗位调整和基于改进灰狼算法优化的多级阈值处理后的结果；

图2为基于多尺度融合的血管增强结果；

图3a为分割得到的升主动脉三维重构结果；

图3b为光流法检测得到的每层最大光流位移距离结果；

图3c为左右冠状动脉起始层中，具有最大光流位移距离的像素坐标和该层质心间的位置关系图；

图3d-e分别为光流法检测得到的左右冠状动脉起始层；

图4a为使用的冠状动脉骨架线示意图；

图4b为传统区域生长法分割的冠状动脉结果；

图4c为结合端点检测的区域生长法分割的冠状动脉结果；

图5为本发明方法流程图；

图6为基于改进灰狼优化算法的多级阈值处理流程图；

图7为基于结合端点检测的区域生长法冠状动脉分割流程图；

具体实施方式

以下结合附图对本发明作进一步说明。

如图5一种基于cta图像的冠状动脉分割方法，包括如下步骤：

步骤一、获取原始cta心脏图像；

步骤二、分割冠状动脉感兴趣区域：

2-1对cta原始图像进行窗宽窗位调整，去除肺部血管；

获取dicom头文件中的窗宽(windowwidth，ww)和窗位(windowlevel，wl)信息，截取ct值在范围内的原始cta图像并归一化到0～255的灰度级，ct值小于的像素点设置为0，ct值大于的像素点设置为255；映射过程如公式(1)所示：

式(1)中y表示归一化后每个像素点的灰度值,x表示原始cta图像中每个像素点的ct值；

2-2对窗宽窗位调整后的cta图像进行基于改进灰狼优化算法的多级阈值处理，用以去除与冠状动脉灰度值相近的心房与心室等组织，获得冠状动脉感兴趣区域如图6，具体是：

2-2-1设定最大迭代次数k；随机初始化n组可行解向量每组可行解向量包含dim个取值范围在0到255之间的灰度阈值，可将0到255的图像灰度范围分为dim+1个阈值区间c0,c1…cdim，其中上述每一组可行解向量都为目标函数可行解；

2-2-2将上述n组可行解向量分别代入到目标函数，将目标函数值按照大小排序，取前三大目函数标值所对应的可行解向量分别赋予最优解次优解第三优解

多级阈值的目标函数：

式(2)中hm是某可行解向量内灰度阈值所构建的第m个阈值区间内的熵函数；

式(3)中h0,h1,…hm为每个阈值区间内的kapur熵，t1,t2,…tm为每组可行解向量中的灰度阈值，p0，p1…pi为阈值区间内每一个灰度值在图像中出现的概率，ω0,ω1…ωm为每一个阈值区间的概率，l＝256；

2-2-3根据公式(4-6)更新参数向量和控制系数a：

a＝2(1-t/k)(6)

其中a根据迭代次数t从2线性递减至0，k为总的迭代次数，是[0,1]内的随机向量,n＝1,2,3；

2-2-4根据当前迭代次数中三个等级依次递减的最优解向量通过公式(7-10)更新其余解向量

式(7)中代表每一组解向量和三个最优解向量之间的距离，由公式(5)得到；

式(8)中是每一组可行解向量根据三个最优解解向量计算得到的候选可行解向量，由公式(4)得到；

根据公式(2)分别计算的目标函数值j1、j2、j3并根据公式(9)计算权重系数w1、w2、w3；

最后利用公式(10)根据权重系数w1、w2、w3更新每一组可行解向量得到

2-2-5将更新后的n组可行解向量分别代入到目标函数，将目标函数值按照大小排序，取前三大目函数标值所对应的可行解向量分别赋予最优解次优解第三优解判断当前迭代次数是否满足最大迭代次数，若满足则结束，输出最优解反之则返回2-2-3继续迭代；

2-2-6根据最优解对窗宽窗位调整后的cta图像保留最高一级灰度阈值区间cdim内的像素，最终得到冠状动脉感兴趣区域；

步骤三、基于多尺度融合的血管特征提取，增强不同直径大小的冠状动脉；首先对cta数据内所有体素点计算不同尺度σ下的hessian矩阵h(a)，如式(11)所示；

式(10)中每一个子项iij(a)表示图像i在点a处的二阶偏导，矩阵h(a)存在三个特征值λ1、λ2、λ3且满足|λ1|≤|λ2|≤|λ3|的关系，以及三个特征值各自对应的特征向量为frangi血管增强函数v(σ,a)可以利用多尺度hessian矩阵计算所得的三个特征值λ1、λ2、λ3对血管状区域进行增强，对非血管区域进行抑制，血管相似性函数v(σ,a)如式(12)所示：

式(11)中，ra、rb用于量度管状、盘状和斑状等不同结构之间的特征差异，s是为对比度系数，用于描述目标和周围背景的亮度差异，s的值就越大，亮度差异越大；frangi血管增强函数v(σ,a)的取值范围在0到1之间，v(σ,a)越接近于1表示该点处的组织与管状结构的相似性最大；α、β和γ为三个调参系数，通过影响ra、rb和s在血管相似性函数中所占的比重来调整式(12)的灵敏度；ra、rb和s定义分别如式(13)、式(14)和式(15)表示：

步骤四、升主动脉分割以及左右冠状动脉起始层识别；结合光流法检测相邻两层升主动脉间光流位移距离最大值的变化和心脏医学解剖结构知识判断左右冠状动脉的起始层；

4-1由于冠状动脉的根节点位于升主动脉根部，升主动脉的分割和提取是识别左右冠状动脉起始层的前提；采用霍夫变换检测类圆形的升主动脉区域，以被检测到的类圆形圆心作为种子点，利用区域生长法分割每一层升主动脉，最后得到完整的升主动脉；

所述的霍夫变换检测类圆形的升主动脉区域属于常规技术，故不详解；

4-2使用lucas-kanade光流法检测相邻两层升主动脉间的最大光流及所有像素中的最大位移距离；

在cta图像序列中，大部份的升主动脉区域都呈现类圆形，但是当某一升主动脉层包含冠状动脉分支时，类圆形升主动脉区域会多一小段凸起的冠状动脉区域，该层相对于前一层在冠状动脉根结点处的光流位移距离会存在巨大的变化，通过对比相邻2层间光流位移距离最大值的变化即可初步筛选出冠状动脉的起始层；光流法假设相邻的2层cta图像强度为i(x,y,t)和i(x+δx,y+δy,t+δt)，并且相邻的2层图像强度不变，可以得到公式(16)

i(x+δx,y+δy,t+δt)＝i(x,y,t)(16)

将公式(16)使用泰勒级数展开后可得：

其中为高次项，大小可以忽略，为此结合公式(16)和(17)可得：

其中和分别为图像相对于x轴和y轴的空间导数，为时间导数，公式(18)可改写为

(ix,iy)(u,v)＝-it(19)

其中，u和v就是图像中每个像素分别沿着x轴和y轴的光流位移距离，但是只有公式(19)一个约束方程无法求解2个未知变量，所以光流法假设目标像素的周围像素也存在相同的光流位移距离，即在目标像素周围m*m的窗口内，每个像素点都具有相似的光流位移距离u和v；简单的来说，就是通过增加约束方程来求解光流矢量u和v，如公式(20)所示；

式(19)中(xn，yn)表示m*m的窗口内第n个像素点；

通过公式(20)求解所得的u和v可以计算出每个像素点的位移距离，如公式(21)所示；

根据每一层distances最大的点，筛选出最大的前10层做为候选冠脉起始层，然后结合三个心脏医学解剖结构先验知识进一步筛选出左右冠状动脉的起始层；

所述的三个心脏医学解剖结构先验知识具体如下：

第一，升主动脉根部只存在左右两个冠状动脉开口，在10层冠脉候选层中，位移距离最大的是左冠状动脉起始层sleft；

第二，左右冠脉起始层间相隔10层以上，且左冠状动脉起始层位于右冠状动脉起始层的z轴上方；

第三，左冠状动脉起始点ol位于该层质心的-90°～45°之间，右冠状动脉起始点or位于该层质心的位于45°～135°之间；

最后在满足上述第二、三条件的剩余候选冠脉起始层中选择层号最小的层作为右侧冠状动脉起始层sright；

步骤五、在步骤三增强后的cta图像上，从步骤四获得的左右冠状动脉起始层开始，使用结合端点检测的区域生长法分割每一层cta图像中的血管区域；具体是从步骤三获得的左右冠状动脉起始层sleft和sright开始，沿着z轴向下或向上利用迭代区域生长法进行冠状动脉的初始分割；

所述的迭代区域生长法由步骤5-1种子点自动选择、步骤5-2迭代区域生长法生长准则、步骤5-3端点检测以及步骤5-4端点再生长四部分组成，如图7具体是：

5-1种子点自动选择

利用血管上下层之间的空间连通性，即当前层si与相邻层si+1血管区域内的二维坐标存在部分重叠作为每层种子点的选择依据；由于三维的冠状动脉血管分叉众多，所以在二维的单层冠状动脉cta图像中可能存在多个独立的血管连通区域，当si层分割完成后，先检测血管连通区域个数c，然后分别在si层的每个独立的血管区域中以一定的步长选取m(m＝15)个种子点作为相邻层si+1血管区域生长的种子点；

举例：血管会有很多分叉，每个分叉在一层图像上就是一个区域，如果有三个分叉，那么这一层就有三个血管区域了，在每个血管区域中取15个种子点，一共会有15*3＝45个种子点；

5-2迭代区域生长法生长准则

利用迭代区域生长法分割每一层的冠状动脉血管区域；

5-2-1首先从5-1中选取的每一层c*15个种子点依次开始对该层血管区域进行分割，每一个种子点会将区域边缘像素点8领域内与种子点灰度值相差在1％以内的像素点不断归并到感兴趣区域中，直到不能生长为止；

5-2-2由于5-1中选取的种子点只有部分会落在血管区域内，所以每个种子点生长完成后，通过计算分割得到的血管感兴趣区域内的像素数量num判断该种子点是否在血管区域内：如果像素数num小于nthresh，则表明该种子点位于血管区域内，保留该种子点分割得到的血管感兴趣区域；如果像素数num大于nthresh，则说明该种子点落在血管区域之外，出现了严重的过分割现象，丢弃该种子点分割得到的区域，根据经验nthresh＝2000；然后将所有保留下来的血管感兴趣区域之和作为该层最终的冠状动脉血管分割结果；

5-2-3判断每一层最终的冠状动脉血管分割结果像素数是否为0，如果等于0，则表明无法再继续分割血管，结束迭代，得到初始冠状动脉感兴趣区域；如果不等于0，则返回5-1继续执行下一层种子点的选择以及5-2迭代区域生长法；

5-3端点检测

对5-2中提取的初始冠状动脉感兴趣区域利用细化方法提取血管骨架线；具体是利用血管骨架线的26邻域单连通性，将骨架线上的体素分为端点p1、普通连接点p2、分叉点p3；从各个端点p1和分叉点p3开始，继续使用结合5-1的种子点选择机制和5-2的迭代区域生长法尝试z轴向上和向下的双向生长，用来分割冠状动脉分支中拓扑结构复杂，上下起伏的血管部分；

所述端点p1在该体素的26邻域点中，只有1个邻域点位于骨架线上；所述普通连接点p2，在该体素的26邻域点中，只有2个邻域点属于骨架线上；所述分叉点p3，在该体素的26邻域点中，至少有3个邻域点位于骨架线上；

5-4端点再生长

判断结合端点检测的区域生长法继续分割后的冠状动脉区域是否有新的端点p1和分叉点p3，如果有，则返回5-4继续分割；如果没有，则结束分割。

为证明上述方法的可行性，下面采用具体的冠状动脉cta图像进行举例说明。本发明验证所采用的cta数据集来自山东大学齐鲁医学院，该数据集由年龄在51-85岁之间的14例男性和6例女性病例组成，均带有医生标注，每例数据都有大约200-300层包含冠脉血管的cta轴向切片，每层轴向切片图像分辨率为512*512，像素间距dx，dy，dz为0.35mm，0.35mm，0.80mm。

实验过程如下：先对原始的cta图像序列(如图1a所示)进行窗宽窗位调整和基于改进灰狼算法优化的多级阈值处理，分割得到的冠状动脉感兴趣区域如图1b所示，可有效抑制肺部血管等非心脏组织并去除大部分与冠状动脉灰度值相近的心房与心室等心脏组织，减少非冠状动脉组织对后续冠脉分割产生的影响；对预处理后的图1b进行基于多尺度融合的血管特征提取，结果如图2所示，不同直径的冠状动脉血管被增强；图3a为利用霍夫变换和区域生长分割得到的升主动脉三维重构图，图中升主动脉的左右凸起分叉就是左右冠状动脉的起始根结点，包含冠状动脉的升主动脉层不是普通的类圆形。本发明首先通过光流法计算每一层所有像素相对于前一层的光流位移距离，并保存每一层中的最大值，结果如图3b所示。在该例cta数据的90层升主动脉中，普通类圆形升主动脉层，不会出现位移激增的现象，对应于图3b曲线中较平缓的区域，但是当遇到包含冠状动脉根节点的不规则升主动脉层时，类圆形的升主动脉层由于多了一段冠状动脉凸起，会出现最大位移距离激增的现象，所以可将图3b中5个最大峰值的层作为冠状动脉候选层，再结合冠状动脉的医学解剖结构先验知识进一步筛选得到左右冠状动脉的起始层。首先，5层冠状动脉候选层中位移距离最大峰即图3b中第53层就是左冠状动脉支的起始层，其次第53层中最大位移像素点的坐标(图3c中的黑圈)位于该层升主脉质心的-90°～45°之间，也符合心脏医学解剖结构的先验构知识，本例第53层左冠状动脉支的起始层结果如图3d所示。右冠状动脉起始层在本数据集中距离左冠状动脉支的起始层在10层以上即8mm左右，所以63层以后，峰值最大的74层最可能是右冠状动脉，其最大位移像素点的坐标(图3e中的黑圈)位于该层升主脉质心的45°～135°之间，也符合心脏医学解剖结构的先验构知识，本例第74层右冠状动脉支的起始层结果如图3e所示。本发明在20例cta数据上对该算法进行了测试，结果如表1所示，左侧冠状动脉起始层检测成功19例，准确率达到95％，右侧冠状动脉起始层检测成功17例，准确率达到85％。

表1光流法识别冠状动脉起始层

最后从53层和74层分别使用结合端点检测的迭代区域生长法分割每一层cta图像中的血管区域。图4b和4c所示分别为传统区域生长法与结合端点检测的区域生长法的冠状动脉分割结果对比，图4c所示结合端点检测的区域生长法对于冠状动脉末端拓扑结构复杂，上下起伏的细小分支有更好的分割能力。为了评估算法的分割结果精度，本发明使用dice、jaccard、msd(meansquaredsurfacedistance)、maxsd(maximumsurfacedistance)作为分割评价指标，结果如表2所示，采用本发明提出的结合端点检测的区域生长法得到的分割结果，dice和jaccard系数相比于传统区域生长法，分别平均提高了4％和8％，达到了0.70和0.57，msd和maxsd相比于传统区域生长法，分别平均降低了2％和4％，达到了0.40和2.66。

表2不同方法结果比较

本发明提出了一种面向cta数据的冠状动脉自动分割方法。首先，通过图像预处理有效抑制了非冠脉组织，提升了冠状动脉和背景的对比度；其次，结合光流法与心脏解剖结构先验知识检测带冠状动脉分叉的不规则升主动脉层，避免了手动初始化左右冠状动脉的起始点；最后，相比于传统区域生长法，本文提出的结合端点检测的自适应区域生长法对于灰度不均且拓扑结构复杂的细小分支有更好的分割能力和准确性。