心脏冠脉血管自动分段命名方法与流程

本发明涉及心脏冠脉血管图像处理领域，具体涉及一种心脏冠脉血管自动分段命名方法。

背景技术：

对于心脏相关疾病的诊断而言，对心脏左右冠脉的血管识别和命名是至关重要的一步。然而心脏冠脉血管不同的个体有各种各样的变异，有的人只有十几根血管，而有的人则有多几十甚至上百根血管，所以对心脏冠脉血管的划分识别一直是一个非常耗时也耗费人工的工作。

技术实现要素：

本发明的目的在于提供一种基于心脏冠脉血管中心线数据的心脏冠脉血管自动分段命名方法。

为实现上述目的，本发明采用以下技术方案：

心脏冠脉血管自动分段命名方法，包括：

s1、提取心脏冠脉3d图像的血管中心线，定义血管中心线中各点的三维坐标；

定义沿着人的两肩方向为x轴，人左手所在方向为x轴正方向；

定义沿着人的胸口方向，垂直于x轴为y轴正方向；

定义沿着人体向下方向为z轴正方向；

s2、基于血管中心线从心脏冠脉血管中识别左冠脉及右冠脉；

s21、从呈两棵树状结构的心脏冠脉血管中分别寻找最长的血管；

s22、从获得的两根最长的血管的起始点开始往其延伸方向寻找第m个像素点，将第m个像素点处的x坐标值较大的血管标记为左冠脉，将第m个像素点处的x坐标值较小的血管标记为右冠脉；m∈[30,80]。

进一步地，m＝50。

进一步地，还包括步骤s3：从识别的右冠脉中识别rca；

步骤s3具体执行以下步骤：

s31、在右冠脉中寻找在y轴方向走向最大的血管，标记为rca_possible；

s32、比较rca_possible的终点的y轴坐标和右冠脉起点的y轴坐标，若rca_possible的终点的y轴坐标大于右冠脉起点的y轴坐标，则将rca_possible认定为rca，否则进入步骤s33；

s33、在右冠脉血管中寻找在z轴走向最大的血管，标记为rca。

进一步地，所述步骤s3还包括识别r-plb以及r-pda的步骤，其具体执行以下步骤：

s34、确定左冠脉的起始点leftroot和右冠脉的起始点rightroot，计算他们的中点坐标midroot[x，y，z]；

s35、在所识别的rca的底部平缓部分寻找所有在rca的x轴正方向或者y轴负方向的分叉，将这些血管按照在rca上分叉点出现的顺序保存进pda_plb_possible_list，并进入步骤s36；若找不到任何一条符合条件的血管组成pda_plb_possible_list，则r-pda和r-plb不存在，步骤s3终止；

s36、计算pda_plb_possible_list中所有血管在xy平面上和midroot点的距离，将距离最小的血管识别为r-pda；

s37、在pda_plb_possible_list中排在r-pda之后的血管中，寻找一根最长的血管，识别为r-plb；若pda_plb_possible_list在r-pda之后没有血管，则r-plb不存在，步骤s3终止。

进一步地，还包括步骤s4：从识别的左冠脉中识别lm、lad及lcx；

s41、定义lad与lcx重叠的部分为lm，限定lm的总长度值的最大值为lmmax；

s42、采用循环策略，从左冠脉中识别lad及lcx：

s421、定义最短分支血管长度值阈值区间[n1,n2]；以阈值n1开始循环，删除左冠脉中长度小于n1的分支血管；则对于后续循环有n'＝(n+n)，式中，n'为当前阈值，n为上次循环阈值，n为循环增量，删除左冠脉中长度小于n'的分支血管；

s422、在每次循环中寻找可能的lad与lcx组合，将整个循环结果中识别的组合可能记录为lcx_lad_list；

可能的lcx为：在lm后的分支中寻找走向为y轴正方向或者x轴负方向，且在z轴上取得最大值的血管；可能的lad为：在lm后的分支中寻找走向为y轴负方向或者x轴正方向，且在z轴上取得最大值的血管；

