CN106023202B - 基于心脏ct图像的冠状动脉血流储备分数计算方法 - Google Patents

Abstract

本发明公开了一种基于心脏CT图像的冠状动脉血流储备分数计算方法，包括通过提取心肌图像，提取冠状动脉口，对冠状动脉精确分割；通过对冠状动脉体数据的边缘检测，生成FFR计算所需要的三角网格模型；最后通过公式计算血流储备分数，其中，P_aorta是主动脉平均压，(P_outlet)_i是三维网格图像中点i的压力值。本发明基于心脏CT图像，能快速准确的全自动得到血流储备分数FFR。

Description

基于心脏CT图像的冠状动脉血流储备分数计算方法

技术领域

本发明涉及冠状动脉影像学评价领域，具体地涉及一种基于心脏CT图像的全自动冠状动脉血流储备分数计算方法。

背景技术

冠状动脉造影及血管内超声均被认为是诊断冠心病的“金标准”，但它们只能对病变狭窄程度进行影像学评价，而狭窄到底对远端血流产生了多大影响却不得而知；血流储备分数(FFR)现已经成为冠脉狭窄功能性评价的公认指标，其最重要的功能是对一个未知影响的冠脉狭窄的功能后果进行准确评价。

血流储备分数(FFR)是指在冠状动脉存在狭窄病变的情况下，目标测量血管所供心肌区域能获得的最大血流量与同一区域理论上正常情况下所能获得的最大血流量之比。FFR主要通过计算冠状动脉狭窄远端压力与主动脉根部压力之比来获得。狭窄远端压力可以通过压力导丝在最大灌注血流(通过冠脉内或静脉内注射罂粟碱或腺苷或ATP)时测得。

FFR＝Pd/Pa(Pd为指引导管测量的主动脉压，Pa为压力导丝测量的冠脉狭窄远端压力)一般来说，FFR是指最大充血状态下，不存在“静息FFR”这个概念。

正常的心外膜冠状动脉对血流的阻力很小，FFR的正常值为1.0；FFR的值将小于1.0时即表明当前心外膜冠脉有狭窄病变的存在。

FFR＜0.75的情况时，所代表的狭窄情况几乎都会导致心肌缺血，FFR≥0.75的情况时，所代表的狭窄则造成心肌缺血的可能性非常小。

冠脉CTA能准确评估冠脉狭窄程度，且能辨别管壁斑块性质，是一种无创、操作简单的诊断冠状动脉病变检查方法，可作为筛查高危人群的首选方法。因此，如果对于冠心病患者的血管进行干预，前期应该对患者冠脉进行CTA的评价。冠状动脉慢性完全闭塞病变(CTO)如果采用CTA进行评价，评价结果肯定有一些有价值的信息。

通过冠状动脉CT血管造影CCTA计算无创获得的FFR(CTFFR)不仅无需额外影像检查或药物，与造影时测定的FFR相关性良好，这种一体化技术能从根本上避免不必要的冠脉血管造影与血运重建治疗。DeFacto试验结果也清楚地表明，在冠状动脉CT中，CTFFR结果的分析提供了那些真正限制血流及增加病人危险性的病变的生理信息。CTFFR结合了冠脉CTA和FFR的优势，可以从结构和功能两方面来评估冠状动脉狭窄，成为一种提供冠脉病变解剖学和功能学信息的崭新无创性检测体系。

发明内容

为了解决上述的技术问题，本发明目的是：提供一种全新的基于心脏CT图像的冠状动脉血流储备分数计算方法，基于心脏CT图像，能快速准确的全自动得到血流储备分数FFR。

本发明的技术方案是：

一种基于心脏CT图像的冠状动脉血流储备分数计算方法，其特征在于，包括以下步骤：

S01：提取心肌图像：对心脏CT图像进行分割，通过形态学操作得到心脏图像，对该心脏图像进行直方图分析得到心室心房图像，通过心脏图像与心室心房图像做差得到心肌图像；

S02：提取冠状动脉口：对主动脉图像的二值化图像进行形态学膨胀，得到全主动脉的二值图像，并通过像素取反得到全主动脉互补图像，根据主动脉中心线上点的平均灰度进行区域生长，得到含有冠状动脉口的主动脉图像，用含有冠状动脉口的主动脉图像与全主动脉互补图像做图像乘法，得到含有冠状动脉口的图像，并确定冠状动脉口；

S03：提取冠状动脉：在心肌图像上以冠状动脉口为种子点，通过区域生长提取冠状动脉，计算冠状动脉的平均灰度和平均方差，根据冠脉灰度分布，沿着冠状动脉方向提取冠脉树；

