# Joint reconstruction and segmentation of noisy velocity images as an inverse
Navier–Stokes problem

Alexandros Kontogiannis1<EMAIL_ADDRESS>Scott V. Elgersma2 Andrew J.
Sederman2  Matthew P. Juniper 1  1Department of Engineering, University of
Cambridge, Trumpington Street, Cambridge CB2 1PZ, UK 2Department of Chemical
Engineering & Biotechnology, University of Cambridge, Philippa Fawcett Drive,
Cambridge CB3 0AS, UK

###### Abstract

We formulate and solve a generalized inverse Navier–Stokes problem for the
joint velocity field reconstruction and boundary segmentation of noisy flow
velocity images. To regularize the problem we use a Bayesian framework with
Gaussian random fields. This allows us to estimate the uncertainties of the
unknowns by approximating their posterior covariance with a quasi-Newton
method. We first test the method for synthetic noisy images of 2D flows and
observe that the method successfully reconstructs and segments the noisy
synthetic images with a signal-to-noise ratio (SNR) of 3. Then we conduct a
magnetic resonance velocimetry (MRV) experiment to acquire images of an
axisymmetric flow for low ($\simeq 6$) and high ($>30$) SNRs. We show that the
method is capable of reconstructing and segmenting the low SNR images,
producing noiseless velocity fields and a smooth segmentation, with negligible
errors compared with the high SNR images. This amounts to a reduction of the
total scanning time by a factor of 27. At the same time, the method provides
additional knowledge about the physics of the flow (e.g. pressure), and
addresses the shortcomings of MRV (low spatial resolution and partial volume
effects) that otherwise hinder the accurate estimation of wall shear stresses.
Although the implementation of the method is restricted to 2D steady planar
and axisymmetric flows, the formulation applies immediately to 3D steady flows
and naturally extends to 3D periodic and unsteady flows.

###### keywords:

##  1 Introduction

Experimental measurements of fluid flows inside or around an object often
produce velocity images that contain noise. These images may be post-processed
in order to either reveal obscured flow patterns or to extract a quantity of
interest (e.g. pressure or wall shear stress). For example, magnetic resonance
velocimetry (MRV) (Fukushima, 1999; Mantle & Sederman, 2003; Elkins & Alley,
2007; Markl et al., 2012; Demirkiran et al., 2021) can measure all three
components of a time varying velocity field but the measurements become
increasingly noisy as the spatial resolution is increased. To achieve an image
of acceptable signal-to-noise ratio (SNR), repeated scans are often averaged,
leading to long signal acquisition times. To address that problem, fast
acquisition protocols (pulse sequences) can be used, but these may be
difficult to implement and can lead to artefacts depending on the magnetic
relaxation properties and the magnetic field homogeneity of the system
studied. Another way to accelerate signal acquisition is by using sparse
sampling techniques in conjunction with a reconstruction algorithm. The latter
approach is an active field of research, commonly referred to as compressed
sensing (Donoho, 2006; Lustig et al., 2007; Benning et al., 2014; Peper et
al., 2019; Corona et al., 2021). Compressed sensing (CS) algorithms exploit a
priori knowledge about the structure of the data, which is encoded in a
regularization norm (e.g. total variation, wavelet bases), but without
considering the physics of the problem. Even though the present study concerns
the reconstruction of fully-sampled, noisy MRV images, the method that we
present here can be applied to sparsely-sampled MRV data.

For images depicting fluid flow, a priori knowledge can come in the form of a
Navier–Stokes problem. The problem of reconstructing and segmenting a flow
image then can be expressed as a generalized inverse Navier–Stokes problem
whose flow domain, boundary conditions, and model parameters have to be
inferred in order for the modeled velocity to approximate the measured
velocity in an appropriate metric space. This approach not only produces a
reconstruction that is an accurate fluid flow inside or around the object (a
solution to a Navier–Stokes problem), but also provides additional physical
knowledge (e.g. pressure), which is otherwise difficult to measure. Inverse
Navier–Stokes problems have been intensively studied during the last decade,
mainly enabled by the increase of available computing power. Recent
applications in fluid mechanics range from the forcing inference problem
(Hoang et al., 2014), to the reconstruction of scalar image velocimetry (SIV)
(Gillissen et al., 2018; Sharma et al., 2019) and particle image velocimetry
(PIV) (Gillissen et al., 2019) signals, and the identification of optimal
sensor arrangements (Mons et al., 2017; Verma et al., 2020). Regularization
methods that can be used for model parameters are reviewed by Stuart (2010)
from a Bayesian perspective and by Benning & Burger (2018) from a variational
perspective. The well-posedness of Bayesian inverse Navier–Stokes problems is
addressed by Cotter et al. (2009).

Recently, Koltukluoğlu & Blanco (2018) treat the reduced inverse Navier–Stokes
problem of finding only the Dirichlet boundary condition for the inlet
velocity that matches the modeled velocity field to MRV data for a steady 3D
flow in a glass replica of the human aorta. They measure the model-data
discrepancy using the $L^{2}$-norm and introduce additional variational
regularization terms for the Dirichlet boundary condition. The same
formulation is extended to periodic flows by Koltukluoğlu (2019); Koltukluoğlu
et al. (2019), using the harmonic balance method for the temporal
discretization of the Navier–Stokes problem. Funke et al. (2019) address the
problem of inferring both the inlet velocity (Dirichlet) boundary condition
and the initial condition, for unsteady blood flows and 4D MRV data, with
applications to cerebral aneurysms. We note that the above studies consider
rigid boundaries and require a priori an accurate, and time-averaged,
geometric representation of the blood vessel.

To find the shape of the flow domain, e.g. the blood vessel boundaries,
computed tomography (CT) or magnetic resonance angiography (MRA) is often
used. The acquired image is then reconstructed, segmented, and smoothed. This
process not only requires substantial effort and the design of an additional
experiment (e.g. CT, MRA), but it also introduces geometric uncertainties
(Morris et al., 2016; Sankaran et al., 2016), which, in turn, affect the
predictive confidence of arterial wall shear stress distributions and their
mappings (Katritsis et al., 2007; Sotelo et al., 2016). For example, Funke et
al. (2019) report discrepancies between the modeled and the measured velocity
fields near the flow boundaries, and they suspect they are caused by geometric
errors that were introduced during the segmentation process. In general, the
assumption of rigid boundaries either implies that a time-averaged geometry
has to be used, or that an additional experiment (e.g. CT, MRA) has to be
conducted to register the moving boundaries to the flow measurements.

A more consistent approach to this problem is to treat the blood vessel
geometry as an unknown when solving the generalized inverse Navier–Stokes
problem. In this way, the inverse Navier–Stokes problem simultaneously
reconstructs and segments the velocity fields and can better adapt to the MRV
experiment by correcting the geometric errors and improving the
reconstruction.

In this study, we address the problem of simultaneous velocity field
reconstruction and boundary segmentation by formulating a generalized inverse
Navier–Stokes problem, whose flow domain, boundary conditions, and model
parameters are all considered unknown. To regularize the problem, we use a
Bayesian framework and Gaussian measures in Hilbert spaces. This further
allows us to estimate the posterior Gaussian distributions of the unknowns
using a quasi-Newton method, which has not yet been addressed for this type of
problem. We provide an algorithm for the solution of this generalized inverse
Navier–Stokes problem, and demonstrate it on synthetic images of 2D steady
flows and real MRV images of a steady axisymmetric flow.

This paper consists of two parts. In section 2, we formulate the generalized
inverse Navier–Stokes problem and an algorithm that solves it. In section 3 we
test the method using both synthetic and real MRV velocity images and describe
the setup of the MRV experiment.

##  2 An inverse Navier–Stokes problem for noisy flow images

In this section, we formulate the generalized inverse Navier–Stokes problem
and provide an algorithm for its solution. 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)$ the space
of square-integrable functions with $k$ square-integrable derivatives in
$\Omega$. For a given covariance operator, $\mathcal{C}$, we also define the
covariance-weighted $L^{2}$ spaces, endowed with the inner product
${\big{\langle}\cdot,\cdot\big{\rangle}_{\mathcal{C}}:=\big{\langle}\cdot,\mathcal{C}^{-1}\cdot\big{\rangle}}$,
which generates the norm $\big{\lVert}\cdot\big{\rVert}_{\mathcal{C}}$. The
Euclidean norm in the space of real numbers $\mathbb{R}^{n}$ is denoted by
$\lvert\cdot\rvert_{\mathbb{R}^{n}}$. We use the superscript $(\cdot)^{\star}$
to denote a measurement, $(\cdot)^{\circ}$ to denote a reconstruction, and
$(\cdot)^{\bullet}$ to denote the ground truth.

###  2.1 The inverse Navier–Stokes problem

A $n$-dimensional velocimetry experiment usually provides noisy flow velocity
images on a domain $I\subset\mathbb{R}^{n}$, depicting the measured flow
velocity $\boldsymbol{u}^{\star}$ inside an object $\Omega\subset I$ with
boundary $\partial\Omega=\Gamma\cup\Gamma_{i}\cup\Gamma_{o}$ (figure 1). An
appropriate model is the Navier–Stokes problem

$\left\\{\begin{aligned}
\boldsymbol{u}\boldsymbol{\cdot}\nabla\boldsymbol{u}-\nu{\rmDelta}\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,$
(1)

where $\boldsymbol{u}$ is the velocity,
$p\mathrel{\reflectbox{$\mapsto$}}p/\rho$ is the reduced 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 denote the data space by $\boldsymbol{D}$ and the model space by
$\boldsymbol{M}$, and assume that both spaces are subspaces of
$\boldsymbol{L}^{2}$. In the 2D case,
$\boldsymbol{u}^{\star}=(u^{\star}_{x},u^{\star}_{y})$, and we introduce the
covariance operator

${{\mathcal{C}_{\boldsymbol{u}}}=\mathrm{diag}\Big{(}\sigma^{2}_{u_{x}}\mathrm{I},~{}\sigma^{2}_{u_{y}}\mathrm{I}\Big{)}}\quad,$
(2)

where $\sigma^{2}_{u_{x}},\sigma^{2}_{u_{y}}$ are the Gaussian noise variances
of $u^{\star}_{x},u^{\star}_{y}$, respectively, and $\mathrm{I}$ is the
identity operator. The discrepancy between the measured velocity field
${\boldsymbol{u}^{\star}}\in\boldsymbol{D}$ and the modeled velocity field
$\boldsymbol{u}\in\boldsymbol{M}$ is measured on the data space
$\boldsymbol{D}$ using the reconstruction error functional

$\mathscr{E}(\boldsymbol{u})\equiv\frac{1}{2}\big{\lVert}\boldsymbol{u}^{\star}-\mathcal{S}\boldsymbol{u}\big{\rVert}^{2}_{\mathcal{C}_{\boldsymbol{u}}}:=\frac{1}{2}\int_{I}\big{(}\boldsymbol{u}^{\star}-\mathcal{S}\boldsymbol{u}\big{)}\mathcal{C}_{\boldsymbol{u}}^{-1}\big{(}\boldsymbol{u}^{\star}-\mathcal{S}\boldsymbol{u}\big{)}\quad,$
(3)

where $\mathcal{S}:\boldsymbol{M}\to\boldsymbol{D}$ is the
$L^{2}$-projection111Since the discretized space consists of bilinear
quadrilateral finite elements (see section 2.7), this projection is a linear
interpolation. from the model space $\boldsymbol{M}$ to the data space
$\boldsymbol{D}$.

Our goal is to infer the unknown parameters of the Navier–Stokes problem (1)
such that the model velocity $\boldsymbol{u}$ approximates the noisy measured
velocity $\boldsymbol{u}^{\star}$ in the covariance-weighted $L^{2}$-metric
defined by $\mathscr{E}$. In the general case, the unknown model parameters of
(1) are the shape of $\Omega$, the kinematic viscosity $\nu$, and the boundary
conditions $\boldsymbol{g}_{i},\boldsymbol{g}_{o}$. This inverse Navier–Stokes
problem leads to the nonlinearly constrained optimization problem

$\text{find}\quad\boldsymbol{u}^{\circ}\equiv\underset{\Omega,\boldsymbol{x}}{\mathrm{argmin}}~{}\mathscr{E}\big{(}\boldsymbol{u}(\Omega;\boldsymbol{x})\big{)},\quad\text{such
that $\boldsymbol{u}$ satisfies}\ \eqref{eq:navierstokes_bvp}\quad,$ (4)

where $\boldsymbol{u}^{\circ}$ is the reconstructed velocity field, and
$\boldsymbol{x}=(\boldsymbol{g}_{i},\boldsymbol{g}_{o},\nu)$. Like most
inverse problems, (4) is ill-posed and hard to solve. To alleviate the ill-
posedness of the problem we need to restrict our search of the unknowns
$(\Omega,\boldsymbol{x})$ to function spaces of sufficient regularity.

Figure 1: Given the images of a measured velocity field
$\boldsymbol{u}^{\star}$, we solve an inverse Navier–Stokes problem to infer
the boundary $\Gamma$ (or $\partial\Omega$), the kinematic viscosity, and the
inlet velocity profile on $\Gamma_{i}$. The solution to this inverse problem
is a reconstructed velocity field $\boldsymbol{u}^{\circ}$, from which the
noise and the artefacts
$(\boldsymbol{u}^{\star}-\mathcal{S}\boldsymbol{u}^{\circ})$ have been
filtered out.

###  2.2 Regularization

If $x(t)\in L^{2}(\mathbb{R})$ is an unknown parameter, one way to regularize
the inverse problem (4) is to search for minimizers of the augmented
functional $\mathscr{J}\equiv\mathscr{E}+\mathscr{R}$, where

$\mathscr{R}(x)=\sum_{j=0}^{k}\int_{\mathbb{R}}\alpha_{j}\big{|}\partial_{x}^{j}(x-\bar{x})\big{|}^{2}$
(5)

is a regularization norm for a given (and fixed) prior assumption
$\bar{x}(t)\in H^{k}(\mathbb{R})$, weights ${\alpha_{j}\in\mathbb{R}}$, and
positive integer $k$. This simple idea can be quite effective because by
minimizing $\mathscr{R}$ we force $x$ to lie in a subspace of $L^{2}$ having
higher regularity, namely $H^{k}$, and as close to the prior value $\bar{x}$
as $\alpha_{j}$ allow222The regularization term, given by (5), can be further
extended to fractional Hilbert spaces by defining the norm
$\big{\lVert}x\big{\rVert}_{H^{s}(\mathbb{R})}:=\big{\lVert}(1+\lvert
t\rvert^{s})\mathcal{F}x\big{\rVert}_{L^{2}(\mathbb{R})}$ for noninteger $s$,
with $0<s<\infty$, and where $\mathcal{F}$ denotes the Fourier transform.
Interestingly, under certain conditions, which are dictated by Sobolev’s
embedding theorem (Evans, 2010, Chapter 5), these Hilbert spaces can be
embedded in the more familiar spaces of continuous functions.. However, as
Stuart (2010) points out, in this setting, the choice of $\alpha_{j}$, and
even the form of $\mathscr{R}$, is arbitrary.

There is a more intuitive approach that recovers the form of the
regularization norm $\mathscr{R}$ from a probabilistic viewpoint. In the
setting of the Hilbert space $L^{2}$, the Gaussian measure
$\gamma\sim\mathcal{N}(m,\mathcal{C})$ has the property that its finite-
dimensional projections are multivariate Gaussian distributions, and it is
uniquely defined by its mean $m\in L^{2}$, and its covariance operator
$\mathcal{C}:L^{2}\to L^{2}$ (appendix A). It can be shown that there is a
natural Hilbert space $H_{\gamma}$ that corresponds to $\gamma$, and that
(Bogachev, 1998; Hairer, 2009)

$H_{\gamma}=\sqrt{\mathcal{C}}\big{(}L^{2}\big{)}\quad.$

In other words, if $x$ is a random function distributed according to $\gamma$,
any realization of $x$ lies in $H_{\gamma}$, which is the image of
$\sqrt{\mathcal{C}}$. Furthermore, the corresponding inner product

$\displaystyle\big{\langle}x,x^{\prime}\big{\rangle}_{\mathcal{C}}=\big{\langle}\mathcal{C}^{-1/2}x,\
\mathcal{C}^{-1/2}x^{\prime}\big{\rangle}$ (6)

is the covariance between $x$ and $x^{\prime}$, and the norm
$\big{\lVert}x\big{\rVert}^{2}_{\mathcal{C}}=\big{\langle}x,x\big{\rangle}_{\mathcal{C}}$
is the variance of $x$. Therefore, if $x$ is an unknown parameter for which a
priori statistical information is available, and if the Gaussian assumption
can be justified, we can choose

$\mathscr{R}(x)=\frac{1}{2}\big{\lVert}x-\bar{x}\big{\rVert}^{2}_{\mathcal{C}}\quad.$
(7)

In this way, $\mathscr{J}\equiv\mathscr{E}+\mathscr{R}$ increases as the
variance of $x$ increases. Consequently, minimizing $\mathscr{J}$ penalizes
improbable realizations.

As mentioned in section 2.1, the unknown model parameters of the Navier–Stokes
problem (1) are the kinematic viscosity $\nu$, the boundary conditions
$\boldsymbol{g}_{i},\boldsymbol{g}_{o}$, and the shape of $\Omega$. Since we
consider the kinematic viscosity $\nu$ to be constant, the regularizing norm
is simply

$\displaystyle\frac{1}{2}\big{|}\nu-\bar{\nu}\big{|}^{2}_{\Sigma_{\nu}}=\frac{1}{2\sigma^{2}_{\nu}}\big{|}\nu-\bar{\nu}\big{|}^{2}_{\mathbb{R}}\quad,$
(8)

where $\bar{\nu}\in\mathbb{R}$ is a prior guess for $\nu$, and
$\sigma^{2}_{\nu}\in\mathbb{R}$ is the variance. For the Dirichlet boundary
condition, $\boldsymbol{g}_{i}\in\boldsymbol{L}^{2}(\Gamma_{i})$, we choose
the exponential covariance function

$\displaystyle
C(x,x^{\prime})=\frac{\sigma_{\boldsymbol{g}_{i}}^{2}}{2\ell}\exp\bigg{(}-\frac{\lvert
x-x^{\prime}\rvert}{\ell}\bigg{)}$ (9)

with variance $\sigma_{\boldsymbol{g}_{i}}^{2}\in\mathbb{R}$ and
characteristic length $\ell\in\mathbb{R}$. For zero-Dirichlet (no-slip) or
zero-Neumann boundary conditions on $\partial\Gamma_{i}$, (9) leads to the
norm (Tarantola, 2005, Chapter 7.21)

$\displaystyle\big{\lVert}\boldsymbol{g}_{i}\big{\rVert}_{\mathcal{C}_{\boldsymbol{g}_{i}}}^{2}\simeq\frac{1}{\sigma_{\boldsymbol{g}_{i}}^{2}}\int_{\Gamma_{i}}\boldsymbol{g}^{2}_{i}+\ell^{2}~{}\big{(}\nabla\boldsymbol{g}_{i}\big{)}^{2}\quad.$
(10)

Using integration by parts we find that the covariance operator is

$\displaystyle\mathcal{C}_{\boldsymbol{g}_{i}}=\sigma_{\boldsymbol{g}_{i}}^{2}\Big{(}\mathrm{I}-\ell^{2}\widetilde{\rmDelta}\Big{)}^{-1}\quad,$
(11)

where $\widetilde{\rmDelta}$ is the $L^{2}$-extension of the Laplacian
$\rmDelta$ that incorporates the boundary condition
$\boldsymbol{g}_{i}=\boldsymbol{0}$ on ${\partial\Gamma_{i}}$. For the natural
boundary condition, $\boldsymbol{g}_{o}\in\boldsymbol{L}^{2}(\Gamma_{o})$, we
can use the same covariance operator, but equip $\widetilde{\rmDelta}$ with
zero-Neumann boundary conditions, i.e.
$\partial_{\boldsymbol{\nu}}\boldsymbol{g}_{o}=0$ on $\partial\Gamma_{o}$.
Lastly, for the shape of $\Omega$, which we implicitly represent with a signed
distance function ${\phi_{\pm}}$ (defined in section 2.4), we choose the norm

$\displaystyle\frac{1}{2}\big{\lVert}\bar{\phi}_{\pm}-{\phi_{\pm}}\big{\rVert}^{2}_{\mathcal{C}_{\phi_{\pm}}}=\frac{1}{2\sigma^{2}_{\phi_{\pm}}}\big{\lVert}\bar{\phi}_{\pm}-{\phi_{\pm}}\big{\rVert}^{2}_{L^{2}(I)},$
(12)

where $\sigma_{\phi_{\pm}}\in\mathbb{R}$ and $\bar{\phi}_{\pm}\in L^{2}(I)$.
Additional regularization for the boundary of $\Omega$ (i.e. the zero level-
set of ${\phi_{\pm}}$) is needed and it is described in section 2.4. Based on
the above results, the regularization norm for the unknown model parameters is

$\displaystyle\mathscr{R}(\boldsymbol{x},{\phi_{\pm}})=~{}$
$\displaystyle\frac{1}{2}\big{|}\nu-\bar{\nu}\big{|}^{2}_{\Sigma_{\nu}}+\frac{1}{2}\big{\lVert}\boldsymbol{g}_{i}-\bar{\boldsymbol{g}}_{i}\big{\rVert}^{2}_{\mathcal{C}_{\boldsymbol{g}_{i}}}$
$\displaystyle+$
$\displaystyle\frac{1}{2}\big{\lVert}\boldsymbol{g}_{o}-\bar{\boldsymbol{g}}_{o}\big{\rVert}^{2}_{\mathcal{C}_{\boldsymbol{g}_{o}}}+\frac{1}{2}\big{\lVert}\bar{\phi}_{\pm}-{\phi_{\pm}}\big{\rVert}^{2}_{\mathcal{C}_{\phi_{\pm}}}\quad.$
(13)

###  2.3 Euler–Lagrange equations for the inverse Navier–Stokes problem

Testing the Navier–Stokes problem (1) with functions
${(\boldsymbol{v},q)\in\boldsymbol{H}^{1}(\Omega)\times L^{2}(\Omega)}$, and
after integrating by parts, we obtain the weak form

$\displaystyle\mathscr{M}(\Omega)(\boldsymbol{u},p,\boldsymbol{v},q;\boldsymbol{x})\equiv\int_{\Omega}\Big{(}\boldsymbol{v}\boldsymbol{\cdot}\big{(}\boldsymbol{u}\boldsymbol{\cdot}\nabla\boldsymbol{u}\big{)}$
$\displaystyle+\nu\nabla\boldsymbol{v}\boldsymbol{:}\nabla\boldsymbol{u}-(\nabla\boldsymbol{\cdot}\boldsymbol{v})p-q(\nabla\boldsymbol{\cdot}\boldsymbol{u})\Big{)}+\int_{\Gamma_{o}}\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{g}_{o}$
$\displaystyle+\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})={0}\quad,$
(14)