s423、访问lcx_lad_list中的每一组lad和lcx组合，获取每组lad和lcx组合的len(lm)和abs(len(lad)–len(lcx))，寻找其中abs(len(lad)–len(lcx))值最小且满足len(lm)＜lmmax的那组lad与lcx组合作为识别结果；式中，len(lm)为lm血管的长度，abs(len(lad)–len(lcx))为lad和lcx长度差的绝对值；

s43：在步骤s42所识别的lcx上查找沿着z轴向下且相对于lcx沿x轴负方向并且沿z轴正方向延伸的血管，将寻找到的血管定义为新的lcx；若寻找不到符合条件的血管，保留步骤s42所识别的lcx结果。

进一步地，lmmax＝120个像素点，n1＝30个像素点，n1＝80个像素点。

进一步地，还包括从识别的lcx中识别om1、om2及om3的步骤s5，步骤s5具体执行以下步骤：

s51、识别lcx上的所有分支血管，组成om1om2list；

s52、定义距离比例drxy＝dx/dy，式中，dx为lad与lcx在xz平面上投影在x轴方向上的平均距离，dy为lad与lcx在yz平面上投影在y轴方向上的平均距离；

s53、给定距离比例阈值drc，进行om1、om2及om3的识别：

当drxy＜drc时，在om1om2list中按照各分支血管在lcx上的出现顺序寻找3根相对于lcx沿y轴负方向走向的分支血管，获得om1、om2及om3的识别结果；

当drxy＞drc时，在om1om2list中按照各分支血管在lcx上的出现顺序寻找3根相对于lcx沿x轴正方向走向的分支血管，获得om1、om2及om3的识别结果。

进一步地，所述距离比例阈值drc＝0.13。

进一步地，还包括从识别的lad中识别d1及d2的步骤s6，步骤s6具体执行以下步骤：

s61、按照在lad上出现的顺序保存在lad上所有分支血管进dlist；

s62、在dlist上依次找出三根相对lad沿着y轴正方向走向的血管，获得d1、d2及d3的识别结果。

进一步地，还包括从识别的lcx及lad中识别ri的步骤s7；

s7：计算om1或者d1距离lm的距离，若om1距离lm的距离不超过6个点，则将步骤s53识别的om1、om2、om3分别修改为ri、om1及om2；若d1距离lm的距离不超过6个点，则将步骤s62识别的d1、d2及d3分别修改为ri、d1及d2。

采用上述技术方案后，本发明与背景技术相比，具有如下优点：

本发明基于心脏冠脉血管的中心线数据，能够适应各类变异状况，可以准确的从几十到上百根血管中自动识别出需要的的十几根目标血管并准确命名，实现了心脏冠脉血管的自动分段命名。

附图说明

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

图2为心脏冠脉3d图像的血管中心线提取流程图；

图3为人体心脏冠脉的中心线示意图；

图4为冠脉rca，r-pda和r-plb的示意图；

图5干扰血管示意图；

图6为lad与lcx示意图；

图7为lcx调整为lcx_updated的示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

实施例

如图1所示，心脏冠脉自动分段命名方法，其包括：

s1、提取心脏冠脉3d图像的血管中心线，定义血管中心线中各点的三维坐标；

s2、从心脏冠脉血管中识别左冠脉及右冠脉；

s3、从识别的右冠脉中识别rca、r-pda以及r-plb；

s4、从识别的左冠脉中识别lm、lad及lcx；

s5、从识别的lcx中识别om1、om2及om3；

s6、从识别的lad中识别d1、d2及d3；

s7、从识别的lad及lcx中识别ri。

如图2所示，心脏冠脉3d图像的血管中心线通过如下方式生成：

s11、获取患者心脏冠脉血管的原始ct图像；

s12、分析所述原始ct图像，三维重构出心脏的冠脉3d图像；

s13、提取心脏冠脉3d图像的血管中心线，定义并保存血管中心线中各点的三维坐标。

三维坐标定义如下：定义沿着人的两肩方向为x轴，人左手所在方向为x轴正方向；定义沿着人的胸口方向，垂直于x轴为y轴正方向；定义沿着人体向下方向为z轴正方向。

