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

TL;DR

This paper introduces a systematic and stable method for reconstructing flow geometry and wall shear stress from MRI data, with error quantification and validation on experimental data.

Contribution

It proposes a new regularization-based reconstruction approach that improves stability and error estimation in flow geometry and wall shear stress from MRI measurements.

Findings

The method provides stable estimates under measurement noise.

Error bounds are derived based on smoothness assumptions.

Experimental data demonstrates the approach's accuracy and robustness.

Abstract

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.

Figures

Tables

Table 1. Table 1: Relative reconstruction errors ‖ R α δ − R † ‖ H 2 ​ ( 0 , 2 ​ π ) / ‖ R † ‖ H 2 ​ ( 0 , 2 ​ π ) subscript norm superscript subscript 𝑅 𝛼 𝛿 superscript 𝑅 † superscript 𝐻 2 0 2 𝜋 subscript norm superscript 𝑅 † superscript 𝐻 2 0 2 𝜋 \|R_{\alpha}^{\delta}-R^{\dagger}\|_{H^{2}(0,2\pi)}/\|R^{\dagger}\|_{H^{2}(0,2\pi)} for different data resolutions h ℎ h and regularization parameters α 𝛼 \alpha ; optimal results in bold.

$h\setminus \alpha$	0.16	0.08	0.04	0.02	0.01
1.00	0.0220	0.0117	0.0098	0.0128	0.0194
0.75	0.0161	0.0111	0.0118	0.0155	0.0215
0.60	0.0132	0.0086	0.0089	0.0129	0.0201
0.43	0.0104	0.0079	0.0090	0.0121	0.0172
0.30	0.0064	0.0053	0.0063	0.0087	0.0127

Table 2. Table 2: Relative reconstruction errors ‖ v ~ α , β δ , ϵ − v r ​ e ​ f ‖ H 2 ​ ( B ) / ‖ v r ​ e ​ f ‖ H 2 ​ ( B ) subscript norm superscript subscript ~ 𝑣 𝛼 𝛽 𝛿 italic-ϵ superscript 𝑣 𝑟 𝑒 𝑓 superscript 𝐻 2 𝐵 subscript norm superscript 𝑣 𝑟 𝑒 𝑓 superscript 𝐻 2 𝐵 \|\widetilde{v}_{\alpha,\beta}^{\delta,\epsilon}-v^{ref}\|_{H^{2}(B)}/\|v^{ref}\|_{H^{2}(B)} , with v r ​ e ​ f = u r ​ e ​ f ∘ ϕ † superscript 𝑣 𝑟 𝑒 𝑓 superscript 𝑢 𝑟 𝑒 𝑓 superscript italic-ϕ † v^{ref}=u^{ref}\circ\phi^{\dagger} denoting the simulated velocity on the reference domain, for different data resolutions h ℎ h and regularization parameters β 𝛽 \beta ; optimal results in bold.

$h\setminus \beta$	$3.2\cdot {10}^{-4}$	$1.6\cdot {10}^{-4}$	$8.0\cdot {10}^{-5}$	$4.0\cdot {10}^{-5}$	$2.0\cdot {10}^{-5}$
1.00	0.3815	0.2661	0.1818	0.2249	0.3985
0.75	0.3317	0.2233	0.1926	0.2912	0.4541
0.60	0.3377	0.2313	0.1813	0.2507	0.3899
0.43	0.3098	0.2011	0.1795	0.2851	0.4478
0.30	0.3245	0.2090	0.1518	0.2106	0.3190

Table 3. Table 3: Relative errors ‖ τ α , β δ , ϵ − τ † ‖ L 2 ​ ( 0 , 2 ​ π ) subscript norm superscript subscript 𝜏 𝛼 𝛽 𝛿 italic-ϵ superscript 𝜏 † superscript 𝐿 2 0 2 𝜋 \|\tau_{\alpha,\beta}^{\delta,\epsilon}-\tau^{\dagger}\|_{L^{2}(0,2\pi)} in the reconstruction of the wall shear stress for different data resolutions resp. voxel size h ℎ h .

h	1.00	0.75	0.60	0.43	0.30
${\Vert {\tau }_{\alpha ,\beta }^{\delta ,ϵ}-{\tau }^{†}\Vert }_{L^2(0,2\pi )}$	0.0789	0.0385	0.0391	0.0264	0.0405

Equations

$${\Omega }_R=\{(r\cos \phi ,r\sin \phi ):0\le r<R(\phi ),0\le \phi \le 2\pi \},$$

$$R^{†}\in H_{per}^k(0,2\pi )\text{mbox}with\text{mbox}{r}_0\le R(\phi )\le r_1$$

