On the stable estimation of flow geometry and wall shear stress from magnetic resonance images
OOn the stable estimation of flow geometry and wallshear stress from magnetic resonance images
H. Egger and G. Teschner
AG Numerik, TU Darmstadt, Dolivostraße 15, 62393 Darmstadt, GermanyE-mail: egger,[email protected]
Abstract.
We consider the stable reconstruction of flow geometry and wall shearstress from measurements obtained by magnetic resonance imaging. As noted in areview article by Petersson, most approaches considered so far in the literature seemnot be satisfactory. We therefore propose a systematic reconstruction procedure thatallows to obtain stable estimates of flow geometry and wall shear stress and we areable to quantify the reconstruction errors in terms of bounds for the measurementerrors under reasonable smoothness assumptions. A full analysis of the approachis given in the framework of regularization methods. In addition, we discuss theefficient implementation of our method and we demonstrate its viability, accuracy,and regularizing properties for experimental data.
Keywords : wall shear stress estimation, magnetic resonance imaging, ill-posed problem,Tikhonov regularization, conditional stability
1. Introduction
Due to its clinical relevance, the estimation of wall shear stress in arteries, i.e., ofthe normal derivative of tangential velocity components at aterial walls, has attractedsignificant interest in the literature; see e.g. [7, 12, 15, 16, 17, 19] and the referencestherein. Unfortunately, most of the approaches utilized so far suffer at least from oneof 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 usedin [15] to fit cubic polynomials to the flow profile. In the presence of a logarithmicturbulent boundary layer [9], this leads to a flattened approximation of the velocitynear the boundary, thus underestimating the velocity derivative and overestimating theradius of the flow region. Due to a kink of the velocity profile at the boundary, a splineinterpolation in the whole domain, as proposed in [19], leads to inaccurate representationof the velocity, in particular, near the boundary, which makes the estimate of wall shear a r X i v : . [ m a t h . NA ] D ec n stable estimation of flow geometry and wall shear stress even in theabsence of noise and for relatively simple velocity profiles, all methods evaluated werefound 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 abovedifficulties and that allows for a rigorous stability and convergence analysis. The overallreconstruction problem will be decomposed into the following three natural steps: • estimation of the flow boundary from magnetic resonance images (geometryidentification); • reconstruction of flow velocity from phase contrast images in a function class thatensures zero velocity at the estimated flow boundary (velocity estimation); • evaluation of the normal derivative of the velocity at the boundary (wall shear stresscomputation).In principle, various methods for the solution of the three sub-problems are available. Wehere consider one particular approach that allows to conduct a full convergence analysisof the overall reconstruction process under reasonable smoothness assumptions.For the geometry identification problem, we utilize a parametric formulationwhich leads to a nonlinear inverse problem with a non-differentiable forward operator.We prove a conditional stability estimate and derive convergence rates for Tikhonovregularization under simple smoothness assumptions on the true geometry. Theparametrization of the flow domain obtained in the first step is then used to formulatethe velocity reconstruction problem on a reference domain, resulting in a linearinverse problem with additional data perturbation stemming from the inexact geometryrepresentation. Again, a full convergence analysis of Tikhonov regularization ispresented for this problem. By choosing sufficiently strong regularization terms in thefirst two steps, we obtain stability of the geometry and velocity reconstruction in strongnorms that allow us to compute the third step in a stable way. Apart from a completetheoretical analysis of our approach, we also discuss the efficient numerical realizationand demonstrate its viability by application to experimental data.
2. Geometry identification
For ease of presentation, we assume in the sequel that the flow geometry is cylindricaland that the flow field has the particular form (0 , , u ( x, y )), which allows us to restrictthe considerations to a two dimensional setting; the extension to three dimensions will n stable estimation of flow geometry and wall shear stress D = ( − , , which we call the field of view .For a given radius function R : [0 , π ] → R + , we define the domainΩ R = { ( r cos ϕ, r sin ϕ ) : 0 ≤ r < R ( ϕ ) , ≤ ϕ ≤ π } , (1)parametrized by R , and we denote by Ω † = Ω R † the true flow geometry. Throughoutthe presentation, we will assume that R † ∈ H kper (0 , π ) with r ≤ R ( ϕ ) ≤ r (2)for some constants 2 ≤ k ≤ < r < r <
1; here H sper (0 , π ) is the subspace ofperiodic functions in the Sobolev space H s (0 , π ), s ≥
0. Let us note that the aboveassumptions imply in particular that Ω † is of class C k − ,α , compactly embedded in D ,and star shaped with respect to the origin.The geometry identification from possibly perturbed magnetic resonance imagescan now be formulated as a nonlinear inverse problem F ( R ) = m δ on D, (3)with forward operator F : D ( F ) ⊂ H per (0 , π ) → L ( D ), R (cid:55)→ χ Ω R defined on thedomain D ( F ) = { R ∈ H per (0 , π ) : r ≤ R ( ϕ ) ≤ r } . We further assume that theperturbations in the data are bounded by (cid:107) χ Ω † − m δ (cid:107) L ( D ) ≤ δ (4)and that knowledge of the noise level δ is available. For the stable approximation ofsolutions to (3), we consider Tikhonov regularization with the functional J δα ( R ) = (cid:107) F ( R ) − m δ (cid:107) L ( D ) + α (cid:107) R (cid:107) H (0 , π ) . (5)We call R δα ∈ D ( F ) an approximate minimizer for regularization parameter α >
0, if J δα ( R δα ) ≤ inf R ∈ D ( F ) J ( R ) + δ . (6)Note that the operator F here is continuous but not differentiable and, therefore,standard convergence rate results about Tikhonov regularization, cf. [5, Chapter 10], donot apply. Instead, we utilize recent results on Tikhonov regularization in Hilbert scalesunder conditional stability [3, 20], which allow us to prove the following assertions. Theorem 2.1.
Assume that R † ∈ H kper (0 , π ) ∩ D ( F ) , ≤ k ≤ . Then anyapproximate minimizer R δα of the Tikhonov functional (5) with α = δ /k satisfies (cid:107) F ( R δα ) − χ Ω † (cid:107) L ( D ) ≤ Cδ and (cid:107) R δα − R † (cid:107) H r (0 , π ) ≤ Cδ − r/k (7) for ≤ r ≤ , with a constant C that only depends on the norm (cid:107) R † (cid:107) H k (0 , π ) of the truesolution. The same estimates hold true, if α is chosen by a discrepancy principle, i.e., α = max { α − n : n ≥ such that (cid:107) F ( R δα ) − m δ (cid:107) L ( D ) ≤ δ } . (8) n stable estimation of flow geometry and wall shear stress X s = H sper (0 , π ), s ∈ R defines a Hilbert scale and Y = L ( D ) is aHilbert space. In view of the results in [3], it thus remains to establish a conditionalstability estimate for the operator F . To do so, let us assume that R, (cid:101) R ∈ D ( F ) anddefine r min ( ϕ ) = min { R ( ϕ ) , (cid:101) R ( ϕ ) } and r max ( ϕ ) = max { R ( ϕ ) , (cid:101) R ( ϕ ) } . Then | F ( R ) − F ( (cid:101) R ) | ( r cos ϕ, r sin ϕ ) = (cid:40) , for r min ( ϕ ) ≤ r ≤ r max ( ϕ ) , . This allows us to express the difference in the observations by (cid:107) F ( R ) − F ( (cid:101) R ) (cid:107) L ( D ) = (cid:90) π (cid:90) r | F ( R ) − F ( (cid:101) R ) | r dr dϕ = (cid:90) π (cid:90) r max ( ϕ ) r min ( ϕ ) r dr dϕ = 12 (cid:90) π r max ( ϕ ) − r min ( ϕ ) dϕ = 12 (cid:90) π | R ( ϕ ) − (cid:101) R ( ϕ ) | | R ( ϕ ) + (cid:101) R ( ϕ ) | dϕ ≥ (cid:90) π | R ( ϕ ) − (cid:101) R ( ϕ ) | r | R ( ϕ ) − (cid:101) R ( ϕ ) || r − r | dϕ = r r − r (cid:107) R − (cid:101) R (cid:107) L (0 , π ) . Hence (cid:107) R − (cid:101) R (cid:107) L (0 , π ) ≤ C (cid:107) F ( R ) − F ( (cid:101) R ) (cid:107) L ( D ) for all R, (cid:101) R ∈ D ( F ), i.e., the operator F satisfies a conditional Lipschitz stability estimate. The assertions of the theorem thenfollow from Theorems 2.1 and 2.2 in [3] with a = 0, s = 2, u = k , and γ = 1. (cid:3) Remark 2.2.
Note that only a simple smoothness condition for the function R † , andthus on the domain Ω † , was required to derive a quantitative convergence result here.If in addition, an upper bound (cid:107) R † (cid:107) H (0 , π ) ≤ C is available, then one obtains (cid:107) R δα − R † (cid:107) H (0 , π ) ≤ δ R , (9)with δ R = Cδ / and we may assume that the bound δ R is known for furthercomputations. The results of [3] even provide more general estimates of the form (cid:107) R δα − R † (cid:107) H r ≤ Cδ − r/k , ≤ r ≤ s, if regularization is performed in the norm of H s (0 , π ) for 2 ≤ s ≤ k/ α = δ s/k . Hence, convergence rates arbitrarily close to onecan, in principle, be obtained in any desired norm if the true solution R † is sufficientlysmooth and the regularization norm is chosen sufficiently strong.
3. Velocity approximation
Let us recall from equation (1) the definition of a domain Ω R parametrized by a radiusfunction R ∈ D ( F ). Using the scaling transformation φ R : ( r cos ϕ, r sin ϕ ) (cid:55)→ ( r r + ( R ( ϕ ) − r ) r η ) · (cos ϕ, sin ϕ ) , (10)with η ≥ k ≥ k as in the previous section, we can express Ω R equivalently asimage Ω R = φ R ( B ) of the unit disk B = { x ∈ R : | x | < } . Some important propertiesof this transformation are derived in Appendix A. In the following, we assume that R † , R δα ∈ D ( F ) , R † ∈ H per (0 , π ) , and (cid:107) R δα − R † (cid:107) H (0 , π ) ≤ δ R , (11) n stable estimation of flow geometry and wall shear stress δ R ≤ C ; see the discussion in Remark 2.2. For ease of notation, wewrite φ † = φ R † , φ δα = φ R δα , and denote by Ω † = φ † ( B ) and Ω δα = φ δα ( B ) the correspondingdomains parametrized by the radius functions R † and R δα , respectively.The velocity reconstruction from noisy velocity data u (cid:15) can then be formulatedcompactly as a linear inverse problem over the reference domain, i.e., T v = u (cid:15) ◦ φ δα in B, (12)with operator T : H ( B ) ∩ H ( B ) → L ( B ) defined by T ( v ) = v , i.e., as simpleembedding between Sobolev spaces. We assume that a bound on the data perturbations (cid:107) u † − u (cid:15) (cid:107) L ( D ) ≤ (cid:15) (13)is available, where u † denotes the true velocity field to be reconstructed, and we requirethat u † is piecewise smooth and vanishes outside Ω † , i.e., u † ∈ H ( D ) , u † | Ω † ∈ H (Ω † ) , and u † | D \ Ω † ≡ . (14)Let us note that these assumptions are reasonable, if the flow domain Ω † is smooth. Remark 3.3.
Observe that information about zero velocity at the boundary of theflow domain is encoded explicitly into the definition of the operator T . The particularformulation (12) over the reference domain B is used for the following reason: Inexactknowledge about the flow domain is shifted to the data u (cid:15) ◦ φ δα and the operator T ,therefore, is independent of the geometry reconstruction. This simplifies the analysisgiven in the following. Let us note that invertibility of the transformation φ δα isguaranteed by Lemma A.4., which allows to associate to any approximate solution v ofequation (12) a velocity field u = v ◦ ( φ δα ) − in the physical domain.For the stable solution of the nonlinear inverse problem (12), we again considerTikhonov regularization and we define the regularized approximations by v δ,(cid:15)α,β = arg min v ∈ D ( T ) (cid:107) T v − u (cid:15) ◦ φ δα (cid:107) L ( B ) + β (cid:107) ∆ v (cid:107) L ( B ) , (15)where D ( T ) = H ( B ) ∩ H ( B ) by definition. For our further analysis, we will requiretwo auxiliary results that we will present next. As a first step, we derive an estimatefor the data error in the reference domain. Lemma 3.4.
Let the assumptions (11), (13) and (14) hold. Then (cid:107) u † ◦ φ † − u (cid:15) ◦ φ δα (cid:107) L ( B ) ≤ δ U (16) with δ U = C ( δ / R (cid:107) u † (cid:107) H (Ω † ) + (cid:15) ) and a constant C that can be computed explicitly. Proof. By the integral transformation theorem we have (cid:107) u † ◦ φ † − u (cid:15) ◦ φ δα (cid:107) L ( B ) ≤ (cid:107) det( J δα ) − (cid:107) L ∞ ( B ) (cid:107) u † ◦ φ † ◦ ( φ δα ) − − u (cid:15) (cid:107) L (Ω δα ) ≤ C ( (cid:107) u † ◦ φ † ◦ ( φ δα ) − − u † (cid:107) L (Ω δα ) + (cid:107) u † − u (cid:15) (cid:107) L (Ω δα ) ) , n stable estimation of flow geometry and wall shear stress J δα of the transformation φ δα inthe second step. The last term can then be readily estimated by the bound (13) on thedata noise. The first term on the right hand side can further be split as (cid:107) u † ◦ φ † ◦ ( φ δα ) − − u † (cid:107) L (Ω δα ) ≤ (cid:107) u † ◦ φ † ◦ ( φ δα ) − − u † (cid:107) L (Ω δα ∩ Ω † ) + (cid:107) u † ◦ φ † ◦ ( φ δα ) − (cid:107) L (Ω δα \ Ω † ) = ( i ) + ( ii ) , where we used u † ≡ δα \ Ω † which follows from assumption (14). From thesmoothness of R δα , R † , and u † , and using the bound on the difference of the inversetransformations ( φ δα ) − and ( φ † ) − provided by Lemma A.5., we then obtain( i ) = (cid:107) u † ◦ φ † ◦ ( φ δα ) − − u † ◦ φ † ◦ ( φ † ) − (cid:107) L (Ω δα ∩ Ω † ) ≤ (cid:107) u † ◦ φ † (cid:107) W , ∞ ( B ) (cid:107) ( φ δα ) − − ( φ † ) − (cid:107) L ∞ (Ω δα ∩ Ω † ) | Ω δα ∩ Ω † | / ≤ C (cid:107) φ † (cid:107) W , ∞ ( B ) (cid:107) u † (cid:107) H (Ω † ) δ R ≤ C (cid:107) u † (cid:107) H (Ω † ) δ R . Observe that (cid:107) φ † (cid:107) W , ∞ ( B ) ≤ C (cid:107) R † (cid:107) W , ∞ (0 , π ) ≤ C (cid:48) (cid:107) R † (cid:107) H (0 , π ) by definition (10) of thetransformation and continuous embedding. Using Lemma A.1. and assumption (11),we can control the area | Ω δα \ Ω † | of the geometry mismatch by δ R , which allows us tobound the remaining term in the above estimate by (cid:107) u † ◦ φ † ◦ ( φ δα ) − (cid:107) L (Ω δα \ Ω † ) ≤ (cid:107) u † (cid:107) L ∞ ( D ) | Ω δα \ Ω † | / ≤ C (cid:107) u † (cid:107) H (Ω † ) δ / R . Note that the constants C in the above estimates are generic and independent of (cid:15) and δ . A combination of the bounds then yields the assertion of the lemma. (cid:3) Remark 3.5.
Following Remark 2.2. and the arguments used in the proof above, thebound δ U in (16) is computable in terms of the data noise levels δ and (cid:15) in (4) and (13),if bounds on (cid:107) R † (cid:107) H (0 , π ) and (cid:107) u † (cid:107) H (Ω † ) 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 ofthe operator T , which will allow us to utilized simple source conditions below. Lemma 3.6.
Let assumptions (11) and (14) hold and define v † := u † ◦ φ † .Then v † ∈ R (( T ∗ T ) µ ) for all µ < / . However, v † (cid:54)∈ R (( T ∗ T ) / ) , in general. Proof. We equip D ( T ) = H ( B ) ∩ H ( B ) with the norm (cid:107) v (cid:107) := (cid:107) ∆ v (cid:107) L ( B ) . Then forarbitrary v ∈ H ( B ) ∩ H ( B ) and f ∈ L ( B ), we have( T ∗ f, v ) H ( B ) ∩ H ( B ) = (∆ T ∗ f, ∆ v ) L ( B ) = ( f, T v ) L ( B ) = ( f, v ) L ( B ) . Thus w = T ∗ f is given as the unique solution of the boundary value problem∆ v = f in B with v = 0 and ∆ v = 0 on ∂B. From standard elliptic regularity theory [6], we can conclude that w = T ∗ f ∈ H ( B )for any f ∈ L ( B ). Using R (( T ∗ T ) / ) = R ( T ∗ ), see [5, Thm. 2.6], we thus arrive at R (( T ∗ T ) ) = { v ∈ H ( B ) | v = 0 on ∂B } , R (( T ∗ T ) / ) = { v ∈ H ( B ) | v = ∆ v = 0 on ∂B } . n stable estimation of flow geometry and wall shear stress u † and R † , we deduce that v † ∈ H ( B ) with v † = 0on ∂B , but ∆ v † (cid:54) = 0 on ∂B , in general. By interpolation of Sobolev spaces [10], we thusdeduce that v † ∈ R (( T ∗ T ) µ ) for µ < /
8, but not for µ ≥ /
8, in general. (cid:3)
Remark 3.7.
A range condition v † ∈ R (( T ∗ T ) µ ) would hold with µ = 1 /
4, if ∆ v † = 0at ∂B would be valid additionally; this can however not be expected in general. Thelimiting factor for the regularity index µ in the range condition, therefore, is themismatch of the higher order boundary conditions. This could be circumvented bychoosing a different equivalent norm on H ( B ) ∩ H ( B ); see [14] for details.From standard results about Tikhonov regularization in Hilbert spaces for linearinverse problems [5], we can now immediately conclude the following results. Theorem 3.8.
Let assumptions (11)–(14) hold and let v δ,(cid:15)α,β be defined by (15) withregularization parameter β = δ / (2 µ +1) U . Then (cid:107) v δ,(cid:15)α,β − v † (cid:107) H ( B ) ≤ Cδ µ µ +1 U for all ≤ µ < / . (17) The same estimates hold for a-posteriori parameter choice by the discrepancy principle β = max { β − n : n ≥ and such that (cid:107) v δ,(cid:15)α,β − v † (cid:107) ≤ δ U } . (18)Corresponding bounds for the error in the velocity u can be obtained, in principle, byback transformation into physical domain and some elementary computations. Let usclose 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 (cid:107) T v − u (cid:15) ◦ φ δα (cid:107) L ( B ) + β (cid:107) ∆ t v (cid:107) L ( B ) , with T : D ( T ) ⊂ H t ( B ) ∩ H ( B ) → L ( B ) defined by T v = v and t ≥ µ/ (2 µ + 1) sufficiently close to one.Further note that the data residual could also be measured in physical space. Theregularized approximate solution is then defined by (cid:101) v δ,(cid:15)α,β = arg min v ∈ D ( (cid:101) T ) (cid:107) (cid:101) T v − u (cid:15) (cid:107) L (Ω δα ) + β (cid:107) ∆ t v (cid:107) L ( B ) , with operator (cid:101) T : H ( B ) ∩ H ( B ) → L (Ω δα ) given by (cid:101) T v = v ◦ ( φ δα ) − . This simplyamounts to a change to an equivalent norm in the data space. A quick inspection of thearguments in the above proof reveals that the assertions of Theorem 3.8. remain validalso for this choice of regularization method, which is more easy to implement and willthus be used in our numerical tests in Section 6.
4. 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, weassume that fluid viscosity is normalized, and define the associated wall shear stress by τ R ( v )( ϕ ) = − n R ( ϕ ) · J R (cos ϕ, sin ϕ ) − · ( ∇ v )(cos ϕ, sin ϕ ) , (19) n stable estimation of flow geometry and wall shear stress J R is the Jacobian of φ R , and n R ( ϕ ) is the outward pointing unit normal vectorat the corresponding point ( R ( ϕ ) cos ϕ, R ( ϕ ) sin ϕ ) ∈ ∂ Ω R in the physical domain. Letus note that n R ( ϕ ) can be expressed explicitly by n R ( ϕ ) = ( R ( ϕ ) cos ϕ + R (cid:48) ( ϕ ) sin ϕ, R ( ϕ ) sin ϕ − R (cid:48) ( ϕ ) cos ϕ ) (cid:112) ( R ( ϕ )) + ( R (cid:48) ( ϕ )) . (20)For ease of notation, we will identify ∂B with the interval (0 , π ) in the sequel. Inaddition, we again define symbols τ δ,(cid:15)α,β = τ R δα ( v δ,(cid:15)α,β ) and τ † = τ R † ( v † ) where v † = u † ◦ φ † .A combination of the results derived so far and some elementary geometric computationsnow leads to the following assertion. Theorem 4.10.
Let the assumptions of Theorem 3.7 hold. Then (cid:107) τ δ,(cid:15)α,β − τ † (cid:107) L (0 , π ) ≤ C ( δ R + δ (2 µ ) / (2 µ +1) U ) for all ≤ µ < / . (21)Proof. We use the definition of τ R ( v ) and decompose the error into the three parts (cid:107) τ δ,(cid:15)α,β − τ † (cid:107) L (0 , π ) = (cid:107) n R δα J − R δα ∇ v δ,(cid:15)α,β − n R † J − R † ∇ v † (cid:107) L ( ∂B ) ≤ (cid:107) ( n R δα − n R † ) J − R δα ∇ v δ,(cid:15)α,β (cid:107) L ( ∂B )) + (cid:107) n R † ( J − R δα − J − R † ) ∇ v δ,(cid:15)α,β (cid:107) L ( ∂B )) + (cid:107) n R † J − R † ( ∇ v δ,(cid:15)α,β − ∇ v † ) (cid:107) L ( ∂B )) = ( i ) + ( ii ) + ( iii ) . Using Lemmas A.2. and A.3., and the embedding of Sobolev spaces, we obtain( i ) ≤ (cid:107) n R δα − n R † (cid:107) L ∞ ( ∂B )) (cid:107) J − R δα (cid:107) L ∞ ( ∂B )) (cid:107)∇ v δ,(cid:15)α,β (cid:107) L ( ∂B )) ≤ C (cid:107) R δα − R † (cid:107) H (0 , π ) (cid:107) R δα (cid:107) H (0 , π ) (cid:107) v δ,(cid:15)α,β (cid:107) H ( B ) ≤ Cδ R . By the geometric arguments of Lemma A.5., the second term can be estimated by( ii ) ≤ (cid:107) n R † (cid:107) L ∞ ( ∂B )) (cid:107) J − R δα − J − R † (cid:107) L ∞ ( ∂B )) (cid:107)∇ v δ,(cid:15)α,β (cid:107) L ( ∂B )) ≤ C (cid:107) R δα − R † (cid:107) H (0 , π ) (cid:107) v δ,(cid:15)α,β (cid:107) H ( B ) ≤ Cδ R . With the help of the results of the previous section, the third term, which measures theamplification of the velocity approximation error, can be bounded by( iii ) ≤ (cid:107) n R † (cid:107) L ∞ ( ∂B )) (cid:107) J − R † (cid:107) L ∞ ( ∂B )) (cid:107)∇ v δ,(cid:15)α,β − ∇ v † (cid:107) L ( ∂B )) ≤ C (cid:107) R † (cid:107) H (0 , π ) (cid:107) v δ,(cid:15)α,β − v † (cid:107) H ( B ) ≤ Cδ (2 µ ) / (2 µ +1) U . The assertion of the theorem then follows by combination of these estimates. (cid:3)
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 trueflow geometry and velocity are sufficiently smooth and the regularization terms in theTikhonov functionals (5) and (15) are chosen sufficiently strong. The main observationof the previous theorem therefore is, that it is possible to obtain a quantitative estimateunder reasonable smoothness assumptions. n stable estimation of flow geometry and wall shear stress
5. Remarks on the extension to three dimensions
For a three-dimensional flow, the wall shear stress is a tensor defined by τ = − µ (cid:0) ∇ u + ( ∇ u ) T (cid:1) · n, where µ is the dynamic viscosity, u is the velocity vector, and n the outer unit normalvector at the vessel wall. If appropriate measurements of the density and of all threevelocity components are available, then the reconstruction approach and the theoreticalresults presented in the previous sections can be generalized almost verbatim to the threedimensional setting; only the computational realization becomes more complicated.Note that the geometry and velocity reconstruction naturally decompose into a sequenceof 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 z -coordinate, one has to employ regularization also in the z -direction, which fully couplesthe problems for the individual cross-sections. The additional computational complexitydue to this coupling can be overcome by a Kaczmarz strategy [8, 11], which allows toreduce the numerical solution to the iterated solution of two-dimensional problems forthe 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.
6. Numerical validation
In order to demonstrate the viability of the proposed approach, we now report about thereconstruction of flow geometry, flow velocity, and wall shear stress from experimentaldata obtained in a clinical whole-body 3 Tesla magnetic resonance imaging scanner(Prisma, Siemens Healthcare, Erlangen) at the University Medical Center Freiburg. As atest case, we consider the flow of water through a cylindrical pipe with constant diameter d = 25 . ◦ C;details about the experimental setup are described in [1]. Let us note that the flowconditions lead to a Reynolds number of about 5300 and thus amount to a turbulent flowwith steep velocity gradients in the boundary layer. For our reconstruction procedures,we utilize magnetic resonance images of density and the axial flow velocity acquired withdifferent resolutions corresponding to in plane voxel sizes of h = 1 . h = 0 . † = φ † ( B ) of the pipe is knownhere. A reference solution u ref for the flow velocity can thus be computed by numericalsimulation [9] and will be used for comparison with our results. Let us emphasizethat this reference solution represents a time averaged velocity field, in which alltemporal 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. Onetherefore cannot expect to get a perfect fit to the simulated reference velocity data. n stable estimation of flow geometry and wall shear stress For the geometry identification, we utilize the standard magnetic resonance images ofthe density. After a simple scaling procedure, see Appendix B for details, the data canbe interpreted as m δ | V i = 1 | V i | (cid:90) V i χ Ω † ( x ) dx + noise , where V i denotes the i th voxel of the measurement array of size h × h . The action of theforward operator for our computational tests is then defined as F γ,h ( R ) | V i = 1 | V i | (cid:90) V i H γ ( R ( ϕ ( ξ )) − | ξ | ) dξ. (22)where H γ ( x ) = π arctan (cid:16) xγ (cid:17) + , γ >
0, is a smooth approximation of the Heavisidefunction. In our computational tests, we thus minimize the Tikhonov functional J δα,γ ( R ) = (cid:107) F γ,h ( R ) − m δ (cid:107) L ( D ) + α (cid:107) R (cid:107) H (0 , π ) over the finite dimensional space of radius functions V N = { R N ∈ H per (0 , π ) | R N ( ϕ ) = b + N (cid:88) k =1 a k sin( kϕ ) + b k cos( kϕ ) } ∩ D ( F ) . Due to the smoothing with γ >
0, the forward operator F h,γ is continuously differentiableand 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 asthe barycenter of the measurement data m δ and the true radius was R † = 12 . (cid:107) R δα − R † (cid:107) H (0 , π ) / (cid:107) R † (cid:107) H (0 , π ) obtained formeasurement acquired with different data resolution h , and for different choices of theregularization parameter α . h \ α Table 1.
Relative reconstruction errors (cid:107) R δα − R † (cid:107) H (0 , π ) / (cid:107) R † (cid:107) H (0 , π ) for differentdata resolutions h and regularization parameters α ; optimal results in bold. As expected intuitively, reconstruction errors are slightly decreasing with increasingdata resolution. Further note that the reconstructions are stable with respect to thechoice 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. n stable estimation of flow geometry and wall shear stress R † and R δα obtained in ournumerical tests. The maximal error in the reconstructed radius is about 0 . Figure 1.
Density measurements ( h = 1mm) and reconstructed geometry (left);corresponding radius functions R † , R δα (right); axis labeling in mm. is substantially less than the voxel size h = 1mm of the data. This illustrates that sub-pixel resolution can be obtained by the proposed geometry reconstruction procedure. We now turn to the reconstruction of the velocity field, for which we utilize phase-contrast 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 parameteraccording 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 u (cid:15) | V i = 1 | V i ∩ Ω δα | (cid:90) V i ∩ Ω δα u † ( x ) dx + noise . The remaining voxels only contain information about the noise and will not be used inthe reconstruction; see Figure 2 and the remarks given in Appendix B.Following Remark 3.9., we define the corresponding forward operator mapping tothe physical domain by (cid:101) T h : H ( B ) ∩ H ( B ) → L (Ω δα ) with (cid:101) T h ( v ) | V i = 1 | V i ∩ Ω δα | (cid:90) V i ∩ Ω δα v ◦ ( φ δα ) − ( x ) dx. Note that only voxels V i with | V i ∩ Ω δα | > (cid:101) T h . Theregularized approximation for the velocity field in the reference domain is then definedas the minimizer of the Tikhonov functional (cid:101) v δ,(cid:15)α,β = arg min v ∈ V M (cid:107) (cid:101) T h v − u (cid:15) (cid:107) L (Ω δα ) + β (cid:107) ∆ v (cid:107) L ( B ) , n stable estimation of flow geometry and wall shear stress Figure 2.
Velocity raw data u (cid:15) with large noise outside the flow domain (left) andtruncated velocity raw data actually used in the reconstruction (right). over the space V M = span { v ∈ H ( B ) | − ∆ v = λv for some λ ≤ M } of eigen functionsof the Dirichlet Laplace operator on the unit disk. The solution (cid:101) v δ,(cid:15)α,β can be computedefficiently by Cholesky factorization or the conjugate gradient method. In Table 2,we display the reconstruction errors obtained with our algorithm for different dataresolutions h and regularization parameters β . h \ β . · − . · − . · − . · − . · − Table 2.
Relative reconstruction errors (cid:107) (cid:101) v δ,(cid:15)α,β − v ref (cid:107) H ( B ) / (cid:107) v ref (cid:107) H ( B ) , with v ref = u ref ◦ φ † denoting the simulated velocity on the reference domain, for different dataresolutions h and regularization parameters β ; optimal results in bold. Recall that u ref , which was obtained by numerical simulation, corresponds to atime averaged velocity field and therefore does not contain temporal fluctuations thatare present in the measurements; see Figure 3. This explains the relatively large errorsin Table 2 and their independence of the data resolution. The optimal regularizationparameters β 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 inFigure 3 shows that the reconstructed velocity field correctly reproduces some turbulentfluctuations which, due to time averaging, are not present in the simulated referencevelocity field. Apart from these differences, the reconstruction is in good agreementwith the reference solution. As a final step in our numerical tests, we now utilize the reconstructed radius function R δα and velocity fields v δ,(cid:15)α,β to compute the approximation τ δ,(cid:15)α,β for the wall shear stress n stable estimation of flow geometry and wall shear stress Figure 3.
Reference velocity v ref = u ref ◦ φ † (left) obtained by simulation andreconstructed velocity (cid:101) v δ,(cid:15)α,β (right) for data resolution h = 1mm and β = 8 . · − . by (19). In Table 3, we display the results obtained for optimal α and β = 1 . · − .Let us note that most part of the error stems from perturbations in higher modes, whichh 1.00 0.75 0.60 0.43 0.30 (cid:107) τ δ,(cid:15)α,β − τ † (cid:107) L (0 , π ) Table 3.
Relative errors (cid:107) τ δ,(cid:15)α,β − τ † (cid:107) L (0 , π ) in the reconstruction of the wall shearstress for different data resolutions resp. voxel size h . can be suppressed efficiently by application of a low-pass filter. In Figure 4, we plot thereconstruction of the wall shear stress τ δ,(cid:15)α,β and its constant approximation τ δ,(cid:15)α,β againstthe reference value τ ref obtained from numerical flow simulation [9]. Figure 4.
Reconstructed wall shear stress τ δ,(cid:15)α,β and constant approximation τ δ,(cid:15)α,β (red)in comparison to the reference wall shear stress τ ref (blue) obtained by simulation.All results are functions of angle ϕ with values in Pa. Note that the average wall shear stress τ δ,(cid:15)α,β is in very good agreement with thereference value τ ref obtained for the simulated data. The local variations in thereconstruction τ δ,(cid:15)α,β can be explained by the turbulent variations of velocity in the data;compare with the plots in Figure 2 and 3. n stable estimation of flow geometry and wall shear stress
7. Discussion
As reviewed by Petersson [16], the stable and accurate estimation of wall shear stressfrom magnetic resonance imaging data, is a delicate issue and most reconstructionapproaches reported in literature seem not to work properly. In this paper, we thereforeconsidered a systematic approach for the estimation of wall shear stress from magneticresonance images of density and velocity, for which stability and convergence could beestablished under simple and realistic smoothness assumptions on the flow geometry andvelocity. The theoretical results were validated by numerical tests for experimental datawhich demonstrate that wall shear stress can be estimated from magnetic resonanceimaging data with relative errors of a few percent and practically independent of thedata 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 empiricalMoody charts [13] or the Clauser plot method [2, 18], which are, however, limitedto axisymmetric geometry or fully developed turbulent flow. Numerical simulationscould also be used, in principle, to compute wall shear stress estimates [9], but preciseknowledge about the rheological properties of the fluid are required. In contrast to theseapproaches, the method considered in this paper is generally applicable and, therefore,seems most appropriate for application in a clinical context.
Acknowledgements
Experimental data were acquired at the the University Medical Center, Freiburg,together with A. Bauer (SLA, TU Darmstadt) and A. Krafft, N. Shokina (Med.Phys.,UMC Freiburg). Funding of the authors by the German Research Foundation (DFG)via grant Eg-331/1-1 and through grant GSC 233 of the “Excellence Initiative” of theGerman Federal and State Governments is gratefully acknowledged.
Appendix A. Geometric results
In the following, we present some auxiliary results concerning geometrical details. Recallthat D ( F ) = { R ∈ H per (0 , π ) | r < R < r } for some 0 < r < r < R = φ R ( B ) is the image of the unit ball B under the transformation φ R ( r cos( ϕ ) , r sin( ϕ )) = ( r + ( R ( ϕ ) − r ) r η ) (cos( ϕ ) , sin( ϕ ))with η ≥
2. The Jacobian matrix of φ R is denoted by J R , and for given functions R , R ,the subscripts transfer to the associated objects, e.g. φ = φ R and J = J R .As a first result, we estimate the differences of domains Ω R = φ R ( B ) in terms ofdifferences of their parameterizing radius functions. Lemma A.1.
Let R , R ∈ D ( F ) . Then | Ω \ Ω | ≤ C (cid:107) R − R (cid:107) H (0 , π ) . (A.1) n stable estimation of flow geometry and wall shear stress R max ( ϕ ) = max { R ( ϕ ) , R ( ϕ ) } . Then | Ω \ Ω | = (cid:90) π (cid:90) R max ( ϕ ) R ( ϕ ) r dr dϕ = 12 (cid:90) π ( R max ( ϕ ) + R ( ϕ ))( R max ( ϕ ) − R ( ϕ )) dϕ ≤ πr (cid:107) R max − R (cid:107) L ∞ (0 , π ) ≤ πr C (cid:107) R − R (cid:107) H (0 , π ) . In the last step, we used the continuous embedding H (0 , π ) (cid:44) → L ∞ (0 , π ). (cid:3) As a next step, we estimate differences in the normal vector n R defined in (20), interms of differences in the radius function. Lemma A.2.
Let R , R ∈ D ( F ) . Then (cid:107) n − n (cid:107) L ∞ ( ∂B ) ≤ C (cid:107) R − R (cid:107) H (0 , π ) . (A.2)Proof. Let us introduce e r = (cos( ϕ ) , sin( ϕ )) and e ϕ = ( − sin( ϕ ) , cos( ϕ )). Omitting theexplicit notion of the dependence on ϕ , we then obtain for any angle ϕ that | n − n | = | λ r e r + λ ϕ e ϕ | = λ r + λ ϕ . The radial component of the difference can be estimated by | λ r | = (cid:12)(cid:12)(cid:12)(cid:12) R √ R + R (cid:48) − R √ R + R (cid:48) (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:12)(cid:12)(cid:12)(cid:12) R − R √ R + R (cid:48) + R √ R + R (cid:48) − √ R + R (cid:48) √ R + R (cid:48) √ R + R (cid:48) (cid:12)(cid:12)(cid:12)(cid:12) ≤ | R − R | r + R √ ( R − R ) +( R (cid:48) − R (cid:48) ) √ R + R (cid:48) √ R + R (cid:48) ≤ √ r (cid:107) R − R (cid:107) W , ∞ (0 , π ) In a similar way, one can estimate λ ϕ , and by Sobolev’s embedding theorem, we obtain (cid:107) n − n (cid:107) L ∞ ( ∂B ) ≤ (cid:107) λ r (cid:107) L ∞ ( ∂B ) + (cid:107) λ ϕ (cid:107) L ∞ ( ∂B ) ≤ C (cid:107) R − R (cid:107) H (0 , π ) , which already yields the desired estimate. (cid:3) The following result ensures smoothness of the transformations φ R whenever theradius function R is sufficiently smooth. Lemma A.3.
Let R , R ∈ D ( F ) ∩ H kper (0 , π ) for k ≤ η with η ≥ as in (10). Then (cid:107) φ (cid:107) W k − , ∞ ( B ) ≤ C (cid:107) R (cid:107) H k (0 , π ) . (A.3) where φ = φ R with φ R defined in (10). Moreover, (cid:107) φ − φ (cid:107) W k − , ∞ ( B ) ≤ C (cid:107) R − R (cid:107) H k (0 , π ) (A.4)Proof. The continuity of φ is obvious and the Jacobian J = J R obtained by derivationwith respect to coordinates x = ( r cos φ, r sin φ ) reads J ( r cos ϕ, r sin ϕ ) = (cid:32) r + η ( R ( ϕ ) − r ) r η − R (cid:48) ( ϕ ) r η − r + ( R ( ϕ ) − r ) r η − (cid:33) . (A.5)Since η ≥
2, the Jacobian can be seen to be continuous also in r = 0, and we have (cid:107) J (cid:107) L ∞ ( B ) ≤ C ( (cid:107) R (cid:107) L ∞ (0 , π ) + (cid:107) R (cid:48) (cid:107) L ∞ (0 , π ) ) ≤ C (cid:107) R (cid:107) H (0 , π ) . This shows the first estimate for k = 1. The assertion for higher order derivatives of φ follow in a similar way. From the formula (10), one can see that φ R is affine linear in R . n stable estimation of flow geometry and wall shear stress (cid:3) As a next step, we show that the transformations φ R defined in (10) are invertible. Lemma A.4.
Let R ∈ D ( F ) . Then the transformation φ R defined in (10) is adiffeomorphism with inverse transformation ψ R = ( φ R ) − and (cid:107) J − R (cid:107) L ∞ ( B ) = (cid:107) J ψ R (cid:107) L ∞ (Ω R ) ≤ C (cid:107) R (cid:107) H (0 , π ) (A.6)Proof. Using (A.5), we can estimate the determinant bydet( J R ( r cos ϕ, r sin ϕ )) = ( r + η ( R ( ϕ ) − r ) r η − ) · ( r + ( R ( ϕ ) − r ) r η − ) ≥ r . The existence of an inverse transformation ψ R = ( φ R ) − then follows form the implicitfunction theorem and the Jacobian of the inverse mapping is J ψ R ◦ φ R = J − R . Thebounds for J ψ R can then be deduced in an elementary way. (cid:3) Using the previous results, we can also bound differences in the inverse mappings.
Lemma A.5.
Let R , R ∈ D ( F ) and assume that R ∈ H per (0 , π ) . Furthermore, let ψ , ψ denote the corresponding inverse transformations with Jacobians J ψ , J ψ . Then (cid:107) ψ − ψ (cid:107) L ∞ (Ω ∩ Ω ) ≤ C (cid:107) R − R (cid:107) H (0 , π ) (A.7) and the difference in the Jacobians can be bounded by (cid:107) J ψ − J ψ (cid:107) L ∞ (Ω ∩ Ω ) ≤ C (cid:107) R − R (cid:107) H (0 , π ) (A.8)Proof. Step 1:
To show (A.7), let x ∈ Ω ∩ Ω . Then there are x , x ∈ B with φ ( x ) = φ ( x ) = x . Since the transformations φ and φ preserve angles, the angularparts ϕ ( x ) = ϕ ( x ) = ϕ ( x ) =: ϕ are equal. φ has an inverse ψ , hence x − x = ψ ( φ ( x )) − ψ ( φ ( x )) = ψ ( x ) − ψ ( x + dx ) , where the defect is given by dx = φ ( x ) − φ ( x ) = ( R ( ϕ ) − R ( ϕ )) | x | η (cos( ϕ ) , sin( ϕ )) . Since x and dx have the same angular coordinate and Ω is star shaped with respect tothe origin, we have { x + t dx | t ∈ [0 , } ⊂ Ω . The mean value theorem yields | ψ ( x ) − ψ ( x ) | = | x − x | = | J ψ ( ξ ) dx | ≤ (cid:107) J ψ (cid:107) W , ∞ (Ω ) (cid:107) R − R (cid:107) H (0 , π ) . The assertion (A.7) follows by the estimate of Lemma A.4.
Step 2: Show (A.8).
Startingat (A.5) elementary calculus yields (cid:107) det J φ − det J φ (cid:107) L ∞ ( B ) ≤ C (cid:107) R − R (cid:107) L ∞ (0 , π ) , where the constant C depends only on (cid:107) R (cid:107) L ∞ (0 , π ) , (cid:107) R (cid:107) L ∞ (0 , π ) , r and η . For theinverse Jacobians, we further conclude that (cid:107) J − φ − J − φ (cid:107) L ∞ ( B ) ≤ (cid:13)(cid:13)(cid:13) J φ J φ − J φ J φ (cid:13)(cid:13)(cid:13) L ∞ ( B ) ≤ (cid:13)(cid:13)(cid:13) det J φ − det J φ det J φ det J φ J φ (cid:13)(cid:13)(cid:13) L ∞ ( B ) + (cid:13)(cid:13)(cid:13) J φ ( J φ − J φ ) (cid:13)(cid:13)(cid:13) L ∞ ( B ) ≤ C (cid:16) r (cid:107) R (cid:107) H (0 , π ) + r (cid:17) (cid:107) R − R (cid:107) H (0 , π ) . n stable estimation of flow geometry and wall shear stress φ has the representation J − φ = ˜ J φ / det J φ , where ˜ J φ isa rearrangement of J φ and all expressions are continuously differentiable, one can verifywithout difficulty that J − φ ∈ W , ∞ ( B ) with the associated norm bounded in terms of η , r , and (cid:107) R (cid:107) H (0 , π ) . Therefore we arrive at (A.8) by (cid:107) J ψ − J ψ (cid:107) L ∞ (Ω ∩ Ω ) = (cid:107) J − φ ◦ ψ − J − φ ◦ ψ (cid:107) L ∞ (Ω ∩ Ω ) ≤ (cid:107) J − φ (cid:107) W , ∞ ( B ) (cid:107) ψ − ψ (cid:107) L ∞ (Ω ∩ Ω ) + (cid:107) J − φ − J − φ (cid:107) L ∞ ( B ) ≤ C (cid:107) R − R (cid:107) H (0 , π ) . This completes the proof of the second estimate of the lemma. (cid:3)
Appendix B. Interpretation of the data
Let us briefly comment on the physical interpretation and the preprocessing of theexperimental 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 overvoxels; this values are additionally perturbed by data noise. From a histogram of themagnitude data, see Figure B1, one can deduce to peak values m and m , which areused for scaling of the data. Figure B1.
Magnitude raw data (left) and magnitude value histogram (right) usedto normalize the magnitude data.
Based on the two values m , m , the normalized magnitude m δ are then defined by m δ | V i = T (cid:18) m δraw | V i − m m − m (cid:19) , with truncation function T ( x ) = max(0 , min(1 , x )). These measurements then representan approximation for the characteristic function of the flow domain. n stable estimation of flow geometry and wall shear stress Appendix B.2. Velocity data
The raw data acquired in phase contrast magnetic resonance velocimetry can beinterpreted as d δ | V i = (cid:90) V i ρ ( x ) e i πu ( x ) /v enc + noise , where ρ is the propton density and v enc the velocity encoding parameter. An averagevalue of the velocity in the voxel V i is then recovered as the phase of this signal, i.e., u δ | V i = v enc · arg( d δ | V i ) . (B.1)Phase unwrapping may be required, if the maximal velocity is larger then v enc . In theabsence of data noise and assuming that the density is given by ρ = cχ Ω , one can useTaylor approximation to deduce that u δ | V i ≈ | V i ∩ Ω | (cid:90) V i ∩ Ω u ( x ) dx, where Ω is the flow domain and u the true flow velocity. This is the measurementmodel used in the numerical tests. Let us note that the phase retrieval in (B.1) isparticularly sensitive to data perturbations, if abs( d δ | V i ) is small, which is the case closeto the boundary and outside the flow domain. Therefore, particularly large noise in thevelocity data is expected for these voxels; compare with Figure 2. Bibliography [1] A. Bauer, S. Wegt, C. Tropea, A. Krafft, N. Shokina, J. Hennig, G. Teschner, and H. Egger.Ground-truth for measuring wall shear stress with magnetic resonance velocimetry. In .[2] F. H. Clauser. The turbulent boundary layer.
Advances in Applied Mechanics , 4:1–51, 1956.[3] H. Egger and B. Hofmann. Tikhonov regularization in Hilbert scales under conditional stabilityassumptions.
Inverse Problems , 34:115015, 2018.[4] H. Egger and A. Leit˜ao. Nonlinear regularization methods for ill-posed problems with piecewiseconstant or strongly varying solutions.
Inverse Problems , 25:115014, 19, 2009.[5] H. W. Engl, M. Hanke, and A. Neubauer.
Regularization of inverse problems , volume 375 of
Mathematics and its Applications . Kluwer Academic Publishers Group, Dordrecht, 1996.[6] D. Gilbarg and N. S. Trudinger.
Elliptic Partial Differential Equations of Second Order , volume224 of
Grundlehren der mathematischen Wissenschaften . Springer, Berlin, 2001.[7] M. A. Gimbrone, T. Nagel, and J. N. Topper. Biomechanical activation: an emerging paradigmin endothelial adhesion biology.
J. Clinical Investigation , 99:1809–1813, 1997.[8] B. Kaltenbacher, A. Neubauer, and O. Scherzer.
Iterative regularization methods for nonlinear ill-posed problems , volume 6 of
Radon Series on Computational and Applied Mathematics . Walterde Gruyter GmbH & Co. KG, Berlin, 2008.[9] G. K. E. Khoury, P. Schlatter, A. Noorani, P. F. Fischer, G. Brethouwer, and A. V. Johansson.Direct numerical simulation of turbulent pipe flow at moderately high Reynolds numbers.
J.Flow Turbul. Combust. , 91:475–495, 2013.[10] A. Lunardi.
Interpolation Theory . Edizione della Normale, Pisa, 2009.[11] F. Margotti, A. Rieder, and A. Leit˜ao. A Kaczmarz version of the reginn -Landweber iterationfor ill-posed problems in Banach spaces.
SIAM J. Numer. Anal. , 52:1439–1465, 2014. n stable estimation of flow geometry and wall shear stress [12] M. Markl, A. Frydrychowicz, S. Kozerke, M. Hope, and O. Wieben. 4D flow MRI. Journal ofMagnetic Resonance Imaging , 36:1015–1036, 2012.[13] L. F. Moody. Friction Factors for Pipe Flow.
Trans. ASME , 66:671–684, 1944.[14] A. Neubauer. When do Sobolev spaces form a Hilbert scale?
Proc. Amer. Math. Soc. , 103:557–562, 1988.[15] S. Oyre, S. Ringgaard, S. Kozerke, W. P. Paaske, M. B. Scheidegger, P. Boesiger, and E. M.Pedersen. Quantitation of circumferential subpixel vessel wall position and wall shear stress bymultiple sectored three-dimensional paraboloid modeling of velocity encoded Cine MR.
Magn.Res. Med. , 40:645–655, 1998.[16] S. Petersson, P. Dyverfeldt, and T. Ebbers. Assessment of the accuracy of MRI wall shear stressestimation using numerical simulations.
Magn. Res. Imag. , 36:128–138, 2012.[17] W. V. Potters, P. Ooij, H. Marquering, E. van Bavel, and A. J. Nederveen. Volumetric arterial wallshear stress calculation based on cine phase contrast MRI.
J. Magn. Res. Imag. , 41:505–516,2015.[18] N. Shokina, W. B. Buchenberg, M. Menza, A. Bauer, G. Teschner, C. Tropea, H. Egger, J. Hennig,and A. J. Krafft. Accurate MR-based wall shear stress measurements in fully developed turbulentflow using the Clauser-plot method. Abstract No 2946 in Joint annual meeting ISMRM-ESMRMB, Paris, France, 16-21 June 2018.[19] A. F. Stalder, M. F. Russe, A. Frydrychowicz, J. Bock, J. Hennig, and M. Markl. Quantitative2D and 3D phase contrast MRI: Optimized analysis of blood flow and vessel wall parameters.
Magn. Res. Imag. , 60:1218–1231, 2008.[20] U. Tautenhahn. On a general regularization scheme for nonlinear ill-posed problems. II.Regularization in Hilbert scales.