# Physics-informed compressed sensing for PC-MRI: an inverse Navier–Stokes
problem

Alexandros Kontogiannis, Matthew P. Juniper Engineering Department, University
of Cambridge, Trumpington Street, Cambridge CB2 1PZ, UK

###### Abstract

We formulate a physics-informed compressed sensing (PICS) method for the
reconstruction of velocity fields from noisy and sparse phase-contrast
magnetic resonance signals. The method solves an inverse Navier–Stokes
boundary value problem, which permits us to jointly reconstruct and segment
the velocity field, and at the same time infer hidden quantities such as the
hydrodynamic pressure and the wall shear stress. Using a Bayesian framework,
we regularize the problem by introducing a priori information about the
unknown parameters in the form of Gaussian random fields. This prior
information is updated using the Navier–Stokes problem, an energy-based
segmentation functional, and by requiring that the reconstruction is
consistent with the $k$-space signals. We create an algorithm that solves this
reconstruction problem, and test it for noisy and sparse $k$-space signals of
the flow through a converging nozzle. We find that the method is capable of
reconstructing and segmenting the velocity fields from sparsely-sampled (15%
$k$-space coverage), low ($\sim$$10$) signal-to-noise ratio (SNR) signals, and
that the reconstructed velocity field compares well with that derived from
fully-sampled (100% $k$-space coverage) high ($>$$40$) SNR signals of the same
flow.

###### Index Terms:

phase contrast magnetic resonance imaging (PC-MRI), physics-informed
compressed sensing, velocity reconstruction and segmentation

##  I Introduction

In phase contrast magnetic resonance imaging one seeks to reconstruct complex
images (${w}$) whose magnitude is proportional to the nuclear spin density
($\rho$) and whose phase difference ($\Delta\varphi$) is proportional to a
flow velocity component ($u_{i}$). For fully-sampled $k$-space signals
(${s}^{\bullet}$) there is a one-to-one correspondence between the $k$-space
and the physical (complex) space, which is given by a Fourier transform
($\mathcal{F}$) such that $s^{\bullet}=\mathcal{F}w$. Then, for sufficiently
high ($>3$) signal-to-noise ratios (SNR), the noise in the magnitude and the
phase of the complex image can be assumed to be white and Gaussian [1]. In
this case, it is reasonable to directly reconstruct (denoise) the complex
image in physical space, using either general purpose image denoising
algorithms [2, 3], or physics-informed algorithms [4, 5, 6, 7, 8, 9]. For
sparsely-sampled $k$-space signals (${s}^{\star}$) there is no longer a one-
to-one correspondence between $k$-space and physical space. A naive solution
to this problem is to zero-fill the subsampled $k$-space signal and then
perform an inverse Fourier transform. This approach, known as the ‘zero-
filling solution’, reconstructs a complex image that is corrupted by artefacts
and interference noise. Incoherent subsampling can eliminate artefacts, but
interference noise remains a problem. In this case, it is better to
reconstruct the complex image using compressed sensing (CS) [10, 11, 12].

Figure 1: A $d$-dimensional velocity field ($\boldsymbol{u}$) that can be
described by a Navier–Stokes problem ($\mathcal{Z}$) has an underlying
($d-$1)-dimensional structure
($\boldsymbol{x}=\mathcal{Z}^{-1}\boldsymbol{u}$), which is the parameter
vector containing the boundary of the object ($\partial\Omega$), the boundary
conditions $(\boldsymbol{g}_{i},\boldsymbol{g}_{o})$, and the kinematic
viscosity ($\nu$). Images of $n^{d}$ voxels depicting the $d$ velocity
components can be compressed/decompressed by solving an inverse/forward N–S
problem.

Compressed sensing relies on the assumption that the complex image, or its
constituent parts (e.g. the magnitude and the phase), have a sparse
representation in a transform domain that is considered to be known based on a
priori knowledge. The reconstructed complex image can then be recovered by
solving a (generally) nonlinear optimization problem of the form [13]

$\textrm{min}\ \big{\lVert}Tw\big{\rVert}_{\ell^{1}}\quad\text{subject
to}\quad\big{\lVert}s^{\star}-\mathcal{P}\mathcal{F}w\big{\rVert}_{\ell^{2}}<N_{s}\varepsilon\quad,$

where $T$ is a sparsity-inducing transformation, $\mathcal{P}$ is the sparse
sampling pattern, $N_{s}$ is the total number of sampled $k$-space points, and
$\varepsilon\propto\sigma$ is a user-selected value proportional to the
standard deviation of noise ($\sigma$). In the above example, the functional
to be minimized acts as a regularizer while the constraint enforces
consistency between the reconstruction and the measurement. It is also
possible to decompose the complex image into real and imaginary [14, 15], or
into magnitude and phase components [16][17, Chapter 9.3], and reconstruct
them separately. With this decomposition, different regularizers can be used
for each component. This is particularly important for PC-MRI because suitable
regularizers for the reconstructed magnitude, e.g. total variation (TV), may
produce large errors or artefacts (e.g. staircasing) in the reconstructed
phase. This is because the magnitude image often has high-contrast features
resembling a piecewise-constant function, which is sparse in the TV-transform
domain [12, 13], while the phase usually has smooth features since it encodes
velocity information, and thus is no longer sparse in the same domain. A
suitable transform for the phases, and for piecewise-smooth functions in
general, is second-order total generalized variation (TGV2) [18, 13].
Complications arise due to the fact that the phase of the complex image alone
does not represent velocity, but the phase difference does. Therefore, if
there is available a priori information on the velocity, the regularization
should be imposed on the phase difference rather than on the individual phases
[17, Chapter 9.3][19].

Even though generic CS methods perform very well in magnitude reconstruction,
accurate velocity reconstruction from subsampled $k$-space signals remains a
challenge. It has been suggested [17, Chapter 11.1] that even sparser PC-MRI
signals could be reconstructed if a regularization method based on the
Navier–Stokes (N–S) equations is used. However, the effect of the boundary
conditions should be equally important in order to capture the velocity
profiles in both the lumen and the near-wall region, and to estimate the wall-
shear stress111Wall shear stress is often sought-after in fluid mechanics
applications. with greater confidence. Based on these observations, we believe
that a general and accurate way of injecting a priori knowledge for velocity
(phase difference) regularization is in the form of a N–S boundary value
problem, and not in the form of the N–S equations alone. Several physics-
informed velocity regularization methods have been proposed in the past [4,
15, 5, 6, 7, 8], but none of them exploits the full structure of an inverse
Navier–Stokes problem where the domain boundary ($\partial\Omega$), the
boundary conditions ($\boldsymbol{g}_{i},\boldsymbol{g}_{o}$), and the
kinematic viscosity ($\nu$) are all considered unknown. A similar approach is
discussed in [9], but the method applies only to fully-sampled PC-MRI signals
and unwrapped velocity images.

Solving an inverse N–S problem amounts to finding its unknown parameters
($\boldsymbol{x}$), which produce a modeled velocity ($\boldsymbol{u}$) that
approximates the measured velocity
($\boldsymbol{u}^{\star}\propto\Delta\varphi$) in an appropriate norm. In this
way, not only do we obtain a regularized, noiseless (SNR$=\infty$) velocity
field, but we also infer hidden flow-related quantities such as the
hydrodynamic pressure ($p$) and the wall shear stresses [20, 21], which cannot
be measured using conventional MRI or PC-MRI methods, but which naturally
arise from the N–S problem. It is also important to note that, because we
jointly reconstruct and segment the velocity field, geometric errors are
minimized [22, 23].

Our approach to physics-informed CS extends the standard notion of sparsity
used in conventional CS methods to a more general notion of a structure [24],
which is dictated by the N–S problem. Instead of enforcing sparsity during the
nonlinear optimization process, we recover a sparse (hidden) structure of the
velocity field by enforcing the Navier–Stokes problem as a constraint. This is
better explained in figure 1, which shows that the velocity exhibits an
underlying low-dimensional structure that is encoded in the unknown parameters
$\boldsymbol{x}$. Based on this, a velocity field can be
compressed/decompressed by solving an inverse/forward Navier–Stokes problem.
In contrast to our approach, a conventional sparsity-promoting and physics-
informed velocity regularization method is described in [25], but it relies on
a pre-existing database (library) containing the dominant eigenmodes of a
Navier–Stokes problem that is defined in a pre-set geometry, with pre-set
inlet boundary condition.

In this paper, we build upon the work of [9] in order to formulate a physics-
informed compressed sensing (PICS) method for the joint reconstruction and
segmentation of velocity fields from sparse and noisy $k$-space signals. We
provide an algorithm that solves the reconstruction problem and demonstrate it
on $k$-space signals of a steady axisymmetric flow. Because the acquired
signals were originally fully-sampled, we sparsify them using two different
sparse sampling patterns and study the effect of sampling density on the
velocity reconstruction error. In section II we formulate the physics-informed
compressed sensing method, and an algorithm that implements it. In section III
we test the method on sparse and noisy $k$-space signals.

##  II Physics-informed compressed sensing as an inverse Navier–Stokes problem

In what follows, $L^{2}(\Omega)$ denotes the space of square-integrable
functions in $\Omega$, with inner product
$\big{\langle}\cdot,\cdot\big{\rangle}$ and norm
$\big{\lVert}\cdot\big{\rVert}_{L^{2}(\Omega)}$, and $H^{k}(\Omega)$ denotes
the space of square-integrable functions with $k$ square-integrable
derivatives. For a given covariance operator, $\mathcal{C}$, we also define
the covariance-weighted $L^{2}(\Omega)$ spaces, endowed with the inner product
${\big{\langle}\cdot,\cdot\big{\rangle}_{\mathcal{C}(\Omega)}:=\big{\langle}\cdot,\mathcal{C}^{-1}\cdot\big{\rangle}}_{\Omega}$,
which generates the norm
$\big{\lVert}\cdot\big{\rVert}_{\mathcal{C}(\Omega)}$. The Euclidean norm in
the space of real numbers $\mathbb{R}^{n}$ is denoted by
$\lvert\cdot\rvert_{\mathbb{R}^{n}}$, and the measure (volume) of the domain
$\Omega$ by $\lvert\Omega\rvert$. The first variation of a functional
$\mathscr{J}:L^{2}\to\mathbb{R}$ with respect to an unknown $z\in L^{2}$ is
defined by

$\delta_{z}\mathscr{J}\equiv\frac{d}{d\tau}\mathscr{J}({z}+\tau{z}^{\prime},\dots)\Big{|}_{\tau=0}=\big{\langle}D_{z}\mathscr{J},z^{\prime}\big{\rangle}\quad,$
(1)

where $\tau\in\mathbb{R}$, $z^{\prime}\in L^{2}$ is an allowed perturbation of
$z$, and $D_{z}\mathscr{J}$ is the generalized gradient. If $z$ is defined on
a covariance-weighted $L^{2}$ space, we furthermore define the steepest
descent direction $\widehat{D}_{z}\mathscr{J}$ such that

$\big{\langle}\widehat{D}_{z}\mathscr{J},z^{\prime}\big{\rangle}_{\mathcal{C}}=\big{\langle}\mathcal{C}^{-1}\widehat{D}_{z}\mathscr{J},z^{\prime}\big{\rangle}=\big{\langle}{D}_{z}\mathscr{J},z^{\prime}\big{\rangle}\quad,$
(2)

