本发明涉及冠状动脉影像学评价领域，具体地涉及一种基于造影图像和流体力学模型的微循环阻力指数计算方法。
背景技术：
冠心病是严重的心脏疾病，经皮冠脉动脉介入(PCI)已成为冠心病治疗的一个重要手段，多年的研究表明，通过PCI治疗冠心病尤其是急性冠脉综合征(acute coronary syndrome，ACS)的病死率得以显著降低，但已有的研究表明，即使通过冠脉介入或冠脉搭桥术治疗，仍有25％患者在心外膜血管成功重建后，并未能实现组织水平心肌再灌注，其主要原因在于冠脉微循环功能障碍。冠脉微循环功能受损是决定急性心肌梗塞预后的独立危险因素，因此对患者冠脉微循环功能状态的评估越来越受到重视。
目前评估冠脉微循环可以通过负荷心电图、放射性核素显像技术、核磁共振成像技术等方法。但以上各种评估技术受检查技术手段本身及检查人员影响较大，不能精确评估冠脉微循环功能，同时这些指标反映的是心外膜血管及微循环共同作用的结果。我们需要一个更精确，更稳定的指标，其只反映冠脉微循环功能状态，不是心外膜血管的影响。目前IMR可能是能满足上述条件的一个理想指标。IMR是反映冠状动脉微循环阻力的指标，被定义为冠脉远端动脉压力(Pd)除以最大充血状态下的平均传导时间(Tmn)的倒数，换言之，即Pd与Tmn的乘积，单位为mmHg·s。IMR＜25为正常，＞30为异常，提示微循环阻力升高，而25-30为灰区。
现有的IMR测量方法是通过0.014英寸软压力导丝同步记录冠脉压力和温度，导丝杆上的两个温度感受器探测到温度变化的时间差就可知道盐水从指引导管到达导丝头端温度感受器运行的平均传导时间(transit mean time，Tmn)，根据定义Pd与Tmn的乘积就可得出IMR值。
但上述压力导丝测量IMR需要介入血管末端，增加手术难度和风险，同时压力导丝昂贵的价格也限制其大规模应用。
技术实现要素：
为了解决上述的技术问题，本发明目的是：提供一种于造影图像和流体力学模型的微循环阻力指数计算方法，基于线冠脉造影图像和流体力学模型计算微循环阻力指数(IMR)，无需通过压力导丝测量，操作简便大大降低手术难度和风险，可在临床上大规模推广应用。
本发明的技术方案是：
一种基于造影图像和流体力学模型的微循环阻力指数计算方法，包括以下步骤：
S01：选取两个体位的冠状动脉造影图像，通过分析定义时间序列的X射线造影不同帧，去除动静态噪声，对冠状动脉进行分割；并对不同角度拍摄的冠脉造影图像进行分割，分别得到冠脉中心线及直径；
S02：将S01步骤生成的多方向的冠脉中心线在三维空间进行投影计算，生成冠脉三维模型；
S03：通过TIMI记帧法获取血流在最大充血状态下的传导时间Tmn及血流速度；
S04：基于X射线重构的冠脉三维模型，创建血管三维网格；并基于X射线重构的冠脉中心线及直径，拉直血管，按对称轴切开，生成轴对称二维平面模型，然后创建血管二维轴对称网格；
S05：测量主动脉平均压Pa；
S06：基于S03步骤得到的血液流速与S04步骤生成的血管三维网格，求解不可压缩流的基本公式，计算沿着血管中心线从入口到下游各点的压力降ΔPi，则冠脉远端动脉压力Pd＝Pa-ΔPi；
S07：计算微循环阻力指数IMR＝Pd*Tmn。
优选的技术方案中，所述步骤S01包括以下步骤：
S11：定义时间序列X射线冠脉造影图像中有导管出现的第一帧为参考图像，有完整冠脉出现的需要分割的图像为目标图像；
S12：提取参考图像中导管的特征点，对目标图像中的冠状动脉进行图像增强，提取冠状动脉的区域图像；
S13：以参考图像中导管的特征点为种子点，将冠状动脉的区域图像进一步处理后进行动态区域生长，获得冠脉造影图像。
优选的技术方案中，所述步骤S12包括：
将参考图像减去目标图像，去除静态噪声，通过均值滤波去除动态噪声；
通过灰度直方图分析，利用阈值进一步去除噪声；
通过多尺度海森矩阵对将去噪后的图像的导管进行图像增强；
二值化图像：将导管对应的像素标记为1，其余像素标记0；对二值化后的图像进行细化，得到一组导管的特征点。
优选的技术方案中，所述步骤S12包括：
将目标图像减去参考图像，去除静态噪声，通过均值滤波去除动态噪声；
通过灰度直方图分析，利用阈值进一步去除噪声；
通过多尺度海森矩阵对将去噪后的图像的冠状动脉进行图像增强，得到增强后的目标图像；
根据增强后的目标图像中各区域与特征点的位置关系，确定冠状动脉的区域；提取该连通区域，获得冠状动脉的区域图像。
优选的技术方案中，所述步骤S13包括：
对冠状动脉的区域图像进行二值化；
进行形态学运算，按照对应位置还原到目标图像中，进行动态区域生长，获得冠脉造影图像。
优选的技术方案中，所述步骤S02包括：
心跳周期内相同时间点、不同角度拍摄的冠脉造影图像，用步骤S01处理后，在中心线上提取特征点；
基于已知的X射线拍照角度，三维还原这些特征点，进而三维还原中心线。
优选的技术方案中，所述步骤S03中计算血液流速的方法使用下述方法中的一种：造影剂运输时间算法、造影剂遍历距离算法、Stewart—Hamilton算法、First—pass分布分析法、光流法或者流体连续法。
优选的技术方案中，所述步骤S04包括：
基于步骤S02中重构的冠脉三维模型，用标准的扫掠法、切分法或混合法进行网格划分，生成结构性三维六面体网格；
基于步骤S02中重构的冠脉三维模型，在血管壁上生成三角形面网格，然后长成四面体网格，在血管中生成非结构性三维四面体网格；
基于步骤S02中重构的冠脉三维模型，拉直血管，按对称轴切开，生成轴对称二维平面模型，用切分法或者扫掠法进行网格划分，生成结构性二维四边形网格。
优选的技术方案中，所述步骤S05包括：
对血管三维网格进行求解，用数值法求解连续性和Navier-Stokes方程：
其中，P，ρ，μ分别为流速、压力、血流密度、血流粘性；
入口边界条件为血流速度，出口边界条件为out-flow边界条件。
优选的技术方案中，所述步骤S05还包括：
对血管二维轴对称网格，数值法求解连续性和Navier-Stokes方程：
其中，ρ表示血液的密度，uz、ur分别表示z向、r方向的流速，μ表示血液的动力粘度，p表示血液的压强；
入口边界条件为血流速度，出口边界条件为out-flow边界条件；
针对各种类型的冠脉弯曲，用三维模型计算从入口到出口的压力差，对照二维轴对称模型计算的结果，建立用于存储各种类型的弯曲对二维轴对称结果的修正系数的数据库；利用三维模型计算结果修正二维模型计算结果得到ΔPi值。
与现有技术相比，本发明的优点是：
1、基于线冠脉造影图像和流体力学模型计算微循环阻力指数(IMR)，速度快，准确性高，不需要额外插入压力导丝进行测量，操作简便，大大降低手术难度和风险，可在临床上大规模推广应用。
2、该方法还有助于临床医师详细了解患者的冠脉功能状态，明晰患者血管及患者病变并优化介入策略，评估微循环功能状态并指导药物治疗，并且对PCI患者的预后提供更为准确的预测。
附图说明
下面结合附图及实施例对本发明作进一步描述：
图1为本发明的方法流程图；
图2为参考图像；
图3为待分割的一目标图像；
图4为待分割的另一目标图像；
图5为增强后的导管图像；
图6为导管的特征点；
图7为增强后的目标图像；
图8为目标图像的冠状动脉图像；
图9为分割的结果图像；
图10为横切面截图；
图11为纵切面截图；
具体实施方式
为使本发明的目的、技术方案和优点更加清楚明了，下面结合具体实施方式并参照附图，对本发明进一步详细说明。应该理解，这些描述只是示例性的，而并非要限制本发明的范围。此外，在以下说明中，省略了对公知结构和技术的描述，以避免不必要地混淆本发明的概念。
如图1所示，本发明的基于造影图像和流体力学模型的微循环阻力指数计算方法，包括以下步骤。
定义时间序列X射线冠脉造影图像中的有导管出现的第一帧为参考图像A，如图2所示，有完整冠脉出现的需要分割的为目标图像B，本发明选取两个体位的冠状动脉造影图像进行分割，如图3、4所示，最后得到的分割结果为结果图像P4。
1：提取冠脉造影图像：
从时间序列图像中可以明显看出，参考图像A(图2)和目标图像B(图3)中的图像像素可以分成三部分：需要分割的结果，称之为导管或者冠状动脉；在时间中静止不变的噪声，如胸腔中的肋骨，称之为静态噪声；在时间中变化的噪声，如部分肺部组织、部分心脏组织，称之为动态噪声。
首先，因为参考图像A(图2)中的导管结构比较简单，容易定位，所以先提取参考图像A中的导管中的特征点。
先用参考图像A-目标图像B，除去部分静态噪声；
做均值滤波，除去部分动态噪声；
通过灰度直方图分析，利用阈值进一步去噪；
将去噪后的图像进行多尺度海森矩阵，对导管进行图像增强，得到增强后的导管图像，如图5所示；
二值化图像，将导管对应的像素标记为1，其余像素标记0；
对二值化后的图像进行细化，得到一组导管的特征点P1，如图6所示。
其次，对目标图像B(图3)中的冠状动脉进行图像增强。
用目标图像B-参考图像A，除去部分静态噪声；
均值滤波，除去部分动态噪声；
通过灰度直方图分析，利用阈值进一步去噪；
通过多尺度海森矩阵对去噪后的图像的冠状动脉进行图像增强，得到图像P2，如图7所示。
通过对目标图像B进行以上操作得到的图像中，除被增强后的冠状动脉之外，还可能有被海森矩阵误认为管状物而进行了增强的噪声。于是需要将图像P2中冠状动脉和噪声进行区分。
特征点P1是由导管细化后得到的，所以P1的位置一定在导管中；目标图像和参考图像中，导管的位置虽然会有一定变化，但是该变化必定只是在一个小范围内的。所以，可以根据图像P2中各个区域与特征点P1的位置关系，确定冠状动脉的区域；
提取该连通区域，获得图像P3，如图8所示。
最后，以特征点P1为种子点，将图像P3进一步处理后进行动态区域生长，获得结果图像，如图9所示。
P3图像二值化；
形态学运算；
按照对应位置还原到目标图像中；
进行动态区域生长，得到结果图像P4，如图9所示。
2：生成冠脉三维模型：
对两张不同角度拍摄的冠脉造影图像按照上述步骤进行分割，分别得到冠脉中心线及直径，将上述步骤生成的两个方向的冠脉中心线在三维空间进行投影计算，生成冠脉三维模型即可获取准确的血管长度和直径；
3：计算冠脉血流速度：
通过TIMI记帧法获取血流在最大充血状态下的传导时间Tmn及血流速度，对冠脉造影图像按时序帧进行扫描计算出造影剂从血管入口流到血管末端需要的时间，再通过血管长度除以时间得到血液流速；
4：冠脉网格划分：
基于重构的冠脉三维模型，用标准的扫掠法进行网格划分，生成结构性三维六面体网格；
基于重构的冠脉三维模型，用其它方法(例如：切分法、混合法)进行网格划分，生成结构性三维六面体网格；
基于重构的冠脉三维模型，在血管壁上生成三角形面网格，然后长成四面体网格，在血管中生成非结构性三维四面体网格；
基于重构的冠脉三维模型，拉直血管，按对称轴切开，生成轴对称二维平面模型，用切分法(或者扫掠法)进行网格划分，生成结构性二维四边形网格，图10为网格横切面截图，图11为网格纵切面截图。
5：测量主动脉平均压Pa，例如通过指引导管测量。
6：冠脉远端动脉压力Pd计算：
1)三维模型
基于X射线造影图像重构的几何结构，划分三维网格，用数值法(如：有限差分、有限元、有限体积法等)求解连续性和Navier-Stokes方程：
其中，P，ρ，μ分别为流速、压力、血流密度、血流粘性。
入口边界条件为血流速度，出口边界条件为out-flow边界条件；
计算沿着血管中心线从入口到下游各点的压力降ΔPi，则冠脉远端动脉压力Pd＝Pa-ΔPi。
2)二维模型
基于X射线造影重构的几何结构，把有狭窄的血管拉直(二维轴对称模型)，划分二维网格，用数值法(如：有限差分、有限元、有限体积法等)求解连续性和Navier-Stokes方程：
其中，ρ表示血液的密度，uz、ur分别表示z向、r方向的流速，μ表示血液的动力粘度，p表示血液的压强。
入口边界条件为血流速度，出口边界条件为out-flow边界条件；
针对各种类型的冠脉弯曲，用三维模型计算从入口到出口的压力差，对照二维轴对称模型计算的结果，建立用于存储各种类型的弯曲对二维轴对称结果的修正系数的数据库；利用三维模型计算结果修正二维模型计算结果得到ΔPi；
则冠脉远端动脉压力Pd＝Pa-ΔPi。7：计算微循环阻力指数IMR＝Pd*Tmn(单位mmHg·s)。
应当理解的是，本发明的上述具体实施方式仅仅用于示例性说明或解释本发明的原理，而不构成对本发明的限制。因此，在不偏离本发明的精神和范围的情况下所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。此外，本发明所附权利要求旨在涵盖落入所附权利要求范围和边界、或者这种范围和边界的等同形式内的全部变化和修改例。