On the stable estimation of flow geometry and wall shear stress from magnetic resonance images

We consider the stable reconstruction of flow geometry and wall shear stress from measurements obtained by magnetic resonance imaging. As noted in a review article by Petersson, most approaches considered so far in the literature seem not be satisfactory. We therefore propose a systematic reconstruction procedure that allows to obtain stable estimates of flow geometry and wall shear stress and we are able to quantify the reconstruction errors in terms of bounds for the measurement errors under reasonable smoothness assumptions. A full analysis of the approach is given in the framework of regularization methods. In addition, we discuss the efficient implementation of our method and we demonstrate its viability, accuracy, and regularizing properties for experimental data.


Introduction
Due to its clinical relevance, the estimation of wall shear stress in arteries, i.e., of the normal derivative of tangential velocity components at aterial walls, has attracted significant interest in the literature; see e.g. [7,12,15,16,17,19] and the references therein. Unfortunately, most of the approaches utilized so far suffer at least from one of the following artifacts: • wall shear stress is systematically underestimated for coarse data resolution; • for increasing data resolution the estimates become increasingly unstable.
Let us briefly discuss some reasons for these problems: Internal measurements are used in [15] to fit cubic polynomials to the flow profile. In the presence of a logarithmic turbulent boundary layer [9], this leads to a flattened approximation of the velocity near the boundary, thus underestimating the velocity derivative and overestimating the radius of the flow region. Due to a kink of the velocity profile at the boundary, a spline interpolation in the whole domain, as proposed in [19], leads to inaccurate representation of the velocity, in particular, near the boundary, which makes the estimate of wall shear stress unreliable. These observations are amplified by the fact that noise in the velocity measurements can be expected to be particularly high close to the boundary; compare with the data in Section 6 and the remarks in the appendix. Let us note that even for exact velocity data the evaluation of the wall shear stress is unstable with respect to the boundary location due to the discontinuity of the normal derivative of velocity at the boundary. These observations made the authors of [16] conclude that even in the absence of noise and for relatively simple velocity profiles, all methods evaluated were found to be impacted by considerable errors.
In this paper, we investigate the stable estimation of flow geometry, velocity, and wall shear stress by a problem adapted approach that overcomes the above difficulties and that allows for a rigorous stability and convergence analysis. The overall reconstruction problem will be decomposed into the following three natural steps: • estimation of the flow boundary from magnetic resonance images (geometry identification); • reconstruction of flow velocity from phase contrast images in a function class that ensures zero velocity at the estimated flow boundary (velocity estimation); • evaluation of the normal derivative of the velocity at the boundary (wall shear stress computation).
In principle, various methods for the solution of the three sub-problems are available. We here consider one particular approach that allows to conduct a full convergence analysis of the overall reconstruction process under reasonable smoothness assumptions. For the geometry identification problem, we utilize a parametric formulation which leads to a nonlinear inverse problem with a non-differentiable forward operator. We prove a conditional stability estimate and derive convergence rates for Tikhonov regularization under simple smoothness assumptions on the true geometry. The parametrization of the flow domain obtained in the first step is then used to formulate the velocity reconstruction problem on a reference domain, resulting in a linear inverse problem with additional data perturbation stemming from the inexact geometry representation. Again, a full convergence analysis of Tikhonov regularization is presented for this problem. By choosing sufficiently strong regularization terms in the first two steps, we obtain stability of the geometry and velocity reconstruction in strong norms that allow us to compute the third step in a stable way. Apart from a complete theoretical analysis of our approach, we also discuss the efficient numerical realization and demonstrate its viability by application to experimental data.