S04：生成冠状动脉网格模型：将冠状动脉图像进行二值化，然后绘制等值面图像，得到冠状动脉三维网格图像；

S05：冠状动脉FFR计算：通过公式计算血流储备分数，其中，P_aorta是主动脉平均压，(P_outlet)_i是三维网格图像中点i的压力值。

优选的，所述步骤S02中确定冠状动脉口通过如下步骤得到：对含有冠状动脉口的图像进行连通域分析，用不同的灰度标签标识各个连通域，确定冠状动脉口。

优选的，所述步骤S02中的主动脉图像为在心脏图像上，利用主动脉截面成圆形的特征，提取升主动脉及中心线，得到主动脉图像。

优选的，所述步骤S04中冠状动脉图像二值化包括：遍历冠状动脉图像V1中体素，如果像素A1等于0，则该像素值不变；如果A1不等于0，则将A1的像素值设为1，得到一个新的数据V2。

与现有技术相比，本发明的优点是：

基于心脏CT图像，能快速准确的全自动得到血流储备分数FFR。

附图说明

下面结合附图及实施例对本发明作进一步描述：

图1为心脏CT图像；

图2为心脏图像；

图3为主动脉图像；

图4为本发明的基于心脏CT图像的冠状动脉血流储备分数(FFR)计算方法的流程图；

图5为心肌图像；

图6为含有冠状动脉口的主动脉图像；

图7为分割的冠状动脉口；

图8为冠状动脉图像；

图9为冠状动脉的三维网格图像。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明了，下面结合具体实施方式并参照附图，对本发明进一步详细说明。应该理解，这些描述只是示例性的，而并非要限制本发明的范围。此外，在以下说明中，省略了对公知结构和技术的描述，以避免不必要地混淆本发明的概念。

给定心脏CT图像如图1，根据逆向方法，提取心脏：以非目标区域的降主动脉、脊椎、肋骨为对象进行处理，通过逐步的去除胸腔壁、肺部、椎骨和降主动脉等非心脏组织来提取得到心脏图像，如图2所示。

在得到的心脏图像上，通过利用主动脉截面成圆形的特征，提取升主动脉及中心线，得到主动脉图像，如图3所示。

如图4所示，本发明的基于心脏CT图像的冠状动脉血流储备分数(FFR)计算方法包括提取心肌图像、提取冠状动脉口、提取冠状动脉、生成冠状动脉网格模型、冠状动脉FFR计算。

1：提取心肌图像：

对心脏CT图像进行分割，通过形态学操作得到心脏图像，对该心脏图像进行直方图分析得到心室心房图像，通过心脏图像与心室心房图像做差得到心肌图像，如图5所示。

2：提取冠状动脉口：

对主动脉图像的二值化图像进行形态学膨胀，得到全主动脉的二值图像，并通过像素取反得到全主动脉互补图像。

根据主动脉中心线上点的平均灰度进行区域生长，得到含有冠状动脉口的主动脉图像，如图6所示。

用含有冠状动脉口的主动脉图像与全主动脉互补图像做图像乘法，得到含有冠状动脉口的图像，对含有冠状动脉口的图像进行连通域分析，用不同的灰度标签标识各个连通域，确定冠状动脉口，如图7所示。

3：提取冠状动脉：

在心肌图像上，以冠状动脉口为种子点，通过区域生长提取冠状动脉，计算冠状动脉的平均灰度和平均方差，根据冠脉灰度分布，沿着冠状动脉方向提取冠脉树，如图8所示。

4：生成冠状动脉网格模型：

通过步骤三，得到冠状动脉图像数据V1，该数据中的体素在空间上构成一个立方体，属于冠状动脉部分的体素像素值不为0(像素值大约在-3000到3000之间)，其余体素像素值都为0。

本步骤需要把数据变成空间三维网格数据V3，以便于步骤五中的FFR计算。

(1)冠状动脉数据二值化

遍历冠状动脉图像数据V1中体素，做简单的像素值判断，如果像素A1等于0，则该像素值不变；如果A1不等于0，则将A1的像素值设为1。

最终会得到一个新的图像数据V2，该图像中，属于冠状动脉部分的体素像素值为1，其余部分为0。

(2)等值面生成

本方法中，等值面生成使用Marching Cubes算法。

该算法中，体素被定义为一个极小的六面体，相邻上下层之间的四个像素组成的立方体上的八个顶点。而等值面就是在空间中所以具有某个相同属性值的点的集合。它可以表示成：

{(x,y,z)│f(x,y,z)＝c}，c是常数 [A1]