$$F(R)=m^{\delta }\text{mbox}onD,$$

$$\Vert {\chi }_{{\Omega }^{†}}-m^{\delta }{\Vert }_{L^2(D)}\le \delta$$

$$J_{\alpha }^{\delta }(R)=\Vert F(R)-m^{\delta }{\Vert }_{L^2(D)}^2+\alpha \Vert R{\Vert }_{H^2(0,2\pi )}^2.$$

$$J_{\alpha }^{\delta }(R_{\alpha }^{\delta })\le \underset{R\in D(F)}{\inf }J(R)+{\delta }^2.$$

$$\Vert F(R_{\alpha }^{\delta })-{\chi }_{{\Omega }^{†}}{\Vert }_{L^2(D)}\le C\delta \text{mbox}and\Vert R_{\alpha }^{\delta }-R^{†}{\Vert }_{H^r(0,2\pi )}\le C{\delta }^{1-r/k}$$

$$\alpha =\max \{{\alpha }_0{2}^{-n}:n\ge 0\text{mbox}suchthat\Vert F(R_{\alpha }^{\delta })-m^{\delta }{\Vert }_{L^2(D)}\le 4\delta \}.$$

$$|F(R)-F(\widetilde{R})|(r\cos\varphi,r\sin\varphi)=\left\{\begin{array}[]{ll}1,&\mbox{for }r_{min}(\varphi)\leq r\leq r_{max}(\varphi),\\ 0&\mbox{else}.\end{array}\right.$$

$$\Vert F(R)-F(\tilde{R}){\Vert }_{L^2(D)}^2={\int }_0^{2\pi }{\int }_0^{r_1}|F(R)-F(\tilde{R}){|}^{2}rdrd\phi ={\int }_0^{2\pi }{\int }_{r_{min}(\phi )}^{r_{max}(\phi )}rdrd\phi$$

$$=\frac{1}{2}{\int }_0^{2\pi }{r}_{max}(\phi {)}^2-r_{min}(\phi {)}^{2}d\phi =\frac{1}{2}{\int }_0^{2\pi }|R(\phi )-\tilde{R}(\phi )||R(\phi )+\tilde{R}(\phi )|d\phi$$

$$\ge \frac{1}{2}{\int }_0^{2\pi }|R(\phi )-\tilde{R}(\phi )|2r_0\frac{|R(\phi )-\tilde{R}(\phi )|}{|r_1-r_0|}d\phi =\frac{r_0}{r_1-r_0}\Vert R-\tilde{R}{\Vert }_{L^2(0,2\pi )}^2.$$

$$\Vert R_{\alpha }^{\delta }-R^{†}{\Vert }_{H^2(0,2\pi )}\le {\delta }_R,$$

$$\Vert R_{\alpha }^{\delta }-R^{†}{\Vert }_{H^r}\le C{\delta }^{1-r/k},0\le r\le s,$$

$${\varphi }_R:(r\cos \phi ,r\sin \phi )$$

$$R^{†},R_{\alpha }^{\delta }\in D(F),R^{†}\in H_{per}^4(0,2\pi ),\text{mbox}and\Vert R_{\alpha }^{\delta }-R^{†}{\Vert }_{H^2(0,2\pi )}\le {\delta }_R,$$

$$Tv=u^{ϵ}\circ {\varphi }_{\alpha }^{\delta }\text{mbox}inB,$$

$$\Vert u^{†}-u^{ϵ}{\Vert }_{L^2(D)}\le ϵ$$

$$u^{†}\in H^1(D),u^{†}{|}_{{\Omega }^{†}}\in H^3({\Omega }^{†}),\text{mbox}and{u}^{†}{|}_{D\setminus {\Omega }^{†}}\equiv 0.$$

$$v_{\alpha ,\beta }^{\delta ,ϵ}=\arg \underset{v\in D(T)}{\min }\Vert Tv-u^{ϵ}\circ {\varphi }_{\alpha }^{\delta }{\Vert }_{L^2(B)}^2+\beta \Vert \Delta v{\Vert }_{L^2(B)}^2,$$

$$\Vert u^{†}\circ {\varphi }^{†}-u^{ϵ}\circ {\varphi }_{\alpha }^{\delta }{\Vert }_{L^2(B)}\le {\delta }_U$$

$$\Vert u^{†}\circ {\varphi }^{†}-u^{ϵ}\circ {\varphi }_{\alpha }^{\delta }{\Vert }_{L^2(B)}$$

$$\Vert u^{†}\circ {\varphi }^{†}\circ ({\varphi }_{\alpha }^{\delta }{)}^{-1}-u^{†}{\Vert }_{L^2({\Omega }_{\alpha }^{\delta })}$$

$$\le \Vert u^{†}\circ {\varphi }^{†}\circ ({\varphi }_{\alpha }^{\delta }{)}^{-1}-u^{†}{\Vert }_{L^2({\Omega }_{\alpha }^{\delta }\cap {\Omega }^{†})}+\Vert u^{†}\circ {\varphi }^{†}\circ ({\varphi }_{\alpha }^{\delta }{)}^{-1}{\Vert }_{L^2({\Omega }_{\alpha }^{\delta }\setminus {\Omega }^{†})}=(i)+(ii),$$

$$(i)$$

$$\Vert u^{†}\circ {\varphi }^{†}\circ ({\varphi }_{\alpha }^{\delta }{)}^{-1}{\Vert }_{L^2({\Omega }_{\alpha }^{\delta }\setminus {\Omega }^{†})}\le \Vert u^{†}{\Vert }_{L^{\infty }(D)}|{\Omega }_{\alpha }^{\delta }\setminus {\Omega }^{†}{|}^{1/2}\le C\Vert u^{†}{\Vert }_{H^3({\Omega }^{†})}{\delta }_R^{1/2}.$$

$$(T^{\ast }f,v{)}_{H^2(B)\cap H_0^1(B)}=(\Delta T^{\ast }f,\Delta v{)}_{L^2(B)}=(f,Tv{)}_{L^2(B)}=(f,v{)}_{L^2(B)}.$$

$${\Delta }^{2}v=f\text{mbox}inB\text{mbox}withv=0\text{mbox}and\Delta v=0\text{mbox}on\partial B.$$

$$\begin{array}[]{rcl}\mathcal{R}((T^{*}T)^{0})&=&\{v\in H^{2}(B)\mid v\,=\,0\mbox{ on }\partial B\},\\ \mathcal{R}((T^{*}T)^{1/2})&=&\{v\in H^{4}(B)\mid v\,=\Delta v\,=\,0\mbox{ on }\partial B\}.\end{array}$$

$$\Vert v_{\alpha ,\beta }^{\delta ,ϵ}-v^{†}{\Vert }_{H^2(B)}\le C{\delta }_U^{\frac{2\mu }{2\mu +1}}\text{mbox}forall0\le \mu <1/8.$$

$$\beta =\max \{{\beta }_0{2}^{-n}:n\ge 0\text{mbox}andsuchthat\Vert v_{\alpha ,\beta }^{\delta ,ϵ}-v^{†}\Vert \le 2{\delta }_U\}.$$

$$\Vert Tv-u^{ϵ}\circ {\varphi }_{\alpha }^{\delta }{\Vert }_{L^2(B)}^2+\beta \Vert {\Delta }^{t}v{\Vert }_{L^2(B)}^2,$$

Full text

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., 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.

2 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\pi]\to\mathbb{R}_{+}$, we define the domain

