基于造影影像及冠状动脉的血流速度计算方法[发明专利]

(12)发明专利申请

(10)申请公布号 CN 110448319 A(43)申请公布日 2019.11.15

(21)申请号 201810430295.4(22)申请日 2018.05.08

(71)申请人 博动医学影像科技（上海）有限公司

地址 200233 上海市闵行区桂平路555号46

栋202室(72)发明人 李莹光　张义敏　刘新峰　梁夫友　

涂圣贤　(74)专利代理机构 北京迎硕知识产权代理事务

所(普通合伙) 11512

代理人 吕良　张群峰(51)Int.Cl.

A61B 6/00(2006.01)

权利要求书1页 说明书6页 附图3页

()发明名称

基于造影影像及冠状动脉的血流速度计算方法

(57)摘要

一种基于造影影像的冠状动脉血流速度计算方法，该方法包括以下步骤：获取冠状动脉造影图像序列；对获取的每一帧冠状动脉造影图像进行分割和二值化处理，提取出冠脉树；从冠脉树中提取目标血管中心线；根据不同帧的目标血管中心线长度的变化情况，拟合计算出充盈阶段的目标血管平均血流速度。

CN 110448319 ACN 110448319 A

权　利　要　求　书

1/1页

1.一种基于造影影像的血流速度计算方法，其特征在于，该方法包括以下步骤：获取目标血管的造影图像序列；

对获取的每一帧血管造影图像进行分割处理，提取出目标血管形状；提取目标血管的中心线；

根据不同图像帧中目标血管中心线长度的变化，拟合计算出血管造影充盈阶段的平均血流速度。

2.根据权利要求1所述的基于造影影像的血流速度计算方法，其特征在于，所述的分割处理，包括先对血管进行特征提取，然后再进行二值化，分割后的图像进行降噪二次处理。

3.一种基于造影影像的冠状动脉血流速度计算方法，其特征在于，该方法包括以下步骤：

获取冠状动脉造影图像序列；

对获取的每一帧冠状动脉造影图像进行分割处理，提取出冠状动脉血管树；从冠状动脉血管树中提取目标血管中心线；

根据不同图像帧的目标血管中心线长度的变化，拟合计算出血管造影充盈阶段的目标血管平均血流速度。

4.根据权利要求3所述的基于造影影像的冠状动脉血流速度计算方法，其特征在于，在对冠状动脉造影图像进行分割前要进行降噪预处理。

5.根据权利要求3所述的基于造影影像的冠状动脉血流速度计算方法，其特征在于，对冠状动脉造影图像分割时，先对冠状动脉血管进行特征提取，然后再进行二值化，分割后的图像进行降噪二次处理。

6.根据权利要求5所述的基于造影影像的冠状动脉血流速度计算方法，其特征在于，在对目标血管进行特征提取时，采用多尺度Hessian矩阵的方式。

7.根据权利要求5所述的基于造影影像的冠状动脉血流速度计算方法，其特征在于，在对目标血管进行特征提取时，采用多尺度Gabor函数滤波器的方式。

8.根据权利要求5所述的基于造影影像的冠状动脉血流速度计算方法，其特征在于，在对冠状动脉造影图像进行冠状动脉血管树进行目标血管特征提取和二值化以后，采取追踪的方式提取目标血管中心线。

9.根据权利要求8所述的基于造影影像的冠状动脉血流速度计算方法，其特征在于，具体步骤包括：

统计每一帧冠状动脉造影图像的目标血管中心线的长度，根据图像中目标血管中心线由多个离散的点构成，计算每一帧图像中目标血管中心线的长度；

得到帧数-中心线长度变化关系后，自动选取多个与血流充盈过程对应的中心线长度点做线性拟合；

通过直线斜率计算出目标血管平均血流速度。

10.根据权利要求3所述的基于造影影像的冠状动脉血流速度计算方法，其特征在于，若关注的血管不是视野中出现的血管时，则确定一个感兴趣区域，通过这个感兴趣区域的圈定，提取出该感兴趣区域内的目标血管中心线。

2

CN 110448319 A

说　明　书

基于造影影像及冠状动脉的血流速度计算方法

1/6页

技术领域

[0001]本发明涉及医疗领域，尤其涉及应用在基于造影影像的血流速度计算方法，以及在冠状动脉血流速度计算中的使用。

背景技术

[0002]近年来，随着人口老龄化的加速，中国心血管病危险因素流行趋势明显，导致了心血管病的发病人数持续增加，已成为中国城乡居民死亡的首要原因，也加重了社会和家庭的负担，成为重大的公共安全问题。心血管疾病主要是缺血性心脏病，通常指冠状动脉粥样硬化性心脏病或冠心病。