本方法中的c是在三维重构过程中给定的像素值1。

Marching Cubes算法提取等值面的流程如下：

(1)将原始数据经过预处理之后，读入特定的数组中；

(2)从网格数据体中提取一个单元体成为当前单元体，同时获取该单元体的所有信息；

(3)将当前单元体8个顶点的函数值与给定等值面值C进行比较，得到该单元体的状态表；

(4)根据当前单元体的状态表索引，找出与等值面相交的单元体棱边，并采用线性插值的方法计算出各个交点的位置坐标；

(5)利用中心差分法求出当前单元体8个顶点的法向量，再采用线性插值的方法得到三角面片各个顶点的法向；

(6)根据各个三角面片顶点的坐标和顶点法向量进行等值面图象的绘制。

最终得到冠状动脉的三维网格图像数据V3，如图9所示。

5：冠状动脉FFR计算：

1)三维模型

基于CT重构的几何结构，划分三维网格，用数值法(如：有限差分、有限元、有限体积法等)求解连续性和Navier-Stokes方程：

其中，P，ρ，μ分别为流速、压力、血流密度、血流粘性，表示梯度，t表示单位时间，T表示转置矩阵。

入口边界条件为P_aorta(主动脉平均压，来自于患者的测量结果)-13毫米汞柱(静脉注入腺苷后的压力降)-P₀(零流压，51毫米汞柱)。

出口阻力边界条件的推导如下：

其中，Q_total为心肌体积乘以心肌血流密度，(Q_outlet)_i为出口i的血流，D_i为出口i的直径，(R_outlet)_i为出口i的阻力，冠状动脉血流储备CFR设为2.7。设想冠状动脉血管正常情况下从入口到出口有1毫米汞柱的压力降。算出三维网格中各点的压力(P_outlet)_i后，通过公式A6得到FFR：

其中，P_aorta是主动脉平均压，(P_outlet)_i是三维网格图像中各点的压力值。

2)二维模型

基于CT重构的几何结构，把有狭窄的血管拉直(二维轴对称模型)，划分二维网格，用数值法(如：有限差分、有限元、有限体积法等)求解连续性和Navier-Stokes方程：

其中，ρ表示血液的密度，u_z、u_r分别表示z向、r方向的流速，μ表示血液的动力粘度，p表示血液的压强。

入口边界条件为P_aorta(主动脉平均压，来自于患者的测量结果)-13毫米汞柱(静脉注入腺苷后的压力降)-P₀(零流压，51毫米汞柱)。

出口阻力边界条件的推导如下：

其中，Q_total为心肌体积乘以心肌血流密度，(Q_outlet)_i为出口i的血流，D_i为出口i的直径，(R_outlet)_i为出口i的阻力，冠状动脉血流储备CFR设为2.7。

针对各种类型的弯曲，用三维模型计算从入口到出口的压力差，对照二维轴对称模型计算，建立用于存储各种类型的弯曲对二维轴对称结果的修正系数的数据库。这样算出压力后对照数据库中的修正系数，精确得到从入口到出口的压力差，然后计算FFR。

应当理解的是，本发明的上述具体实施方式仅仅用于示例性说明或解释本发明的原理，而不构成对本发明的限制。因此，在不偏离本发明的精神和范围的情况下所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。此外，本发明所附权利要求旨在涵盖落入所附权利要求范围和边界、或者这种范围和边界的等同形式内的全部变化和修改例。

Claims (5)

1.一种基于心脏CT图像的冠状动脉血流储备分数计算方法，其特征在于，包括以下步骤：

S01：提取心肌图像：对心脏CT图像进行分割，通过形态学操作得到心脏图像，对该心脏图像进行直方图分析得到心室心房图像，通过心脏图像与心室心房图像做差得到心肌图像；

S02：提取冠状动脉口：对主动脉图像的二值化图像进行形态学膨胀，得到全主动脉的二值图像，并通过像素取反得到全主动脉互补图像，根据主动脉中心线上点的平均灰度进行区域生长，得到含有冠状动脉口的主动脉图像，用含有冠状动脉口的主动脉图像与全主动脉互补图像做图像乘法，得到含有冠状动脉口的图像，并确定冠状动脉口；

S03：提取冠状动脉：在心肌图像上以冠状动脉口为种子点，通过区域生长提取冠状动脉，计算冠状动脉的平均灰度和平均方差，根据冠脉灰度分布，沿着冠状动脉方向提取冠脉树；

S04：生成冠状动脉网格模型：将冠状动脉图像进行二值化，然后绘制等值面图像，得到冠状动脉三维网格图像；