therefore $\widehat{D}_{z}\mathscr{J}=\mathcal{C}{D}_{z}\mathscr{J}$. When the
covariance takes the form $\mathcal{C}=\sigma^{2}\mathrm{I}$, where
$\sigma^{2}\in\mathbb{R}$ is the variance and $\mathrm{I}$ is the identity
operator, we write $\sigma^{-2}\big{\lVert}\cdot\big{\rVert}_{L^{2}}$ instead
of $\big{\lVert}\cdot\big{\rVert}_{\mathcal{C}}$ for simplicity. We often use
the superscript $(\cdot)^{\star}$ to denote a measurement, $(\cdot)^{\circ}$
to denote a reconstruction, and $(\cdot)^{\bullet}$ to denote the ground
truth. Note that, for the velocity, $\boldsymbol{u}^{\star}$ denotes the
velocity obtained from the phase differences, and $\boldsymbol{u}$ denotes the
velocity obtained from the Navier–Stokes problem.

###  II-A Phase contrast magnetic resonance imaging

Phase contrast magnetic resonance imaging can measure $d$-dimensional velocity
fields $\boldsymbol{u}$ inside or around an object $\Omega$, by requiring a
minimum of $d$ sets of $k$-space signals $\\{\boldsymbol{s}\\}_{i=1}^{d}$, one
for each velocity component. The velocity component $u_{i}$ can be recovered
from the signal set $\boldsymbol{s}_{i}=\\{s\\}_{j=1}^{4}$ by computing the
phase difference

$u_{i}=c_{i}(\varphi_{1}-\varphi_{2}-\varphi_{3}+\varphi_{4})\equiv
c_{i}(\Delta\varphi)_{i}\quad,$ (3)

where $c_{i}$ is a known constant that depends on the gyromagnetic ratio of
hydrogen and the gradient pulse properties, and
$\varphi_{j}\equiv\textrm{arg}\big{(}\mathcal{F}^{-1}s_{j}\big{)}$ when
$s_{j}$ is a fully-sampled $k$-space signal. Note that the last two phases in
the above phase difference correspond to zero-flow experiments in order to
remove any phase shift contributions that are not caused by the flow.

For sparsely-sampled $k$-space signals we define the sparse sampling operator
$\mathcal{P}:I\to I_{s}$, which projects from the full image space
$I\in\mathbb{C}^{d}$ to the sparse image space $I_{s}\subset I$. Then, if
$s_{j}$ is sparse, the corresponding complex image is given by

$w_{j}\equiv\rho_{j}e^{i\varphi_{j}}=\mathcal{F}^{-1}\mathcal{P}^{-1}s_{j}\quad.$
(4)

Since the backprojection $\mathcal{P}^{-1}$ is ill-posed, we define
$\mathcal{P}^{-1}s_{j}$ as the zero-filled $k$-space signal in $I$. We also
obtain the zero-filling solution for the velocity by using equation (3) with
$\varphi_{j}=\textrm{arg}(w_{j})$. As was mentioned in section I, the zero-
filling velocity solution is corrupted by interference noise and artefacts
that strongly depend on the sparse sampling pattern $\mathcal{P}$. We
therefore need to reconstruct the complex images $w_{j}$ using a different
approach.

###  II-B Physics-informed compressed sensing (PICS) formulation

####  II-B1 Phase regularization

A priori knowledge of the velocity, i.e. the phase difference, comes in the
form of a Navier–Stokes boundary value problem. Assuming steady,
incompressible flow and a Newtonian fluid, the fluid dynamics are governed by
the following problem

$\left\\{\begin{aligned}
\boldsymbol{u}\boldsymbol{\cdot}\nabla\boldsymbol{u}-\nu{\Delta}\boldsymbol{u}+\nabla
p&=\boldsymbol{0}\quad&&\textrm{in}\quad\Omega\\\
\nabla\boldsymbol{\cdot}\boldsymbol{u}&=0\quad&&\textrm{in}\quad\Omega\\\
\boldsymbol{u}&=\boldsymbol{0}\quad&&\textrm{on}\quad\Gamma\\\
\boldsymbol{u}&=\boldsymbol{g}_{i}\quad&&\textrm{on}\quad\Gamma_{i}\\\
-\nu\partial_{\boldsymbol{\nu}}\boldsymbol{u}+p\boldsymbol{\nu}&=\boldsymbol{g}_{o}\quad&&\textrm{on}\quad\Gamma_{o}\end{aligned}\right.\quad,$
(5)

where $\boldsymbol{u}$ is the velocity,
$p\mathrel{\reflectbox{$\mapsto$}}p/\rho$ is the reduced hydrodynamic
pressure222From now on we refer to $p$ simply as the pressure., $\rho$ is the
density, $\nu$ is the kinematic viscosity, $\boldsymbol{g}_{i}$ is the
Dirichlet boundary condition at the inlet $\Gamma_{i}$, $\boldsymbol{g}_{o}$
is the natural boundary condition at the outlet $\Gamma_{o}$,
$\boldsymbol{\nu}$ is the unit normal vector on $\partial\Omega$, and
${\partial_{\boldsymbol{\nu}}\equiv\boldsymbol{\nu}\boldsymbol{\cdot}\nabla}$
is the normal derivative. We furthermore define the projection operator
$\mathcal{S}:M\to I$ that projects from the model space $M$ to the image space
$I$. The model space $M$ is a fixed space that always contains the deformable
domain $\Omega$ ($\Omega\subset M$), such that
$\Gamma_{i},\Gamma_{o}\subset\partial M$, and

$\displaystyle\partial M\cap\textrm{Int}\ \Gamma=\varnothing,\ \partial
M\cap\Gamma_{i}=\partial M_{i},\ \partial M\cap\Gamma_{o}=\partial M_{o},$ (6)

where Int denotes the interior of a set, i.e. $\textrm{Int}\
\Gamma\equiv\big{(}\Gamma-\partial\Gamma\big{)}$. The domain $\Omega$ does not
have to be path-connected. In general, it can be defined as the union of
disjoint, path-connected sets $\Omega_{j}$, i.e.
$\Omega=\Omega_{1}\cup\Omega_{2}\cup\dots$, provided that every $\Omega_{j}$
includes at least one inlet and one outlet.

The discrepancy between the measured ($\boldsymbol{u}^{\star}$) and the
modeled ($\boldsymbol{u}$) velocity fields is given by

$\mathscr{U}(\boldsymbol{\varphi},\boldsymbol{u})\equiv\frac{1}{2}\sum^{d}_{i=1}\big{\lVert}{u}^{\star}_{i}-\mathcal{S}u_{i}\big{\rVert}^{2}_{\mathcal{C}_{u_{i}}(I_{w})}\quad,$
(7)

where $\boldsymbol{\varphi}$ contains the phases needed to compute
$\boldsymbol{u}^{\star}=\boldsymbol{u}^{\star}(\boldsymbol{\varphi})$ using
(3), $I_{w}\subseteq I$ is a user-selected area (window) of interest, and
$\mathcal{C}_{u_{i}}$ is the covariance operator for the velocity discrepancy
of the $i$-th component. In [9] we had assumed a diagonal covariance operator
for the velocity discrepancy because the noise in the phase images can be
assumed to be white and Gaussian for fully-sampled $k$-space signals with
$\text{SNR}>3$ [1]. Here, we model the interference (correlated) noise using
an exponential covariance operator such that

$\mathcal{C}_{u_{i}}g\equiv\sigma^{2}_{u_{i}}\big{(}C*g\big{)}\quad\text{for
any}\quad g\in L^{2}(I)\quad,$ (8)

where $\sigma_{u_{i}}$ is the standard deviation of noise in the $i$-th
velocity component, ‘$*$’ denotes convolution, and $C$ is the exponential
kernel

$C(r)=\frac{1}{d!V_{d}(\ell)}\ e^{-\lvert r\rvert/\ell}\quad\text{for
any}\quad r\in\mathbb{R}^{d}\quad,$ (9)

where $d$ is the dimension of the velocity field $\boldsymbol{u}$, $\ell$ is
the characteristic length, $V_{d}(\ell)=\pi^{d/2}\ell^{d}/\Gamma(d/2+1)$ is
the volume of the $d$-dimensional Euclidean ball of radius $\ell$, and
$\Gamma$ is the gamma function. Adopting a Bayesian inference setting similar
to [9], we introduce a $2\pi$-periodic prior for the phase [16]. The combined
phase regularization functional is then given by

$\displaystyle\mathscr{R}_{\varphi}(\boldsymbol{\varphi},\boldsymbol{u})\equiv\mathscr{U}+\frac{1}{2}\sum^{4d}_{j=1}\big{\lVert}e^{i\varphi_{j}}-e^{i\bar{\varphi}_{j}}\big{\rVert}^{2}_{\mathcal{C}_{\varphi_{j}}(I)}\quad,$
(10)

where $\bar{\varphi}_{j}$ is the prior assumption of $\varphi_{j}$,
$\mathcal{C}_{\varphi_{j}}=\widetilde{\sigma}_{{\varphi_{j}}}^{2}(C*\cdot)$,
$\widetilde{\sigma}_{{\varphi_{j}}}=\xi_{\varphi}\sigma_{{\varphi_{j}}}$, and
$\xi_{\varphi}\in\mathbb{R}$ is a user-selected parameter that determines the
confidence in the prior assumption.

####  II-B2 Navier–Stokes parameters regularization

Additional regularization is required to penalize improbable solutions of the
Navier–Stokes parameters $\boldsymbol{x}$. The parameters $\boldsymbol{x}$
consist of the shape of the domain $\partial\Omega$, the Dirichlet boundary
condition at the inlet(s) $\boldsymbol{g}_{i}$, the natural boundary condition
at the outlet(s) $\boldsymbol{g}_{o}$, and the kinematic viscosity $\nu$. Note
that we identify the object $\Omega$ and its boundary with a signed distance
function ${\phi_{\pm}}\in L^{2}(M)$ such that

$\displaystyle\Omega=\big{\\{}y\in\Omega:\ {\phi_{\pm}}(y)<0\big{\\}}\ ,\
\partial\Omega=\big{\\{}y\in\Omega:\ {\phi_{\pm}}(y)=0\big{\\}}.$

The regularization functional for $\boldsymbol{x}$ is given by [9]

$\displaystyle\mathscr{R}_{x}(\boldsymbol{x})=$
$\displaystyle\phantom{+}\frac{1}{2\sigma_{\phi_{\pm}}^{2}}\big{\lVert}\bar{\phi}_{\pm}-{\phi_{\pm}}\big{\rVert}^{2}_{L^{2}(M)}+\frac{1}{2}\big{\lVert}\boldsymbol{g}_{i}-\bar{\boldsymbol{g}}_{i}\big{\rVert}^{2}_{\mathcal{C}_{\boldsymbol{g}_{i}}(\partial
M_{i})}$
$\displaystyle+\frac{1}{2}\big{\lVert}\boldsymbol{g}_{o}-\bar{\boldsymbol{g}}_{o}\big{\rVert}^{2}_{\mathcal{C}_{\boldsymbol{g}_{o}}(\partial
M_{o})}+\frac{1}{2\sigma_{\nu}^{2}}\big{|}\nu-\bar{\nu}\big{|}^{2}_{\mathbb{R}}\quad,$
(11)

where $(\bar{\cdot})$ is a prior assumption, $\sigma_{(\cdot)}$ is the
standard deviation of the prior, and
$\mathcal{C}_{\boldsymbol{g}_{i}},\mathcal{C}_{\boldsymbol{g}_{o}}$ are
covariance operators that are based on the exponential covariance kernel (9),
but with additional boundary conditions, which are described in [9].

####  II-B3 Magnitude segmentation

To exploit information about $\partial\Omega$ from the nuclear spin density
(magnitude) images we use an energy-based segmentation functional
$\mathscr{S}$, which assumes that the image consists of two regions with
approximately uniform magnitudes (e.g. stars in a night sky) [26, 27]. The
functional is given by

$\displaystyle\mathscr{S}(\boldsymbol{\rho},{\phi_{\pm}},\alpha,\beta)$
$\displaystyle\equiv\frac{1}{8d}\sum^{4d}_{j=1}\frac{1}{\sigma^{2}_{\rho_{j}}}\bigg{(}\big{\lVert}\big{(}\rho_{j}-\alpha\big{)}\mathcal{S}\mathcal{H}({\phi_{\pm}})\big{\rVert}^{2}_{L^{2}(I_{w})}$
$\displaystyle+$
$\displaystyle\big{\lVert}\big{(}\rho_{j}-\beta\big{)}\big{(}\mathcal{S}\mathcal{H}({\phi_{\pm}})-1\big{)}\big{\rVert}^{2}_{L^{2}(I_{w})}\bigg{)}\quad,$
(12)

where the mean value of the average magnitude inside $\Omega$ is
$\alpha\in\mathbb{R}$ and outside $\Omega$ is $\beta\in\mathbb{R}$. In
addition, $\sigma_{\rho_{j}}$ is the standard deviation of noise in the
magnitude image, and $\mathcal{H}$ denotes the (modified) Heaviside function
such that

$\mathcal{H}(t)=1\quad\text{if}\quad t<0\quad\text{else}\quad 0\quad.$ (13)

Taking into account the prior contribution, the combined functional is given
by

$\displaystyle\mathscr{R}_{\rho}\big{(}\boldsymbol{\rho},{\phi_{\pm}},\alpha,\beta\big{)}\equiv\mathscr{S}+\frac{1}{2}\sum^{4d}_{j=1}\big{\lVert}\rho_{j}-\bar{\rho}_{j}\big{\rVert}^{2}_{\mathcal{C}_{\rho_{j}}(I)}\quad,$
(14)

where $\mathcal{C}_{\rho_{j}}=\widetilde{\sigma}_{{\rho_{j}}}^{2}(C*\cdot)$,
$\widetilde{\sigma}_{{\rho_{j}}}=\xi_{\rho}\sigma_{{\rho_{j}}}$, and
$\xi_{\rho}\in\mathbb{R}$ is a user-selected parameter that determines the
confidence in the prior assumption.

####  II-B4 Consistency with measured $k$-space signal

We consider the noise in the sparse $k$-space signals $s^{\star}_{j}$ to be
white and Gaussian, with zero mean and standard deviation $\sigma_{j}$. The
discrepancy between the measured sparse signals $s^{\star}_{j}$ and a signal
$s_{j}\equiv\mathcal{F}(\rho_{j}e^{i\varphi_{j}})$ is thus given by

$\displaystyle\mathscr{E}(\boldsymbol{\varphi},\boldsymbol{\rho})\equiv\sum^{4d}_{j=1}\frac{1}{2\sigma_{j}^{2}}\big{\lVert}s^{\star}_{j}-\mathcal{P}\mathcal{F}(\rho_{j}e^{i\varphi_{j}})\big{\rVert}^{2}_{L^{2}(I_{s})}\quad,$
(15)

with $\rho_{j}e^{i\varphi_{j}}\equiv w$ being the complex image.

Figure 2: A set of sparse $k$-space signals ($\boldsymbol{s}^{\star}$) can be
zero-filled and Fourier-transformed to the physical space in order to produce
a set of complex images ($\boldsymbol{w}^{\star}$), which provide a first
estimate for the measured velocity field ($\boldsymbol{u}^{\star}$) and the 1H
spin density inside an object. A Navier–Stokes problem is then used to
reconstruct the measured velocity field and the boundary of the object. The
reconstructed complex image ($\boldsymbol{w}^{\circ}$) is then assembled and
transformed back to $k$-space ($\boldsymbol{s}^{\circ}$), where it can be
compared with the original sparse signal.

####  II-B5 The inverse problem

We now collect the above results to formulate an inverse problem where the
reconstructed phases $\boldsymbol{\varphi}^{\circ}$, N–S parameters
$\boldsymbol{x}^{\circ}$, magnitudes $\boldsymbol{\rho}^{\circ}$, and
segmentation constants
$\boldsymbol{c}^{\circ}=(\alpha^{\circ},\beta^{\circ})$, are the minimizers of
the following nonlinearly-constrained optimization problem

$\displaystyle\big{(}\boldsymbol{\varphi}^{\circ},\boldsymbol{x}^{\circ},\boldsymbol{\rho}^{\circ},$
$\displaystyle\
\boldsymbol{c}^{\circ}\big{)}=\underset{\boldsymbol{\varphi},\boldsymbol{x},\boldsymbol{\rho},\boldsymbol{c}}{\textrm{argmin}}\
\big{(}\mathscr{R}_{\varphi}+\mathscr{R}_{x}+\mathscr{R}_{\rho}\big{)}$ (16)
subject to
$\displaystyle\text{i)}\quad\big{\lVert}s^{\star}_{j}-\mathcal{P}\mathcal{F}\big{(}\rho_{j}e^{i\varphi_{j}}\big{)}\big{\rVert}_{L^{2}(I_{s})}\leq\lvert
I_{s}\rvert\varepsilon_{j}\ \text{for all}\ j$
$\displaystyle\text{ii)}\quad\boldsymbol{u}=\mathcal{Z}\ \boldsymbol{x}\quad,$

