Fast reconstruction of 3D blood flows from Doppler ultrasound images and reduced models

TL;DR

This paper introduces improved reduced models for fast 3D blood flow reconstruction from Doppler ultrasound, enhancing accuracy over traditional linear methods through task-adapted, potentially nonlinear mappings.

Contribution

It proposes new reduced space construction methods tailored to blood flow reconstruction, outperforming classical linear PBDW approaches.

Findings

Task-adapted reduced models improve reconstruction accuracy.

Nonlinear mappings can outperform linear approaches.

Numerical experiments validate the proposed methods.

Abstract

This paper deals with the problem of building fast and reliable 3D reconstruction methods for blood flows for which partial information is given by Doppler ultrasound measurements. This task is of interest in medicine since it could enrich the available information used in the diagnosis of certain diseases which is currently based essentially on the measurements coming from ultrasound devices. The fast reconstruction of the full flow can be performed with state estimation methods that have been introduced in recent years and that involve reduced order models. One simple and efficient strategy is the so-called Parametrized Background Data-Weak approach (PBDW). It is a linear mapping that consists in a least squares fit between the measurement data and a linear reduced model to which a certain correction term is added. However, in the original approach, the reduced model is built a priori and independently of the reconstruction task (typically with a proper orthogonal decomposition or a greedy algorithm). In this paper, we investigate the construction of other reduced spaces which are built to be better adapted to the reconstruction task and which result in mappings that are sometimes nonlinear. We compare the performance of the different algorithms on numerical experiments involving synthetic Doppler measurements. The results illustrate the superiority of the proposed alternatives to the classical linear PBDW approach.

Tables

Table 1. Table 1: Computational cost of each method.

	POD-lin	P-POD-aff	P-Greedy-aff	P-DB-aff