Geometry identification
For ease of presentation, we assume in the sequel that the flow geometry is cylindrical and that the flow field has the particular form (0, 0, u(x, y)), which allows us to restrict the considerations to a two dimensional setting; the extension to three dimensions will be discussed briefly in Section 5. Without loss of generality, we further assume that measurements are available in the domain D = (−1, 1) 2 , which we call the field of view.
For a given radius function R : [0, 2π] → R + , we define the domain parametrized by R, and we denote by Ω † = Ω R † the true flow geometry. Throughout the presentation, we will assume that for some constants 2 ≤ k ≤ 4 and 0 < r 0 < r 1 < 1; here H s per (0, 2π) is the subspace of periodic functions in the Sobolev space H s (0, 2π), s ≥ 0. Let us note that the above assumptions imply in particular that Ω † is of class C k−1,α , compactly embedded in D, and star shaped with respect to the origin.
The geometry identification from possibly perturbed magnetic resonance images can now be formulated as a nonlinear inverse problem with forward operator F : We further assume that the perturbations in the data are bounded by and that knowledge of the noise level δ is available. For the stable approximation of solutions to (3), we consider Tikhonov regularization with the functional We call R δ α ∈ D(F ) an approximate minimizer for regularization parameter α > 0, if Note that the operator F here is continuous but not differentiable and, therefore, standard convergence rate results about Tikhonov regularization, cf. [5, Chapter 10], do not apply. Instead, we utilize recent results on Tikhonov regularization in Hilbert scales under conditional stability [3,20], which allow us to prove the following assertions.
Then any approximate minimizer R δ α of the Tikhonov functional (5) with α = δ 4/k satisfies for 0 ≤ r ≤ 2, with a constant C that only depends on the norm R † H k (0,2π) of the true solution. The same estimates hold true, if α is chosen by a discrepancy principle, i.e., Proof. Note that X s = H s per (0, 2π), s ∈ R defines a Hilbert scale and Y = L 2 (D) is a Hilbert space. In view of the results in [3], it thus remains to establish a conditional stability estimate for the operator F . To do so, let us assume that R, R ∈ D(F ) and define r min (ϕ) = min{R(ϕ), R(ϕ)} and r max (ϕ) = max{R(ϕ), R(ϕ)}. Then This allows us to express the difference in the observations by i.e., the operator F satisfies a conditional Lipschitz stability estimate. The assertions of the theorem then follow from Theorems 2.1 and 2.2 in [3] with a = 0, s = 2, u = k, and γ = 1.
Remark 2.2. Note that only a simple smoothness condition for the function R † , and thus on the domain Ω † , was required to derive a quantitative convergence result here. If in addition, an upper bound R † H 4 (0,2π) ≤ C is available, then one obtains with δ R = Cδ 1/2 and we may assume that the bound δ R is known for further computations. The results of [3] even provide more general estimates of the form if regularization is performed in the norm of H s (0, 2π) for 2 ≤ s ≤ k/2 and with regularization parameter α = δ 2s/k . Hence, convergence rates arbitrarily close to one can, in principle, be obtained in any desired norm if the true solution R † is sufficiently smooth and the regularization norm is chosen sufficiently strong.

Velocity approximation
Let us recall from equation (1) the definition of a domain Ω R parametrized by a radius function R ∈ D(F ). Using the scaling transformation φ R : (r cos ϕ, r sin ϕ) → (r 0 r + (R(ϕ) − r 0 ) r η ) · (cos ϕ, sin ϕ), with η ≥ k ≥ 2 and k as in the previous section, we can express Ω R equivalently as image Ω R = φ R (B) of the unit disk B = {x ∈ R 2 : |x| < 1}. Some important properties of this transformation are derived in Appendix A. In the following, we assume that with known bound δ R ≤ C; see the discussion in Remark 2.2. For ease of notation, we write φ † = φ R † , φ δ α = φ R δ α , and denote by Ω † = φ † (B) and Ω δ α = φ δ α (B) the corresponding domains parametrized by the radius functions R † and R δ α , respectively. The velocity reconstruction from noisy velocity data u can then be formulated compactly as a linear inverse problem over the reference domain, i.e., with operator T : i.e., as simple embedding between Sobolev spaces. We assume that a bound on the data perturbations is available, where u † denotes the true velocity field to be reconstructed, and we require that u † is piecewise smooth and vanishes outside Ω † , i.e., Let us note that these assumptions are reasonable, if the flow domain Ω † is smooth.
Observe that information about zero velocity at the boundary of the flow domain is encoded explicitly into the definition of the operator T . The particular formulation (12) over the reference domain B is used for the following reason: Inexact knowledge about the flow domain is shifted to the data u • φ δ α and the operator T , therefore, is independent of the geometry reconstruction. This simplifies the analysis given in the following. Let us note that invertibility of the transformation φ δ α is guaranteed by Lemma A.4., which allows to associate to any approximate solution v of equation (12) a velocity field u = v • (φ δ α ) −1 in the physical domain. For the stable solution of the nonlinear inverse problem (12), we again consider Tikhonov regularization and we define the regularized approximations by v δ, α,β = arg min where For our further analysis, we will require two auxiliary results that we will present next. As a first step, we derive an estimate for the data error in the reference domain.
Lemma 3.4. Let the assumptions (11), (13) and (14) hold. Then and a constant C that can be computed explicitly. Proof. By the integral transformation theorem we have where we used Lemma A.4. to estimate the Jacobian J δ α of the transformation φ δ α in the second step. The last term can then be readily estimated by the bound (13) on the data noise. The first term on the right hand side can further be split as (14). From the smoothness of R δ α , R † , and u † , and using the bound on the difference of the inverse transformations (φ δ α ) −1 and (φ † ) −1 provided by Lemma A.5., we then obtain (10) of the transformation and continuous embedding. Using Lemma A.1. and assumption (11), we can control the area |Ω δ α \ Ω † | of the geometry mismatch by δ R , which allows us to bound the remaining term in the above estimate by R . Note that the constants C in the above estimates are generic and independent of and δ. A combination of the bounds then yields the assertion of the lemma.
Remark 3.5. Following Remark 2.2. and the arguments used in the proof above, the bound δ U in (16) is computable in terms of the data noise levels δ and in (4) and (13), if bounds on R † H 4 (0,2π) and u † H 3 (Ω † ) are available. For our further computations, we may thus assume that δ U is known.
As a next step, we interpret standard smoothness assumptions on u † in terms of the operator T , which will allow us to utilized simple source conditions below. Lemma 3.6. Let assumptions (11) and (14) hold and define v † : . Thus w = T * f is given as the unique solution of the boundary value problem with v = 0 and ∆v = 0 on ∂B.
at ∂B would be valid additionally; this can however not be expected in general. The limiting factor for the regularity index µ in the range condition, therefore, is the mismatch of the higher order boundary conditions. This could be circumvented by choosing a different equivalent norm on H 2 (B) ∩ H 1 0 (B); see [14] for details. From standard results about Tikhonov regularization in Hilbert spaces for linear inverse problems [5], we can now immediately conclude the following results.
Theorem 3.8. Let assumptions (11)-(14) hold and let v δ, α,β be defined by (15) with The same estimates hold for a-posteriori parameter choice by the discrepancy principle Corresponding bounds for the error in the velocity u can be obtained, in principle, by back transformation into physical domain and some elementary computations. Let us close this section with a remark indicating some natural generalizations.
Remark 3.9. If u † is sufficiently smooth and the functional in (15) is replaced by defined by T v = v and t ≥ 1 sufficiently large, one could, in principle, obtain any rate 2µ/(2µ + 1) sufficiently close to one. Further note that the data residual could also be measured in physical space. The regularized approximate solution is then defined by v δ, α,β = arg min This simply amounts to a change to an equivalent norm in the data space. A quick inspection of the arguments in the above proof reveals that the assertions of Theorem 3.8. remain valid also for this choice of regularization method, which is more easy to implement and will thus be used in our numerical tests in Section 6.

