Pressure drop estimation

Abstract

Embodiments and aspects described herein provide methods and systems for determining pressure difference across a tube arising from fluid flow within the tube, comprising: obtaining three-dimensional time dependent fluid velocity data at a plurality of points along the tube; processing the three-dimensional time dependent fluid velocity data to determine: i) a flow rate (Q) of the fluid through the tube; ii) the kinetic energy (K) of the fluid flow through the tube; iii) an advective energy rate (A) of the fluid flow through the tube; and iv) a viscous dissipation rate (V) pertaining to the fluid flow; and calculating the pressure difference in dependence on all of the flow rate (Q), kinetic energy (K), advective energy rate (A), and viscous dissipation rate (V). Further embodiments are also described.

Claims

1. A method of determining pressure difference across a tube arising from fluid flow within the tube, comprising: obtaining fluid velocity data at the inlet and outlet planes of the tube; and processing the fluid velocity data to determine: i) a flow rate (Q) of the fluid through the tube; and ii) an advective energy rate (A) of the fluid flow through the tube; the method further comprising calculating the pressure difference in dependence on both of the flow rate (Q), and the advective energy rate (A), without considering a viscous dissipation rate (V) and the kinetic energy (K) of the fluid flow through the tube, wherein: a) the flow rate Q is dependent on a surface integral of the fluid velocity data across either the inlet or outlet plane of the tube, or across any other plane dividing the tube in two sub-sections; and/or b) the advective energy rate A is dependent on a sum of the surface integrals of the time dependent fluid velocity data, or data derived therefrom, across the inlet and outlet planes and the fluid density.

2. A method according to claim 1, wherein the pressure difference across the tube is given by one of:
Δp.sub.cAW=ρ/2(∫.sub.Γ.sub.OUTLET|v|.sup.2(v.Math.n)dx+∫.sub.Γ.sub.INLET|v|.sup.2(v.Math.n)dx)/∫.sub.Γ.sub.OUTLETv.Math.ndx  A) where v is a velocity field at a generic voxel, n is the normal direction on inlet and outlet planes, rho is the fluid density respectively, and Γ.sub.INLET and Γ.sub.OUTLET are the respective inlet and outlet planes of the tube; or
Δp.sub.sAW=−ρ/2∫.sub.Γ.sub.OUTLET|v|.sup.2(v.Math.n)dx/∫.sub.Γ.sub.OUTLETv.Math.ndx  B) where v is a velocity field at a generic voxel, rho is the fluid density respectively, Γ.sub.OUTLET is a selected plane of the tube, and n is the normal direction on the selected plane, wherein the selected plane is one of the planes at any end of the tube.

3. A method according to claim 1, wherein fluid velocity data is not obtained for anywhere else in the tube other than the inlet and outlet planes.

4. A method according to claim 1, wherein fluid velocity data is not obtained for anywhere else in the tube other than the selected plane.

5. A method according to claim 1, wherein the tube is a blood vessel, and the fluid flow is in vivo blood flow, and wherein the fluid velocity data is obtained from one of: i) 4D phase contrast magnetic resonance imaging data; and/or ii) 2D phase contrast magnetic resonance imaging data; and/or iii) Doppler echocardiographic data; and/or iv) Ultrafast echocardiographic data.

6. A method of determining pressure difference across a tube arising from fluid flow within the tube, comprising: obtaining partial fluid velocity data from the outlet of the tube; estimating what a full velocity profile is from the partial data; and calculating the pressure difference in dependence on the fluid velocity data and a known fluid density; wherein the tube is a blood vessel, and the fluid flow is in vivo or ex-vivo blood flow, and where the fluid velocity data is obtained by i) 2D Doppler echocardiography; and/or ii) ultrafast 2D echocardiography, where the observation of data along a line is directly taken as the best estimation of the full velocity profile, and as such the pressure difference is given by: $\Delta p_{\mathrm{sAW}}=-\frac{\rho }{2}{\int }_{\lambda }{.Math.v.Math.}^{3}d\lambda \text{/}{\int }_{\lambda }.Math.v.Math.d\lambda$ where v is the fluid velocity along the line λ, and rho is the fluid density.

7. A method of determining pressure difference across a tube arising from fluid flow within the tube, comprising: obtaining partial fluid velocity data from the outlet of the tube; estimating what a full velocity profile is from the partial data; and calculating the pressure difference in dependence on the fluid velocity data and a known fluid density; wherein the tube is a blood vessel, and the fluid flow is in vivo or ex-vivo blood flow, and where the fluid velocity data is obtained by i) 2D Doppler echocardiography; and/or ii) ultrafast 2D echocardiography, where the full velocity profile is estimated as a circular profile divided in two axisymmetric parts, each of them as the extrapolation in the circumferential direction of the data along the observed line, and as such the pressure difference is given by: $\Delta p_{\mathrm{sAW}}=-\frac{\rho }{2}\underset{\lambda =0}{\overset{\lambda =L}{\int }}{.Math.v.Math.}^3.Math.\lambda -L\text{/}2.Math.d\lambda \text{/}\underset{\lambda =0}{\overset{\lambda =L}{\int }}.Math.v.Math..Math.\lambda -L\text{/}2.Math.d\lambda$ where v is the fluid velocity along the line λ of length L, and rho is the fluid density.

8. A method of determining pressure difference across a tube arising from fluid flow within the tube, comprising: obtaining partial fluid velocity data from the outlet of the tube, wherein either the partial data is acquired along several lines at the outlet plane of the tube, each line being defined by intersecting an imaging plane with the outlet plane of the tube or the partial data is acquired at a sparse set of points or lines around the outlet plane, in such a manner that the velocity at the outlet plane can be inferred by interpolation or extrapolation; estimating what a full velocity profile is from the partial data; and calculating the pressure difference in dependence on the fluid velocity data and a known fluid density; wherein the tube is a blood vessel, and the fluid flow is in vivo or ex-vivo blood flow, and where the fluid velocity data is obtained by i) 2D Doppler echocardiography; and/or ii) ultrafast 2D echocardiography, where the pressure difference is given by: $\Delta p_{\mathrm{sAW}}=-\frac{\rho }{2}{\int }_{{\Gamma }_{\mathrm{OUTLET}}}{.Math.v_R.Math.}^2\left(v_R.Math.n\right)\mathrm{dx}\text{/}{\int }_{{\Gamma }_{\mathrm{OUTLET}}}{v}_R.Math.\mathrm{ndx}$ where v.sub.R is the fluid velocity vector reconstructed from the obtained partial fluid velocity, data and rho is the fluid density.

9. A system configured to perform: obtaining fluid velocity data at inlet and outlet planes of a tube; and processing the fluid velocity data to determine: i) a flow rate (Q) of the fluid through the tube; and ii) an advective energy rate (A) of the fluid flow through the tube; the method further comprising calculating the pressure difference in dependence on both of the flow rate (Q), and the advective energy rate (A), without considering a viscous dissipation rate (V) and the kinetic energy (K) of the fluid flow through the tube, wherein: a) the flow rate Q is dependent on a surface integral of the fluid velocity data across either the inlet or outlet plane of the tube, or across any other plane dividing the tube in two sub-sections; and/or b) the advective energy rate A is dependent on a sum of the surface integrals of the time dependent fluid velocity data, or data derived therefrom, across the inlet and outlet planes and the fluid density.

10. The system of claim 9, wherein the pressure difference across the tube is given by one of: $\begin{array}{cc}\Delta p_{\mathrm{cAW}}=-\frac{\rho }{2}\left({\int }_{{\Gamma }_{\mathrm{OUTLET}}}{.Math.v.Math.}^2\left(v.Math.n\right)\mathrm{dx}+{\int }_{{\Gamma }_{\mathrm{INLET}}}{.Math.v.Math.}^2\left(v.Math.n\right)\mathrm{dx}\right)\text{/}{\int }_{{\Gamma }_{\mathrm{OUTLET}}}v.Math.\mathrm{ndx}& A)\end{array}$ where v is a velocity field at a generic voxel, n is the normal direction on inlet and outlet planes, rho is the fluid density respectively, and Γ.sub.INLET and Γ.sub.OUTLET are the respective inlet and outlet planes of the tube; or $\begin{array}{cc}\Delta p_{\mathrm{sAW}}=-\frac{\rho }{2}{\int }_{{\Gamma }_{\mathrm{OUTLET}}}{.Math.v.Math.}^2\left(v.Math.n\right)\mathrm{dx}\text{/}{\int }_{{\Gamma }_{\mathrm{OUTLET}}}v.Math.\mathrm{ndx}& B)\end{array}$ where v is a velocity field at a generic voxel, rho is the fluid density respectively, Γ.sub.OUTLET is a selected plane of the tube, and n is the normal direction on the selected plane, wherein the selected plane is one of the planes at any end of the tube.

11. The system of claim 9, wherein fluid velocity data is not obtained for anywhere else in the tube other than the inlet and outlet planes.