[TABLE]

parametrized by $R$, and we denote by $\Omega^{\dagger}=\Omega_{R^{\dagger}}$ the true flow geometry. Throughout the presentation, we will assume that

[TABLE]

for some constants $2\leq k\leq 4$ and $0<r_{0}<r_{1}<1$; here $H^{s}_{per}(0,2\pi)$ is the subspace of periodic functions in the Sobolev space $H^{s}(0,2\pi)$, $s\geq 0$. Let us note that the above assumptions imply in particular that $\Omega^{\dagger}$ is of class $C^{k-1,\alpha}$, 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

[TABLE]

with forward operator $F:D(F)\subset H^{2}_{per}(0,2\pi)\to L^{2}(D)$, $R\mapsto\chi_{\Omega_{R}}$ defined on the domain $D(F)=\{R\in H^{2}_{per}(0,2\pi):r_{0}\leq R(\varphi)\leq r_{1}\}$. We further assume that the perturbations in the data are bounded by

[TABLE]

and that knowledge of the noise level $\delta$ is available. For the stable approximation of solutions to (3), we consider Tikhonov regularization with the functional

[TABLE]

We call $R_{\alpha}^{\delta}\in D(F)$ an approximate minimizer for regularization parameter $\alpha>0$, if

[TABLE]

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.

Theorem 2.1.

Assume that $R^{\dagger}\in H^{k}_{per}(0,2\pi)\cap D(F)$, $2\leq k\leq 4$. Then any approximate minimizer $R_{\alpha}^{\delta}$ of the Tikhonov functional (5) with $\alpha=\delta^{4/k}$ satisfies

[TABLE]