where $\mathscr{N}$ is the Nitsche (1971) penalty term

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

which weakly imposes the Dirichlet boundary condition
$\boldsymbol{z}\in\boldsymbol{L}^{2}(T)$ on a boundary $T$, given a
penalization constant $\eta$333The penalization $\eta$ is a numerical
parameter with no physical significance (see section 2.7).. We define the
augmented reconstruction error functional

$\displaystyle\mathscr{J}(\Omega)(\boldsymbol{u},p,\boldsymbol{v},q;\boldsymbol{x})\equiv\
$
$\displaystyle\mathscr{E}(\boldsymbol{u})+\mathscr{R}(\boldsymbol{x},{\phi_{\pm}})+\mathscr{M}(\Omega)(\boldsymbol{u},p,\boldsymbol{v},q;\boldsymbol{x})\quad,$
(16)

which contains the regularization terms $\mathscr{R}$ and the model constraint
$\mathscr{M}$, such that $\boldsymbol{u}$ weakly satisfies (1). To reconstruct
the measured velocity field $\boldsymbol{u}^{\star}$ and find the unknowns
$(\Omega,\boldsymbol{x})$, we minimize $\mathscr{J}$ by solving its associated
Euler–Lagrange system.

####  2.3.1 Adjoint Navier–Stokes problem

In order to derive the Euler–Lagrange equations for $\mathscr{J}$, we first
define

$\displaystyle\boldsymbol{\mathcal{U}}^{\prime}=\Big{\\{}\boldsymbol{u}^{\prime}\in\boldsymbol{H}^{1}(\Omega):\boldsymbol{u}^{\prime}\big{|}_{\Gamma\cup\Gamma_{i}}\equiv\boldsymbol{0}\Big{\\}}$
(17)

to be the space of admissible velocity perturbations
$\boldsymbol{u}^{\prime}$, and $\mathcal{P}^{\prime}\subset L^{2}(\Omega)$ to
be the space of admissible pressure perturbations $p^{\prime}$, such that
$(-\partial_{\boldsymbol{\nu}}\boldsymbol{u}^{\prime}+p^{\prime}\boldsymbol{\nu})\big{|}_{\Gamma_{o}}\equiv\boldsymbol{0}$.
We start with

$\displaystyle\delta_{\boldsymbol{u}}\mathscr{E}\equiv\frac{d}{d\tau}\mathscr{E}(\boldsymbol{u}+\tau\boldsymbol{u}^{\prime})\Big{|}_{\tau=0}$
$\displaystyle=\int_{\Omega}-{\mathcal{C}^{-1}_{\boldsymbol{u}}}\big{(}\boldsymbol{u}^{\star}-\mathcal{S}\boldsymbol{u}\big{)}\boldsymbol{\cdot}\mathcal{S}\boldsymbol{u}^{\prime}$
$\displaystyle=\int_{\Omega}-\mathcal{S}^{\dagger}{\mathcal{C}^{-1}_{\boldsymbol{u}}}\big{(}\boldsymbol{u}^{\star}-\mathcal{S}\boldsymbol{u}\big{)}\boldsymbol{\cdot}\boldsymbol{u}^{\prime}\equiv\Big{\langle}D_{\boldsymbol{u}}\mathscr{E},\boldsymbol{u}^{\prime}\Big{\rangle}_{\Omega}\quad.$
(18)

Adding together the first variations of $\mathscr{M}$ with respect to
$(\boldsymbol{u},p)$,

$\displaystyle\delta_{\boldsymbol{u}}\mathscr{M}\equiv\frac{d}{d\tau}\mathscr{M}(\cdot)(\boldsymbol{u}+\tau\boldsymbol{u}^{\prime},\dots)\Big{|}_{\tau=0}\quad,\quad\delta_{{p}}\mathscr{M}\equiv\frac{d}{d\tau}\mathscr{M}(\cdot)(\dots,p+\tau
p^{\prime},\dots)\Big{|}_{\tau=0},$

and after integrating by parts, we find

$\displaystyle\delta_{\boldsymbol{u}}\mathscr{M}+\delta_{p}\mathscr{M}=$
$\displaystyle\int_{\Omega}\Big{(}-\boldsymbol{u}\boldsymbol{\cdot}\big{(}\nabla\boldsymbol{v}+(\nabla\boldsymbol{v})^{\dagger}\big{)}-\nu\rmDelta\boldsymbol{v}+\nabla
q\Big{)}\boldsymbol{\cdot}\boldsymbol{u}^{\prime}+\int_{\Omega}(\nabla\boldsymbol{\cdot}\boldsymbol{v})p^{\prime}$
$\displaystyle+$
$\displaystyle\int_{\partial\Omega}\big{(}{(\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}\big{)}\boldsymbol{\cdot}\boldsymbol{u}^{\prime}$
$\displaystyle+$
$\displaystyle\int_{\Gamma\cup\Gamma_{i}}\boldsymbol{v}\boldsymbol{\cdot}(-\nu\partial_{\boldsymbol{\nu}}\boldsymbol{u}^{\prime}+p^{\prime}\boldsymbol{\nu})+\mathscr{N}_{\Gamma_{i}\cup\Gamma}(\boldsymbol{v},q,\boldsymbol{u}^{\prime};\boldsymbol{0})\quad.$
(19)

Since $\mathscr{R}$ does not depend on $(\boldsymbol{u},p)$, we can use (18)
and (19) to assemble the optimality conditions of $\mathscr{J}$ for
$(\boldsymbol{u},p)$

$\displaystyle\Big{\langle}D_{\boldsymbol{u}}\mathscr{J},\boldsymbol{u}^{\prime}\Big{\rangle}_{\Omega}={0}\quad,\quad\Big{\langle}D_{{p}}\mathscr{J},{p}^{\prime}\Big{\rangle}_{\Omega}=0\quad.$
(20)

For (20) to hold true for all perturbations
${(\boldsymbol{u}^{\prime},p^{\prime})\in\boldsymbol{\mathcal{U}}^{\prime}\times\mathcal{P}^{\prime}}$,
we deduce that $(\boldsymbol{v},q)$ must satisfy the following adjoint
Navier–Stokes problem

$\displaystyle\left\\{\begin{aligned}
-\boldsymbol{u}\boldsymbol{\cdot}\big{(}\nabla\boldsymbol{v}+(\nabla\boldsymbol{v})^{\dagger}\big{)}-\nu\rmDelta\boldsymbol{v}+\nabla
q&=-D_{\boldsymbol{u}}\mathscr{E}\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.\quad.$
(21)

In this context, $\boldsymbol{v}$ is the adjoint velocity and $q$ is the
adjoint pressure, which both vanish when
$\boldsymbol{u}^{\star}\equiv\mathcal{S}\boldsymbol{u}$. Note also that we
choose boundary conditions for the adjoint problem (21) that make the boundary
terms of (19) vanish, and that these boundary conditions are subject to the
choice of $\boldsymbol{\mathcal{U}}^{\prime}$, which, in turn, depends on the
boundary conditions of the (primal) Navier–Stokes problem.

####  2.3.2 Shape derivatives for the Navier–Stokes problem

To find the shape derivative of an integral defined in $\Omega$, when the
boundary $\partial\Omega$ deforms with speed $\mathscr{V}$, we use Reynold’s
transport theorem. For the bulk integral of $f:\Omega\to\mathbb{R}$, we find

$\displaystyle\frac{d}{d\tau}\bigg{(}\int_{\Omega(\tau)}f\bigg{)}\bigg{|}_{\tau=0}=\int_{\Omega}f^{\prime}+\int_{\partial\Omega}f~{}(\mathscr{V}\boldsymbol{\cdot}\boldsymbol{\nu})\quad,$
(22)

while for the boundary integral of $f$ we find (Walker, 2015, Chapter 5.6)

$\displaystyle\frac{d}{d\tau}\bigg{(}\int_{\partial\Omega(\tau)}f\bigg{)}\bigg{|}_{\tau=0}=\int_{\partial\Omega}f^{\prime}+(\partial_{\boldsymbol{\nu}}+\kappa)f~{}(\mathscr{V}\boldsymbol{\cdot}\boldsymbol{\nu})\quad,$
(23)

where $f^{\prime}$ is the shape derivative of $f$ (due to $\mathscr{V}$),
$\kappa$ is the summed curvature of $\partial\Omega$, and
$\mathscr{V}\equiv\zeta\boldsymbol{\nu}$, with $\zeta\in
L^{2}(\partial\Omega)$, is the Hadamard parameterization of the speed field.
Any boundary that is a subset of $\partial I$, i.e. the edge of the image $I$,
is non-deforming and therefore the second term of the above integrals
vanishes. The only boundary that deforms is $\Gamma\subset\partial\Omega$. For
brevity, let $\delta_{\mathscr{V}}I$ denote the shape perturbation of an
integral $I$. Using (22) on $\mathscr{E}$, we compute

$\displaystyle\delta_{\mathscr{V}}\mathscr{E}=\Big{\langle}D_{\boldsymbol{u}}\mathscr{E},\boldsymbol{u}^{\prime}\Big{\rangle}_{\Omega}\quad,$
(24)

where $D_{\boldsymbol{u}}\mathscr{E}$ is given by (18). Using (22) and (23) on
$\mathscr{M}$, we obtain the shape derivatives problem for
$(\boldsymbol{u}^{\prime},p^{\prime})$

$\left\\{\begin{aligned}
\boldsymbol{u}^{\prime}\boldsymbol{\cdot}\nabla\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot}\nabla\boldsymbol{u}^{\prime}-\nu{\rmDelta}\boldsymbol{u}^{\prime}+\nabla
p^{\prime}&=\boldsymbol{0}\quad&&\textrm{in}\quad\Omega\\\
\nabla\boldsymbol{\cdot}\boldsymbol{u}^{\prime}&=0\quad&&\textrm{in}\quad\Omega\\\
\boldsymbol{u}^{\prime}&=-\partial_{\boldsymbol{\nu}}\boldsymbol{u}(\mathscr{V}\boldsymbol{\cdot}\boldsymbol{\nu})\quad&&\textrm{on}\quad\Gamma\\\
\boldsymbol{u}^{\prime}&=\boldsymbol{0}\quad&&\textrm{on}\quad\Gamma_{i}\\\
-\nu\partial_{\boldsymbol{\nu}}\boldsymbol{u}^{\prime}+p^{\prime}\boldsymbol{\nu}&=\boldsymbol{0}\quad&&\textrm{on}\quad\Gamma_{o}\end{aligned}\right.\quad,$
(25)

which can be used directly to compute the velocity and pressure perturbations
for a given speed field $\mathscr{V}$. We observe that
$(\boldsymbol{u}^{\prime},p^{\prime})\equiv\boldsymbol{0}$ when
$\zeta\equiv\mathscr{V}\boldsymbol{\cdot}\boldsymbol{\nu}\equiv{0}$. Testing
the shape derivatives problem (25) with $(\boldsymbol{v},q)$, and adding the
appropriate Nitsche terms for the weakly enforced Dirichlet boundary
conditions, we obtain

$\displaystyle\delta_{\mathscr{V}}\mathscr{M}=$
$\displaystyle\int_{\Omega}\Big{(}\boldsymbol{v}\boldsymbol{\cdot}\big{(}\boldsymbol{u}^{\prime}\boldsymbol{\cdot}\nabla\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot}\nabla\boldsymbol{u}^{\prime}\big{)}+\nu\nabla\boldsymbol{v}\boldsymbol{:}\nabla\boldsymbol{u}^{\prime}-(\nabla\boldsymbol{\cdot}\boldsymbol{v})p^{\prime}-q(\nabla\boldsymbol{\cdot}\boldsymbol{u}^{\prime})\Big{)}$
$\displaystyle+\int_{\Gamma\cup\Gamma_{i}}$
$\displaystyle\boldsymbol{v}\boldsymbol{\cdot}(-\nu\partial_{\boldsymbol{\nu}}\boldsymbol{u}^{\prime}+p^{\prime}\boldsymbol{\nu})+\mathscr{N}_{\Gamma_{i}}(\boldsymbol{v},q,\boldsymbol{u}^{\prime};\boldsymbol{0})+\mathscr{N}_{\Gamma}(\boldsymbol{v},q,\boldsymbol{u}^{\prime};-\zeta\partial_{\boldsymbol{\nu}}\boldsymbol{u})={0}\quad.$
(26)

If we define $I_{i}$ to be the first four integrals in (19), integrating (26)
by parts yields

$\displaystyle\delta_{\mathscr{V}}\mathscr{M}=$
$\displaystyle\sum_{i=1}^{4}I_{i}+\mathscr{N}_{\Gamma_{i}}(\boldsymbol{v},q,\boldsymbol{u}^{\prime};\boldsymbol{0})+\mathscr{N}_{\Gamma}(\boldsymbol{v},q,\boldsymbol{u}^{\prime};-\zeta\partial_{\boldsymbol{\nu}}\boldsymbol{u})={0}\quad,$
(27)

and, due to the adjoint problem (21), we find

$\displaystyle\delta_{\mathscr{V}}\mathscr{M}=$
$\displaystyle-\delta_{\mathscr{V}}\mathscr{E}+\int_{\Gamma}(-\nu\partial_{\boldsymbol{\nu}}\boldsymbol{v}+q\boldsymbol{\nu})\boldsymbol{\cdot}\zeta\partial_{\boldsymbol{\nu}}\boldsymbol{u}={0}\quad,$
(28)

since $\boldsymbol{u}^{\prime}\big{|}_{\Gamma_{i}}\equiv\boldsymbol{0}$, and
$\boldsymbol{u}\big{|}_{\Gamma}\equiv\boldsymbol{0}$. Therefore, the shape
perturbation of $\mathscr{J}$ is

$\displaystyle\delta_{\mathscr{V}}\mathscr{J}\equiv\Big{\langle}D_{\mathscr{V}}\mathscr{J},\mathscr{V}\boldsymbol{\cdot}\boldsymbol{\nu}\Big{\rangle}_{\Gamma}\equiv\Big{\langle}D_{\zeta}\mathscr{J},\zeta\Big{\rangle}_{\Gamma}=\delta_{\mathscr{V}}\mathscr{E}+\delta_{\mathscr{V}}\mathscr{M}+\delta_{\mathscr{V}}\mathscr{R}=0\quad,$
(29)

which, due to (28) and $\delta_{\mathscr{V}}\mathscr{R}\equiv 0$, takes the
form

$\Big{\langle}D_{\zeta}\mathscr{J},\zeta\Big{\rangle}_{\Gamma}=\Big{\langle}\partial_{\boldsymbol{\nu}}\boldsymbol{u}\boldsymbol{\cdot}\big{(}-\nu\partial_{\boldsymbol{\nu}}\boldsymbol{v}+q\boldsymbol{\nu}\big{)},\zeta\Big{\rangle}_{\Gamma}\quad,$
(30)

where $D_{\zeta}\mathscr{J}$ is the shape gradient. Note that the shape
gradient depends on the normal gradient of the (primal) velocity field and the
pseudotraction,
$\big{(}-\nu\nabla\boldsymbol{v}+q\mathrm{I}\big{)}\boldsymbol{\cdot}\boldsymbol{\nu}$,
that the adjoint flow exerts on $\Gamma$.

####  2.3.3 Generalized gradients for the unknown model parameters
$\boldsymbol{x}$

The unknown model parameters $\boldsymbol{x}$ have an explicit effect on
$\mathscr{M}$ and $\mathscr{R}$, and can therefore be obtained by taking their
first variations. For the Dirichlet-type boundary condition at the inlet we
find

$\displaystyle\Big{\langle}D_{\boldsymbol{g}_{i}}\mathscr{J},~{}\boldsymbol{g}_{i}^{\prime}\Big{\rangle}_{\Gamma_{i}}$
$\displaystyle=\Big{\langle}\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{)},~{}\boldsymbol{g}_{i}^{\prime}\Big{\rangle}_{\Gamma_{i}}$
(31)
$\displaystyle=\Big{\langle}{\mathcal{C}_{\boldsymbol{g}_{i}}}\big{(}\nu\partial_{\boldsymbol{\nu}}\boldsymbol{v}-q\boldsymbol{\nu}-\eta\boldsymbol{v}\big{)}+\boldsymbol{g}_{i}-\bar{\boldsymbol{g}}_{i},~{}\boldsymbol{g}_{i}^{\prime}\Big{\rangle}_{{\mathcal{C}_{\boldsymbol{g}_{i}}}}=\Big{\langle}\widehat{D}_{\boldsymbol{g}_{i}}\mathscr{J},~{}\boldsymbol{g}_{i}^{\prime}\Big{\rangle}_{{\mathcal{C}_{\boldsymbol{g}_{i}}}},$