12. The system of claim 9, wherein fluid velocity data is not obtained for anywhere else in the tube other than the selected plane.

13. The system of claim 9, wherein the tube is a blood vessel, and the fluid flow is in vivo blood flow, and wherein the fluid velocity data is obtained from one of: i) 4D phase contrast magnetic resonance imaging data; and/or ii) 2D phase contrast magnetic resonance imaging data; and/or iii) Doppler echocardiographic data; and/or iv) Ultrafast echocardiographic data.

14. A system configured to perform: obtaining partial fluid velocity data from an outlet of a tube; estimating what a full velocity profile is from the partial data; and calculating the pressure difference in dependence on the fluid velocity data and a known fluid density; wherein the tube is a blood vessel, and the fluid flow is in vivo or ex-vivo blood flow, and where the fluid velocity data is obtained by i) 2D Doppler echocardiography; and/or ii) ultrafast 2D echocardiography.

15. The system of claim 14, wherein the partial data is acquired along a line at an outlet plane of the tube, the line being defined by intersecting an imaging plane with the outlet plane of the tube.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) Embodiments of the invention will now be further described by way of example only and with reference to the accompanying drawings, wherein like reference numerals refer to like parts, and wherein:

(2) FIG. 1 is a set of diagrams illustrating various uses of embodiments of the invention;

(3) FIG. 2 is a set of diagrams illustrating measurements collected by embodiments of the invention, and associated mathematical formula for the computation of the pressure drop;

(4) FIG. 3 is a drawing illustrating the comparison of the process of one of the embodiment of the invention with the existing alternatives in the state of the art;

(5) FIG. 4 is a drawing illustrating an aspect of the theoretical basis of the first embodiment of the invention;

(6) FIG. 5 is a block diagram of a first embodiment of the invention;

(7) FIG. 6 is a flow diagram relating to the operation of the first embodiment of the invention;

(8) FIG. 7 is a flow diagram relating to the operation of a second embodiment of the invention;

(9) FIG. 8 is a flow diagram relating to the operation of a third embodiment of the invention;

(10) FIG. 9 is a block diagram of a second embodiment of the invention;

(11) FIG. 10 is a block diagram of a second and third embodiment of the invention; and

(12) FIGS. 11 to 17(a) and (b) are various graphs illustrating results and effectiveness of embodiments of the invention.

DESCRIPTION OF THE EMBODIMENTS

(13) Embodiments of the invention will now be described. Three embodiments will be described, the first embodiment relating to pressure drop estimation in the heart using the full WERP approach which models the pressure drop in detail over a whole length of the blood vessel of interest, the second embodiment relating to the so-called “complete Advective WERP” approach which models the pressure drop via measurements at the inlet and outlet of the vessel only, and the third embodiment being a so-called “simplified advective WERP” method, which makes use of measurements at the outlet of the vessel only. FIG. 2 illustrates the measurements and the computations required for each approach, whereas FIG. 12 gives further details of the anatomical regions of interest (particularly used in the second and third embodiments) (see FIG. 1(A)), and the different data that can be obtained from different imaging modalities (FIG. 1(B)).

(14) In FIG. 1 (A) the definition of anatomical regions where the pressure drops are computed are as follows: Trans Valvular Region (TVR): from the left ventricle outflow tract (LVOT) (Plane 1) to the vena contract (VC) (Plane 2); Ascending Aorta Region (AAR): from the VC to the brachiocephalic artery (Plane 3); Descending Aorta Region (DAR) from the left subclavian artery (Plane 4) to a plane at the same height of the aortic valve plane (Plane 5).

(15) In FIG. 1(B) is shown schematics of the velocity field at the VC acquired at peak systole with different techniques and increasing level of data availability, as follows. Top left: A single velocity value is obtained from continuous 1D Doppler. Top right: A 1D velocity profile along a plane of insonation is available from a 2D color Doppler. Bottom left: A complete 2D velocity profile is available from a 3D color Doppler. Bottom right: The complete velocity vector field over the entire plane, and not simply its projection along the insonation line, is available from 4D PC-MRI.

(16) With respect to the first embodiment, referred to as the full WERP approach, this can be summarised as follows: Given a fluid through a tubular structure, the computation of the difference of pressure between the inlet and outlet is given by the work-energy relative pressure (WERP) method, which is derived from the work-energy principle obtained from the Navier Stokes equation, and computes the pressure difference as an addition of three rates of energy transfer (kinetic, advective and viscous) divided by the net flow through the tube. The kinetic term integrates the rate of temporal change in kinetic energy of all particles within the tubular domain. The second term describes the rate of energy transfer due to the physical movement of a fluid in and out of the domain. And the third term describes the energy dissipated due to viscous friction of the fluid.

(17) With respect to the second embodiment, referred to as the complete advective WERP approach, this can be summarised thus: Given a fluid through a tubular structure, the computation of the difference of pressure between the inlet and outlet is given by the advective WERP equation, which is a simplification of the WERP equation by only taking one of the three additive terms of energy, the advective term. In more detail, this component of energy neglects the contribution of the advected energy from the lateral walls, as velocities perpendicular to the wall are small in the near-wall regions, and the computation of the pressure difference reduces to the integral of the advected energy rate at the inlet and outlet planes of the tubular structure.

(18) Finally, with respect to the third embodiment, being the simplified advective WERP approach, the following is a summary thereof: Given a fluid through a tubular structure, the computation of the difference of pressure between the inlet and outlet is given by the simplified advective WERP equation, which is a simplification of the complete advective WERP equation by only taking the contribution of either the inlet or outlet planes of the tubular domain. The assumption here is that the integral of advected energy rate at the other end of the tubular domain is negligible.

First Embodiment: Full WERP

(19) 2. Methods

(20) Starting from the work-energy principle, we derive the formula for the pressure difference over a vascular segment (Section 2.1). Subsequently, we detail its discrete formulation (Section 2.2) and pre-processing steps (Section 2.3) required to work with 4D PC-MRI data.

(21) 2.1. Pressure Difference from Fluid Work Energy

(22) Pressure differences in a fluid system are related to the kinematics of the flow field. This relationship is described by the well-known Navier-Stokes equations where, in the absence of gravity, variations in pressure are balanced by fluid accelerations and viscous stresses. Using the conservation of mass and momentum for closed systems, the work-energy for an incompressible isothermal Newtonian fluid over a Region Of Interest (ROI) (Ω) with boundary Γ yields,

(23) $\underset{\underset{\frac{\partial }{\partial t}{K}_e}{︸}}{\frac{\rho }{2}\frac{\partial }{\partial t}{\int }_{\Omega }\left(v.Math.v\right)\mathrm{dx}}+\underset{\underset{A_e}{︸}}{\frac{\rho }{2}{\int }_{\Gamma }{.Math.v.Math.}^2\left(v.Math.n\right)\mathrm{dx}}+\underset{\underset{H\left(p\right)}{︸}}{{\int }_{\Gamma }\mathrm{pv}.Math.n\mathrm{dx}}-\underset{\underset{S_e}{︸}}{{\int }_{\Gamma }\mu \left[D\left(v\right).Math.n\right].Math.v\mathrm{dx}}+\underset{\underset{V_e}{︸}}{\frac{\mu }{2}{\int }_{\Omega }D\left(v\right)\text{:}D\left(v\right)\mathrm{dx}}=0,$

(24) where v represents the velocity, ρ the pressure, n is the normal vector on Γ, D(⋅)=[∇(⋅)+∇T(⋅)], and ρ and μ as the fluid density and dynamic viscosity. Here, ∂/∂t K.sub.e is the temporal derivative of the kinetic energy within Ω, A.sub.e the advected energy rate describing the energy transfer due to the physical movement of a fluid in and out of Ω and V.sub.e is the rate of viscous dissipation. H(p) and S.sub.e represent energy inputs to the fluid system, the hydraulic power and the shear energy rate, respectively. Here we assume that the boundary of the Ω can be written as Γ=Γ.sub.i∪Γ.sub.o∪Γ.sub.w, where i, o and w indicate contributions from the vessel inlet, outlet and walls surface. The mathematical details behind the work-energy principle derivation can be found in our published paper Donati et al Non-invasive pressure difference estimation from PC-MRI using the work-energy equation Medical Image Analysis, Vol. 26, pp. 159-1′72, December 2015, the entire contents of which necessary for understanding the present embodiments being incorporated herein by reference.

(25) Starting from this work-energy balance, as a first approximation, we ignore the contribution to the advected energy A.sub.e from the lateral walls Γ.sub.w, as velocities are small in the near-wall regions compared to the core blood flow. Consequently, computations are limited to the inlet and outlet cross-sections, e.g.

(26) $\begin{array}{cc}{A}_e=\frac{\rho }{2}{\int }_{{\Gamma }_l.Math.{\Gamma }_v}{.Math.v.Math.}^2\left(v.Math.n\right)\mathrm{dx}& \mathrm{Equation}\left(1\right)\end{array}$