where $\varepsilon_{j}\propto\sigma_{j}$ is a user-selected value proportional
to the standard deviation of noise. We observe that the N–S problem (5),
encoded in the operator $\mathcal{Z}$, acts as a nonlinear equality constraint
in order to ensure that $\boldsymbol{u}$ is always a N–S solution, which is
uniquely defined by the parameters $\boldsymbol{x}$
($\boldsymbol{u}^{\circ}=\mathcal{Z}\boldsymbol{x}^{\circ}$).

####  II-B6 Relaxing the constraints

We recast problem (16) into a variational form by defining the Lagrangian

$\displaystyle\mathscr{J}\equiv\mathscr{R}_{\varphi}+\mathscr{R}_{x}+\mathscr{R}_{\rho}+\mathscr{E}+\mathscr{M}\quad,$
(17)

where $\mathscr{M}$ is the weak formulation of the N–S problem (5), for test
functions $(\boldsymbol{v},q)\in\boldsymbol{H}^{1}(\Omega)\times
L^{2}(\Omega)$, given by

$\displaystyle\mathscr{M}(\boldsymbol{u},p$
$\displaystyle,\boldsymbol{v},q;\boldsymbol{x})\equiv\int_{\Omega}\Big{(}\boldsymbol{v}\boldsymbol{\cdot}\big{(}\boldsymbol{u}\boldsymbol{\cdot}\nabla\boldsymbol{u}\big{)}+\nu\nabla\boldsymbol{v}\boldsymbol{:}\nabla\boldsymbol{u}-(\nabla\boldsymbol{\cdot}\boldsymbol{v})p$
$\displaystyle-q(\nabla\boldsymbol{\cdot}\boldsymbol{u})\Big{)}+\int_{\Gamma_{o}}\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{g}_{o}+\int_{\Gamma\cup\Gamma_{i}}\boldsymbol{v}\boldsymbol{\cdot}(-\nu\partial_{\boldsymbol{\nu}}\boldsymbol{u}+p\boldsymbol{\nu})$
$\displaystyle+\mathscr{N}_{\Gamma_{i}}(\boldsymbol{v},q,\boldsymbol{u};\boldsymbol{g}_{i})+\mathscr{N}_{\Gamma}(\boldsymbol{v},q,\boldsymbol{u};\boldsymbol{0})\quad,$
(18)

and $\mathscr{N}$ is the Nitsche penalty term [28]

$\displaystyle\mathscr{N}_{\Gamma}(\boldsymbol{v},q,\boldsymbol{u};\boldsymbol{g})\equiv\int_{\Gamma}(-\nu\partial_{\boldsymbol{\nu}}\boldsymbol{v}+q\boldsymbol{\nu}+\eta\boldsymbol{v})\boldsymbol{\cdot}(\boldsymbol{u}-\boldsymbol{g})\quad,$
(19)

which we use to weakly impose a Dirichlet boundary condition
$\boldsymbol{g}\in\boldsymbol{L}^{2}(\Gamma)$ on a boundary $\Gamma$, for
fixed penalty $\eta\in\mathbb{R}$ [9, Section 2.7].

####  II-B7 Euler–Lagrange equations

Instead of solving (16) as is, it is preferable to search for critical points
of (17) by solving the nonlinear Euler–Lagrange (E–L) system for the unknowns

$\boldsymbol{v},\ q,\ \boldsymbol{u},\ p,\ \boldsymbol{\varphi},\
\boldsymbol{x},\ \boldsymbol{\rho},\ \boldsymbol{c}\ ,$

where $\boldsymbol{v}$ is the adjoint velocity and $q$ is the adjoint
pressure, acting as Lagrange multipliers in (17) via $\mathscr{M}$. We find
that the E–L equations

$\delta_{\boldsymbol{v}}\mathscr{J}=\delta_{\boldsymbol{v}}\mathscr{M}=\boldsymbol{0}\quad,\quad\delta_{q}\mathscr{J}=\delta_{q}\mathscr{M}=0\quad,$
(20)

are in fact the Navier–Stokes equations and the weakly-enforced boundary
conditions, i.e. problem (5). The next set of equations, namely

$\delta_{\boldsymbol{u}}\mathscr{J}=\delta_{\boldsymbol{u}}\big{(}\mathscr{U}+\mathscr{M}\big{)}=\boldsymbol{0}\quad,\quad\delta_{p}\mathscr{J}=\delta_{p}\mathscr{M}=0\quad,$
(21)

gives rise to the adjoint Navier–Stokes problem [9]

$\displaystyle\left\\{\begin{aligned}
-\boldsymbol{u}\boldsymbol{\cdot}\big{(}\nabla\boldsymbol{v}+(\nabla\boldsymbol{v})^{\dagger}\big{)}-\nu\Delta\boldsymbol{v}+\nabla
q&=-D_{\boldsymbol{u}}\mathscr{J}\quad&&\textrm{in}\quad\Omega\\\
\nabla\boldsymbol{\cdot}\boldsymbol{v}&=0\quad&&\textrm{in}\quad\Omega\\\
\boldsymbol{v}&=\boldsymbol{0}\quad&&\textrm{on}\quad\Gamma\cup\Gamma_{i}\\\
{(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nu})\boldsymbol{v}}+{(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{v})\boldsymbol{\nu}}+\nu\partial_{\boldsymbol{\nu}}\boldsymbol{v}-q\boldsymbol{\nu}&=\boldsymbol{0}\quad&&\textrm{on}\quad\Gamma_{o}\end{aligned}\right.$
(22)

