阅读说明：本技术 一种利用cta和dsa测量冠状动脉血流储备分数的方法 (Method for measuring coronary artery flow reserve fraction by using CTA (computed tomography angiography) and DSA (digital signal amplification) ) 是由 陈荣民 于 2020-11-20 设计创作，主要内容包括：本发明公开了一种利用CTA和DSA测量冠状动脉血流储备分数的方法,包括以下步骤：利用CTA冠脉造影进行三维重建,得到获得冠脉几何模型G1、心肌体积V以及目标血管的长度L；根据几何模型G1和目标血管的近端入口位置,截取用于仿真模拟的分支几何模型G2；利用连续时间内的DSA造影图得到血流流经目标血管的时间t,根据速度公式结合L和t计算出目标血管近端的血流速度Vin；结合心肌体积V以及血压,计算出几何模型G1的远端阻力Rtotal；结合Rtotal、G1和G2的分支个数,计算G2的远端阻力R；以Vin作为G2的入口边界条件,R作为G2的远端微循环阻力,通过求解流体力学方程计算出FFR值。本发明结合了两种技术的优势,并依靠流体力学公式计算出所需的FFR值。(The invention discloses a method for measuring coronary artery blood flow reserve fraction by using CTA and DSA, which comprises the following steps: performing three-dimensional reconstruction by using CTA coronary angiography to obtain a coronary geometric model G1, a myocardial volume V and a length L of a target blood vessel; intercepting a branch geometric model G2 for simulation according to the geometric model G1 and the proximal entrance position of the target blood vessel; obtaining the time t of the blood flow flowing through the target blood vessel by using DSA (digital angiography) contrast maps in continuous time, and calculating the blood flow velocity Vin at the near end of the target blood vessel by combining L and t according to a velocity formula; calculating the far-end resistance Rtotal of the geometric model G1 by combining the myocardial volume V and the blood pressure; calculating the far-end resistance R of G2 by combining the branch numbers of Rtotal, G1 and G2; the FFR value is calculated by solving the fluid mechanics equation with Vin as the inlet boundary condition for G2 and R as the distal microcirculation resistance for G2. The invention combines the advantages of the two technologies and calculates the required FFR value by relying on a fluid mechanics formula.)

1. A method for measuring coronary fractional flow reserve using CTA and DSA, comprising the steps of:

performing three-dimensional reconstruction by using CTA coronary angiography to obtain a coronary geometric model G1, a myocardial volume V and a length L of a target blood vessel;

intercepting a branch geometric model G2 for simulation according to the geometric model G1 and the proximal entrance position of the target blood vessel;

obtaining the time t of the blood flow flowing through the target blood vessel by using DSA (digital angiography) contrast maps in continuous time, and calculating the blood flow velocity Vin at the near end of the target blood vessel by combining L and t according to a velocity formula;

calculating the far-end resistance Rtotal of the geometric model G1 by combining the myocardial volume V and the blood pressure;

calculating the far-end resistance R of G2 by combining the branch numbers of Rtotal, G1 and G2;

the FFR value is calculated by solving the fluid mechanics equation with Vin as the inlet boundary condition for G2 and R as the distal microcirculation resistance for G2.

2. The method of claim 1, wherein the acquisition of time t for blood flow through a target vessel comprises: DSA contrast obtains image data, each frame of image is arranged in sequence, in the image sequence of the sequence, the image of the contrast agent flowing to the appointed initial position is used as the first frame of image, the image of the contrast agent flowing to the appointed end position is used as the last frame of image, and then according to the frame rate of the DSA image, the time of the contrast agent from the initial position to the end point, namely the blood flow time t, is calculated.

3. The method of claim 1A method for measuring coronary flow reserve fraction by CTA and DSA, wherein the obtaining of the distal resistance Rtotal of G1 comprises: according to the formula Q ═ V^βObtaining flow Q, wherein V is the myocardial volume, beta is a coefficient, and beta is 1.5; calculating an average pressure P (1/3 Pa +2/3 Pb) according to the systolic pressure Pa and the diastolic pressure Pb of the case; rtotal is calculated according to ohm's law Rtotal ═ P/Q.