where $-\widehat{D}_{\boldsymbol{g}_{i}}\mathscr{J}$ is the steepest descent
direction that corresponds to the covariance-weighted norm. For the natural
boundary condition at the outlet we find

$\displaystyle\Big{\langle}D_{\boldsymbol{g}_{o}}\mathscr{J},~{}\boldsymbol{g}_{o}^{\prime}\Big{\rangle}_{\Gamma_{o}}$
$\displaystyle=\Big{\langle}\boldsymbol{v}+{\mathcal{C}^{-1}_{\boldsymbol{g}_{o}}}\big{(}\boldsymbol{g}_{o}-\bar{\boldsymbol{g}}_{o}\big{)},~{}\boldsymbol{g}_{o}^{\prime}\Big{\rangle}_{\Gamma_{o}}$
(32)
$\displaystyle=\Big{\langle}{\mathcal{C}_{\boldsymbol{g}_{o}}}\boldsymbol{v}+\boldsymbol{g}_{o}-\bar{\boldsymbol{g}}_{o},~{}\boldsymbol{g}_{o}^{\prime}\Big{\rangle}_{{\mathcal{C}_{\boldsymbol{g}_{o}}}}=\Big{\langle}\widehat{D}_{\boldsymbol{g}_{o}}\mathscr{J},~{}\boldsymbol{g}_{o}^{\prime}\Big{\rangle}_{{\mathcal{C}_{\boldsymbol{g}_{o}}}}.$

Lastly, since the kinematic viscosity is considered to be constant within
$\Omega$ its generalized gradient is

$\displaystyle\Big{\langle}D_{\nu}\mathscr{J},~{}\nu^{\prime}\Big{\rangle}_{\mathbb{R}}$
$\displaystyle=\Big{\langle}\int_{\Omega}\nabla\boldsymbol{v}\boldsymbol{:}\nabla\boldsymbol{u}+{\Sigma^{-1}_{\nu}}\big{(}\nu-\bar{\nu}\big{)},~{}\nu^{\prime}\Big{\rangle}_{\mathbb{R}}$
(33)
$\displaystyle=\Big{\langle}{\Sigma_{\nu}}\int_{\Omega}\nabla\boldsymbol{v}\boldsymbol{:}\nabla\boldsymbol{u}+\nu-\bar{\nu},~{}\nu^{\prime}\Big{\rangle}_{{\Sigma_{\nu}}}=\Big{\langle}\widehat{D}_{\nu}\mathscr{J},~{}\nu^{\prime}\Big{\rangle}_{\Sigma_{\nu}}.$

For a given step size $\mathbb{R}\ni\tau>0$, the steepest descent directions
(31)–(33) can be used either to update an unknown parameter $x$ through

$\displaystyle x_{k+1}=x_{k}+\tau s_{k},$ (34)

with $s_{k}=-\widehat{D}_{x}\mathscr{J}$, or to reconstruct an approximation
$\widetilde{H}$ of the inverse Hessian matrix, in the context of a quasi-
Newton method, and thereby to compute
$s_{k}=-\widetilde{H}\widehat{D}_{x}\mathscr{J}$. We adopt the latter
approach, which is discussed in section 2.5.

###  2.4 Geometric flow

To deform the boundary $\partial\Omega$ using the simple update formula (34)
we need a parametric surface representation. Here we choose to implicitly
represent $\partial\Omega$ using signed distance functions ${\phi_{\pm}}$. The
object $\Omega$ and its boundary $\partial\Omega$ are then identified with a
particular function ${\phi_{\pm}}$ so that

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

####  2.4.1 Implicit representation of $\Omega$ using signed distance
functions

A signed distance function ${\phi_{\pm}}$ for $\Omega$ can be obtained by
solving the Eikonal equation

$\displaystyle\lvert\nabla{\phi_{\pm}}(x)\rvert=1\quad\text{subject
to}\quad{\phi_{\pm}}\big{|}_{\partial\Omega}=0\quad,\quad x\in I.$ (35)

One way to solve this problem is with level-set methods (Osher & Sethian,
1988; Sethian, 1996; Burger, 2001, 2003; Burger & Osher, 2005; Yu et al.,
2019). There is, however, a different approach, which relies on the heat
equation (Varadhan, 1967b, a; Crane et al., 2017). The main result that we
draw from Varadhan (1967b), in order to justify the use of the heat equation
for the approximation of ${\phi_{\pm}}$, states that

$\displaystyle d(x,\partial\Omega)=\lim_{\tau_{1}\to
0}\big{(}-\frac{\sqrt{\tau_{1}}}{2}\log u(x,\tau_{1})\big{)}\quad,\quad x\in
I\quad,$ (36)

where $d(x,\partial\Omega)$ is the Euclidean distance between any point $x\in
I$ and $\partial\Omega$, and $u$ is the solution of heat propagation away from
$\partial\Omega$