如此，最终生成的血管中心线被保存在vtk格式的文件中，这些vtk文件可以保存心脏冠脉中每根血管从起点到终点每个点的三维坐标，因而，最终得到的中心线信息是每条血管的点的信息，这里的每个点就是对应的ct图象上的一个截面。后续步骤s2-s7实质上是根据血管中心线而进行的分段及命名。

在对步骤s2-s7进行论述之前，首先对人体心脏冠脉的大体概况进行介绍。如图3所示的是人体心脏冠脉的中心线示意图，图中，水平面为xy平面，竖直且平行于纸面的平面为xz平面，竖直且垂直于纸面的面为yz平面，从左往右是x轴正方向，从上往下是z轴正方向，从里往外时y轴正方向。图中rca为右冠脉、r-plb为右冠起源后侧支、r-pda为右后降支、lm为左主干、d1为第一对角支、d2为第二对角支、lad为左前降支、om1为第一钝缘支、om2为第二钝缘支、lcx为左回旋支。

可以看出，左右冠脉实际上是两颗树状结构，其中左冠脉所有血管的起点都是从左心房主动脉发出，右冠脉所有血管都从右心房主动脉发出。因此，对于这两颗冠脉树如何区分那一颗树是左冠脉那一颗是右冠脉是划分工作的第一步。由于左冠脉在心脏左侧，右冠脉在心脏右侧，理论上我们只需要比较两个root点的x轴大小就可以判断，可是由于某些变异情况，有些人的左右冠脉血管会发生交叉。因此步骤s2具体执行以下步骤来判断左冠脉和右冠脉：

s21、从呈两棵树状结构的心脏冠脉血管中分别寻找最长的血管，分别记为tree1vessel及tree2vessel；

s22、从tree1vessel及tree2vessel的起始点(root)开始往其延伸方向寻找第m个像素点，记录为tree1point和tree2point，比较tree1point和tree2point的x坐标值，将其x轴坐标值较大的血管标记为左冠脉，将其坐标值较小的血管标记为右冠脉，m∈[30,80]。由于任何冠脉的最长血管都不可能少于50个像素点，因此，本实施例中优选m＝50。

区别完左右冠脉，既可以针对左右冠脉划分其内部的血管。对于右冠脉，所需识别的血管主要有rca，r-pda，r-plb。其中rca是一定存在的，而r-pda和r-plb则不一定存在。如图4所示的是一个心脏有冠脉rca，r-pda和r-plb的示意图。

对于rca，申请人通过分析大量冠脉血管中心线数据得出结论，从大概率上(90％)以上，rca是右冠脉中沿着患者胸口正面延伸最远(y轴方向取得最大值)的那根血管。但是由于存在一些变异，因此，不能单纯地依靠这一结论进行rca的判断。

因此，步骤s3具体执行以下步骤：

s31、在右冠脉中寻找在y轴方向走向最大的血管，标记为rca_possible；

s32、比较rca_possible的终点的y轴坐标和右冠脉起点的y轴坐标，若rca_possible的终点的y轴坐标大于右冠脉起点的y轴坐标，则将rca_possible认定为rca，否则进入步骤s33；

s33、在右冠脉血管中寻找在z轴走向最大的血管，标记为rca。

rca从外观上看很像一个反写的字母c，在rca的开始部分，血管的走向通常是沿着x轴负方向平缓的下降(c字母的上端)；然后开始迅速的接近垂直的沿着z轴正方向快速下降(字母c的中段)；当下降到心脏底部以后，又会沿着y轴正方向平缓的前行(字母c的下端)。

在识别了rca以后，下一步就是识别r-plb和r-pda，应当注意，这两根血管有可能不存在，因此，步骤s3还执行以下步骤：

s34、确定左冠脉的起始点leftroot和右冠脉的起始点rightroot，计算他们的中点坐标midroot[x，y，z]；