4. The method of claim 1, wherein the distal resistance R of G2 is calculated as: r is Rtotal m/n, wherein m is the number of coronary outlets of G1, and n is the number of coronary outlets of G2.

5. A method of measuring coronary flow reserve fraction using CTA and DSA as claimed in claim 3, wherein the calculation of the FFR value comprises: obtaining the blood flow velocity of the G2 inlet according to the blood flow velocity Vin as an inlet boundary condition, obtaining the microcirculation resistance of the far-end blood vessel of the G2 outlet according to R by adopting a three-unit impedance model at the outlet boundary, calculating the far-end blood pressure Pd by solving a fluid mechanics equation, and obtaining the target coronary artery FFR value according to the FFR (fractional flow rate) Pd/Pa.

6. The method of claim 5, wherein the acquisition of the blood flow velocity at the G2 inlet comprises: with Vin as the maximum value, a piecewise function v (t) is constructed, where v (t) Vin x sin (pi t) when 0< t <0.25s, and v (t) Vin when 0.25s ≦ t ≦ 0.5 s.

7. The method of claim 6, wherein the obtaining of the distal vascular microcirculation resistance comprises: let G2 have n outlets, the distal microcirculation resistance of the ith outlet is R_iThe resistance of the microcirculation of the blood vessels at the far ends of each branch is equal according to ohm's law R_i＝n*R(i＝1,2...n)。

8. The method of claim 7, wherein the obtaining of the distal blood pressure Pd comprises:

u(t)＝Vin(t) on Γ^inlet

wherein, the inlet adopts a speed boundary condition, and the outlet adopts a three-unit impedance model; respectively using resistorsAndrepresenting the resistance of the proximal and downstream vessels, wherein the resistance of the microcirculation distal to the ith outletBy means of a capacitor C_iIndicates the compliance of the downstream blood vessel, C_iObtained by flow calculation; gamma-shaped^inletAnd Γ^outletCalculating an entry boundary and an exit boundary, p, of a domain for a fluid, respectively_iIs the ith outlet pressure, Q_iIs the flow rate of the ith outlet and,n_iis the normal vector of the ith outlet; p is a radical of_i(0) And Q_i(0) Initial pressure and flow rate of the ith outlet respectively; since p is defined throughout the fluid calculation region, solving the above equation results in p including the blood pressure value Pd at the distal end of the blood vessel.

Technical Field

The invention relates to the field of biomedical engineering, in particular to a method for measuring coronary artery blood flow reserve fraction by using CTA and DSA.

Background

The fractional flow reserve, FFR for short, is an important index for judging the smoothness of blood flow in coronary artery, and refers to the ratio of the maximum blood flow obtained from the myocardial region of the blood vessel to the maximum blood flow obtained from the same region under the theoretical normal condition, i.e. the ratio of the mean pressure (Pd) in the narrow distal coronary artery to the mean pressure (Pa) at the coronary artery entrance under the maximum hyperemia state of the myocardium. At present, the FFR value is mainly obtained by measuring the coronary artery of a human body by a clinician through intervention of a pressure guide wire, and the measuring method has higher cost and certain risk.

Several non-invasive FFR measurement methods have been developed in recent years, such as nuclide imaging, magnetic resonance perfusion, CT non-invasive fractional flow reserve measurement (i.e., FFR-CT), DSA-based FFR measurement (i.e., QFR). The former two methods are also traditional methods, and have strong dependence on equipment, for example, a nuclide imaging method requires a SPECT (single photon emission computed tomography) device, a magnetic resonance perfusion method requires a magnetic resonance device, and the detection cost is high. The latter two methods have been recently developed, and QFR was first approved for use in china in the last year. Other techniques are also included, such as a method of determining the FFR of the coronary arteries as disclosed in the invention of publication No. CN 108992057A.

FFR-CT was the first FDA certificate obtained in 2014 by the united states, which is based on CT coronary angiography (CTA) images, simulates coronary blood flow, simulates major cardiovascular indexes of human body, such as blood pressure, blood flow velocity, coronary blood flow reserve fraction, and the like, and is used for post-processing of clinical quantitative analysis. The flow method comprises the steps of firstly obtaining a CTA image in a CT room, then processing the CTA image, reconstructing a three-dimensional structure of coronary artery, then carrying out grid and mathematical modeling on the three-dimensional structure, and finally submitting some boundary conditions to a super computer for calculation to obtain the FFR value of each coronary artery. The method obtains the three-dimensional structure of the coronary artery based on the image of CTA, but cannot accurately obtain boundary conditions, such as blood pressure or blood flow velocity at the outlet. In addition, the method needs to carry out hemodynamic numerical simulation on the global coronary blood flow region and a part of ascending aorta region, and the calculation time is long.