$\displaystyle\left\\{\begin{aligned}
\big{(}I-\tau_{1}\rmDelta\big{)}u&=0\quad&&\textrm{in}\ I\\\
u&=1\quad&&\textrm{on}\ \partial\Omega\end{aligned}\right.\quad.$ (37)

Crane et al. (2017) used the above result to implement a smoothed distance
function computation method which they called the ‘heat method’. Here, we
slightly adapt this method to compute signed distance functions ${\phi_{\pm}}$
in truncated domains (figure 2(b)). To compute ${\phi_{\pm}}$ we therefore
solve (37) for $\tau_{1}\ll 1$, and then obtain ${\phi_{\pm}}$ by solving

$\displaystyle\left\\{\begin{aligned}
\nabla\boldsymbol{\cdot}\nabla{\phi_{\pm}}&=\nabla\boldsymbol{\cdot}X\quad&&\textrm{in}\
I\\\
\partial_{\boldsymbol{\nu}}{\phi_{\pm}}&=X\boldsymbol{\cdot}\boldsymbol{\nu}\quad&&\textrm{on}\
\partial I\\\ {\phi_{\pm}}&=0\quad&&\textrm{on}\
\partial\Omega\end{aligned}\right.\quad,\quad X=-\text{sgn}(\psi)\frac{\nabla
u}{\lvert\nabla u\rvert}\quad,$ (38)

with $X$ being the normalized heat flux and $\psi$ being a signed function
such that $\psi(x)$ is negative for points $x$ in $\Omega$ and positive for
points $x$ outside ${\Omega}$. This intermediate step (the solution of two
Poisson problems (37)-(38) instead of one) is taken to ensure that
$\lvert\nabla{\phi_{\pm}}\rvert=1$.

(a) Shape $\partial\Omega$

(b) Level-sets of ${\phi_{\pm}}$

(c) Magnitude of ${\phi_{\pm}}$ and ${\mathring{\boldsymbol{\nu}}}$

(d) Shape gradient $\zeta$ on $\partial\Omega$

(e) ${\mathring{\zeta}}$ in $I$ ($\Rey_{\zeta}=1$)

(f) ${\mathring{\zeta}}$ in $I$ ($\Rey_{\zeta}=0.01$)

Figure 2: The geometric flow of $\partial\Omega$ (figure 2(a)) relies on the
computation of its signed distance field ${\phi_{\pm}}$ (figure 2(b)) and its
normal vector extension ${\mathring{\boldsymbol{\nu}}}$ (figure 2(c)). The
shape gradient $\zeta$ (figure 2(d)), which is initially defined on
$\partial\Omega$, is extended to the whole image $I$ (${\mathring{\zeta}}$ in
figures 2(e), 2(f)). Shape regularization is achieved by increasing the
diffusion coefficient $\epsilon_{\zeta}$ in order to mitigate small scale
perturbations when assimilating noisy velocity fields
$\boldsymbol{u}^{\star}$. Figure 2(f) shows results at a lower value of
$\Rey_{\zeta}$ than figure 2(e).

####  2.4.2 Propagating the boundary of $\Omega$

To deform the boundary $\partial\Omega$ we transport ${\phi_{\pm}}$ under the
speed field ${\mathscr{V}\equiv\zeta\boldsymbol{\nu}}$. The convection-
diffusion problem for ${\phi_{\pm}}(x,t)$ reads

$\displaystyle\left\\{\begin{aligned}
\partial_{t}{\phi_{\pm}}+{\mathring{\mathscr{V}}}\boldsymbol{\cdot}\nabla{\phi_{\pm}}-\epsilon_{\phi_{\pm}}\rmDelta{\phi_{\pm}}&=0\quad&&\textrm{in}\
I\times(0,\tau]\\\ {\phi_{\pm}}&={{\phi_{\pm}}}_{0}\quad&&\textrm{in}\
I\times\\{t=0\\}\end{aligned}\right.\quad,\quad\epsilon_{\phi_{\pm}}=\frac{\lvert\mathscr{V}\rvert_{\infty}\iota}{\Rey_{\phi_{\pm}}}\quad,$
(39)

where ${{\phi_{\pm}}}_{0}$ denotes the signed distance function of the current
domain $\Omega$, $\epsilon_{\phi_{\pm}}$ is the diffusion coefficient, $\iota$
is a length scale, $\Rey_{\phi_{\pm}}$ is a Reynolds number, and
${{\mathring{\mathscr{V}}}:I\to\mathbb{R}\times\mathbb{R}}$ is an extension of
$\mathscr{V}:\partial\Omega\to\mathbb{R}\times\mathbb{R}$. If we solve (39)
for ${\phi_{\pm}}(x,\tau)$ we obtain the implicit representation of the
perturbed domain $\Omega_{\tau}$, at time $t=\tau$ (the step size), but to do
so we first need to extend $\mathscr{V}$ to the whole space of the image $I$.

To extend $\mathscr{V}$ to $I$ we extend the normal vector $\boldsymbol{\nu}$
and the scalar function $\zeta$, which are both initially defined on
$\partial\Omega$. The normal vector extension (figure 2(c)) is easily obtained
by

$\displaystyle{\mathring{\boldsymbol{\nu}}}(x)=\frac{\nabla{\phi_{\pm}}}{\lvert\nabla{\phi_{\pm}}\rvert}=\nabla{\phi_{\pm}}\quad,\quad
x\in I\quad,$ (40)

since $\lvert\nabla{\phi_{\pm}}\rvert=1$, and an outward-facing extension is
given by

$\displaystyle{\mathring{\boldsymbol{\nu}}}_{o}=\text{sgn}({\phi_{\pm}})~{}{\mathring{\boldsymbol{\nu}}}.$
(41)

We then use the extended normal vector ${\mathring{\boldsymbol{\nu}}}_{o}$ to
extend $\zeta\in L^{2}(\partial\Omega)$ to ${\mathring{\zeta}}\in L^{2}(I)$,
using the convection-diffusion problem

$\displaystyle\left\\{\begin{aligned}
\partial_{t}{\mathring{\zeta}}+{{\mathring{\boldsymbol{\nu}}}_{o}}\boldsymbol{\cdot}\nabla{\mathring{\zeta}}-\epsilon_{\zeta}\rmDelta{\mathring{\zeta}}&=0\quad&&\textrm{in}\quad
I\times(0,\tau_{\zeta}]\\\
{\mathring{\zeta}}&=\zeta\quad&&\textrm{on}\quad\partial\Omega\times(0,\tau_{\zeta}]\\\
{\mathring{\zeta}}&\equiv 0\quad&&\textrm{in}\quad
I\times\\{t=0\\}\end{aligned}\right.\quad,\quad\epsilon_{\zeta}=\frac{\lvert{\mathring{\boldsymbol{\nu}}}_{o}\rvert_{\infty}\iota}{\Rey_{\zeta}}\quad.$
(42)

In other words, we convect $\zeta$ along the predefined
${\mathring{\boldsymbol{\nu}}}_{o}$-streamlines and add isotropic diffusion
for regularization (figures 2(e), 2(f)). The choice of $\epsilon_{\phi_{\pm}}$
in (39) and $\epsilon_{\zeta}$ in (42) has been made in order for the shape
regularization to depend only on the length scale $\iota$ and the Reynolds
numbers $\Rey_{\phi_{\pm}},\Rey_{\zeta}$. More precisely, the shape
regularization depends only on $\Rey_{\phi_{\pm}}$ and $\Rey_{\zeta}$ because
we fix the length scale $\iota$ to equal the smallest possible length scale of
the modelled flow, which is the numerical grid spacing $h$ for a uniform
cartesian grid. For illustration, if we consider $\zeta$ to be the
concentration of a dye on $\partial\Omega$ (figure 2(d)), using a simplified
scaling argument similar to the growth of a boundary layer on a flat plate, we
observe that the diffusing dye at distance $d$ from $\partial\Omega$ will
extend over a width $\delta$ such that

$\displaystyle\delta\sim\sqrt{\frac{\epsilon_{\zeta}d}{\lvert{\mathring{\boldsymbol{\nu}}}_{o}\rvert_{\infty}}}=\sqrt{\frac{\iota
d}{\Rey_{\zeta}}}\quad,\quad\text{or}\quad\frac{\delta}{\iota}\sim\sqrt{\frac{\alpha}{\Rey_{\zeta}}},\quad\text{when}\quad
d=\alpha\iota\quad.$ (43)

The above scaling approximation describes the dissipation rate of small-scale
features such as roughness away from $\partial\Omega$. This is therefore how
$\Rey_{\phi_{\pm}}$ and $\Rey_{\zeta}$ control the regularity of the boundary
$\partial\Omega_{\tau}$ at time $t=\tau$, which is given by (39). We take
$\tau_{\zeta}$ to be large enough to find a steady-state for (42). We recast
the linear initial value problems (39) and (42) into their corresponding
boundary value problems using backward-Euler temporal discretization because
the time dependent solution does not interest us here.

The extended shape gradient (30), after taking into account the regularizing
term for ${\phi_{\pm}}$, is therefore given by

$\displaystyle\Big{\langle}D_{\mathring{\zeta}}\mathscr{J},~{}{\mathring{\zeta}}^{\prime}\Big{\rangle}_{I}$
$\displaystyle=\Big{\langle}{\mathring{\zeta}}+\mathcal{C}^{-1}_{\phi_{\pm}}\big{(}\bar{\phi}_{\pm}-{\phi_{\pm}}\big{)},~{}{\mathring{\zeta}}^{\prime}\Big{\rangle}_{I}$
(44)
$\displaystyle=\Big{\langle}\mathcal{C}_{\phi_{\pm}}{\mathring{\zeta}}+\bar{\phi}_{\pm}-{\phi_{\pm}},~{}{\mathring{\zeta}}^{\prime}\Big{\rangle}_{\mathcal{C}_{\phi_{\pm}}}=\Big{\langle}\widehat{D}_{\mathring{\zeta}}\mathscr{J},~{}{\mathring{\zeta}}^{\prime}\Big{\rangle}_{\mathcal{C}_{\phi_{\pm}}}\quad,$

where ${\mathring{\zeta}}$ is the extension of the shape gradient
$\zeta(x)=\partial_{\boldsymbol{\nu}}\boldsymbol{u}\boldsymbol{\cdot}\big{(}-\nu\partial_{\boldsymbol{\nu}}\boldsymbol{v}+q\boldsymbol{\nu}\big{)}$,
for $x$ on $\Gamma$.

###  2.5 Segregated approach for the Euler–Lagrange system

The inverse Navier–Stokes problem for the reconstruction and the segmentation
of noisy velocity images $\boldsymbol{u}^{\star}$ can be written as the saddle
point problem (Benzi et al., 2005)

$\text{find}\quad\boldsymbol{u}^{\circ}\equiv\mathrm{arg}\
\underset{\Omega,\boldsymbol{x}}{\mathrm{min}}\
\underset{\boldsymbol{v},q}{\mathrm{max}}~{}\mathscr{J}(\Omega)(\boldsymbol{u},p,\boldsymbol{v},q;\boldsymbol{x}),$
(45)

where $\mathscr{J}$ is given by (16). The above optimization problem leads to
an Euler–Lagrange system whose optimality conditions were formulated in
section 2.3. We briefly describe our segregated approach to solve this
Euler–Lagrange system in algorithm 1.

Input: $\boldsymbol{u}^{\star}$, initial guesses for the unknowns
$(\Omega_{0},x_{0})$, regularization parameters.

begin

$k\leftarrow 0$

$\makebox[0.0pt][l]{$({\phi_{\pm}})_{k}$}\phantom{banana}\leftarrow$ signed
distance field (eq. (37)-(38))

$\makebox[0.0pt][l]{$(\boldsymbol{u},p)_{k}$}\phantom{banana}\leftarrow$
Navier–Stokes problem for $({\phi_{\pm}},\boldsymbol{x})_{k}$ (eq. (1))

while  _convergence criterion is not met_ do

$\makebox[0.0pt][l]{($\boldsymbol{v},q$)}\phantom{avocadoss}\leftarrow$
adjoint Navier–Stokes problem with $\boldsymbol{u}_{k}$ (eq. (21))

$\makebox[0.0pt][l]{$\widehat{D}_{(\cdot)}\mathscr{J}$}\phantom{avocadoss}\leftarrow$
steepest ascent directions (eq. (31)–(33) and (44))

$\makebox[0.0pt][l]{$\boldsymbol{s},\tau$}\phantom{avocadoss}\leftarrow$
search directions and step-size (eq. (46))

$\makebox[0.0pt][l]{$({\phi_{\pm}},\boldsymbol{x})_{k+1}$}\phantom{avocadoss}\leftarrow$
perturb ${\phi_{\pm}}$ (eq. (39)) and model parameters $\boldsymbol{x}$ (eq.
(34))

$\makebox[0.0pt][l]{$(\boldsymbol{u},p)_{k+1}$}\phantom{avocadoss}\leftarrow$
linearized Navier–Stokes problem for $({\phi_{\pm}},\boldsymbol{x})_{k+1}$
(eq. (47))

$k\leftarrow k+1$

$(\boldsymbol{u}^{\circ},p^{\circ})$avocados$\leftarrow(\boldsymbol{u},p)_{k}$

$(\Omega^{\circ},\boldsymbol{x}^{\circ})$avocados$\leftarrow({\phi_{\pm}},\boldsymbol{x})_{k}$

Output: reconstruction $(\boldsymbol{u}^{\circ},p^{\circ})$ and inferred model
parameters $(\Omega^{\circ},\boldsymbol{x}^{\circ})$.

Optional output: wall shear rate $\gamma_{w}^{\circ}$ from
$\boldsymbol{u}^{\circ}$ and $\partial\Omega^{\circ}$.

Algorithm 1 Reconstruction and segmentation of noisy flow velocity images.

To precondition the steepest descent directions (31)–(33) and (44), we
reconstruct the approximated inverse Hessian $\widetilde{H}$ of each unknown
using the BFGS quasi-Newton method (Fletcher, 2000) with damping (Nocedal &
Wright, 2006). Due to the large scale of the problem, it is only possible to
work with the matrix-vector product representation of $\widetilde{H}$.
Consequently, the search directions are given by

$\displaystyle\boldsymbol{s}=-\begin{pmatrix}\widetilde{H}_{\mathring{\zeta}}&\cdot&\cdot&\cdot\\\
\cdot&\widetilde{H}_{\boldsymbol{g}_{i}}&\cdot&\cdot\\\
\cdot&\cdot&\widetilde{H}_{\boldsymbol{g}_{o}}&\cdot\\\
\cdot&\cdot&\cdot&\widetilde{H}_{\nu}\end{pmatrix}\begin{pmatrix}\widehat{D}_{\mathring{\zeta}}\mathscr{J}\\\
\widehat{D}_{\boldsymbol{g}_{i}}\mathscr{J}\\\
\widehat{D}_{\boldsymbol{g}_{o}}\mathscr{J}\\\
\widehat{D}_{\nu}\mathscr{J}\end{pmatrix}\quad,$ (46)

and the unknown variables $\boldsymbol{x}$ are updated according to (34). The
signed distance function ${\phi_{\pm}}$ is perturbed according to (39), with
${\mathring{\mathscr{V}}}\equiv-\big{(}\widetilde{H}_{\mathring{\zeta}}\widehat{D}_{\mathring{\zeta}}\mathscr{J}\big{)}{\mathring{\boldsymbol{\nu}}}$.
We start every line search with a global step size $\tau=1$, and halve the
step size until
${\mathscr{J}(({\phi_{\pm}},\boldsymbol{x})_{k+1})<\mathscr{J}(({\phi_{\pm}},\boldsymbol{x})_{k})}$.
To update the flowfield $\boldsymbol{u}_{k}$ to $\boldsymbol{u}_{k+1}$ we
solve the Oseen problem for the updated parameters
$({\phi_{\pm}},\boldsymbol{x})_{k+1}$

$\displaystyle\boldsymbol{u}_{k}\boldsymbol{\cdot}\nabla\boldsymbol{u}_{k+1}-\nu{\rmDelta}\boldsymbol{u}_{k+1}+\nabla
p_{k+1}=\boldsymbol{0},\quad\nabla\boldsymbol{\cdot}\boldsymbol{u}_{k+1}=0\quad,$
(47)

with the boundary conditions given by (1). Algorithm 1 terminates if either
the covariance-weighted norm for the perturbations of the model parameters is
below the user-specified tolerance, or the line search fails to reduce
$\mathscr{J}$.

###  2.6 Uncertainty estimation

We now briefly describe how the reconstructed inverse Hessian $\widetilde{H}$
can provide estimates for the uncertainties of the model parameters. To
simplify the description, let $x$ denote an unknown parameter distributed
according to $\mathcal{N}(x_{k},\mathcal{C}_{x})$. The linear approximation to
the data $\boldsymbol{u}^{\star}$ is given by

$\displaystyle\boldsymbol{u}^{\star}=\mathcal{Z}x+\varepsilon\quad,\quad\varepsilon\sim\mathcal{N}(0,{\mathcal{C}_{\boldsymbol{u}}})\quad,$
(48)

where $\boldsymbol{u}=\mathcal{Z}x$, where $\mathcal{Z}$ is the operator that
encodes the linearized Navier–Stokes problem around the solution
$\boldsymbol{u}_{k}$. To solve (45), we update $x$ as

$\displaystyle
x_{k+1}=x_{k}+\mathcal{C}\mathcal{Z}^{\dagger}{\mathcal{C}^{-1}_{\boldsymbol{u}}}\big{(}\boldsymbol{u}^{\star}-\mathcal{Z}x_{k}\big{)}\quad,\quad\text{with}\quad\mathcal{C}=\big{(}\mathcal{Z}^{\dagger}{\mathcal{C}^{-1}_{\boldsymbol{u}}}\mathcal{Z}+\mathcal{C}^{-1}_{x}\big{)}^{-1}\quad,$
(49)

where $\mathcal{Z}^{\dagger}$ is the operator that encodes the adjoint
Navier–Stokes problem, and $\mathcal{C}$ is the posterior covariance operator.
It can be shown that (Tarantola, 2005, Chapter 6.22.8)

$\displaystyle\mathcal{C}=\big{(}\mathcal{Z}^{\dagger}{\mathcal{C}^{-1}_{\boldsymbol{u}}}\mathcal{Z}+\mathcal{C}^{-1}_{x}\big{)}^{-1}=\big{(}\mathcal{C}_{x}\mathcal{Z}^{\dagger}{\mathcal{C}^{-1}_{\boldsymbol{u}}}\mathcal{Z}+\mathrm{I}\big{)}^{-1}\mathcal{C}_{x}\simeq\widetilde{H}_{x}\mathcal{C}_{x}\quad,$
(50)

where $\widetilde{H}_{x}$ is the reconstructed inverse Hessian for $x$. Note
that $\widetilde{H}$ by itself approximates
$\big{(}\mathcal{C}_{x}\mathcal{Z}^{\dagger}{\mathcal{C}^{-1}_{\boldsymbol{u}}}\mathcal{Z}+\mathrm{I}\big{)}^{-1}$,
and not $\mathcal{C}$, because we use the steepest ascent directions
$\widehat{D}_{(\cdot)}\mathscr{J}$ (prior-preconditioned gradients), instead
of the gradients ${D}_{(\cdot)}\mathscr{J}$, in the BFGS formula. Therefore,
if $\widetilde{\mathcal{C}}\equiv\widetilde{H}_{x}\mathcal{C}_{x}$ is the
approximated covariance matrix, then samples $x^{s}_{k+1}$ from the posterior
distribution can be drawn using the Karhunen–Loève expansion

$\displaystyle
x^{s}_{k+1}=x_{k}+\sum_{k}\eta_{k}~{}\sqrt{\lambda_{k}}\varphi_{k}\quad,\quad\text{with}\quad\eta_{k}\sim\mathcal{N}(0,1)\quad,$
(51)

where $(\lambda,\varphi)_{k}$ is the eigenvalue/eigenvector pair of
$\widetilde{\mathcal{C}}$. The variance of $x_{k+1}$ can then be directly
computed from the samples.

###  2.7 Numerics

To solve the above boundary value problems numerically, we use an immersed
boundary finite element method. In particular, we implement the fictitious
domain cut-cell finite element method (FEM), introduced by Burman (2010);
Burman & Hansbo (2012) for the Poisson problem, and later on extended to the
Stokes and the Oseen problems (Schott & Wall, 2014; Massing et al., 2014;
Burman et al., 2015; Massing et al., 2018). We define $\mathcal{T}_{h}$ to be
a tessellation of $I$ produced by square cells (pixels) $K\in\mathcal{T}_{h}$,
having sides of length $h$. We also define the set of cut-cells
$\mathcal{T}_{h}^{\triangleright}$ consisting of the cells that are cut by the
boundary $\partial\Omega$, and $\mathcal{T}_{h}^{\scaleobj{0.7}{\Box}}$ the
set of cells that are found inside $\Omega$ and which remain intact (not cut)
(see figure 1). We assume that the boundary $\partial\Omega$ is well-resolved,
i.e. $\ell_{\partial\Omega}/h\gg 1$ where $\ell_{\partial\Omega}$ is the
smallest length scale of $\partial\Omega$. For the detailed assumptions on
$\partial\Omega$ we cite Burman & Hansbo (2012). The discretized space is
generated by assigning a bilinear quadrilateral finite element
$\mathcal{Q}_{1}$ to every cell $K$. To compute the integrals we use standard
Gaussian quadrature for cells $K\in\mathcal{T}_{h}^{\scaleobj{0.7}{\Box}}$,
while for cut-cells $K\in\mathcal{T}_{h}^{\triangleright}$, where integration
must be considered only for the intersection $K\cap\Omega$, we use the
approach of Mirtich (1996), which relies on the divergence theorem and simply
replaces the integral over $K\cap\Omega$ with an integral over
$\partial\big{(}K\cap\Omega\big{)}$. The boundary integral on
$\partial\big{(}K\cap\Omega\big{)}$ is then easily computed using one-
dimensional Gaussian quadrature (Massing et al., 2013). Since we use an inf-
sup unstable finite element pair ($\mathcal{Q}_{1}$-$\mathcal{Q}_{1}$)
(Brenner, Susanne, Scott, 2008) we use a pressure-stabilizing Petrov-Galerkin
formulation (Tezduyar, 1991; Codina, 2002) and $\nabla$-div stabilization for
preconditioning (Benzi & Olshanskii, 2006; Heister & Rapin, 2013). Typical
values and formulas for numerical parameters, e.g. Nitsche’s penalization
$\eta$, are given by Massing et al. (2014, 2018). Here, we take
$\eta=\gamma\nu/h$ (Massing et al., 2018), with $\gamma=100$. To solve the
Navier–Stokes problem we use fixed-point iteration (Oseen linearization), and
at each iteration we solve the coupled system using the Schur complement; with
an iterative solver (LGMRES) for the outer loops, and a direct sparse solver
(UMFPACK) for the inner loops. The immersed FEM solver, and all the necessary
numerical operations of algorithm 1, are implemented in Python, using its
standard libraries for scientific computing, namely SciPy (Virtanen et al.,
2020) and NumPy (Harris et al., 2020). Computationally intensive functions are
accelerated using Numba (Lam et al., 2015) and CuPy (Okuta et al., 2017).

##  3 Reconstruction and segmentation of flow images

In this section we reconstruct and segment noisy flow images by solving the
inverse Navier–Stokes problem (45) using algorithm 1. We then use the
reconstructed velocity field to estimate the wall shear rate on the
reconstructed boundary. First, we apply this to three test cases with known
solutions by generating synthetic 2D Navier–Stokes data. Next, we perform a
magnetic resonance velocimetry experiment in order to acquire images of a 3D
axisymmetric Navier–Stokes flow, and apply algorithm 1 to these images.

We define the signal-to-noise ratio (SNR) of the $u^{\star}_{x}$ image as

$\displaystyle\text{SNR}_{x}=\frac{\mu_{x}}{\sigma_{u_{x}}}\quad,\quad\mu_{x}\equiv\frac{1}{\lvert\Omega^{\bullet}\rvert}\int_{\Omega^{\bullet}}\lvert
u_{x}^{\bullet}\rvert\quad,$ (52)

where $\sigma_{u_{x}}$ is the standard deviation, $\Omega^{\bullet}$ is the
ground truth domain, $\lvert\Omega^{\bullet}\rvert$ is the volume of this
domain, and $\lvert u_{x}^{\bullet}\rvert$ is the magnitude of the ground
truth $x$-velocity component in $\Omega^{\bullet}$. We also define the
componentwise averaged, noise relative reconstruction error
$\mathscr{E}^{\bullet}_{x}$, and the total relative reconstruction error
$\mathcal{E}^{\bullet}$ by

$\displaystyle\mathscr{E}^{\bullet}_{x}\equiv\log\bigg{(}\frac{1}{\lvert\Omega\rvert}\int_{\Omega}\frac{\lvert
u^{\bullet}_{x}-\mathcal{S}u^{\circ}_{x}\rvert}{\sigma_{u_{x}}}\bigg{)}\quad\text{and}\quad\mathcal{E}^{\bullet}\equiv\frac{\big{\lVert}\boldsymbol{u}^{\bullet}-\mathcal{S}\boldsymbol{u}^{\circ}\big{\rVert}_{L^{1}(I)}}{\big{\lVert}\boldsymbol{u}^{\bullet}\big{\rVert}_{L^{1}(I)}}\quad.$
(53)

respectively. Similar measures also apply for the $u^{\star}_{y}$ image.

We define the volumetric flow rate $Q$, the cross-section area at the inlet
$A$, and the diameter at the inlet $D$. The Reynolds number is based on the
reference velocity ${U\equiv Q/A}$, and the reference length $D$.

###  3.1 Synthetic data for 2D flow in a converging channel

We start by testing algorithm 1 on a flow through a symmetric converging
channel having a taper ratio of $0.67$. To generate synthetic 2D Navier–Stokes
data we solve the Navier–Stokes problem (1) for a parabolic inlet velocity
profile ($\boldsymbol{g}_{i}$), zero-pseudotraction boundary conditions at the
outlet ($\boldsymbol{g}_{o}\equiv\boldsymbol{0}$), and $\Rey\simeq 534$, in
order to obtain the ground truth velocity $\boldsymbol{u}^{\bullet}$. We then
generate the synthetic data $\boldsymbol{u}^{\star}$ by corrupting the
components of $\boldsymbol{u}^{\bullet}$ with white Gaussian noise such that
$\text{SNR}_{x}=\text{SNR}_{y}=3$. For this test case, we are only trying to
infer $\Omega$ and $\boldsymbol{g}_{i}$. Note that, in our method, the initial
guess $x_{0}$ of an unknown $x$ equals the mean of its prior distribution
$\bar{x}$, i.e. $x_{0}\equiv\bar{x}$. We start the algorithm using bad initial
guesses (high uncertainty in priors) for both the unknown parameters (see
table 1). The initial guess for $\Omega$, labelled $\Omega_{0}$, is a
rectangular domain with height equal to $0.7D$, centered in the image domain.
For ${\boldsymbol{g}_{i}}_{0}$ we take a parabolic velocity profile with a
peak velocity of approximately $2U$ that fits the inlet of $\Omega_{0}$. For
comparison, $\boldsymbol{g}_{i}^{\bullet}$ has a peak velocity of $1.5U$,
while it is also defined on a different domain, namely $\Omega^{\bullet}$.

|  | image dimension | model dimension | $\sigma_{u_{x}}/U$ | $\sigma_{u_{y}}/U$
---|---|---|---|---|---
converging channel | (2D) | $192^{2}$ | $200^{2}$ | $3.97\times 10^{-1}$ | $8.97\times 10^{-3}$
Regularization | $\sigma_{\phi_{\pm}}/D$ | $\sigma_{\boldsymbol{g}_{i}}/U$ | $\sigma_{\nu}/UD$ | $\Rey_{\phi_{\pm}}$ | $\Rey_{\zeta}$ | $\ell/h$
converging channel | (2D) | $1.0$ | $2.0$ | $\cdot$ | 0.025 | 0.025 | 3
Table 1: Input parameters for the inverse 2D Navier–Stokes problem.

(a) Synthetic image $u^{\star}_{x}$

(b) Our reconstruction $u_{x}^{\circ}$

(c) Discrepancy
$\sigma_{u_{x}}^{-1}\big{(}u^{\star}_{x}-\mathcal{S}u_{x}^{\circ}\big{)}$

(d) Synthetic image $u^{\star}_{y}$

(e) Our reconstruction $u_{y}^{\circ}$

(f) Discrepancy
$\sigma_{u_{y}}^{-1}\big{(}u^{\star}_{y}-\mathcal{S}u_{y}^{\circ}\big{)}$

Figure 3: Reconstruction (algorithm 1) of synthetic noisy velocity images
depicting the flow (from left to right) in a converging channel. Figures
3(a)-3(b) and 3(d)-3(e) show the horizontal, $u_{x}$, and vertical, $u_{y}$,
velocities and share the same colormap (colorbar not shown). Figures 3(c) and
3(f) show the discrepancy between the noisy velocity images and the
reconstruction (colorbars apply only to figures 3(c) and 3(f)).

(a) Velocity magnitude $\lvert\boldsymbol{u}^{\circ}\rvert$ and shape
$\partial\Omega^{\circ}$

(b) Reconstruction error history

(c) Noisy data (grey) and reconstruction

(d) Ground truth velocity distributions

Figure 4: Reconstruction (algorithm 1) of synthetic images depicting the flow
(from left to right) in a converging channel. Figure 4(a) depicts the
reconstructed boundary $\partial\Omega^{\circ}$ (cyan line), the $2\sigma$
confidence region computed from the approximated posterior covariance
$\widetilde{\mathcal{C}}_{\mathring{\zeta}}\equiv\widetilde{H}_{\mathring{\zeta}}\mathcal{C}_{\phi_{\pm}}$
(blue region), the ground truth boundary $\partial\Omega^{\bullet}$ (yellow
line), and the initial guess $\partial\Omega_{0}$ (white line). Figure 4(b)
shows the reconstruction error as a function of iteration number. Velocity
slices are drawn for 10 equidistant cross-sections (labelled with the letters
A to J) for both the reconstructed images (figure 4(c)) and the ground truth
(figure 4(d)), colored red for positive values and blue for negative.

(a) Our reconstruction $p^{\circ}$

(b) Ground truth $p^{\bullet}$

Figure 5: (a) Reconstructed and (b) ground truth reduced hydrodynamic pressure
($p$) for the flow (from left to right) in the converging channel in figure 4.

The algorithm manages to reconstruct and segment the noisy flow images in 39
iterations, with total reconstruction error $\mathcal{E}^{\bullet}\simeq
1.44\%$. The results are presented in figures 3 and 4. We observe that the
inverse Navier–Stokes problem performs very well in filtering the noise
$(\boldsymbol{u}^{\star}-\mathcal{S}\boldsymbol{u}^{\circ})$ (figures 3(c),
3(f)), providing noiseless images for each component of the velocity (figures
3(b), 3(e)). As we expect, the discrepancies
${\mathcal{C}^{-1}_{\boldsymbol{u}}}(\boldsymbol{u}^{\star}-\mathcal{S}\boldsymbol{u}^{\circ})$
(figures 3(c), 3(f)) consist mainly of Gaussian white noise, except at the
corners of the outlet (figure 3(f)), where there is a weak correlation. For a
more detailed presentation of the denoising effect, we plot slices of the
reconstructed velocity (figure 4(c)) and the ground truth velocity (figure
4(d)). The reconstructed pressure $p^{\circ}$, which is consistent with the
reconstructed velocity $\boldsymbol{u}^{\circ}$ to machine precision accuracy,
is, in effect, indistinguishable from the ground truth $p^{\bullet}$ (figure
5).

Having obtained the reconstructed velocity $\boldsymbol{u}^{\circ}$, we can
compute the wall shear rate $\gamma_{w}^{\circ}$ on the reconstructed boundary
$\partial\Omega^{\circ}$, which we compare with the ground truth
$\gamma_{w}^{\bullet}$ in figure 6. Using the upper
($\partial\Omega^{\circ}_{+}$) and lower ($\partial\Omega^{\circ}_{-}$) limits
of the $2\sigma$ confidence region for $\partial\Omega^{\circ}$ (figure 4(a))
we estimate a confidence region for $\gamma_{w}^{\circ}$; although this has to
be interpreted carefully. Note that, for example, $\partial\Omega^{\circ}_{+}$
and $\partial\Omega^{\circ}_{-}$ can be smoother than the mean
$\partial\Omega^{\circ}$, and, therefore, $\partial\Omega^{\circ}$ may be
found outside this confidence region. A better estimate of the confidence
region could be obtained by sampling the posterior distribution of
$\partial\Omega^{\circ}$ in order to solve a Navier–Stokes problem for each
sample $\partial\Omega_{k}$ and find the distribution of $\gamma^{\circ}_{w}$.
Since the latter approach would be computationally intensive, we only provide
our estimate, which requires the solution of only two Navier–Stokes problems.

(a) Lower boundary

(b) Upper boundary

Figure 6: Wall shear rate
$\gamma_{w}\equiv\boldsymbol{\tau}\boldsymbol{\cdot}\partial_{\boldsymbol{\nu}}\boldsymbol{u}$,
where $\boldsymbol{\tau}$ is the unit tangent vector of $\partial\Omega$, for
the converging channel flow in figure 4. The wall shear stress is found by
multiplying this by the viscosity. The reconstructed wall shear rate
($\gamma_{w}^{\circ}$) is calculated on $\partial\Omega^{\circ}$ and for
$\boldsymbol{u}^{\circ}$, while the ground truth ($\gamma_{w}^{\bullet}$) is
calculated on $\partial\Omega^{\bullet}$ and for $\boldsymbol{u}^{\bullet}$.
The blue region is bounded by the two wall shear rate distributions for
$\boldsymbol{u}^{\circ}$, calculated on the upper
($\partial\Omega^{\circ}_{+}$) and lower ($\partial\Omega^{\circ}_{-}$) limits
of the $2\sigma$ confidence region of $\partial\Omega^{\circ}$. Note that the
reconstructed solution can sometimes be found outside the blue region because
the reconstructed shape $\partial\Omega^{\circ}$ may be less regular than
$\partial\Omega^{\circ}_{+}$ or $\partial\Omega^{\circ}_{-}$.

###  3.2 Synthetic data for 2D flow in a simulated abdominal aortic aneurysm

Next, we test algorithm 1 in a channel that resembles the cross-section of a
small abdominal aortic aneurysm, with $D_{\text{max}}/D\simeq 1.5$, where
$D_{\text{max}}$ is the maximum diameter at the midsection. We generate
synthetic images for $\boldsymbol{u}^{\star}$ as in section 3.1, again for
${\text{SNR}_{x}=\text{SNR}_{y}=3}$, but now for $\Rey=153$. The ground truth
domain $\Omega^{\bullet}$ has horizontal symmetry but the inlet velocity
profile deliberately breaks this symmetry. The inverse problem is the same as
that in section 3.1 but with different input parameters (see table 2). The
initial guess $\Omega_{0}$ is a rectangular domain with height equal to
$0.85D$, centered in the image domain. For ${\boldsymbol{g}_{i}}_{0}$ we take
a skewed parabolic velocity profile with a peak velocity of approximately $2U$
that fits the inlet of $\Omega_{0}$.

|  | image dim. | model dim. | $\sigma_{u_{x}}/U$ | $\sigma_{u_{y}}/U$
---|---|---|---|---|---
simul. abd. aortic aneurysm | (2D) | $192^{2}$ | $200^{2}$ | $2.80\times 10^{-1}$ | $5.26\times 10^{-3}$
Regularization | $\sigma_{\phi_{\pm}}/D$ | $\sigma_{\boldsymbol{g}_{i}}/U$ | $\sigma_{\nu}/UD$ | $\Rey_{\phi_{\pm}}$ | $\Rey_{\zeta}$ | $\ell/h$
simul. abd. aortic aneurysm | (2D) | $1.0$ | $2.0$ | $\cdot$ | 0.1 | 0.1 | 3
Table 2: Input parameters for the inverse 2D Navier–Stokes problem.

(a) Synthetic image $u^{\star}_{x}$

(b) Our reconstruction $u_{x}^{\circ}$

(c) Discrepancy
$\sigma_{u_{x}}^{-1}\big{(}u^{\star}_{x}-\mathcal{S}u_{x}^{\circ}\big{)}$

(d) Synthetic image $u^{\star}_{y}$

(e) Our reconstruction $u_{y}^{\circ}$

(f) Discrepancy
$\sigma_{u_{y}}^{-1}\big{(}u^{\star}_{y}-\mathcal{S}u_{y}^{\circ}\big{)}$

Figure 7: As for figure 3, but for the synthetic images depicting the flow
(from left to right) in the simulated 2D model of an abdominal aortic
aneurysm.

(a) Velocity magnitude $\lvert\boldsymbol{u}^{\circ}\rvert$ and shape
$\partial\Omega^{\circ}$

(b) Reconstruction error history

(c) Noisy data (grey) and reconstruction

(d) Ground truth velocity distributions

Figure 8: As for figure 4, but for the synthetic images depicting the flow
(from left to right) in the simulated 2D model of an abdominal aortic
aneurysm.

(a) Our reconstruction $p^{\circ}$

(b) Ground truth $p^{\bullet}$

Figure 9: (a) Reconstructed and (b) ground truth reduced hydrodynamic pressure
($p$) for the flow (from left to right) in the simulated 2D model of an
abdominal aortic aneurysm in figure 8.

The algorithm manages to reconstruct and segment the noisy flow images in 39
iterations, with total reconstruction error $\mathcal{E}^{\bullet}\simeq
2.87\%$. The results are presented in figures 7 and 8. We observe that the
discrepancy (figures 7(c), 7(f)) consists mainly of Gaussian white noise.
Again, some correlations are visible in the discrepancy of the $y$-velocity
component at the upper inlet corner and the upper boundary of the simulated
abdominal aortic aneurysm. The latter correlations (figure 7(f)) can be
explained by the associated uncertainty in the predicted shape
$\partial\Omega^{\circ}$ (figure 8(a)), which is well estimated for the upper
boundary but slightly underestimated for the upper inlet corner. It is
interesting to note that the upward skewed velocity profile at the inlet
creates a region of low velocity magnitude on the lower boundary. The velocity
profiles in this region produce low wall shear stresses, as seen in figures
8(c), 8(d), and 10(b). These conditions are particularly challenging when one
tries to infer the true boundary $\partial\Omega^{\bullet}$ because the local
SNR is low ($\text{SNR}\ll 1$), meaning that there is considerable information
loss there. Despite the above difficulties, algorithm 1 manages to approximate
the posterior distribution of $\partial\Omega^{\circ}$ well, and successfully
predicts extra uncertainty in this region (figure 8(a)). Again, the
reconstructed pressure $p^{\circ}$ is indistinguishable from the ground truth
$p^{\bullet}$ (figure 9).

Using the reconstructions $\boldsymbol{u}^{\circ}$ and
$\partial\Omega^{\circ}$ we compute the wall shear rate and we compare it with
the ground truth in figure 10(b). We observe that the reconstructed solution
approximates the ground truth well, even for very low signal-to-noise ratios
($\text{SNR}=3$). Note that the waviness of the ground truth
$\gamma_{w}^{\bullet}$ is due to the relatively poor resolution of the level
set function that we intentionally used to implicitly define this domain.

(a) Lower boundary

(b) Upper boundary

Figure 10: As for figure 6 but for the synthetic images depicting the flow
(from left to right) in the simulated 2D model of an abdominal aortic aneurysm
in figure 8.

###  3.3 Synthetic data for 2D flow in a simulated aortic aneurysm

Next, we test algorithm 1 in a channel that resembles the cross-section of an
aorta that has an aneurysm in its ascending part. This test case is designed
to demonstrate that the algorithm is applicable to realistic geometries with
multiple inlets/outlets and for abnormal flow conditions (e.g. separation and
recirculation zones). We generate synthetic images for
$\boldsymbol{u}^{\star}$ as in section 3.1, but for
${\text{SNR}_{x}=\text{SNR}_{y}=2.5}$, and for $\Rey=500$. For increased
Reynolds numbers ($\Rey=1000,1500$), we observed vortex shedding within the
aneurysm and we could not find a steady flow solution to generate synthetic
images of steady flow. The inverse problem is the same as that in section 3.1
but with different input parameters (see table 3). The initial guess for the
boundary of $\Omega_{0}$ (figure 11(a)) is generated by using the Chan–Vese
segmentation method (Chan & Vese, 2001; Getreuer, 2012a; van der Walt et al.,
2014) on the noisy mask of the ground truth domain $\Omega^{\bullet}$ (figure
11(b)). The prior standard deviation $\sigma_{\phi_{\pm}}$ corresponds to the
length of approximately $7$ pixels of the noisy mask. The initial guess for
the inlet velocity profile ${\boldsymbol{g}_{i}}_{0}$ is also shown in figure
11(a). Using the prior information of the boundary and the inlet velocity
profile, algorithm 1 generates an initial guess for the Navier–Stokes velocity
field (figures 12(a), 12(b)) during its zeroth iteration.

|  | image dim. | model dim. | $\sigma_{u_{x}}/U$ | $\sigma_{u_{y}}/U$
---|---|---|---|---|---
simul. aortic aneurysm | (2D) | $300^{2}$ | $325^{2}$ | $1.17\times 10^{-1}$ | $2.62\times 10^{-1}$
Regularization | $\sigma_{\phi_{\pm}}/D$ | $\sigma_{{g_{i}}_{x}}/U$ | $\sigma_{{g_{i}}_{y}}/U$ | $\sigma_{\nu}/UD$ | $\Rey_{\phi_{\pm}}$ | $\Rey_{\zeta}$ | $\ell/h$
simul. aortic aneurysm | (2D) | $0.025$ | $0.5$ | $0.4$ | . | $1$ | $1$ | $5$
Table 3: Input parameters for the inverse 2D Navier–Stokes problem.

(a) Initial guesses (priors) for $\partial\Omega$ and $\boldsymbol{g}_{i}$

(b) Noisy mask of $\Omega^{\bullet}$

Figure 11: Initial guesses (input for algorithm 1) for the geometry
($\partial\Omega_{0}$) and the inlet velocity profile
(${\boldsymbol{g}_{i}}_{0}$) versus their corresponding ground truth (figure
11(a)), for the flow in the simulated 2D model of an aortic aneurysm. The
initial guess $\partial\Omega_{0}$ (figure 11(a)) is generated by segmenting
the noisy mask (figure 11(b)) of the ground truth domain $\Omega^{\bullet}$.

(a) $(u_{x})_{0}$

(b) $(u_{y})_{0}$

(c) $\boldsymbol{u}_{0}$

(d) $\sigma_{u_{x}}^{-1}\big{(}u^{\star}_{x}-\mathcal{S}(u_{x})_{0}\big{)}$

(e) $\sigma_{u_{y}}^{-1}\big{(}u^{\star}_{y}-\mathcal{S}(u_{y})_{0}\big{)}$

(f)
${\mathcal{C}^{-1}_{\boldsymbol{u}}}(\boldsymbol{u}^{\star}-\mathcal{S}\boldsymbol{u}_{0})$

Figure 12: Zeroth iteration (N–S solution for the initial guesses in figure
11) velocity images (figures 12(a), 12(b)), streamlines (figure 12(c)), and
discrepancies with the data (figures 12(d)-12(f)), for the flow in the
simulated 2D model of an aortic aneurysm. Streamlines are plotted on top of
the velocity/discrepancy magnitude image, and streamline thickness increases
as the velocity magnitude increases. Figures 12(a)-12(c) (colorbar not shown)
and 12(d)-12(f) (colorbar shown on the right) share the same colormap.

(a) Synthetic image $u^{\star}_{x}$

(b) Our reconstruction $u_{x}^{\circ}$

(c) Discrepancy
$\sigma_{u_{x}}^{-1}\big{(}u^{\star}_{x}-\mathcal{S}u_{x}^{\circ}\big{)}$

(d) Synthetic image $u^{\star}_{y}$

(e) Our reconstruction $u_{y}^{\circ}$

(f) Discrepancy
$\sigma_{u_{y}}^{-1}\big{(}u^{\star}_{y}-\mathcal{S}u_{y}^{\circ}\big{)}$

Figure 13: Reconstruction (final iteration of algorithm 1) of synthetic noisy
velocity images depicting the flow in the simulated 2D model of an aortic
aneurysm. Figures 13(a)-13(b) and 13(d)-13(e) show the horizontal, $u_{x}$,
and vertical, $u_{y}$, velocities and share the same colormap (colorbar not
shown). Figures 13(c) and 13(f) show the discrepancy between the noisy
velocity images and the reconstruction (colorbars apply only to figures 13(c)
and 13(f)).

(a) Velocity magnitude $\lvert\boldsymbol{u}^{\circ}\rvert$ and shape
$\partial\Omega^{\circ}$

(b) Reconstruction error history

Figure 14: Reconstruction (final iteration of algorithm 1) of synthetic images
depicting the flow in the simulated 2D model of an aortic aneurysm. Figure
14(a) depicts the reconstructed boundary $\partial\Omega^{\circ}$ (cyan line),
the $2\sigma$ confidence region computed from the approximated posterior
covariance
$\widetilde{\mathcal{C}}_{\mathring{\zeta}}\equiv\widetilde{H}_{\mathring{\zeta}}\mathcal{C}_{\phi_{\pm}}$
(blue region), and the ground truth boundary $\partial\Omega^{\bullet}$
(yellow line). Figure 14(b) shows the reconstruction error as a function of
iteration number.

(a) Zeroth iteration $p_{0}$

(b) Our reconstruction $p^{\circ}$

(c) Ground truth $p^{\bullet}$

Figure 15: (a) Zeroth iteration (N–S solution for the initial guesses in
figure 11), (b) reconstructed (final iteration of algorithm 1), and (c) ground
truth reduced hydrodynamic pressure for the simulated 2D model of an aortic
aneurysm in figure 14. All subfigures share the same colormap (symmetric
logarithmic scale) and the same colorbar.

(a) Synthetic data $\boldsymbol{u}^{\star}$

(b) Our reconstruction $\boldsymbol{u}^{\circ}$

(c) Ground truth $\boldsymbol{u}^{\bullet}$

(d) TV-B $\lambda/\lambda_{0}=0.1$

(e) TV-B $\lambda/\lambda_{0}=0.01$

(f) TV-B $\lambda/\lambda_{0}=0.001$

Figure 16: Streamlines for the flow in the simulated 2D model of an aortic
aneurysm (figures 13 and 14), and comparison with total variation denoising
using Bregman iteration (TV-B) with different weights $\lambda$. Streamlines
are plotted on top of the velocity magnitude image, and streamline thickness
increases as the velocity magnitude increases.

The algorithm manages to reconstruct and segment the noisy flow images in 15
iterations, with total reconstruction error $\mathcal{E}^{\bullet}\simeq
5.73\%$. The results are presented in figures 13 and 14. We observe that the
discrepancy of the last iteration (figures 13(c), 13(f)) consists mainly of
Gaussian white noise. Some correlations are visible in the discrepancy of the
$x$-velocity component near the stagnation points of the upper branches, but
these correlations are explained by the extra uncertainty in the predicted
shape $\partial\Omega^{\circ}$ (figure 14(a)). By comparing figures 13(c) and
13(f) with figures 12(d) and 12(e), we confirm that the algorithm has
successfully assimilated the remaining information from the noisy velocity
measurements. Figure 15 shows the pressure of the zeroth iteration (figure
15(a)), and the reconstructed pressure $p^{\circ}$ (figure 15(b)), which
compares well to the ground truth pressure $p^{\bullet}$ (figure 15(c)).

We further compare the performance of algorithm 1 with a state-of-the-art
image denoising algorithm, namely total variation denoising using Bregman
iteration (TV-B) (Getreuer, 2012b; van der Walt et al., 2014), in figure 16.
We first observe that algorithm 1 denoises the velocity field without losing
contrast near the walls of the aorta, and accurately identifies the low-speed
vortical structure within the aneurysm, which is obscured by noise. We then
test three different values of the TV-B parameter $\lambda/\lambda_{0}$444The
parameter $\lambda_{0}=\lambda_{0}(\sigma)$, where $\sigma$ is the noise
standard deviation in the image, is given by Getreuer (2012b) as an optimal
value for $\lambda$., which controls the total variation regularization, and
observe that, even though TV-B manages to denoise the velocity field and
reveal certain large scale vortices, there is considerable loss of contrast
near the walls of the aorta and a systematic error (e.g. decreasing peak
velocity) that increases as $\lambda$ decreases.

Using the reconstructions $\boldsymbol{u}^{\circ}$ and
$\partial\Omega^{\circ}$ we compute the reconstructed wall shear rate
($\gamma_{w}^{\circ}$) and compare it with the ground truth
($\gamma_{w}^{\bullet}$) (figure 17). We observe that $\gamma_{w}^{\circ}$
approximates $\gamma_{w}^{\bullet}$ well, and that discrepancies are well
accounted for by the $\gamma_{w}^{\circ}\pm 2\sigma$-bounds.

(a) Zeroth iteration $(\gamma_{w})_{0}$

(b) Our reconstruction $\gamma^{\circ}_{w}$

(c) Ground truth $\gamma^{\bullet}_{w}$

(d) Lower conf. bound $\gamma^{\circ}_{w}-2\sigma$

(e) Upper conf. bound $\gamma^{\circ}_{w}+2\sigma$

Figure 17: Wall shear rate
$\gamma_{w}\equiv\boldsymbol{\tau}\boldsymbol{\cdot}\partial_{\boldsymbol{\nu}}\boldsymbol{u}$,
where $\boldsymbol{\tau}$ is the unit tangent vector of $\partial\Omega$, for
the flow in the simulated 2D model of an aortic aneurysm in figure 14. The
wall shear stress is found by multiplying this by the viscosity. The
reconstructed wall shear rate ($\gamma_{w}^{\circ}$) is calculated on
$\partial\Omega^{\circ}$ and for $\boldsymbol{u}^{\circ}$, while the ground
truth ($\gamma_{w}^{\bullet}$) is calculated on $\partial\Omega^{\bullet}$ and
for $\boldsymbol{u}^{\bullet}$. The $\pm 2\sigma$-bounds are calculated on the
upper ($\partial\Omega^{\circ}_{+}$) and lower ($\partial\Omega^{\circ}_{-}$)
limits of the confidence region of $\partial\Omega^{\circ}$. All subfigures
share the same colormap (symmetric logarithmic scale) and the same colorbar.

###  3.4 Magnetic resonance velocimetry experiment

We measured the flow through a converging nozzle using magnetic resonance
velocimetry (Fukushima, 1999; Mantle & Sederman, 2003; Elkins & Alley, 2007).
The nozzle converges from an inner diameter of 25mm to an inner diameter of
13mm, over a length of 40mm (figure 18(b)). On either side of the converging
section, the entrance-to-exit length equals 10 times the local diameter
(figure 18(b)) in order to ensure the absence of entrance/exit effects. We
acquired velocity images for a Reynolds number of 162 (defined at the nozzle
outlet). We used a 40 wt% glycerol in water solution (Cheng, 2008; Volk &
Kähler, 2018) as the working fluid in order to increase the viscosity and
minimize the effect of thermal convection in the resulting velocity field due
to the temperature difference between the magnet bore and the working fluid.
The nozzle is made of polyoxymethylene to minimize magnetic susceptibility
differences between the nozzle wall and the working fluid (Wapler et al.,
2014). Figure 18(a) depicts the schematic of the flow loop of the MRV
experiment. To pump the water/glycerol solution we used a Watson Marlow 505S
peristaltic pump (Watson Marlow, Falmouth UK) with a 2L dampening vessel at
its outlet to dampen flow oscillations introduced by the peristaltic pump. To
make the flow uniform, we installed porous polyethylene distributor plates
(SPC technologies, Fakenham UK) at the entrance and the exit of the nozzle.

We acquired the velocity images on a Bruker Spectrospin DMX200 with a $4.7$T
superconducting magnet, which is equipped with a gradient set providing
magnetic field gradients of a maximum strength of 13.1Gcm-1 in three
orthogonal directions, and a birdcage radiofrequency coil tuned to a
${}^{1}\mathrm{H}$ frequency of 199.7 MHz with a diameter and a length of
6.3cm. To acquire 2D velocity images we used slice-selective spin-echo imaging
(Edelstein et al., 1980) combined with pulsed gradient spin-echo (PGSE)
(Stejskal & Tanner, 1965) for motion encoding (figure 18(c)). We measured each
of the three orthogonal velocity components in a 1mm thick transverse slice
through the converging section of the nozzle, which is centered along the
nozzle centerline. The flow images we acquired have a field of view of
84.2$\times$28.6mm at 512$\times$128 pixels, giving an in-plane resolution of
165$\times$223$\mu$m. For velocity measurements in the net flow direction, we
used a gradient pulse duration, $\delta$, of 0.3 to 0.5ms and flow observation
times, ${\rmDelta}$, of 9 to 12ms. For velocity measurements in the
perpendicular to the net flow direction, we used an increased gradient pulse
duration, $\delta$, of 1.0ms and an increased observation time, ${\rmDelta}$,
of 25 to 30ms, due to the lower velocity magnitudes in this direction. We set
the amplitude, $g$, of the flow encoding gradient pulses to $\pm$3Gcm-1 for
the direction parallel to the net flow and to $\pm$1.5Gcm-1 for the direction
perpendicular to the net flow, in order to maximize phase contrast whilst
avoiding velocity aliasing by phase wrapping. To obtain an image for each
velocity component, we took the phase difference between two images acquired
with flow encoding gradients having equal magnitude $g$ but opposite signs. To
remove any phase shift contributions that are not caused by the flow, we
corrected the measured phase shift of each voxel by subtracting the phase
shift measured under zero-flow conditions. The gradient stabilization time
that we used is 1ms and we acquired the signal with a sweep width of 100kHz.
We used hard 90∘ excitation pulses with a duration of 85$\mu$s, and a
512$\mu$s Gaussian-shaped soft 180∘ pulse for slice selection and spin-echo
refocusing. We found the $T_{1}$ relaxation time of the glycerol solution to
be 702ms, as measured by an inversion recovery pulse sequence. To allow for
magnetization recovery between the acquisitions, we used a repetition time of
1.0s. To eliminate unwanted coherences and common signal artefacts, such as DC
offset, we used a four step phase cycle.

To be consistent with the standard definition used in MRI/MRV, we define the
SNR of each MRV image using (52), but with $\mu_{x}$ replaced by the mean
signal intensity (images of the ${}^{1}\mathrm{H}$ spin density) over the
nozzle domain ($\mu_{I}$), and $\sigma_{u_{x}}$ replaced by the standard
deviation of the Rayleigh distributed noise in a region with no signal
($\sigma_{I}$) (Gudbjartsson & Patz, 1995). The standard deviation for the
phase is therefore ${\sigma_{\varphi}=1/\text{SNR}}$. The MRV images are
acquired by taking the sum/difference of four phase images, and then
multiplying by the constant factor $1/2\gamma g\delta\Delta$, where $\gamma$
is the gyromagnetic ratio of ${}^{1}\mathrm{H}$ (linear relation between the
image phase and the velocity). The error in the MRV measured velocity is
therefore $\sigma_{u}=\sigma_{\varphi}/\gamma g\delta\Delta$. To acquire high
SNR images (figure 19), we averaged 32 scans, resulting in a total acquisition
time of 137 minutes per velocity image ($\sim 4.6$ hours for both velocity
components). To evaluate the denoising capability of the algorithm we acquired
poor SNR images by averaging only 4 scans (the minimum requirement for a full
phase cycle) and decreasing the repetition time to 300ms, resulting in a total
acquisition time of 5.1 minutes per velocity image ($10.2$ minutes for both
velocity components).

To verify the quantitative nature of the MRV experiment we compared the
volumetric flow rates calculated from the MRV images (using 2D slice-selective
velocity imaging in planes normal to the direction of net flow) with the
volumetric flow rates measured from the pump outlet. The results agree with an
average error of $\pm$1.8%.

(a) Magnet and flow loop

(b) Converging nozzle

(c) Spin-echo pulse sequence with slice selective refocusing and flow encoding

Figure 18: Schematic of the rig that we use to conduct the MRV experiment
consisting of: (1) 20L holding tank, (2) peristaltic pump, (3) 2L vessel, (4)
clamp valves, (5) porous polyethylene distributor, (6) radiofrequency probe,
(7) converging nozzle, (8) $4.7$T superconducting magnet, (9) volumetric
cylinder for flow measurements. Figure 18(b) shows a sketch of the converging
nozzle with the active area of the spectrometer shown by a red box. The pulse
sequence that we use for 2D velocity imaging is shown in figure 18(c).

(a) $u^{\bullet}_{z}$

(b) $u^{\bullet}_{r}$

Figure 19: High SNR images (average of 32 scans with $\text{SNR}_{z}\simeq
44,\text{SNR}_{r}\simeq 34$) that we acquired for the flow through the
converging nozzle using MRV (units in [cm/s]).

###  3.5 Magnetic resonance velocimetry data in a converging nozzle

We now use algorithm 1 to reconstruct and segment the low SNR images
($\boldsymbol{u}^{\star}$) that we acquired during the MRV experiment (section
3.4), and compare them with the high SNR images of the same flow
($\boldsymbol{u}^{\bullet}$ in figure 19). The flow is axisymmetric with zero
swirl. The subscript ‘$x$’ is replaced by ‘$z$’, which denotes the axial
component of velocity, and the subscript ‘$y$’ is replaced by ‘$r$’, which
denotes the radial component of velocity. The low SNR images
($\text{SNR}_{z}=6.7$, $\text{SNR}_{r}=5.8$) required a total scanning time of
5.1 minutes per velocity image (axial and radial components), and the high SNR
images ($\text{SNR}_{z}=44.2$, $\text{SNR}_{r}=34.4$) required a total
scanning time of 137 minutes per velocity image. Since the signal intensity of
an MRV experiment corresponds to the ${}^{1}\mathrm{H}$ spin density, we
segment the spin density image using a thresholding algorithm (Otsu, 1979) in
order to obtain a mask $\psi$, such that $\psi=1$ inside $\Omega$ (the nozzle)
and $\psi=0$ outside $\Omega$. We consider $\psi$ to be the prior information
for the geometry of the nozzle, which also serves as an initial guess for
$\Omega$ $(\Omega_{0})$. For ${\boldsymbol{g}_{i}}_{0}$ we take a parabolic
velocity profile with a peak velocity of $0.6U$, where $U\simeq 5$ cm/s is the
characteristic velocity for this problem. In this case we treat the kinematic
viscosity as an unknown, with a prior distribution
$\mathcal{N}\big{(}\bar{\nu},(0.1\bar{\nu})^{2}\big{)}$, and
$\bar{\nu}=4\times 10^{-6}$m2/s. Note that the axis of the nozzle is not
precisely known beforehand, and since we only solve an axisymmetic
Navier–Stokes problem on the $z-r$ half-plane, we also introduce an unknown
variable for the vertical position of the axis (see appendix C).

|  | image dimension | model dimension | $\sigma_{u_{z}}/U$ | $\sigma_{u_{r}}/U$
---|---|---|---|---|---
nozzle | (3D) | $255\times 128$ |  $300\times 130$ (half-plane) | $1.4168\times 10^{-1}$ | $3.0679\times 10^{-2}$
Regularization | $\sigma_{\phi_{\pm}}/D$ | $\sigma_{\boldsymbol{g}_{i}}/U$ | $\sigma_{\nu}/UD$ | $\Rey_{\phi_{\pm}}$ | $\Rey_{\zeta}$ | $\ell/h$
nozzle | (3D) | 0.25 | 0.5 | $6.2\times 10^{-4}$ | 0.025 | 0.025 | 3
Table 4: Input parameters for the inverse 3D axisymmetric Navier–Stokes
problem.

(a) Low SNR MRV image $u^{\star}_{z}$

(b) Our reconstruction $u_{z}^{\circ}$

(c) Discrepancy
$\sigma_{u_{z}}^{-1}\big{(}u^{\star}_{z}-\mathcal{S}u_{z}^{\circ}\big{)}$

(d) Low SNR MRV image $u^{\star}_{r}$

(e) Our reconstruction $u_{r}^{\circ}$

(f) Discrepancy
$\sigma_{u_{r}}^{-1}\big{(}u^{\star}_{r}-\mathcal{S}u_{r}^{\circ}\big{)}$

Figure 20: Reconstruction (algorithm 1) of low SNR MRV velocity images
depicting the axisymmetric flow (from left to right) in the converging nozzle
(figure 18(b)). Figures 20(a)-20(b) and 20(d)-20(e) show the horizontal,
$u_{x}$, and vertical, $u_{y}$, velocities and share the same colormap
(colorbar not shown). Figures 20(c) and 20(f) show the discrepancy between the
noisy velocity images and the reconstruction (colorbars apply only to figures
20(c) and 20(f)) . The reconstructed flow $\boldsymbol{u}^{\circ}$ is
axisymmetric by construction, therefore $u^{\circ}_{z}$ depicts an even
reflection and $u^{\circ}_{r}$ depicts an odd reflection, so that they can be
compared with the MRV images (see appendix C).

(a) Velocity magnitude $\lvert\boldsymbol{u}^{\circ}\rvert$ and shape
$\partial\Omega^{\circ}$

(b) Reconstruction error history

(c) Low SNR data (grey) and reconstruction

(d) High SNR velocity data

Figure 21: Reconstruction (algorithm 1) of synthetic images depicting the
axisymmetric flow (from left to right) in the converging nozzle (figure
18(b)). Figure 21(a) depicts the reconstructed boundary
$\partial\Omega^{\circ}$ (cyan line), the $2\sigma$ confidence region computed
from the approximated posterior covariance
$\widetilde{\mathcal{C}}_{\mathring{\zeta}}\equiv\widetilde{H}_{\mathring{\zeta}}\mathcal{C}_{\phi_{\pm}}$
(blue region), the ground truth boundary $\partial\Omega^{\bullet}$ (yellow
line), and the initial guess $\partial\Omega_{0}$ (white line). Figure 21(b)
shows the reconstruction error as a function of iteration number. Velocity
slices are drawn for 10 equidistant cross-sections (labelled with the letters
A to J) for both the reconstructed images (figure 21(c)) and the high SNR
images (figure 21(d)), colored red for positive values and blue for negative.

(a) $u_{z}$, TV-B $\lambda/\lambda_{0}=0.1$

(b) $u_{z}$, TV-B $\lambda/\lambda_{0}=0.01$

(c) $u_{z}$, TV-B $\lambda/\lambda_{0}=0.001$

(d) $u_{r}$, TV-B $\lambda/\lambda_{0}=0.1$

(e) $u_{r}$, TV-B $\lambda/\lambda_{0}=0.01$

(f) $u_{r}$, TV-B $\lambda/\lambda_{0}=0.001$

Figure 22: Total variation denoising using Bregman iteration with different
weights $\lambda$ for the low SNR MRV images (figures 20(a) and 20(d))
depicting the axisymmetric flow (from left to right) in the converging nozzle.

(a) Wall shear rates $\gamma^{\circ}_{w}$ and $\gamma^{\bullet}_{w}$

(b) Our pressure reconstruction $p^{\circ}$

Figure 23: (a) Wall shear rates (as for figure 6) and (b) reduced hydrodynamic
pressure inferred from the MRV images depicting the axisymmetric flow in the
converging nozzle.

Using the input parameters of table 4, the algorithm manages to reconstruct
the noisy velocity image and reduce segmentation errors in just 6 iterations,
with total reconstruction error $\mathcal{E}^{\bullet}\simeq 5.94\%$. The
results are presented in figures 20 and 21. We observe that algorithm 1
manages to filter out the noise, the outliers, and the acquisition artefacts
of the low SNR MRV images depicting the axial $u^{\star}_{z}$ (figure 20(a))
and the radial $u^{\star}_{r}$ (figure 20(d)) component of velocity. A notable
difference between these real MRV images and the synthetic MRV images in
sections 3.1 and 3.2, is that the real MRV images display artefacts and
contain outliers. We have not pre-processed the MRV images for example by
removing outliers. The estimated posterior uncertainty of
$\partial\Omega^{\circ}$ is depicted in figure 21(a), in which we observe that
regions with gaps in the data coincide with regions of higher uncertainty.
Although we treat the kinematic viscosity $\nu$ as an unknown parameter, the
posterior distribution of $\nu$ remains effectively unchanged. More precisely,
we infer a kinematic viscosity of ${\nu^{\circ}=3.995\times
10^{-6}\simeq\bar{\nu}}$, with a posterior variance of
$(0.1005\bar{\nu})^{2}$. This is because we use a Bayesian approach to this
inverse problem, where the prior information for $\nu$ is already rich enough.
Technically, the reconstruction functional $\mathscr{E}$ is insensitive to
small changes of $\nu$ (or $1/\Rey$), and, as a result, the prior term in the
gradient of $\nu$ (equation (33)) dominates; i.e. the model $\mathscr{M}$ is
not informative. Physically, it is not possible to infer $\nu$ (with
reasonable certainty) for this particular flow without additional information
on pressure.

As in section 3.3, we compare the denoising performance of algorithm 1 (figure
20) with TV-B (Getreuer, 2012b; van der Walt et al., 2014) (figure 22). We
again observe that algorithm 1 has managed to filter out both the noise and
the artefacts, while the TV-B-denoised images present artefacts, loss of
contrast, and a systematic error that depends on the parameter $\lambda$.

Figure 23(a) shows the reconstructed wall shear rate $\gamma^{\circ}_{w}$,
computed for the reconstructed velocity field $\boldsymbol{u}^{\circ}$ on the
segmented shape $\partial\Omega^{\circ}$, and compares it with the ground
truth wall shear rate $\gamma^{\bullet}_{w}$, computed for the high SNR
velocity field $\boldsymbol{u}^{\bullet}$ (figure 19) on the high SNR shape
$\partial\Omega^{\bullet}$ (1H spin density). We observe that the ground truth
wall shear rate is particularly noisy, as MRV suffers from low resolution and
partial volume effects (Bouillot et al., 2018; Saito et al., 2020) near the
boundaries $\partial\Omega$. Certainly, it is possible to smooth the boundary
$\partial\Omega^{\bullet}$ (which we obtained using the method of Otsu (1979)
for the ${}^{1}\mathrm{H}$ spin density) using conventional image processing
algorithms. However, the velocity field $\boldsymbol{u}^{\bullet}$ will not be
consistent with the new smoothed boundary (the no-slip boundary condition will
not be satisfied). The method that we propose here for the reconstruction and
segmentation of MRV images tackles exactly this problem: it infers the most
likely shape of the boundary ($\partial\Omega^{\circ}$) from the velocity
field itself, without requiring an additional experiment (e.g. CT, MRA) or
manual segmentation using another software. Furthermore, in this Bayesian
setting we can use the ${}^{1}\mathrm{H}$ spin density to introduce a priori
knowledge of $\partial\Omega$ in the form of a prior, which would prove useful
in areas of low velocity magnitudes where the velocity field itself does not
provide enough information in order to segment the boundaries, e.g. flow
within an aneurysm or a heart ventricle (Demirkiran et al., 2021). As a
result, algorithm 1 performs very well in estimating the posterior
distribution of wall shear rate, a quantity which depends both on the velocity
field and the boundary shape, and which is hard to measure otherwise.

###  3.6 Choosing the regularization parameters

Regularization is crucial in order to successfully reconstruct the velocity
field and segment the geometry of the nozzle in the presence of noise,
artefacts, and outliers. Regularization comes from the Navier–Stokes problems
(primal and adjoint) ($\mathscr{M}$), and the regularization of the model
parameters ($\mathscr{R}$).

####  3.6.1 Notes on prior information for the Navier–Stokes unknowns

By adopting a Bayesian inference framework, we assume that the prior
information of an unknown $x$ is a Gaussian random field with mean $\bar{x}$
and covariance $\mathcal{C}_{x}$, i.e.
${x\sim\mathcal{N}(\bar{x},\mathcal{C}_{x})}$ (see section 2.2 and appendix
A). We, therefore, need to provide algorithm 1 with a prior mean and a prior
covariance for every N–S unknown. For the inlet velocity boundary condition,
$\bar{\boldsymbol{g}}_{i}$ can be a smooth approximation to the noisy velocity
data at the inlet, and then $\sigma_{\boldsymbol{g}_{i}}$ is the prior
standard deviation around this mean. For the outlet natural boundary
condition, $\bar{\boldsymbol{g}}_{o}$ can be $0$, and then
$\sigma_{\boldsymbol{g}_{o}}$ determines the confidence of the user regarding
whether or not the outlet is a pseudotraction-free boundary. For both the
inlet and the outlet boundary conditions, the parameter $\ell$, which can be
different for each boundary condition, controls the regularity of the
functions ${\boldsymbol{g}_{i}}$ and ${\boldsymbol{g}_{o}}$, i.e. length
scales smaller than $\ell$ are suppressed. For the shape, $\bar{\phi}_{\pm}$
can be a rough segmentation of the original geometry, and then
$\sigma_{\phi_{\pm}}$ is the prior standard deviation around this mean. For
example, in section 3.3, we set $\sigma_{\phi_{\pm}}$ approximately equal to a
length of 7 pixels by visually inspecting the noisy mask (figure 11(b)). The
same methodology applies to the determination of prior information regarding
the kinematic viscosity $\nu$.

The advantage of this probabilistic framework is that when prior information
is available it can be readily exploited in order to regularize the inverse
problem and facilitate its numerical solution. On the other hand, if there is
no prior information available regarding an unknown, we can assume that this
unknown is distributed according to a zero-mean Gaussian distribution with a
sufficiently large standard deviation $\sigma$.

####  3.6.2 Notes on shape regularization and the choice of
$\Rey_{{\phi_{\pm}}},\Rey_{\zeta}$

For the axisymmetric nozzle (see section 3.5), we avoid overfitting the shape
$\partial\Omega$ by choosing the Reynolds numbers for the geometric flow to be
${\Rey_{{\phi_{\pm}}}=\Rey_{\zeta}=0.025}$. Increasing these Reynolds numbers
to around $1.0$, we start noticing that the assimilated boundary becomes more
susceptible to noise in the image. However, for the simulated aortic aneurysm
(see section 3.3) we chose ${\Rey_{{\phi_{\pm}}}=\Rey_{\zeta}=1.0}$ in order
to preserve high curvature regions. From numerical experiments we have
observed that typical successful values for the Reynolds numbers
$\Rey_{{\phi_{\pm}}},\Rey_{\zeta}$ lie in the interval $(0.01,0.1)$ for low
SNR images ($\text{SNR}<10$) with relatively flat boundaries, in $(0.1,1.0)$
for higher SNR images ($\text{SNR}\geq 10$) with relatively flat boundaries,
and ${\geq 1}$ for geometries with regions of high curvature. Physical
intuition that justifies the use of $\Rey_{{\phi_{\pm}}},\Rey_{\zeta}$ as the
preferred shape regularization parameters is provided in section 2.4.2.

##  4 Conclusions

We have formulated a generalized inverse Navier–Stokes problem for the joint
reconstruction and segmentation of noisy velocity images of steady
incompressible flow. To regularize the inverse problem, we adopt a Bayesian
framework by assuming Gaussian prior distributions for the unknown model
parameters. Although the inverse problem is formulated using variational
methods, every iteration of the nonlinear problem is actually equivalent to a
Gaussian process in Hilbert spaces. We implicitly define the boundaries of the
flow domain in terms of signed distance functions and use Nitsche’s method to
weakly enforce the Dirichlet boundary condition on the moving front. The
moving of the boundaries is expressed by a convection-diffusion equation for
the signed distance function, which allows us to control the regularity of the
boundary by tuning an artificial diffusion coefficient. We use the steepest
ascent directions of the model parameters in conjunction with a quasi-Newton
method (BFGS), and we show how the posterior Gaussian distribution of a model
parameter can be estimated from the reconstructed inverse Hessian.

We devise an algorithm that solves this inverse Navier–Stokes problem and test
it for noisy ($\text{SNR}=2.5,3$) 2D synthetic images of Navier–Stokes flows.
The algorithm successfully reconstructs the velocity images, infers the most
likely boundaries of the flow and estimates their posterior uncertainty. We
then design a magnetic resonance velocimetry (MRV) experiment to obtain images
of a 3D axisymmetric Navier–Stokes flow in a converging nozzle. We acquire MRV
images of poor quality ($\text{SNR}\simeq 6$), intended for
reconstruction/segmentation, and images of higher quality ($\text{SNR}>30$)
that serve as the ground truth. We show that the algorithm performs very well
in reconstructing and segmenting the poor MRV images, which were obtained in
just $10.2$ minutes, and that the reconstruction compares well to the high SNR
images, which required a total acquisition time of $\sim 4.6$ hours. Lastly,
we use the reconstructed images and the segmented (smoothed) domain to
estimate the posterior distribution of the wall shear rate and compare it with
the ground truth. Since the wall shear rate depends on both the shape and the
velocity field, we note that our algorithm provides a consistent treatment to
this problem by jointly reconstructing and segmenting the flow images,
avoiding the design of an additional experiment (e.g. CT, MRA) for the
measurement of the geometry, or the use of external (non physics-informed)
segmentation software.

The present method has several advantages over general image reconstruction
and segmentation algorithms, which do not respect the underlying physics and
the boundary conditions, and, at the same time, provides additional knowledge
of the flow physics (e.g. pressure field and wall shear stress), which is
otherwise difficult to measure. It can be used to substantially decrease
signal acquisition times and provides additional knowledge of the physical
system being imaged. Although our current implementation is restricted to 2D
planar and axisymmetric flows, the method naturally extends to periodic and
unsteady Navier–Stokes problems in complicated 3D geometries.

Declaration of Interests. The authors report no conflict of interest.

##  Appendix A Gaussian measures in Hilbert spaces

The mean of a Gaussian measure $\gamma\sim\mathcal{N}(m,\mathcal{C})$ in
$L^{2}$ is given by

$\displaystyle m\equiv\mathbb{E}h:=\int_{L^{2}}h\ \gamma(dh)\quad.$ (54)

The covariance operator $\mathcal{C}:L^{2}\to L^{2}$ and the covariance
$C:L^{2}\times L^{2}\to\mathbb{R}$ are defined by

$\displaystyle\mathcal{C}x:=\int_{L^{2}}h\big{\langle}x,h\big{\rangle}\
\gamma(dh)\quad,\quad
C(x,x^{\prime}):=\int_{L^{2}}\big{\langle}x,h\big{\rangle}\big{\langle}x^{\prime},h\big{\rangle}\
\gamma(dh)\quad,$ (55)

noting that
$\big{\langle}\mathcal{C}x,x^{\prime}\big{\rangle}=C(x,x^{\prime})$. The above
(Bochner) integrals define integration over the function space $L^{2}$, and
under the measure $\gamma$, and are well defined due to Fernique’s theorem
(Hairer, 2009). These integrals can be directly computed by sampling the
Gaussian measure $\gamma$ with Karhunen–Loève expansion (as in section 2.6).

##  Appendix B Euler–Lagrange system

The integration by parts formulae for the nonlinear term (equation (19)) are

$\displaystyle\int_{\Omega}\big{(}\boldsymbol{u}^{\prime}\boldsymbol{\cdot}\nabla\boldsymbol{u}\big{)}\boldsymbol{\cdot}\boldsymbol{v}$
$\displaystyle=\int_{\Omega}u_{j}^{\prime}\partial_{j}u_{i}\
v_{i}=-\int_{\Omega}u_{j}^{\prime}u_{i}\
\partial_{j}v_{i}+\partial_{j}u_{j}^{\prime}u_{i}\
v_{i}+\int_{\partial\Omega}u^{\prime}_{j}\nu_{j}\ u_{i}v_{i}$
$\displaystyle=-\int_{\Omega}\big{(}\boldsymbol{u}\boldsymbol{\cdot}(\nabla\boldsymbol{v})^{\dagger})\boldsymbol{\cdot}\boldsymbol{u}^{\prime}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\int_{\partial\Omega}(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{v})(\boldsymbol{\nu}\boldsymbol{\cdot}\boldsymbol{u}^{\prime})}\quad,$
(56)
$\displaystyle\int_{\Omega}\big{(}\boldsymbol{u}\boldsymbol{\cdot}\nabla\boldsymbol{u}^{\prime}\big{)}\boldsymbol{\cdot}\boldsymbol{v}$
$\displaystyle=\int_{\Omega}u_{j}\partial_{j}u^{\prime}_{i}\
v_{i}=-\int_{\Omega}\partial_{j}u_{j}u^{\prime}_{i}\
v_{i}+u_{j}u^{\prime}_{i}\
\partial_{j}v_{i}+\int_{\partial\Omega}u_{j}u^{\prime}_{i}\ \nu_{j}v_{i}$
$\displaystyle=-\int_{\Omega}\big{(}\boldsymbol{u}\boldsymbol{\cdot}\nabla\boldsymbol{v}\big{)}\boldsymbol{\cdot}\boldsymbol{u}^{\prime}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\int_{\partial\Omega}(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nu})(\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{u}^{\prime})}\quad.$
(57)

##  Appendix C Axisymmetric inverse Navier–Stokes problem

The axisymmetric Navier–Stokes problem is

$\displaystyle\boldsymbol{u}\boldsymbol{\cdot}\nabla\boldsymbol{u}-\nu\rmDelta\boldsymbol{u}+\nabla
p+\boldsymbol{f}=\boldsymbol{0},\quad\nabla\boldsymbol{\cdot}\boldsymbol{u}=0\quad,$
(58)

where

$\displaystyle\boldsymbol{u}=u_{z}\hat{\boldsymbol{z}}+u_{r}\hat{\boldsymbol{r}}\quad,\quad\nabla\boldsymbol{u}=(\partial_{z}\boldsymbol{u},~{}\partial_{r}\boldsymbol{u})\quad,\quad\rmDelta\boldsymbol{u}=\partial_{z}^{2}u_{z}+\partial_{r}^{2}u_{r}+\frac{1}{r}~{}\partial_{r}u_{r}\quad,$
$\displaystyle\nabla\boldsymbol{\cdot}\boldsymbol{u}=\partial_{z}u_{z}+\partial_{r}u_{r}+\frac{u_{r}}{r}\quad,\quad\boldsymbol{f}=\Big{(}0,\frac{\nu
u_{r}}{r^{2}}\Big{)}\quad,$

and the nonlinear term $\boldsymbol{u}\boldsymbol{\cdot}\nabla\boldsymbol{u}$
retains the same form as in the Cartesian frame.

In order to compare the axisymmetric modeled velocity field with the MRV
images, we introduce two new operators: i) the reflection operator
${\mathcal{R}:\mathbb{R}^{+}\times\mathbb{R}\to\mathbb{R}\times\mathbb{R}}$,
and ii) a rigid transformation
${\mathcal{T}:\mathbb{R}^{2}\to\mathbb{R}^{2}}$. The reconstruction error is
then expressed by

