# POD method for acceleration of blood flow reconstruction in a vessel with contrast enhanced X-ray CT

S. Huang, Monica Sigovan, Bruno Sixou

<figure><img src="image_1.png" /></figure>

## To cite this version:

S. Huang, Monica Sigovan, Bruno Sixou. POD method for acceleration of blood flow reconstruction in a vessel with contrast enhanced X-ray CT. Computer Methods in Biomechanics and Biomedical Engineering: Imaging & Visualization, 2022, 11 (4), pp.1478-1490. 10.1080/21681163.2022.2146316. hal-04177980

HAL Id: hal-04177980

https://hal.science/hal-04177980

Submitted on 7 Aug 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
POD method for acceleration of blood flow reconstruction in a vessel
with contrast enhanced X-ray CT.

S.Huangª, M.Sigovan ª and B.Sixou ª

ª CREATIS, CNRS UMR5220; Inserm U630; INSA-Lyon; Université Lyon 1; Université de Lyon, F-69621 Villeurbanne Cedex, France.

ARTICLE HISTORY

Compiled August 6, 2023

ABSTRACT

The reconstruction of blood velocity in a vessel from contrast enhanced X-ray CT projections is a complex inverse problem. In this work, we investigate an inversion approach based on the adjoint method with a partial differential equation as constraint. The proper orthogonal decomposition method is used to accelerate the reconstruction with reduced basis for the velocity, the density and the adjoint method. The effectiveness of the reconstruction method is shown with simulated realistic stationary and non-stationary blood flows in a vessel.

KEYWORDS

Tomography, inverse problems, hemodynamics, proper orthogonal
decomposition(POD).

1. Introduction

The prognosis of partially blocked vessels is an extremely important step in the treat-
ment planning for patients with cardiovascular diseases. It relies on the understanding
of blood flow in the cardiovascular system which is a major challenge ( Shi 2011;
Crosetto 2011; Quarteroni 2017). The use of measurements to reliably infer fluid ve-
locity and pressure/stress fields is not a straightforward task and quantifying the blood
flow in the vascular network remains very difficult. Magnetic resonance and ultrasound
allow the measurement of blood flow velocity but with a low spatial resolution ( Pelc
1994; MacDonald 2015; Markl 2016). There are very few studies about the measure-
ment of blood velocity with X-ray CT ( Prevrhal 2011; Korporaal 2016; Bouillot 2019;
Barfett et al. 2014; Daly 2013). A time of flight study is presented by Barfett et al.
(Barfett et al. 2014). A method to use the change of the density gradients to evalu-
ate the velocity field is investigated in (Korporaal 2016). The spatial evolution of the
time-intensity curves is examined in (Bouillot 2019). Recently, we have proposed a
new method to reconstruct the velocity field in a vessel with contrast enhanced spec-
tral CT (Huang 2021). It is estimated using the measurement of 2D Radon projections
perpendicular to the main vessel axis direction. We assume that the propagation of the
tracer in the vessel can be modeled by a partial differential equation (PDE) based on
the unknown velocity field. We leverage this tranport equation and the PDE is used as

CONTACT B.Sixou. Email: bruno.sixou@insa-lyon.fr
a constraint in the optimization algorithm which relies on the adjoint method (Herzog 2010; Plessix 2006). Using a comparison with ground truth quantities, simulations experiments have shown that good reconstructions can be obtained with a better spatial resolution and that all the components of the velocity can be estimated with a good accuracy. Yet, the computational time remains very high and from a computational point of view, the challenge remains to reach a good resolution for the velocity field within acceptable computation cost.

It is thus necessary to derive a computationally efficient surrogate model to reduce these computational costs since we have to evaluate repeatedly the output of expensive forward and adjoint partial differential equations to solve the inverse problem. With reduced-order modeling, the structure of the physical problem and the governing equations are taken into account. The reduced basis method represents PDE solutions in terms of known basis functions that contain information about the solution structure (Quarteroni 2015). Model reduction can be obtained with Proper Orthogonal Decomposition (POD) combined with Galerkin projection for time dependent and parametrized problems (Ghattas 2021; Volkwein 2013). The behaviour of the full order simulations can be recovered with a large reduction of the computational time. The principle of the method is to obtain a low-dimensional approximation with projection on a low-dimensional subspace. POD reduced order methods are among the most successful tools for studying fluid flows (Rowley 2005; Du 2013; Lu 2019). The reduced model can be used over a range of physical or geometrical parameter values which is crucial for inverse problems, and uncertainty quantification (Ghattas 2021; Volkwein 2013). This approximation methodology has been considered for parametrized optimal control problems and inverse problems governed by parametrized partial differential equations that play a major role in several applications, like biomedical applications (Ballarin 2017; Ravindran 2000).

In this work, inspired by these approaches and in order to improve the convergence properties of our adjoint method to solve the inverse problem of blood velocity reconstruction, we investigate the approximation of the solution with proper orthogonal decomposition methods. The quantitites of interest, the velocity field, the density and the adjoints variables are expanded in terms of appropriate sets of basis functions obtained from snapshots. In order to take into account the uncertainty of several parameters of the forward model, the approach is applied for parametrically varying conditions. The efficiency of the framework is demonstrated with numerical examples with stationary and non-stationary velocity fields with realistic simulations. We also present comparisons with full order approximation simulations.

The paper is organized as follows. After the introduction, in section 2, the inverse problem will be presented and formulated as an optimal problem. In section 3, we will present our approach to approximate the unknown velocity field combining the POD decompositions applied to the velocity, density and adjoint variables. In section 4, the simulation results achieved with the POD reduced method to solve the inverse problem are detailed and discussed.

## 2. Inverse problem formulation

In this section, we present the inverse problem set up and the solution methodology of our former work and we extend it to parametrically varying conditions. The inverse problem considered in the following is parametrized by a parameter $\mu$ belonging to some parameter space $\mathcal{M}$ and we reformulate the inverse problem set up and solution
methodology with this new parameter. The velocity reconstruction of the blood flow can be considered as an inverse problem with a partial differential equation constraint. The inverse problem can be formulated as an optimal control problem. The new inversion method was presented in detail in (Huang 2021). The main flow direction is the z axis and two-dimensional Radon projections are acquired perpendicularly to this axis, in the xy plane. The space-time box considered is denoted $Q = \Omega \times [0, T]$ where $\Omega$ is a bounded spatial domain and $[0, T]$ the time interval. We will consider the Hilbert space $L_2([0, T], L_2(\Omega))$ of square integrable functions, $t \to f(\mathbf{x}, t)$, $t \in [0, T]$, with values in the Hilbert space $L_2(\Omega)$ of square integrable functions on $\Omega$, embedded with the norm:

$$ \|f\|_{\mathcal{L}} = \left( \int_{0}^{T} \|f\|_{L_2(\Omega)}^2 dt \right)^{1/2}. \qquad (1) $$