(27) Furthermore, we assume the pressure to be nearly constant on the inlet and outlet planes, making
H(p)=p.sub.i∫.sub.Γ.sub.iv.Math.ndx+p.sub.o∫.sub.Γ.sub.ov.Math.ndx+∫.sub.Γ.sub.wpv.Math.ndx.  Equation (2)

(28) When little or no compliance is present, |v.Math.n|<<1 on the wall, the global mass balance compatibility condition yields,
ƒ.sub.Γv.Math.ndx=∫.sub.Γ.sub.ov.Math.ndx+∫.sub.Γ.sub.ov.Math.ndx=0,  Equation (3)
letting,
H(p)=ΔpΛ,  equation (4)

(29) where Δp=p.sub.o−p.sub.i is the pressure difference between the outlet and inlet and Λ=∫.sub.Γ.sub.ov.Math.n dx is a term accounting for the flux through surfaces, a term that can be expressed as a function of the inlet surface only by means of Eq. 3.

(30) Regarding the shear energy S.sub.e, we consider the contribution over each boundary segment—inlet, outlet and wall—to be effectively zero. On inlet/outlet planes, this term contributes if there are significant gradients in the direction of the boundary normal. While these gradients can occur particularly in bending or tapering vessels—they are extremely mild and effectively scaled away by the low viscosity of blood. This argument on the flow gradients cannot be assumed near the vessel walls, where a significant wall shear stress is induced. However, as this shear stress is principally orthogonal to the wall velocity (which predominantly dilates in the boundary normal direction), the contribution of these shear stresses to S.sub.e is assumed negligible.

(31) With the assumptions above, the Work-Energy Relative Pressure (WERP) formulation to estimate the pressure difference based on energy contributions yields,

(32) $\begin{array}{cc}\Delta p=-\frac{1}{\Lambda }\left(\frac{\partial }{\partial t}{K}_e+A_e+V_e\right).& \mathrm{Equation}\left(5\right)\end{array}$

(33) From this equation, we observe that all RHS terms are directly derived from flow data, enabling the computation of the pressure difference. However, we also observe that this computation requires that lΛl>0 (e.g. that flow is observed through the vascular segment).

(34) 2.2. Computation from 4D PC-MRI

(35) Let Vt represent the velocity image acquired at time t, Vt(i, j, k) the velocity field evaluated at time t at the voxel (i, j, k) and Δt the discrete time step between two consecutive acquisitions. We discretize derivatives in Eq. 5 using a central finite difference method and estimate the pressure difference between inlet/outlet planes at time t+1/2 as

(36) $\begin{array}{cc}\Delta p^{t+\frac{1}{2}}=-\frac{1}{\Lambda \left(V^{t+\frac{1}{2}}\right)}\left(\frac{K_e\left(V^{t+1}\right)-K_e\left(V^t\right)}{\Delta t}+A_e\left(V^{t+\frac{1}{2}}\right)+V_e\left(V^{t+\frac{1}{2}}\right)\right),& \mathrm{Equation}\left(6\right)\end{array}$

(37) where velocities at t+1/2 are approximated to second order accuracy O(Δt.sup.2) by equation (7)

(38) $V^{t+\frac{1}{2}}=\frac{1}{2}\left(V^t+V^{t+1}\right).$

(39) Computation of the WERP formulation terms is performed by integrating over a voxelized version of Ω, I.sub.ROI. Surface integrals are evaluated on the planes obtained by clipping the 3D mask to define inlet I.sub.in.sup.2D and outlet I.sub.out.sup.2D cross-sections (see FIG. 3) and the normal vectors N.sub.2D(i, j). The discrete terms are then estimated from the image-based velocity field as,

(40) $\begin{array}{cc}\Lambda \left(V\right)=\mathrm{dS}\underset{\left(i,j\right)\in I_{\mathrm{out}}^{\mathrm{ID}}}{.Math.}{ℳ}_{2D}\left(V\right)\left(i,j\right).Math.N_{2D}\left(i,j\right),K_e\left(V\right)=\rho \mathrm{dV}\underset{\left(i,j,k\right)\in I_{\mathrm{ROI}}}{.Math.}{.Math.ℳ\left(V\right)\left(i,j,k\right).Math.}^2,A_e\left(V\right)=\frac{\rho }{2}\mathrm{dS}\underset{\left(i,j\right)\in I_m^{2D}.Math.I_{\mathrm{int}}^{2D}}{.Math.}{.Math.{ℳ}_{2D}\left(V\right)\left(i,j\right).Math.}^2.Math.\left({ℳ}_{2D}\left(V\right)\left(i,j\right).Math.N_{2D}\left(i,j\right)\right),V_e\left(V\right)=\frac{\mu \mathrm{dV}}{2}\underset{\left(i,j,k\right)\in I_{\mathrm{ROI}}}{.Math.}𝒟\left(V\right)\left(i,j,k\right)\text{:}𝒟\left(V\right)\left(i,j,k\right).& \mathrm{Equation}\left(8\right)\end{array}$

(41) where dS=Δx.sup.2 and dV=Δx.sup.3 are the pixel surface and voxel volume, respectively, based on the voxel length Δx. The discrete evaluation of all the contributions relies on the definition of the approximated velocity fields M(V) and M.sub.2D(V), obtained through averaging over the 3D mask and on the 2D planes defined above,

(42) $\begin{array}{cc}\left(V\right)\left(i,j,k\right)=\frac{1}{2\left(1+q\right)}\underset{n=o}{\overset{q}{.Math.}}(V\left(i+{\delta }_{n1},j+{\delta }_{n2},k+{\delta }_{n3}\right)+V\hspace{1em}\left(i-{\delta }_{n1},j-{\delta }_{n2},k-{\delta }_{n3}\right)),{2D}_{}\left(V\right)\left(i,j\right)=\frac{1}{2.Math.\max \left(1,q\right)}\underset{n=o}{\overset{\max \left(0,q-1\right)}{.Math.}}\left(V\left(i+{\delta }_{n1},j+{\delta }_{n2}\right)+V\left(i-{\delta }_{n1},j-{\delta }_{n2}\right)\right).& \mathrm{Equation}\left(9\right)\end{array}$

(43) In the above, δ.sub.ij is the Kronecker delta and q is a parameter used to smooth the underlying data based on O(Δx.sup.2) approximations to the velocity value (see FIG. 4, which shows a schematic representation of the finite differences central stencils. Velocity field V evaluation at P=(i,j,k) through operator M(V) using standard (a) and filtered approach (b) and velocity field V derivative evaluation at P=(i,j,k) through operator D(V) using standard (c) and filtered central differences approach (d)).

(44) In the above if q=0,M(V)(i,j,k)=V(i,j,k) and M.sub.2D(V)(i,j)=V(i,j) return the velocity measured at the voxel (i, j, k) and (i, j), respectively. Alternatively, if q=3, the measurement of the velocity field is taken as a weighted sum of O(Δx.sup.2) approximations based on neighboring voxel measurements, effectively averaging out potential artefacts due to noise.

(45) Similarly, in Eq. 8, the discrete tensor D(V) is calculated as,
D(V)(i,j,k)=(G(V)(i,j,k)+G(V)T(i,j,k)),  Equation (10)

(46) where G(V) is a velocity gradient tensor defined as,

(47) $\begin{array}{cc}{\mathrm{mn}}_{}\left(V\right)\left(i,j,k\right)=\frac{1}{2.Math.\max \left(1,q\right)}\underset{x=0}{\overset{q}{.Math.}}\left(1-{\delta }_{\mathrm{rn}}\right){\tilde{D}}_n^r{V}_m\left(i,j,k\right),\mathrm{where}& \mathrm{Equation}\left(11\right)\\ {\tilde{D}}_n^r{V}_m\left(i,j,k\right)=\frac{1}{2}\left(D_n{V}_m\left(i+{\delta }_{r1},j+{\delta }_{r2},k+{\delta }_{r3}\right)+D_n{V}_m\left(i-{\delta }_{r1},j-{\delta }_{r2},k-{\delta }_{r3}\right)\right),\mathrm{and}& \mathrm{equation}\left(12\right)\\ D_n{V}_m\left(i,j,k\right)=\frac{\begin{array}{c}{V}_m\left(i+{\delta }_{n1},j+{\delta }_{n2},k+{\delta }_{n3}\right)-\\ V_m\left(i-{\delta }_{n1},j-{\delta }_{n2},k-{\delta }_{n3}\right)\end{array}}{2\Delta x}.& \mathrm{Equation}\left(13\right)\end{array}$

(48) Again, if q=0 velocity gradients are approximated by second order central differences centred at the voxel (i, j, k). Imposing q=3, a filtered approach is adopted, where the velocity derivative is approximated using weighted average of derivatives computed with second order central differences at neighboring voxels, therefore reducing noise contamination (see FIG. 2).

(49) 2.3. Required Pre-Processing

