Assessment of 4D flow MRI's quality by verifying its Navier–Stokes compatibility

Abstract 4D Flow Magnetic Resonance Imaging (MRI) is the state‐of‐the‐art technique to comprehensively measure the complex spatio‐temporal and multidirectional patterns of blood flow. However, it is subject to artifacts such as noise and aliasing, which due to the 3D and dynamic structure is difficult to detect in clinical practice. In this work, a new mathematical and computational model to determine the quality of 4D Flow MRI is presented. The model is derived by assuming the true velocity satisfies the incompressible Navier–Stokes equations and that can be decomposed by the measurements u→meas plus an extra field w→. Therefore, a non‐linear problem with w→ as unknown arises, which serves as a measure of data quality. A stabilized finite element formulation tailored to this problem is proposed and analyzed. Then, extensive numerical examples—using synthetic 4D Flow MRI data as well as real measurements on experimental phantom and subjects—illustrate the ability to use w→ for assessing the quality of 4D Flow MRI measurements over space and time.


| INTRODUCTION
Time-resolved 3D flow magnetic resonance imaging, known as 4D Flow MRI, has shown in the last years potential in assessing cardiovascular diseases since it offers full coverage of the region of interest, therefore allowing for its analysis after the scan. 1,2 Additionally, it allows the computation of several hemodynamic parameters, which can be used as new biomarkers. 2 However, high-quality 4D Flow in subjects involves long time scans (>20 min) even with coarse spatio-temporal resolutions making it challenging for everyday clinical use. In order to accelerate the acquisition time, several strategies have been proposed, such as parallel imaging 3,4 which accelerates the acquisition by exploiting the sensitivity of multiple receivers, and k-space undersampling [5][6][7][8] which exploits data redundancies in frequency and time. This scan time reduction comes at the price of reducing the signal-to-noise ratio (SNR). SNR also decreases when reducing the image's voxel size. Moreover, the velocity field can only be obtained under a certain predefined range which depends on the magnetic gradient setup, therefore being potentially subject of velocity aliasing. Other artifacts may also appear due to subject's respiration and motion during the scan.
To the best of the author's knowledge, quality control of 4D Flow in clinics is based on calculation of peak/mean flows, mean velocities, flow patterns, and stroke volumes. 9,10,2 A more systematic approach is to compute the divergence field of the data: assuming the blood flow is incompressible, jumps in the divergence field may indicate the presence of artifacts. Indeed, the incompressibility assumption has been used for denoising 11,12,9,13 and as a regularization term during the reconstruction process. 14,15 However, an important limitation of the divergence is that any measured velocity will have "infinitely" large divergence compared to its reference value which would be zero. Therefore, to the best of our understanding, only the spatial distribution of the divergence may be used as an error indicator but not in absolute terms for image quality check. Moreover, it will be shown later in the article, measurement artifacts may lead to no important changes in the divergence field. Therefore, this work introduces an alternative quantitative approach for assessing 4D Flow quality by verifying the compatibility with the linear momentum conservation part of the Navier-Stokes, which when being written appropriately, includes also angular momentum and mass conservation.
The rest of this article is structured as follows. In Section 2, the mathematical model will be introduced, and a numerical method will be developed and analyzed. In Section 3, a set of relevant examples, using synthetic data, will be detailed. Numerical computations of the new model for several types of artifacts and its comparison against the divergence of the data are also shown. Results with phantom and subject's 4D Flow MRI are shown in Section 4. Finally, in Section 5, we discuss potential applications of this metric in the context of reconstruction and processing 4D Flow MRI.

