CONVERGENCE ANALYSIS OF PRESSURE RECONSTRUCTION METHODS FROM DISCRETE VELOCITIES

. Magnetic resonance imaging allows the measurement of the three-dimensional velocity field in blood flows. Therefore, several methods have been proposed to reconstruct the pressure field from such measurements using the incompressible Navier–Stokes equations, thereby avoiding the use of invasive technologies. However, those measurements are obtained at limited spatial resolution given by the voxel sizes in the image. In this paper, we propose a strategy for the convergence analysis of state-of-the-art pressure reconstruction methods. The methods analyzed are the so-called Pressure Poisson Estimator (PPE) and Stokes Estimator (STE). In both methods, the right-hand side corresponds to the terms that involving the field velocity in the Navier–Stokes equations, with a piecewise linear interpolation of the exact velocity. In the theoretical error analysis, we show that many terms of different order of convergence appear. These are certainly dominated by the lowest-order term, which in most cases stems from the interpolation of the velocity field. However, the numerical results in academic examples indicate that only the PPE may profit of increasing the polynomial order, and that the STE presents a higher accuracy than the PPE, but the interpolation order of the velocity field always prevails. Furthermore, we compare the pressure estimation methods on real MRI data, assessing the impact of different spatial resolutions and polynomial degrees on each method. Here, the results are consistent with the academic test cases in terms of sensitivity to polynomial order as well as the STE showing to be potentially more accurate when compared to reference pressure measurements.


Introduction
The intra-arterial spatial distribution of the blood pressure can be measured by means of catheterization [2].This technique consists in inserting a catheter equipped with a pressure transducer into the vasculature of the patient and manoeuvring it, under local anaesthesia and guided by fluoroscopy, to the location of interest.Although it is the "clinical gold standard" for pressure quantification, the invasive nature of the method is Figure 1: Pressure map estimation from in an experimental phantom of a thoraxic aorta [22], adapted and reprinted from [16]: Left: The 4D Flow MRI velocity measurements.Right: Relative pressure map computed from velocity data.Shown: cuts through roughly the center of the vessels.associated with a risk of complications [17,25,27].Given that there are recommendations to avoid its use [24], to compute the pressure difference from measured flow fields is strongly preferred.
Time-resolved 3D velocity encoded magnetic resonance imaging, or 4D flow MRI, offers measuring the complete 3D velocity field within a region of interest [14,20], Additionally, 4D flow MRI allows the computation of several hemodynamic parameters, which can be used as new biomarkers [20].Among those hemodynamic outputs, relative pressures are one of the most popular ones due to its potential to replace the invasive catheterisation procedures.
When the full velocity field is measured as in 4D flow MRI -naturally at a finite spatio-temporal resolutionit can be inserted in the linear momentum balance of the incompressible Navier-Stokes equations (NSE).Then, the velocity terms laid in the right-hand-side while the pressure holds as an unknown, which needs to be found using appropriate discretization approaches.An example of pressure map estimation from real 4D Flow MRI data is shown in Figure 1.
In practice, and as it can be appreciated in Figure 1-Left, those measurements are obtained at limited spatial resolution -given by the voxel size in the image -and therefore the velocity entering to the right-hand-side corresponds to an interpolated version of the exact velocity.Therefore, there is not a unique numerical approach to compute the reconstructed pressures.A review and preliminary numerical comparison of methods can be found in [3].Among those methods, only a few can compute pressure fields and not just averaged pressure differences between two locations.
The first one is the so-called Pressure Poisson Estimator (PPE) [8,13] and it consists of applying the divergence to the NSE obtaining a pressure Poisson equation, similarly as it is used in projection methods [11].However, the original PPE method cannot include the viscous contribution to the pressure gradient at the level of accuracy of the measured data.Therefore, recently in [18] the PPE method was modified by adding a boundary term with the viscous contribution.
Another more modern method corresponds to the Stokes Estimator (STE) was reported in [23].The STE consists in adding to the NSE the Laplacian of an artificial incompressible velocity field with null trace leading to a linear Stokes problem for both pressure and artificial velocity fields.Such artificial velocity is supposed to be zero for perfect velocity measurements.The STE has shown more accurate results than the PPE in numerically simulated data [3,23] and in real phantom and patient data [16].However, the STE method is considerably more expensive computationally than the PPE.
To the best of the authors' knowledge, neither a mathematical convergence analysis of both PPE and STE methods or a comparison among discretization schemes for each of the methods has been reported.
Therefore, the purpose of this work is to propose a strategy for performing a priori error analysis and applied it to the PPE and STE methods.The strategy is based on the splitting of the solution in two components and adding their contributions to the overall error.Moreover, for both methods we studied different discretization strategies in order to verify the theoretical analysis and give insights on the cost-effectiveness of each approach.In order to assess the impact of discretizations on each of the methods, calculations of pressure fields based on experimental MRI data are also included.
The remainder of this work is organized as follows.In Section 2 we present and analyze the PPE method in the standard and modified variants using Continuos Galerkin approaches.Section 3 introduces the STE and analyzes the classical Taylor-Hood and a tailored PSPG discretization.Then, in Section 4 we show numerical results using three known analytical solutions for the NSE, confirming the a priori error analysis.In Section 5 the methods are assessed under different spatial resolutions and polynomial degree on experimental MRI data.Finally, in Section 6 we draw some interpretation of the results and give recommendations for the use of these methods.