Computation of the wall shear stress
Let R ∈ D(F ) be a given radius function with associated transformation φ R and let v be a given velocity field defined on the reference domain B. For ease of notation, we assume that fluid viscosity is normalized, and define the associated wall shear stress by where J R is the Jacobian of φ R , and n R (ϕ) is the outward pointing unit normal vector at the corresponding point (R(ϕ) cos ϕ, R(ϕ) sin ϕ) ∈ ∂Ω R in the physical domain. Let us note that n R (ϕ) can be expressed explicitly by For ease of notation, we will identify ∂B with the interval (0, 2π) in the sequel. In addition, we again define symbols τ δ, Proof. We use the definition of τ R (v) and decompose the error into the three parts . Using Lemmas A.2. and A.3., and the embedding of Sobolev spaces, we obtain By the geometric arguments of Lemma A.5., the second term can be estimated by With the help of the results of the previous section, the third term, which measures the amplification of the velocity approximation error, can be bounded by The assertion of the theorem then follows by combination of these estimates.
Remark 4.11. Let us note that, following the considerations of Remark 2.2. and 3.9., one could in principle again obtain convergence rates arbitrarily close to one, if the true flow geometry and velocity are sufficiently smooth and the regularization terms in the Tikhonov functionals (5) and (15) are chosen sufficiently strong. The main observation of the previous theorem therefore is, that it is possible to obtain a quantitative estimate under reasonable smoothness assumptions.

Remarks on the extension to three dimensions
For a three-dimensional flow, the wall shear stress is a tensor defined by where µ is the dynamic viscosity, u is the velocity vector, and n the outer unit normal vector at the vessel wall. If appropriate measurements of the density and of all three velocity components are available, then the reconstruction approach and the theoretical results presented in the previous sections can be generalized almost verbatim to the three dimensional setting; only the computational realization becomes more complicated. Note that the geometry and velocity reconstruction naturally decompose into a sequence of two-dimensional inverse problems for the individual cross-sections for fixed coordinate z in the flow direction. To ensure continuity of the reconstructions with respect to the zcoordinate, one has to employ regularization also in the z-direction, which fully couples the problems for the individual cross-sections. The additional computational complexity due to this coupling can be overcome by a Kaczmarz strategy [8,11], which allows to reduce the numerical solution to the iterated solution of two-dimensional problems for the individual cross-sections. A detailed investigation of these computational aspects is, however, not in the scope of the current paper and will be given elsewhere.