[0003]常规的心血管成像技术，比如X射线冠脉造影，可以显示血管的损伤部位与损伤范围，但是对于确定狭窄是否造成缺血以及定位犯罪血管存在一定的局限性。血流储备分数(FFR)可以弥补影像技术的不足，通过测量冠脉狭窄远端与狭窄近端的压力比值来反应狭窄病变对最大血流量的程度，以此判断狭窄部位是否引起缺血与是否需要进行支架术来重建血运。目前，FFR已成为临床上评价缺血性心脏病的金标准，被欧洲心脏病学会(ESC)指南规定为IA级临床证据与被美国心脏病学会(ACC指南规定为IIa级临床依据。[0004]尽管压力导丝引导的血流储备分数(FFR)测量是检测冠心病的金标准，但由于其局限性，在临床推广依旧不是很理想。这些局限因素包括，FFR检测是有创的，且有些病人在注入腺苷时反应不适，FFR检测的价格也比较高，增加了医疗成本等。

[0005]而定量血流分数(QFR)是一种新型的用于评估冠脉狭窄功能学意义的方法，由冠状动脉造影的三维重建结果和流体力学方法得到。FAVOR Pilot研究显示核心实验室的QFR分析能很好地评估冠脉狭窄的功能学意义，FAVOR II China研究显示导管室在线实时QFR分析也有很好的可行性和准确性。

[0006]现有的血流计算方法主要还包括多普勒导丝法和温度稀释法。多普勒导丝可以持续测量血管内压力和速度，是测量速度的金标准。但这种测量方式是有创的，测量值也会因为时刻和导丝位置的改变而改变，重复性较差，且导丝造价高昂，也增加了医疗成本。温度稀释法是用冷生理盐水作为指示剂，具有热敏电阻的漂浮导管作为心导管，检测血流温度的变化情况，求得血液温度随时间变化的曲线，血流速度与指示剂的平均转运时间成反比。这种方法仍然是有创的，会受冷水注入速度影响，得到的也是温度随时间变化的曲线，并没有直接得到速度值，导丝费用高昂也是一个问题。

[0007]现有技术公开了一种血管血流速度的计算方法，该方法包括：确定血管的感兴趣区域；计算并拟合出所述感兴趣区域内的灰度拟合曲线；确定一预定时间区间内的最大灰度值曲线或最小灰度值曲线；计算预定时间区间内最大灰度值曲线或最小灰度值曲线与灰度拟合曲线所包围的区域面积值；基于该区域面积值，获取该面积值所对应的单位时间血流量；基于所述单位时间血流量以及血管管腔面积，获得所述血管的血流速度。优选的，所述感兴趣区域包括注入造影剂的主支血管及其分支。优选的，通过目标图像跟踪检测感兴趣区域在不同心跳时刻下的位置的变化，从而获得最佳的感兴趣区域。优选的，所述方法进

3

CN 110448319 A

说　明　书

2/6页

一步包括：接受血管的X射线造影影像序列，选择感兴趣区域；选择起始时间为造影剂充盈之前，提取每帧造影中感兴趣区域内灰度直方图，通过所述灰度直方图计算出每帧下感兴趣区域内的灰度值，并根据灰度值拟合出灰度随时间变化的灰度拟合曲线。优选的，所述方法进一步包括：确定第一时间点，以及以第一时间点为中心的一预定时间区间内灰度拟合曲线的最大值和最小值。该技术方案，由于获取的造影图像收到心脏跳动等的因素影响，图像噪声不容易抑制，造成计算结果不够稳定。发明内容

[0008]本发明的实施例的目的是提供一种基于造影影像的血流速度计算方法，以及基于造影影像的冠状动脉血流速度计算方法，以优化对冠状动脉血流速度的计算。[0009]本发明的实施例之一，一种基于造影影像的血流速度计算方法，该方法包括以下步骤：

[0010]获取目标血管的造影图像序列；

[0011]对获取的每一帧血管造影图像进行分割处理，提取出目标血管形状；[0012]提取目标血管的中心线；

[0013]根据不同图像帧中目标血管中心线长度的变化，拟合计算出血管造影充盈阶段的平均血流速度。

[0014]本发明的实施例之一，一种基于造影影像的冠状动脉血流速度计算方法，该方法包括以下步骤：

[0015]获取冠状动脉造影图像序列；

[0016]对获取的每一帧冠状动脉造影图像进行分割和二值化处理，提取出冠状动脉树；[0017]从冠状动脉树中提取目标血管中心线；

[0018]根据不同帧的目标血管中心线长度的变化情况，拟合计算出充盈阶段的目标血管平均血流速度。冠状动脉树即是冠状动脉血管树。

[0019]根据本发明实施例通过对于造影图像序列的处理，对于血管血流速度计算，相对于现有技术，获得的有益效果包括：[0020](1)全程自动处理，减少了人为因素的影响，具有很强的可重复性；[0021](2)该方法操作简便，便于医生学习使用；[0022](3)相对于现有技术，计算方法更加合理，计算精度高；[0023](4)对于造影图像中的血管评估均可使用，应用更加广泛。

[0024]本发明实施例提供的基于造影影像的冠状动脉血流速度的定量计算方法，弥补现有血流计算方法的不足，实现对冠状动脉血流的定量评估，包括但不限于将计算结果用于QFR值的计算。

附图说明

[0025]通过参考附图阅读下文的详细描述，本发明示例性实施方式的上述以及其他目的、特征和优点将变得易于理解。在附图中，以示例性而非性的方式示出了本发明的若干实施方式，其中：

[0026]图1本发明实施例中多尺度Hessian矩阵对血管特征提取效果图。

4

CN 110448319 A[0027][0028][0029][0030][0031][0032]

说　明　书

3/6页

图2本发明实施例中多尺度Gabor滤波器对血管特征提取效果图。图3本发明实施例中用Fast Marching算法提取中心线效果图。图4本发明实施例中计算出连续图像序列的中心线变化示意图。图5本发明实施例中中心线坐标示意图。图6本发明实施例中血流速度拟合示意图。图7本发明实施例中感兴趣区域示意图。

具体实施方式

[0033]本发明实施例中的冠脉动脉血管说明如下：[0034]冠状动脉分为左右两支，其开口分别位于左、右冠状窦内，左冠状动脉又分为左主干(LM)、左前降支(LAD)、回旋支(LCX)。左主干为一自主动脉根部发出的长约1～3cm的主干，左前降支走行于室间沟内，回旋支走行于左房室沟内。右冠状动脉(RCA)主干走行于右房室沟内。

[0035]根据一个或多个实施例，一种基于造影影像的冠状动脉血流速度计算方法，该方法包括以下步骤：通过冠状动脉造影技术获取冠状动脉造影图像序列。该图像序列反映了冠状动脉树的血流灌注情况。对获取的每一帧冠状动脉造影图像进行分割处理，提取出冠状动脉血管树。从冠状动脉血管树中提取目标血管中心线。目标血管可以是冠脉树的主支或者分支。主支可能是左前降支(LAD)、回旋支(LCX)或者右冠状动脉(RCA)。图1、图2或图3中，分别体现了上述主支的形态。

[0036]根据不同图像帧的目标血管中心线长度的变化，拟合计算出血管造影充盈阶段的目标血管平均血流速度。

[0037]在对冠状动脉造影图像进行分割前要进行降噪预处理。降噪预处理包括抑制背景噪声，突出血管结构，增强前景和背景的对比度，同时保证血管边缘清晰。所述背景噪声包括，随着心脏跳动而移动的膈肌噪声和不随心脏跳动而移动的脊柱噪声。[0038]对冠状动脉造影图像分割时，先对冠状动脉血管进行特征提取，然后再进行二值化，分割后的图像进行降噪二次处理。该降噪二次处理的目的是突出血管结构，减少背景噪声，保证血管分割结果的连续性。[0039]根据一个或多个实施例，在对血管进行特征提取时，采用多尺度Hessian矩阵的方式。Hessian矩阵被广泛用于检测和分析特定形状，在点P＝(x,y)处，该曲面的曲率可以用Hessian矩阵表示为：

[0040][0041]

假设两个特征值|λ1|<|λ2|，λ1和其特征向量v1表示曲率小的强度和方向，λ2和其

特征向量v1表示曲率大的强度和方向，v1平行与血管，v2垂直与血管轴，血管点的特征是λ1≈0，λ2≥0，针对二维图像中斑状与线状几何结构，相应的Hessian矩阵特征值大小及符号不同，可构建血管增强滤波器。

[0042]Frangi等定义二维的血管相似性相应函数为

5

CN 110448319 A

说　明　书

4/6页

[0043]

[0044]式中，RB＝λλ1/2用于在二维图像中区分球状结构和管状结构；

为Hessian矩阵的F范数；β和c分别为控制线性滤波器在RB和S下敏感度的比例因子。

[0045]多尺度Hessian矩阵滤波时，由图像和一定尺度下高斯滤波器g(P,σ)的二阶偏微分卷积，可以求得图像的二阶导数，即：[0046]Ix，σ)＝gx，σ)*I(P)y(P，y(P，[0047]其中，高斯函数表示为：

[0048][0049]

高斯函数作为卷积核是实现尺度变换的唯一线性核，使用高斯滤波器与冠脉图像进行卷积滤波，可以平滑图像去除噪声。通过改变σ值可以实现多尺度的计算，将不同尺度下的最大相应值作为输出值，实现血管特征增强结果。在实际计算中，需要确定血管的不同尺度大小，确定Frangi函数中的两个常数λλ即可得到血管增强后的效果。如图1所示多1和2，尺度Hessian矩阵对血管特征提取效果图。[0050]可选择的，可以采用形态学Bottem-Hat变换来进一步强调血管。Bottem-Hat变换是一种常见的用于提取图像暗结构的形态学算子，定义如下：[0051]g＝f-(f·b)

[0052]

其中，f是输入图像，b是结构元素函数，代表闭操作，由膨胀操

作和腐蚀操作组成，闭操作可以去除由结构元素b表示的图像中较暗的感兴趣结构，并保留

较亮像素。在造影图像中，目标物体(血管)的灰度值较背景低，因此，f·b可以看作是背景，由原图减去背景图像可以得到血管部位。[0053]根据一个或多个实施例，在对血管进行特征提取时，采用多尺度Gabor函数滤波器的方式。Gabor函数实际上是由复数正弦函数调制而得到的高斯函数，其表达式如下：[00]h(x,y)＝g(x′，y′)exp(j2πFx′)

[0055]

其中(x′，y′)＝(x cosθ+y sinθ，-x sinθ+y cosθ)，θ为滤波器的方位角，可以通过在x-y平面旋转获得任意的θ。F为滤波器的中心频率，即滤波器的带通中心所在的频域中的位置。σσ分别表示x方向和y方向的常量值，并由σx、y是高斯包络的空间常量，x由频率带宽确定，σy由方位带宽所确定：

[0056][0057]

6

CN 110448319 A

说　明　书

5/6页

[0058]

其中，BF是频率带宽，它体现了空域及频域上滤波器的局部变化大小，Bθ是方位角

带宽，它体现了滤波器对不同方位角的敏感程度。如图2所示多尺度Gabor滤波器对血管特征提取效果图。

[0060]根据一个或多个实施例，在对冠状动脉造影图像进行冠状动脉血管树进行血管特征提取和二值化以后，采取追踪的方式提取主支血管中心线。[0061]中心线的提取方法主要有：基于细化/骨骼化的方法；基于追踪的方法。[0062]基于细化的中心线提取算法的基本思想是：先分割出目标结构，再对分割得到的结果进行细化操作，从而得到目标结构的中心线(骨架)。血管中心线可以看作血管的“骨架”，即血管在其真实中轴附近的一个近似。在满足拓补不变和几何约束条件的前提下，血管被细化算法层层剥离，最终成为在中心线附近的单像素细线，在视觉上仍保留原有血管树拓补结构。[0063]这里，所谓细化是把一个具有一定面积的区域用一条(或一组)曲线(或细线)来表示它。从广义角度来说，细化操作属于连接成份的变形操作。如果将连接的成份(集合)用符号“1”来表示，背景用符号“0”来表示。则细化操作就是用改变连接成份的形状使符号“1”中某些像素由“1”变成“0”，迭代这一过程，直到最后由一组单个像素组成的一组曲线或细线来代表这个区域。这组曲线(或细线)应保留连接成份的连通性和轮廓的几何特征。[00]基于追踪的算法的基本思想是：首先给定初始点，在通过确定其与邻域点之间的关系，自动迭代跟踪出中心线上所有的点集。跟踪通常从初始点开始探测血管中心线或通过分析正交于跟踪方向的像素来检测边缘。实现探测血管中心线或轮廓的方法有很多，最直接的便是先进行边缘检测，而后利用边缘连接信息引导跟踪算法进行跟踪。[0065]当采取追踪的方式提取血管中心线，可用Fast Marching算法提取中心线。Eikonal方程：

T＝0 on Γ描述了闭曲线在法向速度F(x,y)下的演化，速度函数F

仅依赖位置，求解这个方程可以得到曲线演化到一个点的时间，从而得到图像的时间图。通过快速行进算法可以对它进行求解。如图3所示Fast Marching算法提取中心线效果图。在图3中，我们找到了冠状动脉树的中心线，但是我们关注的是前降支的情况，所以我们需要确定一个起始点，寻找出最长中心线。[0066]根据一个或多个实施例，通过图像分割和中心线提取，我们能计算出连续图像序列的中心线变化情况，中心线由箭头标出，如图4所示。图4-1，图4-2，图4-3，图4-4，图4-5，图4-6表示出一个血管充盈过程中心线的变化。图中的箭头表示血管中心线。[0067]根据一个或多个实施例，在提取每一图像帧的主支血管中心线，统计每一帧冠状动脉造影图像的主支血管中心线的长度，统计方法如图5所示。图像中主支血管中心线由多个离散的点构成，假设两个点x1坐标(m1,n1),x2坐标(m2,n2),则这两个点之间的距离

所以中心线的长度，即x1到xn之间的长度为

由

[0059]

此可以得到每帧中心线的长度。

[0068]计算每一帧图像中主支血管中心线的长度，可以得到如图6-1所示的帧数-中心线

7

CN 110448319 A

说　明　书

6/6页

长度变化关系，自动选取多个与血流充盈过程对应的中心线长度点做线性拟合，通过直线斜率计算出主支血管平均血流速度，如图6-2所示。关注血流充盈过程，即对中心线长度均匀增长的阶段，进行线性拟合。[0069]根据一个或多个实施例，如图7所示。对于冠脉血管来说，随着心脏的跳动，冠脉血管会随着产生移动和变形。当关注的血管不是视野中最长的血管时，就需要确定一个感兴趣区域，通过这个感兴趣区域的圈定，我们可以找到在该区域内的最长血管中心线，也就是主支血管中心线。

[0070]根据一个或多个实施例，基于造影影像的冠状动脉血流速度计算系统，包括图像处理模块、中心线提取模块和速度计算模块。系统对冠脉造影图像序列进行分割和二值化，提取冠脉树，然后从冠脉树中提取主支中心线，根据不同帧血管中心线长度的变化情况，可以拟合出充盈阶段的平均血流速度。

[0071]图像处理模块的功能是对图像预处理和图像分割。血管造影图像背景复杂，并且有膈肌和脊椎图像的影像，需要再分割前对血管进行合适的预处理，抑制背景噪声，突出血管结构，增强前景和背景的对比度，同时还要保证血管边缘清晰。图像分割时可以先对血管进行特征提取，然后再进行二值化，需要尽可能突出血管结构，尽可能减少背景噪声的分割，且需要保证血管分割结果的连续性。[0072]中心线提取模块：在对冠脉造影进行特征提取和二值化以后，采取追踪的方式提取血管中心线。可用Fast Marching算法提取中心线。Eikonal方程：

T＝0 on Γ

描述了闭曲线在法向速度F(x,y)下的演化，速度函数F仅依赖位置，求解这个方程可以得到曲线演化到一个点的时间，从而得到图像的时间图。通过快速行进算法可以对它进行求解。在图3中，我们找到了冠脉树的中心线，但是我们关注的是前降支的情况，所以我们需要确定一个起始点，寻找出最长中心线。[0073]速度计算模块：通过图像分割和中心线提取，计算出连续序列的中心线变化情况，中心线由箭头标出，如图4所示。

[0074]统计每一帧的中心线的长度，统计方法如图5所示：中心线由多个离散的点构成，假设两个点x1坐标(m1,n1),x2坐标(m2,n2),则这两个点之间的距离所以中心线的长度，即x1到xn之间的长度为

[0075]

由此可以得到每帧中心线的长度。

计算每一帧中心线的长度，可以得到如图6-1的帧数-中心线长度变化图，关注血流充盈过程，即中心线长度均匀增长的阶段，进行线性拟合，通过直线斜率可以计算出平均血流速度。如图6-2所示是线性拟合图。[0076]值得说明的是，虽然前述内容已经参考若干具体实施方式描述了本发明创造的精神和原理，但是应该理解，本发明并不限于所公开的具体实施方式，对各方面的划分也不意味着这些方面中的特征不能组合，这种划分仅是为了表述的方便。本发明旨在涵盖所附权利要求的精神和范围内所包括的各种修改和等同布置。

8

CN 110448319 A

说　明　书　附　图

1/3页

图1

图2

图3

9

CN 110448319 A

说　明　书　附　图

2/3页

图4

图5

图6

10

CN 110448319 A

说　明　书　附　图

3/3页

图7

11