The FFR measurement method based on DSA (QFR) is that DSA plane images of two angles of coronary artery are firstly obtained under a DSA device through DSA images, then the diameter of the blood vessel is measured according to the images, then a three-dimensional blood vessel is reconstructed according to the diameter, then the blood flow speed is determined according to the flowing condition of contrast agent in the blood vessel, and a blood pressure sensor designed at a contrast agent catheter is prepared in advance to obtain the blood pressure at the entrance of the coronary artery, so that the FFR value of the coronary artery can be calculated, the calculation amount is small, and a single computer can calculate a case in a few minutes. QFR, the calculation speed is high, but there is a big disadvantage that the DSA image is a plane image, the three-dimensional structure of the blood vessel can not be obtained accurately, especially the three-dimensional structure at the narrow part, which is decisive for the calculation result.

The above prior arts all have respective defects, which results in that the parameters of blood flow and the blood vessel model cannot be well combined, so that the calculation result is distorted, and the accuracy and efficiency of calculation cannot be obtained at the same time.

Disclosure of Invention

Aiming at the problem that the blood flow parameters and the model cannot be well combined in the prior art, the invention provides a method for measuring the fractional flow reserve of the coronary artery by using CTA (computed tomography angiography) and DSA (computed angiography-specific angiography), which combines the advantages of the two technologies and calculates the required FFR (fractional flow rate) value by depending on a fluid mechanics formula.

The technical scheme of the invention is as follows.

A method for measuring coronary flow reserve fraction using CTA and DSA, comprising the steps of: performing three-dimensional reconstruction by using CTA coronary angiography to obtain a coronary geometric model G1, a myocardial volume V and a length L of a target blood vessel; intercepting a branch geometric model G2 for simulation according to the geometric model G1 and the proximal entrance position of the target blood vessel; obtaining the time t of the blood flow flowing through the target blood vessel by using DSA (digital angiography) contrast maps in continuous time, and calculating the blood flow velocity Vin at the near end of the target blood vessel by combining L and t according to a velocity formula; calculating the far-end resistance Rtotal of the geometric model G1 by combining the myocardial volume V and the blood pressure; calculating the far-end resistance R of G2 by combining the branch numbers of Rtotal, G1 and G2; the FFR value is calculated by solving the fluid mechanics equation with Vin as the inlet boundary condition for G2 and R as the distal microcirculation resistance for G2.

According to the scheme, CTA can be used for accurately modeling and obtaining the length of the blood vessel, but the blood volume velocity cannot be accurately measured, DSA can be used for obtaining the blood flow time but cannot be used for accurately measuring the length of the blood vessel, and after the two methods are combined, the blood flow velocity can be calculated by utilizing respective advantage parts and respective accurate parameters, so that the FFR value can be calculated. Meanwhile, by adopting the CTA image and the DSA image, the three-dimensional reconstruction can be accurately carried out, and the boundary condition can be accurately obtained, so that the FFR value can be more accurately calculated.

The acquisition of the myocardial volume V is based on CTA images, the segmentation processing is carried out on the myocardium through medical image software, and then the volume V is obtained by measuring. The acquisition process of the target blood vessel length L comprises the following steps: acquiring a central line of G2; determining the main branch vessel of G2; the length of the main branch vessel centerline is measured as L.

Preferably, the obtaining of the time t when the blood flow flows through the target blood vessel includes: DSA contrast obtains image data, each frame of image is arranged in sequence, in the image sequence of the sequence, the image of the contrast agent flowing to the appointed initial position is used as an initial frame image, the image of the contrast agent flowing to the appointed end position is used as an end frame image, and then the time of the contrast agent from the initial position to the end point, namely the blood flow time t, is calculated according to the frame rate of the DSA images. The images are arranged according to the time sequence of the frames, so that the time can be segmented, the flowing condition of the contrast agent can be observed conveniently, and the time calculation is facilitated.