The Poisson Pressure Estimator
Given measurements of the velocity , with  being free divergence [18], the term   has the form   := ( • ∇) − ∆, and  corresponds to the outward normal vector of Ω.
We will make use of the function spaces Assuming  ∈   , the weak formulation of the ( 1) is given by: find  ∈  such that where (, ) = (∇, ∇) Ω and where   is the Kronecker delta.We refer to Standard-PPE if  = 1 and Modified-PPE if  = 2 [18].The uniqueness of the solution of Problem (2) follows from the Lax-Milgram Lemma [9].Indeed, the coercivity of the left-side is straightforward.The continuity of    follows from In spite of the last term of    , for  = 2 is given as in (3), in the practice we will use the identity

Continuous Galerkin discretizations
Let { ℎ } ℎ>0 be a shape-regular family of partitions of the polygonal domain Ω, composed by elements  of diameter ℎ  with mesh size ℎ = max ∈ ℎ ℎ  .ℰ ℎ denotes the set of all edges of  ℎ and  the edges that compose it.In addition,   () denotes the polynomial function spaces defined on  of total degree less than or equal to .
The finite element spaces for the pressure approximation and velocity interpolation are: We will also consider the interpolation operators  ℎ :  ∩  +1 (Ω) →  ℎ and ℒ ℎ : where  ℎ and ℒ ℎ corresponds to a modified Lagrange interpolator with average zero and the classical Lagrange interpolator, respectively.Note that it is possible to prove that  ℎ satisfy the same error approximation properties than the classical Lagrange interpolator.Thus, the Galerkin scheme associated with the continuous variational formulation (2) reads as follows: Find with According to discrete Lax-Milgram Theorem, problem (6) has a unique solution  ℎ ∈  ℎ .
Remark 2.1.Note that from the definitions (3) and (7) we can assure that the problem (6) is not a Galerkin scheme of the continuous problem (2).Indeed, the scheme is not consistent.
The strategy to prove convergence is to use the known Strang's lemma for conformal and non-consistent cases.
In order to prove the convergence theorems, let us state the next result.
Lemma 2.2.Let us assume that  ∈  2 (Ω).Then, there exists  independent of ℎ such that Proof.In this proof we will assume that  1 ,  2 are the error interpolation constants, C is an injection constant and   an inverse inequality constant.
Using properties of interpolation given by (5) and Young inequality we obtain Then, there exists  2 and  3 independent of ℎ such that Proof. sup For the first term in the above inequality, we use Lemma 2.2 and for the third term we have Now, thanks to [7, Lemma 1.46], we have where C1 = (︁ )︁ 1/2 (see [26,Theorem 3]).Besides, from [21, Lemma 10.8] we get and then, from (10) to (12) we arrive to which allows us to arrive at the desired result, wherefrom Lemma 2.2, Finally, the next theorem holds.
In addition, we assume that ]  , for  = 1 and  = 2 respectively.Then, there exists constants  1 ,  2 and  3 such that Proof.Thanks to the Strang's Lemma (see [9,Lemma 2.27]) we have that The bound for the first term on the right-hand side follows directly from Lemma 2.3.For the second term, we will consider the interpolation operator, and then inf Theorem 2.5.Let the hypothesis of Theorem 2.4 hold with Ω is a convex polygonal domain.Then, with   is the Poincaré inequality constant and  ≥ 1.
Proof.The proof starts taking qℎ ∈  ℎ such that satisfies the equation where  is the continuous vector function representing the measured velocity field.It follows from the triangle inequality that Given that  − qℎ ∈  ∩  +1 (Ω), there exists a unique  ∈  2 (Ω) such that (see [4]) In addition, thanks to the convexity of Ω, by elliptic regularity we have there exists  Ω > 0 such that Notice, the weak formulation for ( 15) is given by Now, replacing  by  − qℎ , using the orthogonality property of ( 17), Cauchy-Schwarz inequality and interpolation properties given in (7), we get Applying the inequality ( 16) we arrive to For the second term of the right-hand side in (14) we proceed as follows: From Lemma 2.3 and Poincaré inequality, we get Hence, the result is a direct consequence of the estimates ( 18) and ( 19).