Numerical validation
In order to demonstrate the viability of the proposed approach, we now report about the reconstruction of flow geometry, flow velocity, and wall shear stress from experimental data obtained in a clinical whole-body 3 Tesla magnetic resonance imaging scanner (Prisma, Siemens Healthcare, Erlangen) at the University Medical Center Freiburg. As a test case, we consider the flow of water through a cylindrical pipe with constant diameter d = 25.885mm at a constant flow rate of about 6l/min and a temperature of about 22 • C; details about the experimental setup are described in [1]. Let us note that the flow conditions lead to a Reynolds number of about 5300 and thus amount to a turbulent flow with steep velocity gradients in the boundary layer. For our reconstruction procedures, we utilize magnetic resonance images of density and the axial flow velocity acquired with different resolutions corresponding to in plane voxel sizes of h = 1.5mm to h = 0.3mm. Due to the axisymmetic geometry, we expect an almost radially symmetric flow, but this a-priori knowledge will not be utilized in our computational tests. Up to a simple translation, the exact geometry Ω † = φ † (B) of the pipe is known here. A reference solution u ref for the flow velocity can thus be computed by numerical simulation [9] and will be used for comparison with our results. Let us emphasize that this reference solution represents a time averaged velocity field, in which all temporal fluctuations are suppressed. The experimental data, on the other hand, contain such turbulent fluctuations; see Figure 2. In addition, also the flow conditions, e.g., temperature and flow rate, do not match exactly with the simulation data. One therefore cannot expect to get a perfect fit to the simulated reference velocity data.

Geometry identification
For the geometry identification, we utilize the standard magnetic resonance images of the density. After a simple scaling procedure, see Appendix B for details, the data can be interpreted as where V i denotes the ith voxel of the measurement array of size h×h. The action of the forward operator for our computational tests is then defined as where H γ (x) = 1 π arctan x γ + 1 2 , γ > 0, is a smooth approximation of the Heaviside function. In our computational tests, we thus minimize the Tikhonov functional over the finite dimensional space of radius functions Due to the smoothing with γ > 0, the forward operator F h,γ is continuously differentiable and a projected Gauss-Newton method [4,8] can be used for the minimization process. For our computational tests, the center of the coordinate system was chosen as the barycenter of the measurement data m δ and the true radius was R † = 12.9425mm. In Table 1 As expected intuitively, reconstruction errors are slightly decreasing with increasing data resolution. Further note that the reconstructions are stable with respect to the choice of the regularization parameter and the optimal regularization parameters, i.e., those for which the error is minimal, are practically independent of the data resolution.
In Figure 1, we display one of the data sets used for our computations, together with the reconstructed geometry and the radius functions R † and R δ α obtained in our numerical tests. The maximal error in the reconstructed radius is about 0.15mm, which is substantially less than the voxel size h = 1mm of the data. This illustrates that subpixel resolution can be obtained by the proposed geometry reconstruction procedure.

Velocity approximation
We now turn to the reconstruction of the velocity field, for which we utilize phasecontrast magnetic resonance imaging data [12,17]; also see Appendix B. Let Ω δ α = φ δ α (B) denote the approximation of the flow domain obtained with the radius function R δ α reconstructed in the first step with α chosen as the best regularization parameter according to Table 1. For voxels V i lying at least partially in the flow domain, i.e., with |V i ∩ Ω δ α | > 0, the measurements can interpreted as The remaining voxels only contain information about the noise and will not be used in the reconstruction; see Figure 2 and the remarks given in Appendix B. Following Remark 3.9., we define the corresponding forward operator mapping to the physical domain by T h : Note that only voxels V i with |V i ∩ Ω δ α | > 0 are required in the definition of T h . The regularized approximation for the velocity field in the reference domain is then defined as the minimizer of the Tikhonov functional v δ, α,β = arg min  Table 2, we display the reconstruction errors obtained with our algorithm for different data resolutions h and regularization parameters β.  Recall that u ref , which was obtained by numerical simulation, corresponds to a time averaged velocity field and therefore does not contain temporal fluctuations that are present in the measurements; see Figure 3. This explains the relatively large errors in Table 2 and their independence of the data resolution. The optimal regularization parameters β are again independent of the voxel size h used in the data acquisition.
A quick comparison with the raw data, depicted in Figure 2, with the plots in Figure 3 shows that the reconstructed velocity field correctly reproduces some turbulent fluctuations which, due to time averaging, are not present in the simulated reference velocity field. Apart from these differences, the reconstruction is in good agreement with the reference solution.