Note that

$D_{\boldsymbol{u}}\mathscr{J}=D_{\boldsymbol{u}}\mathscr{U}=-\mathcal{S}^{\dagger}\mathcal{C}_{\boldsymbol{u}}^{-1}\big{(}\boldsymbol{u}^{\star}-\mathcal{S}\boldsymbol{u}\big{)}\quad,$
(23)

and that both $\boldsymbol{v}$ and $q$ vanish when
$D_{\boldsymbol{u}}\mathscr{U}\equiv\boldsymbol{0}$ in $I_{w}$, i.e. when the
measured velocity is identical to the modeled velocity in the user-selected
area of interest $I_{w}$. The rest of the equations comprising the E–L system
are

$\displaystyle\delta_{\boldsymbol{\varphi}}\mathscr{J}$
$\displaystyle=\delta_{\boldsymbol{\varphi}}\big{(}\mathscr{R}_{\varphi}+\mathscr{E}\big{)}$
$\displaystyle=\boldsymbol{0}$ (24)
$\displaystyle\delta_{\boldsymbol{x}}\mathscr{J}$
$\displaystyle=\delta_{\boldsymbol{x}}\big{(}\mathscr{R}_{x}+\mathscr{S}+\mathscr{M}\big{)}$
$\displaystyle=\boldsymbol{0}$ (25)
$\displaystyle\delta_{\boldsymbol{\rho}}\mathscr{J}$
$\displaystyle=\delta_{\boldsymbol{\rho}}\big{(}\mathscr{R}_{\rho}+\mathscr{E}\big{)}$
$\displaystyle=\boldsymbol{0}$ (26)
$\displaystyle\delta_{\boldsymbol{c}}\mathscr{J}$
$\displaystyle=\delta_{\boldsymbol{c}}\mathscr{S}$
$\displaystyle=\boldsymbol{0}\quad,$ (27)

which we expand in Appendix A to obtain an explicit formula for the steepest
ascent direction of each unknown. We then use these directions in order to
find a critical point of $\mathscr{J}$ using an optimization algorithm.

####  II-B8 Solving the E–L system

For large scale problems such as (16) it is prohibitively expensive to solve
the full, implicit E–L system. Instead, we devise a segregated method that
solves the E–L system, which is briefly explained in algorithm 1. Each
iteration in the main loop of algorithm 1 consists of four main stages. First,
since
$\mathscr{J}=\mathscr{J}(\boldsymbol{\varphi},\boldsymbol{x},\boldsymbol{\rho},\boldsymbol{c})$,
we fix the phases $\boldsymbol{\varphi}$, the magnitudes $\boldsymbol{\rho}$,
and the segmentation constants $\boldsymbol{c}$, and update the unknown N–S
parameters $\boldsymbol{x}$ by solving one iteration of the inverse N–S
problem. This stage is actually identical to the algorithm presented in [9],
which was devised for the reconstruction of noisy velocity images (fully-
sampled PC-MRI signals), where $\boldsymbol{\varphi},\boldsymbol{\rho}$ are
always constant, and the segmentation functional $\mathscr{S}$ is absent, i.e.
$\boldsymbol{c}$ is not involved. The updated N–S parameters
$\boldsymbol{x}_{k+1}$ are then used to update the modeled velocity field to
$\boldsymbol{u}_{k+1}=\mathcal{Z}^{\boldsymbol{u}_{k}}\boldsymbol{x}_{k+1}$,
where $\mathcal{Z}^{\boldsymbol{u}_{k}}$ denotes the Oseen linearization of
the N–S problem (5) around $\boldsymbol{u}_{k}$. The Oseen linearization
around $\boldsymbol{u}_{k}$ replaces the nonlinear convective term
$\boldsymbol{u}\boldsymbol{\cdot}\nabla\boldsymbol{u}$ in (5) with the linear
convective term $\boldsymbol{u}_{k}\boldsymbol{\cdot}\nabla\boldsymbol{u}$.
Next, we fix $\boldsymbol{x},\boldsymbol{\rho}$ and $\boldsymbol{c}$ and
reconstruct the phases, drawing information from the modeled velocity field
$\boldsymbol{u}_{k+1}$, a priori information, and the $k$-space signals.
Finally, we fix $\boldsymbol{\varphi},\boldsymbol{x}$ and $\boldsymbol{\rho}$
in order to compute the segmentation constants $\alpha,\beta$, given by the
explicit formulas (43),(44), and then fix
$\boldsymbol{\varphi},\boldsymbol{x}$ and $\boldsymbol{c}$ in order to
reconstruct the magnitudes, drawing information from the energy-based
segmentation functional and the $k$-space signals. During the first stage of
each iteration, we further update an inverse Hessian approximation of the
unknown parameters $\boldsymbol{x}$ using the damped BFGS quasi-Newton method
[29, 30], and this allows us to estimate the uncertainty in the predicted
shape, the boundary conditions, and the kinematic viscosity [9]. Algorithm 1
terminates if either the update for the unknowns
$\boldsymbol{\varphi},\boldsymbol{x},\boldsymbol{\rho},\boldsymbol{c}$ is
below the user-specified tolerance, or the line search fails to reduce
$\mathscr{J}$. Every line search starts with a step size of $\tau=1$, and the
step size is halved until $\mathscr{J}_{k+1}<\mathscr{J}_{k}$.

Input: sparse PC-MRI signals $\boldsymbol{s}^{\star}$, noise variances,
initial guesses (priors) for the unknowns
$\bar{\boldsymbol{\varphi}},\bar{\boldsymbol{x}},\bar{\boldsymbol{\rho}}$

begin

$k\leftarrow 0$

Initialization

$\boldsymbol{u}^{\star}_{0}$banana $\leftarrow$ zero-filled velocity from
$\bar{\boldsymbol{\varphi}}$ \- eq. (3)

$\bar{\phi}_{\pm}$banana $\leftarrow$ signed distance field from
$\sum_{j}\bar{\rho}_{j}$ [9]

$(\boldsymbol{u},p)_{0}$banana $\leftarrow$ $\mathcal{Z}\
\bar{\boldsymbol{x}}$ \- N–S problem (5)

_ _ while  _convergence_criterion_is_not_met_

 _[1.5pt]_ do

inverse N–S - min
$\mathscr{J}({\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\boldsymbol{\varphi}_{k}},\boldsymbol{x},{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\boldsymbol{\rho}_{k}},{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\boldsymbol{c}_{k}})$:

$(\boldsymbol{v},q)_{k}$avocados $\leftarrow$ adjoint N–S problem (22)

$\widehat{D}_{\boldsymbol{x}}\mathscr{J}$avocados $\leftarrow$ s. a.
directions (45)-(51)

$\widetilde{H}^{k}_{x},\tau$avocados $\leftarrow$ approx. inv. Hessian and l.
s. [9]

$\boldsymbol{x}_{k+1}$avocados $\leftarrow$ $\boldsymbol{x}_{k}-\tau\
\widetilde{H}^{k}_{x}\widehat{D}_{\boldsymbol{x}}\mathscr{J}$

$(\boldsymbol{u},p)_{k+1}$avocados $\leftarrow$
$\mathcal{Z}^{\boldsymbol{u}_{k}}\boldsymbol{x}_{k+1}$ \- lin. N–S problem (5)

[3pt]  phases - min
$\mathscr{J}(\boldsymbol{\varphi},{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\boldsymbol{x}_{k+1}},{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\boldsymbol{\rho}_{k}},{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\boldsymbol{c}_{k}})$:

$\widehat{D}_{\boldsymbol{\varphi}}\mathscr{J},\tau$avocados $\leftarrow$ s.
a. directions (39) and l. s.

$\boldsymbol{\varphi}_{k+1}$avocados $\leftarrow$
$\boldsymbol{\varphi}_{k}-\tau\ \widehat{D}_{\boldsymbol{\varphi}}\mathscr{J}$

$\boldsymbol{u}^{\star}_{k+1}$avocados $\leftarrow$ measured vel. from
$\boldsymbol{\varphi}_{k+1}$ \- eq. (3)

[3pt]  segmentation - min
$\mathscr{J}({\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\boldsymbol{\varphi}_{k+1}},{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\boldsymbol{x}_{k+1}},{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\boldsymbol{\rho}_{k}},\boldsymbol{c})$:

$\boldsymbol{c}_{k+1}$avocados $\leftarrow$ compute $(\alpha,\beta)_{k+1}$ \-
eq. (43),(44)

[3pt]  magnitudes - min
$\mathscr{J}({\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\boldsymbol{\varphi}_{k+1}},{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\boldsymbol{x}_{k+1}},\boldsymbol{\rho},{\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\boldsymbol{c}_{k+1}})$:

$\widehat{D}_{\boldsymbol{\rho}}\mathscr{J},\tau$avocados $\leftarrow$ s. a.
directions (40) and l. s.

$\boldsymbol{\rho}_{k+1}$avocados $\leftarrow$ $\boldsymbol{\rho}_{k}-\tau\
\widehat{D}_{\boldsymbol{\rho}}\mathscr{J}$

[3pt]  $k\leftarrow k+1$

Output:

$\boldsymbol{x}^{\circ}$avocados$\leftarrow$ $\boldsymbol{x}_{k}$bananas
(inferred N–S parameters)

$(\boldsymbol{u}^{\circ},p^{\circ})$avocados$\leftarrow$
$(\boldsymbol{u},p)_{k}$bananas (N–S velocity and pressure)

$(\boldsymbol{\varphi}^{\circ},\boldsymbol{\rho}^{\circ})$avocados$\leftarrow$
$(\boldsymbol{\varphi},\boldsymbol{\rho})_{k}$bananas (reconstructed
phases/mag.)

[3pt] Optional output:

$\gamma_{w}\leftarrow$ wall shear rate from $\boldsymbol{u}^{\circ}$,
$\boldsymbol{x}^{\circ}$

s. a.: steepest ascent, l. s.: line search, lin.: linearized

Algorithm 1 PICS for sparse PC-MRI signals.

####  II-B9 Numerics

To solve boundary value problems, problem (5) and (22), for example, we use an
immersed boundary finite element method, the details of which can be found in
[9, Section 2.7]. We implement algorithm 1, and all the necessary numerical
methods, in Python modules, using Python’s standard libraries for scientific
computing: NumPy [31] and SciPy [32].

##  III Sparse PC-MRI reconstruction using PICS

We now apply algorithm 1 to sparse and noisy (low SNR) PC-MRI signals of an
axisymmetric flow of water/glycerol through a converging nozzle [9]. The
description of the PC-MRI experiment can be found in Appendix B. Since the
$k$-space was fully-sampled, we sparsify the $k$-space signals using two
different sparse sampling patterns $\mathcal{P}$: i) a sampling pattern
$\mathcal{P}_{\odot}$, created from a two-dimensional normal distribution, and
ii) a sampling pattern $\mathcal{P}_{\parallel}$, created from a one-
dimensional normal distribution of lines. The sampling pattern
$\mathcal{P}_{\parallel}$ is commonly used in MRI because conventional
acquisition protocols sample $k$-space lines for each system excitation. For
these sparse sampling patterns, we further investigate the effect of
subsampling on the velocity reconstruction error. The reconstruction error of
a modeled velocity component $u^{\circ}_{i}$ is measured by