The continuous problem
The STE consists then in adding the Laplacian of an incompressible auxiliary velocity  ∈ V and to the left-hand side of the Navier-Stokes equations, i.e.: Let us define the space  =  2 0 (Ω).Hence, we can define the weak problem of (20) as: where For the analysis, we will use the following norm and if  is an linear functional operator we use the norm Note that the problem ( 21) is well posed thanks to the V-ellipticity of the bilinear form (∇, ∇) Ω , and (, ∇ • ) Ω satisfy an inf-sup condition (cf.[9, Prop.2.36]).Lemma 3.1.There exists a positive constant  ℬ such that ‖ℬ‖ ≤  ℬ .
Proof.Using triangular inequality, Cauchy-Schwarz inequality, together to the inequalities 1 ≤ √  and and the result is obtained straightforwardly, with  ℬ = 2 √ .

Discrete spaces
Now, for each ℎ let W ℎ and  ℎ be finite-dimensional spaces such that: For our error analysis, we will need to make use of some known results.

Taylor-Hood discretization
For the discrete STE we set V ℎ = W ℎ ∩ V where inf-sup stable pairs of finite elements require the use of different spaces for velocity and pressure and for this reason, we take Taylor-Hood, where  =  + 1.Otherwise, it is not possible to use conforming spaces of the lowest order for the discrete velocity.Furthermore, we will consider the property of the interpolation operator Thereby, the discrete version of the problem (21) reads as follows: where the bilinear form ℬ is like in the continuous case, and with ℒ ℎ : [ 2 (Ω)]  −→  ℎ a Lagrange interpolant.
For the first term of the right-hand side, we use Lemma 2.2.
For the second term of the right-hand-side, we have from the interpolation bounds: Finally, using the above inequalities and Poincaré inequality we get where   is the Poincaré constant.Thereby, we arrive straight at the result of the lemma.
Finally, we can derive the first main convergence result.