| The continuous problem
We assume a physical velocity field u ! , which satisfies the conservation of linear momentum of the incompressible Navier-Stokes equations in the vessel lumen Ω: where p is the fluid pressure field, and ρ and μ are the density and dynamic viscosity of the fluid, respectively. We recall that both the convective representation of the advection term and Laplacian representation of the viscous term are written in a simplified form using that r Á u ! ¼ 0. Moreover, the constitutive model (incompressible Newtonian model) represented by the stress terms (viscous plus pressure) is derived by enforcing conservation of angular momentum. Let us denote u ! meas the 4D Flow measurement field. We assume that there exist a field w ! , that satisfies: with ∂Ω being the whole boundary of Ω. By writing (1) in weak form, and using relations (2)-(4), we can formulate the following weak problem: Find The following remarks are in order.
Remark 1. The left-hand-side of Problem (5) resembles the incompressible Navier-Stokes equation, up to two additional terms. Unfortunately, none of these terms are (semi-)positive-definite. Therefore, the analysis of the well posedness of this continuous problem (existence, uniqueness, time-stability) becomes challenging. However, those properties can be ensured at the discrete level by including adequate stabilization terms and constraints on the physical constants.
Remark 2. Equation (1) uses the so-called convective form of the advective term. While alternative forms of Equation (5) could be derived starting from other forms for the advection (e.g., conservative), the resulting discrete problem will need to be stabilized to become a solvable problem, leading to the same expression for the bilinear form. There will be, however, a difference in right-hand-side terms. There is, however, no particular reason to choose one above the other since all formulations are consistent with perfect (i.e., divergence-free) measurements.
where the second term in the right-hand-side is needed to enforce the compatibility with respect to the boundary condition (4). However, as will be show later on, this leads w ! to have in general larger values than using the divergence-free model, even in places with no or little measurement errors, reducing its potential utility.
Remark 4. The regularity assumptions required to solve for w ! are stronger than the usual perturbations in real data, which are often not additive in nature and present strong jumps in space and time. Therefore, we do not aim that w ! is capable of "correcting" the data (e.g., by computing u ! meas þ w ! ) as it may be the purpose of, for example, data assimilation approaches. The field w ! rather aims to provide a complementary metric to detect defects in the measurements.
Remark 5. The choice of homogeneous Dirichlet boundary conditions is based on its consistency with the case of perfect measurements. In contrast, natural boundary conditions are not consistent: if w ! ¼ 0 ! everywhere the pressure p cannot become the physical pressure on boundaries where homogenous Neumann boundary condition is enforced.
Remark 6. Note that the so-called "Stokes estimator" (STE) method for pressure reconstruction 16 is recovered dropping the first three terms of the left-hand-side.

| Stabilized finite element formulation
The tetrahedral mesh with characteristic element size h obtained from the segmented medical image is denoted by Ω h , which is the discrete domain over we define the following functions spaces In order to aim for the clinical applicability, it is crucial to use fast and robust numerical schemes. For the spatial discretization, we adopt V h and Q h as spaces for w ! and p, respectively, using stabilized finite elements to ensure solvability.
For the time discretization, we consider a backward Euler method with fixed time step τ to avoid GCL-type conditions. In order to avoid a root-finding problem at each time step the non-linear term on w ! will be treated semi-implicitly.
The resulting fully discrete stabilized formulation reads as follows. Given The stabilized bilinear form is defined as: being the bilinear form associated to the non-stabilized weak form of (5), while the convection stabilization term is given by and the pressure stabilization term as with δ > 0 some user-defined parameter. Finally, the right-hand-side is given by meas . Additional theoretical properties of Problem (5) in terms of its solvability and time stability of the solution are shown in Appendix A.
Remark 7. The stabilization term S k conv is in general not consistent with the solution of (5). However, it is consistent when the measurements are perfect since, in such case w ! ¼ 0 and r Á u ! meas ¼ 0. The stabilization term S k press is weakly consistent due to the inclusion of h 2 .
3 | NUMERICAL EXPERIMENTS USING SYNTHETIC DATA