$\mathcal{E}^{\bullet}_{u^{\circ}_{i}}=(\sqrt{\lvert
M\rvert}\sigma_{u^{\bullet}_{i}})^{-2}\big{\lVert}u^{\bullet}_{i}-u^{\circ}_{i}\big{\rVert}_{L^{2}(M)}\quad,$
(28)

where $u^{\bullet}_{i}$ is the corresponding ground truth velocity image,
which, in this paper, is a high signal-to-noise ratio image of the same flow
(figure 3), and $\sigma_{u^{\bullet}_{i}}^{2}$ is the estimated variance of
Gaussian white noise in the ground truth image.

Figure 3: Ground truth for the axial velocity component $u^{\bullet}_{z}$
($\text{SNR}_{z}\simeq 52$), and the radial velocity component
$u^{\bullet}_{r}$ ($\text{SNR}_{r}\simeq 43$), for the left-to-right flow
through the converging nozzle (colormap scale shown in figure 5). The images
were obtained using phase-contrast MRI (see Appendix B).

###  III-A Noise in the signal, the phase, and the magnitude

The noise in each $k$-space signal $s^{\star}_{j}$ is assumed to be Gaussian,
white, and additive, with zero mean and variance $\sigma^{2}_{j}$. It is
common practice in MRI to estimate the SNR from regions of the magnitude
images with no signal. For fully-sampled signals, the magnitude of the signal
$s^{\star}_{j}\in\mathbb{C}^{d}$, given by
$\rho_{j}\equiv\lvert\mathcal{F}^{-1}s^{\star}_{j}\rvert\in\mathbb{R}^{d}$, is
distributed according to a Rayleigh distribution [1]. The noise variance in
$s^{\star}_{j}$ is thus given by

$\sigma^{2}_{j}=2(4-\pi)^{-1}\ \sigma^{2}_{\rho_{j}}\quad,$ (29)

where $\sigma^{2}_{\rho_{j}}$ is the noise variance of the magnitude image
$\rho_{j}$, estimated from a region with no signal. The SNR of the $j$-th
signal is then computed from

$\text{SNR}_{j}=\frac{\mu_{\rho_{j}}}{\sigma_{\rho_{j}}}\quad,$ (30)

where $\mu_{\rho_{j}}$ is the mean of $\rho_{j}$, taking into account only the
active regions (non-zero signal). When $\text{SNR}_{j}>3$, the noise in the
phase images $\varphi_{j}$ approximates a zero-mean Gaussian distribution with
variance [1]

$\sigma^{2}_{\varphi_{j}}=\text{SNR}_{j}^{-2}\quad.$ (31)

Considering that the measured velocity, given by equation (3), is computed
from the phase differences, the noise variance of the $k$-th velocity
component is estimated by

$\sigma^{2}_{u_{k}}=c^{2}_{k}\sum_{j}\sigma^{2}_{\varphi_{j}}\quad,$ (32)

for $j=1,\dots,4d$ and $k\equiv(j-1)\textrm{div}4+1$, where div denotes
integer division.

Based on the above, for fully-sampled $k$-space signals the noise
$\varepsilon_{u_{k}}$ in the velocity image of the $k$-th component is
distributed according to $\mathcal{N}(0,\sigma_{u_{k}}^{2}\mathrm{I})$. For
sparsely-sampled $k$-space signals we observe that $\varepsilon_{u_{k}}$ is
correlated. Here, we model this correlation using an exponential covariance
function (see section II-B1). Figure 4 shows the velocity discrepancy
($u^{\star}_{z}-\mathcal{S}u_{z}$) between the measured ($u^{\star}_{z}$) and
the modeled ($u_{z}$) axial velocity components for fully-sampled signals, and
15% sparsely-sampled signals using $\mathcal{P}_{\odot}$. For fully-sampled
$k$-spaces, we observe that

$u^{\star}_{k}-\mathcal{S}u_{k}\sim\mathcal{N}(0,\sigma^{2}_{u_{k}}\mathrm{I})\quad,$
(33)

and this justifies the use of the norm
$\sigma^{-2}_{u_{k}}\big{\lVert}u^{\star}_{k}-\mathcal{S}u_{k}\big{\rVert}^{2}_{L^{2}}$
in [9], where fully-sampled signals were considered. For sparsely-sampled
$k$-spaces, we observe that the noise is correlated and we assume that

$u^{\star}_{k}-\mathcal{S}u_{k}\sim\mathcal{N}(0,\mathcal{C}_{u_{k}})\quad,$
(34)

where $\mathcal{C}_{u_{k}}$ is given by equation (8), with characteristic
length $\ell$ equal to the smallest resolved length scale in the image $I$.

Figure 4: Left frames: velocity discrepancy for fully-sampled (100% $k$-space
coverage) signals. The noise $\varepsilon$ in the measured velocity can be
considered to be white and Gaussian when the signals are fully-sampled, i.e.
$u^{\star}_{k}-\mathcal{S}u_{k}\sim\mathcal{N}(0,\sigma^{2}\mathrm{I})$. Right
frames: velocity discrepancy for sparsely-sampled signals (15% $k$-space
coverage with $\mathcal{P}_{\odot}$). For sparsely-sampled signals we observe
correlations that depend on the sampling pattern $\mathcal{P}$. To model these
correlations, we introduce a covariance operator $\mathcal{C}$ that is
constructed using an exponential kernel, i.e.
$u^{\star}_{k}-\mathcal{S}u_{k}\sim\mathcal{N}(0,\mathcal{C})$. TABLE I: Noise
level for the acquired noisy signals and the resulting velocity images.

$k$-space signals $s^{\star}_{j}$
---
$\sigma_{1}$ | $\sigma_{2}$ | $\sigma_{3}$ | $\sigma_{4}$ | $\sigma_{5}$ | $\sigma_{6}$ | $\sigma_{7}$ | $\sigma_{8}$
38.4 | 38.5 | 14.7 | 14.6 | 43.1 | 42.9 | 19.2 | 19.2

vel. component |  $c_{k}$ [cm/s] |  $\sigma_{u_{(\cdot)}}$ [cm/s] | SNRj
---|---|---|---
axial $z$ ($j=1,\dots 4$) | 1.988 | 0.348 | 11.4
radial $r$ ($j=5,\dots 8$) | 0.529 | 0.113 | 9.3

###  III-B Generating the sparse sampling patterns $\mathcal{P}_{\odot}$ and
$\mathcal{P}_{\parallel}$

To generate the sampling pattern $\mathcal{P}_{\odot}$ for a discretized two-
dimensional $k$-space, we draw random samples $\boldsymbol{y}^{s}_{k}$ from a
normal distribution $\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$, whose
mean and covariance are given by

$\boldsymbol{\mu}=\Big{(}\frac{n_{1}}{2},\frac{n_{2}}{2}\Big{)}\quad,\quad\boldsymbol{\Sigma}=\textrm{diag}\Big{(}\frac{\omega_{1}n_{1}}{4},\frac{\omega_{2}n_{2}}{4}\Big{)}\quad,$
(35)

where $n_{i}$ denotes the number of points (or pixels), and $\omega_{i}$
denotes the coverage, in the sense of a $2\sigma$ interval, along the $i$-th
direction. We further round the samples
$\boldsymbol{y}^{s}_{k}=(y^{s}_{k1},y^{s}_{k2})\mapsto\lfloor\boldsymbol{y}^{s}_{k}\rceil=(\lfloor
y^{s}_{k1}\rceil,\lfloor y^{s}_{k2}\rceil)$ to the closest integer, so that
they correspond to the pixel index. The degree of subsampling is controlled by
the sampling density, which is defined by $N_{s}/N$, where $N_{s}$ is the
number of admissible sampled $k$-space points (we reject samples that either
repeat on the same pixel or lie outside the $k$-space domain) and
$N=n_{1}n_{2}$ is the total number of $k$-space points. The sampling pattern
$\mathcal{P}_{\odot}:\mathbb{Z}^{n_{1}\times n_{2}}\to\mathbb{R}$ is then
given by

$\mathcal{P}_{\odot}(\boldsymbol{y})=1\quad\text{if}\quad\boldsymbol{y}=\lfloor\boldsymbol{y}^{s}_{k}\rceil\quad\text{else}\quad
0\ ,$ (36)

for $k=1,\dots,N_{s}$. The sampling pattern
$\mathcal{P}_{\parallel}:\mathbb{Z}^{n_{1}\times n_{2}}\to\mathbb{R}$ consists
of random samples of vertical lines, and is given by

$\mathcal{P}_{\parallel}(y_{1},\cdot)=1\quad\text{if}\quad
y_{1}=\lfloor{y}^{s}_{k1}\rceil\quad\text{else}\quad 0\ ,$ (37)

for $k=1,\dots,N_{s}$.

###  III-C Reconstructing sparse PC-MRI signals

Using the sparse sampling patterns $\mathcal{P}_{\odot}$ and
$\mathcal{P}_{\parallel}$, with $\omega_{1}=\omega_{2}=0.35$, and for various
sampling densities $N_{s}/N$, we sparsify the originally full $k$-space
signals that we have acquired for the flow through the nozzle. We use
algorithm 1 to reconstruct the sparse, noisy (low SNR) signals
$\boldsymbol{s}^{\star}$ (see table I for the noise level), and compare the
results with the ground truth (high SNR) images that we have acquired for the
same flow (figure 3). The low SNR images required a total scanning time of 2.6
minutes per velocity image (axial and radial components), and the high SNR
images required a total scanning time of 68 minutes per velocity image. These
scanning times correspond to 100% $k$-space coverage. Thus, all other things
being equal, the 15%-sampled low SNR signals would have been acquired in $\sim
23$ seconds, assuming that acquisition time linearly scales with sampling
density333The validity of this assumption depends on both the sparse sampling
pattern and the pulse sequence of the magnetic resonance experiment..

####  III-C1 Input and initialization

The input set of algorithm 1 consists of the $k$-space signals
$\boldsymbol{s}^{\star}$, the estimated noise variances (table I), and the
prior assumption for $\boldsymbol{\varphi},\boldsymbol{\rho}$ and
$\boldsymbol{x}$. The prior assumptions are in fact Gaussian random fields,
and we therefore need to define both their mean and their variance (confidence
level of a priori knowledge). We choose the prior mean of the phase
($\bar{\boldsymbol{\varphi}}$) and the magnitude ($\bar{\boldsymbol{\rho}}$)
images to be their respective zero-filling solution. For the N–S unknowns, we
choose the prior mean of the signed distance function (SDF)
($\bar{\phi}_{\pm}$), to be the SDF that corresponds to the Chan–Vese
segmentation [26, 27] of the averaged zero-filling magnitude image (we show
how to compute the SDF from a segmentation in [9]). We set the prior mean for
the inlet velocity boundary condition to
$\bar{\boldsymbol{g}}_{i}=(\bar{g}_{iz},0)$, where $\bar{g}_{iz}$ is a
parabolic velocity profile that satisfies the zero-velocity boundary condition
on the boundary $\partial\Gamma_{i}$ of the inlet $\Gamma_{i}$, and has peak
velocity equal to $2.5$ cm/s. The prior mean of the outlet boundary condition
is $\bar{\boldsymbol{g}}_{o}\equiv\boldsymbol{0}$ (traction-free boundary
condition), and of the kinematic viscosity is $\bar{\nu}\simeq 2.54\times
10^{-5}$ m2/s [33, 34], which corresponds to a 70% glycerol in water mixture
(see Appendix B). We explain how we choose the prior variances in Appendix C.
Lastly, the initial guess for the measured velocity is obtained using equation
(3) with $\bar{\boldsymbol{\varphi}}$, and the initial guess for the modeled
velocity by solving the Navier–Stokes problem $\mathcal{Z}\
\bar{\boldsymbol{x}}$.