(50) Prior to application in a clinical setting, a number of pre-processing steps are required. Field inhomogeneities and eddy currents may be corrected (1) using the pre-processing tools outlined in Bock et al. (2011). Subsequently, a binary mask I.sub.ROI is defined (2), based on a thresholding of the velocity magnitude calibrated by the maximum velocity V.sub.max including voxels with a velocity magnitude greater than V.sub.max, with S being the segmentation thresholding parameter). Inlet and outlet points are manually selected by the user (3) depending on the clinical problem under investigation. A skeletonisation of the binary mask is then used to define the inlet and outlet planes perpendicular to the vessel (4). As a result of this process, the binary masks of the raw 3D image and the inlet/outlet planes needed for the WERP computation are defined. Within this work, the image acquisition process was mimicked in silico for the validations tests presented in Results 3.1 and 3.2 and 3.3 in Donati et al Non-invasive pressure difference estimation from PC-MRI using the work-energy equation Medical Image Analysis, Vol. 26, pp. 159-172, December 2015. Simulated PC MRI images were subsequently processed following steps (2)-(4) prior to application of the WERP method.

(51) 2.4 System Description

(52) FIG. 5 illustrates an example imaging system that is able to apply the above described techniques to calculate blood flow through a blood vessel using the full WERP method. Here, a 4D phase contrast magnetic resonance imaging (4D PC MRI) system 60 is provided, comprising MR imaging coils within which the subject is located, controlled by an MRI imaging control system 58, including an MRI controller processor 50. The MRI imaging control system 58 including the MRI controller processor 50 function in a conventional manner to allow 4 dimensional phase contrast magnetic resonance image data to be obtained, for example of an internal blood vessel of the subject the pressure difference across it is desired to know. For example, it may be desirable to measure transvalvular pressure drops (TPD) along the transvalvular region in the heart, between the left ventricular outflow tract and the vena contracta. Of course, other blood vessels may also be monitored, as desired. For example, it may be possible to take into account secondary branches.

(53) The 4D PC MRI system 58, 60, collects 4D PC MRI data 52 of the imaging subject, which is saved for further processing and analysis. As explained previously with respect to FIG. 1, the complete velocity vector field over the entire imaged plane, and not simply its projection along the insonation line, is available from 4D PC-MRI.

(54) The 4D PC MRI system 58, 60, also includes a WERP computation program 54, which acts to process the 4D PC MRI data 52 as described above in accordance with the WERP processing method described, to obtain WERP pressure estimation data 56. The WERP processing method computes the pressure difference along the vessel as an addition of three rates of energy transfer (kinetic, advective and viscous) divided by the net flow through the vessel, as described above. The resulting WERP pressure estimation data 56 gives an estimate of pressure at each point along the vessel between the measured inlet and outlet planes, taking into effect the vessel walls between the inlet and outlet. It provides a complete solution for every point within the vessel between the inlet and outlet.

(55) FIG. 6 is a flow diagram summarising some of the important elements of the full WERP process. In step 6.2 the 4D phase contrast MRI data relating to the blood vessel of the subject for which pressure is to be found is obtained. Note that this may be obtained “live” i.e. just before the calculation is performed, or may be recorded or stored data from a previous scan. In this respect, there is no need for the human subject to be present in order for the WERP process to be performed to find the pressure drops, as the process can be run on recorded or stored MRI data.

(56) At step 6.4, the 4D PC MRI data is segmented, as shown in FIG. 3, to define the region of interest corresponding to the blood vessel for which the pressure drop is to be found. Then, at step 6.6, and again as shown in FIG. 3, the vessel inlet and outlet planes are defined.

(57) Next, at step 6.8, for each WERP formulation term the image based velocity field from the MRI data is integrated over the defined region of interest. The precise calculations performed are detailed in Equation 8 and the related equations above. The result are a number of individual terms that can then be combined to find the whole pressure drop over the defined region of the vessel, using Equation 5 above, as shown at 6.12. The calculated WERP pressure drop data is then saved, and output to the user, at step 6.14, for example by being displayed on a screen (not shown).

(58) Second and Third Embodiments: Complete Advective WERP and Simplifed Advective WERP

(59) The second and third embodiments relate to simplifications of the full WERP process, that mean that less complex imaging modalities which collect less information may be used. In particular, as described below, the complete advective WERP approach may be performed using any modality that renders velocity data in two anatomical planes, for example 2D PC MRI or 3D Doppler ECG, whereas the simplified advective approach may be used with any modality that renders velocity data in a single anatomical plane, for example 2D PC MRI, 3D Doppler ECG or 2D Doppler ECG data. Accessing this partial information (i.e. velocity data in a single anatomical plane) is feasible by Doppler ECG imaging apparatus that are very common in many clinical environments.

(60) There follows a discussion of the approaches, and presentation of various results and comparisons with the Bernoulli based approaches of the prior art.

(61) As shown in the first embodiment of the invention, using the WERP formulation the contributions to the pressure drop due to temporal acceleration (kinetic drops), spatial transport of momentum (advective drops) and friction losses (viscous drops) are revealed to better understand the hemodynamics of the human aorta towards improved stratification of patients.sup.22.

(62) Building on this, we conduct a thorough comparison between drops computed with WERP and Bernoulli formulation on 3D PCMRI (or 4D flow MRI) data, for the assessment of severity of AS on a cohort of 32 patients with BAV. Studying the correlation between non-invasive drop estimates obtained with the different approaches, a consistent overestimation was observed in the Bernoulli formulation. Further investigation demonstrated that this discrepancy is the result of the fact that a Bernoulli approach does not take into account information about the velocity profile. Examining the WERP prediction of the TPD, a clear dependence on the advective drop was observed, particularly in stenosed patients. Utilizing this fact, we propose a novel simplification of the WERP formulation, which improves the accuracy of computed TPD but limits the dependence on comprehensive flow data, enabling estimation from either 2D PCMRI or 3D Doppler echocardiography data. The applicability of the method to 2D Doppler echocardiographic data is also tested.

(63) Material and Methods

(64) Patient Data

(65) For these embodiments, a cohort of 32 subjects with different degree of aortic valve stenosis is selected. Each subject underwent two cardiovascular magnetic resonance scans—one on a 1.5 T system (Avanto, Siemens, Erlangen, Germany) for anatomical imaging, the second on a 3 T system (Trio, Siemens, Erlangen, Germany) for 4D PCMRI assessment, both using a 32-channel cardiac coil. All images were electrocardiogram (ECG)-gated.

(66) Subjects were divided between Group I (n=20) and Group II (n=12) according to an assessment of the mean TPD following current clinical guidelines.sup.23, diagnosing mild stenosis if mean TPD>20 mmHg. The Bernoulli principle using peak velocity values across the valve was used for the computation of these pressure values. Aortic dimensions and hemodynamics data are shown in Table 1.

(67) TABLE-US-00001 TABLE 1 Aortic dimensions and flow hemodynamics in n = 32 patients divided in two groups based on the mean systolic pressure drop: Group I (Δp ≤ 20 mmHg, n = 20) and Group II (Δp > 20 mmHg, n = 12). Group I Group II Male (%) 35% 91% Age 28.2 ± 14.1 38.8 ± 20.2 Aortic Diameters/BSA [mm/m.sup.2] Left ventricle outflow tract 13.4 ± 2.6  15.4 ± 3.8  Aortic valve 14.5 ± 1.7  17.6 ± 4.0  Brachiocephalic artery 14.5 ± 1.5  20.7 ± 3.9  Left subclavian artery 12.2 ± 1.2  13.8 ± 1.8  Mid descending aorta 11.1 ± 1.1  12.0 ± 1.4  Values are mean ± std

(68) Pre-Processing and Definition of Anatomical Regions

(69) 4D PC-MRI images field inhomogeneities and eddy currents.sup.24-26 were corrected using the pre-processing tools developed by Bock et al..sup.27. The lumen of the left ventricle and aorta were identified using a thresholding criterion based on the peak velocity magnitude, defining a binary mask. A skeletonization algorithm.sup.28 was used to extract the centerline of the aorta and its perpendicular planes.

(70) Pressure drops were computed over three anatomical regions illustrated in FIG. 1A, whereas FIG. 2 illustrates the measurements and the computations required for each approach, whereas FIG. 12A. The assessment of AS was then based on the TPD computed along the TransValvular Region (TVR), between the Left-Ventricular Outflow Tract (LVOT) (Plane 1) and the Vena Contracta (VC) (Plane 2), which are reported in Results, below. The LVOT plane was located 12 mm before the aortic valve plane, following the definition used by Garcia et al..sup.29, and the VC was detected from the medical image as the plane containing the peak velocity magnitude, which is also the plane of maximum narrowing of the aortic valve jet.

(71) For completeness, pressure drops evaluated in the Ascending Aorta Region (AAR)—from the VC (Plane 2) to the brachiocephalic artery (Plane 3)—and the Descending Aorta Region (DAR)—from the left subclavian artery (Plane 4) to a plane at the same height of the aortic valve plane (Plane 5)—are also included in Supplemental Material B.