Preferably, the obtaining process of the distal resistance Rtotal of G1 includes: according to the formula Q ═ V^βObtaining flow Q, wherein V is the myocardial volume, beta is a coefficient, and beta is 1.5; calculating an average pressure P (1/3 Pa +2/3 Pb) according to the systolic pressure Pa and the diastolic pressure Pb of the case; rtotal is calculated according to ohm's law Rtotal ═ P/Q.

Preferably, the calculation formula of the distal resistance R of G2 is: r is Rtotal m/n, wherein m is the number of coronary outlets of G1, and n is the number of coronary outlets of G2.

Preferably, the calculation process of the FFR value includes: obtaining the blood flow velocity of the G2 inlet according to the blood flow velocity Vin as an inlet boundary condition, obtaining the microcirculation resistance of the far-end blood vessel of the G2 outlet according to R by adopting a three-unit impedance model at the outlet boundary, calculating the far-end blood pressure Pd by solving a fluid mechanics equation, and obtaining the target coronary artery FFR value according to the FFR (fractional flow rate) Pd/Pa.

Preferably, the process of acquiring the blood flow velocity of the G2 inlet comprises: constructing a piecewise function V (t) with Vin as the maximum value, wherein V (t) is Vin sin (pi t) when t is more than 0 and less than 0.25s, and V (t) is Vin when t is more than 0.25s and less than 0.5 s. In the numerical simulation stage, a quasi-steady state model is adopted, the calculation time is set to be 0.5s, the unit time step length is set to be 0.01s, and a piecewise function is constructed.

Preferably, the process of obtaining the microcirculation resistance of the distal blood vessel comprises the following steps: let G2 have n outlets, the distal microcirculation resistance of the ith outlet is R_iThe resistance of the microcirculation of the blood vessels at the far ends of each branch is equal according to ohm's law R_i＝n*R(i＝1，2…n)。

The fluid mechanics governing equation for this example is as follows:

wherein the content of the first and second substances,which is the stress tensor of the flow field, u denotes the blood flow velocity,is a gradient operator, t is a time variable, I is an identity matrix,to representF is the source term, here gravity, p is blood pressure, ρ is blood density, μ is viscosity coefficient of blood, and Ω is the fluid calculation region.

The boundary conditions are set as follows:

u(t)＝Vin(t)onΓ^inlet

wherein, the inlet adopts a speed boundary condition, and the outlet adopts a three-unit impedance model; respectively using resistorsAndrepresenting the resistance of the proximal and downstream vessels, wherein the resistance of the microcirculation distal to the ith outletBy means of a capacitor C_iIndicates the compliance of the downstream blood vessel, C_iObtained by flow calculation; gamma-shaped^inletAnd Γ^outletCalculating an entry boundary and an exit boundary, p, of a domain for a fluid, respectively_iIs the ith outlet pressure, Q_iIs the flow rate of the ith outlet and,n_iis the normal vector of the ith outlet; p is a radical of_i(0) And Q_i(0) Initial pressure and flow rate of the ith outlet respectively; since p is defined throughout the fluid calculation region, solving the above equation results in p including the blood pressure value Pd at the distal end of the blood vessel.

The substantial effects of the invention include: meanwhile, accurate data in CTA and DSA are adopted for calculation; in the CTA reconstructed image, the length required by calculating the blood flow velocity by adopting a midline method is higher than that directly calculated from DSA, and the far-end blood flow velocity is accurate, so that the calculated far-end pressure Pd is accurate, and the final FFR is accurate naturally; in DSA image data, accurate flow time information is obtained according to a frame image corresponding to the position of a contrast agent, and the far-end blood flow velocity is more accurate; the three-dimensional structure of the blood vessel is extracted by CTA to obtain the three-dimensional structure of the branch blood vessel, and the component of the blood flow of the branch blood vessel can be calculated, so that the final calculation result can be calculated more accurately; on the whole, the mathematical model, the blood vessel length and the calculation parameter setting required by the invention can be completed before the DSA is carried out after the CTA is finished, and the calculation can be rapidly submitted as long as the blood flow time is obtained after the DSA radiography is carried out. In addition, the area of the simulation solution is limited to the lesion blood vessel and the downstream branch blood vessel, so the whole time can be greatly shortened, the calculation cost is reduced, enough time is provided for quality control, and the calculation result can be more reliably grasped.