TABLE II: Reconstruction error $\mathcal{E}_{z}^{\bullet}$ / $\mathcal{E}_{r}^{\bullet}$ for each sampling pattern $\mathcal{P}$. |  $k$-space coverage
---|---
sampling | 5% | 10% | 15% | 25%
$\mathcal{P}_{\odot}$ | 1.13 / 0.47 | 0.66 / 0.42 | 0.55 / 0.30 | $\cdot$
$\mathcal{P}_{\parallel}$ | $\cdot$ | 2.41 / 1.78 | 0.63 / 0.39 | 0.61 / 0.29
full sampling | 0.56 / 0.31

####  III-C2 Reconstruction results

We first test the algorithm on sparse $k$-space signals that we generate using
$\mathcal{P}_{\odot}$ for 5%, 10% and 15% sampling. The generated sampling
patterns, the reconstructed measured velocity $\boldsymbol{u}^{\star}$
(obtained from the phase difference), and the reconstructed modeled velocity
$\boldsymbol{u}$ (N–S solution), are shown in figure 5, alongside their
corresponding zero-filling solution. The fourth column of figure 5 depicts the
result of the image reconstruction algorithm [9], which we use when the
signals are fully-sampled. For fully-sampled $k$-space signals, there is a
one-to-one correspondence between the $k$-space and the physical space, and,
consequently, we can directly reconstruct the velocity field in physical space
($\boldsymbol{u}^{\star}$ is known, but noisy, and fixed during the
reconstruction).

In the experiments we selected velocity encoding parameters that cause
aliasing due to phase wrapping, in order to increase the SNR of the velocity
images. For the 100% sampling case with this flow, we use basic knowledge of
fluid mechanics to unwrap the image by adding or subtracting multiples of
$2\pi$ to aliased pixels of the phase difference image. We can do this because
we already know that the velocity components have well-defined regions with
fixed signs. We use this simple method of unwrapping only for the 100%
sampling case, in order to be able to reconstruct it using the algorithm in
[9]. For 5%, 10% and 15% sampling, the present algorithm (PICS, algorithm 1)
unwraps and reconstructs the measured velocity simultaneously and autonomously
(figure 7). This is possible with algorithm 1 because it treats
$\boldsymbol{u}^{\star}$ as a function of the unknown phases
$\boldsymbol{\varphi}$, i.e. with $\boldsymbol{u}^{\star}$ no longer being
fixed.

The velocity reconstruction errors, defined by equation (28), are presented in
table II for $\mathcal{P}_{\odot}$-subsampled and
$\mathcal{P}_{\parallel}$-subsampled signals. We observe that for 15%
$\mathcal{P}_{\odot}$-sampling, the reconstruction errors are
$\mathcal{E}^{\bullet}_{z}=0.55$ and $\mathcal{E}^{\bullet}_{r}=0.30$, and
compare this result to the 100%-sampling reconstruction errors of
$\mathcal{E}^{\bullet}_{z}=0.56$ and $\mathcal{E}^{\bullet}_{r}=0.31$. The
fact that the $\mathcal{P}_{\odot}$-subsampled signal has a lower
reconstruction error than the 100%-sampled signal may seem counterintuitive,
but in reality the difference in the errors is negligible, and it can be
attributed to either the technical implementation details of this large-scale
problem, or in the outliers, the artefacts and the noise in the remaining 85%
of the $k$-space. For the sampling pattern $\mathcal{P}_{\odot}$, both the
reconstruction error values (table II) and a visual inspection of the velocity
images (figures 5,3) show that 15% sampling is sufficient to accurately
reconstruct the flow through the converging nozzle. Next, we test the
algorithm on sparse $k$-space signals that we generate using
$\mathcal{P}_{\parallel}$ for 10%, 15% and 25% sampling. This sampling pattern
is more coherent than $\mathcal{P}_{\odot}$, and, therefore, it produces
strong artefacts in the phase and magnitude images for low sampling densities
(figure 6), making the reconstruction problem harder to solve. For the
sampling pattern $\mathcal{P}_{\parallel}$, both the reconstruction error
values (table II) and a visual inspection of the velocity images (figures 6,3)
show that 25% sampling is sufficient in order to accurately reconstruct the
flow through the converging nozzle.

Although we treat the kinematic viscosity $\nu$ as an unknown parameter, for
all reconstructions we find that the posterior distribution of $\nu$ remains
effectively unchanged, i.e. $\nu^{\circ}\simeq\bar{\nu}$ and
$\sigma^{\circ}_{\nu}\simeq{\sigma}_{\nu}$, where $\sigma^{\circ}_{\nu}$ is
the posterior standard deviation. This is because the prior information is
accurate enough (see Appendix C, where the prior variance $\sigma^{2}_{\nu}$
is small) and the model $\mathcal{M}$ cannot further improve this prediction
because the velocity reconstruction functional $\mathcal{U}$ is insensitive to
such small changes in $\nu$ (the prior term in equation (47) dominates) [9,
Section 3.5].

Figure 5: We reconstruct sparse $k$-space signals with 5%, 10% and 15%
$\mathcal{P}_{\odot}$-sampling using algorithm 1 (first three columns), and
compare the results with the reconstruction of the 100%-sampled $k$-space
signal using the algorithm in [9] (last column). For the reconstructed sparse
signals (first three columns), we show the zero-filled velocity (rows 2,5),
the reconstructed measured velocity (rows 3,6), and the reconstructed modeled
velocity (rows 4,7). For the 100%-sampled $k$-space signal (last column), we
show the resulting velocity using Fourier inversion (row 2,5), the
postprocessed (unwrapped and masked) Fourier inversion (row 3,5) and the
reconstructed velocity (row 4,7), obtained using the algorithm in [9] (flow is
from left to right). Figure 6: As for figure 5, but for sparse $k$-space
signals with 10%, 15% and 25% $\mathcal{P}_{\parallel}$-sampling. Figure 7:
Iterates $k$ showing the synchronous velocity reconstruction and unwrapping
for the measured axial velocity component $u^{\star}_{z}$ (left), and the
measured radial velocity component $u^{\star}_{r}$ (right), starting from
their respective zero-filled solution (colormap scale show in figure 6).

###  III-D Further reducing sampling density

####  III-D1 Targeted sampling for the N–S unknowns $\boldsymbol{x}$

At this point, it is important to mention that the sampling density can be
further decreased (for fixed reconstruction error) if the a priori information
on $\boldsymbol{x}$ becomes more accurate. The present method indicates that
when the Navier–Stokes problem is an accurate and descriptive model, only the
boundary $\partial\Omega$, the boundary conditions
$\boldsymbol{g}_{i},\boldsymbol{g}_{o}$, and the kinematic viscosity are
needed in order to find the velocity field (see figure 1). Therefore, to
reduce signal acquisition time further, we can obtain more accurate spin
density (magnitude) images separately (without encoding velocity) and use them
as priors for $\partial\Omega$ ($\bar{\phi}_{\pm}$). Similarly, we can obtain
$d-1$ dimensional scans for the inlet velocity boundary condition, and use
them as priors for $\bar{\boldsymbol{g}}_{i}$. Conventional MRI or PC-MRI
cannot measure $\boldsymbol{g}_{o}$ and $\nu$, but $\nu$ is usually known with
high certainty compared to the other parameters. Targeted sampling using a
Navier–Stokes problem is an interesting generalization [24] of sparsity in CS,
and closely follows the concept of compressed sampling.

####  III-D2 Design of optimal sampling patterns $\mathcal{P}$

Another way to further decrease the sampling density is to design optimal
sampling patterns [35], also taking into account the constraints of the PC-MRI
signal acquisition protocols. For example, in section III-C2 we found that the
generated 15% $\mathcal{P}_{\odot}$-sampling pattern produces slightly better
results that the generated 25% $\mathcal{P}_{\parallel}$-sampling pattern, but
note that the former adequately covers the center of $k$-space (figure 5),
while the latter leaves gaps (figure 6). This can be corrected by fully-
sampling the central region of $k$-space and then sparsely-sampling regions
that are not as important. To identify the important regions, we can use our
algorithm to backpropagate the errors from the velocity images to the
$k$-space. The physics-informed design of optimal sparse sampling patterns is
left for future work.

###  III-E Wall shear stress estimation

Finally we compute the reconstructed wall shear rate $\gamma^{\circ}_{w}$
(figure 8) using the reconstructed (noiseless, $\text{SNR}=\infty$) velocity
field $\boldsymbol{u}^{\circ}$ ($=\mathcal{Z}\boldsymbol{x}^{\circ}$) and the
inferred shape $\partial\Omega^{\circ}$, and compare it with the ground truth
wall shear rate $\gamma^{\bullet}_{w}$, computed for the ground truth (high
SNR) velocity measurement $\boldsymbol{u}^{\bullet}$ (figure 3) on the ground
truth (high SNR) shape $\partial\Omega^{\bullet}$. We further compute the
uncertainty in $\gamma^{\circ}_{w}$ by propagating the uncertainty of the
reconstructed N–S unknowns $\boldsymbol{x}^{\circ}$ through the N–S problem.
We observe that the ground truth wall shear rate is particularly noisy because
PC-MRI suffers from low resolution and partial volume effects [36, 37] near
the boundaries $\partial\Omega$. Using standard image processing algorithms,
it is possible to either denoise $\gamma^{\bullet}_{w}$ directly, or compute
$\gamma^{\bullet}_{w}$ for a denoised velocity $\boldsymbol{u}^{\bullet}$ and
a smoothed boundary $\partial\Omega^{\bullet}$. This denoised velocity field
will not, however, be consistent with the new smoothed boundary (the no-slip
wall boundary condition will not be satisfied), and the steep, near-wall
velocity gradients that dictate the wall shear rate magnitude will not be
accurately recovered. On the other hand, our inverse Navier–Stokes boundary
value problem (algorithm 1), jointly reconstructs the velocity field and
infers the boundary. This accurately reconstructs the near-wall velocity
distribution, and, if needed, can be made more accurate by increasing the N–S
model resolution. This method infers the boundary by drawing information from
both the velocity field and the spin density (magnitude) images so that
geometric errors that affect the estimated wall shear stress distribution can
be minimized.