Computation of the wall shear stress
As a final step in our numerical tests, we now utilize the reconstructed radius function R δ α and velocity fields v δ, α,β to compute the approximation τ δ, α,β for the wall shear stress by (19). In Table 3, we display the results obtained for optimal α and β = 1.6 · 10 −4 . Let us note that most part of the error stems from perturbations in higher modes, which h 1.00 0.75 0.60 0.43 0.30 0.0789 0.0385 0.0391 0.0264 0.0405 Table 3. Relative errors τ δ, α,β − τ † L 2 (0,2π) in the reconstruction of the wall shear stress for different data resolutions resp. voxel size h.
can be suppressed efficiently by application of a low-pass filter. In Figure 4, we plot the reconstruction of the wall shear stress τ δ, α,β and its constant approximation τ δ, α,β against the reference value τ ref obtained from numerical flow simulation [9]. Note that the average wall shear stress τ δ, α,β is in very good agreement with the reference value τ ref obtained for the simulated data. The local variations in the reconstruction τ δ, α,β can be explained by the turbulent variations of velocity in the data; compare with the plots in Figure 2 and 3.

Discussion
As reviewed by Petersson [16], the stable and accurate estimation of wall shear stress from magnetic resonance imaging data, is a delicate issue and most reconstruction approaches reported in literature seem not to work properly. In this paper, we therefore considered a systematic approach for the estimation of wall shear stress from magnetic resonance images of density and velocity, for which stability and convergence could be established under simple and realistic smoothness assumptions on the flow geometry and velocity. The theoretical results were validated by numerical tests for experimental data which demonstrate that wall shear stress can be estimated from magnetic resonance imaging data with relative errors of a few percent and practically independent of the data resolution; this is in stark contrast to the results of Stalder [19].
Let us note that stable wall shear stress estimates can also be obtained via empirical Moody charts [13] or the Clauser plot method [2,18], which are, however, limited to axisymmetric geometry or fully developed turbulent flow. Numerical simulations could also be used, in principle, to compute wall shear stress estimates [9], but precise knowledge about the rheological properties of the fluid are required. In contrast to these approaches, the method considered in this paper is generally applicable and, therefore, seems most appropriate for application in a clinical context.
As a next step, we estimate differences in the normal vector n R defined in (20), in terms of differences in the radius function.
The following result ensures smoothness of the transformations φ R whenever the radius function R is sufficiently smooth. (10). Then (10). Moreover, Proof. The continuity of φ 1 is obvious and the Jacobian J 1 = J R 1 obtained by derivation with respect to coordinates x = (r cos φ, r sin φ) reads Since η ≥ 2, the Jacobian can be seen to be continuous also in r = 0, and we have J 1 L ∞ (B) ≤ C ( R 1 L ∞ (0,2π) + R 1 L ∞ (0,2π) ) ≤ C R 1 H 2 (0,2π) .
This shows the first estimate for k = 1. The assertion for higher order derivatives of φ 1 follow in a similar way. From the formula (10), one can see that φ R is affine linear in R.
Hence the second estimate follows directly from the first.
As a next step, we show that the transformations φ R defined in (10) are invertible.
Using the previous results, we can also bound differences in the inverse mappings.
Lemma A.5. Let R 1 , R 2 ∈ D(F ) and assume that R 1 ∈ H 3 per (0, 2π). Furthermore, let ψ 1 , ψ 2 denote the corresponding inverse transformations with Jacobians J ψ 1 , J ψ 2 . Then and the difference in the Jacobians can be bounded by Since the inverse Jacobian of φ 1 has the representation J −1 φ 1 =J φ 1 / det J φ 1 , whereJ φ 1 is a rearrangement of J φ 1 and all expressions are continuously differentiable, one can verify without difficulty that J −1 φ 1 ∈ W 1,∞ (B) with the associated norm bounded in terms of η, r 0 , and R 1 H 3 (0,2π) . Therefore we arrive at (A.8) by This completes the proof of the second estimate of the lemma.

Appendix B. Interpretation of the data
Let us briefly comment on the physical interpretation and the preprocessing of the experimental data used for the reconstructions in Section 6.

Appendix B.1. Magnitude data
The magnitude raw data m δ raw represent integral means of the proton density over voxels; this values are additionally perturbed by data noise. From a histogram of the magnitude data, see Figure B1, one can deduce to peak values m 0 and m 1 , which are used for scaling of the data.  (1, x)). These measurements then represent an approximation for the characteristic function of the flow domain.