(72) Simulated 3D, 2D and continuous (M) Doppler Echocardiography

(73) In order to remove any observer and inter-modality dependence, echocardiography data was derived by sampling 4D PC-MRI in each subject. Idealized conditions were taken: a perfect alignment between the blood flow jet and the line of insonation, with no shadowing. These conditions enable an optimal comparison between the different mathematical formulations introduced in next section.

(74) Echocardiographic data was simulated in Planes 1 and Plane 2 for the computation of the TPD (see FIG. 2 illustrates the measurements and the computations required for each approach, whereas FIG. 12B, whereas FIG. 2 illustrates the measurements and the computations required for each approach, whereas FIG. 2). Original PC-MRI data was linearly interpolated onto a grid of 1 mm×1 mm sample points in each plane. 3D Doppler echocardiographic acquisitions were defined by projection of the velocity along the direction of insonation, taking into account the funneling effect of the probe. To achieve this, the probe location was simulated 10 cm upstream of the VC in the direction of the aortic jet flow, where the aortic valve jet direction was defined as the direction of the velocity vector at the pixel with maximum velocity magnitude, thus defining an idealized 2D velocity profile by color Doppler. Similarly, a set of 1D velocity profiles from 2D color Doppler acquisition were defined by the intersection of the previously projected velocity field with hypothetical insonation planes—since velocity profiles are non-axisymmetric, a total of 12 profiles containing the peak velocity and with arbitrarily oriented lines (with increments of 15°) were generated in each case (see Supplemental Material D for an illustration of the 2D and 1D profiles in each case). Finally, continuous (1D) Doppler echocardiographic data was simulated using the magnitude of the peak velocity pixel projected in the aortic jet flow direction at each time point.

(75) Non-Invasive Pressure Drop Estimates

(76) To assess the patient-specific severity of AS, we compared TPD obtained using a range of approaches ranging from the most complete (WERP) to the simplest (Bernoulli) formulation. Taking a number of simplifications, the derived approaches can be all obtained from the Navier-Stokes equation for Newtonian isothermal fluids, as extensively reported in Supplemental Material A. For clarification purposes, all the methods compared and the assumptions made to obtain the different formulations are schematically presented in FIG. 2.

(77) The Bernoulli principle is widely accepted in the assessment of AS. Obtained from the Navier-Stokes momentum equation by considering the flow along a streamline, it neglects any contribution from the kinetic and viscous pressure components.sup.30,31. A series of additional assumptions yields the definition of an approximate pressure drop Δp.sub.sB in mmHg as,
Δp.sub.sB=4v.sub.MAX.sup.2,  (Equation 101)

(78) where v.sub.MAX is the peak velocity at the VC estimated from the 4D PC-MRI image. Factor 4 allows conversion of pressure from Pa to mmHg (133 Pa=1 mmHg), if blood density ρ=1060 kg/m.sup.3 is assumed. In this work this was referred to as the simplified Bernoulli (sB) formulation.sup.32.

(79) Building on this approach and using the same imaging data, we also accounted for the proximal velocity acquired at the LVOT v.sub.PROX, thus defining the corrected Bernoulli (cB) pressure drop Δp.sub.cB=4(v.sub.MAX.sup.2−v.sub.PRX.sup.2).sup.33-35. Clinical guidelines make use of the corrected formulation in patients evaluation when v.sub.PROX>1.5 m/s or v.sub.MAX3<0.0 m/s. It follows that with Bernoulli-based formulations just one single velocity value or two are needed to estimate the pressure drops. This prompted applicability to continuous 1D Doppler echocardiography, thus enabling definition of pressure drops Δp.sub.sB.sup.D1D and Δp.sub.cB.sup.D1D).

(80) The availability of comprehensive velocity fields from 4D PC-MRI also enables computation of the pressure drop using fewer assumptions on the flow. Assuming negligible compliant effects and nearly constant pressure at the inlet and outlet planes of a tubular region, the WERP method.sup.21 approximates the Navier-Stokes equation to compute the pressure drop from 4D PC-MRI data Δp.sub.w as,

(81) $\begin{array}{cc}\Delta p_w=-\frac{1}{Q}\left(\frac{\partial K}{\partial t}+A+V\right),& \left(\mathrm{Equation}102\right)\end{array}$

(82) where Q is the flow rate computed at the outlet, ∂K/∂t is the temporal derivative of the kinetic energy within the vascular region, A is the advected energy rate describing the energy transfer due to the physical movement of a fluid in and out of the domain, and V is the rate of viscous dissipation describing energy losses due to friction. These quantities were estimated directly from the image as described by Donati et al..sup.21.

(83) The sum of energy contributions of Equation 102 allows separation of the components, yielding to the definition of the kinetic pressure drop Δp.sub.KW=−1/Q(∂K/∂t) the advective pressure drop Δp.sub.AW=−A/Q and the viscous pressure drop Δp.sub.AW=−V/Q. As the computation of K and V requires the whole velocity field (volume integrals), contrariwise to A and Q (surface integrals), we looked at approximations eliminating these to reduce the data dependence. This is the reason why the rest of this section focuses on the advective term of the TPD.

(84) The complete advective WERP (cAW) approach estimates the pressure drop Δp.sub.cAW from the velocity field extracted at the inlet and outlet planes (This is analogous to the previously defined Δp.sub.AW. The notation Δp.sub.cAW stresses that drops in the complete advective WERP formulation are evaluated using both inlet and outlet planes.). It takes into account the three-dimensionally encoded velocity field over the 2D region instead of the single peak velocity value of sB.

(85) Analogously to the sB approach, the assumption of outlet velocities much larger than inlet velocities triggers the definition of the simplified advective WERP (sAW) pressure drop from the 4D PC-MRI data Δp.sub.sAW, which only accounts for the contribution to the advective energy at the outlet plane.

(86) However, as the availability of 4D PC-MRI velocity fields is limited in clinical practice, the above mentioned WERP-based formulations can also be applied to the more widely used echocardiographic data. Consequently, by using the advective WERP approaches to 3D Doppler echocardiographic images, we estimated drops Δp.sub.sAW.sup.D3D and Δp.sub.cAW.sup.D3D based on the velocity field acquired along the direction of the beam in planes perpendicular to the blood jet at the VC and LVOT.

(87) The sAW approach in particular can be further simplified to compute the advected energy rate A by only relying on velocity values at the VC sampled along a line, therefore yielding the pressure drop Δp.sub.sAW,1D.sup.D2D from 2D color Doppler echocardiographic data.

(88) Within this work we mainly focused on mean pressure drops, as clinically accepted guidelines for the echocardiographic assessment of the AS are based on these. Trivially, the mean pressure drops (or mean pressure drop components) Δp were estimated as the temporal average of the instantaneous drops during systole Δp(t). For completeness, the peak pressure drop (or peak pressure drop component) Δp.sub.MAX was instead computed as the maximum absolute value of the instantaneous drops.

(89) Results

(90) Analysis of the Pressure Components of the TPD

(91) The impact of the kinetic, advective and viscous pressure components to the TPD was analyzed for both groups of patients by using the WERP formulation. We computed the temporal transient of the drops over systole, together with their peak values (see FIG. 11), with peak values of the kinetic pressure drops evaluated on acceleration and deceleration phases separately. Clear differentiation between groups is revealed for the peak advective (P<0.001) and viscous drops (P<0.001), unlike for the kinetic drops during acceleration (P<0.01).

(92) Subjects in Group II have an average advective TPD (absolute peak of 28.74±6.01 mmHg, accounting for 96.9% of the total TPD on average, range 91.5% to 98.5%) highly dominant over the kinetic term by almost one order of magnitude (absolute peak of 4.76±3.35 mmHg during acceleration and of 4.63±2.96 mmHg during deceleration), and over the viscous term by over two orders of magnitude (absolute peak of 0.14±0.07 mmHg). Prevalence of the advective component is also shown in Group I (accounts for 86.2% of the total TPD on average, range 60.5% to 98.7%) although to a lesser extent (absolute peak means of 4.83±3.10 mmHg, 2.51±1.67/0.97±0.08 mmHg and 0.02±0.01 mmHg for the advective, kinetic during acceleration/deceleration and viscous components, respectively).

(93) Comparison of the TPD by WERP and Bernoulli Formulations

(94) We compare the TPD computed from the 4D PC-MRI images with the cAW formulation against results obtained with sAW, sB and cB formulations. We compared the mean pressure drops estimated in the TVR with the different approaches, reporting the linear regression and correlation coefficient on Group I and II, separately. FIG. 12 reveals a consistent overestimation of the drops computed with sB (linear regression slopes 0.492 for Group I and 0.367 for Group II, with average overestimation of 99.1%, range of 49.2% to 145.1%), cB (0.573 for Group I and 0.291 for Group II, with average overestimation of 71.3%, range of 14.6% to 141.1%) and sAW (0.887 for Group I and 0.817 for Group II) formulations. Higher correlations were reported between the cAW approach and each of the compared methods for patients in Group I.