| Procedure for synthetic data generation
All testcases were created in the following way: • A reference flow u ! is created using a finite element solver setup on an unstructured mesh of an aorta with a coarctation using high spatial and temporal resolutions. Physiologically relevant boundary conditions were set. The geometry with boundary conditions and resulting velocity field are shown in Figure 1. The simulation details are given in Appendix B. • The results of the simulation where interpolated in time with a time step of 0.03 s, hence downsampled by a factor of 30. • In order to simulate the 4D Flow MRI acquisition, a complex magnetization field was computed for every component of the reference velocity as: , with u 1 , …, u 3 the time undersampled velocity fields. The value of ϕ 0 was assumed to be constant in space and time and equal to 7:5 Á 10 À2 rad, and m 0 assumed constant and equal to 0:5 on the volume. The venc parameter setup will be explained in short. The background magnetization field was computed without any dependence on the velocity as M 0 ¼ m 0 exp iϕ 0 ð Þ. • Each magnetization field M 0 , …, M 3 was interpolated into a fine voxel-mesh with an element size of h b ¼ 1 mm. Then, a smoothing filter on space via convolution with a Gaussian kernel with a standard deviation of 5h b was applied, in order to take into account partial volume effects. The filtered magnetization was interpolated into a grid mesh with F I G U R E 1 Boundary conditions for the simulation with the reference solution obtained for generate synthetic data sets. In (A) boundary conditions, (B) Velocity vector field in the reference mesh at peak systole the desired image resolution of 2 Â 2 Â 2 mm 3 around the reference mesh. Such interpolation was performed using a piecewise linear Lagrangian interpolator. Then, the four components were arranged as a multidimensional array, which will serve to create the images, denoted as M 0 , M 1 , M 2 , M 3 . • We perturbed the magnetizations adding noise and aliasing by: The venc parameter was set lower than the peak true velocity. Therefore, when the velocity (in absolute value) exceeds the venc, the reconstructed velocity will be wrapped. Two venc parameters as the 150% and 70% of the maximum reference velocity were chosen, resulting in the values of 123 and 57 cm/s, respectively. The magnetization M j , j ¼ 0, …, 3 includes Gaussian noise, having real and imaginary parts independent noise realizations with standard deviation 0.25. This results in a velocity noise with a variance of 17.73% of the maximum reference velocity, in the high venc case, and a variance of 11.41% in the lower venc case. • Then, the 4D Flow measurements are given by with = representing an element-wise division of the arrays and f M j corresponding to the magnetization perturbed with the noise. Here, the time step index has been omitted in the description for the sake of readability.
• An image mask is created from the reference simulation on the uniform-grid mesh. Then, a new semi-structured tetrahedral mesh following the aortic shape is created using the algorithm reported in Reference 17. The velocity is defined as a P 1 finite element field on such mesh to visualize the results and quantify the errors to the reference solution interpolated to the same mesh. The different velocity measurements generated are shown in Figure 2, first column from left to right.