$\mathscr{E}(\boldsymbol{u})\equiv\frac{1}{2}\big{\lVert}\boldsymbol{u}^{\star}-\mathcal{S}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{T}\mathcal{R}}\boldsymbol{u}\big{\rVert}^{2}_{{\mathcal{C}_{\boldsymbol{u}}}}:=\frac{1}{2}\int_{I}\big{(}\boldsymbol{u}^{\star}-\mathcal{S}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{T}\mathcal{R}}\boldsymbol{u}\big{)}{\mathcal{C}^{-1}_{\boldsymbol{u}}}\big{(}\boldsymbol{u}^{\star}-\mathcal{S}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{T}\mathcal{R}}\boldsymbol{u}\big{)}~{}\mathrm{d}x\mathrm{d}y\quad.$
(59)

We introduce an unknown variable for the vertical position of the axisymmetry
axis by letting $\mathcal{T}u=u(x,y+y_{0})$, for $y_{0}=\text{const}$. Then,
the generalized gradient for $y_{0}$ is

$\displaystyle\Big{\langle}D_{y_{0}}\mathscr{J},y^{\prime}_{0}\Big{\rangle}_{\mathbb{R}}=\Big{\langle}-\int_{I}{\mathcal{C}^{-1}_{\boldsymbol{u}}}\big{(}\boldsymbol{u}^{\star}-\mathcal{S}\mathcal{T}\mathcal{R}\boldsymbol{u}\big{)}\big{(}\mathcal{S}\mathcal{T}\mathcal{R}~{}\partial_{y}\boldsymbol{u}\big{)},~{}y^{\prime}_{0}\Big{\rangle}_{\mathbb{R}}\quad,$
(60)