Drawings

FIG. 1 is a schematic diagram of the acquisition process of G1 and G2 according to the embodiment of the present invention;

fig. 2 is an image of a DSA start frame and an end frame according to an embodiment of the present invention.

Detailed Description

The technical solution of the present application will be described with reference to the following examples. In addition, numerous specific details are set forth below in order to provide a better understanding of the present invention. It will be understood by those skilled in the art that the present invention may be practiced without some of these specific details. In some instances, methods, means, elements and circuits that are well known to those skilled in the art have not been described in detail so as not to obscure the present invention.

Example (b):

a method for measuring coronary flow reserve fraction using CTA and DSA, comprising the steps of: as shown in fig. 1, a three-dimensional reconstruction is performed by using CTA coronary angiography to obtain a coronary geometric model G1, a myocardial volume V, and a length L of a target blood vessel; intercepting a branch geometric model G2 for simulation according to the geometric model G1 and the proximal entrance position of the target blood vessel; obtaining the time t of the blood flow flowing through the target blood vessel by using DSA (digital angiography) contrast maps in continuous time, and calculating the blood flow velocity Vin at the near end of the target blood vessel by combining L and t according to a velocity formula; calculating the far-end resistance Rtotal of the geometric model G1 by combining the myocardial volume V and the blood pressure; calculating the far-end resistance R of G2 by combining the branch numbers of Rtotal, G1 and G2; the FFR value is calculated by solving the fluid mechanics equation with Vin as the inlet boundary condition for G2 and R as the distal microcirculation resistance for G2.

According to the scheme, CTA can be used for accurately modeling and obtaining the length of the blood vessel, but the blood volume velocity cannot be accurately measured, DSA can be used for obtaining the blood flow time but cannot be used for accurately measuring the length of the blood vessel, and after the two methods are combined, the blood flow velocity can be calculated by utilizing respective advantage parts and respective accurate parameters, so that the FFR value can be calculated. Meanwhile, by adopting the CTA image and the DSA image, the three-dimensional reconstruction can be accurately carried out, and the boundary condition can be accurately obtained, so that the FFR value can be more accurately calculated.

The acquisition of the myocardial volume V is based on CTA images, the segmentation processing is carried out on the myocardium through medical image software, and then the volume V is obtained by measuring. The acquisition process of the target blood vessel length L comprises the following steps: acquiring a central line of G2; determining the main branch vessel of G2; the length of the main branch vessel centerline is measured as L.

As shown in fig. 2, the acquisition process of the time t when the blood flow flows through the target blood vessel includes: DSA contrast obtains image data, each frame of image is arranged in sequence, in the image sequence of the sequence, the image of the contrast agent flowing to the appointed initial position is used as an initial frame image, the image of the contrast agent flowing to the appointed end position is used as an end frame image, and then the time of the contrast agent from the initial position to the end point, namely the blood flow time t, is calculated according to the frame rate of the DSA images. The images are arranged according to the time sequence of the frames, so that the time can be segmented, the flowing condition of the contrast agent can be observed conveniently, and the time calculation is facilitated.

Wherein the obtaining process of the far-end resistance Rtotal of G1 comprises the following steps: according to the formula Q ═ V^βObtaining flow Q, wherein V is the myocardial volume, beta is a coefficient, and beta is 1.5; systolic blood pressure Pa andcalculating the mean pressure P (1/3 Pa +2/3 Pb) of the diastolic pressure Pb; rtotal is calculated according to ohm's law Rtotal ═ P/Q.

Wherein the distal resistance R of G2 is calculated as: r is Rtotal m/n, wherein m is the number of coronary outlets of G1, and n is the number of coronary outlets of G2.

The calculation process of the FFR value comprises the following steps: obtaining the blood flow velocity of the G2 inlet according to the blood flow velocity Vin as an inlet boundary condition, obtaining the microcirculation resistance of the far-end blood vessel of the G2 outlet according to R by adopting a three-unit impedance model at the outlet boundary, calculating the far-end blood pressure Pd by solving a fluid mechanics equation, and obtaining the target coronary artery FFR value according to the FFR (fractional flow rate) Pd/Pa.

