1 Introduction In recent years, cardiac dysfunction has become the most common cause of death in the western world. Advances in imaging science and, more recently, computational physiology provide significant potential to circumvent many of the current limitations in diagnosis and therapy planning. Key to the application of these technologies is the extraction of clinically relevant information from patient data. It is within this context that the focus of this work is on an automated computational model that computes relative pressure fields based on dynamic flow velocity information. Ongoing clinical research has established phase contrast magnetic resonance velocity mapping as a useful tool to gain non-invasive insight into dynamic cardiovascular blood flow in a wide range of contexts (Chai and Mohiaddin, 2005). Multi-directional intra-cardiac flow has been studied, among others, by Kilner et al. (2000) relative to selected planes and by Wigström et al. (1999), Hope et al. (2007), Markl et al. (2007, 2011a), and Pitcher et al. (2011) relative to volumes of acquisition. To build on the success of these relatively comprehensive non-invasive measurements of flow, we set out to derive, directly from the velocity data, intra-vascular differences of pressure. To date, fluid mechanics models have been widely used to analyse cardiovascular blood flow and more recently been integrated with tissue mechanics to understand coupled cardiac function. The fluid domain behaviour can be evaluated based on numerical discretisation techniques (McQueen and Peskin, 2000; Oertel and Krittian, 2011). The flow domain can be coupled numerically to solid mechanics models based on monolithic (Lemmon and Yoganathan, 2000; Nordsletten et al., 2008, 2010) or even mixed approaches (Krittian et al., 2010a,). In these models, both the intra-ventricular pressure and velocity fields are a direct consequence of the continuum mechanics principles of mass and momentum conservation as well as the imposed boundary conditions. An alternative to calculating flow fields from pressure boundary conditions is to determine pressure from known flow fields. In this case the so-called Pressure Poisson Equation (PPE) can be derived directly from the well-known Navier–Stokes equations as shown among others by Gresho and Sani (1987). This pressure information is in turn very valuable for the formulation of more realistic boundary conditions for the models described above. Song et al. (1994) used the PPE to determine relative pressure fields from a sequence of ultrafast cardiac CT images. Yang et al. (1996) and, more recently, Ebbers et al. (2001, 2002), Ebbers and Farnebäck (2009) applied a similar approach for the computation of flow pressure fields from MR velocity mapping. Their applied mathematical formulation is based on the assumption that the contribution of viscous terms to the pressure calculation can be neglected which holds true only for high Reynolds number flow. Furthermore, the underlying numerical discretisation requires an iterative solution in order to determine unknown boundary conditions. The need for applying these boundary conditions on the fluid domain further complicates the direct use of the actual imaging space as computational domain. Therefore, Bock et al. (2011) have shown that a more accurate pressure field estimation is achieved by extracting the actual fluid domain from its surrounding imaging space. With this work we provide a numerical model implemented in OpenCMISS 1 (Bradley et al., 2011) that allows the direct computation of cardiovascular pressure based on 4D flow data, and relative to a given reference point. The novelty of this work consists in the numerical discretisation through finite elements which turned out to enable a number of significant advantages both from the computational and methodological point of view (Hassanzadeh et al., 1994). Although the imaging space is meshed automatically based on the actual MR resolution, data projection in general allows modified element sizes and interpolation orders. By applying the functional projection of the given velocity data, we bypass the sensitive determination of boundary conditions (Gresho, 1991, 1998a,b) and establish an overall volume source field as the sole driving force of the PPE problem. This important step now allows a modular composition of the actual computational fluid domain as a subspace of the original imaging space. Degrees of freedom that lie outside the fluid domain can be eliminated based on the masking information provided. Furthermore, a mathematical isolation of separated fluid domains inside the same imaging space becomes possible. This embedded-yet-isolated flow field approach enables a quick, robust and accurate calculation of the cardiovascular pressure field not only on predefined integration paths but at any point in space where masked velocity data is present. 2 Materials and methods In order to mathematically describe the relative pressure caused by the dynamic, three-dimensional, viscous and often highly complex character of cardiovascular blood flow, we need to consider the complete formulation of the well-known Navier–Stokes equations (Oertel and Krittian, 2011). Having enforced a divergence-free state of the given velocity data, the Navier–Stokes equations can be further transformed into the so-called Pressure Poisson Equation in order to compute the corresponding relative pressure field. In the following sections we provide the mathematical and numerical foundations for the pressure estimation process, and highlight the advantages of the finite element technique for discretising the PPE problem within the cardiovascular pressure estimation process. 2.1 Governing equations The pressure estimation process is based on the continuum mechanics principles of mass and momentum conservation. The underlying equations can be used to derive the PPE foundations needed for the pressure estimation process presented. Computing the pressure distribution p that corresponds to a given incompressible flow field u , we expect u to satisfy the divergence-free condition ∇ · u = ! 0 . Following Newton’s second law, the relative pressure distribution can be seen as a consequence of transient and convective momentum, viscous resistance and volume forces. On an Arbitrary Lagrangian Eulerian (ALE) reference frame (Nordsletten et al., 2008), this condition is formulated by the Navier–Stokes equations, (1) ρ ∂ u ∂ t + ( ( u - w ) · ∇ ) u = f - ∇ p + μ △ u , where u represents the velocity, w the reference velocity, t the time, f is a volume force, p the pressure and ρ and μ the fluid density and viscosity, respectively. Obtaining a pressure distribution from its gradient given in Eq. (1) is not straightforward. Most approaches identify spatial integration paths which are often significantly sensitive when applied to noisy input data (Ebbers et al., 2002). In order to include smoothing options and to avoid boundary condition sensitivities, we have chosen to start with a higher-order pressure derivative which yields the PPE problem as the basis of our finite element approach. Rearranging Eq. (1) for ∇p yields, (2) ∇ p = b , ∇ p = b , where the right-hand side vector b is a function of given velocity data and depends on the constitutive properties of blood: (3) b = f + μ △ u - ρ ∂ u ∂ t + ( ( u - w ) · ∇ ) u . The PPE is now defined as the divergence of Eq. (2) and gives us a higher order derivative of the pressure field p: (4) △ p = ∇· b . △ p = ∇ · b . 2.2 Numerical discretisation In order to solve for cardiovascular pressure fields, especially finite difference approaches have been used in the past (Ebbers et al., 2001, 2002; Ebbers and Farnebäck, 2009; Yang et al., 1996) driven by so-called Neumann boundary conditions, which are often sensitive and difficult to determine (see Appendix A). In this section we introduce a finite element based approach to the field of cardiovascular pressure estimation, driven by volume sources rather than surface fluxes (Hassanzadeh et al., 1994). This not only avoids the use of gradient boundary conditions but also allows a reduction of the computational domain at a later stage. We transform Eq. (4) into the weak form by multiplying by a test function q and integrating over the computational domain Ω. We then look for a p ∈ H 1 (Ω) such that, (5) ∫ Ω (∇·∇ p)qd Ω = ∫ Ω (∇· b )qd Ω, ∀ q ∈ H 1(Ω) ∫ Ω ( ∇ · ∇ p ) qd Ω = ∫ Ω ( ∇ · b ) qd Ω , ∀ q ∈ H 1 ( Ω ) An advantage for the pressure estimation approach is based on applying integration by parts to both the left-hand and right-hand side of Eq. (5), which yields the surface integral ∫ Γ ((∇ p - b )· n )qd Ω = 0, ∫ Γ ( ( ∇ p - b ) · n ) qd Ω = 0 , leaving only the volume integrals (6) ∫ Ω ∇ p·∇ qd Ω = ∫ Ω b ·∇ qd Ω. ∫ Ω ∇ p · ∇ qd Ω = ∫ Ω b · ∇ qd Ω . Following a standard Galerkin finite element discretisation we get the matrix system K mn p n = s m for Eq. (6) for m, n = 1, …, N with the number of degrees of freedom (DOF) N and (7) K mn = ∫ Ω ∇ φ n · ∇ φ m d Ω = ∫ Ω ∂ φ n ∂ x k · ∂ φ m ∂ x k d Ω . The right-hand side source is discretised as: (8) s m = ∫ Ω b · ∇ φ m d Ω = ∫ Ω b k · ∂ φ m ∂ x k d Ω , with (9) b k = - ρ ∂ u k ∂ t + ( u i - w i ) ∂ u k ∂ x i ︸ acceleration terms + μ ∂ 2 u k ∂ x i ∂ x i ︸ viscous terms . The scalar function φ m represents the finite element test functions, φ n the basis functions for the resulting pressure field with m, n = 1, …, N. The vector x contains global coordinates in D-dimensional space and i, k = 1, …, D. It should be noted that Eq. (9) accounts for both viscous and inertia terms. It can therefore be applied to the whole range of laminar, low- and high-Reynolds number flows (see Appendix A). In order to capture the second derivative terms correctly, a tri-cubic Lagrangian or hermite basis function is suited best for the finite-element implementation. The latter shows further potential to improve real data quality by applying additional projection and data smoothing methods. 2.3 Embedded pressure Poisson approach For the general purpose of cardiovascular pressure estimation, we characterise the embedded velocity fields by introducing the element-based labelling factor κ into Eq. (6) which yields: (10) K ∼ mn p n = κ s m , where (11) K ∼ mn = ∫ Ω κ ( ∇ φ m · ∇ φ n ) d Ω , and s m is defined in Eq. (8). Assuming a velocity screen procedure that results in a discretised domain Ω containing both the fluid domain of interest Ω int and the surrounding area Ω ext , we can now use κ to perform the PPE computation without extra segmentation or mesh adaptation where κ = 1 on Ω int and κ = 0 on Ω ext . Elements e are treated as boundary elements if one degree of freedom (DOF) of Ω e is labelled as an outside voxel corresponding to Ω ext . Masking information may be treated as piecewise constant or, alternatively, evaluated and scaled with 0 ⩽ κ ⩽ 1. This avoids the propagation of Ω int -source signals to Ω ext and any external influence from Ω ext on Ω int . Within this work, the matrix system given by Eq. (11) is solved using PETSc for iterative Krylov sub-space linear system solvers. 2.4 Verification and validation In order to test the discretised PPE problem, we perform the following verification process using a self-adjoint analytic solution of the form: (12) △ p = - 12 π 2 L 2 p , which represents a complex spatial pressure field (13) p ( x , y , z ) = sin 2 π x L sin 2 π y L sin 2 π z L + p 0 inside a regular cube with side length L = 1 and reference pressure p 0 = 0. We shall consider the numerical solution under mesh refinement and for different orders of interpolation. The results for mesh refinements by a factor of 2 are shown in Fig. 1 where r is the error measured as L2 norm, h represents element size and logarithmus dualis is defined as ld(r) = log 2(r). For linear, quadratic and cubic Lagrange interpolation order the corresponding graph exhibits the expected error convergence rate. In addition to demonstrating the ability of our approach to solve the matrix systems accurately, we want to provide the proof of concept when embedding given flow fields in a data acquisition space. Two different examples can be seen in Figs. 2 and 3 where Ω ext (black) represents parts of Ω that do not contain any flow information, and Ω int (coloured) represents the embedded field. According to Section 2.3, Ω ext must be excluded from the actual computational domain. The left-hand side of Fig. 2 shows the velocity magnitude fields in Ω int of two fully developed channel flows surrounded by Ω ext . Arrows indicate the respective flow direction resulting in one positive and one negative pressure gradient. The analytic values for velocity and pressure distribution are given by the following formulas: (14) u ( y ) = ( D 2 - 4 y 2 ) u max D 2 , p ( x ) = p 0 - 8 μ x u max D 2 , respectively, where u max is the maximum velocity inside a channel with width D, flow direction x, centre line y = 0 and reference pressure p 0. Due to the elimination of DOFs outside the fluid domain, the two respective pipe flows (u 1 = u(y) = −u 2; p 1 = p(x) = −p 2) are completely separated and isolated from each other. Whereas Fig. 2 represents a friction-driven test case, Fig. 3 shows the results of a third, inertia-driven test problem where velocity magnitude fields are generated by four pairs of cylinders (inner radius R i , outer radius R o ) rotating at different speeds (inner velocity ω i , outer velocity ω o ). The analytic values for velocity and pressure distribution are given by the following formula (15) u ( r ) = C r , p ( r ) = p 0 - C 2 2 r 2 , with R = R i < r < R o and C = ( 9 + c / 10 ) R 2 · 1 s . The left-hand side of Fig. 3 shows this linearly increasing velocity magnitude from case c = 1 to case c = 4; the right-hand side shows the expected increase of pressure on the outer wall of the cylinder system due to the increasing centrifugal forces. The calculated pressure level for all four cases is in very good agreement with the analytical values given by Eq. (15). As in the previous example, all four velocity fields are mathematically isolated as needed for the following cardiovascular application. Parameters and R for Figs. 2 and 3 have been chosen such that the vertical geometric extension of each flow domain is 40% of the overall domain. With L = 1, this means D = 0.4 and R = 0.2 in Figs. 2 and 3, respectively. Initial and constitutive parameters have been set to 1 (mu and u m ax in Fig. 2), reference pressures to mid range (p 0 = 25 in Fig. 2 and p 0 = 0 in Fig. 3) to allow for a pressure scale starting from 0. Fig. 1 shows the effect of finite-element order and discretisation density on the convergence test example demonstrating the expected decreasing error with an increasing number of elements. The test example in Fig. 2 illustrates how our FEM formulation solves a simple Pouiseulle flow, what could not be done by traditional finite difference methods (see Appendix A). The test example in Fig. 3 verifies the correct solution of the equation as compared to previous methods Song et al. (1994). 3 Modelling and application In the following section, we apply the data processing approach described above in order to estimate the velocity-based pressure field for datasets of four healthy human subjects. Fig. 4 explains the following workflow from data acquisition to pressure evaluation. 3.1 Data acquisition and processing The velocity input for this study is provided by phase-contrast MR imaging, a technique that allows blood flow velocity to be measured and post-processed non-invasively. Measurements were performed on a 3 T system (Magnetom Trio, Siemens AG, Erlangen, Germany) with a standard 8-channel phased-array coil. 4D flow data with three-directional velocity encoding and covering the whole heart fluid domain were acquired using a navigator respiration controlled and ECG-gated rf-spoiled gradient echo sequence (Markl et al., 2007) (spatial resolution: 2.95 × 2.50 × 2.90 mm3, temporal resolution: 38.4 ms, velocity encoding: 150 cm/s, time frames per cardiac cycle: 17). Initial raw data normally contains magnitude and three-dimensional phase information for each voxel of the initial imaging space. Voxel-based phase shifts can be directly transformed into velocity vectors which marks the starting point for our cardiovascular pressure estimation. Data-processing was established to further enhance quality (e.g. eddy-current elimination, velocity aliasing or noise filtering) and to allow for fluid domain representation (i.e. MR segmentation and flow field masking). Although the pressure estimation approach presented in this work can be independently applied to any spatially distributed velocity field, we follow the data-processing steps defined by Velomap (Bock et al., 2007), a Matlab (The MathWorks, Inc.) based pre-processing tool developed by the Diagnostic Radiology Department, University of Freiburg, Germany. Velomap has already proven to supply reliable flow fields in previous in vivo flow studies (Markl et al., 2011a,b; Stalder et al., 2008). Noise masking was performed by thresholding of the signal deviation of the magnitude data in order to exclude regions with low signal intensity. Further noise reduction and separation of static tissue and vessels is achieved by comparing the standard deviation of the velocity–time course for each pixel in the flow data set (Bock et al., 2007). MR segmentation output was used to mask geometric entities based on an averaged fluid domain representation (Bock et al., 2010). The following procedures are typically applied before the 4D flow data enters our pressure estimation workflow: Anti aliasing, Noise masking, Eddy current reduction. Based on a ’speed sum squares’ iso-surface representation over all time-steps T (16) σ = ∑ j = 1 T ∑ i = 1 3 u i 2 ( t j ) = const . , an averaged segmentation of the cardiovascular geometry was created. Since this iso-value represents the basis for fluid domain masking (Ω int /Ω ext ), careful distinction of adjacent cavities or vessels was taken into account. 3.2 Pressure estimation workflow In order to allow a smooth and straightforward representation of the cardiovascular geometry of interest, our approach follows the mean fluid domain approach often used in 4D flow analysis based on MR segmentation information. 3.2.1 Geometrical representation The left-hand side of Fig. 5 shows an isosurface (highlighted in red) of σ, defined in Eq. (16), clearly indicating the left and right ventricle, the aorta and adjacent large vessels. Valid velocity information is available in the entire imaging space Ω (square box) following an Eulerian description approach. However, we will focus on the internal volume Ω int – with σ = 0.2 m2/s2 set as lower bound – assuming Ω int is a sub-domain of the true fluid domain over the whole cardiac cycle. This value has been chosen according to Bock et al. (2011) in order to separate the aortic arch allowing for the best mean representation of the cardiovascular velocity field, and must be optimised for the respective cases under consideration. Following this mean geometrical representation, the right-hand side of Fig. 5 now visualises the fluid domain masking and the subsequent separation of the whole imaging domain Ω into Ω int und Ω ext . This information will be used in a similar way as for the cases in Figs. 2 and 3; they are also the basis for the numerical approach and its parameter settings for κ in Eq. (10). 3.2.2 Pressure field estimation In order to illustrate clinically relevant results, computation of the pressure differences in the rest of the article is performed in the isolated aorta, which is manually segmented from the automatic identification of the fluid domain represented in Fig. 5. This isolation removes the contamination from the rest of vascular structures (mainly the pulmonary arteries). Fig. 6 shows the discretised imaging domain Ω, where each voxel of the image is a control point of a tri-cubic Lagrangian hexahedra. The masked DOFs corresponding to the aorta are used to distinguish between the fluid domain of interest Ω int and Ω ext . Whereas nodes which belong to Ω int carry velocity information, all nodes of Ω ext get eliminated from the initial computational domain, effectively decreasing the system size and, hence, increasing computational efficiency. 3.3 Visualisation and post-processing 3.3.1 Velocity field and estimated relative pressure The left-hand side of Fig. 7 shows 6 different snapshots at different points in time during systole chosen for best visualisation results. Starting with almost steady-state conditions (snapshot 1), one can clearly identify the increasing blood velocity magnitude in early systole (snapshots 2 and 3) and the blood momentum transported further downstream the aorta (snapshots 4 and 5) until the blood velocity magnitude decreases to a minimum (snapshot 6). The right-hand side shows the corresponding relative pressure fields directly computed from the respective time-frames. Since there is almost no blood flow present at the very beginning of systole (snapshot 1) no differences of relative pressure can be identified. However, as soon as the early systole begins, a pressure drop over the aortic valve plane can be seen (snapshot 2) followed by an increase of both the magnitude and differential values of the relative pressure field (snapshot 3). After the main blood flow has passed the aortic arch, the highest relative pressure value develops due to the centrifugal forces (snapshot 4). Finally, the relative pressure field returns to its initial state (snapshot 5 and 6). The reference point for the relative pressure field has been located at the end of the descending aorta. Pressure values are measured relative to the pressure at this reference point. Consequently, the relative pressure at this point must be equal to zero as shown in Figs. 7–9, respectively. 3.3.2 Spatial and temporal correlation In order to understand the velocity/pressure interdependence, we plot in Fig. 8 velocity and pressure magnitudes for three different time-frames during the systolic phase. Two different cutting planes (one through the aorta and one through the ventricle) are used to project and visualise both velocity magnitude and pressure. Relative to the reference point at the end of the descending aorta, one can clearly identify a pressure drop during early systole followed by the highest pressure in the aortic arch as a consequence of the stagnation pressure and the change of direction of the blood column. Whereas Fig. 8 gives a first indication of the spatial relation between velocity and pressure, Fig. 9 allows a temporal analysis of the cardiac cycle. Point 1 is placed inside the left ventricle, followed downstream by points 2 and 3 until the end of the descending aorta is reached at point 4; close to the pressure reference point. It is interesting to note that the early systole causes positive and negative relative pressure changes at points 1 and 2, respectively, before the pressure level returns to its force-free state. The effect of the directional change due to the aortic arch causes an earlier pressure increase accompanied by a higher magnitude. As expected, control point 4 shows no significant pressure change, being located at the end of the descending aorta. Fig. 10 illustrates the acceleration of blood at early systole (positive values), followed by the deceleration during late systole (negative values) and a relatively flat profile during diastole. 4 Discussion In this study we have presented a methodology based on the Pressure-Poisson equation to enable direct access to relative pressure fields derived from 4D flow measurements. Whereas 4D flow information are increasingly available from non-invasive imaging to date, there has been no computationally efficient approach established to determine the corresponding pressure field due to both inertial and viscous effects. The proposed framework provides a platform for a wide range of applications, both in clinically relevant diagnosis and in computational cardiac modelling. Our approach exhibits several novel features that overcome many of the deficiencies associated with existing methods based on alternative approaches. Firstly, conventional methods are typically driven by Neumann conditions in the form of pressure gradients defined on the relevant boundaries of the region on interest. However, such boundary data is generally not readily available. To address this issue conventional methods resort to costly iterative schemes to computationally approximate this boundary data, but these iterative schemes often lack robustness and may converge slowly or not at all. In contrast, the approach outlined in this paper removes the need for such boundary conditions and iterative schemes. Using the finite-element formulation enables integration by parts of the source term in the Pressure-Poisson equation, effectively eliminating unknown boundary integrals, as previously demonstrated by Hassanzadeh et al. (1994). It is this source term which drives the underlying pressure estimation process. Secondly, the method eliminates the no-flow domain associated with a zero source term, which, in combination with eliminating the need to determine a solution via iteration, reduces computational cost relative to current schemes. A final benefit is that our approach, by accounting for pressure changes due to both acceleration and viscous resistance (see Appendix A), is valid for both low- and high-Reynolds number laminar flows. While the magnitude of these viscous terms was relatively small for the healthy volunteer data, for pathological cases these terms are likely to be significantly more important. For these cases taking into account simple turbulent assumptions may also be relevant, however, this extension was beyond the scope of the current study. To verify the correctness and demonstrate the potential of this technique, we have presented various numerical results in Sections 2 and 3. A convergence study has shown that the proposed method yields the expected rate of convergence under mesh refinement. Furthermore, the Poiseuille-flow model in two channels reproduced the correct linear pressure profile. This result is in contrast to results determined using methods employing iterative determination of boundary conditions which typically fail for this problem due to the absence of a source term from which to derive the boundary conditions. For the third test problem of four cylinders rotating at different speeds, our computed solution compares well with the analytic solution as well as with results published in literature. While the above results serve to verify the method, in Section 3 we also assessed the potential of our method on 4D flow measurements acquired on a human aorta. We presented relative pressures on the surface of the aorta as well as on cutting planes, and variations of relative pressures over a heart cycle at specific points of interest in the aortic branch. These results are in good agreement with results published in Yang et al. (1996), Ebbers and Farnebäck (2009) and Bock et al. (2011). This demonstrates the viability and potential of the proposed methodology for the non-invasive estimation of relative pressures to complement 4D flow velocity data and assist in the assessment of healthy and diseased conditions. While these results illustrate the potential of the FEM system to solve the PPE problem, there remains further work to understand the behaviour of pressure estimation in vivo. The influence of MR image acquisition factors (limited spatial and temporal resolution, average through different heart cycles, Eddy currents and patient movement among others) on physical correctness is of particular interest in terms of interpretation of cardiovascular Reynolds numbers. The uncertainty is less influenced by physical blood flow properties, but more by the uncertainty of velocity magnitude and its gradients. While traditional FEM techniques for measuring accuracy rely on supposed smoothness, understanding the influence of interpolation and mesh resolution on the computed vector b , and consequently the pressure, is unknown. As far as geometrical limitations are concerned, a mean representation of the fluid domain is taken into account. However, domain discontinuity like through a closed cardiac valve must be treated carefully, as the pressure field determined is relative to some reference pressure that is unique within a connected fluid domain, but generally distinct for separated regions. This means that changes of the pressure reference value through discontinuities, such as during the valve opening and closing process, require more sophisticated measurement input. Future work will address data quality and enhanced cleaning of data with noise. Stipulating that the flow is governed by the incompressible Navier–Stokes equations, the velocity field should theoretically be divergence-free. Consequently, the degree to which the measured velocity field violates the divergence-free condition provides a useful quality measure of the acquired data. Accordingly, to improve the quality of the data, the error components corresponding to a violation of the divergence-free condition can be eliminated by projecting the velocity data onto a divergence-free finite-element approximation before using it as input to the pressure estimation. Further validation will then be carried out by comparing our pressure estimation results to data from combined 4D flow and pressure catheter measurements.