(95) Improved TPD Computation: Understanding the Velocity Field

(96) Finally, we explored the reasons behind the discrepancies highlighted in FIG. 12 between the different methods, at the same time investigating the applicability of the WERP-based formulations to different imaging data. To this end, the advective TPD from the original 4D PC-MRI velocity fields (Δp.sub.sAW) were compared against those computed using simulated and idealized 3D Doppler (Δp.sub.sAW.sup.D3D) and 2D Doppler (Δp.sub.sAW.sup.D2D) echocardiographic data. Additionally, we reported TPD estimates from the sB approach to simulated continuous Doppler echocardiographic images (Δp.sub.sB.sup.DID).

(97) FIG. 13 illustrates the overestimation obtained with the sB formulation (regression slope of 0.496) and with the sAW formulation applied to 1D velocity profiles obtained from an idealized color Doppler acquisition (regression slope of 0.755 when comparing against averaging results from 12 velocity profiles in each case). Variability is introduced by the specific choice of the profile, determined by the insonation plane. A higher correlation is achieved with an idealized color Doppler 3D acquisition using the sAW formulation (regression slope of 1.076).

(98) This section is complemented in Supplemental Material C with an in-silico analysis of three velocity profiles in an idealized stenosis, which demonstrates that the overestimation observed with the Bernoulli formulation is uniquely driven by the shape of the velocity profile.

(99) Discussion

(100) Herein we report the existence of an overestimation of the TPD by the Bernoulli principle, explain its fundamental cause, and propose the formulation to correct it by accounting for the velocity profile in the cross section of the blood jet at the point of the vena contracta.

(101) The TPD is driven by the forces that spatially accelerate the flow ejected by the ventricle through a narrow orifice regardless of the other two components of the pressure drops, as illustrated in FIG. 11. Referring to the Navier-Stokes additive terms that contribute to pressure, the transport of momentum (i.e. the advective forces) is the principal contributor to the TPD, with a secondary role played by temporal accelerations (or kinetic drops), and minor impact of viscous dissipative stresses. This finding confirms the sensible choice of the Bernoulli principle to quantify the pressure drop from continuous Doppler recordings adopted in clinical guidelines.sup.2,13,36,37, since the Bernoulli principle is the simplification of the physics of flow through a pipe accounting only to the advective forces.

(102) However, the results described herereveal an overestimation of Bernoulli-based drop computations compared to results obtained with a more complete formulation of the advective drop through the WERP method. The cause of this fundamental bias is that a Bernoulli principle relies on the simplification of the vessel as a single streamline.sup.38-40, as if the vessel was a one-dimensional pipe. But blood flows through 3D vessels, and therefore the computation of the advective energy, and the advective drop, should account for the blood velocity in a cross-section of the vessel, not for only a single peak velocity value. Only if blood velocity profiles were flat, with all particles of the cross section of a vessel having the same velocity, the Bernoulli principle could be adopted without any loss of accuracy (see Supplemental Material C below for a quantitative description of the error in three idealized velocity profiles).

(103) The analysis of our cohort of 32 subjects reveals a large variability in the morphology of the velocity profiles, as illustrated in FIG. 17. As a consequence, the amount of overestimation of the Bernoulli principle is heterogeneous: in the group of stenosed cases Bernoulli overestimated the TPD of 99.1% on average (range of 49.2% to 145.1%) with the simplified and of 71.3% (range of 14.6% to 141.1%) with the corrected formulation, leading to a poor linear regression fitting (see FIG. 12, panels A & B). As an additional reference, in-silico results of Supplemental Material C reveal how the shape of a perfect paraboloid with β=2 introduces an overestimation of 100.9%. The magnitude of proposed correction is therefore quite large and variable, and justifies therefore the need for further research to quantify its impact in the stratification of subjects with valve stenosis.

(104) Our mathematical analysis reveals that, as in the widely adopted Bernoulli formulations, the adoption of sensible assumptions simplify the computation of the advective TPD to only account for the velocity at the inlet and outlet of the vascular domain. This finding is quite relevant for the clinical translation of our findings, since velocity data is only needed at two planes of the vascular anatomy, and not in the entire lumen or ventricular blood pool. Even more, the need for data in the proximal (inlet) region is acknowledged in guidelines only for cases with small level of stenosis, where the corrected Bernoulli formulation should be used.sup.23,34. Our results confirm the improved accuracy of this strategy for the control cases, but illustrate an increment in the discrepancy with the adoption of a corrected Bernoulli formulation in the mild-to-severe disease cases (drop of correlation and slope in Group II against Group I with correction, compare panels A and B in FIG. 12). These findings suggest that, as currently done in most practical clinical situations, accounting only for the outlet velocity profile (at the point of the vena contracta), neglecting the proximal velocities, is a sensible methodological strategy. The goal of a non-invasive method to stratify VS is then set in the robust and accurate assessment of the advective energy at the point of the vena contracta.

(105) Our results then reveal how this computation of the advective drop depends on the amount, and quality, of data available about the velocity at the point of the VC. We provide the mathematical formulation adapted to different sources of blood velocity data, and analyze the existence of fundamental biases between them. Our results reveal that idealized 3D echocardiographic data, offering a 2D velocity profile at the VC with artifacts from the funneling effect and from the projection of the velocity alongside the line of insonation, introduce a small bias (slope of the regression line of 1.076, see FIG. 13). Accounting for the same artifacts, but with only the information of a part of the 2D velocity profile along one line of insonation as with a 2D echo probe, advective energy is overestimated (slope of the linear regression of 0.755, see FIG. 13). And if only a peak value is available, and the formulation then reduces to the Bernoulli equation, the overestimation is the largest (slope of 0.496, see FIG. 13). The translation of the theoretical finding is thus feasible to a wide range of imaging acquisition protocols and modalities, and further research is needed to define the optimal strategy to control the location of the vena contracta to be imaged (point of maximum advective drop), to identify the direction of the jet and its perpendicular plane, and to maximize the amount and quality of velocity data at that plane.

(106) Our findings also offer a plausible explanation for the overestimation of the TPD using continuous Doppler echocardiography compared to catheter recordings, currently attributed to the pressure recovery after the narrowing of the jet across the valve.sup.5,6,9. The fundamental cause is the intrinsic assumption of the Bernoulli principle of the reduction of the vessel flow to a single streamline. This work demonstrates that, accounting for the physical principles that govern human hemodynamics (the Navier-Stokes momentum equation) in order to compute the pressure drop that best explains the changes observed in a dense velocity field, the simplifications made by the Bernoulli equation introduce an overestimation of the TPD.

(107) Supplemental Material

(108) A. Mathematical Details Behind the WERP Pressure Drops Estimates

(109) The original form of the WERP formulation is defined from the Navier-Stokes equation based on the work-energy principle—as presented by Donati et al.—and estimates the pressure drop Δp.sub.W across the vascular region c with inlet plane Γ.sub.INLET and outlet planer OUTLET defined from the lumen segmented from the 4D PCMRI image as,

(110) 0$\begin{array}{cc}\Delta p_w=-\frac{1}{Q}\left(\frac{\partial K}{\partial t}+A+V\right).& \left(\mathrm{Equation}A\mathrm{.1}\right)\end{array}$

(111) The blood flow rate Q, the kinetic energy K, the advective energy rate A and the viscous dissipation rate V can be evaluated by solving numerical surface and volume integrals as,

(112) $\begin{array}{cc}Q={\int }_{{\Gamma }_{\mathrm{OUTLET}}}v.Math.\mathrm{ndx}=-{\int }_{{\Gamma }_{\mathrm{INLET}}}v.Math.\mathrm{ndx},K=\frac{\rho }{2}{\int }_{\Omega }\left(v.Math.v\right)\mathrm{dx},A=\frac{\rho }{2}\left({\int }_{{\Gamma }_{\mathrm{INLET}}}{.Math.v.Math.}^2\left(v.Math.n\right)\mathrm{dx}+{\int }_{{\Gamma }_{\mathrm{OUTLET}}}{.Math.v.Math.}^2\left(v.Math.n\right)\mathrm{dx}\right),V=\frac{\mu }{2}{\int }_{\Omega }D\left(V\right)\text{:}D\left(v\right)\mathrm{dx},& \left(\mathrm{Equation}A\mathrm{.2}\right)\end{array}$

(113) where v is the three-dimensional time-dependent velocity field at the generic voxel, n is the normal direction on the inlet/outlet plane, ρ=1060 kg/m.sup.3 and μ=0.004 Pa.Math.s are the blood density and dynamic viscosity, respectively, and D(⋅)=[∇(⋅)+∇(∩).sup.T]. Using separation of the pressure components the complete advective pressure drop evaluated using the WERP method Δp.sub.cAW=−A/Q yields from Equation A.2,