Figure 8: Inferred shape $\partial\Omega^{\circ}$ (first row) for 5%, 10% and
15% $\mathcal{P}_{\odot}$-sampling using algorithm 1, and for 100%-sampling
using the algorithm in [9]. The inferred shape is shown in cyan color and the
2$\sigma$ confidence interval in blue color. The yellow line denotes the
ground truth of $\partial\Omega$, computed from the segmentation of the high
SNR magnitude images. The inferred (nondimensional) reduced hydrodynamic
pressure (second row) $p^{\circ}/U^{2}$, where $U=4.39$ cm/s, is immediately
obtained from the Navier–Stokes problem without any additional computation
(see algorithm 1). The predicted wall shear rate
$\gamma^{\circ}_{w}\equiv\boldsymbol{\tau}\boldsymbol{\cdot}\partial_{\boldsymbol{\nu}}\boldsymbol{u}^{\circ}$
(third row), where $\boldsymbol{\tau}$ is the unit tangent vector of
$\partial\Omega^{\circ}$, is shown in cyan, and the 2$\sigma$ confidence
interval, estimated by propagating the uncertainties of
$\boldsymbol{x}^{\circ}$ through the Navier–Stokes problem, is shown in blue.
The red line denotes the wall shear rate computed from the high SNR velocity
data ($\gamma^{\bullet}_{w}$). The wall shear stress $\tau_{w}$ is found by
multiplying the wall shear rate by the dynamic viscosity $\mu$ ($=\rho\nu$,
where $\rho=1183.6$ kg/m3), i.e. $\tau_{w}=\mu\gamma_{w}$.

##  IV Conclusion

We propose a physics-informed compressed sensing (PICS) reconstruction method
that optimally combines computational fluid dynamics with phase-contrast
magnetic resonance imaging. By adopting a Bayesian framework, we use a priori
knowledge to regularize the inverse problem and find the most likely measured
velocity field, segmentation, and Navier–Stokes twin (maximum a posteriori
estimates) for the flow through a converging nozzle. We show that the physics-
informed reconstruction of the sparse (15% $\mathcal{P}_{\odot}$-sampling),
low SNR signals, compares well to fully-sampled, high SNR images, of the same
flow.

PICS extends the capabilities of conventional compressed sensing methods for
phase-contrast MRI. In particular, it integrates an optimal filter (the
regularity of the velocity field is controled by the Navier–Stokes problem),
it infers the hydrodynamic pressure without any additional computation, and
allows us to estimate the wall shear stresses at negligible additional
numerical cost. At the same time it provides us with a digital twin, which can
subsequently be used to simulate different flow conditions, and paves the way
for the development of new applications in patient-specific cardiovascular
modeling. Our formulation naturally extends to three-dimensional and periodic
flows in complicated geometries, and the Navier–Stokes problem can be updated
to model complex (e.g. non-Newtonian) fluid flows.

## Acknowledgments

The authors would like to thank Scott V. Elgersma and Andrew J. Sederman for
generating the PC-MRI dataset, which is the same as that of [9], and for their
useful input and comments during the preparation of this manuscript. Author
A.K. is financially supported by the W.D. Armstrong scholarship from the
Cambridge Trust.

##  Appendix A E–L equations for
$\boldsymbol{\varphi},\boldsymbol{\rho},\boldsymbol{c}$ and $\boldsymbol{x}$

Here we expand (24)-(27) in order to obtain explicit relations for the
corresponding generalized gradients ${D}_{(\cdot)}\mathscr{J}$, which are
defined by (1). Note that the steepest ascent directions, which are defined by
(2), are given by

$\widehat{D}_{(\cdot)}\mathscr{J}=\mathcal{C}_{(\cdot)}{D}_{(\cdot)}\mathscr{J}\quad.$
(38)

For the phases and magnitudes we find

$\displaystyle{D}_{{\varphi}_{j}}$
$\displaystyle\mathscr{J}=\chi_{I_{w}}(-1)^{j\textrm{div}2}c_{k}~{}\mathcal{C}_{u_{k}}^{-1}\big{(}u^{\star}_{k}-\mathcal{S}u_{k}\big{)}$
$\displaystyle+\mathcal{C}^{-1}_{\varphi_{j}}\Re\Big{(}ie^{i\varphi_{j}}\big{(}e^{i\varphi_{j}}-e^{i\bar{\varphi}_{j}}\big{)}^{\dagger}\Big{)}$
$\displaystyle+\sigma_{j}^{-2}\rho_{j}~{}\Im\Big{(}e^{i\varphi_{j}}\big{(}\mathcal{F}^{-1}\mathcal{P}^{\dagger}(s^{\star}_{j}-\mathcal{P}\mathcal{F}\rho_{j}e^{i\varphi_{j}})\big{)}^{\dagger}\Big{)}\
,$ (39) $\displaystyle{D}_{{\rho}_{j}}$
$\displaystyle\mathscr{J}=\chi_{I_{w}}(\sqrt{4d}\
\sigma_{\rho_{j}})^{-2}\Big{(}\big{(}\rho_{j}-\alpha\big{)}\mathcal{S}\mathcal{H}({\phi_{\pm}})$
$\displaystyle+\big{(}\rho_{j}-\beta\big{)}\big{(}\mathcal{S}\mathcal{H}({\phi_{\pm}})-1\big{)}\Big{)}+\mathcal{C}^{-1}_{\rho_{j}}\big{(}\rho_{j}-\bar{\rho}_{j})$
$\displaystyle-\sigma_{j}^{-2}~{}\Re\Big{(}e^{i\varphi_{j}}\big{(}\mathcal{F}^{-1}\mathcal{P}^{\dagger}(s^{\star}_{j}-\mathcal{P}\mathcal{F}\rho_{j}e^{i\varphi_{j}})\big{)}^{\dagger}\Big{)}\
,$ (40)

for $j=1,\dots,4d$, where $k\equiv(j-1)\textrm{div}4+1$, $\chi_{I_{w}}$ is the
characteristic function of $I_{w}\subseteq I$. Note that
${D}_{{\boldsymbol{\varphi}}}\mathscr{J}(\boldsymbol{\varphi},\boldsymbol{u},\boldsymbol{\rho})\in
L^{2}(I)$, and
${D}_{{\boldsymbol{\rho}}}\mathscr{J}(\boldsymbol{\rho},{\phi_{\pm}},\boldsymbol{\varphi},\alpha,\beta)\in
L^{2}(I)$.

For the segmentation constants we find

$\displaystyle{D}_{\alpha}\mathscr{J}$
$\displaystyle=-\sum_{j=1}^{4d}(\sqrt{4d}\
\sigma_{\rho_{j}})^{-2}\int_{\Omega}\rho_{j}-\alpha\quad,$ (41)
$\displaystyle{D}_{\beta}\mathscr{J}$
$\displaystyle=-\sum_{j=1}^{4d}(\sqrt{4d}\
\sigma_{\rho_{j}})^{-2}\int_{I_{w}-\Omega}\rho_{j}-\beta\quad,$ (42)

where $I_{w}-\Omega$ denotes the complement of $\Omega$ in $I_{w}$. Note that
${D}_{{\alpha}}\mathscr{J}(\boldsymbol{\rho},{\phi_{\pm}},\alpha)\in\mathbb{R}$,
and
${D}_{\beta}\mathscr{J}(\boldsymbol{\rho},{\phi_{\pm}},\beta)\in\mathbb{R}$.
Fixing $\boldsymbol{\rho}$ and $\Omega$ (or, equivalently, ${\phi_{\pm}}$),
and setting (41) and (42) to zero, we find the following explicit relations
for the constants $\alpha,\beta$ that minimize $\mathscr{J}$

$\displaystyle\alpha$
$\displaystyle=\sum_{j=1}^{4d}(\sqrt{4d\lvert\Omega\rvert}\sigma_{\rho_{j}})^{-2}\int_{\Omega}\rho_{j}\quad,$
(43) $\displaystyle\beta$ $\displaystyle=\sum_{j=1}^{4d}(\sqrt{4d\lvert
I_{w}-\Omega\rvert}\sigma_{\rho_{j}})^{-2}\int_{I_{w}-\Omega}\rho_{j}\quad.$
(44)

The gradients for the N–S unknowns $\boldsymbol{x}$ are drawn directly from
[9]

$\displaystyle{D}_{\boldsymbol{g}_{i}}\mathscr{J}$
$\displaystyle=\nu\partial_{\boldsymbol{\nu}}\boldsymbol{v}-q\boldsymbol{\nu}-\eta\boldsymbol{v}+{\mathcal{C}^{-1}_{\boldsymbol{g}_{i}}}\big{(}\boldsymbol{g}_{i}-{\bar{\boldsymbol{g}}_{i}}\big{)}\
,$ (45) $\displaystyle{D}_{\boldsymbol{g}_{o}}\mathscr{J}$
$\displaystyle=\boldsymbol{v}+{\mathcal{C}^{-1}_{\boldsymbol{g}_{o}}}\big{(}\boldsymbol{g}_{o}-\bar{\boldsymbol{g}}_{o}\big{)}\
,$ (46) $\displaystyle{D}_{\nu}\mathscr{J}$
$\displaystyle=\nabla\boldsymbol{v}\boldsymbol{:}\nabla\boldsymbol{u}+\sigma_{\nu}^{-2}\big{(}\nu-\bar{\nu}\big{)}\
,$ (47) $\displaystyle{D}_{{\phi_{\pm}}}\mathscr{(J-\mathscr{S})}$
$\displaystyle={\mathring{\zeta}}+\mathcal{C}^{-1}_{\phi_{\pm}}\big{(}\bar{\phi}_{\pm}-{\phi_{\pm}}\big{)}\
,$ (48)

where ${\mathring{\zeta}}\in L^{2}(M)$ is an extension of the shape gradient
$\zeta\in L^{2}(\Gamma)$ that is generated by an advection-diffusion problem
[9, Section 2.4]. The shape gradient is given by [9]

$\displaystyle\zeta=\partial_{\boldsymbol{\nu}}\boldsymbol{u}\boldsymbol{\cdot}\big{(}-\nu\partial_{\boldsymbol{\nu}}\boldsymbol{v}+q\boldsymbol{\nu}\big{)}\quad.$
(49)

The segmentation functional $\mathscr{S}$, given by (12), was not included in
[9] because the boundary $\partial\Omega$ was inferred only from the velocity
images. Here, $\mathscr{S}$ allows us to draw information from the spin
density (magnitude) images in order to better infer the boundary. For example,
near regions of low velocity gradients (e.g. nearly separated flow) the
velocity field alone is not informative. Consequently, taking into account
$\mathscr{S}$, we add the following term in (48)

$\displaystyle D_{{\phi_{\pm}}}\mathscr{S}={\mathring{\zeta}}_{\mathscr{S}}\in
L^{2}(M)\quad,$ (50)

where ${\mathring{\zeta}}_{\mathscr{S}}$ is an advection-diffusion extension
of $\zeta_{\mathscr{S}}\in L^{2}(\Gamma)$, similar to $\zeta^{\circ}$ [9], and

$\displaystyle\zeta_{\mathscr{S}}=\delta_{\Gamma}\sum_{j=1}^{4d}(\sqrt{8d}\sigma_{j})^{-2}\mathcal{S}^{\dagger}\Big{(}(\rho_{j}-\alpha)^{2}-(\rho_{j}-\beta)^{2}\Big{)},$
(51)

where $\delta_{\Gamma}$ is the Dirac measure on $\Gamma$, or, equivalently,
the characteristic function of $\Gamma$ ($\delta_{\Gamma}=\chi_{\Gamma}$).

##  Appendix B Phase-Contrast MRI experiment