s35、在所识别的rca底部平缓部分(具体为z轴最大的部分，该部分z轴坐标变化相对平缓)寻找所有在rca的x轴正方向或者y轴负方向的分叉，将这些血管按照在rca上分叉点出现的顺序保存进pda_plb_possible_list，并进入步骤s36；若找不到任何一条符合条件的血管组成pda_plb_possible_list，则r-pda和r-plb不存在，步骤s3终止；

s36、计算pda_plb_possible_list中所有血管在xy平面上和midroot点的距离，将距离最小的血管识别为r-pda；

s37、在pda_plb_possible_list中排在r-pda之后的血管中，寻找一根最长的血管，识别为r-plb；若pda_plb_possible_list在r-pda之后没有血管，则r-plb不存在，步骤s3终止。

对于左冠脉，所需识别的血管有：lad，lcx，ri，d1，d2，om1，om2。

首先，进行lm、lad和lcx的识别。理论上，左冠脉所有血管从左心房出发然后分成两个血管，其中向y轴正方向或则x轴负方向延伸的是lcx，向y轴负方向或者x轴正方向延伸的是lad。因此，lad和lcx重叠的部分定义为lm，在从lm的终点之后分叉的血管中寻找lcx和lad。

理想状态下，找到左冠脉所有血管重合的部分，就是lm；在lm之后分叉的血管中找一根沿着y轴正方向，或者x轴方向延伸，并且在z轴能够取得最大值的血管，定义为lcx；在lm之后分叉的血管中找一根沿着x轴正方向，或者y轴负方向延伸，并且在z轴取得最大值的血管，定义为lad。

该方法可以实现60％左右的lm、lad和lcx的快速识别，可是由于冠脉血管中的变异，如图5所示，往往在lm上会有长短不一的干扰血管，很容易被误认为lad和lcx。如图5中出现的血管就很容易被算法误判为lad，导致划分失败。

为了避免被这类干扰血管所误导，最有效的做法就是在寻找lad和lcx前把他们删除。然而有时候这些看起来比较短的干扰血管就是真实的lad或者lcx，直接进行删除极易导致误删，同时，如何定义干扰血管的长度，也是所需解决的问题。

因此，本实施例中步骤s4具体执行以下分步骤：

s41、明确lm即是lad与lcx重叠的部分，同时lm的长度有限，限定lm的总长度值的最大值lmmax不超过120个像素点。

s42、采用循环策略，从左冠脉中识别lad及lcx。

s421、定义最短分支血管长度值阈值区间[n1,n2]；以阈值n1开始循环，删除左冠脉中长度小于n1的分支血管；则对于后续循环有n'＝(n+n)，式中，n'为当前阈值，n为上次循环阈值，n为循环增量，删除左冠脉中长度小于n'的分支血管。本实施例中，n1＝30，n2＝80，n＝10。

s422、在每次循环中寻找可能的lad与lcx组合，将整个循环结果中识别的组合可能记录为lcx_lad_list；

可能的lcx为：在lm后的分支中寻找走向为y轴正方向或者x轴负方向，且在z轴上取得最大值的血管；可能的lad为：在lm后的分支中寻找走向为y轴负方向或者x轴正方向，且在z轴上取得最大值的血管；

s423、访问lcx_lad_list中的每一组lad和lcx组合，获取每组lad和lcx组合的len(lm)和abs(len(lad)–len(lcx))，寻找其中abs(len(lad)–len(lcx))值最小且满足len(lm)＜lmmax的那组lad与lcx组合作为识别结果；式中，len(lm)为lm血管的长度，abs(len(lad)–len(lcx))为lad和lcx长度差的绝对值。

如图6所示，有三对疑似lcx和lad，分别是(lcx，lad1)，(lcx，lad2)，(lcx，lad3)。其中(lcx，lad3)这个组合的lm长度超过了lmmax，所以排除；对于(lcx，lad1)和(lcx，lad2)，我们可以得出abs(len(lad)–len(lad2))<abs(len(lad)–len(lad1))。所以，lcx和lad2可以作为最终确定的lcx和lad。

然而，步骤s42识别的lcx并不一定总是准确的，有时候识别出来的lcx有可能是om1或者om2，因此还需要执行步骤s43。