The propagation of the contrast agent with the density $f(\mathbf{x}, t)$ is modeled with a partial differential equation based on a realistic velocity field obtained with the Navier-Stokes equation (Glowinski 1992; Galdi 1994; Temam 2001). We assume that the density of the other materials can be neglected and the evolution of the Radon projections with time is related to its propagation of the tracer. The propagation of the contrast agent in the vessel is modeled by a convection-diffusion equation with a diffusion parameter $D$ and the velocity $V = (u, v, w) : [0, T] \times \Omega \to \mathbb{R}^3$. In this work, stationary and non-stationary velocity field are considered, obtained with the stationary and non-stationary Navier-Stokes equations. The velocity field $V$ is parametrized the first component $\mu_1$ of $\mu$ related to the boundary/initial conditions. The parameter $\mu_1$ will be specified in the following. The contrast agent enters the vessel through a disk ($\mathcal{D}$) in the plane $z = 0$. We assume that the value of the function $f$ on the disk ($\mathcal{D}$) is associated to the parameter $\mu_2$ of the parametrized inverse problem. The convection-diffusion equation with initial and boundary conditions can be written as (Huang 2021):

$$ \left\{ \begin{array}{l} e(f, \mathbf{V}, \mu) = \frac{\partial f(\mathbf{x},t)}{\partial t} + \mathbf{V} \cdot \nabla f(\mathbf{x}, t) - D \Delta f(\mathbf{x}, t) = 0 \\ f(\mathbf{x}, t) = \mu_2 \quad \forall x \in (\mathcal{D}) \quad \forall t \in [0, T] \\ f(\mathbf{x}, 0) = 0 \quad \forall x \in \Omega - (\mathcal{D}) \end{array} \right. \qquad (2) $$

where the convective term is given by:

$$ \mathbf{V} \cdot \nabla f(\mathbf{x}, t) = u \frac{\partial f(\mathbf{x}, t)}{\partial x} + v \frac{\partial f(\mathbf{x}, t)}{\partial y} + w \frac{\partial f(\mathbf{x}, t)}{\partial z} \qquad (3) $$

and where $\mu_2$ is the second component of $\mu$. It should be noted that the dynamics of the scalar field $f$ is thus considered as a prior knowledge for our inverse problem.

Let $\Sigma \in \mathbb{R}^2$ be a bounded Lipschitz domain, the Radon transform of the density of the tracer can be written as:

$$ Rf(\theta, s) = \int_{\Sigma \cap L(\theta, s)} f(\mathbf{x}) dl(\mathbf{x}) \qquad (4) $$

As displayed on Figure1, the line $L(\theta, s)$ is determined by the angle $\theta \in [0, \pi]$ and by the distance $s \in [-a, a]$ from the origin (Natterer 1986). The scanner is rotating
<figure><img src="image_3.png" /><figcaption>Figure 1. Principle of integral Radon transform.</figcaption></figure>

with time and the Radon projections are acquired for $N_{\theta}$ projection angles regularly distributed in the angular range $\Delta\theta(t)$ which is changing with time. For each measurement angle, $N_p$ projections are measured and the 2D Radon projections are obtained for $N_z$ values along the z axis, with $z \in [0, z_{max}]$. Generally speaking, the aim of the inverse problem is to evaluate the density $f(\mathbf{x}, t)$ from the linear equations $Rf(\theta, s, z) = p^{\delta}(\theta, s, z)$, where $p^{\delta}$ are the measured projections with noise, s and z are coordinates of the two-dimensional detector. A variational approach is investigated based on a cost functional written as:

$$J(f, \mathbf{V}, \mu) = \frac{1}{2} \int_{0}^{T} \|Rf(\theta, s, z) - p^{\delta}(t)\|_{2}^{2} dt + \frac{\alpha_{s}}{2} \|\nabla \mathbf{V}\|_{2}^{2} \quad (5)$$

with $\mu = (\mu_1, \mu_2)$. The first term is a data term with a $L_2$ norm defined on $[0, \pi] \times [-a, a] \times [0, z_{max}]$. The second one is a spatial regularization term. In the regularization functional, $\|\nabla \mathbf{V}\|_2$ is based on the Frobenius norm of the Jacobian matrix of the velocity $\mathbf{V}$ combined with the $L_2(\Omega)$ norm of the square of each entry of this matrix. The parameters $\alpha_s$ is a regularization parameter.

The constraint $e(f(\mathbf{V}), \mathbf{V}, \mu)$ is replaced by an operator $f(\mathbf{V}, \mu)$ mapping the velocity field to the density. The inverse problem can be rewritten with a reduced cost functional $\hat{J}$ depending on the velocity $\mathbf{V}$ and on $\mu$ and formulated as an optimal control problem:

$$\min_{\mathbf{V}, f \in U \times Y} J(f, \mathbf{V}, \mu) = \min_{\mathbf{V} \in U} J(f(\mathbf{V}), \mathbf{V}, \mu) = \min_{\mathbf{V} \in U} \hat{J}(\mathbf{V}, \mu) \quad (6)$$

In order to obtain the gradient of the reduced functional, we have to consider a Lagrangian which combines the regularization functional with the weak form of the partial differential equation. The Lagrangian $L$ is obtained with a dual variables $p$:
$[0, T] \to L_2(\Omega)$ and written as:

$$L(f, \mathbf{V}, p, \mu) = J(f, \mathbf{V}, \mu) + \beta \int_{0}^{T} <e(f, \mathbf{V}, \mu), p(t)> dt \quad (7)$$

where $\beta$ is a Lagrangian parameter and $<.,>.$ denotes the $L_2(\Omega)$ scalar product.

This inverse problem for a given value of $\mu$ has been studied in (Huang 2021). The minimizer can be obtained with the adjoint method by finding the solution of a system of three equations. The first equation is the convection-diffusion Eq.(2). As a necessary condition, we obtain the following equation for the adjoint variable $p$:

$$\beta\left(\frac{\partial p}{\partial t} + \mathbf{V} \cdot \nabla p + D \nabla \cdot \nabla p\right) = R^t(Rf - g) \quad p(\cdot, T) = 0 \quad (8)$$

This adjoint equation is a backward equation and $R^t$ is the transpose of the Radon operator. The gradient of the reduced functional is obtained with the variations of the Lagrangian with respect to the velocity field. It is given by:

$$\nabla \hat{J} = \int_{0}^{T} (\beta p \nabla f - \alpha_s \Delta \mathbf{V}) dt \quad (9)$$

The algorithmic scheme and the numerical implementation of the adjoint method are detailed in (Huang 2021) coupling Eq.(2), Eq.(8), Eq.(9). It has been shown that the inverse problem is well-posed for any value of $\mu$. The proof is based on the implicit function theorem, on the fact that $J$ and $e$ are Fréchet differentiable, and that for every value $\mu \in \mathcal{M}$ and for every control variable $\mathbf{V}$, there exists an unique solution $f(\mathbf{V}, \mu)$. The reduced order method proposed in this work is based on the assumption that the solution manifold

$$\mathcal{M} = \{f(\mu), \mathbf{V}(\mu), p(\mu); \mu \in \mathcal{M}\} \quad (10)$$

can be approximated by a low dimensional subspace in an offline phase that may be computationally expensive. The set $\mathcal{M}$ is sampled with a discrete subset of parameters $\mathcal{M}_s \subset \mathcal{M}$. This finite dimensional parametric set has cardinality $N_s = |\mathcal{M}_s|$. A solution manifold can be defined for the snapshots, $\{f(\mu), \mathbf{V}(\mu), p(\mu); \mu \in \mathcal{M}_s\}$. If the finite parameter space $\mathcal{M}_s$ is large enough, this manifold is a good representation of the manifold $\{f(\mu), \mathbf{V}(\mu), p(\mu); \mu \in \mathcal{M}_s\}$. The full order discretized finite element solutions will be denoted as $\{f_{\mathcal{F}}(\mu), \mathbf{V}_{\mathcal{F}}(\mu), p_{\mathcal{F}}(\mu), \mu \in \mathcal{M}_s\}$. After finite element discretization, they are vectors of $\mathbb{R}^m$ or $\mathbb{R}^{3m}$ for the density, adjoint and velocity fields respectively, where $m$ is the number of coefficients in the full order FE approximation. They will be projected on low dimensional spaces to obtain the reduced basis $(\Psi_{f,i})_{1 \le i \le N_f}$, $(\Psi_{V,j})_{1 \le j \le N_v}$, $(\Psi_{p,k})_{1 \le k \le N_p}$ with dimensions $N_f, N_v$ and $N_p$.

### 3. Inversion method with proper orthogonal decompositions

In subsection 3.1, we give a description of the POD method with the construction of the basis functions from the collection of evaluations of the velocity, density and adjoint variables. We continue in subsection 3.2 with the methodology to solve the inverse problem.
### 3.1. The offline step:building POD basis of density, velocity and adjoint variables

In this work, we have performed separate PODs on each of the solution manifolds {$f_F(\mu_1), \dots, f_F(\mu_{N_s})$}, {$V_F(\mu_1), \dots, V_F(\mu_{N_s})$}, {$p_F(\mu_1), \dots, p_F(\mu_{N_s})$} to obtain low dimensional spaces $Y_{N_f}, U_{N_v}$ and $P_{N_p}$ of dimension $N_f, N_v^3$ and $N_p$ respectively (Strazzullo 2017). For any finite-dimensional space $V$, denote $\Pi_V$ the orthogonal projector onto this subspace. The goal of the POD is to obtain a few basis vector giving the essential information in some data and the algorithm provides finite dimensional spaces that minimizes the empirical projection errors:

$$ \frac{1}{N_s} \sum_{\mu \in \mathcal{M}_s} \| f_F(\mu) - \Pi_{Y_{N_f}} f_F(\mu) \|_{\mathcal{L}}^2 \quad (11) $$

$$ \frac{1}{N_s} \sum_{\mu \in \mathcal{M}_s} \| V_F(\mu) - \Pi_{U_{N_v}} V_F(\mu) \|_{\mathcal{L}}^2 \quad (12) $$

$$ \frac{1}{N_s} \sum_{\mu \in \mathcal{M}_s} \| p_F(\mu) - \Pi_{P_{N_p}} p_F(\mu) \|_{\mathcal{L}}^2 \quad (13) $$

for the density, velocity field and adjoint variable respectively.

We will detail the construction of the POD basis for the density. Similar methods have been used for the velocity components and the adjoint variable. The POD basis {$\psi_{f,j}$}$_{1 \le j \le N_f}$ of a rank $N_f$ is obtained with the snapshots {$y_j$}$_{1 \le j \le N_s}$ = {$f(\mu_j)$}$_{1 \le j \le N_s} \in \mathbb{R}^m$. The POD problem can be formulated as a constrained optimization problem and solved with a Lagrangian framework. The first-order necessary optimality conditions are related to the SVD decomposition of the matrix $Y = [y_1, \dots, y_{N_s}] \in \mathbb{R}^{m \times N_s}$ whose columns are given by the snapshots $y_j$ (Volkwein 2013; Ghattas 2021). The singular values decomposition of the matrix $Y$ can be written:

$$ \Psi_f^t Y \Phi_f = \begin{pmatrix} D & 0 \\ 0 & 0 \end{pmatrix} = \Sigma \in \mathbb{R}^{m \times N_s} \quad (14) $$

where $D = \text{diag}(\sigma_1, \dots, \sigma_d) \in \mathbb{R}^{d \times d}$ with eigenvalues in decreasing order $\sigma_j \ge \sigma_{j+1} \ge 0$, and $\Sigma$ is a matrix of rank $d$. The matrices $\Phi_f \in \mathbb{R}^{N_s \times N_s}$ and $\Psi_f \in \mathbb{R}^{m \times m}$ are orthogonal matrices with columns {$\phi_{f,j}$}$_{1 \le j \le N_s}$ and {$\psi_{f,j}$}$_{1 \le j \le m}$ respectively. The vectors $\psi_{f,j}$ and $\phi_{f,j}$ are eigenvectors of the empirical non centered covariance matrices $YY^t$ and $Y^tY$ respectively with eigenvalues $\lambda_j = \sigma_j^2$ for $j = 1, \dots, d$ and they satisfy:

$$ Y \phi_{f,j} = \sigma_j \psi_{f,j} \quad Y^t \psi_{f,j} = \sigma_j \phi_{f,j} \quad \text{for } j = 1, \dots, d \quad (15) $$

In the following, the POD basis of rank $N_f$ is determined as follows: We solve the $N_s \times N_s$ eigenvalue problem

$$ Y^t Y \phi_{f,i} = \lambda_i \phi_{f,i} \quad \text{for } i = 1, \dots, N_f \quad (16) $$
and compute the $N_f$ eigenvectors $\phi_{f,1} ... \phi_{f,N_f} \in \mathbb{R}^{N_s}$ and then set:

$$ \psi_{f,i} = \frac{1}{\sqrt{\lambda_i}} Y \phi_{f,i} \quad \text{for} \quad i = 1, \dots, N_f \qquad (17) $$

Several reduced basis matrix of rank $N_f$, $N_v$ and $N_p$ are obtained for the density, the velocity components and the adjoint variables with this methodology. They will be denoted $\Psi_f, \Psi_{v_x}, \Psi_{v_y}, \Psi_{v_z}$ and $\Psi_p$ respectively in the following.

### 3.2. Online step: POD reduced model for the inverse problem

During the online phase, we perform a Galerkin projection of $f(\mu), V(\mu), p(\mu)$ onto $Y_{N_f}, U_{N_v}$ and $P_{N_p}$ respectively. The reduced model with the reduced operators can be computed once the basis for the density, the velocity and the adjoint variable are defined. Model reduction is used to find a low-dimensional approximation space and then a system of ordinary differential equations is obtained in the low dimensional space. To form the reduced model, the density $f$ is approximated in the POD basis

$$ f = \Psi_f \hat{f} \qquad (18) $$

where $\Psi_f \in \mathbb{R}^{m \times N_f}$ is the basis matrix for the density. The reduced state $\hat{f}$ represents the coefficients of expansion in the basis. Given a linear operator $A$ and using Galerkin projection, we obtain the reduced model

$$ \frac{d\hat{f}}{dt} = \hat{A}\hat{f} \qquad (19) $$

with $\hat{f}(0) = \Psi_f^t f(0)$. The reduced-order operator is defined by projection of the full-order operator $\hat{A} = \Psi_f^t A \Psi_f$. This equation is solved with finite difference discretization in time. After finite difference discretization in time, the PDE constraint $e(f, V, \mu) = 0$ can be rewritten as

$$ A^n f = f^n \qquad (20) $$

with

$$ A^n f = f + (u^n \frac{\partial f}{\partial x} + v^n \frac{\partial f}{\partial y} + w^n \frac{\partial f}{\partial z} - D \Delta f) \Delta t \qquad (21) $$

where $\Delta t$ denotes the time step of the discretization and $f^n, u^n, v^n$, and $w^n$ refer to the density and velocity components at time step $n$ projected on the low dimensional subspaces $Y_{N_f}$ and $U_{N_v}$. After projection onto the density space, the $N_f$ coefficients of the expansion of the density at time step $n+1$ are calculated as the solution of a linear system with a matrix in $\mathbb{R}^{N_f \times N_f}$.

Similar model reductions can be applied to the adjoint equation and to the gradient of the functional. For the adjoint equation, the reduced ordinary differential equation is obtained with the matrix $\Psi_p$ build with the adjoint snapshots. The velocity gradient is projected on the velocity basis. For instance, let us denote $\nabla \hat{J}_x$ the component along $x$
of the gradient obtained with Eq.9, then the component along $x$ of projected gradient is calculated with:

$$ \nabla \tilde{J}_x = \Psi_{V_x} \Psi_{V_x}^t \nabla \hat{J}_x \quad (22) $$

The component along $y$ and $z$ of the reduced gradient are derived similarly with the basis matrices $\Psi_{V_y}$ and $\Psi_{V_z}$.

# 4. Numerical Experiments

In this section, we detail some numerical experiments to test the former approach on a realistic numerical phantom with complex flows described by Navier-Stokes equations.

## 4.1. Finite elements simulation details

The numerical experiments have been performed with a realistic phantom displayed in Figure 2. It consists of a 10 cm long vessel with a bifurcation. The average transverse diameter is around 10 millimeters. The finite elements simulations have been implemented with the Fenics software (Goda 1979; Logg 2012). For the discretization of the time domain $[0, T] = [0s, 0.5s]$, we have used a regular grid of size $N_t = 50$ with a small time step $dt = 0.01s$ in agreement with the Courant-Friedrich-Lax (CFL) condition. The velocity fields to be reconstructed $V^*(x)$ or $V^*(x, t)$ are calculated with the Navier-Stokes equations. These velocity fields are used with the transport equation Eq.2 to obtain the ground truth density of the contrast agent, $f^*$. The approach presented in Section 3 combining the adjoint method and POD decompositions is then applied to reconstruct the velocity field, starting from a first guess solution. The PDE transport equation Eq.2 is solved with finite elements for the space dimension and with the Euler explicit discretization for the time dimension (Morton 2014). The simulation domain is discretized with $P1$ finite elements and 56845 nodes. The term inverse crime is used when the finite difference or finite element grid used in the inversion process is the same as the one used to compute the synthetic data, reducing the problem to a well-posed finite-dimensional one that behaves fundamentally different from the original infinite-dimensional one (Kaipio 2005). In order to avoid this, the data for the forward model have been generated with a discretization with 102500 nodes. Stability conditions $Ddt/dx^2 < 1$ and CFL conditions are fulfilled with the chosen space and time discretizations (Morton 2014). The values for the diffusion coefficient and the kinematic viscosity are $D = 0.001cm^2/s$ and $\nu = 0.01cm^2/s$ respectively. The proposed reconstruction method with PODs has been tested on velocity fields obtained with the stationary and non-stationary Navier-Stokes equations detailed in (Huang 2021). It is well-known that stable spatial discretizations are required for these type of simulations. In this work, we have used the Taylor-Hood ($P2-P1$) finite elements (Hood 1973). The algorithm used for the nonstationary Navier-Stokes velocity field is the Incremental Pressure Correction scheme in Fenics(Goda 1979; Logg 2012). The Radon projections perpendicular to the main vessel direction are estimated with the Scipy python library (Virtanen et al. 2020) for $N_z = 40$ values regularly spaced. The Radon projections are measured for $N_p = 100$ values, for $N_\theta = 180$ projection angles in the angular range $\Delta\theta(t) = [p\delta\theta, (p+1)\delta\theta]$, with a measurement interval $\delta\theta=180^\circ$ and for each time step $t_p = pT/N_t$ in the time domain, $0 \le t_p \le T$. A Gaussian white
noise with peak-to-peak signal-to-noise ratio (PPSNR) between 0 and 20 dB is used to corrupt the projections. The peak-to-peak signal to noise ratio is defined as:

$$ PPSNR = 20\log\left(\frac{S_{max}}{n_{max}}\right) \qquad (23) $$

where $S_{max}$ is the maximum signal amplitude and $n_{max}$ the maximum noise amplitude. The simulation results presented in the following subsections are obtained with the initial velocity field $V_0$ such that $V_{0,x} = 0.75V_x^*$, $V_{0,y} = 0.75V_y^*$, $V_{0,z} = 0.75V_z^*$ where $V^*$ is the ground truth velocity. The optimal decrease of the discrepancy term or of the Lagrangian is used as a criteria for the choice of the regularization parameter $\alpha_s$ and of the Lagrangian parameter $\beta$.

In order to have quantitative error metrics for the reconstruction and to study the convergence characteristics of the POD approach, we have also estimated the evolution of the relative root mean square errors (RMSE) as a function of the iteration number $k$ for several quantities of interest using $L_2$ norms. These quantities include the velocity field $V_k$, its components and the reconstructed density $f_k$ for the iteration $k$. After discretization in time, the norm given in Eq.1 is equivalent to $N(f) = (\sum_{p=1}^{N_t} \|f(\cdot, t_p)\|_{L_2(\Omega)}^2)^{1/2}$. The error metrics are defined as follows:

$$ RMSE(f_k) = N(f_k - f^*)/N(f^*) \qquad (24) $$

$$ RMSE(V_k) = \frac{(N(u-u_k)^2+N(v-v_k)^2+N(w-w_k)^2)^{1/2}}{(N(u)^2+N(v)^2+N(w)^2)^{1/2}} \qquad (25) $$

$$ RMSE(w_k) = N(w - w_k)/N(w) \qquad (26) $$

Similar definitions will be used for the components of the velocity along the x and y axis.

### 4.1.1. POD bases construction

In this section, we detail the POD bases construction. The inverse problem depends on a vector $\mu$ of parameters. The snapshots are sampled both in time and parameter space. The POD has been applied for time varying and parametrically varying conditions and we have used regular grid sampling over the desired range of conditions. The parameter $\mu_1$ corresponds to the initial/boundary conditions for the velocity profile in the whole volume. To obtain realistic velocity snapshots, stationary and non-stationary Navier-Stokes equations have been considered (Glowinski 1992; Temam 2001). The unsteady NS equation in an incompressible viscous flow is:

$$ \left\{ \begin{array}{l} \frac{\partial \mathbf{V}}{\partial t} + \mathbf{V} \cdot \nabla \mathbf{V} + \nabla p - \nu \Delta \mathbf{V} = 0 \\ \nabla \cdot \mathbf{V} = 0 \\ \mathbf{V} = \mathbf{V}_{\partial\Omega} \quad \text{on } \partial\Omega \end{array} \right. \qquad (27) $$
<figure><img src="image_4.png" /><figcaption>Figure 2. The finite element mesh used for the simulations.</figcaption></figure>

and the steady-state equation is given by :

$$ \left\{ \begin{array}{l} \mathbf{V} \cdot \nabla \mathbf{V} + \nabla p - \nu \Delta \mathbf{V} = 0 \\ \nabla \cdot \mathbf{V} = 0 \\ \mathbf{V} = \mathbf{V}_{\partial \Omega} \quad \text{on } \partial \Omega \end{array} \right. \qquad (28) $$

where $\nu$ is the kinematic viscosity, $p(\mathbf{x})$ the pressure field and $\mathbf{V}_{\partial\Omega}$ the velocity field on the boundary $\partial\Omega$. The boundary $\partial\Omega = \Gamma_{in} \cup \Gamma_{out} \cup \Gamma_s$ of the spatial domain $\Omega$ consists of an inflow part $\Gamma_{in}$ on the plane $z = 0$, of an outflow boundary $\Gamma_{out}$ on the output plane and of the lateral side of the vessel $\Gamma_s$. We will assume that the velocity on the boundary $\Gamma_s$ vanishes and we will use the velocity $\mathbf{V}_{\Gamma_{in}}$ on the inflow part $\Gamma_{in}$ and the pressure $p_{\Gamma_{out}}$ on the outflow part $\Gamma_{out}$ as control variables. We will assume that the velocity on the boundary $\Gamma_s$ vanishes and we will use the velocity $\mathbf{V}_{\Gamma_{in}}$ on the inflow part $\Gamma_{in}$ and the pressure $p_{\Gamma_{out}}$ on the outflow part $\Gamma_{out}$ as control variables. The Dirichlet boundary condition at the inlet is thus $\mathbf{V} = \mathbf{V}_{\Gamma_{in}}$ on $\Gamma_{in}$. For the outlet, we use $\nu\partial_{\mathbf{N}}\mathbf{V} + p\mathbf{N} = 0$ where $\mathbf{N}$ the unit normal vector and $\partial_{\mathbf{N}} = \mathbf{N} \cdot \nabla$ the normal derivative. A parabolic profile has been chosen with a maximum value $\mathbf{V}_{\Gamma_{in,max}}$ in the middle of the inlet boundary. The stationary velocity fields will thus be parametrized with $\mu_1 = (\mathbf{V}_{\Gamma_{in,max}}, p_{\Gamma_{out}})$. The velocity snapshots have been obtained with regular grid sampling of the inlet velocity and of the outlet pressure. For the inlet velocity, 10 values have been chosen between 0.2 and 0.6 m/s and for the outlet pressure 10 values have been selected between 0 and 1000Pa. For the non-stationary velocity field, the inlet velocity is defined as:

$$ \left\{ \begin{array}{l} \mathbf{V}_{\Gamma_{in,max}} + \mathbf{V}_{\Gamma_{in,max}} \sin(\pi t) \quad \text{for } t \le 0.5 \\ \mathbf{V}_{\Gamma_{in,max}} + \mathbf{V}_{\Gamma_{in,max}} (3/2 - 0.5 \cos(2\pi(t - 0.5))) \quad \text{for } t \le 0.5 \end{array} \right. \qquad (29) $$
The unsteady velocity fields will thus be parametrized with $\mu_1 = (V_{\Gamma_{in,max}}, p_{\Gamma_{out}}, t_i)$, where $(t_i)_{1\le i\le 50}$ are the time steps of the simulation time interval. To compute the POD basis for density we consider snapshots which are state solutions computed with the convection-diffusion Eq.2 at different times or for different velocity field and inlet density. The parameter $\mu_2$ corresponds to the inlet boundary condition for the density. Three different inlet densities between 0.5 and 1.5 kg/m³ have been used. A similar method has been used for the adjoint snapshots. The reduced states are for the density, the velocity and the adjoint variable are obtained from these snapshots as explained in Section 3.1.

The POD bases solve the minimization of the projection errors given by Eq.11, 12 and 13 over the set of $N_f$, $N_v$ or $N_p$ dimensional subspaces. It is well-known that the empirical projection error is the sum of squares of the singular values for the singular values not included in the pod basis. For the density for instance, the projection error is given by $\sum_{i=N_f+1}^{N_s} \sigma_i^2$. The quality of the POD approximation is thus controlled by rate of decay of the spectrum of the empirical covariance matrix. For the density, we chose the basis number with a tolerance $\kappa$ for the singular values such that:

$$ \frac{\sum_{i=1}^{N_f} \sigma_i^2}{\sum_{i=1}^{N_s} \sigma_i^2} > \kappa \qquad (30) $$

with $\kappa = 1\%$. The same choice criteria have been used for the velocity and adjoint variables. In order to check that the bases considered are sufficient and they well represent the full order approximation solutions, we have calculated the average relative error for the density, velocity and adjoint variable for several values of $\mu$. For a given value of $\mu$, and a time step $t_k$, the relative error for the density $f_{\mathcal{F}}$ can be written as:

$$ e_f(t_k) = \| \sum_{i=1}^{N_f} \langle f_{\mathcal{F}}(t_k), \psi_i^f \rangle_{\mathbb{R}^m} \psi_i^f - f_{\mathcal{F}}(t_k) \| / \| f_{\mathcal{F}}(t_k) \| \qquad (31) $$

where $N_f$ is the number of basis function used for the density. Monte-carlo estimates of the relative errors have obtained with 100 testing examples collected with a uniform sampling of the $\mu$ parameter.

For the stationary velocity field, the basis functions for the velocity, the density and the adjoint variables are obtained by retaining the first 2, 75 or 6 eigenvectors of the snapshots matrix respectively. The total reduced dimension is thus $2^3 * 75 * 6$. For the non-stationary velocity field, we have selected $N_f = 6$, $N_v = 40$ and $N_p = 2$. With these values, the relative projection errors are below 1% and with these truncations the PODs reach a high accuracy. We have thus obtained a very low dimensional spaces which remains very accurate. In the next section, we show that the use of these bases is an effective method for the solution of the inverse problem.

## 4.2. Results and discussion

In this section, we present numerical results and we compare the proposed POD method against finite elements simulations with full order approximation. We will test the proposed algorithm on two examples, with a stationary velocity field and a non-stationary one. We will compare the results for different noise levels, visualize the recovered solutions and analyze the convergence properties of the various approaches on two benchmark problems.
### 4.2.1. Stationary velocity field

In order to illustrate the results with a stationary velocity field, we display the reconstruction results for the velocity and density fields obtained with $\mu_1 = (0.55\,m/s, 500\,Pa)$ and an inlet density of $\mu_2 = 1.1\,kg/m^3$. Cross-sections of the reconstructed density and of the true density are displayed on Figure 3. Some examples of reconstructed velocity fields $w$ along the $z$ axis obtained with the POD based approach are displayed in Figure 4 together with the ground truth velocity for a low noise level, for selected simulation times and cross-sections. Similar figures are displayed for the component along the $y$ and $x$ components in Figure 5 and Figure 6 respectively. These figures show that good reconstruction results are obtained with the proposed approach. Figure 7 and 8 compare the error maps for the velocity and the density for the POD and full order approximation. Figures 9, 10, 11, 12 and 13 display the evolution with the iterations of the RMSE for the density, for the velocity components along different directions and for the velocity norm, for different noise levels. These figures compare the errors as a function of the iterations for POD and full order approximation simulations. The reconstruction errors are summarized in Table 1, 2 and 3 for the density, the velocity components and the velocity norm for the two noise levels considered and the two reconstruction strategies. It should be noted that similar reconstruction errors are obtained for the range of parameters $\mu$ investigated.

With the POD approach, the reconstruction errors for the density and the velocity show large and fast reductions and the method is able to give a good approximation of these unknowns fields. The reconstruction accuracy for the density and the $z$ component of the velocity decreases when the noise level on the projections increases. We obtain also good reconstruction results for the transverse components of the velocity. These components are less sensitive to the noise corrupting the Radon projection. With the proposed approach and the partial differential equation constraint, the simultaneous reconstruction of the axial and transverse components of the velocity can be achieved. When more basis functions are included in the POD approach, the errors stagnates and the basis can be considered as optimal for the inversion for the range of $\mu$ values studied.

The proposed approach significantly outperforms the full order approximation method. The errors for the density and velocity field are much decreased with respect to the ones obtained with the classical finite element simulations. It can be seen on figures 8 and 9 that the POD based inversion strategy suppresses reconstruction errors on the boundaries of the vessel. As illustrated in these figures, the full order simulation has an non-optimal behavior and the optimization problem is stuck in local minima away from the optimal solution. Increasing the number of parameters for the unknown field does not improve the error due to the increased complexity of the optimization problem. It is thus beneficial to use reduced bases to reduce the complexity of the optimization and to obtain a better approximation of the solution.

The POD inversion method has a high offline cost. The offline cost for building a reduced model consists in obtaining the full model snapshot solutions, plus the computational cost of forming the POD bases. To enable large scale simulations, the offline basis construction for the model reduction could be implemented for parallel computation on distributed-memory clusters. Without this parallelel implementation, the repeated solution of the full-scale system of equations requires 20h on a modern HPC cluster with ten processors and each of the 15000 snapshots is around 5M.

Yet, it has a much lower online cost than the full order simulation. The number of degrees of freedom and the average simulation times required for the Eq.2 and
one iteration of the whole adjoint method are summarized in Table 4. The speed up index is defined as the ratio of the computation time of the FOA solution divided by the computation time of the reduced basis approximation. It can be estimated to 90 for the forward step and to 15 for the inversion method with the large decrease of the number of degrees of freedom. With the basis numbers considered and our implementation, the runtime of the method is much reduced with the POD approach. This makes it the most suitable approach for applications where the velocity field has to be estimated several times and for real-time applications. The results of this section thus clearly demonstrate the effectiveness of the proposed methodology to solve the inverse problem for a stationary velocity field.

### 4.2.2. Non-stationary velocity field

In this section, we summarize the results obtained with an non-stationary velocity field. It must be mentioned that we are not considering the regimes where the Navier-Stokes equations become chaotic and turbulent. We display the reconstruction results for the velocity and density corresponding to $\mu_1 = (0.55\,\text{m/s}, 500\,\text{Pa})$ and an inlet density of $\mu_2 = 1.1\,\text{kg/m}^3$. The velocity field is calculated with Eq. 29. The evolution with the iterations of the relevant quantitative metrics, RMSE for the density, for the velocity norm and for the velocity components are displayed in figure 14 and figure 15 for different noise levels and the two reconstruction approaches. The reconstruction errors are summarized in Table 4 and 5. The reconstruction errors for the density and velocity fields are much reduced with the proposed approach and it clearly outperforms the full order approximation method. The largest decrease is obtained for the $z$ component of the velocity. The accuracy of the reconstruction of the $z$ component of the velocity and of the density is weakly dependent on the noise level. These results are similar to the ones obtained for the stationary velocity and a simultaneous reconstruction of the different components of the velocity field can be obtained for the noise levels studied. The POD methodology presented can be used for more challenging inverse problems such as those arising with time-dependent velocity fields.

### 4.2.3. Discussion

The main contribution of this study was to demonstrate, through several realistic examples, the ability of the algorithm to infer the velocity field from the Radon projections by leveraging the known underlying transport equation. The approach is able to make accurate predictions for steady or time dependent flows. We have investigated relatively high noise levels but the dominant flow features can be extracted through dimensionality reduction. We have shown that the POD methodology shows good numerical performance for stationary and non-stationary velocity fields. The surrogate model obtained with POD decompositions combined with Galerkin projection is computationally efficient and it achieves lower reconstruction errors compared to approaches based on full order meshes. It is also a suitable approach to be used in order to reduce the computational cost. With the conventional full order approximation, one needs to solve significantly more expensive and high-dimensional optimization problems. The information on the density and its gradients in the simulation domain and near its boundaries is not sufficient to infer a good solution for the velocity field. The POD algorithm is capable of accurately reconstruct the velocity fields without having access to full observations of the density field, and without the knowledge of initial and boundary conditions. Moreover, the boundary conditions for the velocity are more
accurately recovered.

The presented approach can be used for any application dealing with incompressible transient or turbulent fluid dynamic problems that has the response depending on parameter changes and demanding real-time computing and visualization. It can be useful in fluids mechanics and hydrology where the underlying physical processes are heterogeneous in both space and time. Our primary goal was to characterize blood flow in vessels for clinical applications. Assessing the hemodynamics of the blood inside the cardiovascular system is very useful for patients with cardiovascular diseases. Yet, applications can be found in many engineering physical problems involving industrial viscous flows.

We can mention several limitations of the proposed approach. The transformation used to obtain the POD basis is an orthogonal linear transformation. Linear projection-based reduced order methods may not be efficient for dimensionality reduction for complex flow patterns because the number of required modes increases significantly. We have assumed the blood to be a Newtonian fluid. Whereas blood can be assumed to be a Newtonian fluid in large vessels, its viscosity varies with flow rate in smaller arteries and capillaries. For medical conditions such as stenosis, a non-Newtonian model may be preferred over a Newtonian model. More complex physics should be included in models to account for realistic conditions, such as non-Newtonian blood behavior and viscoelastic properties of the arterial wall. We have not considered a variation of the domain/structure and fluid-structure interaction problems. More generally, we have assumed that the PDE and operators governing the fluid dynamics and tracer propagation are known during the projection step. Another possible limitation of the current study is the use of synthetically generated data.

A simple generalization of the method can be obtained with replacing the Navier-Stokes equation with a nonlinear non-Newtonian equation. The presented approach based on POD can be used for higher Reynolds numbers but with more POD basis elements. We could also generalize the approach with an efficient geometrical parametrization for moving domains since our method exploits a parametrized formulation of the problem. The method could thus be extended to varying domain with fluid-structure interactions. This work leaves many interesting directions open for future research. Many concepts from deep learning can be used with classical methods to solve complex inverse problems involving dynamical systems. It is possible to combine model reduction and neural networks and deep learning approaches for parametric PDEs (Bhattacharya 2021). The reduced bases used in this work seems sufficiently rich and the projected model can give accurate solution of the inverse problem. Yet, we have mentioned that POD may not the optimal choice for dimension reduction. The development of deep learning methods and autoencoders on function spaces may outperform the principal component analysis type method proposed here.

Neural network-based autoencoders have been investigated as an alternative for nonlinear approximation because they can address some of the limitations of linear projection techniques (Gonzalez 2018). They provide a greater flexibility for the dimensionality reduction.

More generally, our geometry is precisely defined but the initial and boundary conditions for the density and the velocity field are not known accurately. Yet, we have assumed that the available data, i.e the density field, respect a given physical law, the transport equation, described by a partial differential equation. The prior knowledge of the underlying physical law introduces a prior that effectively regularizes the minimization procedure and help us to solve the inverse problem and to infer the velocity field. In future works, we will investigate non intrusive and purely data driven
<table><caption>**Table 1.** Comparison of density reconstruction errors for POD and FOA for two noise levels. Initial test velocity 0.55m/s, outlet pressure 500Pa.</caption><thead><tr><td>Noise level(dB)</td><td>POD</td><td>FOA</td></tr></thead><tbody><tr><td>0</td><td>0.024</td><td>0.053</td></tr><tr><td>20</td><td>0.014</td><td>0.046</td></tr></tbody></table>

<table><caption>**Table 2.** Comparison of velocity reconstruction errors for POD and FOA for two noise levels. Initial test velocity 0.55m/s, outlet pressure 500Pa.</caption><thead><tr><td>Noise level(dB)</td><td>POD</td><td>FOA</td></tr></thead><tbody><tr><td>0</td><td>0.190</td><td>0.475</td></tr><tr><td>20</td><td>0.098</td><td>0.442</td></tr></tbody></table>

methods.

## 5. Conclusion

In this work, we have considered a new POD approach for the reconstruction of the velocity field of the blood inside a vessel. The method is based on the propagation of a contrast agent and on the acquisition of Radon projections perpendicularly to the main propagation direction. The PDE constraint modeling the transport of the contrast agent is used to eliminate the density in the optimization problem. We leverage our former solution of the inverse problem with the adjoint method. The full order approximation is replaced by reduced order approximations with POD for the density, the velocity and the adjoint fields. We demonstrate the effectiveness of this new methodology and its numerical performance on a vessel with a bifurcation with stationary and non-stationary velocity flows based on the Navier-Stokes equations. Good reconstructions results are obtained for the velocity and the density. The proposed algorithms produced better results than the full order approximation. The computational cost is much reduced and the accuracy is improved.

<table><caption>**Table 3.** Comparison of velocity components reconstruction errors for POD and FOA for two noise levels. Initial test velocity 0.55m/s, outlet pressure 500Pa.</caption><thead><tr><td>Direction</td><td>Noise level(dB)</td><td>POD</td><td>FOA</td></tr></thead><tbody><tr><td>X</td><td>0</td><td>0.132</td><td>0.176</td></tr><tr><td>X</td><td>20</td><td>0.054</td><td>0.166</td></tr><tr><td>Y</td><td>0</td><td>0.013</td><td>0.173</td></tr><tr><td>Y</td><td>20</td><td>0.044</td><td>0.160</td></tr><tr><td>Z</td><td>0</td><td>0.0442</td><td>0.115</td></tr><tr><td>Z</td><td>20</td><td>0.032</td><td>0.116</td></tr></tbody></table>
<table><caption>**Table 4.** Comparison of the number of degrees of freedom, average simulation time for the transport equation (Eq.2), and for one iteration of the online step of the whole adjoint method implemented with POD and FOA and a stationary velocity field.</caption><thead><tr><th></th><th>POD</th><th>FOA</th></tr></thead><tbody><tr><td>Mesh size</td><td>-</td><td>56845</td></tr><tr><td>Degrees of freedom</td><td>3600</td><td>56845<sup>5</sup></td></tr><tr><td>Simulation time for Eq.2 (s)</td><td>10<sup>-2</sup> (online)</td><td>0.9</td></tr><tr><td>Simulation time for adjoint method (s)</td><td>35</td><td>450</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Comparison of horizontal sections of the ground truth and reconstructed density for z=0.125 and for a PPSNR=20 dB.</figcaption></figure>

## 6. Acknowledgment

The authors acknowledge financial support of the China Scholarship Council(CSC).

## References

* Adams RA. 1975. Sobolev Spaces. New York: Academic Press.

* Ballarin F, Faggiano A, Manzoni A, Quarteroni G, Rozza G, Ipollito S, Antona C, Scrofani R. 2017. Numerical modeling of hemodynamics scenarios of patient-specific coronary artery bypass grafts. Biomechanics and modeling in mechanobiology. 16:1373–1390.

* Barfett JJ, Velauthapillai N, Fierstra J. 2014. Intra-vascular blood velocity and volumetric flow rate calculated from dynamic 4D CT angiography using a time of flight technique. The international journal of cardiovascular imaging. 30:1383–1392.

* Bhattacharya K, Hosseini B, Kovachki NB, Stuart AM. 2021. Model reduction and neural networks for parametric PDEs. The SMAI journal of computational mathematics. 7:121–

<table><caption>**Table 5.** Comparison of velocity reconstruction errors for POD and FOA for two noise levels. Case 3: initial test velocity 0.55m/s, outlet pressure 500Pa.</caption><thead><tr><th>Noise level(dB)</th><th>POD</th><th>FOA</th></tr></thead><tbody><tr><td>0</td><td>0.374</td><td>0.404</td></tr><tr><td>20</td><td>0.378</td><td>0.409</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Comparison of horizontal sections of the ground truth and reconstructed velocity along x for z=0.125 and for a PPSNR=20 dB.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparison of horizontal sections of the ground truth and reconstructed velocity along y for z=0.125 and for a PPSNR=20 dB.</figcaption></figure>

157.

Bouillot P, Brina 0, Chnafa C, Cancelliere NM, Vargas MI, Radovanic I, Krings T, Steinman DA, Pereira VM. 2019. Robust cerebrovascular blood velocity and flow rate estimation form 4D-CTA. Medical physics. 46:2126–2136.

Carere G, Strazzullo M, Ballarin F, Rozza G, Stevenson R. 2021. A weighted POD-reduction approach for parametrized PDE-constrained optimal control problems with random inputs and application to environmental science. SIAM Journal Scientific Computing. 37:1111–1138.

Courant R. 1946. Calcul of Variations. New York.

Carere G, Strazzullo M, Ballarin F, Rozza G, Stevenson R. 2021. A weighted POD-reduction approach for parametrized PDE-constrained optimal control problems with random inputs and application to environmental science. SIAM Journal Scientific Computing. 37:1111–1138.

Crosetto P, Reymond P, Deparis S, Kontaxakis D, Stergiopulos N, Quarteroni A. 2011. Fluid-structure interaction simulation of aortic blood flow. Computers and fluids. 43:46–57.

Daly SM, Leahy MJ. 2013. Go with the flow: a review of methods and advancements in blood flow imaging. Journal of biophotonics. 6:217–255.

Du J, Fang F, Pain CC, Navo IM, Zhu J, Ham DA. 2013. Go with the flow: a review of methods and advancements in blood flow imaging. Computers and Mathematics with Applications.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Comparison of horizontal sections of the ground truth and reconstructed velocity along z for z=0.125 and for a PPSNR=20 dB.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Error map of densities reconstructed with POD and full order approximation for z=0.125 and for a PPSNR=20 dB.</figcaption></figure>

65:362–379.

Galdi GP. 1994. An introduction to the mathematical theory of the Navier-Stokes equations. Springer Tracts in Natural Philosophy. 38.

Ghattas O and Willcox K. 2021. Learning physics-based models from data: perspectives from inverse problems and model reduction. Acta Numerica. 445–554.

Girault V, Raviart PA. 1986. Finite element methods for Navier-Stokes equations. Springer Series in Computational Mathematics, Springer-Verlag, Berlin.

Glowinski P and Pironneau O 1992. Finite element methods for Navier-Stokes equations. Annual review of fluid mechanics. 24:167–204.

Goda K. 1979. A multistep technique with implicit difference schemes for calculating two or three dimensional cavity flows. Journal of Computational Physics. 30:76–95.

Gonzalez F and Balajewicz M. 2018. Deep convolutional recurrent autoencoders for learning low-dimensional feature dynamics of fluid systems arXiv preprint arXiv:1808.01346

Herzog M and Kunisch K. 2010. Algorithms for PDE-constrained optimization. GAMM-Mitteilungen. 33:163–176.

Hinze M, Pinnau R, Ulbrich M, Ulbrich S. 2009. Optimization with PDE constraints. Mathematical Modelling:theory and applications, Springer Netherlands.

Hood P, Taylor C. 1973. A numerical solution of the Navier-Stokes equations using the finite element technique. Comp.and Fluids. 1:73–100.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Error map of reconstructed velocity along z with POD or full order approximation for z=0.125 and for a PPSNR=20 dB.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Evolution of the density error as a function of the iteration number for different noise levels for the stationary velocity: plain line, full order approximation; dashed line, POD method; (△, 20 dB), (○, 0 dB).</figcaption></figure>

Huang S, Sigovan M, Sixou B. 2021. Reconstruction of vascular blood flow in a vessel from tomographic projections. Biomedical Physics and Engineering Express. 7:065032.

Kaipio J and Somersalo E 2005. Statistical and computational inverse problems. New York: Springer.

Korporaal JG, Benz MR, Schindera ST, Flohr TG, Schmidt B. 2016. Contrast gradient-based blood velocimetry with computed tomography. Investigative Radiology. 51:1–9.

Logg A, Margal KA, Wells GN. 2012. Automated solution of differential equations by the finite element method, The Fenics book. Lecture notes in computational science and engineering, Springer-Verlag Berlin Heidelberg.

Lu K, Jin Y, Chen Y, Yang Y, Hou L, Zhang Z, Li Z, Fu C. 2019. Review for order reduction based on proper orthogonal decomposition and outlooks of applications in mechanical systems. Mechanical Systems and Signal Processing. 123: 264–297.

MacDonald ME, Frayne R. 2015. Cerebrovascular MRI: a review of state-of-the-art approaches, methods and techniques. NMR in Biomedicine. 28:767–791.

Markl M, Schnell S, Wu C, Bollache E, Jarvis K, Barker AJ, Robinson JD, Rigsby CK.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Evolution of the velocity error as a function of the iteration number for different noise levels for the stationary velocity: plain line, full order approximation; dashed line, POD method; (△, 20 dB), (○, 0 dB).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Evolution of the error for the velocity along x as a function of the iteration number for different noise levels for the stationary velocity: plain line, full order approximation; dashed line, POD method; (△, 20 dB), (○, 0 dB).</figcaption></figure>

2016. Cerebrovascular MRI: a review of state-of-the-art approaches, methods and techniques. Clinical Radiology. 71:779-795.

Morton KW, Mayers DF. 2014. Numerical Solution of Partial Differential Equations. Cambridge University Press, Cambridge.

Natterer F. 1986. The Mathematics of Computerized Tomography. Wiley, Chichester, UK.

Pelc NJ, Sommer FG, Li KCP, Brosnan TJ, Herfkens RJ, Enzmann DR. 1994. Quantitative magnetic resonance flow imaging. J.Magn.Reson.Q. 10:125-147.

Plessix RE. 2006. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications. Geophys.J.Int. 167:495-503.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Evolution of the error for the velocity along *y* as a function of the iteration number for different noise levels for the stationary velocity: plain line, full order approximation; dashed line, POD method; (△, 20 dB), (o, 0 dB).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Evolution of the error for the velocity along *z* as a function of the iteration number for different noise levels for the stationary velocity: plain line, full order approximation; dashed line, POD method; (△, 20 dB), (o, 0 dB).</figcaption></figure>

Prevrhal S, Forsythe CH, Harnish RJ, Saeed M, Yeh BM. 2011. CT angiographic measurement of vascular blood flow velocity by using projection data. Radiology.261:923-929.

Quarteroni A, Manzoni A, Negri F. 2015. Reduced basis methods for partial differential equations: an introduction. Springer, New York.

Quarteroni A, Manzoni A, Vergara C. 2017. The cardiovascular system: mathematical modeling, numerical algorithms, clinical application. Acta Numerica. 001-202.

Ravindran SS, Sivaguru S. 2000. A reduced-order approach for optimal control of fluids using proper orthogonal decomposition. International journal for numerical methods in fluids.
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Evolution of the error for the density and the velocity as a function of the iteration number for different noise levels for the non-stationary velocity: plain line, full order approximation; dashed line, POD method; (△, 20 dB), (○, 0 dB).</figcaption></figure>

34:425–446.

Roubicek T. 2013. Nonlinear partial differential equations with applications. International Series of Numerical Mathematics, Birkhauser.

Rowley CW. 2005. Model reduction for fluids, using balanced proper orthogonal decomposition. International Journal of Bifurcation and Chaos. 15:997–1013.

Shi Y, Lawford P, Hose R. 2011. Review of zero-d and 1-d models of blood flow in the cardiovascular system. Biomedical engineering on line. 10:1–38.

Strazzullo M, Ballarin F, Mossetti R, Rozza G. 2017. Model reduction for parametrized optimal control problems in environmental marine sciences and engineering. SIAM Journal Scientific Computing. 40:1055–1079.

Temam R. 2001. Navier-Stokes equations. AMS Chelsea Publishing, Providence, RI.

Virtanen P. 2020. SciPy: Fundamental algorithms for scientific computing in python. Nature methods. 17:261–272.

Volkwein S. 2013. Proper orthogonal decomposition: theory and reduced-order modelling. Lecture notes, University of Konstanz. 4:1–29.
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Evolution of the error for the velocity along y and z as a function of the iteration number for different noise levels for the non-stationary velocity: plain line, full order approximation; dashed line, POD method; (△, 20 dB), (○, 0 dB).</figcaption></figure>