The PC-MRI experiment for the flow through the converging nozzle is described
in detail in [9, Section 3.3]. In this paper, however, in order to demonstrate
the de-aliasing capability of algorithm 1, we use $k$-space signals (of the
same dataset) that produced severely aliased velocity images. These signals
correspond to slightly different flow conditions than those described in [9].
They were acquired for a 70 wt%, instead of a 40 wt%, glycerol in water
solution, and the Reynolds number was 22.4. We found the $T_{1}$ relaxation
time of this glycerol solution to be $313$ms. For the high SNR images, the
repetition time was $0.5$s, which resulted in a total acquisition time of 68
minutes per velocity image. For the low SNR images, the repetition time was
$150$ms, which resulted in a total acquisition time of 2.6 minutes per
velocity image. The total acquisition times correspond to fully-sampled
$k$-spaces.

##  Appendix C Prior information for the N–S unknowns

The prior variances of the Navier–Stokes unknowns are given in table III. Note
that, for 5% $\mathcal{P}_{\odot}$/$\mathcal{P}_{\parallel}$-sampling, instead
of $\sigma_{\phi_{\pm}}=2h_{z}$ we set ${\sigma_{\phi_{\pm}}=10h_{z}}$ as the
prior information on the geometry becomes less accurate due to subsampling
artefacts.

TABLE III: Confidence level of prior information. | shape $\partial\Omega$ | inlet $\boldsymbol{g}_{i}$ | viscosity $\nu$ | resolution
---|---|---|---|---
$\mathcal{P}$ | $(\sigma_{\phi_{\pm}},\text{Re}_{{\phi_{\pm}}},\text{Re}_{\zeta})$ |  $(\sigma_{\boldsymbol{g}_{i}}$ [cm/s], $\ell)$ |  $\sigma_{\nu}$ [m2/s] |  $h_{z}\times h_{r}$ [$\mu$m]
$\mathcal{P}_{\odot}$ | (2$h_{z}$, 0.025, 0.025) | (1.0, 3$h_{r}$) | 2.54$\times$106 | 165$\times$223
$\mathcal{P}_{\parallel}$ | (2$h_{z}$, 0.010, 0.010) | (1.0, 3$h_{r}$) | 2.54$\times$106 | 165$\times$223

## References

  * [1] H. Gudbjartsson and S. Patz, “The Rician Distribution of Noisy MRI Data,” _Magnetic Resonance in Medicine_ , vol. 34, no. 6, pp. 910–914, 1995.
  * [2] S. G. Chang, B. Yu, and M. Vetterli, “Adaptive wavelet thresholding for image denoising and compression,” _IEEE Transactions on Image Processing_ , vol. 9, no. 9, pp. 1532–1546, 2000.
  * [3] G. Pascal, “Rudin–Osher–Fatemi Total Variation Denoising using Split Bregman,” _Image Processing On Line_ , vol. 2, no. 1, pp. 74–95, 2012.
  * [4] N. Fatouraee and A. A. Amini, “Regularization of flow streamlines in multislice phase-contrast MR imaging,” _IEEE Transactions on Medical Imaging_ , vol. 22, no. 6, pp. 699–709, 2003.
  * [5] F. Ong, M. Uecker, U. Tariq, A. Hsiao, M. T. Alley, S. S. Vasanawala, and M. Lustig, “Robust 4D flow denoising using divergence-free wavelet transform,” _Magnetic Resonance in Medicine_ , vol. 73, no. 2, pp. 828–842, 2015.
  * [6] J. Mura, A. M. Pino, J. Sotelo, I. Valverde, C. Tejos, M. E. Andia, P. Irarrázaval, and S. Uribe, “Enhancing the Velocity Data From 4D Flow MR Images by Reducing its Divergence,” _IEEE Transactions on Medical Imaging_ , vol. 35, no. 10, pp. 2353–2364, 2016.
  * [7] T. S. Koltukluoğlu and P. J. Blanco, “Boundary control in computational haemodynamics,” _Journal of Fluid Mechanics_ , vol. 847, pp. 329–364, 2018\.
  * [8] S. W. Funke, M. Nordaas, Ø. Evju, M. S. Alnæs, and K. A. Mardal, “Variational data assimilation for transient blood flow simulations: Cerebral aneurysms as an illustrative example,” _International Journal for Numerical Methods in Biomedical Engineering_ , vol. 35, no. 1, pp. 1–27, 2019\.
  * [9] A. Kontogiannis, S. V. Elgersma, A. J. Sederman, and M. P. Juniper, “Joint reconstruction and segmentation of noisy velocity images as an inverse Navier–Stokes problem,” _Submitted to the Journal of Fluid Mechanics_ , 2021.
  * [10] E. J. Candès, J. Romberg, and T. Tao, “Robust Uncertainty Principles : Exact Signal Frequency Information,” _IEEE Transactions on Information Theory_ , vol. 52, no. 2, pp. 489–509, 2006.
  * [11] D. L. Donoho, “Compressed sensing,” _IEEE Transactions on Information Theory_ , vol. 52, no. 4, pp. 1289–1306, 2006.
  * [12] M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: The application of compressed sensing for rapid MR imaging,” _Magnetic Resonance in Medicine_ , vol. 58, no. 6, pp. 1182–1195, 2007.
  * [13] M. Benning, L. Gladden, D. Holland, C. B. Schönlieb, and T. Valkonen, “Phase reconstruction from velocity-encoded MRI measurements - A survey of sparsity-promoting variational approaches,” _Journal of Magnetic Resonance_ , vol. 238, pp. 26–43, 2014. [Online]. Available: http://dx.doi.org/10.1016/j.jmr.2013.10.003
  * [14] D. J. Holland, D. M. Malioutov, A. Blake, A. J. Sederman, and L. F. Gladden, “Reducing data acquisition times in phase-encoded velocity imaging using compressed sensing,” _Journal of Magnetic Resonance_ , vol. 203, no. 2, pp. 236–246, 2010. [Online]. Available: http://dx.doi.org/10.1016/j.jmr.2010.01.001
  * [15] T. Roberts, N. Kingsbury, and D. J. Holland, “Sparse recovery of complex phase-encoded velocity images using iterative thresholding,” _2013 IEEE International Conference on Image Processing, ICIP 2013 - Proceedings_ , pp. 350–354, 2013.
  * [16] F. Zhao, D. C. Noll, J. F. Nielsen, and J. A. Fessler, “Separate magnitude and phase regularization via compressed sensing,” _IEEE Transactions on Medical Imaging_ , vol. 31, no. 9, pp. 1713–1723, 2012.
  * [17] A. Reci, “Signal sampling and processing in magnetic resonance applications,” Ph.D. dissertation, University of Cambridge, 2019.
  * [18] K. Bredies, K. Kunisch, and T. Pock, “Total generalized variation,” _SIAM Journal on Imaging Sciences_ , vol. 3, no. 3, pp. 492–526, 2010.
  * [19] V. Corona, M. Benning, L. F. Gladden, A. Reci, A. J. Sederman, and C.-B. Schönlieb, “Joint Phase Reconstruction and Magnitude Segmentation from Velocity-Encoded MRI Data,” _Time-dependent Problems in Imaging and Parameter Identification_ , pp. 1–24, 2021. [Online]. Available: http://arxiv.org/abs/1908.05285
  * [20] J. Sotelo, J. Urbina, I. Valverde, C. Tejos, P. Irarrazaval, M. E. Andia, S. Uribe, and D. E. Hurtado, “3D Quantification of Wall Shear Stress and Oscillatory Shear Index Using a Finite-Element Method in 3D CINE PC-MRI Data of the Thoracic Aorta,” _IEEE Transactions on Medical Imaging_ , vol. 35, no. 6, pp. 1475–1487, 2016.
  * [21] J. Zhang, M. C. Brindise, S. Rothenberger, S. Schnell, M. Markl, D. Saloner, V. L. Rayz, and P. P. Vlachos, “4D Flow MRI Pressure Estimation Using Velocity Measurement-Error-Based Weighted Least-Squares,” _IEEE Transactions on Medical Imaging_ , vol. 39, no. 5, pp. 1668–1680, 2020.
  * [22] P. Berg, D. Stucht, G. Janiga, O. Beuing, O. Speck, and D. Thévenin, “Cerebral blood flow in a healthy circle of willis and two intracranial aneurysms: Computational fluid dynamics versus four-dimensional phase-contrast magnetic resonance imaging,” _Journal of Biomechanical Engineering_ , vol. 136, no. 4, pp. 1–9, 2014.
  * [23] J. Id and P. Berg, “2018 ( MATCH )— Phase Ib : Effect of morphology on hemodynamics,” vol. 2018, 2019.
  * [24] M. F. Duarte and Y. C. Eldar, “Structured compressed sensing: From theory to applications,” _IEEE Transactions on Signal Processing_ , vol. 59, no. 9, pp. 4053–4085, 2011.
  * [25] I. Bright, G. Lin, and J. N. Kutz, “Compressive sensing based machine learning strategy for characterizing the flow around a cylinder with limited pressure measurements,” _Physics of Fluids_ , vol. 25, no. 12, 2013.
  * [26] T. F. Chan and L. A. Vese, “Active contours without edges,” _IEEE Transactions on Image Processing_ , vol. 10, no. 2, pp. 266–277, 2001.
  * [27] P. Getreuer, “Chan–Vese Segmentation,” _Image Processing On Line_ , vol. 2, pp. 214–224, 2012.
  * [28] J. Nitsche, “Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind,” _Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg_ , vol. 36, no. 1, pp. 9–15, 1971.
  * [29] R. Fletcher, _Practical Methods of Optimization_. John Wiley & Sons, 2000.
  * [30] J. Nocedal and S. J. Wright, _Numerical optimization_ , 2006.
  * [31] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, “Array programming with NumPy,” _Nature_ , vol. 585, no. 7825, pp. 357–362, Sep. 2020. [Online]. Available: https://doi.org/10.1038/s41586-020-2649-2
  * [32] P. Virtanen _et al._ , “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” _Nature Methods_ , vol. 17, pp. 261–272, 2020.
  * [33] N.-S. Cheng, “Formula for the Viscosity of a Glycerol–Water Mixture,” _Industrial & Engineering Chemistry Research_, vol. 47, no. 9, pp. 3285–3288, may 2008. [Online]. Available: https://doi.org/10.1021/ie071349z
  * [34] A. Volk and C. J. Kähler, “Density model for aqueous glycerol solutions,” _Experiments in Fluids_ , vol. 59, no. 5, p. 75, 2018. [Online]. Available: https://doi.org/10.1007/s00348-018-2527-y
  * [35] L. F. Gladden and A. J. Sederman, “Magnetic Resonance Imaging and Velocity Mapping in Chemical Engineering Applications,” _Annual Review of Chemical and Biomolecular Engineering_ , vol. 8, no. 1, pp. 227–247, 2017.
  * [36] P. Bouillot, B. M. Delattre, O. Brina, R. Ouared, M. Farhat, C. Chnafa, D. A. Steinman, K. O. Lovblad, V. M. Pereira, and M. I. Vargas, “3D phase contrast MRI: Partial volume correction for robust blood flow quantification in small intracranial vessels,” _Magnetic Resonance in Medicine_ , vol. 79, no. 1, pp. 129–140, 2018.
  * [37] K. Saito, S. Abe, M. Kumamoto, Y. Uchihara, A. Tanaka, K. Sugie, M. Ihara, M. Koga, and H. Yamagami, “Blood Flow Visualization and Wall Shear Stress Measurement of Carotid Arteries Using Vascular Vector Flow Mapping,” _Ultrasound in Medicine and Biology_ , vol. 46, no. 10, pp. 2692–2699, 2020\.