for $0\leq r\leq 2$, with a constant $C$ that only depends on the norm $\|R^{\dagger}\|_{H^{k}(0,2\pi)}$ of the true solution. The same estimates hold true, if $\alpha$ is chosen by a discrepancy principle, i.e.,

[TABLE]

Proof. Note that $X_{s}=H^{s}_{per}(0,2\pi)$, $s\in\mathbb{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,\widetilde{R}\in D(F)$ and define $r_{min}(\varphi)=\min\{R(\varphi),\widetilde{R}(\varphi)\}$ and $r_{max}(\varphi)=\max\{R(\varphi),\widetilde{R}(\varphi)\}$. Then

[TABLE]

This allows us to express the difference in the observations by

[TABLE]

Hence $\|R-\widetilde{R}\|_{L^{2}(0,2\pi)}\leq C\|F(R)-F(\widetilde{R})\|_{L^{2}(D)}$ for all $R,\widetilde{R}\in D(F)$, 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 $\gamma=1$. $\square$

Remark 2.2.

Note that only a simple smoothness condition for the function $R^{\dagger}$, and thus on the domain $\Omega^{\dagger}$, was required to derive a quantitative convergence result here. If in addition, an upper bound $\|R^{\dagger}\|_{H^{4}(0,2\pi)}\leq C$ is available, then one obtains

[TABLE]

with $\delta_{R}=C\delta^{1/2}$ and we may assume that the bound $\delta_{R}$ is known for further computations. The results of [3] even provide more general estimates of the form

[TABLE]

if regularization is performed in the norm of $H^{s}(0,2\pi)$ for $2\leq s\leq k/2$ and with regularization parameter $\alpha=\delta^{2s/k}$. Hence, convergence rates arbitrarily close to one can, in principle, be obtained in any desired norm if the true solution $R^{\dagger}$ is sufficiently smooth and the regularization norm is chosen sufficiently strong.

3 Velocity approximation

Let us recall from equation (1) the definition of a domain $\Omega_{R}$ parametrized by a radius function $R\in D(F)$. Using the scaling transformation

[TABLE]

with $\eta\geq k\geq 2$ and $k$ as in the previous section, we can express $\Omega_{R}$ equivalently as image $\Omega_{R}=\phi_{R}(B)$ of the unit disk $B=\{x\in\mathbb{R}^{2}:|x|<1\}$. Some important properties of this transformation are derived in Appendix A. In the following, we assume that

[TABLE]

with known bound $\delta_{R}\leq C$; see the discussion in Remark 2.2. For ease of notation, we write $\phi^{\dagger}=\phi_{R^{\dagger}}$, $\phi_{\alpha}^{\delta}=\phi_{R_{\alpha}^{\delta}}$, and denote by $\Omega^{\dagger}=\phi^{\dagger}(B)$ and $\Omega_{\alpha}^{\delta}=\phi_{\alpha}^{\delta}(B)$ the corresponding domains parametrized by the radius functions $R^{\dagger}$ and $R_{\alpha}^{\delta}$, respectively.

The velocity reconstruction from noisy velocity data $u^{\epsilon}$ can then be formulated compactly as a linear inverse problem over the reference domain, i.e.,

[TABLE]

with operator $T:H^{1}_{0}(B)\cap H^{2}(B)\to L^{2}(B)$ defined by $T(v)=v$, i.e., as simple embedding between Sobolev spaces. We assume that a bound on the data perturbations

[TABLE]

is available, where $u^{\dagger}$ denotes the true velocity field to be reconstructed, and we require that $u^{\dagger}$ is piecewise smooth and vanishes outside $\Omega^{\dagger}$, i.e.,

[TABLE]

Let us note that these assumptions are reasonable, if the flow domain $\Omega^{\dagger}$ is smooth.

Remark 3.3.

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^{\epsilon}\circ\phi_{\alpha}^{\delta}$ 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 $\phi_{\alpha}^{\delta}$ is guaranteed by Lemma A.4., which allows to associate to any approximate solution $v$ of equation (12) a velocity field $u=v\circ(\phi_{\alpha}^{\delta})^{-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

[TABLE]

where $D(T)=H^{2}(B)\cap H_{0}^{1}(B)$ by definition. 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

[TABLE]

with $\delta_{U}=C(\delta_{R}^{1/2}\|u^{\dagger}\|_{H^{3}(\Omega^{\dagger})}+\epsilon)$ and a constant $C$ that can be computed explicitly.

Proof. By the integral transformation theorem we have

[TABLE]

where we used Lemma A.4. to estimate the Jacobian $J_{\alpha}^{\delta}$ of the transformation $\phi_{\alpha}^{\delta}$ 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

[TABLE]

where we used $u^{\dagger}\equiv 0$ in $\Omega_{\alpha}^{\delta}\setminus\Omega^{\dagger}$ which follows from assumption (14). From the smoothness of $R_{\alpha}^{\delta}$, $R^{\dagger}$, and $u^{\dagger}$, and using the bound on the difference of the inverse transformations $(\phi_{\alpha}^{\delta})^{-1}$ and $(\phi^{\dagger})^{-1}$ provided by Lemma A.5., we then obtain

[TABLE]

Observe that $\|\phi^{\dagger}\|_{W^{1,\infty}(B)}\leq C\|R^{\dagger}\|_{W^{1,\infty}(0,2\pi)}\leq C^{\prime}\|R^{\dagger}\|_{H^{2}(0,2\pi)}$ by definition (10) of the transformation and continuous embedding. Using Lemma A.1. and assumption (11), we can control the area $|\Omega_{\alpha}^{\delta}\setminus\Omega^{\dagger}|$ of the geometry mismatch by $\delta_{R}$, which allows us to bound the remaining term in the above estimate by

[TABLE]

Note that the constants $C$ in the above estimates are generic and independent of $\epsilon$ and $\delta$. A combination of the bounds then yields the assertion of the lemma. $\square$

Remark 3.5.

Following Remark 2.2. and the arguments used in the proof above, the bound $\delta_{U}$ in (16) is computable in terms of the data noise levels $\delta$ and $\epsilon$ in (4) and (13), if bounds on $\|R^{\dagger}\|_{H^{4}(0,2\pi)}$ and $\|u^{\dagger}\|_{H^{3}(\Omega^{\dagger})}$ are available. For our further computations, we may thus assume that $\delta_{U}$ is known.

As a next step, we interpret standard smoothness assumptions on $u^{\dagger}$ 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^{\dagger}:=u^{\dagger}\circ\phi^{\dagger}$.

Then $v^{\dagger}\in R((T^{*}T)^{\mu})$ for all $\mu<1/8$. However, $v^{\dagger}\not\in R((T^{*}T)^{1/8})$, in general.*

Proof. We equip $D(T)=H^{2}(B)\cap H^{1}_{0}(B)$ with the norm $\|v\|:=\|\Delta v\|_{L^{2}(B)}$. Then for arbitrary $v\in H^{2}(B)\cap H^{1}_{0}(B)$ and $f\in L^{2}(B)$, we have

[TABLE]

Thus $w=T^{*}f$ is given as the unique solution of the boundary value problem

[TABLE]

From standard elliptic regularity theory [6], we can conclude that $w=T^{*}f\in H^{4}(B)$ for any $f\in L^{2}(B)$. Using $\mathcal{R}((T^{*}T)^{1/2})=\mathcal{R}(T^{*})$, see [5, Thm. 2.6], we thus arrive at

[TABLE]

From the regularity assumptions on $u^{\dagger}$ and $R^{\dagger}$, we deduce that $v^{\dagger}\in H^{3}(B)$ with $v^{\dagger}=0$ on $\partial B$, but $\Delta v^{\dagger}\neq 0$ on $\partial B$, in general. By interpolation of Sobolev spaces [10], we thus deduce that $v^{\dagger}\in R((T^{*}T)^{\mu})$ for $\mu<1/8$, but not for $\mu\geq 1/8$, in general. $\square$

Remark 3.7.

A range condition $v^{\dagger}\in\mathcal{R}((T^{*}T)^{\mu})$ would hold with $\mu=1/4$, if $\Delta v^{\dagger}=0$ at $\partial B$ would be valid additionally; this can however not be expected in general. The limiting factor for the regularity index $\mu$ 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)\cap H_{0}^{1}(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_{\alpha,\beta}^{\delta,\epsilon}$ be defined by (15) with regularization parameter $\beta=\delta_{U}^{2/(2\mu+1)}$. Then

[TABLE]

The same estimates hold for a-posteriori parameter choice by the discrepancy principle

[TABLE]

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^{\dagger}$ is sufficiently smooth and the functional in (15) is replaced by

[TABLE]

with $T:D(T)\subset H^{2t}(B)\cap H_{0}^{1}(B)\to L^{2}(B)$ defined by $Tv=v$ and $t\geq 1$ sufficiently large, one could, in principle, obtain any rate $2\mu/(2\mu+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

[TABLE]

with operator $\widetilde{T}:H^{2}(B)\cap H_{0}^{1}(B)\to L^{2}(\Omega_{\alpha}^{\delta})$ given by $\widetilde{T}v=v\circ(\phi_{\alpha}^{\delta})^{-1}$. 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.

4 Computation of the wall shear stress

Let $R\in\mathcal{D}(F)$ be a given radius function with associated transformation $\phi_{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

[TABLE]

where $J_{R}$ is the Jacobian of $\phi_{R}$, and $n_{R}(\varphi)$ is the outward pointing unit normal vector at the corresponding point $(R(\varphi)\cos\varphi,R(\varphi)\sin\varphi)\in\partial\Omega_{R}$ in the physical domain. Let us note that $n_{R}(\varphi)$ can be expressed explicitly by

[TABLE]

For ease of notation, we will identify $\partial B$ with the interval $(0,2\pi)$ in the sequel. In addition, we again define symbols $\tau_{\alpha,\beta}^{\delta,\epsilon}=\tau_{R_{\alpha}^{\delta}}(v_{\alpha,\beta}^{\delta,\epsilon})$ and $\tau^{\dagger}=\tau_{R^{\dagger}}(v^{\dagger})$ where $v^{\dagger}=u^{\dagger}\circ\phi^{\dagger}$. A combination of the results derived so far and some elementary geometric computations now leads to the following assertion.

Theorem 4.10.

Let the assumptions of Theorem 3.7 hold. Then

[TABLE]

Proof. We use the definition of $\tau_{R}(v)$ and decompose the error into the three parts

[TABLE]

Using Lemmas A.2. and A.3., and the embedding of Sobolev spaces, we obtain

[TABLE]

By the geometric arguments of Lemma A.5., the second term can be estimated by

[TABLE]

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

[TABLE]

The assertion of the theorem then follows by combination of these estimates. $\square$

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.

5 Remarks on the extension to three dimensions

For a three-dimensional flow, the wall shear stress is a tensor defined by

[TABLE]

where $\mu$ 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 $z$-coordinate, 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.

6 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.885$mm at a constant flow rate of about $6$l/min and a temperature of about $22^{\circ}$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.5$mm to $h=0.3$mm. 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 $\Omega^{\dagger}=\phi^{\dagger}(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.

6.1 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

[TABLE]

where $V_{i}$ denotes the $i$th voxel of the measurement array of size $h\!\times\!h$. The action of the forward operator for our computational tests is then defined as

[TABLE]

where $H_{\gamma}(x)=\frac{1}{\pi}\arctan\left(\frac{x}{\gamma}\right)+\frac{1}{2}$, $\gamma>0$, is a smooth approximation of the Heaviside function. In our computational tests, we thus minimize the Tikhonov functional

[TABLE]

over the finite dimensional space of radius functions

[TABLE]

Due to the smoothing with $\gamma>0$, the forward operator $F_{h,\gamma}$ 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^{\delta}$ and the true radius was $R^{\dagger}=12.9425$mm. In Table 1, we display the relative errors $\|R_{\alpha}^{\delta}-R^{\dagger}\|_{H^{2}(0,2\pi)}/\|R^{\dagger}\|_{H^{2}(0,2\pi)}$ obtained for measurement acquired with different data resolution $h$, and for different choices of the regularization parameter $\alpha$.

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^{\dagger}$ and $R_{\alpha}^{\delta}$ obtained in our numerical tests.

The maximal error in the reconstructed radius is about $0.15$mm, which is substantially less than the voxel size $h=1$mm of the data. This illustrates that sub-pixel resolution can be obtained by the proposed geometry reconstruction procedure.

6.2 Velocity approximation

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 $\Omega_{\alpha}^{\delta}=\phi_{\alpha}^{\delta}(B)$ denote the approximation of the flow domain obtained with the radius function $R_{\alpha}^{\delta}$ reconstructed in the first step with $\alpha$ 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}\cap\Omega_{\alpha}^{\delta}|>0$, the measurements can interpreted as

[TABLE]

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 $\widetilde{T}_{h}:H^{2}(B)\cap H^{1}_{0}(B)\rightarrow L^{2}(\Omega_{\alpha}^{\delta})$ with

[TABLE]

Note that only voxels $V_{i}$ with $|V_{i}\cap\,\Omega_{\alpha}^{\delta}|>0$ are required in the definition of $\widetilde{T}_{h}$. The regularized approximation for the velocity field in the reference domain is then defined as the minimizer of the Tikhonov functional

[TABLE]

over the space $V_{M}=\mbox{span}\{v\in H^{1}_{0}(B)\;|\;-\Delta v=\lambda v\mbox{ for some }\lambda\leq M\}$ of eigen functions of the Dirichlet Laplace operator on the unit disk. The solution $\widetilde{v}_{\alpha,\beta}^{\delta,\epsilon}$ can be computed efficiently by Cholesky factorization or the conjugate gradient method. In Table 2, we display the reconstruction errors obtained with our algorithm for different data resolutions $h$ and regularization parameters $\beta$.

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 $\beta$ 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.

6.3 Computation of the wall shear stress

As a final step in our numerical tests, we now utilize the reconstructed radius function $R_{\alpha}^{\delta}$ and velocity fields $v_{\alpha,\beta}^{\delta,\epsilon}$ to compute the approximation $\tau_{\alpha,\beta}^{\delta,\epsilon}$ for the wall shear stress by (19). In Table 3, we display the results obtained for optimal $\alpha$ and $\beta=1.6\cdot 10^{-4}$.

Let us note that most part of the error stems from perturbations in higher modes, which can be suppressed efficiently by application of a low-pass filter. In Figure 4, we plot the reconstruction of the wall shear stress $\tau_{\alpha,\beta}^{\delta,\epsilon}$ and its constant approximation $\overline{\tau}_{\alpha,\beta}^{\delta,\epsilon}$ against the reference value $\tau^{ref}$ obtained from numerical flow simulation [9].

Note that the average wall shear stress $\overline{\tau}_{\alpha,\beta}^{\delta,\epsilon}$ is in very good agreement with the reference value $\tau^{ref}$ obtained for the simulated data. The local variations in the reconstruction $\tau_{\alpha,\beta}^{\delta,\epsilon}$ can be explained by the turbulent variations of velocity in the data; compare with the plots in Figure 2 and 3.

7 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.

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 the German Federal and State Governments is gratefully acknowledged.

Appendix A Geometric results

In the following, we present some auxiliary results concerning geometrical details. Recall that $D(F)=\{R\in H^{2}_{per}(0,2\pi)\mid r_{0}<R<r_{1}\}$ for some $0<r_{0}<r_{1}<1$ and that $\Omega_{R}=\phi_{R}(B)$ is the image of the unit ball $B$ under the transformation

[TABLE]

with $\eta\geq 2$. The Jacobian matrix of $\phi_{R}$ is denoted by $J_{R}$, and for given functions $R_{1},R_{2}$, the subscripts transfer to the associated objects, e.g. $\phi_{1}=\phi_{R_{1}}$ and $J_{2}=J_{R_{2}}$.

As a first result, we estimate the differences of domains $\Omega_{R}=\phi_{R}(B)$ in terms of differences of their parameterizing radius functions.

Lemma A.1.

Let $R_{1},R_{2}\in D(F)$. Then

[TABLE]

Proof. Let us define $R_{max}(\varphi)=\max\{R_{1}(\varphi),R_{2}(\varphi)\}$. Then

[TABLE]

In the last step, we used the continuous embedding $H^{1}(0,2\pi)\hookrightarrow L^{\infty}(0,2\pi)$. $\square$

As a next step, we estimate differences in the normal vector $n_{R}$ defined in (20), in terms of differences in the radius function.

Lemma A.2.

Let $R_{1},R_{2}\in D(F)$. Then

[TABLE]

Proof. Let us introduce $e_{r}=(\cos(\varphi),\sin(\varphi))$ and $e_{\varphi}=(-\sin(\varphi),\cos(\varphi))$. Omitting the explicit notion of the dependence on $\varphi$, we then obtain for any angle $\varphi$ that

[TABLE]

The radial component of the difference can be estimated by

[TABLE]

In a similar way, one can estimate $\lambda_{\varphi}$, and by Sobolev’s embedding theorem, we obtain

[TABLE]

which already yields the desired estimate. $\square$

The following result ensures smoothness of the transformations $\phi_{R}$ whenever the radius function $R$ is sufficiently smooth.

Lemma A.3.

Let $R_{1},R_{2}\in D(F)\cap H^{k}_{per}(0,2\pi)$ for $k\leq\eta$ with $\eta\geq 2$ as in (10). Then

[TABLE]

where $\phi_{1}=\phi_{R_{1}}$ with $\phi_{R}$ defined in (10). Moreover,

[TABLE]

Proof. The continuity of $\phi_{1}$ is obvious and the Jacobian $J_{1}=J_{R_{1}}$ obtained by derivation with respect to coordinates $x=(r\cos\phi,r\sin\phi)$ reads

[TABLE]

Since $\eta\geq 2$, the Jacobian can be seen to be continuous also in $r=0$, and we have

[TABLE]

This shows the first estimate for $k=1$. The assertion for higher order derivatives of $\phi_{1}$ follow in a similar way. From the formula (10), one can see that $\phi_{R}$ is affine linear in $R$. Hence the second estimate follows directly from the first. $\square$

As a next step, we show that the transformations $\phi_{R}$ defined in (10) are invertible.

Lemma A.4.

Let $R\in D(F)$. Then the transformation $\phi_{R}$ defined in (10) is a diffeomorphism with inverse transformation $\psi_{R}=(\phi_{R})^{-1}$ and

[TABLE]

Proof. Using (27), we can estimate the determinant by

[TABLE]

The existence of an inverse transformation $\psi_{R}=(\phi_{R})^{-1}$ then follows form the implicit function theorem and the Jacobian of the inverse mapping is $J_{\psi_{R}}\circ\phi_{R}=J_{R}^{-1}$. The bounds for $J_{\psi_{R}}$ can then be deduced in an elementary way. $\square$

Using the previous results, we can also bound differences in the inverse mappings.

Lemma A.5.

Let $R_{1},R_{2}\in D(F)$ and assume that $R_{1}\in H^{3}_{per}(0,2\pi)$. Furthermore, let $\psi_{1}$, $\psi_{2}$ denote the corresponding inverse transformations with Jacobians $J_{\psi_{1}}$, $J_{\psi_{2}}$. Then

[TABLE]

and the difference in the Jacobians can be bounded by

[TABLE]

Proof. Step 1: To show (29), let $x\in\Omega_{1}\cap\,\Omega_{2}$. Then there are $x_{1},x_{2}\in B$ with $\phi_{1}(x_{1})=\phi_{2}(x_{2})=x$. Since the transformations $\phi_{1}$ and $\phi_{2}$ preserve angles, the angular parts $\varphi(x_{1})=\varphi(x_{2})=\varphi(x)=:\varphi$ are equal. $\phi_{2}$ has an inverse $\psi_{2}$, hence

[TABLE]

where the defect is given by

[TABLE]

Since $x$ and $dx$ have the same angular coordinate and $\Omega_{2}$ is star shaped with respect to the origin, we have $\{x+t\;dx\;|\;t\in[0,1]\}\subset\Omega_{2}$. The mean value theorem yields

[TABLE]

The assertion (29) follows by the estimate of Lemma A.4. Step 2: Show (30). Starting at (27) elementary calculus yields

[TABLE]

where the constant $C$ depends only on $\|R_{1}\|_{L^{\infty}(0,2\pi)}$, $\|R_{2}\|_{L^{\infty}(0,2\pi)}$, $r_{0}$ and $\eta$. For the inverse Jacobians, we further conclude that

[TABLE]

Since the inverse Jacobian of $\phi_{1}$ has the representation $J_{\phi_{1}}^{-1}=\tilde{J}_{\phi_{1}}/\det J_{\phi_{1}}\;$, where $\tilde{J}_{\phi_{1}}$ is a rearrangement of $J_{\phi_{1}}$ and all expressions are continuously differentiable, one can verify without difficulty that $J^{-1}_{\phi_{1}}\in W^{1,\infty}(B)$ with the associated norm bounded in terms of $\eta$, $r_{0}$, and $\|R_{1}\|_{H^{3}(0,2\pi)}$. Therefore we arrive at (30) by

[TABLE]

This completes the proof of the second estimate of the lemma. $\square$

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.

B.1 Magnitude data

The magnitude raw data $m^{\delta}_{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.

Based on the two values $m_{0}$, $m_{1}$, the normalized magnitude $m^{\delta}$ are then defined by

[TABLE]

with truncation function $T(x)=\max(0,\min(1,x))$. These measurements then represent an approximation for the characteristic function of the flow domain.

B.2 Velocity data

The raw data acquired in phase contrast magnetic resonance velocimetry can be interpreted as

[TABLE]

where $\rho$ is the propton density and $v_{enc}$ the velocity encoding parameter. An average value of the velocity in the voxel $V_{i}$ is then recovered as the phase of this signal, i.e.,

[TABLE]

Phase unwrapping may be required, if the maximal velocity is larger then $v_{enc}$. In the absence of data noise and assuming that the density is given by $\rho=c\chi_{\Omega}$, one can use Taylor approximation to deduce that

[TABLE]

where $\Omega$ is the flow domain and $u$ the true flow velocity. This is the measurement model used in the numerical tests. Let us note that the phase retrieval in (31) is particularly sensitive to data perturbations, if $\textrm{abs}(d^{\delta}|_{V_{i}})$ is small, which is the case close to the boundary and outside the flow domain. Therefore, particularly large noise in the velocity data is expected for these voxels; compare with Figure 2.