Stabilized PSPG discretization
Let us consider again the Stokes problem given as in (20) and its respective variational formulation (21).We will now analyze the PSPG Stabilization [12] with the end of comparing the error of convergence between the pressure obtained with both schemes.
We want to use spaces of finite element of order  for the velocity and the pressure, i.e.,  = , by means of the following stabilized formulation. where and   ℎ (•, •) defined as in (27), and Remark 3.8.Note that the term ∆ ℎ is not included in the stabilization.This is possible to do while keeping strong consistency since  = 0. Our choice allows also us to avoid conditional well-posedness of the discrete solution as in standard PSPG stabilized formulations.
Proof.Using the inequalities ( 24) and (32), Theorem 3.1 and Cauchy-Schwarz inequality we obtain that and then the result follows.
In the next lemmas we will consider the pair ( wℎ , qℎ ) ∈ V ℎ ×  ℎ which are solution of the equation where We highlight that the solvability of the problem (34) has been guaranteed in [12].
Proof.We note that (, ) and ( wℎ , qℎ ) satisfy the orthogonality property Indeed, thanks to the consistency of bilinear form ℬ we get By the triangle inequality, we can get, For the second term of the right-hand side, we must consider the result earned in [12] from where we get and so, from this inequality and (35) we obtain and thereby we arrive to Lemma 3.12.Let ( wℎ , qℎ ) and ( ℎ ,  ℎ ) be solutions of (34) and (29), respectively.Additionally, we assume that  ∈ [ 2 (Ω)]  .Then, the following bound is satisfied: Proof.Let   ℎ := wℎ −  ℎ and   ℎ := qℎ −  ℎ .Then, thanks to the stability of ℬ  given in (31) we have We take the term within the sum, making use of the Cauchy-Schwarz inequality and proceeding similarly as in (9), we obtain As a main result of this section, by employing the approximation properties and a priori estimates, we obtain the next result.Theorem 3.13 (Main Result III).Assume that the hypothesis of Theorem 3.7 holds.Then, Proof.The proof follows from combining the results of Lemmas 3.6, 3.11 and 3.12.

Numerical results for the convergence analysis
In this section, we present some numerical examples to illustrate the theoretical results previously described.The legends in the plots follow the notation: In addition, for Modified-PPE and Standard-PPE we will use the legend PPEvisc and PPE respectively.For the STE computed using Taylor Hood spaces and PSPS we will use the legend STE (TH) and STE (PSPG) respectively.
Every numerical routine has been sorted out using the open-source finite element libraries FEniCS [1].
Example 4.1.For the first example, we consider the exact solution of the two dimensional Kovasznay flow where Ω = (︀ − 1 2 , 3

2
)︀ × (0, 2) and the parameter  is given by For this illustration we have taken the Reynold number as in [19] which is given by Re = 1  .
The convergence results for Example 4.1 are shown in Figure 2 and examples of pressure and velocity fields in Figures 3-6.
First, it can be appreciated the lack of convergence of the PPE, while adding the viscous terms recovers it.Also, the STE appears to be more accurate than the PPE (visc) and it seems not to profit from the increase of polynomial order.Moreover, the STE-PSPG appears to deliver more accurate results than the STE-TH.Finally, it is worth saying that the sensitivity of all methods with respect to the polynomial order decreases when increasing the Reynolds number.
Example 4.2.Next we turn to the testing the scheme, where the computational domain is the rectangle Ω = (0, 1) 2 and we consider the exact solution of the Navier-Stokes equation given by The convergence results for Example 4.2 are shown in Figure 7 and the examples of pressure and velocity fields in  Here, the same remarks given about the results in Example 4.1 apply, except that for higher Reynolds numbers the STE methods appear to keep the sensitivity (though worsening) when increasing the polynomial order.