| Results
Problem (7) was solved using synthetic velocity measurements. In all the cases, δ ¼ 1 was used for the pressure stabilization term, chosen as the smallest possible value regarding a local neighborhood insensitive to δ.
The results are shown in Figure 2 at the time of peak systole. For every data set, the two alternative formulations defined by Equations (3) and (6), referred now as Model A and Model B were tested, resulting in the fields w values below 30% of u ! meas . This could be attributed to the possibility that larger errors in the gradient of the measurements are more challenging to handle with the proposed stabilized finite element scheme. The divergence field shows similar behavior. In contrast, w ! B presents very large values in a few spatial locations in spite that the measurements are perturbed only by the spatial blurring. In the pure aliased case (second row, (e)-(h)), the divergence appears to be larger after the coarctation, where aliasing is present. However, in the ascending aorta it appears to detect aliasing only in a few places. This could be related by the fact that in the case of a "straight" flow, a wrapped velocity field does maintain the divergence. In contrast, w ! on both models considerably grows where aliasing occurs, compared to the nonperturbed case, but not as much as δu ! . Also, a coupling among the components of the field appears, therefore w ! on both models, does not point exactly as δu ! . Note that also not all points with aliasing are detected, most likely due to the homogenous Dirichlet boundary condition. In the case of w ! B , it turns out that its magnitude grows considerably in the regions where no aliasing is present, while the growth of w ! A does concentrate around the aliasing. For the pure noise case, w ! B presents higher values than w ! A . Nevertheless, both present a similar "random" behavior as δu but with lower magnitude in the case of w ! A due to its higher regularity. For the case with noise and aliasing, a combination of the two aforementioned behaviors is obtained. Note that since the velocity-to-noise ratio is proportional to the venc, in this case it is smaller than what is shown in the "pure noise scenario" with a larger venc. The divergence field shows again a small sensitivity to aliasing and exhibits an overall increment in the presence of noise.
Regarding the differences between both corrector fields, in summary w ! B presents larger values in artifact-free zones. w ! A , on the other hand, can better detect localized perturbations in the measurements. Moreover, as mentioned in Remark 5 and as it can be appreciated on Figure 2, the field w ! does not match δu due to its higher regularity than the perturbed measured velocities. Consequently, since the purpose of w ! is to simply detect faults in the data, we recommend and adopt Model A for the rest of the manuscript. Additionally, a correlation study of the fields k w ! k and j r Á u ! j for every case with respect to k δu ! k was performed. Only the results obtained by the Model A were used. Pearson and Spearman correlation coefficients were computed for each field at the moment of peak systole. Moreover, a rank of significant association was computed between the values of the fields using the Maximal Information Coefficient estimator (MIC e ) from Reference 18. All these values are shown in Tables 1-4. The k w ! k field shows better correlation with k δu ! k than j r Á u ! meas j in cases where aliasing is present in the measurements. Without aliasing, for both fields the correlation is much lower, having k w ! k a slightly better performance.
The MIC e coefficient shows higher level of significance between the fields when aliasing is present as well. orifice coarctation made of Technyl was placed in the descending aorta (for further details of the setup and the phantom see 19,20 ). A blood mimicking fluid made with 60% water and 40% glycerol (Orica Chemicals, Watkins, CO) was used in the system. The fluid has a density of 1:119 g=cm 3 , dynamic viscosity of 0:0483 P and T1 value of 900 ms, which are in the range of values for human blood. The acquisition was performed with a venc of 350 cm/s and using a Cartesian sampling sequence with no k-space undersampling involved. In MRI, the noise level of the image increases when decreasing the voxel size. Therefore, three isotropic voxel sizes (coarse: 2.5 mm, mid: 2.0 mm and fine: 1.5 mm) were acquired in order to investigate the results of w ! in front of different resolutions and SNR levels. Figure 3 shows the 4D Flow measurements together with their result for w ! at the moment of peak systole. Also the divergence of the measurements are included. First, note that for all three resolutions w ! tends to grow in zones where the flows becomes convection dominated, and not necessarily where the velocities are high, see Figure 3A and B. For instance, in the three regions of high velocity, the inflow tube and the descending aorta post stenosis, the inflow tube presents a negligible w ! value, where we F I G U R E 3 Results for aortic phantom data, for three spatial resolutions of the measurements: measurements (left), (mid) w ! , and (right) divergence of the measurements expect the flow to be closer to a Stokes flow. This observation agrees with the results obtained in the synthetic case. However, real measurements have an additional source of error, namely the measurement technique assumes that the velocity field is constant within the time window of observation, in this case 40 ms. This induces to larger errors in the data that are not present in the synthetic case, possibly explaining the higher w ! values observed in the coarse voxel case. When decreasing the voxel size, as expected from the reduction of the signal-to-noise ratio, w ! increases with the appearance of artifacts in the measurements. This occurs in the ascending aorta, in the stenosis, and in the inferior part of the arch. For all voxel sizes the divergence shows an increment post stenosis, with much lower values in the rest of the domain. Therefore, at least in this experiment the divergence of the measurements seems not to be capable of detecting spatial regions with decreased measurement quality, it only captures the regions where the velocity is larger.