S05：冠状动脉FFR计算：用数值法求解连续性和Navier-Stokes方程：

$\▿\·\stackrel{\→}{V}=0---\[A2\]$

$\mathrm{\ρ}\frac{\∂\stackrel{\→}{V}}{\∂t}+\mathrm{\ρ}\stackrel{\→}{V}\·\▿\stackrel{\→}{V}=-\▿P+\▿\·\mathrm{\μ}(\▿\stackrel{\→}{V}+{\left(\▿\stackrel{\→}{V}\right)}^T)---\[A3\]$

其中，P，ρ，μ分别为流速、压力、血流密度、血流粘性，▽表示梯度，t表示单位时间，T表示矩阵的转置；

入口边界条件为：P_aorta-13毫米汞柱-P₀，其中，P₀为零流压；

出口阻力边界条件的推导如下：

${\left(Q_{outlet}\right)}_i=Q_{total}\×\frac{D_i^{7/3}}{{\mathrm{\Σ}}_{j=1}^N{D}_j^{7/3}}\×CFR---\[A4\]$

${\left(R_{outlet}\right)}_i=\frac{P_{aorta}-63mmHg}{{\left(Q_{outlet}\right)}_i}---\[A5\]$

其中，Q_total为心肌体积乘以心肌血流密度，(Q_outlet)_i为出口i的血流，D_i为出口i的直径，(R_outlet)_i为出口i的阻力，冠状动脉血流储备CFR设为2.7；计算出三维网格中各点的压力(P_outlet）_i，通过公式A6得到FFR：

$FFR=\frac{P_{aorta}-{\left(P_{outlet}\right)}_i}{P_{aorta}}---\[A6\]$

其中，P_aorta是主动脉平均压，(P_outlet)_i是三维网格图像中各点的压力值。

2.根据权利要求1所述的基于心脏CT图像的冠状动脉血流储备分数计算方法，其特征在于，所述步骤S02中确定冠状动脉口通过如下步骤得到：对含有冠状动脉口的图像进行连通域分析，用不同的灰度标签标识各个连通域，确定冠状动脉口。

3.根据权利要求1所述的基于心脏CT图像的冠状动脉血流储备分数计算方法，其特征在于，所述步骤S02中的主动脉图像为在心脏图像上，利用主动脉截面成圆形的特征，提取升主动脉及中心线，得到主动脉图像。

4.根据权利要求1所述的基于心脏CT图像的冠状动脉血流储备分数计算方法，其特征在于，所述步骤S04中冠状动脉图像二值化包括：遍历冠状动脉图像V1中体素，如果像素A1等于0，则该像素值不变；如果A1不等于0，则将A1的像素值设为1，得到一个新的数据V2。

5.根据权利要求1所述的基于心脏CT图像的冠状动脉血流储备分数计算方法，其特征在于，所述步骤S05还包括：

基于CT重构的几何结构，将有狭窄的血管拉直，即构建二维轴对称模型，划分二维网格，用数值法求解连续性和Navier-Stokes方程：

$\frac{1}{r}\frac{\∂}{\∂r}\left({\mathrm{ru}}_r\right)+\frac{\∂u_z}{\∂z}=0---\[A7\]$

$\mathrm{\ρ}(\frac{\∂u_r}{\∂t}+u_r\frac{\∂u_r}{\∂r}+u_z\frac{\∂u_r}{\∂z})=-\frac{\∂p}{\∂r}+\mathrm{\μ}\[\frac{1}{r}\frac{\∂}{\∂r}\left(r\frac{\∂u_r}{\∂r}\right)+\frac{{\∂}^2{u}_r}{\∂z^2}-\frac{u_r}{r^2}\]---\[A8\]$

$\mathrm{\ρ}(\frac{\∂u_z}{\∂t}+u_r\frac{\∂u_z}{\∂r}+u_z\frac{\∂u_z}{\∂z})=-\frac{\∂p}{\∂r}+\mathrm{\μ}\[\frac{1}{r}\frac{\∂}{\∂r}\left(r\frac{\∂u_z}{\∂r}\right)+\frac{{\∂}^2{u}_z}{\∂z^2}\]---\[A9\]$

其中，u_z、u_r分别表示z向、r方向的流速；

针对不同类型的弯曲，用三维模型计算从入口到出口的压力差，对照二维轴对称模型计算，建立用于存储各种类型的弯曲对二维轴对称结果的修正系数的数据库；

算出压力后对照数据库中的修正系数，精确得到从入口到出口的压力差，然后计算FFR。