Computations using experimental MRI data
Experimental MRI data was used to assess the impact of discretization in the pressure estimation methods in realistic data and flow regimes.The setup consisted of a 3D printed, MR compatible phantom of the thoracic aorta with 60% of obstruction in order to produce a typical obstruction.Blood mimicking fluid was pumped into the phantom obtaining physiological velocities.The phantom was equipped with a catheterization unit to measure invasively and simultaneously the pressure gradient across the obstruction.4D Flow MRI was acquired with an isotropic voxel size of 0.9 mm and 25 time instants along the emulated cardiac cycle.We refer to [15,16,22] for the technical details of the experiment.The 4D Flow data is shown in Figure 1.
Two tetrahedral meshes for the pressure computations were constructed.The first one was created using the original 0.9 mm resolution where the nodes of the mesh correspond to the voxels center.The second mesh has 2 mm resolution created using linear interpolation on the first mesh.
Pressure maps were computed from all 4D flow data sets with the PPE, PPEvisc and STE methods.Due to the pulsatile nature of the experiment, the term −(  , ∇) Ω and − (  , ) Ω with   the backward finite difference operator between two measured time instants, were added to the righthand-side of the PPE and STE methods, respectively.This implies that the convergence analysis does not fully apply to the this experimental setup, however, the goal is merely to give an idea on how the discretization setup and methods compare when using real 4D flow MRI data.
The pressure differences, to be compared with the corresponding catheter values were defined as differences of the pressure averages over two spheres with a radius of 4 mm at locations proximally (ascending aorta) and distally to the obstruction.
For the PPE and PPEvisc continuous Galerkin finite elements with  = 1, 2, 3 were considered.For the STE, both Taylor-Hood (TH) and PSPG cases were computed, the latter with stabilization parameter  = 0.01 as in the previous section for the convergence analysis.In the 2 mm element size mesh,  = 1, 2 was tested for both TH and PSPG.In the 0.9 mm element size mesh, only  = 1 was used for TH (due to the very           high computational cost of higher-order) and  = 1, 2 was used for PSPG.These methods were implemented using the FEM library FEniCS [1].
Figure 12 shows the results of the pressure estimation, where the catheter pressure values show that the 4D flow-based pressure estimation delivers reasonable values.However, note that the catheter measurements cannot be considered as ground truth, since the precision of the pressure measurements can be considered within a few mmHg [5,6,16].It can be noted that: -PPE and PPEvisc deliver visually the same results, which may occur due to the fact that viscous effects are negligible in this type of (patho-)physiological flows.-PPE and PPEvisc allow for a larger pressure gradient when increasing  in the coarse mesh.
-PPE and PPEvisc are less sensitive to  for the finest mesh.
-STE methods allow to recover larger pressure differences than PPE methods.
-STE-PSPG delivers equal or better results than STE-TH for  = 1, what is consistent with the convergence results of the numerical tests in the previous section.-STE-PSPG seems not to profit from increasing polynomial order, what is consistent with convergence results for high Reynolds numbers.-If one may take the catheter measurements as a ground truth, STE-PSPG would deliver the most accurate results.

Conclusions
In this article, we have analyzed theoretically and numerically some strategies employed to recover pressure fields from discrete velocities using the incompressible Navier-Stokes equations.
Two main methods were analyzed, the STE and PPE.While the STE is implemented using the classical Taylor-Hood finite element spaces and the Pressure-Stabilizing Petrov-Galerkin (PSPG), the PPE is imple-mented with the traditional Continuous Galerkin Method.For the PPE, two versions have been studied, the standard one without the viscous term and a modified one that includes it.
The error analysis shows that all methods, except the standard PPE, converge to the exact solution when decreasing the element size of the image mesh ℎ.With respect to the convergence rate, terms of several orders appear in the error analysis.The numerical results determine that the PPEvisc linear order dominates in the test cases presented.For the STE, the convergence in the numerical examples varies depending on the test case and the polynomial order.
Numerical results in academic test cases show that as the Reynolds number increases, the results seem to lose sensitivity to an increase in polynomial order, in particular for the STE, while the PPEvisc shows some improvements.In many of the cases, the error also appears to decrease faster with ℎ for the STE than for the PPEvisc.Among both STE discretizations, the PSPG appears to be equally or sometimes more accurate than Taylor-Hood approximations.
The computations with real MRI data are consistent with these observations.Therefore, it seems that STE-PSPG is likely to be the method of choice then it comes to the highest accuracy and a reasonable computational cost.

2. 1 .
The continuous problemThe Poisson Pressure Estimator (PPE) consists in obtaining the pressure from the classical Navier-Stokes equation by mean a Poisson equation.That is, by applying the divergence operator on the Navier-Stokes equations one obtains ⎧

Figure 9 :
Figure 9:  1 -interpolated reference velocity and pressure fields (top) and reconstructed pressure fields with order  = 1 (bottom) for  =  in Example 4.2.

Figure 12 :
Figure 12: Pressure difference in the obstruction over time error computed from 4D flow and compared with the catheter values.