| Subjects
Two healthy volunteers were scanned in the previously described 1:5 Tesla scanner using a 4-channel torso coil. The local committee approved the study, and informed written consent was obtained from the participants. The acquisition parameters were: FOV in the range of 192 Â 192 Â 162 and 224 Â 224 Â 162 mm 3 , voxel resolution of 2:0 Â 2:0 Â 2:5 mm 3 , temporal resolution dt ¼ 34 ms, venc ¼ 150 cm=s, 25 cardiac phases, flip angle of 6 hand TR=TE ¼ 4:9 ms=2:9 ms. From the resulting data, only the aorta was segmented in order to apply our in-house mesh generation algorithm. Afterwards, the problem was solved assuming a blood density of ρ ¼ 1:2 gr=cm 3 and a dynamic F I G U R E 4 4D Flow measurements, the corresponding w ! , velocity and divergence fields for the volunteers at the time of peak systole viscosity of μ ¼ 0:035 P, same values taken for the simulation with synthetic data. A study of the impact of these parameters on the solution w is presented in Appendix C. Figure 4 shows the 4D Flow measurements, their divergence and the resulting w ! at peak systole. In all the cases w ! grows when flow artifacts appear in the measurements. Volunteer 1 shows a highly regular blood flow which results in a small w ! . The divergence presents no large peaks. This suggests that the measurements are mostly perturbed by noise. Only in the ascending aorta w ! takes slightly higher values, probably due to larger convective effects in that region. The divergence field shows a concentrated spike on the top of the aortic arch likely to be caused by a boundary artifact, which is not detected by w ! due to the homogenous boundary condition.
Volunteer 2 shows a case where the velocity measurements are more perturbed, in particular on three different locations. In the ascending aorta, closer to the heart, the measurements presents discontinuities at certain locations. This is confirmed by the high values of w ! , while the divergence seems almost not to be perturbed, compared to other regions.
In the aortic arch, all three fields appear to be perturbed. A last small region with larger w ! and divergence values appear in the distal part of the descending aorta.

| CONCLUSION
We presented a new mathematical model-including a tailor-made discretization-to detect imperfections in full-field velocity measurements as 4D Flow MRI can obtain it. The derived model uses data consistency with the incompressible Navier-Stokes equations, leading to a new vector field w ! , used as an indicator to check the compatibility of the data with the physics induced by Navier-Stokes.
Synthetic data experiments show the ability of the proposed approach to detect typical artifacts in 4D Flow MRI images, such as noise and aliasing, more robustly than computing the divergence of the measurements. With real 4D Flow MRI, we showed how w ! increases as the data quality decreases. That indicator mark regions where the data could be misrepresenting the blood flow, which is a valuable information when further flow quantification needs to be performed. In our experiments, we observed how some of these errors were not detected using the divergence as indicator.
At representative spatial and temporal resolutions of the data, computations of the proposed finite element discretization results take about 50 s without any parallelization using standard personal computers. That makes its clinical application feasible, since most of the burden lies in the blood lumen segmentation, where several flow markers are quantified. Our work still presents some limitations. Concerning the numerical scheme, more advanced discretizations could be investigated to reduce the vector values arising from the discretization itself, particularly in convection-dominated regions. Moreover, the proposed discretization of w ! cannot fully capture strong discontinuities in the measurements, for example, in aliasingcontaminated measurements, mainly due to the regularity of the solution imposed by the discrete spaces used in this work. Therefore future work could consist in investigating the application of Discontinuous Galerkin approaches. Another limitation of this study concerns the measurement generation in the synthetic data example. The analysis did not consider variability, due to a time blurring effect after several cardiac cycles, which occurs in the MRI measurement process due to sequential filling of the frequency space, respiratory motion artifacts, among others.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available on request from the corresponding author. The data are not publicly available due to privacy or ethical restrictions. First, note that u ! k meas H 1 Ω h ð Þ ½ 3 for all k since the velocity data is bounded and interpolated to P 1 finite elements.
Denote the space of the whole solution vector W ¼ V h Â Q h , then we can prove the following results.
Proof. The first term defines a norm in V h and the second term is a seminorm in W , so k Á k W is a seminorm in W . It remains to proof that k v ! , q Indeed, the first two terms have to be zero and therefore v ! ¼ 0 ! . Therefore, rq is a constant and also zero-valued since q Q h .