and $y_{0}$ is treated in the same way as the inverse Navier–Stokes problem
unknowns $\boldsymbol{x}$.

## References

  * Benning & Burger (2018) Benning, Martin & Burger, Martin 2018 Modern regularization methods for inverse problems. Acta Numerica 27, 1–111, arXiv: 1801.09922.
  * Benning et al. (2014) Benning, Martin, Gladden, Lynn, Holland, Daniel, Schönlieb, Carola Bibiane & Valkonen, Tuomo 2014 Phase reconstruction from velocity-encoded MRI measurements - A survey of sparsity-promoting variational approaches. Journal of Magnetic Resonance 238, 26–43.
  * Benzi et al. (2005) Benzi, Michele, Golubt, Gene H. & Liesen, Jörg 2005 Numerical solution of saddle point problems. Acta Numerica 14, 1–137.
  * Benzi & Olshanskii (2006) Benzi, Michele & Olshanskii, Maxim A. 2006 An augmented Lagrangian-based approach to the Oseen problem. SIAM Journal on Scientific Computing 28 (6), 2095–2113.
  * Bogachev (1998) Bogachev, Vladimir I. 1998 Gaussian measures. American Mathematical Society.
  * Bouillot et al. (2018) Bouillot, Pierre, Delattre, Bénédicte M.A., Brina, Olivier, Ouared, Rafik, Farhat, Mohamed, Chnafa, Christophe, Steinman, David A., Lovblad, Karl Olof, Pereira, Vitor M. & Vargas, Maria I. 2018 3D phase contrast MRI: Partial volume correction for robust blood flow quantification in small intracranial vessels. Magnetic Resonance in Medicine 79 (1), 129–140.
  * Brenner, Susanne, Scott (2008) Brenner, Susanne, Scott, Ridgway 2008 The Mathematical theory of finite element methods, , vol. 15. Texts in Applied Mathematics, Springer.
  * Burger (2001) Burger, Martin 2001 A level set method for inverse problems. Inverse Problems 17 (5), 1327–1355.
  * Burger (2003) Burger, Martin 2003 A framework for the construction of level set methods for shape optimization and reconstruction. Interfaces and Free Boundaries 5 (3), 301–329.
  * Burger & Osher (2005) Burger, Martin & Osher, Stanley J. 2005 A Survey in Mathematics for Industry: A survey on level set methods for inverse problems and optimal design. European Journal of Applied Mathematics 16 (2), 263–301.
  * Burman (2010) Burman, Erik 2010 La pénalisation fantôme. Comptes Rendus Mathematique 348 (21-22), 1217–1220.
  * Burman et al. (2015) Burman, Erik, Claus, Susanne & Massing, André 2015 A stabilized cut finite element method for the three field stokes problem. SIAM Journal on Scientific Computing 37 (4), A1705–A1726.
  * Burman & Hansbo (2012) Burman, Erik & Hansbo, Peter 2012 Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche method. Applied Numerical Mathematics 62 (4), 328–341.
  * Chan & Vese (2001) Chan, Tony F. & Vese, Luminita A. 2001 Active contours without edges. IEEE Transactions on Image Processing 10 (2), 266–277.
  * Cheng (2008) Cheng, Nian-Sheng 2008 Formula for the Viscosity of a Glycerol-Water Mixture. Industrial & Engineering Chemistry Research 47 (9), 3285–3288.
  * Codina (2002) Codina, Ramon 2002 Stabilized finite element approximation of transient incompressible flows using orthogonal subscales. Computer Methods in Applied Mechanics and Engineering 191 (39-40), 4295–4321.
  * Corona et al. (2021) Corona, Veronica, Benning, Martin, Gladden, Lynn F., Reci, Andi, Sederman, Andrew J. & Schönlieb, Carola-Bibiane 2021 Joint Phase Reconstruction and Magnitude Segmentation from Velocity-Encoded MRI Data, pp. 1–24. Cham: Springer International Publishing.
  * Cotter et al. (2009) Cotter, S. L., Dashti, M., Robinson, J. C. & Stuart, A. M. 2009 Bayesian inverse problems for functions and applications to fluid mechanics. Inverse Problems 25 (11).
  * Crane et al. (2017) Crane, Keenan, Weischedel, Clarisse & Wardetzky, Max 2017 The heat method for distance computation. Communications of the ACM 60 (11), 90–99.
  * Demirkiran et al. (2021) Demirkiran, Ahmet, van Ooij, Pim, Westenberg, Jos J M, Hofman, Mark B M, van Assen, Hans C, Schoonmade, Linda J, Asim, Usman, Blanken, Carmen P S, Nederveen, Aart J, van Rossum, Albert C & Götte, Marco J W 2021 Clinical intra-cardiac 4D flow CMR: acquisition, analysis, and clinical applications. European Heart Journal - Cardiovascular Imaging 23 (2), 154–165, arXiv: https://academic.oup.com/ehjcimaging/article-pdf/23/2/154/42828964/jeab112.pdf.
  * Donoho (2006) Donoho, David L. 2006 Compressed sensing. IEEE Transactions on Information Theory 52 (4), 1289–1306.
  * Edelstein et al. (1980) Edelstein, W A, Hutchison, J M S, Johnson, G & Redpath, T 1980 Spin warp NMR imaging and applications to human whole-body imaging. Physics in Medicine and Biology 25 (4), 751–756.
  * Elkins & Alley (2007) Elkins, Christopher J & Alley, Marcus T 2007 Magnetic resonance velocimetry: applications of magnetic resonance imaging in the measurement of fluid motion. Experiments in Fluids 43 (6), 823–858.
  * Evans (2010) Evans, Lawrence C. 2010 Partial Differential Equations, 2nd edn. American Mathematical Society.
  * Fletcher (2000) Fletcher, R. 2000 Practical Methods of Optimization. John Wiley & Sons.
  * Fukushima (1999) Fukushima, Eiichi 1999 Nuclear magnetic resonance as a tool to study flow. Annual Review of Fluid Mechanics 31, 95–123.
  * Funke et al. (2019) Funke, Simon Wolfgang, Nordaas, Magne, Evju, Øyvind, Alnæs, Martin Sandve & Mardal, Kent Andre 2019 Variational data assimilation for transient blood flow simulations: Cerebral aneurysms as an illustrative example. International Journal for Numerical Methods in Biomedical Engineering 35 (1), 1–27.
  * Getreuer (2012a) Getreuer, Pascal 2012a Chan-Vese Segmentation. Image Processing On Line 2, 214–224.
  * Getreuer (2012b) Getreuer, Pascal 2012b Rudin–Osher–Fatemi Total Variation Denoising using Split Bregman. Image Processing On Line 2 (1), 74–95.
  * Gillissen et al. (2018) Gillissen, Jurriaan J.J., Vilquin, Alexandre, Kellay, Hamid, Bouffanais, Roland & Yue, Dick K.P. 2018 A space-time integral minimisation method for the reconstruction of velocity fields from measured scalar fields. Journal of Fluid Mechanics 854, 348–366.
  * Gillissen et al. (2019) Gillissen, Jurriaan J. J., Bouffanais, Roland & Yue, Dick K. P. 2019 Data assimilation method to de-noise and de-filter particle image velocimetry data. Journal of Fluid Mechanics 877, 196–213.
  * Gudbjartsson & Patz (1995) Gudbjartsson, HáKon & Patz, Samuel 1995 The rician distribution of noisy mri data. Magnetic Resonance in Medicine 34 (6), 910–914.
  * Hairer (2009) Hairer, Martin 2009 An Introduction to Stochastic PDEs. Lecture Notes., arXiv: 0907.4178.
  * Harris et al. (2020) Harris, Charles R., Millman, K. Jarrod, van der Walt, Stéfan J., Gommers, Ralf, Virtanen, Pauli, Cournapeau, David, Wieser, Eric, Taylor, Julian, Berg, Sebastian, Smith, Nathaniel J., Kern, Robert, Picus, Matti, Hoyer, Stephan, van Kerkwijk, Marten H., Brett, Matthew, Haldane, Allan, del Río, Jaime Fernández, Wiebe, Mark, Peterson, Pearu, Gérard-Marchant, Pierre, Sheppard, Kevin, Reddy, Tyler, Weckesser, Warren, Abbasi, Hameer, Gohlke, Christoph & Oliphant, Travis E. 2020 Array programming with NumPy. Nature 585 (7825), 357–362.
  * Heister & Rapin (2013) Heister, Timo & Rapin, Gerd 2013 Efficient augmented Lagrangian-type preconditioning for the Oseen problem using Grad-Div stabilization. International Journal for Numerical Methods in Fluids 71 (1), 118–134.
  * Hoang et al. (2014) Hoang, Viet Ha, Law, Kody J.H. & Stuart, Andrew M. 2014 Determining white noise forcing from Eulerian observations in the Navier–Stokes equation. Stochastics and Partial Differential Equations: Analysis and Computations 2 (2), 233–261.
  * Katritsis et al. (2007) Katritsis, Demosthenes, Kaiktsis, Lambros, Chaniotis, Andreas, Pantos, John, Efstathopoulos, Efstathios P. & Marmarelis, Vasilios 2007 Wall Shear Stress: Theoretical Considerations and Methods of Measurement. Progress in Cardiovascular Diseases 49 (5), 307–329.
  * Koltukluoğlu (2019) Koltukluoğlu, Taha Sabri 2019 Fourier spectral dynamic data assimilation: Interlacing cfd with 4d flow mri. In Medical Image Computing and Computer Assisted Intervention – MICCAI 2019 (ed. Dinggang Shen, Tianming Liu, Terry M. Peters, Lawrence H. Staib, Caroline Essert, Sean Zhou, Pew-Thian Yap & Ali Khan), pp. 741–749. Cham: Springer International Publishing.
  * Koltukluoğlu & Blanco (2018) Koltukluoğlu, Taha S. & Blanco, Pablo J. 2018 Boundary control in computational haemodynamics. Journal of Fluid Mechanics 847, 329–364.
  * Koltukluoğlu et al. (2019) Koltukluoğlu, Taha S., Cvijetić, Gregor & Hiptmair, Ralf 2019 Harmonic balance techniques in cardiovascular fluid mechanics. In Medical Image Computing and Computer Assisted Intervention – MICCAI 2019, pp. 486–494. Cham: Springer International Publishing.
  * Lam et al. (2015) Lam, Siu Kwan, Pitrou, Antoine & Seibert, Stanley 2015 Numba: a LLVM-based Python JIT compiler. In Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC, pp. 1–6.
  * Lustig et al. (2007) Lustig, Michael, Donoho, David & Pauly, John M. 2007 Sparse MRI: The application of compressed sensing for rapid MR imaging. Magnetic Resonance in Medicine 58 (6), 1182–1195.
  * Mantle & Sederman (2003) Mantle, M D & Sederman, A J 2003 Dynamic MRI in chemical process and reaction engineering. Progress in Nuclear Magnetic Resonance Spectroscopy 43 (1), 3–60.
  * Markl et al. (2012) Markl, Michael, Frydrychowicz, Alex, Kozerke, Sebastian, Hope, Mike & Wieben, Oliver 2012 4D flow MRI. Journal of Magnetic Resonance Imaging 36 (5), 1015–1036.
  * Massing et al. (2013) Massing, André, Larson, Mats G. & Logg, Anders 2013 Efficient implementation of finite element methods on nonmatching and overlapping meshes in three dimensions. SIAM Journal on Scientific Computing 35 (1), C23–C47, arXiv: https://doi.org/10.1137/11085949X.
  * Massing et al. (2014) Massing, André, Larson, Mats G., Logg, Anders & Rognes, Marie E. 2014 A Stabilized Nitsche Fictitious Domain Method for the Stokes Problem. Journal of Scientific Computing 61 (3), 604–628, arXiv: 1206.1933.
  * Massing et al. (2018) Massing, A., Schott, B. & Wall, W. A. 2018 A stabilized Nitsche cut finite element method for the Oseen problem. Computer Methods in Applied Mechanics and Engineering 328, 262–300, arXiv: 1611.02895.
  * Mirtich (1996) Mirtich, Brian 1996 Fast and Accurate Computation of Polyhedral Mass Properties. Journal of Graphics Tools 1 (2), 31–50.
  * Mons et al. (2017) Mons, Vincent, Chassaing, Jean Camille & Sagaut, Pierre 2017 Optimal sensor placement for variational data assimilation of unsteady flows past a rotationally oscillating cylinder. Journal of Fluid Mechanics 823, 230–277.
  * Morris et al. (2016) Morris, Paul D., Narracott, Andrew, Von Tengg-Kobligk, Hendrik, Soto, Daniel Alejandro Silva, Hsiao, Sarah, Lungu, Angela, Evans, Paul, Bressloff, Neil W., Lawford, Patricia V., Rodney Hose, D. & Gunn, Julian P. 2016 Computational fluid dynamics modelling in cardiovascular medicine. Heart 102 (1), 18–28.
  * Nitsche (1971) Nitsche, J. 1971 Ü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 36 (1), 9–15.
  * Nocedal & Wright (2006) Nocedal, Jorge & Wright, Stephen J. 2006 Numerical Optimization, 2nd edn. Springer New York, NY.
  * Okuta et al. (2017) Okuta, Ryosuke, Unno, Yuya, Nishino, Daisuke, Hido, Shohei & Loomis, Crissman 2017 CuPy: A NumPy-Compatible Library for NVIDIA GPU Calculations. In Proceedings of Workshop on Machine Learning Systems (LearningSys) in The Thirty-first Annual Conference on Neural Information Processing Systems (NIPS).
  * Osher & Sethian (1988) Osher, Stanley & Sethian, James A. 1988 Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics 79 (1), 12–49.
  * Otsu (1979) Otsu, Nobuyuki 1979 A threshold selection method from gray-level histograms. IEEE Transactions on Systems, Man, and Cybernetics 9 (1), 62–66.
  * Peper et al. (2019) Peper, Eva S, Gottwald, Lukas M, Zhang, Quinwei, Coolen, Bram F, van Ooij, Pim, Nederveen, Aart J & Strijkers, Gustav J 2019 Highly accelerated carotid 4D flow MRI using Pseudo-Spiral Cartesian acquisition and a Total Variation constrained Compressed Sensing reconstruction. Journal of Cardiovascular Magnetic Resonance p. IN PRESS.
  * Saito et al. (2020) Saito, Kozue, Abe, Soichiro, Kumamoto, Masaya, Uchihara, Yuto, Tanaka, Akito, Sugie, Kazuma, Ihara, Masafumi, Koga, Masatoshi & Yamagami, Hiroshi 2020 Blood Flow Visualization and Wall Shear Stress Measurement of Carotid Arteries Using Vascular Vector Flow Mapping. Ultrasound in Medicine and Biology 46 (10), 2692–2699.
  * Sankaran et al. (2016) Sankaran, Sethuraman, Kim, Hyun Jin, Choi, Gilwoo & Taylor, Charles A. 2016 Uncertainty quantification in coronary blood flow simulations: Impact of geometry, boundary conditions and blood viscosity. Journal of Biomechanics 49 (12), 2540–2547.
  * Schott & Wall (2014) Schott, B. & Wall, W. A. 2014 A new face-oriented stabilized XFEM approach for 2D and 3D incompressible Navier-Stokes equations. Computer Methods in Applied Mechanics and Engineering 276, 233–265.
  * Sethian (1996) Sethian, J. A. 1996 A fast marching level set method for monotonically advancing fronts. Proceedings of the National Academy of Sciences of the United States of America 93 (4), 1591–1595.
  * Sharma et al. (2019) Sharma, Arjun, Rypina, Irina I., Musgrave, Ruth & Haller, George 2019 Analytic reconstruction of a two-dimensional velocity field from an observed diffusive scalar. Journal of Fluid Mechanics 871, 755–774, arXiv: 1904.04919.
  * Sotelo et al. (2016) Sotelo, Julio, Urbina, Jesus, Valverde, Israel, Tejos, Cristian, Irarrazaval, Pablo, Andia, Marcelo E., Uribe, Sergio & Hurtado, Daniel E. 2016 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 35 (6), 1475–1487.
  * Stejskal & Tanner (1965) Stejskal, E. O. & Tanner, J. E. 1965 Spin diffusion measurements: Spin echoes in the presence of a time-dependent field gradient. The Journal of Chemical Physics 42 (1), 288–292.
  * Stuart (2010) Stuart, A. M. 2010 Inverse problems: A Bayesian perspective. Acta Numerica 19 (2010), 451–459.
  * Tarantola (2005) Tarantola, Albert 2005 Inverse Problem Theory and Methods for Model Parameter Estimation. SIAM.
  * Tezduyar (1991) Tezduyar, T.E. 1991 Stabilized finite element formulations for incompressible flow computations. Advances in Applied Mechanics, vol. 28, pp. 1–44. Elsevier.
  * Varadhan (1967a) Varadhan, S. R.S. 1967a Diffusion processes in a small time interval. Communications on Pure and Applied Mathematics 20 (4), 659–685.
  * Varadhan (1967b) Varadhan, S. R.S. 1967b On the behavior of the fundamental solution of the heat equation with variable coefficients. Communications on Pure and Applied Mathematics 20 (2), 431–455.
  * Verma et al. (2020) Verma, Siddhartha, Papadimitriou, Costas, Lüthen, Nora, Arampatzis, Georgios & Koumoutsakos, Petros 2020 Optimal sensor placement for artificial swimmers. Journal of Fluid Mechanics 884, A24.
  * Virtanen et al. (2020) Virtanen, Pauli & others 2020 SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods 17, 261–272.
  * Volk & Kähler (2018) Volk, Andreas & Kähler, Christian J 2018 Density model for aqueous glycerol solutions. Experiments in Fluids 59 (5), 75.
  * Walker (2015) Walker, Shawn W. 2015 The Shapes of Things: A Practical Guide to Differential Geometry and the Shape Derivative. Advances in Design and Control, SIAM, Philadelphia, PA.
  * van der Walt et al. (2014) van der Walt, Stéfan, Schönberger, Johannes L., Nunez-Iglesias, Juan, Boulogne, François, Warner, Joshua D., Yager, Neil, Gouillart, Emmanuelle, Yu, Tony & the scikit-image contributors 2014 scikit-image: image processing in Python. PeerJ 2, e453.
  * Wapler et al. (2014) Wapler, Matthias C, Leupold, Jochen, Dragonu, Iulius, von Elverfeld, Dominik, Zaitsev, Maxim & Wallrabe, Ulrike 2014 Magnetic properties of materials for MR engineering, micro-MR and beyond. Journal of Magnetic Resonance 242, 233–242.
  * Yu et al. (2019) Yu, Hans, Juniper, Matthew P. & Magri, Luca 2019 Combined state and parameter estimation in level-set methods. Journal of Computational Physics 399, 108950, arXiv: 1903.00321.