技术摘要：
本发明公开了一种基于CT图像的冠状动脉中心线提取方法。现有技术由于多种条件限制，不能较完整的、自动的提取出CT图像中的冠状动脉中心线；本发明如下：一、消除被处理的CT图像的散粒噪声。二、对CT图像内的所有像素点进行血管区域增强。三、在CT图像的血管区域内使用  全部
背景技术：
心脏作为人体中，最重要的器官之一，负责向循环系统提供动力，以保证身体各处 正常的新陈代谢，而冠状动脉负责为心脏提供所需养料，带走代谢废物，是保证心脏正常生 理功能重要器官。近年来，由于人们日益重视维护心脑血管健康，冠心病(又称冠状动脉粥 样硬化)作为一种高发的心血管疾病，得到国际社会的广泛关注。《中，国心血管病报告  2018》中，指出经推算我国现有冠心病患者1100万人，且死亡率自2012年有明显的上升趋势。 冠状动脉附着于心脏表面，为心脏提供足量的营养与能量，它源起升主动脉根部， 向左、向右各发散出左冠状动脉(Left  Coronary  Artery，LCA)和右冠状动脉(Right  Coronary Artery，RCA)，两支动脉紧紧包裹着心脏，形似一顶皇冠，倒置与心脏之上。LCA保 证心脏前壁和侧壁的心肌供血，在其分支处10-25mm后分出左前降支(Left  Anterior  Descending Artery，LAD)和回旋支(Left  Circumflex  Artery，LCX)，前者在左右心室间延 伸至心尖，LCA  分支之前的部分被称为左主干(Left  Main  Artery，LM)；RCA负责向右室、右 房、左室部分、后壁、窦房结、房室结供血，从中，又发散出锐缘支(Right  Marginal  Artery， MA)和后降支  (Posterior  Descending  Artery，PDA)，其中，LCA与RCA为冠状动脉的主要血 管，其余的分支都是由它们发散出去。 CTA图像是将血管特征利用计算机重构出来，可以在血管的横截面上提供三维血 管的原始层切图像，可根据此图像分析冠状动脉的病变情况。现有技术中，不能自动提取出 CT  层切图像中的三维血管；近年来基于冠心病的计算机辅助诊断系统(Computer  Aided Diagnosis，CAD)成为研究的热点，而冠状动脉中心线的提取又是CAD研究中的一个难点。在 CAD中冠状动脉中心线是一个定量分析冠心病等心血管疾病重要参数，是对血管狭窄分析、 血管钙化的定量分析、血管斑块类型的判定等多种心血管疾病定量分析的起点。但是，目前 现有的冠状动脉中心线提取方法不能很好的适应冠状动脉树上所有分支的形态及走向。因 此，获得一种能够在各种情况下精确提取冠状动脉中心线的方法对于CAD来说尤为重要。此 外，冠状动脉中心线对多种冠状动脉的可视化重建也有着不可替代的作用，尤其是在曲面 重建技术和多平面重建技术中，发挥了巨大的作用。同时，冠状动脉的中心线还可以为冠状 动脉手术治疗提供术前的路径规划和术中，的路径引导，如经皮冠脉介入治疗和冠脉旁路 移植术。因此，获取精准的冠状动脉中心线对于临床治疗有极大的意义。
技术实现要素：
本发明的目的在于提供一种基于CT图像的冠状动脉中心线提取方法。 本发明的具体步骤如下： 6 CN 111612743 A 说　明　书 2/12 页 步骤一、消除被处理的CT层切图像的散粒噪声。 步骤二、对CT层切图像内的所有体素点进行血管区域增强。 步骤三、经过血管增强的图像血管区域的轮廓被增强，在血管区域内使用区域生 长法，对CT层切图像中的每一层图像中的血管区域进行分割；使得CT层切图像变为二值图像。 步骤四、对经步骤三分割后的CT层切图像的每层图像分别根据原始数据中的空间 分辨率进行三维重构，得到三维空间中的血管数据。 步骤五、对血管数据进行细化，得到血管的初步的中心线；血管数据上各个体素点 的值为1，作为目标体素点；对所有目标体素点进行如下操作：提取出目标体素点的六个6- 邻接点、十二个18-邻接点和八个26-邻接点；对步骤四重构得到的血管数据采用十二方向 的迭代细化处理；每个迭代步骤均对血管数据的所有边界点依次进行十二个方向的细化处 理；十二个方向分别为目标体素点指向十二个18-邻接点的方向。 步骤六、将步骤五中提取出的中心线转换为有向图。在有向图中，仅有一个邻接体 素的体素点视为端点，有且仅有两个邻接体素的体素点视为普通连接点，存在多于两个邻 接体素的体素点视为节点；相邻的两个节点或端点之间的各普通连接点形成一条有向边。 根据各有向边上的节点数量筛选出中心线和分支。 作为优选，步骤五中，若一轮细化处理中删除一个或多个体素点，则再进行一轮迭 代，直至迭代到一轮细化处理中未有能够被删除的体素点。 作为优选，步骤五中，用U、D、W、E、S、N表示点目标体素点p的6-邻接点；用  US、NE、 WD、ES、UW、ND、SW、UN、ED、NW、UE、SD分别表示十二个18-邻接点；组成18-邻接点的两个字母 表示该18-邻接点与那两个6-邻接点相邻。用点A、B、C、F、G、H、I、J分别表示八个26-邻接点； 点A为点U、S、W所在平面的交点；点B为点U、W、N所在平面的交点；点C为点U、N、E所在平面的 交点；点F为点U、S、E  所在平面的交点；点G为点D、S、W所在平面的交点；点H为点D、W、N所在 平面的交点；点I为点D、N、E所在平面的交点；点J为点D、S、E所在平面的交点；十二个方向的 细化处理顺序为：US、NE、WD、ES、UW、ND、SW、UN、ED、NW、UE、SD。十二方向的细化处理各自对应 一组模板；一组模板内有14个模板。在一个方向的细化处理中，若目标体素点的3×3×3邻 域与对应的14个模板中任意一个或多个相符合，则将目标体素点删除。 US方向的十四个模板如下： 第1个模板中，目标体素点p的3×3×3邻域中顶层的九个体素点均标记为“○”，D 点标记为“●”，其余各点均标记为“X”。 第2个模板中，目标体素点p的3×3×3邻域中的点A、US、F、SW、S、ES、G、SD、  J均标 记为“○”，N点标记为“●”，其余各点均标记为“X”。 第3个模板中，目标体素点p的3×3×3邻域中顶层的九个体素点和点SW、S、SE、G、  SD、J均标记为“○”，N点和D点标记为“·”，ND点标记为“●”，其余各点均标记为“X”。 第4个模板中，目标体素点p的3×3×3邻域中点UW和SW标记为“V”，UE点和ES  点标 记为“W”，N点和D点标记为“●”，点U、A、US、F、S标记为“○”，其余各点均标记为“·”。 第5个模板中，目标体素点p的3×3×3邻域中点UW和SW标记为“V”，点UE和ES  标记 为“Z”，点F、N、D标记为“●”，点U、A、US、S标记为“○”，其余各点均标记为“·”。 第6个模板中，目标体素点p的3×3×3邻域中UW点和SW点标记为“Z”，UE点和ES  点 标记为“V”，点A、N、D标记为“●”，点U、F、US、S标记为“○”，其余各点均标记为“·”。 7 CN 111612743 A 说　明　书 3/12 页 第7个模板中，目标体素点p的3×3×3邻域中点UW和SW标记为“V”，点E、N、D  标记 为“●”，点U、A、US、S标记为“○”，其余各点均标记为“·”。 第8个模板中，目标体素点p的3×3×3邻域中UE点和ES点标记为“V”，点W、N、D  标 记为“●”，点U、F、US、S标记为“○”，其余各点均标记为“·”。 第9个模板中，目标体素点p的3×3×3邻域中UW点和SW点标记为“Z”，点A、N、E、  D 标记为“●”，点U、US、S标记为“○”，其余各点均标记为“·”。 第10个模板中，目标体素点p的3×3×3邻域中点UE和ES标记为“Z”，点F、N、W、  D标 记为“●”，点U、US、S标记为“○”，其余各点均标记为“·”。 第11个模板中，目标体素点p的3×3×3邻域中顶层的九个体素点和点SW、S、G、SD  均标记为“○”，ND点和ED点标记为“●”，其余各点均标记为“·”。 第12个模板中，目标体素点p的3×3×3邻域中顶层的九个体素点和点ES、S、SD、J  均标记为“○”，ND点和WD点标记为“●”，其余各点均标记为“·”。 第13个模板中，目标体素点p的3×3×3邻域中前侧的九个体素点和点B、UN、UW、  U、A、US均标记为“○”，NE点和ND点标记为“●”，其余各点均标记为“·”。 第14个模板中，目标体素点p的3×3×3邻域中前侧的九个体素点和点UN、C、U、UE、  US、F均标记为“○”，，NW点和ND点标记为“●”，其余各点均标记为“·”。 符号“○”表示该位置为背景点；符号“●”表示该位置为目标点；符号“·”表示该 位置上既可以为目标点也可以为背景点；同一个模板中的字母“X”表示在这些位置上至少 有一个为目标点；同一个模板中的字母“V”和“W”各自表示这些位置上至少有一个为背景 点；同一个模板中有两个不同的位置标注了字母“Z”，则表示这两个位置上必须一个为目标 点，另一个为背景点。 剩余十一个方向对应的模板均通过US方向对应的十四个模板通过旋转、镜像对称 得到；具体如下：以目标体素点p为坐标原点，点p到点E为X轴正方向，p到点N为Y轴正方向，  p到点U为Z轴正方向建立空间直角坐标系； UN方向的模板由US方向的模板根据AGJF平面镜像对称得到；UW方向的模板由US  方向的模板根据AGIC平面镜像对称得到；UE方向的模板由US方向的模板根据BHJF平面镜像 对称得到；SW方向的模板由US方向的模板根据AJIB平面镜像对称得到；ES方向的模板由US 方向的模板根据FGHC平面镜像对称得到；NW方向的模板由UN方向的模板根据  AJIB平面镜 像对称得到；NE方向的模板由UN方向的模板根据FGHC平面镜像对称得到；  WD方向的模板由 SW方向的模板根据BGJC平面镜像对称得到；SD方向的模板由US方向的模板根据HGJI平面镜 像对称得到；ND方向的模板由UN方向的模板根据HGJI平面镜像对称得到；ED方向的模板由 UE方向的模板根据HGJI平面镜像对称得到。 作为优选，步骤六中，设置一个长度阈值，将小于长度阈值的分支作为需要被删除 的细小分支，大于长度阈值的分支为需要被保留的中心线起点和终点所在的有向边。 作为优选，步骤一的具体过程如下：采用异性扩散滤波器进行散粒噪声的消除。CT 层切图像上的任意点x的各向异性扩散方程的表达式如式(1)所示： 式(1)中， 表示图像在点x的灰度， 表示CT层切图像中，点x的坐标所对应 8 CN 111612743 A 说　明　书 4/12 页 的向量，t表示迭代次数； 表示梯度算子， 表示 的梯度，div(·)表示散度 算子； 表示扩散函数，其表达式如式(2)所示： 式(2)中，k为与噪声级和边界强度相关的常数；求解式(2)的偏微分方程，得到消 除散粒噪声后的滤波结果 其表达式如式(3)所示： u(t Δt)＝u(t) Δt(Φeast-Φwest Φnorth-Φsouth)#(3) 式(3)中，Φeast、Φwest、Φnorth、Φsouth分别为点x的东侧、西侧、上侧、下侧邻接像素 点的扩散通量。Δt表示迭代步长。 作为优选，步骤二中，以CT层切图像内所有体素点分别作为被分析点A进行如下操作： 2-1、本发明利用泰勒级数的展开式分析体素点的局部特征；被分析点A在所在图 像I中的泰勒级数展开式I(A ΔA)如式(4)所示： 式(4)中，I(A)表示图像I中的被分析点A； 表示图像I在被分析点A处的梯度， Δ(·)表示拉普拉斯算子；AT表示被分析点A的转置；H(A)表示被分析点A的Hessian矩阵， 其表达式如式(5)所示。 式(5)中，矩阵H(A)内各个元素均为图像I在被分析点A处的二阶偏导；矩阵H(A)存 在三个特征值λ1、λ2和λ3，以及三个特征值各自对应的特征向量为 和 三个特征 值满足|λ1|≤|λ2|≤|λ3|的关系。 2-2、利用基于Hessian矩阵的Frangi血管增强函数对血管区域进行增强，对非血 管区域进行抑制，血管相似性函数vo(λ)如式(6)所示： 式(6)中，RA、RB为椭球体的两个定义参数；S为对比度系数；α、β和γ为三个调参系 数。  RA、RB和S分别用式(7)、式(8)和式(9)表示： 9 CN 111612743 A 说　明　书 5/12 页 本发明具有的有益效果是： 本发明提出了一种在三维空间中，的针对冠状动脉的拓扑细化算法，该算法在提 取中心线的同时尽可能的保持了冠状动脉原有的拓扑结构；同时，针对细化过程中，所产生 的细小分支，本发明采用Dijkstra算法予以消除。本发明提高了基于CT图像的冠状动脉中 心线的提取精度与效率，优化效果较为明显。 附图说明 图1为本发明验证时使用的原始CTA图像； 图2a为本发明验证过程中，对原始CTA图像进行各项异性扩散滤波后的结果； 图2b为本发明验证过程中，进行血管增强之后的结果； 图2c为本发明验证过程中，对分割出的冠状动脉血管进行三维重构后的结果； 图3a为本发明提出的十二方向细化算法的十二个方向的示意图； 图3b为本发明提出的十二方向细化算法的18-邻接点和26-邻接点的示意图； 图4a-4n为本发明所用进行十二方向拓扑细化算法时所需的细化模板； 图5为本发明验证过程中，对冠状动脉中心线提取的初步结果； 图6为本发明进行细小分支删除后的冠状动脉中心线提取的最终结果。