Offline: Database generation	$𝒪\left({𝒩}^{\alpha }\#\stackrel{~}{}{ℳ}_{\text{train}}\right)$	$𝒪\left({𝒩}^{\alpha }\#\stackrel{~}{}{ℳ}_{\text{train}}\right)$	$𝒪\left({𝒩}^{\alpha }\#\stackrel{~}{}{ℳ}_{\text{train}}\right)$	$𝒪\left({𝒩}^{\alpha }\#\stackrel{~}{}{𝒟}_{\text{train}}\right)$
Offline: Construction reduced model	$𝒪\left(𝒩\#\stackrel{~}{}{ℳ}_{\text{train}}^2\right)$	$𝒪\left(𝒩{\left(\#\stackrel{~}{}{ℳ}_{\text{train}}^{\left(k\right)}\right)}^2\right)$	$𝒪\left(n𝒩\#\stackrel{~}{}{ℳ}_{\text{train}}^{\left(k\right)}\right)$	$𝒪\left(𝒩\#\stackrel{~}{}{𝒟}_{\text{train}}^{\left(k\right)}+m^3\right)$
Online: Reconstruction	$𝒪\left(n^2+n𝒩\right)$	$𝒪\left(n^2+n𝒩\right)$	$𝒪\left(n^2+n𝒩\right)$	$𝒪\left(n^2+nm\#\stackrel{~}{}{𝒟}_{\text{train}}^{\left(k\right)}+n𝒩\right)$

Equations

$${\ell }_i(u),i=1,\dots ,m,$$

$$W_m=\mathrm{span}\{{\omega }_1,\dots ,{\omega }_m\}\subset V.$$

$$\omega =P_{W_m}u.$$

$$A:W_m\mapsto V$$

$$\mathcal{P}(u,y)=0,$$

$$\mathcal{M}:=\{u(y)\in V:y\in \mathrm{Y}\},$$

$$E_{\mathrm{wc}}(A,\mathcal{M})=\underset{u\in \mathcal{M}}{\sup }\Vert u-A(P_{W_m}u)\Vert .$$

$$E_{\mathrm{ms}}^2(A,\mathcal{M})=\mathbb{E}\left(\Vert u-A(P_{W_m}u){\Vert }^2\right)={\int }_V\Vert u-A(P_{W_m}u){\Vert }^{2}dp(u),$$

$${\epsilon }_n:=\underset{u\in \mathcal{M}}{\sup }\mathrm{dist}(u,V_n),\text{or}{\delta }_n^2:=\mathbb{E}\left(\mathrm{dist}(u,V_n{)}^2\right)$$

$$A_{m,n}^{\text{(pbdw)}}:W_m\to V$$

$$\underset{u\in \omega +W^{\perp }}{\min }\mathrm{dist}(u,V_n).$$

$$\beta (X,Y):=\underset{x\in X}{\inf }\underset{y\in Y}{\sup }\frac{⟨x,y⟩}{\Vert x\Vert \Vert y\Vert }=\underset{x\in X}{\inf }\frac{\Vert P_{Y}x\Vert }{\Vert x\Vert }\in [0,1]$$

$$A_{m,n}^{\text{(pbdw)}}(\omega )=u_{m,n}^{\ast }(\omega ):=\underset{u=\omega +W^{\perp }}{\mathrm{argmin}}\mathrm{dist}(u,V_n).$$

$$u_{m,n}^{\ast }(\omega )=v_{m,n}^{\ast }(\omega )+\omega -P_{W_m}{v}_{m,n}^{\ast }(\omega )$$

$$v_{m,n}^{\ast }(\omega )={\left(P_{V_n|W_m}{P}_{W_m|V_n}\right)}^{-1}{P}_{V_n|W_m}(\omega ),$$

$$e_{m,n}^{\text{(wc, pbdw)}}=E_{\mathrm{wc}}(A_{m,n}^{\text{(pbdw)}},\mathcal{M})\le {\beta }^{-1}(V_n,W_m)\underset{u\in \mathcal{M}}{\max }\mathrm{dist}(u,V_n\oplus (V_n^{\perp }\cap W_m))\le {\beta }^{-1}(V_n,W_m){\epsilon }_n,$$

$$e_{m,n}^{\text{(ms, pbdw)}}=E_{\mathrm{ms}}(A_{m,n}^{\text{(pbdw)}},\mathcal{M})$$

$$\le {\beta }^{-1}(V_n,W_m)\mathbb{E}{\left(\mathrm{dist}(u,V_n\oplus (V_n^{\perp }\cap W_m{)}^2\right)}^{1/2}$$

$$\le {\beta }^{-1}(V_n,W_m){\delta }_n,$$

$$e_{m,n_{\mathrm{wc}}^{\ast }}^{\text{(wc, pbdw)}}=\underset{1\le n\le m}{\min }{e}_{m,n}^{\text{(wc, pbdw)}},\text{or}{e}_{m,n_{\mathrm{ms}}^{\ast }}^{\text{(ms, pbdw)}}=\underset{1\le n\le m}{\min }{e}_{m,n}^{\text{(ms, pbdw)}}.$$

$$A_{m,n}^{\text{(aff)}}(\omega ):=\underset{u\in \omega +W^{\perp }}{\mathrm{argmin}}\mathrm{dist}(u,\bar{u}+V_n),$$

$$A_{m,n}^{\text{(aff)}}(\omega )=\bar{u}+A_{m,n}^{\text{(pbdw)}}(\omega -\bar{\omega }),$$

$$\Vert u-A_{m,n}^{\text{(aff)}}(\omega )\Vert$$

$$\le {\beta }^{-1}(V_n,W)\mathrm{dist}(u,\bar{u}+V_n)$$

$$e_{m,n}^{\text{(wc, aff)}}:=E_{\mathrm{wc}}(A_{m,n}^{\text{(aff)}},\mathcal{M})\le {\beta }^{-1}(V_n,W_m){\epsilon }_n^{\text{(aff)}},$$

$$e_{m,n}^{\text{(ms, aff)}}:=E_{\mathrm{ms}}(A_{m,n}^{\text{(aff)}},\mathcal{M}):=\mathbb{E}{\left(\Vert u-A_{m,n}^{\text{(aff)}}(P_{W_m}u){\Vert }^2\right)}^{1/2}\le {\beta }^{-1}(V_n,W_m){\delta }_n^{\text{(aff)}},$$

$${\epsilon }_n^{\text{(aff)}}:=\underset{u\in \mathcal{M}}{\sup }\mathrm{dist}(u,\bar{u}+V_n),\text{or}({\delta }_n^{\text{(aff)}}{)}^2:=\mathbb{E}\left(\mathrm{dist}(u,\bar{u}+V_n{)}^2\right)$$

$$e_{m,n}^{\text{(wc, aff, k)}}=E_{\mathrm{wc}}(A_{m,n}^{\text{(aff)}},{\mathcal{M}}^{(k)})\le {\beta }^{-1}(V_n^{(k)},W_m){\epsilon }_n^{(k)},$$

$$e_{m,n}^{\text{(ms, aff, k)}}=E_{\mathrm{ms}}(A_{m,n}^{\text{(aff)}},{\mathcal{M}}^{(k)}):=\mathbb{E}{\left(\Vert u-A_{m,n}^{\text{(aff)}}(P_{W_m}u){\Vert }^2\right)}^{1/2}\le {\beta }^{-1}(V_n^{(k)},W_m){\delta }_n^{(k)}.$$

$$e_{m,n_{\mathrm{wc}}^{\ast }(k)}^{\text{(wc, aff, k)}}=\underset{1\le n\le m}{\min }{e}_{m,n}^{\text{(wc, aff, k)}},\text{or}{e}_{m,n_{\mathrm{ms}}^{\ast }(k)}^{\text{(ms, aff, k)}}=\underset{1\le n\le m}{\min }{e}_{m,n}^{\text{(ms, aff, k)}}.$$

$$e_m^{\text{(wc, aff)}}=\underset{1\le n\le m}{\max }{e}_{m,n_{\mathrm{ms}}^{\ast }(k)}^{\text{(wc, aff, k)}},\text{or}{e}_{m,n_{\mathrm{ms}}^{\ast }(k)}^{\text{(ms, aff, k)}}=\sum _{k=1}^K{\omega }_k{e}_{m,n_{\mathrm{ms}}^{\ast }(k)}^{\text{(ms, aff, k)}},$$

$$A_{m,n}^{\text{(dd)}}(\omega ):=\underset{u\in \omega +W^{\perp }}{\mathrm{argmin}}\mathrm{dist}(u,\bar{u}+V_n(\omega )),$$

Full text

1 Introduction

Developing artificial intelligence-driven tools to assist doctors in medical decisions and diagnosis requires to solve efficiently and in a reliable manner data assimilation and inverse problems from biomedical applications. These problems share the following general common features, which will guide our subsequent developments:

The available data is often corrupted by noise and obtained with medical imaging techniques, which have the advantage of being non-invasive.

• 2.

There may be morphological constraints that prevent from measuring at specific locations.

• 3.

In some cases, the device may not be able to measure directly the desired quantity of interest (QoI), and a complex post-processing may be required to obtain an estimate of it.

• 4.

Sometimes, the desired QoI is a prediction or a forecast ahead in time, which requires a specific treatment to be inferred.

In this paper, we develop state estimation techniques which aim at taking the above features into account, and which involve reduced modeling of parametrized Partial Differential Equations (PDEs). As an illustrative example of how the reconstruction techniques that we deploy address some of the above points, in this paper, we consider a haemodynamics inverse problem: the real-time reconstruction of the full 3D blood velocity field in an artery from Doppler ultrasound images taken on a restricted portion of the artery. This application is of interest in its own right since a real-time reconstruction of the full flow would enrich the available information for medical diagnosis.

Haemodynamics forward and inverse problems is a wide multidisciplinary topic with a long history, and a complete overview of the numerous existing contributions would go beyond the scope of the present paper. However, as a brief summary of the state of the art for our application on haemodynamics reconstructions using medical images, we could start by citing [2], where a method based on quadratures is used to estimate the velocity from Doppler ultrasound. In [3, 4, 5, 6, 7] sequential data estimation methods are used to provide a description of the haemodynamics in different portions of the vascular tree. A patient specific fluid-structure interaction simulation of the left ventricle is proposed in [8] based on Magnetic Resonance Imaging data. The calibration of computational fluid mechanics simulations was also investigated in [9]. A multi-fidelity approach and a Bayesian framework are detailed in [10]. Recently, some works have explored kernel and deep learning based approaches to perform data assimilation [11, 12, 13].

In this work, we use methods based on reduced-order modelling to obtain reconstructions in close to real-time. This approach has attracted considerable attention and numerous mathematical developments and applications have been developed, see [14, 15, 16, 17, 18, 19, 20, 21, 22, 23]. These methods share connections with other well-known approaches such as 3D and 4D-Var (see [24, 25]) or the Partial Spline Model ([26, Chapter 9]). We refer to [27, 22] for further details on this point. Regarding inverse haemodynamics problems, prior contributions involving reduced-order modeling are, e.g., [28, 29].

Our starting point is the Parametrized Background Data-Weak originally introduced in [1]. The method consists in a linear mapping which is built from a least squares fit between the measurements and a given reduced space and an additional term that is capable of correcting model bias to some extent. In the original approach of PBDW, the reduced model is built a priori and independently of the reconstruction task (typically with a Proper Orthogonal Decomposition or a greedy algorithm). In this paper, we investigate the construction of other reduced spaces which are built to be better adapted to the reconstruction task and which result in mappings that are sometimes nonlinear. We compare the performance of the different algorithms in the context of the fast reconstruction of the velocity field in an artery. The results illustrate the superiority of the proposed alternatives to the classical approach involving linear reduced bases.

The paper is organized as follows, in section 2 we describe mathematically what we understand by state estimation problems and reconstruction algorithms (section 2.1). We next introduce several ways of building a reduced order space, some of which take the measurements into account, henceforth leading to non-linear reconstruction methods (section 2.2). We give details on the practical implementation and computational costs in section 3. We next apply the methodology to the reconstruction of 3D blood flows from Doppler images (section 4). The example is semi-realistic, meant to be a first (still idealized) step towards realistic applications. The main assumptions are the following. First, the Doppler images are noiseless and synthetically computed but they have been generated with a realistic model with respect to the real involved physics and measurement devices. This assumption is mainly due to our lack of real Doppler measurements and also because of the very involved space-time structure of the noise in Doppler images (which is a research topic in itself [30, 31, 32]). Second, we assume that there is no model misfit even if the method can correct it to some extent. However, we use a rather realistic model for flows in large arteries based on incompressible Navier-Stokes equations, with boundary conditions used in haemodynamics simulations (see [33]). The precise parametric PDE model upon which we build our reduced models is explained in section 4.1. Section 4.2 gives some details on Doppler ultrasound imaging and the way in which we have incorporated it to our methodology. Finally, sections 4.3 and 4.4 present numerical results in two test cases. They illustrate the superiority of the proposed reconstruction algorithms with respect to the classical linear PBDW. In addition to this, our second example also shows that the method can be used to estimate quantities of interest which could be helpful in the detection of an arterial blockage.

We conclude this introduction by emphasizing that the assumptions in our computations with no model misfit and synthetic, noiseless measurements are not inherent limitations of the present approach. The focus lies rather on the reduced space construction and the haemodynamics application. The model misfit is automatically corrected and measurement noise can be added to the recostruction pipeline (see [27, 21, 34] for works on measurement noise). We also emphasize the novelty that the current approach represents for the field of blood flow reconstruction. Indeed, estimations performed routinely by medical doctors usually do not involve full reconstructions (perhaps due to the lack of efficient techniques). They often reduce to simple pointwise estimations of the peak and average velocities of the imaged anatomical part. The development of more refined estimations of the velocity field are actually an active research topic, especially for flows in the heart cavities. The only fast technique that we are aware of is based solely on the divergence equation ($\nabla\cdot u=0$) on the imaging plane [35, 36]. Despite its rapidity and simplicity, the contribution of the out-of-plane velocity and the momentum conservation are neglected, leading to very non-physical approximations. Another existing alternative are classical inverse problems in which a joint state-parameter estimation is performed. Their main drawback is the large computational times, which is incompatible to the clinical practice.

2 Reconstruction methods

We start by introducing the state reconstruction methods that we use in the present work. We focus on explaining alternative ways of building reduced basis that are better tailored for the task of state estimation than the usual reduced bases. Our presentation is done for noiseless measurements and assumes that there is no model error. As already brought up, this is due to the fact that we lack from real measurements for our targeted application.

2.1 State estimation and recovery algorithms

Let $\Omega$ be a domain of $\mathbb{R}^{d}$ for a given dimension $d\geq 1$ and let $V$ be a Hilbert space defined over $\Omega$, with inner product $\langle\cdot,\cdot\rangle$ and norm $\|\cdot\|$. Our goal is to recover an unknown function $u\in V$ from $m$ measurement observations

[TABLE]

where the $\ell_{i}$ are linearly independent linear forms over $V$. In many applications, each $\ell_{i}$ models a sensor device which is used to collect the measurement data $\ell_{i}(u)$. If the observations come in the form of an image as in the application of this paper, each $\ell_{i}$ may represent the response of the system in a given pixel. The Riesz representers of the $\ell_{i}$ are denoted by $\omega_{i}$ and span an $m$-dimensional space

[TABLE]

The observations $\ell_{1}(u),\dots,\ell_{m}(u)$ are thus equivalent to knowing the orthogonal projection

[TABLE]

In this setting, the task of recovering $u$ from the measurement observation $\omega$ can be viewed as building a recovery algorithm

[TABLE]

such that $A(P_{{W_{m}}}u)$ is a good approximation of $u$ in the sense that $\|u-A(P_{{W_{m}}}u)\|$ is small.

Recovering $u$ from the measurements $P_{W_{m}}u$ is a very ill-posed problem since there are infinitely many $v\in V$ such that $P_{W_{m}}v=\omega$. It is thus necessary to add some a priori information on $u$ in order to recover the state up to a guaranteed accuracy. We are motivated by the setting where $u$ is a solution to some parameter-dependent PDE of the general form

[TABLE]

where $\mathcal{P}$ is a differential operator and $y$ is a vector of parameters that describes some physical property and lives in a given set $\mathrm{Y}\subset\mathbb{R}^{p}$. Therefore, our prior on $u$ is that it belongs to the set

[TABLE]

which is sometimes referred to as the solution manifold. The performance of a recovery mapping $A$ is usually quantified in two ways:

If the sole prior information is that $u$ belongs to the manifold $\mathcal{M}$, the performance is usually measured by the worst case reconstruction error

[TABLE]

• 2.

In some cases $u$ is described by a probability distribution $p$ on $V$ supported on $\mathcal{M}$. This distribution is itself induced by a probability distribution on $\mathrm{Y}$ that is assumed to be known. In this Bayesian-type setting, the performance is usually measured in an average sense through the mean-square error

[TABLE]

and it naturally follows that $E_{\mathrm{ms}}(A,\mathcal{M})\leq E_{\mathrm{wc}}(A,\mathcal{M})$.

2.2 Linear and nonlinear algorithms using reduced modeling

Reduced models are a family of methods that produce each a hierarchy of spaces $(V_{n})_{n\geq 1}$ that approximate the solution manifold well in the sense that

[TABLE]

decays rapidly as $n$ grows for certain classes of PDEs. Several methods exist to build these spaces among which stand the reduced basis method (see [37]), the (Generalized) Empirical Interpolation Method (see [38, 18, 39]), Proper Orthogonal Decomposition (POD, [40, 41]) and low-rank methods (see [42, 43]).

Linear reconstruction algorithms that make use of reduced spaces $V_{n}$ are the Generalized Empirical Interpolation Method (GEIM) introduced in [18] and further analyzed in [19, 39] and the Parametrized Background Data-Weak Approach (PBDW) introduced in [1] and further analyzed in [44]. Note that some modified versions have been proposed to address measurement noise (see, e.g., [21, 27]) and other recovery algorithms involving reduced modelling have also been recently proposed (see [22]).

2.2.1 PBDW, a linear recovery algorithm

Given a measurement space ${W_{m}}$ and a reduced model $V_{n}$ with $1\leq n\leq m$, the PBDW algorithm

[TABLE]

gives for any $\omega\in{W_{m}}$ a solution of

[TABLE]

Denoting

[TABLE]

for any pair of closed subspaces $(X,Y)$ of $V$, the above optimization problem has a unique minimizer

[TABLE]

as soon as $n\leq m$ and $\beta({V_{n}},{W_{m}})>0$. We adhere to these two assumptions in the following.

As proven in section 3.2, an explicit expression of $u_{m,n}^{*}(\omega)$ is

[TABLE]

with

[TABLE]

where, for any pair of closed subspaces $(X,Y)$ of $V$, $P_{X|Y}:Y\to X$ is the orthogonal projection into $X$ restricted to $Y$. The invertibility of the operator $P_{{V_{n}}|{W_{m}}}P_{{W_{m}}|{V_{n}}}$ is guaranteed under the above conditions.

Formula (6) shows that $A^{\text{(pbdw)}}_{n}$ is a bounded linear map from ${W_{m}}$ to ${V_{n}}\oplus(W_{m}\cap V_{n}^{\perp})$. Depending on whether ${V_{n}}$ is built to address the worst case or mean square error, the reconstruction performance is bounded by

[TABLE]

or

[TABLE]

Note that $\beta({V_{n}},{W_{m}})$ can be understood as a stability constant. It can also be interpreted as the cosine of the angle between ${V_{n}}$ and ${W_{m}}$. The error bounds involve the distance of $u$ to the space $V_{n}\oplus(V_{n}^{\perp}\cap{W_{m}})$ which provides slightly more accuracy than the reduced model ${V_{n}}$ alone. In the following, to ease the reading we will write errors only with the second type of bounds that do not involve this correction part on $V_{n}^{\perp}\cap{W_{m}}$.

An important observation is that for a fixed measurement space $W_{m}$ (which is the setting in our numerical tests), the errors $e_{m,n}^{\text{(wc, pbdw)}}$ and $e_{m,n}^{\text{(ms, pbdw)}}$ reach a minimal value $e_{m,n^{*}_{\mathrm{wc}}}^{\text{(wc, pbdw)}}$ and $e_{m,n^{*}_{\mathrm{ms}}}^{\text{(ms, pbdw)}}$ as the dimension $n$ varies from 1 to $m$. This behavior is due to the trade-off between the increase of the approximation properties of $V_{n}$ as $n$ grows and the degradation of the stability of the algorithm, given here by the decrease of $\beta(V_{n},{W_{m}})$ to 0 as $n\to m$. As a result, the best reconstruction performance with PBDW is given by

[TABLE]

We finish this section with several remarks:

We do not consider measurement noise and model error. However, the PBDW algorithm can correct to some extent the model misfit (through the term $\eta^{*}$, see 3.2). Also, some extensions of the current setting have been proposed to address measurement noise (see [27, 21, 34]). 2. 2.

In our application, the manifold $\mathcal{M}$ will be a family of incompressible fluid flow solutions of a parametric incompressible Navier-Stokes equation. In this case, $V_{n}$ is usually built such that all functions $v\in V_{n}$ satisfy the divergence-free condition $\nabla\cdot v=0$. As a result, a reconstruction with only $v^{*}_{n}$ (see (7)) will yield a divergence-free approximation of the flux. Note however that the full PBDW reconstruction $u^{*}_{n}$ from equation (6) does not guarantee this property since the component $\omega-P_{W_{m}}v^{*}_{m,n}(\omega)$ may not be divergence free. The mass conservation of the reconstruction could be enforced, for instance, by considering its projection on divergence-free fields similar to the one used in fractional methods (like Chorin-Temam) for fluid flows. 3. 3.

Note that in the present setting the measurement space ${W_{m}}$ is fixed and we will adhere to this assumption in the rest of the paper. This is reasonable since usually the nature and location of the sensors is fixed in our application following the experience of the medical doctors. A different, yet related problem, would be to optimize the choice of the measurement space ${W_{m}}$ and we refer to [39, 45, 46, 47, 48, 49] for works on this topic.

2.2.2 Alternative reconstruction algorithms

Taking PBDW as a starting point, we next describe two different strategies to build recovery algorithms. They all incorporate an affine extension of PBDW that we explain next. One motivation to introduce it is because it can give more accuracy when the solution manifold $\mathcal{M}$ is not localized near the origin. This typically happens when the state $u$ is a perturbation of a nominal state $\bar{u}\in V$.

Affine PBDW: Given an average or nominal state $\bar{u}\in V$, we can formulate the equivalent affine version of PBDW, which reads

[TABLE]

As shown in 3.2, this algorithm can also be written as

[TABLE]

where $\bar{\omega}=P_{W_{m}}\bar{u}$ and $A_{m,n}^{\text{(pbdw)}}$ is the linear PBDW algorithm of the previous section. We can easily see through the last expression that $A_{m,n}^{\text{(aff)}}$ is an affine mapping from ${W_{m}}$ to $\bar{u}+V_{n}\oplus({W_{m}}\cap V_{n}^{\perp})$. An interesting feature of this affine extension is that if $\omega=0$, the algorithm yields a nonzero reconstruction in ${W_{m}}^{\perp}$ based on the nominal state $\bar{u}$, $A_{m,n}^{\text{(aff)}}(0)=\bar{u}-A_{m,n}^{\text{(pbdw)}}(\bar{\omega})$. This is in contrast to the linear version where $A_{m,n}^{\text{(pbdw)}}(0)=0$.

Proceeding similarly as before, for a given $u\in\mathcal{M}$, the error is bounded by (see [23])

[TABLE]

and

[TABLE]

or

[TABLE]

where

[TABLE]

In the following, to simplify notation, we will use $\varepsilon_{n}$ and $\delta_{n}$ to denote either the error in the linear PBDW or its affine version since the reasoning and the estimates that will be derived next have the same form.

Partition of $\mathcal{M}$: Our first strategy stems from the fact that we may have very few observations ($m$ small) and since we must have $n\leq m$ for reconstruction, we may not have enough approximation power in the interval where $n$ can range. Also, the approximation errors $(\varepsilon_{n})_{n}$ or $(\delta_{n})_{n}$ provided by reduced basis may not always decrease rapidly to zero (the Kolmogorov $n$-width of $\mathcal{M}$ may decrease slowly). However, the physical structure of the problem could give a natural decomposition of the manifold $\mathcal{M}$ into different subdomains $\mathcal{M}^{(k)}$ that are better adapted for model reduction in the sense that the errors $(\varepsilon^{(k)}_{n})_{n}$ or $(\delta^{(k)}_{n})_{n}$ may decrease faster. Works based on that partition principle have been developed in the context of model reduction of forward problems [50, 51, 52, 53]) and we propose to adapt it for data assimilation. A simple setting where we can decompose the manifold in our context is when we work in an application for which it is possible to know exactly a subset of entries in y for any target $u(y)$ during the online phase, say $\overline{y}\in\mathbb{R}^{\overline{p}}$ (with $\overline{p}<p$). Given $\overline{y}$ and a data-base generated from the governing PDE we may produce a disjoint union of $K$ subsets $\mathrm{Y}^{(k)}$ which yields a decomposition of $\mathcal{M}$ into subsets $\mathcal{M}^{(k)}=u(\mathrm{Y}^{(k)})$. We can thus build reduced models $(V_{n}^{(k)})_{n=1}$ for each subset $\mathcal{M}^{(k)}$ and then reconstruct with the linear or affine PBDW.

Proceeding similarly as in the previous section, the reconstruction performance on subset $\mathcal{M}^{(k)}$ is

[TABLE]

or

[TABLE]

The best reconstruction performance for $\mathcal{M}^{(k)}$ is thus

[TABLE]

It follows that the performance in $\mathcal{M}=\cup_{k=1}^{K}\mathcal{M}^{(k)}$ is

[TABLE]

where $\omega_{k}=p(u\in\mathcal{M}^{(k)})$.

Data-based reduced models: The second strategy is motivated by the observation that the reduced models $V_{n}$ or $V_{n}^{(k)}$ of the previous approaches are built independently of the given measurement space $W_{m}$ and also of the given measurement observation $\omega$. For a given $\omega\in{W_{m}}$, we can build a data-driven $V_{n}(\omega)$ with an orthogonal matching pursuit greedy algorithm (OMP) that we explain next. Once this reduced model has been computed, we reconstruct with the data-driven affine version of PBDW,

[TABLE]

where the difference with respect to (10) is that now ${V_{n}}$ depends on $\omega$. The reconstruction performance of this algorithm is bounded by

[TABLE]

in the worst case setting. Similarly as before, $e_{m,n}^{\text{(wc, dd)}}$ reaches a minimum $e_{m,n^{*}_{\text{(dd)}}}^{\text{(wc, dd)}}$ when $n$ varies from 1 to $m$. Since $V_{n}$ is now adapted to the measurement observations, we expect that the current algorithm performs better than its classical linear counterpart.

We now present the OMP greedy algorithm that we propose. Let

[TABLE]

be the set of normalized functions from $\mathcal{M}$. If $\bar{u}=0$, the first element $\varphi_{1}$ is chosen as

[TABLE]

For $n>1$, given $V_{n}=\mathrm{span}\{\varphi_{1},\ldots,\varphi_{n}\}$, we select

[TABLE]

where $P_{{W_{m}}}V_{n}=\mathrm{span}\{P_{{W_{m}}}\varphi_{1},\ldots,P_{{W_{m}}}\varphi_{n}\}$. We set $V_{n+1}=\mathrm{span}\{V_{n},\varphi_{n+1}\}$.

Note that all operations in this algorithm are done in the space ${W_{m}}$. Hence we can do all calculations in $\mathbb{R}^{m}$, which makes this algorithm be very fast since it does not involve computations with functions from the whole space $V$ (see 3.3 for further details on its implementation).

In the case $\bar{u}\neq 0$, we introduce $\bar{\omega}=P_{W_{m}}\bar{u}$ and the shifted set

[TABLE]

Now it suffices to apply the previous greedy algorithm to the target function $\omega-\bar{\omega}$ instead of $w$ and do the search over $\delta_{\bar{u}}\mathcal{D}$ instead of $\mathcal{D}$.

Remark 1

Note that both nonlinear approaches are very different in nature: the first is based on localization and the second is purely data-driven. Their performance will thus be very dependent on the physical problem under consideration and the measurement setting: the nature of ${W_{m}}$, the nature of the actual measurements $\omega=P_{W_{m}}u$, the nature of the manifold and its alignment with ${W_{m}}$ locally at $\omega$. Note also that both approaches can be combined and we can consider a piece-wise data-driven algorithm (its performance has actually been tested in the numerical tests below).

3 Computational details of the methods

3.1 Discretizations

In the previous section and for the sake of clarity, we have given an idealized description of PBDW and the other reconstruction methods where we do not mention certain discretization aspects that inevitably come into play. To start with, the snapshot functions $u\in\mathcal{M}$ cannot be computed exactly but only up to some tolerance. One typical instance (which is used in our work) is the finite element method which gives an approximation $u_{h}\in V_{h}$ of $u\in V$ where $V_{h}$ is a finite element space with $\mathcal{N}$ degrees of freedom. Therefore the computable states are elements of the perturbed manifold

[TABLE]

The computation of reduced models $V_{n}$ involves a finite training subset $\widetilde{}\mathcal{M}_{\text{train}}\subset\mathcal{M}_{h}$ of snapshots, and we denote by $\#\widetilde{}\mathcal{M}_{\text{train}}$ its cardinality. As a consequence, all the reduced models are low dimensional subspaces of $V_{h}$. Note that the fact that the true states do not belong to $\mathcal{M}_{h}$ can be interpreted as a model bias.

The methodology also requires the computation of the Riesz lifts $\omega_{i}$ in order to define the measurement space $W_{m}$. Since we work in the space $V_{h}$, these can be defined as elements of this space satisfying

[TABLE]

thus resulting in a measurement space $W\subset V_{h}$.

In what follows, we keep our idealized discussion and do not systematically recall that all computations take place in a background discretization space $V_{h}$.

3.2 Algebraic formulation of PBDW

Here we explicitly derive the algebraic formulation of $u_{m,n}^{*}(\omega)$, the function given by the linear PBDW algorithm. Let $X$ and $Y$ be two finite dimensional subspaces of $V$ and let

[TABLE]

be the orthogonal projection into $X$ restricted to $Y$. That is, for any $y\in Y$, $P_{X|Y}(y)$ is the unique element $x\in X$ such that

[TABLE]

Lemma 1

Let ${W_{m}}$ and ${V_{n}}$ be an observation space and a reduced basis of dimension $n\leq m$ such that $\beta({V_{n}},{W_{m}})>0$. Then the linear PBDW algorithm is given by

[TABLE]

with

[TABLE]

Proof 1

By formula (5), $u^{*}_{n}(\omega)$ is a minimizer of

[TABLE]

The last minimization problem is a classical least squares optimization. Any minimizer $v^{*}_{n}\in V_{n}$ satisfies the normal equations

[TABLE]

where $P^{*}_{{W_{m}}|V_{n}}:V_{n}\to{W_{m}}$ is the adjoint operator of $P_{{W_{m}}|V_{n}}$. Note that $P^{*}_{{W_{m}}|V_{n}}$ is well defined since $\beta(V_{n},{W_{m}})=\min_{v\in V_{n}}\|P_{{W_{m}}|V_{n}}v\|/\|v\|>0$, which implies that $P_{{W_{m}}|V_{n}}$ is injective and thus admits an adjoint. Furthermore, since for any $w\in{W_{m}}$ and $v\in V_{n}$, $\langle v,w\rangle=\langle P_{{W_{m}}|V_{n}}v,w\rangle=\langle v,P_{V_{n}|{W_{m}}}w\rangle$, it follows that $P^{*}_{{W_{m}}|V_{n}}=P_{V_{n}|{W_{m}}}$, which finally yields that the unique solution of the least squares problem is

[TABLE]

Therefore $u^{*}_{n}=w+\eta^{*}_{n}=w+v^{*}_{n}-P_{W_{m}}v^{*}_{n}$.

As a direct consequence of this Lemma, we derive the following formulation for the affine PBDW algorithm.

Corollary 2

Let ${W_{m}}$ and ${V_{n}}$ be an observation space and a reduced basis of dimension $n\leq m$ such that $\beta({V_{n}},{W_{m}})>0$. Then the affine PBDW algorithm with respect to a nominal state $\bar{u}$ is given by

[TABLE]

where $\bar{w}=P_{W_{m}}\bar{u}$ and $u^{*}_{m,n}$ is the reconstruction with the linear PBDW method.

Proof 2

The result follows by following the same lines as in the linear case for the shifted minimization problem

[TABLE]

We next derive the symmetric linear system of equations to be solved in order to compute $v^{*}_{n}$ in expression (18). Let $F$ and $H$ be two finite-dimensional spaces of $V$ of dimensions $n$ and $m$ respectively in the Hilbert space $V$ and let $\mathcal{F}=\{f_{i}\}_{i=1}^{n}$ and $\mathcal{H}=\{h_{i}\}_{i=1}^{m}$ be a basis for each subspace respectively. The Gram matrix associated to $\mathcal{F}$ and $\mathcal{H}$ is

[TABLE]

These matrices are useful to express the orthogonal projection $P_{F|H}:H\mapsto F$ in the bases $\mathcal{F}$ and $\mathcal{H}$ in terms of the matrix

[TABLE]

As a consequence, if $\mathcal{V}_{n}=\{v_{i}\}_{i=1}^{n}$ is a basis of the space $V_{n}$ and $\mathcal{W}_{m}=\{\omega_{i}\}_{i=1}^{m}$ is the basis of $W_{m}$ formed by the Riesz representers of the linear functionals $\{\ell_{i}\}_{i=1}^{m}$, the coefficients $\textbf{v}^{*}_{n}$ of the function $v^{*}_{n}$ in the basis $\mathcal{V}_{n}$ are the solution to the normal equations

[TABLE]

where

[TABLE]

since $P^{*}_{{W_{m}}|V_{n}}=P_{{W_{m}}|V_{n}}$ and w is the vector of measurement observations

[TABLE]

Usually $\textbf{v}^{*}_{n}$ is computed with a QR decomposition or any other suitable method. Once $\textbf{v}^{*}_{n}$ is found, the vector of coefficients $\textbf{u}_{n}^{*}$ of $u^{*}_{n}$ easily follows.

3.3 Practical implementation of OMP

Introducing the set $\mathcal{W}=\{P_{{W_{m}}}v\,:\,v\in\mathcal{D}\}$, the recursive step of the OMP algorithm (16) can equivalently written in terms of functions of ${W_{m}}$ as follows. For $n>1$, given $V_{n}=\mathrm{span}\{\varphi_{1},\ldots,\varphi_{n}\}$, we select (see equation (16))

[TABLE]

For the chosen $z_{n+1}$, we take one of the corresponding functions $\varphi_{n+1}$ from $\mathcal{D}$ that satisfy $P_{{W_{m}}}\varphi_{n+1}=z_{n+1}$. We then set $V_{n+1}=\mathrm{span}\{V_{n},\varphi_{n+1}\}$. The computation of $P_{P_{{W_{m}}}V_{n}}\omega$ is a straightforward solving of a linear problem. We look for the coefficients $c=(c_{i})_{i=1}^{n}$ such that $P_{P_{{W_{m}}}(V_{n})}\omega=\sum_{i=1}^{n}c_{i}P_{{W_{m}}}\varphi_{i}$. Since

[TABLE]

it follows that

[TABLE]

which is a linear system of the form

[TABLE]

Denoting $\{\rm w_{1},\dots,w_{m}\}$ an orthonormal basis of ${W_{m}}$, we have

[TABLE]

and

[TABLE]

3.4 Summary of the methods investigated and their computational cost.

Methods:

The methods that we have implemented in our numerical tests and their computational cost are:

Linear PBDW (labelled POD-lin): We build the linear spaces $V_{n}$ from a classical Proper Orthogonal Decomposition of the whole training set $\widetilde{}\mathcal{M}_{\text{train}}$. This is the classical PBDW approach explained in section 2.2.1 and we consider it as the reference that our proposed methods should outperform. 2. 2.

Nonlinear algorithm with manifold partitioning: We partition the training set $\widetilde{}\mathcal{M}_{\text{train}}$ in a number of non-intersecting subsets (see section 2.2.2). For each subset $\widetilde{}\mathcal{M}^{(k)}$ of the partition, we build in an offline phase a reduced model. Two constructions have been tested:

• (a)

P-POD-aff: A Proper Orthogonal Decomposition on each partition $\widetilde{}\mathcal{M}^{(k)}$.

• (b)

P-Greedy-aff: A greedy algorithm: for $n=1$, we set $V^{(k)}_{1}=\operatorname*{span}\{u^{(k)}_{1}\}$ with

[TABLE]

For $n>1$, we select

[TABLE]

and set $V^{(k)}_{n}=\operatorname*{span}\{V^{(k)}_{n-1},u^{(k)}_{n}\}$. 3. 3.

P-DB-aff, Data-driven nonlinear algorithm with manifold partitioning: Since each ultrasound image can be seen as an observation $\omega\in{W_{m}}$, we run the OMP algorithm of section 2.2.2 to build $V_{n}(\omega)$ and do the reconstruction. Note that the greedy search has to be done online since we need the knowledge of the measurement. To speed-up computations, instead of searching in the whole training set

[TABLE]

we restrict the search to the partition

[TABLE]

Computational costs:

We outline the computational cost of the methods in Table 1. Note that all of them can be decomposed into an online and an offline phase. The cost of generating the snapshots database is of the same order in all methods, namely $\mathcal{O}(\mathcal{N}^{\alpha}\#\widetilde{}\mathcal{M}_{\text{train}})$ where $\alpha$ is the scaling behavior of the linear solver (usually $\alpha=2$ for iterative methods). The cost of building the reduced model in POD-lin and P-POD-aff is the one to compute the SVD of the correlation matrix of the snapshots. For POD-lin, the matrix is of size $\#\widetilde{}\mathcal{M}_{\text{train}}^{2}$. For P-POD-aff, we have to compute the SVD of $K$ matrices of size $(\#\widetilde{}\mathcal{M}^{(k)})^{2}$ for $k=1\dots,K$. If these SVD computations are done in parallel, we can obtain important time reductions compared to POD-lin since in general $\#\widetilde{}\mathcal{M}^{(k)}\ll\#\widetilde{}\mathcal{M}_{\text{train}}$. Regarding the cost of the online phase, if we store the QR decomposition of normal equations (20) (at a cost of $\mathcal{O}(n^{3})$ operations, not reported on the table), the first three methods only need to solve the associated linear system, which costs $\mathcal{O}(n^{2})$ operations. In the case of P-DB-aff, we additionally have to run a greedy algorithm in $\mathbb{R}^{m}$ as explained in section 3.3, hence the additional cost reported in the table. Note however that the time required for this extra step is not significantly slowing down computations thanks to the fact that we work in $\mathbb{R}^{m}$ and $m$ is usually moderate. This is in contrast to the greedy algorithm of the offline phase of P-Greedy-aff, which takes place in $\mathbb{R}^{\mathcal{N}}$.

4 Application: Reconstruction of 3D blood velocity fields from Doppler ultrasound images

We apply the above described methodology to reconstruct a 3D blood velocity field on a human carotid artery from Doppler ultrasound images. The images are synthetically generated and the use of data from real patients is deferred to a future work. The main goal of the tests is twofold:

The first goal is to compare the different strategies to construct the space $V_{n}$ and the recovery algorithms. The method that can be consider a sort of a baseline to be compared to is the POD. 2. 2.

Second, we aim at assessing the ability of these recovery strategies to estimate the velocity in a semi-realistic idealised setting. Methods failing in achieving a $10\%$ accuracy on the peak velocity in this setting should not be used in more realistic scenarios.

The section is organized as follows. First, we present the parameter-dependent model that will define the manifold $\mathcal{M}$ on which we will rely to compute different reduced models. Second, we explain how to define a measurement space ${W_{m}}$ from a Doppler velocity image. Finally, we present results on the comparison of the different methods based on their reconstruction performance.

4.1 The model: incompressible Navier-Stokes equations

Let $\Omega$ be a spatial bounded domain of $\mathbb{R}^{3}$ with the shape of a human carotid artery as given in Figure 1(a). The boundary $\Gamma\coloneqq\partial\Omega$ is the union of the inlet part $\Gamma_{i}$ where the blood in entering the domain, the outlets $\Gamma_{o,1}$ and $\Gamma_{o,2}$ where the blood is exiting the domain after a bifurcation, and the walls $\Gamma_{w}$.

We consider the following incompressible Navier-Stokes equations on $\Omega$ and over the time interval $[0,T]$ for $T>0$. For a fluid with density $\rho\in\mathbb{R}^{+}$ and dynamic viscosity $\mu\in\mathbb{R}^{+}$, we search for all $t\in[0,T]$ the couple $(u(t),p(t))\in[H^{1}(\Omega)]^{3}\times L^{2}(\Omega)$ of velocity and pressure such that

[TABLE]

These equations are closed by adding a zero initial condition and the following boundary conditions:

Boundary conditions:

• (a)

No-slip condition for the vessel wall, that is, $u=(0,0,0)^{T}$ on $\Gamma_{w}$.

• (b)

The inlet boundary $\Gamma_{i}$ lies in the $xz$ plane and we apply a Dirichlet condition $u=[0,u_{\text{in}},0]^{T}$. The component b is a function $u_{\text{in}}(t,x,z)=u_{0}~{}g(t)f(x,z)$, where:

• i.

$u_{0}\in\mathbb{R}^{+}$ is an scaling factor. The function $g(t)$ is built by interpolating experimental flow data in the common carotid area taken in [54]. Its behavior is given in Figure 1(b).

• ii.

The function $f$ is a 2D logit-normal distribution

[TABLE]

where the parameter, $s\in\mathbb{R}^{+}$, controls the axial symmetry of the inlet flow.

• (c)

For the outlet boundaries $\Gamma_{o,1}$ and $\Gamma_{o,2}$, we use the so-called Windkessel model (see [33]), which gives the average pressure over each $\Gamma_{o,k}$,

[TABLE]

where $p_{d}^{k}\in\mathbb{R}$ is called distal pressure and is the solution to the ordinary differential equation:

[TABLE]

This model aims to represent the cardiovascular system behavior beyond the boundaries of the working domain with a minimal increase in the computational cost. It is based on an analogy between flow and pressure with current and voltage in electricity. This is the reason why $C_{d}^{k}$ is called distant capacitance and $R_{p}^{k}$ and $R_{d}^{k}$ are respectively called proximal and distant resistances. These three parameters are positive real numbers.

For the numerical solution of (25), we use FeLiSCe, a finite elements method (FEM) software developed at Inria. Spatial discretization of the carotid geometry leads to a tetrahedron mesh of 42659 vertices and 225196 tetrahedra. The meshing process is in charge of MMG [55]. Time is discretized with a semi-implicit backward Euler scheme with time-step $\delta t=2\cdot 10^{-3}s$, which means that the convective term in the Navier-Stokes equations is written with the velocity explicitly and the velocity gradient implicitly in order to circumvent the non-linearity of the system. An explicit scheme is used to numerically solve the ODE on the distal pressure in the Windkessel model. In addition, a backflow stabilization is added in order to address potential instabilities in the outlet boundaries (see, e.g., [56]). Data visualization is done by using ParaView [57] and Vizir [58].

The coupled problem for velocity and pressure is discretized using $\mathbb{P}_{1}-\mathbb{P}_{1}$ Lagrange elements, with a monolithic approach. In order to avoid the inf-sup constraint imposed by the saddle point nature of the problem we use the Brezzi-Pitkäranta stabilization approach (for further details, we refer to [59] and [60]), that basically perturbs the weak form of the mass conservation equation with a stiffness-like term that scales with the square of the elements size. Each forward simulation has $\mathcal{N}=127977$ degrees of freedom at each time step. The linear systems are solved by means of a GMRES [61] method with additive Schwarz preconditioning.

Remark 2

Note that one could use more sophisticated models involving, for instance, fluid-structure interactions or more refined Winkessel models for the pressure. Our present model is a trade-off between its degree of realism and the difficulty and time to solve it. We refer to [33] for a detailed overview of cardiovascular modeling.

Let us define the manifold of solutions that we consider in our numerical experiments. We set the following coefficients to a fixed value

[TABLE]

The ratio of the distal resistances for the Windkessel model at the outlets is introduced:

[TABLE]

This parameter plays an important role since it impacts on how the blood flow splits between the two branches. When $\eta\to 0$ or $\infty$, one branch is obstructed and the blood tends to flow through the other branch. In the following, we call this situation an arterial blockage. When $\eta\approx 1$, the flow splits more or less equally and there is no blockage.

We define the heart rate as the number of cardiac cycles per minute, that is,

[TABLE]

where $T_{c}>0$ is the cardiac cycle duration expressed in seconds. We have $T_{c}=T_{sys}+T_{dia}$, where $T_{sys}$ and $T_{dia}$ are the duration of the systole and diastole respectively.

The manifold $\mathcal{M}$ is generated by the variations of the following six parameters

[TABLE]

Note that the time $t$ is also seen as a parameter. The parameter set is thus

[TABLE]

and the manifold of solutions is

[TABLE]

At this point, several comments are in order:

Note that we only consider velocity fields aince in the present work we are only concerned by the reconstruction of the blood flow velocity. The reconstruction of other quantities of interest, such as the pressure, will be addressed in a forthcoming work.

• 2.

For each $y\in\mathrm{Y}$, the velocity $u(y)$ is a function of $[H^{1}(\Omega)]^{3}$. In the following, we will view it as a function from

[TABLE]

which, endowed with the inner product,

[TABLE]

defines a Hilbert space.

• 3.

Since the time variable has been included as a parameter, a simple way to build nonlinear reduced models is to set a window parameter $\tau>0$ and consider the subset $\mathcal{M}^{(k)}=\mathcal{M}_{[t_{k}-\tau,t_{k}+\tau]}\subset\mathcal{M}$, where $t$ is restricted to the interval $[t_{k}-\tau,t_{k}+\tau]$ of size $2\tau$ centered around a given time $t_{k}$. We can then build reduced models to reconstruct this specific time interval. As we will see in the numerical experiments, this strategy is very effective in our problem because the velocity presents two regimes given by the systole and diastole periods.

• 4.

The computation of reduced models involves a discrete training subset $\widetilde{}\mathcal{M}_{\text{train}}\subset\mathcal{M}_{h}$ which, in the experiments below, involves $\#\widetilde{}\mathcal{M}_{\text{train}}=78528$ snapshots $u(y)$. The parameters are chosen from a uniform random distribution and we only save the solutions during the second cardiac cycle of each simulation.

• 5.

For the purposes of illustrating the potential of the method for diagnoses, we introduce a notion of sickness in terms of the arterial blockage in the following way.

Definition 1** (Sick patient)**

We say that the output of the simulation corresponds to a healthy patient when $\eta\in[0.5,1.5]$. Outside of this range, simulations correspond to sick patients.

We thus have the partition

[TABLE]

with $\mathcal{M}_{\text{healthy}}\coloneqq\mathcal{M}_{\eta\in[0.5,1.5]}$ and $\mathcal{M}_{\text{sick}}=\mathcal{M}_{\eta\in[0.05,0.2]\cup[5,20]}$.

4.2 Measurements

In a regular partition of $[0,T]$, we are given Doppler ultrasound images that contain information on the blood velocity on a subdomain of the carotid. For typical ultrasound machines with acquisition time of 0.1 milliseconds, a CFI is built with, for instance, 32 consecutive B-mode frames, which leads to a measurement sampling of 32 milliseconds. From the image, the observations $\ell_{i}(u)$ are extracted and used to build a complete time-dependent 3D reconstruction of the blood velocity in the whole carotid $\Omega$.

Depending on the technology of the ultrasound device, there are two different types of velocity images. In most cases, ultrasound machines give a scalar mapping which is the projection of the velocity along the direction $b$ of the ultrasound probe. This mapping is called color flow image (CFI, see Figure 2(a)). In more modern devices, it is possible to get a 2D vector flow image corresponding to the projection of the velocity into the plane. This mapping is called vector flow image (VFI, see Figure 2(b)). For both imaging modes, the velocity is estimated by time averaging techniques (we refer to [62] and [63] for further details).

In the following, we work with an idealized version of CFI images. For each time $t$, a given image is a local average in space of the velocity projected into the direction in which the ultrasound probe is steered. More specifically, we consider a partition of $\Omega=\cup_{i=1}^{m}\Omega_{i}$ into $m$ disjoint subdomains (voxels) $\Omega_{i}$. Then, from each CFI image we collect

[TABLE]

where $b\in\mathbb{R}^{3}$ is a unitary vector giving the direction of the ultrasound beam. From (30), it follows that the Riesz representers of the $\ell_{i}$ in $V$ are simply

[TABLE]

where $\chi_{i}$ denotes the characteristic function of the set $\Omega_{i}$. Thus the measurement space is

[TABLE]

Since the voxels $\Omega_{i}$ are disjoint from each other, the functions $\{\omega_{i}\}_{i=1}^{m}$ are orthogonal and therefore having a CFI image is equivalent to having

[TABLE]

Remark 3

The case of VFI images can be treated similarly. This imaging mode gives $2m$ measurements

[TABLE]

and

[TABLE]

where $b$ is the unitary vector giving the direction of the ultrasound beam and $b_{\perp}$ is the unitary vector perpendicular to it and contained in the image plane. Therefore, ${W_{m}}^{(VFI)}={W_{m}}^{(CFI)}\oplus\operatorname*{span}\{\chi_{\Omega_{i}}b_{\perp}\}_{i=1}^{m}$ which is a space of dimension $2m$. This clearly shows that the additional direction $b_{\perp}$ enriches the quality of the measurements in the sense that for any $u\in V$, the approximation error $\|u-P_{{W_{m}}}u\|$ will be smaller with the VFI mode than with the CFI one. As a result, the CFI mode which we consider in our examples is a more challenging case since the measurements contain less information.

4.3 Reconstruction on a first example with healthy patients

To validate our method, we first consider a simple example where we only work with healthy patients, so the manifold is $\mathcal{M}_{\text{healthy}}$. One color flow image contains information of the velocity averaged over 552 voxels. As explained in section 4.2, this sets the dimension of the observations space to $m=552$ (whereas for VFI, $m=1104$). As a result, each image can be seen as an observation $\omega\in{W_{m}}$ (see Figure 3). The dimension $m=552$ may seem quite large but it is representative of the one provided by modern pulsated ultrasound devices. 297 healthy patients are simulated in order to build the training set. For each one of them a number of snapshots containing the second cardiac cycle is stored so there are multiple snaptshots per patient due to time marching. This leads to a training set composed of $\#\widetilde{}\mathcal{M}_{\text{train}}=56383$ snapshots. The performance of the algorithms is tested on a test set $\widetilde{}\mathcal{M}_{\text{test}}$ of 32 healthy patients each one with a parameter configuration in range but different from those inside $\widetilde{}\mathcal{M}_{\text{train}}$.

For the comparison, we implement the algorithms outlined in section 3.4. We recall them below and give specific details related to the current test case:

Linear PBDW: The space $V_{n}$ is built from a classical POD of the whole training set $\widetilde{}\mathcal{M}_{\text{train}}$. 2. 2.

Nonlinear algorithm with manifold partitioning: In real medical examinations, the heart rate HR of the patient and the time $t$ in which the ultrasound image is taken are known. We exploit this fact to decompose $\widetilde{}\mathcal{M}_{\text{train}}$ into $K=IJ$ subsets

[TABLE]

where, for each $k=(i,j)$ in $\{1,\dots,I\}\times\{1,\dots,J\}$,

[TABLE]

Among the possibilities to build the manifold partitioning, our approach is driven by the reconstruction error of the elements inside $\widetilde{}\mathcal{M}_{\text{test}}$. That is to say, we have computed the reconstruction $u^{*}_{i}$ of $u_{i}$ for $i=1,\ldots,\#\widetilde{}\mathcal{M}_{\text{test}}$ for a set of couples $\left(\delta_{\text{HR}},\tau\right)$ with a fixed dimension for the reduced model ($n=30$). Hence, we can pick the optimal values in this sense:

[TABLE]

which has led in our context to the sizes $\tau=T_{c}/10$ and $\delta_{\text{HR}}=5$ beats per minute.

For each subset $\widetilde{}\mathcal{M}^{(k)}$, a reduced space is built. Two constructions have been tested:

• (a)

POD of $\widetilde{}\mathcal{M}^{(k)}$. We call this approach P-POD-aff, Partitioned-POD-affine, in the plots.

• (b)

A greedy algorithm. We call this approach P-Greedy-aff, Partition-Greedy-affine, in the plots

During the online reconstruction, given $t$ and HR, we select the appropriate subset $\widetilde{}\mathcal{M}^{(k)}$ that includes $t$ and HR and reconstruct with a linear PBDW with the reduced model corresponding to $\widetilde{}\mathcal{M}^{(k)}$. 3. 3.

Data-driven nonlinear algorithm with manifold partitioning: This approach is labeled in the plots as P-DB-aff, Partitioned-Data-Based-affine. Since each CFI image can be seen as an observation $\omega\in{W_{m}}$, we run the Orthogonal Matching Pursuit algorithm of section 2.2.2 to build $V_{n}(\omega)$ and do the reconstruction. Note that the greedy search has to be done online since we need the knowledge of the measurement.

For each state $u\in\widetilde{}\mathcal{M}_{\text{test}}$, we compute the relative error:

[TABLE]

where $A_{m,n}(P_{{W_{m}}}u)$ denotes any of the above four reconstruction algorithms. A patient in the training set is represented by a sequence of simulated states during the second cardiac cycle. This raises the interest in evaluating the reconstruction quality by looking at the following relative error in time:

[TABLE]

where we normalize by the total energy in the cardiac cycle $\left(\int_{T_{c}}^{2T_{c}}\|u(t)\|^{2}\mathrm{d}t\right)^{1/2}$.

Figures 4(a) and 4(b) give the average and worst case performance of the four methods,

[TABLE]

Note that as $n$ increases, the error decreases for all methods, except perhaps for the P-DB-aff  approach where the error tends to stagnate for large values of $n$. This could be due to the fact that P-DB-aff heavily relies on the measurement information, which, in the present application, might not deliver enough to learn reduced models $V_{n}(\omega)$ that improve the accuracy as $n$ grows. For instance, the Doppler measurements are close to zero in the diastole of the cardiac cycle so the information that they provide may be insufficient to build a good $V_{n}(\omega)$. We also see from the figure that the nonlinear method P-POD-aff outperforms the rest in the sense that it delivers a given target accuracy with a smaller dimension $n$ of the reduced model. For instance, if we fix a target accuracy on the average performance to $10^{-2}$, we see that the POD-lin requires 40 modes to achieve it, P-Greedy-aff requires 20, P-DB-aff requires 17 and P-POD-aff requires only 10 (see Figure 4(a)).

We next fix $n=30$ and study the error in time $e(t,u,A_{m,30})$ on Figure 5. We observe that the reconstruction tends to be better during the late diastole phase of the cardiac cycle.

As discussed in section 2.2, the inf-sup constant $\beta({V_{n}},{W_{m}})$ might yield to stability issues when $n\to m$ since its value tends to zero (see equations (8) and (9)). Figures 7(a) and 7(b) show its behavior for the four methods during the systole and diastole period. We observe that the four methods perform similarly in terms of stability for the peak systole reconstruction. For the diastole phase, we observe that the inf-sup constant in P-DB-aff performs slightly worse than the rest. We think that this could be due to the fact that the measurement space ${W_{m}}$ is not rich enough to allow P-DB-aff to properly learn reduced models when $n$ becomes large.

4.4 Application to arterial blockage detection

In this example we illustrate that even when the Doppler images do not give information on the whole carotid, we can nevertheless reconstruct the velocity field in the whole domain with our methodology. This is important for actual practice since doctors do not have images in the whole carotid due to morphological constraints. We also show in our example that the method has potential to efficiently estimate in real time relevant quantities of interest.

We illustrate these ideas in the following example: we consider the same setting as before but now the Doppler image does not provide information about the flow in the carotid bifurcations. Therefore, the image does not see the flow split in the common carotid downstream (see Figure 13(a)). In this example we have tested the impact of working with CFI ($m=233$) or VFI images ($m=466$).

We train our reconstruction methods on a set $\widetilde{}\mathcal{M}_{\text{train}}$ containing sick and healthy patients. Here, we only work with two of our previous nonlinear algorithms:

P-POD-aff 2. 2.

P-DB-aff

Figure 8 shows the average and worst case errors

[TABLE]

as a function of the dimension $n$ of $V_{n}$. Like in the previous example, both methods are delivering a very satisfactory accuracy: the average error is below $5\cdot 10^{-2}$ for both methods for all values of $n$. The method P-POD-aff consisting in a partition of the manifold outperforms the data-based one, P-DB-aff.

We next show that the method is efficient to assist in the detection of arterial blockages that may cause severe health problems like a stroke. Since a blockage alters the distribution of the velocity field after the bifurcation, a quantity of interest that could serve as a clinical index is the ratio

[TABLE]

where

[TABLE]

is the blood flow at the outlet $\Gamma_{o}^{i}$, $i=1,\,2$, and $t_{peak}$ is the peak systole instant. Figure 9 shows the evolution of $Q_{i}(t)$ in time for a sick patient and its approximation with our two reconstruction methods. We observe that, regardless of the image format (CFI or VFI), both methods deliver very satisfactory predictions of the flow.

In Figure 10, we compare the value of the exact ratio $r$ with the reconstructed one for sick and healthy patients $u\in\widetilde{}\mathcal{M}_{\text{test}}$.

To define a threshold ratio $r^{*}$ to decide whether the patient has arterial blockage or not, we can take the average of the flow ratios between the healthier of the patients in the sick group and the sicker of the patients in the healthy group, namely,

[TABLE]

where $r(u)$ denotes the flow ratio associated to the velocity field $u$, as defined in (36). In our data-base, we obtain, $r^{*}=1.25$, so any patient for which $r>1.25$ will be considered as presenting high blockage risk. Note that the approximation is very close to the real value for moderate values of $r$ regardless of the image modality. However, we tend to overestimate the value for $r>1.7$. In presence of a blockage, $r$ becomes significantly larger than one so the overestimation is by far more preferable than an underestimation. Indeed, the overestimation makes our method conservative and, in the worst case, we will conclude with a false positive. However, the method will not lead to a false negative diagnosis, which would leave a sick patient without treatment/surgery.

To illustrate the reconstruction for both CFI and VFI, let us consider the forward simulation of a peak systole blood flow as depicted in Figure 11. Assuming this field as a ground truth, we can see the reconstruction with P-POD-aff and P-DB-affin Figures 12 and 13.

5 Conclusions and perspectives

We have investigated the construction of state estimation techniques with different reduced models. We have compared their accuracy for the problem of reconstructing the blood flow velocity field for which partial information is given by Doppler ultrasound images. The numerical tests were synthetically generated and mimic a real context in medical applications. The results show that the classical linear PBDW method with the classical POD reduced model can be outperformed by other models that are built more specifically for the reconstruction task. In particular, the nonlinear reconstruction method built by manifold partitioning and a local POD basis on each partition outperforms all the other bases choices. It presents a good trade-off between simplicity and efficiency and its accuracy is in some cases even ten times better than the classical PBDW. The data-driven approach gave less performant results than the one with local POD bases but superior to the classical approach. We think that the specific application may have had an important impact on the observed performance. Since we build a space $V_{n}(\omega)$ from the Doppler measurements $\omega$, they may not deliver enough information to learn a good reduced model, especially in the diastole phase where $\omega$ is close to zero. As a result, it seems important to investigate further the performance of this method in other reconstruction problems in order to gain better knowledge on its range of applicability.

Acknowledgments

This research was partially supported by the Emergence Project “Models and Measurements” of the Paris City Council. Also, authors thankfully acknowledge the financial support of the ANID Ph. D. Scholarship 72180473.