Proposition 1.
For v ! , w ! ℝ d , C ℝ dÂs the following relation holds: with k k m denoting the ℓ m -norm in ℝ d .
Proof. Directly from the proposition's definition.
Lemma 2. There exists α > 0 such that: Proof. Using standard arguments, Poincaré's inequality and Lemma 1 the following relation holds: under condition (A2). By adding the remaining term S k press w ! , p; w ! , p , relation (A1) follows directly from the definition of the norm. Proof. Since There exists a unique solution in W of Problem (7) under condition (A2) for all k > 0.
, the bilinear and linear forms fulfill the requirements of the Lax-Milgram Theorem (i.e., Lemmas 2-4) for k ≥ 1.

A.2 | TIME STABILITY OF THE DISCRETE SOLUTION
We can furthermore prove the following energy balance: under the condition Proof. Testing (7) with v ! ¼ w ! k and q ¼ p k and using similar arguments as in Lemma 2, it is obtained where the two first terms in the right-hand-side come from the continuous problem, and the two last terms are dissipative due to the numerical scheme. Bounding the former using Poincaré's inequality and the later by zero we obtain which combined with condition (A4) leads to relation (A3).
Remark 8. Condition (A4) is very conservative and in the test cases we observe that the condition is not satisfied. However, in our computations, we have never obtained instabilities, even in situations with temporal jumps in the data as when velocity aliasing appears.
The initial conditions are set as u ! 0 ¼ 0 and π 0 ℓ ¼ 85 mmHg for ℓ ¼ 1, …, K, which correspond to approximately the periodic state of the 3D-0D system. The simulation was performed with a total time of T f ¼ 0:8 s with a timestep of dt ¼ 0:001 s.

A P P END I X C: VISCOSITY SENTIVITY STUDY
In order to know how sensitive is the corrector field w ! to variations of the dynamic viscosity parameter μ, we computed two fields using the values μ À ¼ 0:0175 P and μ þ ¼ 0:070 P with the 4D Flow data set acquired for the healthy volunteers. These values are in a 50% range of reference value for human blood used in this work (μ 0 ¼ 0:035 P), more concretely, we computed the H 1 -norm time curves of the corrector field over time, which are shown in Figure C1. Note that, as the second volunteer has relatively higher velocities, the effect of the viscosity will be "hidden" by the dominance of the convective forces within the fluid. It can be seen that higher differences start to appear during diastole, when the Reynolds number are relatively small compared with the starting of the cycle. This can be explained by the dominance of drag forces in a low velocity regime. In Volunteer 1, the relative difference of the H 1 -norm is no more than a 15% of the values obtained using μ 0 . The latter can be confirmed from the Figure C2, in which the fields are shown at t ¼ 0:6 s, corresponding to the time instant where differences are maximal. In conclusion, within a reasonable interval of viscosity values no major differences in the corrector fields are observed. Therefore, we can safely adopt the reference value μ 0 ¼ 0:035 P.
F I G U R E C 1 Comparison between norms of the corrector field under three different dynamic viscosity values F I G U R E C 2 Corrector fields at t ¼ 0:6 s. The arrow scale have been remained the same with respect to the scaling used in the rest of this work for a better comparison