(114) $\begin{array}{cc}\Delta p_{\mathrm{cAW}}=-\frac{\rho }{2}\left({\int }_{{\Gamma }_{\mathrm{OUTLET}}}{.Math.v.Math.}^2\left(v.Math.n\right)\mathrm{dx}+{\int }_{{\Gamma }_{\mathrm{OUTLET}}}{.Math.v.Math.}^2\left(v.Math.n\right)\mathrm{dx}\right)\text{/}{\int }_{{\Gamma }_{\mathrm{OUTLET}}}v.Math.\mathrm{ndx},& \left(\mathrm{Equation}A\mathrm{.3}\right)\end{array}$

(115) therefore reducing the drops computation to surface integrals on the inlet and outlet planes and making it applicable to 2D PCMRI or 3D Doppler echocardiographic data.

(116) Equation A.3 can be further simplified by assuming outlet velocities much larger than inlet velocities (which is likely to hold in the transvalvular region defined from the LVOT to the VC, especially in stenosed cases) as,

(117) $\begin{array}{cc}\Delta p_{\mathrm{sAW}}=-\frac{\rho }{2}{\int }_{{\Gamma }_{\mathrm{OUTLET}}}{.Math.v.Math.}^2\left(v.Math.n\right)\mathrm{dx}\text{/}{\int }_{{\Gamma }_{\mathrm{OUTLET}}}v.Math.\mathrm{ndx}.& \left(\mathrm{Equation}A\mathrm{.4}\right)\end{array}$

(118) The sAW approach can be further reduced to estimate the advective energy rate from velocity values along of the VC along a single line, not in the complete perpendicular plane, thus enabling applicability to 2D color Doppler echocardiographic images. Equation A.4, by replacing the surface integrals at the outlet plane for line integrals along the line λ defined by intersecting the hypothetical insonation plane with the outlet plane of the aortic lumen plane, and by considering the fact that velocity values are already projected in the direction of the line of insonation, can be rearranged as,

(119) $\begin{array}{cc}\Delta p_{\mathrm{sAW}}=-\frac{\rho }{2}\underset{\lambda }{\int }{.Math.v.Math.}^{2}d\lambda \text{/}\underset{\lambda }{\int }.Math.v.Math.d\lambda & \left(\mathrm{Equation}A\mathrm{.5}\right)\end{array}$

(120) This enables the use of the WERP formulation to compute advective drops based on 2D color Doppler echocardiography.

(121) It is worth noting that advective WERP and Bernoulli formulations are similar—as they both characterize the pressure drop using advective effects—and the mathematical link between them is here explained. In the WERP approach, the blood flow rate Q can be indifferently estimated at the inlet or outlet planes defined from the image data as,
Q=∫.sub.Tv.Math.ndx=v.sub.MAXΨ  (Equation A.6)

(122) Here, v.sub.MAX is the maximum velocity at the inlet/outlet plane and Ψ=∫.sub.ΓΦdx, where Φ is the normalized shape function in the normal direction for the inlet/outlet velocity profile. Using the cAW formulation in Equation A.3 leads to,

(123) $\begin{array}{cc}-\Delta p_{\mathrm{cAW}}Q+\frac{\rho }{2}{\int }_{{\Gamma }_{\mathrm{INLET}}}{v}_{\mathrm{MAX},\mathrm{INLET}}^3\mathrm{dx}+{\int }_{{\Gamma }_{\mathrm{OUTLET}}}{v}_{\mathrm{MAX},\mathrm{OUTLET}}^3\mathrm{dx}=0.& \left(\mathrm{Equation}A\mathrm{.7}\right)\end{array}$

(124) If we assume velocity at the planes mainly aligned to the planes normal n, substitution of Equation A.6 (selectively evaluated at the inlet/outlet planes) into Equation A.7 yields,

(125) $\begin{array}{cc}\Delta p_{\mathrm{cAW}}=\frac{\rho }{2}\left(v_{\mathrm{MAX},\mathrm{OUTLET}}^2{Z}_{\mathrm{OUTLET}}-v_{\mathrm{MAX},\mathrm{INLET}}^2{Z}_{\mathrm{INLET}}\right),& \left(\mathrm{Equation}A\mathrm{.8}\right)\end{array}$
where Z.sub.INLET=∫.sub.Γ.sub.INLET Φ.sub.INLET.sup.3dx/Ψ.sub.INLET and Z.sub.OUTLET=∫.sub.Γ.sub.OUTLET Φ.sub.OUTLET.sup.3dx/Ψ.sub.OUTLET are functions depending on the normalized profile shape only. Consequently, in the hypothesis of a flat velocity profile (i.e. v.sub.OUTLET=V.sub.MAX,OUTLET and v.sub.INLET=v.sub.MAX,INLET as in Bernoulli based formulations) and with blood density ρ=1060 kg/m.sup.3, Equation A.8 simplifies to the cB formulation Δp.sub.cB=4(v.sub.OUTLET.sup.2−v.sub.INLET.sup.2). From Equation A.7, in the hypothesis of v.sub.OUTLET>>v.sub.INLET the sAW formulation in Equation A.4 can be expressed as,

(126) $\begin{array}{cc}\Delta p_{\mathrm{sAW}}=\frac{\rho }{2}{v}_{\mathrm{MAX},\mathrm{OUTLET}}^2{Z}_{\mathrm{OUTLET}},& \left(\mathrm{Equation}A\mathrm{.9}\right)\end{array}$

(127) that assuming a flat velocity profile yields again the sB formulation Δp.sub.sB=4v.sub.OUTLET.sup.2.

(128) B. Aortic Pressure Drops Downstream of the Aortic Valve

(129) FIG. 14 and FIG. 15 show the temporal transients of the averaged pressure drop components (total, kinetic, advective and viscous) during systole for the two groups of patients respectively obtained in the AAR and DAR defined above. Tabulated results report how the widening of the aortic flow jet downstream of the VC captured in the AAR causes a recovery of the total TPD in Group II subjects, with drop magnitudes comparable with that observed in the TVR (see FIG. 11) but with opposite sign (absolute peak means of 28.74 mmHg in AAR and 31.84 mmHg in TVR). Once again, in Group II the spatial acceleration effects due to aortic morphological changes are dominant over dissipative losses by approximately two orders of magnitude (absolute peaks mean of 0.15 mmHg), with an increased impact of the kinetic term to the total pressure drop (absolute peaks mean of 16.30 mmHg during acceleration and 14.14 mmHg during deceleration). Clear differentiation is shown between patients in Group I and II for all the pressure terms.

(130) In the DAR, the total pressure drops decrease in both groups, with the kinetic component prevalent over others and the advective component sensibly reduced, due to the absence of abrupt variations in the aortic geometry or obstacles in the blood flow hemodynamics. These values are similar to previously reported results in healthy controls.sup.22.

(131) C. The Impact of the Velocity Profile on the Pressure Drop Estimation

(132) The observed difference between WERP and Bernoulli formulations is experimentally verified with a in silico study. A steady flow on a straight tube with a change of diameter is considered. Inlet and outlet velocity fields v(x, y) were imposed analytically using the generic formula for poweroids,

(133) $\begin{array}{cc}v\left(x,y\right)=v_{\mathrm{MAX}}\left(1-\frac{{\left({\left(x-x_C\right)}^2+{\left(y-y_C\right)}^2\right)}^{\beta \text{/}2}}{R^{\beta }}\right),& \left(\mathrm{Equation}\mathrm{III}\mathrm{.1}\right)\end{array}$

(134) where v.sub.MAX is the peak velocity, x.sub.c and y.sub.c are the coordinates of the center, R is the radius and β is a coefficient accounting for the shape of the profile. We defined a reference case, by choosing the pipe dimensions and flow properties such as the cardiac output CO=5 L/min the ratio between outlet and inlet radii R.sub.OUTLET/R.sub.INLET=0.25, the density ρ=1060 kg/m.sup.3 and viscosity μ=0.004 Pa.Math.s to be representative of those in the human thoracic aorta in the presence of AS. Additionally, we selected a spatial discretization dx=0.5 mm and a velocity shape coefficient β=4 to reproduce a quasi-paraboloidal profile. We thus compared the pressure gradient ratio PGR=Δp.sub.cB/Δp.sub.cAW estimated with cB and cAW formulations, selectively testing: (1) the impact of the cardiac output (CO=4 L/min and CO=6 L/min), (2) the stenosis level in terms of the ratio between radii (R.sub.OUTLET/R.sub.INLET=0125 and R.sub.OUTLET/R.sub.INLET=0.5), (3) the spatial discretization (dx=0.25 mm and dx=1 mm) and the (4) shape of the velocity profile in terms of the shape coefficient (β=2 and β=10), in order to reproduce configurations that are likely to be found in the human aorta, spanning from paraboloidal (β=2) to blunt profiles (β=10) see FIG. 8.

