[发明专利]基于冠状动脉造影图像分割的SYNTAX自动评分方法

摘要

基于冠状动脉造影图像分割的SYNTAX自动评分方法，所述方法依据患者冠状动脉造影图像，采用血管边支分离和再连接及K主曲线的方法，分割和提取冠状动脉血管树和中心线，生成冠状动脉血管树二值图像；基于SYNTAX评分标准，采用沿中心线方向的追踪方法，自动匹配血管节段编号，选择冠脉左(右)优势类型，结合血管横截面宽度衰减程度和中心线曲率变化，识别和评分冠脉病变不良特征，并计算上述评分，给出合理评分结果，辅助医生确定最优的手术方式。本发明规避了临床上现有的人工计算SYNTAX评分的计算量大、计算步骤繁琐、评分标准复杂、计算时间长等问题。

主权项

基于冠状动脉造影图像分割的SYNTAX自动评分方法，其特征在于：该方法包括肝位和右前斜位两种体位冠状动脉造影图像的读取与显示、冠脉造影图像血管的增强、冠脉造影血管及其分支的识别、血管中心线的提取、血管直径的计算、冠脉树节段的标记、病变部位与病变特征的识别，及SYNTAX自动评分；详细实施流程如下：步骤1，根据冠状动脉造影投照的两种体位，同时读取和显示肝位和右前斜位两种体位的冠状动脉造影图像；步骤2，同时对选定的两个体位的目标冠脉造影图像进行血管增强；步骤3，对增强后的血管进行主干和分支的分离，并将分离后的各血管片段进行标记；步骤4，将分离后的血管主干与分支进行再连接，生成完整的冠脉树二值图像；步骤5，基于完整的冠脉树二值图像，对血管的中心线即血管骨架进行提取；步骤6，基于血管中心线，对图像所有像素位置的血管进行直径的计算和存储；步骤7，基于步骤3已标记的各血管片段和步骤5中所计算的血管方向，依据SYNTAX评分中冠脉16节段法，进行左或右优势型的识别和追踪标记血管节段；步骤8，基于上述已标记的冠脉树图像，结合步骤6所存储的所有像素位置的血管直径，依据SYNTAX评分中病变特征和病变部位的直径变化，识别出对应的血管病变；步骤9，依据SYNTAX评分标准中冠脉树各节段所占的权重因数和各病变特征的评分，自动给出合理的SYNTAX评分，从而辅助医生确定最优的手术方式；所述步骤1读取和显示右前斜位和肝位两种体位的冠脉造影图像方法如下，使用者通过人工交互界面，自主选择导入右前斜位和肝位的目标造影图像；其中，肝位的冠脉造影图像用于右冠状动脉血管信息的提取，右前斜位的冠脉造影图像用于左冠状动脉血管信息的提取；所述步骤2对选定目标冠脉造影图像进行血管增强的方法如下：基于Hessian矩阵的Frangi血管滤波方法，有效地增强造影图像中的目标冠脉树；该方法利用Hessian矩阵的特征值来计算局部血管出现的概率；在Frangi方法中充分考虑了所有的特征值,并对血管特征测度给出了直观的几何解释；在Frangi的方法中,将血管增强看作是一种寻找管状几何结构的滤波过程；由于血管具有不同的直径,因此引进了一个在一定范围内变化的测量刻度；设图像为f(x,y):$H=\left[\begin{array}{cc}{f}_{\mathrm{xx}}& f_{\mathrm{xy}}\\ f_{\mathrm{yx}}& f_{\mathrm{yy}}\end{array}\right]---\left(1\right)$其中，f(x,y)表示图像内像素点(x,y)的灰度值；f_xx、f_xy、f_yx、f_yy分别表示二维图像f(x,y)的四个二阶偏微分，即像素点(x,y)与其相邻的八个像素点的灰度梯度差值；由图像内每个像素点(x,y)的灰度值的二阶梯度值构建Hessian矩阵H(x,y)；X方向上的二阶偏微分：$f_{\mathrm{xx}}=\frac{{\&PartialD;}^{2}f}{{\&PartialD;x}^2}=f(x-1,y)+f(x+1,y)-2f(x,y)---\left(2\right)$Y方向上的二阶偏微分：$f_{\mathrm{yy}}=\frac{{\&PartialD;}^{2}f}{{\&PartialD;y}^2}=f(x,y-1)+f(x,y+1)-2f(x,y)---\left(3\right)$X，Y方向上的混合偏微分：$f_{\mathrm{xy}}=f_{\mathrm{yx}}=\frac{{\&PartialD;}^{2}f}{\&PartialD;x\&PartialD;y}=f(x+1,y+1)+f(x,y)-f(x+1,y)-f(x,y+1)---\left(4\right)$由于f_xy＝f_yx，H是实对称矩阵，因此可用两个特征值λ₁、λ₂来构造增强滤波；在二维图像空间中，Hessian矩阵的两个特征值λ₁、λ₂可以由下面公式计算出：${\mathrm{\&lambda;}}_1=K-\sqrt{K^2-Q^2},---\left(5\right)$${\mathrm{\&lambda;}}_2=K-\sqrt{K^2-Q^2},---\left(6\right)$其中，K＝(f_xx+f_yy)/2，$Q=\sqrt{f_{\mathrm{xx}}{f}_{\mathrm{yy}}-f_{\mathrm{xy}}{f}_{\mathrm{yx}}};$Hessian矩阵的特征值λ₁和λ₂用于判断图像上的点是否为角点，角点是指图像中密度变化剧烈的点；定义|λ₁|≤|λ₂|，血管方向由绝对值最小的特征值λ₁所对应的特征向量给出；由于在造影图像中血管是高亮的，沿血管方向的局部灰度变化理想状况下为零，而血管的剖面方向灰度变化剧烈；因此，对于二维血管结构像素应满足如下条件：|λ₁|≈0,|λ₁|＜＜|λ₂|   (7)在冠状动脉造影拍摄过程中，由于冠状动脉管壁阻力处处不同而引起的各血管中血流量不同，以及血管重叠、肋骨和椎骨的投影成像所造成的造影图像非均匀染色等问题，在造影图像中的血管远端的细小血管分支中极为明显；同时，由于冠状动脉的直径存在变化，不适合使用单一尺度的增强效果；为增强目标小血管的对比度，消除投影所产生的非均匀染色，基于血管直径存在变化，因此引入高斯函数构造多尺度滤波器，采用不同尺度进行增强滤波，补偿各血管间的像素色度差异；将Hessian矩阵的差分运算与高斯函数结合，通过改变高斯函数的标准偏移量来获得不同尺度σ下的线形增强滤波；根据高斯函数的卷积性质，尺度空间倒数g由输入冠脉造影图像与高斯滤波器的二阶导的卷积得到：$g=I\&CircleTimes;\frac{{\&PartialD;}^{2}G(x,y;\mathrm{\&sigma;})}{\&PartialD;a\&PartialD;b}---\left(8\right)$其中，I是给定的冠脉造影图像，G是二维图像中的高斯函数，高斯函数G的表达式为：$G(x,y;\mathrm{\&sigma;})=\frac{1}{2{\mathrm{\&pi;\&sigma;}}^2}\exp (-\frac{x^2+y^2}{2{\mathrm{\&sigma;}}^2})---\left(9\right)$其中，x是水平轴上距初始点的距离；y是垂直轴上距初始点的距离；σ是高斯分布的标准偏差，为空间尺度因子；在满足式(7)的前提下，某尺度σ下的Frangi二维管状滤波公式如下所示：$R_B=\frac{{\mathrm{\&lambda;}}_2}{{\mathrm{\&lambda;}}_1}---\left(11\right)$$S={\left|\right|H\left|\right|}_F=\sqrt{{\mathrm{\&lambda;}}_1^2+{\mathrm{\&lambda;}}_2^2}---\left(12\right)$其中σ,β和c是常数因子，用于调控R_B和S的敏感度，使得满足式(7)的体素对滤波器的反馈最大；在二维管状滤波中，Frangi引入了两个参数：R_B和S；其中R_B用于识别球状结构：当体素位于球状结构中时，R_B取得最大值；同时，由于背景像素具有其导数值非常小的特征，为区分背景像素，故引入测度S去除固有噪声；冠脉造影图像的增强的步骤如下:(1)输入灰度图像，生成像素矩阵I；(2)对I的每一个元素(像素I_ij)，执行(3)～(9)；(3)初始化空间尺度σ；(4)σ若满足停止条件，则跳转(9)；(5)在当前尺度下，计算元素I_ij与高斯函数二阶微分的卷积；(6)生成Hessian矩阵H，并计算特征值λ₁和λ₂；(7)将特征值代入Frangi滤波公式式(4)并计算；(8)所有当前尺度下的元素计算完成后，迭代σ(step)，跳转(4)，计算下一尺度；(9)遍历所有尺度后，尺度迭代结束，记录每个元素在各尺度下的最大值，得到最终结果；(10)图像像素遍历结束，输出增强图像；所述步骤3对增强后的血管进行血管主干和分支的分离的方法如下：首先，对增强后的目标血管进行预处理，采用均值滤波的方法，对血管边缘进行平滑，消除下一步血管片段连接过程中血管边缘可能产生的刺突；其次，在经过图像增强后得到的冠脉二值图像，由于其背景强度低于血管强度，设定一个阈值过滤背景信息；所有信号的像素即目标冠脉树被设置为白色，其它部分为黑色；再次，依据形态学运算得到粗略的中心线；其中，血管的交叉点被认为是具有四个或更多连接的相邻像素，删除作为交叉点的像素，得到没有任何分支的血管片段；分离后的血管片段几乎能完整保留细小血管的相关信息，为得到更为完整的冠状动脉树打下基础；最后，对每个独立且无分支的血管片段进行顺序的标记；基于前三步得到的仅由“1”像素和“0”像素组成的二值图像，其相互邻接的“1”值像素目标血管像素组合成区域，且相互连通，并用边界信息描述每个连通区域即血管片段；对二值图像执行两次扫描：第一次逐行逐列执行扫描像素，判断像素间的相邻关系，对属于同一连通区域的像素赋予相同的连通标号；第二次执行扫描消除重复标记的各连通区域的子区域，合并属于同一连通区域但具有不同标记号的子区域；从而实现对每个独立且无分支的血管片段的顺序标记；所述步骤4连接分离后的冠脉血管片段的方法如下：将各无分支的血管片段端点的连接看作一个聚类问题解决；首先，需要提取各血管片段中心线端点的三个特征：(1)位置x＝(x,y)：血管片段的端点的像素坐标；(2)血管端点的方向：采用Gabor滤波器确定血管片段的方向，因为每个已标记的血管片段滤波的过程是相对独立的，所以在交叉点可能产生的干扰能够被有效地抑制；$g(x,y;\mathrm{\&theta;},\mathrm{\&psi;},\mathrm{\&sigma;},\mathrm{\&gamma;})=\exp (-\frac{x^{\&prime;2}+{\mathrm{\&gamma;}}^2{y}^{\&prime;2}}{2{\mathrm{\&delta;}}^2})\exp \left(i(2\mathrm{\&pi;}\frac{x^{\&prime;}}{\mathrm{\&lambda;}}+\mathrm{\&psi;})\right)---\left(13\right)$x′＝xcosθ+ysinθ,y′＝‑xsinθ+ycosθ式中，λ是正弦因子的波长；θ是Gabor函数法线上平行条纹的方向；ψ是相位偏移；σ是高斯函数的标准偏差；γ是空间方面比，是Gabor核高斯函数的纵横比；(3)已标记血管片段端点的血管宽度，即端点直径：端点的法线方向决定血管的宽度，从端点沿两个法线方向上的像素总和为其宽度；其次，利用测量的两个血管片段间的相似性，引入相似性函数，来决定每个血管片段末端血管的连续性；它包含三部分信息：(1)相邻血管片段端点间方向的连续性：血管片段的相邻端点具有方向上的连续性，是判断是否是同一支血管的关键因素，根据方向的连续性，可以合适匹配一支血管进入和传出的部分；同时，由于存在部分差异，很难完全准确匹配两部分的方向，需要设定一个范围，在此变化范围内认为血管方向一致；(2)端点间的距离；血管片段端点间的距离是另一个相对重要的参考因素；引入血管片段端点的距离和方位角的能量分布场；即，基于Gabor滤波器的角度范围设计一个扇形，其半径就是端点间的距离；因此，定义一个如下的衰减函数：$\mathrm{DK}(x_i,y_j)={\exp }^{-(\left(d^2(x_i,y_j)\right)/{\mathrm{\&delta;}}^2)}---\left(14\right)$式中，d²(x_i,y_i)是两点间的欧氏距离，δ是尺度参数，用于控制随距离增加的衰减速度；(3)相邻血管片段间的血管宽度的差异性：较短距离内，血管的宽度变化率很小；两支血管片段的血管宽度的差异是判断是否属于同一支血管的另一个因素；然而，上述三个因素的重要性不同，其重要性的由大到小依次是方向、距离和宽度；设定加权值，即设p和q代表端点的三维特征空间中的两个点，有相似函数如下所示：$S(x_i,y_j,\mathrm{\&epsiv;})=\left\{\begin{array}{cc}\frac{{\mathrm{\&omega;}}_1}{\mathrm{\&Delta;\&theta;}(x_i,y_j)}+{\mathrm{\&omega;}}_2\mathrm{DK}(x_i,y_j)+\frac{{\mathrm{\&omega;}}_3}{\mathrm{DW}(x_i,y_j)}& \mathrm{\&epsiv;in\; range}\\ 0& \mathrm{otherwise}\end{array}\right.---\left(15\right)$式中，ω₁，ω₂和ω₃用于调整方向、距离和宽度差异性的重要程度；x_i是第i个端点；ε是控制方向的角度合理变化范围；DK是距离衰减函数的定义；DW是血管宽度的差异性；因此，该函数能很好的处理小角度、高曲率或长距离的情况；所述步骤5基于完整的冠脉树二值图像，提取血管中心线即血管骨架的方法如下：根据步骤2、3、4所述工作的基础上，采用K主曲线法可以继续减薄已提取的较为粗略的中心线；首先，对于数据点集合K，曲线f^*称为长度L的主曲线，如果在所有长度小于或者等于L的曲线簇上，f^*最小化距离函数Δ(f)；其中，$\begin{array}{c}\mathrm{\&Delta;}\left(f\right)=E\left[\mathrm{\&Delta;}(X,f)\right]\\ =E\left[\underset{t}{\inf {\left|\right|X-f\left(t\right)\left|\right|}^2}\right]=E\left[{\left|\right|X-f\left(t_f\left(X\right)\right)\left|\right|}^2\right]\end{array}---\left(16\right)$式中，t是投影指标；f(·)表示投影指标与数据点之间的m维矢量函数，是与内在分布变量有非线性关系的描述，且只要数据分布满足有限二阶矩，则K主曲线始终存在；其次，以最短的第一主成分直线作为初始曲线f_1,n，然后迭代算法，每次迭代时都在前次获得的曲线f_i‑1,n线上增加一个新的顶点，以增加曲线的段数；在加入新顶点后，以前的顶点位置在内循环加以更新，使得惩罚距离函数最小化，从而形成新的曲线f_1,n；算法在k超过阈值时停止；内循环由投影步和最优化步组成；最后，在投影步中，数据点按照其投影的顶点或线段被分成不同区域；在最优化步，通过最小均方距离来确定每个顶点新的位置，其中均方距离函数受局部曲率惩罚；投影步与最优化步两个步骤反复执行直到收敛，曲线生成，最终得到冠状动脉树中心线即血管骨架；所述步骤7识别冠状动脉树左或右优势类型和追踪标记血管16节段的方法如下：首先，基于步骤3已标记的各血管片段和步骤5中所计算的血管方向，依据SYNTAX评分中冠脉16节段法，对已标号的“1”像素连通区域进行扫描；对肝位的冠脉造影图像从直径值最大的右冠血管片段端点的像素坐标开始扫描，对步骤4中判断为同一支的血管，从血管片段端点沿中心线方向进行追踪，对已标号的连通区域根据SYNTAX评分中冠脉16节段编号进行顺序标号，判断冠脉树是否具有后降支和后侧支，包括：4.右冠‑后降支，16.右冠‑后降支，16a.右冠‑后侧支第一分支，16b.右冠‑后侧支第二分支，和16c.右冠‑后侧支第三分支，若具有，则为右优势型冠脉树；否则，则为左优势型冠脉树；对右前斜位的冠脉造影图像进行如上方法追踪，判断冠脉树是否具有15.回旋支‑后降支，若具有，则为左优势型冠脉树；否则，则为右优势型冠脉树；所述步骤8识别血管病变的方法如下：在SYNTAX冠状动脉评分体系中，对血管狭窄、完全闭塞、三叉病变、分叉病变、开口病变、严重扭曲、血管病变长度、钙化、血栓、弥漫病变等血管病变进行评分；第一，判断血管狭窄程度并进行评分；根据步骤4获得的直径≥1.5mm的血管宽度衰减函数来判断血管是否狭窄；若某一位置血管片段横截面宽度由≥1.5mm突然衰减至0，即≥1.5mm的血管片段横截面沿中心线方向的下一个相邻的血管横截面为0，则判定该像素位置的血管为完全闭塞，标记该位置，同时设置弹窗询问其完全闭塞时间是否大于3个月，并进行打分；若某一位置血管片段横截面宽度沿该血管方向的相邻血管横截面宽度衰减率≥50％，即判定该像素位置的血管为狭窄，标记该位置，并进行打分；第二，识别分叉病变包括三分叉病变和双分叉病变，并进行评分；三分叉病变的判定包括3步：第1步，已知步骤4所得血管片段的端点像素坐标，判定是否同时有3个直径≥1.5mm血管片段端点的像素位置在同一个步骤4中所设计的基于Gabor滤波器的半径为阈值的圆内，其中包括一根主干血管片段和2根边支血管片段，并且判断3根血管片段的端点方向两两所成角度是否小于70°；第2步，若有符合上述要求的三个血管片段的端点在同一个半径为阈值的圆内，则继续判断其血管边支节段编号是否为下述边支：3/4/16/16a、5/6/11/12、11/12a/12b/13、6/7/9/9a和7/8/10/10a；第3步，若边支血管符合上述节段编号，则继续判定已标记的血管狭窄位置是否在上述半径阈值的圆内，如果存在，统计其狭窄数量如果相邻病变间距<3个参考直径，则作为一个病变计分，否则认为是两个病变，并综合第1、2步进行评分；双分叉的病变包括3步：第1步，已知步骤4所得血管片段的端点像素坐标，判定是否同时有2个直径≥1.5mm血管片段端点的像素位置在同一个步骤4中所设计的基于Gabor滤波器的半径为阈值圆内，其中包括一根主干血管片段和1根边支血管片段，并且判断2根血管片段的端点方向所呈角度是否小于70°；第2步，若有符合上述要求的2个血管片段在同一个半径为阈值的圆内，则继续判断其血管边支节段编号是否为下述边支：5/6/11、6/7/9、7/8/10、11/13/12a、13/14/14a、3/4/16和13/14/15；第3步，若边支血管符合上述节段编号，则继续判定已标记的血管狭窄位置是否在上述半径阈值的圆内，如果存在，统计其狭窄数量(如果相邻病变间距<3个参考直径，则作为一个病变计分，否则认为是两个病变)，并综合第1、2步进行评分；第三，识别开口病变，并进行评分；分别对肝位冠脉造影的第1节段和右前斜位冠脉造影的第5、6、11节段的血管片段的近端，即片段端点所在横截面至沿中心线方向3mm处像素点所在的横截面之间，是否存在50％以上的血管狭窄；如果第1和第5节段存在，则为冠脉的主动脉开口病变；如果第6和第11节段存在，则为回旋支双开口病变；并且，根据SYNTAX评分标准对其进行评分；第四，识别血管扭曲，并进行评分；首先，基于步骤5中提取的血管中心线，对每个标号的血管片段的中心线上的像素点，通过该点前后两点间距离与该点到整个连线距离比作为近似曲率，提取中心线上所有像素点的曲率；其次，通过计算出轮廓上每一点的曲率，进而分别在每个血管片段中找出各自的局部最大值点，并将其作为特征点；最后，如果某段血管存在1或多个≥90°的特征点或者3个或更多45‑90°的特征点，则该片段存在血管扭曲，并进行评分；第五，识别血管钙化点，并进行评分；首先，对肝位和右前斜位两张处理后的造影二值图像中已标记的血管片段像素位置分别对应到原始灰度图像中，并标记为相同的标号；其次，采用分段拉伸函数，设原始图像总灰度级数为M，而其中大部分的灰度分布在[a,b]区间，只有少部分的灰度在此区间外，则可以在[a,b]区间内做线性变换如下$g(x,y)=\left\{\begin{array}{cc}d& 0\&le;f(x,y)\&lt;a\\ \frac{c-d}{b-a}[f(x,y)-a]+d& a\&le;f(x,y)\&lt;b\\ c& b\&lt;f(x,y)\&lt;M\end{array}\right.---\left(17\right)$其中，g(x,y)为变换后的灰度值，[d,c]为对[a,b]段灰度进行拉伸之后的灰度值；根据冠状动脉造影的平均灰度直方图可得冠状动脉部分的灰度变化范围，其余为背景部分，灰度值为0，对原始图像进行灰度拉伸，然后采用Otsu法分割，可得到钙化点；最后，如果在该血管片段区域存在除钙化点外均为背景，则可判定该血管片段已发生钙化，并对应到分割后的二值图像中，根据钙化位置进行评分；第六，识别大于20mm的长病变，并进行评分；判断某血管片段中是否存在横截面宽度衰减≥50％的狭窄，且该狭窄沿其中心线方向长度≥20mm，若存在，则为≥20mm的长病变，根据SYNTAX评分标准进行评分；第七，识别小血管病变，并进行评分；首先，判断各血管片段中是否存在横截面宽度介于1.5mm和2mm的片段区域；其次，如果存在该区域，判断该区域内的中心线像素和在整个血管片段长中心线像素和的占比是否≥75％，如果是，则其发生了小血管病变；最后，根据SYNTAX评分标准进行评分；所述步骤9，SYNTAX自动评分的方法如下：将步骤8中对血管狭窄、完全闭塞、三叉病变、分叉病变、开口病变、严重扭曲、血管病变长度、钙化、血栓、弥漫病变等血管病变特征每项的评分结果，根据SYNTAX评分标准中冠脉树各节段所占的权重因数和病变程度自动叠加或相乘，最终给出SYNTAX评分结果，从而辅助医生确定最优的手术方式。表1为所述SYNTAX评分标准的冠脉16节段各节段的权重因数冠脉病变不良特征评分血管狭窄 ‑完全闭塞*5‑严重病变(50‑99％)*2完全闭塞(T0) ‑时间大于3个月或未知+1‑钝型残端+1‑自身桥侧枝+1‑完全闭塞病变以及第一个可见节段+1/每个不可见部分‑边支‑边支小于1.5mm+1三叉病变 ‑1个病变节段+3‑2个病变节段+4‑3个病变节段+5‑4个病变节段+6分叉病变 ‑A、B、C型病变+1‑E、D、F、G型病变+2‑角度小于70度+1开口病变+1严重扭曲+2长度大于20mm+1严重钙化+2血栓+1弥漫病变/小血管病变+1/每个段号表2为所述SYNTAX评分标准的病变不良特征评分冠脉病变不良特征评分血管狭窄 ‑完全闭塞*5‑严重病变(50‑99％)*2完全闭塞(T0) ‑时间大于3个月或未知+1‑钝型残端+1‑自身桥侧枝+1‑完全闭塞病变以及第一个可见节段+1/每个不可见部分‑边支‑边支小于1.5mm+1三叉病变 ‑1个病变节段+3‑2个病变节段+4‑3个病变节段+5‑4个病变节段+6分叉病变 ‑A、B、C型病变+1‑E、D、F、G型病变+2‑角度小于70度+1开口病变+1严重扭曲+2长度大于20mm+1严重钙化+2血栓+1弥漫病变/小血管病变+1/每个段号