s43：在步骤s42所识别的lcx上查找沿着z轴向下且相对于lcx沿x轴负方向并且沿z轴正方向延伸的血管，将寻找到的血管定义为新的lcx；若寻找不到符合条件的血管，保留步骤s42所识别的lcx结果。

识别出lm、lad及lcx后，在lcx上指向lcx和lad夹角内侧的血管叫做om1和om2。在lad上指向lcx和lad夹角内侧的血管叫做d1和d2。

常规情况下，lad和lcx这两根血管在x轴方向上的距离比较远，y轴方向距离比较近(如果按yz平面投影，这两根血管比较接近；按xz平面投影，这两根血管比较远)。但是在变异情况下，lad和lcx在x轴方向上距离很近，在y轴方向上比较远(如果按yz平面投影，这两根血管比较接远；按xz平面投影，这两根血管比较近)。

针对这两种情况，需要定义变量drxy(dr为distanceratio的缩写)来判断这两根血管在x，y轴方向上的相对位置，这个变量返回的是lad和lcx在x轴和y轴方向上距离的比例。

因此，步骤s5具体执行以下步骤：

s51、识别lcx上的所有分支血管，组成om1om2list。

s52、定义距离比例drxy＝dx/dy，式中，dx为lad与lcx在xz平面上投影在x轴方向上的平均距离，dy为lad与lcx在yz平面上投影在y轴方向上的平均距离。

所述dx及dy通过以下方式计算：

首先，获得lad及lcx分别在xz平面与yz平面上的投影；

其次，在xz平面上，对于lad和lcx分别取z轴相同的点，计算其在x轴上的坐标差值，该坐标差值的平均值即为dx；在yz平面上，对于lad和lcx分别取z轴相同的点，计算其在y轴的坐标差值，该坐标差值的平均值即为dy。

s53、给定距离比例阈值drc(distanceratiocriterion的缩写)，进行om1、om2及om3的识别，本实施例中，drc＝0.13。

则，当drxy≤0.13时(这意味着lad和lcx在x轴上距离很近)，则om的定义就是：相对lcx，沿着y轴负方向走向的血管。因此在om1om2list中按照各分支血管在lcx上的出现顺序寻找3根相对于lcx沿y轴负方向走向的分支血管，即可获得om1、om2及om3的识别结果。

当drxy＞drc时(这意味着lad和lcx在x轴方向上距离比较远)，这个时候om血管的定义就是：相对lcx，沿着x轴的正方向走的血管。因此，在om1om2list中按照各分支血管在lcx上的出现顺序寻找3根相对于lcx沿x轴正方向走向的分支血管，即可获得om1、om2及om3的识别结果。

d1，d2及d3是lad上三根分支血管，他们对应着lcx上om1，om2，及om3。特征是：这三根分支血管都应该是走向都是指向lad和lcx夹角内侧的。因此，步骤s6执行以下步骤：

s61、按照在lad上出现的顺序保存在lad上所有分支血管进dlist；

s62、在dlist上依次找出三根相对lad沿着y轴正方向走向的血管，获得d1、d2及d3的识别结果。

前述步骤在lcx上找了om1、om2及om3；在lad上找了d1、d2及d3，但实际上所需的只有om1、om2和d1、d2。寻找om3和d3只是为了寻找ri。ri的定义是夹在lad和lcx正中间的一根血管，但实际上由于从ct图像中提取血管中心线中产生的误差，这个ri通常并不是正好夹在lad和lcx正中间，而是移动到了lad或者lcx上。

因此，步骤s7执行以下步骤：计算om1或者d1距离lm的距离，若om1距离lm的距离不超过6个点，则将步骤s53识别的om1、om2、om3分别修改为ri、om1及om2；若d1距离lm的距离不超过6个点，则将步骤s62识别的d1、d2及d3分别修改为ri、d1及d2。

本发明在江苏省自然科学基金(项目编号:bk20160971)和中国国家自然科学基金(项目编号:61703212)资助下完成。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。