(135) Results show a global overestimation obtained with the Bernoulli approach, independent of the spatial discretization, the outlet/inlet radii ratio or imposed flow rate. On the contrary, the difference between WERP and Bernoulli estimates is highly dependent on the shape of the 3D velocity profile, with minimal gap obtained with blunt profiles (ratio of PGR=1.18).

(136) D. Velocity Profiles at the Vena Contracta from PCMRI Data

(137) Simulated 3D color Doppler echocardiography velocity profiles from PCMRI data at the VC are shown in FIG. 17, with surface plots of the 2D velocity field and 12 curves for each case showing the 1D velocity profile obtained as explained previously. Velocity profiles are generally blunt, with v.sub.MAX<2.5 m/s for subjects in Group I, clearly showing lower variability of the velocity profiles when compared to subjects in Group II, where peak velocities are v.sub.MAX≥2.5 m/s. The velocity profile is highly non-axisymmetric for Group II patients, therefore motivating the larger variability obtained in this group.

2nd Embodiment: System Description

(138) The second embodiment relates to performing only a subset of the WERP processing, by considering only the advective term of the full WERP equation. As such, less information is required, and simpler imaging systems may be used, particularly 2D PC MRI systems, and 3D Doppler ECG systems. FIG. 9 is a block diagram of a 2D MRI system that may be used to perform the processing to calculate the pressure drop.

(139) Here, a 2D phase contrast magnetic resonance imaging (2D PC MRI) system 902 is provided, comprising MR imaging coils within which the subject is located, controlled by an MRI imaging control system 98, including an MRI controller processor 90. The MRI imaging control system 98 including the MRI controller processor 90 function in a conventional manner to allow 3 dimensional phase contrast magnetic resonance image data to be obtained, for example of an internal blood vessel of the subject the pressure difference across it is desired to know. For example, it may be desirable to measure transvalvular pressure drops (TPD) along the transvalvular region in the heart, between the left ventricular outflow tract and the vena contracta. Of course, other blood vessels may also be monitored, as desired.

(140) The 2D PC MRI system 98, 902, collects 2D PC MRI data 92 of the imaging subject, which is saved for further processing and analysis.

(141) The 2D PC MRI system 98, 902, also includes a complete advective WERP computation program 94, which acts to process the 2D PC MRI data 92 as described above in accordance with the complete advective WERP processing method described, to obtain WERP pressure estimation data 96. The complete advective WERP processing method is a simplification of the WERP equation by only taking one of the three additive terms of energy, the advective term. The resulting pressure estimation data 96 gives an estimate of pressure at the measured inlet and outlet planes, but neglects the effects of the vessel walls between the inlet and outlet. As such, it does not provide a complete solution for every point within the vessel between the inlet and outlet, but instead provides solutions for the inlet and outlet only, but from which pressure drop across the vessel can then be found.

(142) FIG. 10(A) illustrates the 3D Doppler ultrasound ECG equivalent system. Here a transducer 1000 obtained 3D Doppler measurements from the subject, under the control of 3D Doppler ECG system 102, comprising 3D Doppler ECG controller 104. The system produces 3D Doppler ECG data 106, and which is then used as an input to the cAW calculation program 108, which performs the necessary processing as described previously to find solutions for the inlet and outlet planes of the vessel being monitored, which is then stored and output as cAW output data 110.

(143) FIG. 7 illustrates the basic method of operation of the cAW approach. Here, at s.7.2 the MRI data or ECG data of the inlet and outlet planes of the subject vessel for which pressure drop across is to be found is obtained. Then, from this data the complete advective WERP values for the inlet and outlet planes are found, for example using the processing represented by equation A.3 above (s.7.4). Finally the complete advective WERP data is stored and/or output to a user, for example on a display screen (s.7.6).

Third Embodiment: System Description

(144) The third embodiment relates to the simplifed advective WERP approach, which can be performed using 2D Doppler ECG data, although of course it may also be applied with other imaging systems that obtain more information, such as 2D PC MRI, 3D Doppler ECG. In this example, however, FIG. 10(B) shows a 2D Doppler ultrasound ECG equivalent system. Here a transducer 1000 obtains 2D Doppler measurements from the subject, under the control of 2D Doppler ECG system 152, comprising 2D Doppler ECG controller 154. The system produces 2D Doppler ECG data 156, and which is then used as an input to the sAW calculation program 158, which performs the necessary processing as described previously to find solutions for the outlet planes of the vessel being monitored, which is then stored and output as sAW output data 160.

(145) FIG. 8 illustrates the basic method of operation of the sAW approach. Here, at s.8.2 the ECG data of the outlet plane of the subject vessel for which pressure drop across is to be found is obtained. Then, from this data the simplified advective WERP values for the inlet or outlet planes are found, for example using the processing represented by equation A.5 above (s.8.4). Finally the simplified advective WERP data is stored and/or output to a user, for example on a display screen (s.8.6).

(146) For completeness, brief further description of the remaining figures is undertaken below.

(147) FIG. 11 shows the TPD estimated from 4D PC-MRI data using the WERP formulation. Top panel: average temporal transients of the TPD during systole of the total (top left), kinetic (top right), advective (bottom left) and viscous (bottom right) components on Group I (light grey filled area) and Group II (dark grey filled area) subjects. Bottom panel: distribution (mean±standard deviation) for peak TPD Δp.sub.MAX(mmHg), including the significance of differences between Group I and Group II (unpaired T-test). Peak values of the kinetic pressure drops are reported on acceleration and deceleration (into brackets, negative peaks) systolic events separately. It must be noted that the pressure components do not add up to the total TPD, as the peak kinetic, advective and viscous drops do not happen simultaneously.

(148) FIG. 12 shows correlations between mean TPD estimated from 4D PC-MRI data with cAW against sB (A), cB (B) and sAW (C) formulations. Case-specific values and regression line for Group I (blue) and Group II (red) subjects. Grey dashed line in (A) shows the mild-to-severe stenosis threshold defined as Δp.sub.sB=20 mmHg (Baumgartner 2009) from the simplified Bernoulli formulation.

(149) FIG. 13 shows the results of computation of WERP advective drops based on 4D PC-MRI and idealized echocardiographic data. (A) Systolic advective TPD transients computed from 4D PCMRI (black solid line), color Doppler 3D (grey solid line) using the sAW formulation, color Doppler 2D using the sAW formulation to 1D profiles obtained at the peak systolic frame (grey circles), and continuous Doppler 1D using the sB formulation (grey dashed line). (B) Linear regression and correlation factor between peak advective TPD computed from 4D PC-MRI and Doppler echocardiographic data: color Doppler 3D using the sAW formulation (solid blue line), color Doppler 2D using the sAW formulation to 1D profiles 2 (solid black line with errorbars for 12 sampled profiles) and Doppler 1D using the sB formulation (solid red line). Case-specific values are shown for each Doppler based acquisition technique.

(150) FIG. 14 illustrates pressure drops in the AAR estimated from 4D PC-MRI data using the WERP formulation. Top panel: average temporal transients of the pressure drops during systole, with total (top left), kinetic (top right), advective (bottom left) and viscous (bottom right) components on Group I (light grey filled area) and Group II (dark grey filled area) subjects. Bottom panel: distribution (mean±standard deviation) for peak drops Δp.sub.MAX(mmHg) including the significance of differences between Group I and Group II (unpaired T-test). Peak values of the kinetic pressure drops are reported on acceleration and deceleration (into brackets, negative peaks) systolic events separately.

(151) FIG. 15 shows pressure drops in the DAR estimated from 4D PC-MRI data using the WERP formulation. Top panel: average temporal transients of the drops during systole, with total (top left), kinetic (top right), advective (bottom left) and viscous (bottom right) components on Group I (light grey filled area) and Group II (dark grey filled area) subjects. Bottom panel: distribution (mean±standard deviation) for peak drops Δp.sub.Max(mmHg), including the significance of differences between Group I and Group II (unpaired T-test). Peak values of the total and kinetic pressure drops are reported on acceleration and deceleration (into brackets) systolic events separately.

(152) FIG. 16 shows the results of an in silico test on a 3D straight pipe with steady velocity field. Top panel: representation of velocity profiles at inlet/outlet planes obtained with different velocity shape coefficients. Bottom panel: pressure gradient ratio (PGR=Δp.sub.cB/Δp.sub.cAW) between drops estimated with cB and cAW formulations as a function of cardiac output outlet/inlet radii ratio R.sub.OUTLET/R.sub.INLET, spatial discretization dx and velocity shape coefficient β.

(153) FIG. 17 shows 2D velocity profiles from mimicked 2D color Doppler echocardiographic images at the outlet plane for Group I (FIG. 17A) and Group II (FIG. 17B).

(154) Various modifications to the above described embodiments, whether by way of addition, deletion or substitution, will be apparent to the skilled person to provide additional embodiments, any and all of which are intended to be encompassed by the appended claims.