Wherein the process of acquiring the blood flow velocity of the G2 inlet comprises the following steps: constructing a piecewise function V (t) with Vin as the maximum value, wherein V (t) is Vin sin (pi t) when t is more than 0 and less than 0.25s, and V (t) is Vin when t is more than 0.25s and less than 0.5 s. In the numerical simulation stage, a quasi-steady state model is adopted, the calculation time is set to be 0.5s, the unit time step length is set to be 0.01s, and a piecewise function is constructed.

The acquisition process of the microcirculation resistance of the far-end blood vessel comprises the following steps: let G2 have n outlets, the distal microcirculation resistance of the ith outlet is R_iThe resistance of the microcirculation of the blood vessels at the far ends of each branch is equal according to ohm's law R_i＝n*R(i＝1，2…n)。

The fluid mechanics governing equation for this example is as follows:

wherein the content of the first and second substances,which is the stress tensor of the flow field, u denotes the blood flow velocity,is a gradient operator, t is a time variable, I is an identity matrix,to representF is the source term, here gravity, p is blood pressure, ρ is blood density, μ is viscosity coefficient of blood, and Ω is the fluid calculation region.

The boundary conditions are set as follows:

u(t)＝Vin(t)onΓ^inlet

wherein, the inlet adopts a speed boundary condition, and the outlet adopts a three-unit impedance model; respectively using resistorsAndrepresenting the resistance of the proximal and downstream vessels, wherein the resistance of the microcirculation distal to the ith outletBy means of a capacitor C_iIndicates the compliance of the downstream blood vessel, C_iObtained by flow calculation; gamma-shaped^inletAnd Γ^outletCalculating an entry boundary and an exit boundary, p, of a domain for a fluid, respectively_iIs the ith outlet pressure, Q_iIs the flow rate of the ith outlet and,n_iof the ith outletA normal vector; p is a radical of_i(0) And Q_i(0) Initial pressure and flow rate of the ith outlet respectively; since p is defined throughout the fluid calculation region, solving the above equation results in p including the blood pressure value Pd at the distal end of the blood vessel.

The substantial effects of the present embodiment include: meanwhile, accurate data in CTA and DSA are adopted for calculation; in the CTA reconstructed image, the length required by calculating the blood flow velocity by adopting a midline method is higher than that directly calculated from DSA, and the far-end blood flow velocity is accurate, so that the calculated far-end pressure Pd is accurate, and the final FFR is accurate naturally; in DSA image data, accurate flow time information is obtained according to a frame image corresponding to the position of a contrast agent, and the far-end blood flow velocity is more accurate; the three-dimensional structure of the blood vessel is extracted by CTA to obtain the three-dimensional structure of the branch blood vessel, and the component of the blood flow of the branch blood vessel can be calculated, so that the final calculation result can be calculated more accurately; on the whole, the mathematical model, the blood vessel length and the calculation parameter setting required by the embodiment can be completed before the DSA is performed after the CTA is completed, and calculation can be rapidly submitted as long as the blood flow time is obtained after the DSA radiography is performed. In addition, the area of the simulation solution is limited to the lesion blood vessel and the downstream branch blood vessel, so the whole time can be greatly shortened, the calculation cost is reduced, enough time is provided for quality control, and the calculation result can be more reliably grasped.

The technical solution of the embodiments of the present application may be essentially or partially contributed to the prior art, or all or part of the technical solution may be embodied in the form of a software product, where the software product is stored in a storage medium, and includes several instructions to enable a device (which may be a single chip, a chip, or the like) or a processor (processor) to execute all or part of the steps of the method of the embodiments of the present application. And the aforementioned storage medium includes: various media capable of storing program codes, such as a usb disk, a removable hard disk, a Read Only Memory (ROM), a Random Access Memory (RAM), a magnetic disk, or an optical disk.

The above description is only for the specific embodiments of the present application, but the scope of the present application is not limited thereto, and any person skilled in the art can easily conceive of the changes or substitutions within the technical scope of the present application, and shall be covered by the scope of the present application. Therefore, the protection scope of the present application shall be subject to the protection scope of the claims.