# Hidden Fluid Mechanics: A Navier-Stokes Informed Deep Learning Framework for
Assimilating Flow Visualization Data

Maziar Raissi, Alireza Yazdani, and George Em Karniadakis  Division of Applied
Mathematics, Brown University,
Providence, RI, 02912, USA

###### Abstract

We present _hidden fluid mechanics_ (HFM), a physics informed deep learning
framework capable of encoding an important class of physical laws governing
fluid motions, namely the Navier-Stokes equations. In particular, we seek to
leverage the underlying conservation laws (i.e., for mass, momentum, and
energy) to infer hidden quantities of interest such as velocity and pressure
fields merely from spatio-temporal visualizations of a passive scaler (e.g.,
dye or smoke), transported in arbitrarily complex domains (e.g., in human
arteries or brain aneurysms). Our approach towards solving the aforementioned
data assimilation problem is unique as we design an algorithm that is agnostic
to the geometry or the initial and boundary conditions. This makes HFM highly
flexible in choosing the spatio-temporal domain of interest for data
acquisition as well as subsequent training and predictions. Consequently, the
predictions made by HFM are among those cases where a pure machine learning
strategy or a mere scientific computing approach simply cannot reproduce. The
proposed algorithm achieves accurate predictions of the pressure and velocity
fields in both two and three dimensional flows for several benchmark problems
motivated by real-world applications. Our results demonstrate that this
relatively simple methodology can be used in physical and biomedical problems
to extract valuable quantitative information (e.g., lift and drag forces or
wall shear stresses in arteries) for which direct measurements may not be
possible.

###### keywords:

data-driven scientific computing , partial differential equations , physics
informed machine learning , inverse problems , data assimilation

††journal: Journal Name

##  1 Introduction

Recent advances in deep learning and computational resources in addition to
new data recordings and sensor technologies have the potential to
revolutionize our understanding of the physical world in modern application
areas such as image recognition [krizhevsky2012imagenet], drug discovery
[gawehn2016deep], and bioinformatics [alipanahi2015predicting]. Moreover, many
concepts from deep learning can be integrated with classical methods in
mathematical physics to help us solve sophisticated data assimilation and
inverse problems involving complex dynamical systems. This integration of
nonlinear dynamics and deep learning opens the door for principled methods for
model construction and predictive modeling. Furthermore, recent developments
in _physics-informed deep learning_ have enabled us to train deep neural
networks with millions of parameters using very few training examples, for
cases where the available data are known to respect a given physical law
described by a system of partial differential equations [raissi2017physics_I,
raissi2017physics_II]. The prior knowledge of an underlying physical law
introduces important structure that effectively regularizes the minimization
procedure in the training of neural networks, and enables them to generalize
well even when only a few training examples are available. An important class
of these physical laws govern the fluid motions in nature, which will be the
focus of the present work.

Quantifying the flow dynamics essentially requires the knowledge of spatio-
temporal fluid velocity and pressure fields, and has been the centerpiece of
experimental and theoretical fluid mechanics for centuries. Traditionally,
this has been achieved by measuring the instantaneous fluid velocity and
pressure experimentally, or by solving the governing equations of fluid
mechanics in a precisely defined geometry with proper initial and boundary
conditions. Despite significant advances in experimental fluid mechanics
(e.g., with the emergence of particle image velocimetry), the use of
measurements to reliably infer fluid velocity and pressure/stress fields is
not a straightforward task. Furthermore, although experimental measurements of
external flows (e.g., flow past a bluff object) are obtained relatively
easily, albeit in small subdomains, leading to two or even three dimensional
vector fields for velocity, the quantification of velocity fields for internal
flows (e.g., blood flow in the vascular networks) could become very difficult
or impractical. From the theoretical standpoint, the governing equations of
fluid mechanics have been derived from conservation laws (i.e., conservation
of mass, momentum and energy) leading to partial differential equations such
as the well-known Navier-Stokes equations for Newtonian fluids
[batchelor2000introduction]. With the increasing computational power, accurate
solutions of such equations in a “forward” setting are now available with
lower cost even in the turbulence regime using Direct Numerical Simulations
(DNS), hence resolving all energetically important scales down to Kolmogorev
dissipation scales. However, the assimilation of experimental fluid mechanics
data into the mathematical models to infer velocity and pressure fields in an
“inverse” setting has not been adequately addressed so far. The main objective
of this article is to extend and incorporate the aforementioned physics-
informed deep learning algorithms [raissi2017physics_I, raissi2017physics_II]
to leverage the hidden physics of fluid mechanics to infer the latent
quantities of interest (e.g., the velocity and pressure fields) through
minimalistic usage of data obtained by the “passive scalar” transport (e.g.,
transport of dye or smoke).

Let us consider the transport of a passive scalar field $c$ by a velocity
field $\bm{u}$. Such a problem arises, for example, while studying the
spreading of smoke or dye advected by a given velocity field and subject to
molecular diffusion. We expect the dynamics of $c$ to be governed by an
effective transport equation [vergassola1997scalar] written in the form of

$\rho\left(\partial_{t}c+\bm{u}\cdot\nabla c\right)=\kappa\nabla^{2}c,$ (1)

where $\kappa$ is the molecular diffusivity and $\rho$ is the density.
Transport of scalar fields in fluid flow has been studied in numerous
applications such as aerodynamics, biofluid mechanics, and non-reactive flow
mixing to name a few. The use of smoke in wind tunnels or dye in water tunnels
for flow visualization and quantification has long been practiced in
experimental fluid mechanics [barlow2015low]. Moreover, recent techniques in
planar laser induced fluorescence imaging combined with particle image
velocimetry have been developed to assess the relationships between scalar and
velocity/vorticity fields [koochesfahani1985laser, crimaldi2008planar]. The
use of scalar transport in conjunction with advanced imaging modalities to
quantify blood flow in the vascular networks is now a common practice. For
example, coronary computed tomography (CT) angiography is typically performed
on multidetector CT systems after the injection of non-diffusible iodine
contrast agent, which allows coronary artery visualization and the detection
of coronary stenoses [voros2011coronary]. Another example is the
quantification of cerebral blood flow, which is detrimental in the prognostic
assessments in stroke patients using a contrast agent and perfusion CT
[wintermark2002prognostic], and in cognitive neuroscience with the use of
functional magnetic resonance imaging that only relies on the blood-oxygen-
level dependent contrast to measure brain activity [aron2007triangulating].

Inspired by recent developments in _physics-informed deep learning_
[raissi2017physics_I, raissi2017physics_II] and _deep hidden physics models_
[raissi2018deep], we propose to leverage the hidden physics of fluid mechanics
(i.e., the Navier-Stokes equations) and infer the latent quantities of
interest (e.g., the velocity and pressure fields) by approximating them using
deep neural networks. This choice is motivated by modern techniques for
solving forward and inverse problems associated with partial differential
equations, where the unknown solution is approximated either by a neural
network [raissi2017physics_I, raissi2017physics_II, raissi2018forwardbackward,
raissi2018multistep] or a Gaussian process [raissi2018numerical,
raissi2018hidden, raissi2017inferring, raissi2017machine,
raissi2017parametric, perdikaris2017nonlinear, raissi2016deep]. Moreover,
placing a prior on the solution itself is fully justified by the similar
approach pursued in the past century using classical methods of solving
partial differential equations such as finite elements and spectral methods,
where one would expand the unknown solution in terms of an appropriate set of
basis functions. Our focus here is the transport of a passive scalar by
incompressible Newtonian flows in both unbounded geometries (i.e., external
flows) and confined internal flows. We demonstrate the success of our Navier-
Stokes informed deep learning algorithm by recovering the flow velocity and
pressure fields solely from time series data collected on the passive scalar
in arbitrary domains.

##  2 Problem Setup and Solution Methodology

Starting with the transport equation (1), we can rewrite it in the following
non-dimensional form

$c_{t}+uc_{x}+vc_{y}+wc_{z}=\text{\footnotesize
Pec}^{-1}(c_{xx}+c_{yy}+c_{zz}),$ (2)

which governs the evolution of the normalized concentration $c(t,x,y,z)$ of a
passive scalar transported by an incompressible Newtonian fluid whose dynamics
are described by the Navier-Stokes and continuity equations (also non-
dimensional) given below

$\begin{array}[]{l}u_{t}+uu_{x}+vu_{y}+wu_{z}=-p_{x}+\text{\footnotesize
Re}^{-1}(u_{xx}+u_{yy}+u_{zz}),\\\
v_{t}+uv_{x}+vv_{y}+wv_{z}=-p_{y}+\text{\footnotesize
Re}^{-1}(v_{xx}+v_{yy}+v_{zz}),\\\
w_{t}+uw_{x}+vw_{y}+ww_{z}=-p_{z}+\text{\footnotesize
Re}^{-1}(w_{xx}+w_{yy}+w_{zz}),\\\ u_{x}+v_{y}+w_{z}=0.\end{array}$ (3)

In equations (2) and (3), $u$, $v$, and $w$ are the $x$, $y$, and $z$
components of the velocity field $\bm{u}=(u,v,w)$, respectively, and $p$
denotes the pressure. A passive scalar is a diffusive field in the fluid flow
that has no dynamical effect (such as the effect of temperature in a buoyancy-
driven flow) on the fluid motion itself. Smoke and dye are two typical
examples of passive scalars. Here, Re represents the Reynolds number, a
dimensionless quantity in fluid mechanics, and Pec denotes the Péclet number,
a dimensionless quantity relevant in the study of transport phenomena in
flowing fluids. In this work, we assume that the only observables are noisy
data $\\{t^{n},x^{n},y^{n},z^{n},c^{n}\\}_{n=1}^{N}$ on the concentration
$c(t,x,y,z)$ of the passive scalar. Given such data, scattered in space and
time, we are interested in inferring the latent (hidden) quantities
$u(t,x,y,z)$, $v(t,x,y,z)$, $w(t,x,y,z)$, and $p(t,x,y,z)$. To solve the
aforementioned data assimilation problem, we would like to design an algorithm
that is agnostic to the geometry as well as the initial and boundary
conditions. This is enabling as it will give us the flexibility to work with
data acquired in arbitrarily complex domains such as human arteries or brain
aneurysms.

One question that would naturally arise is whether the information on the
passive scalar in the training domain and near its boundaries is sufficient to
result in a unique velocity field. The answer is that normally there are no
guarantees for unique solutions unless some form of boundary conditions are
explicitly imposed on the domain boundaries. However, as shown later for the
benchmark problems studied in the current work, an informed selection of the
training boundaries in the regions where there are sufficient gradients in the
concentration of the passive scalar would eliminate the requirement of
imposing velocity and pressure boundary conditions. In addition to the proper
design of the training domain and the use of the passive scalar $c$, we can
improve the model predictions further by introducing an auxiliary variable
$d:=1-c$ (essentially the complement of $c$) that satisfies the transport
equation

$d_{t}+ud_{x}+vd_{y}+wd_{z}=\text{\footnotesize
Pec}^{-1}(d_{xx}+d_{yy}+d_{zz}).$ (4)

This can be clearly seen by a change of variable from $c$ to $d=1-c$ and using
equation (2). The complementary nature of the auxiliary variable $d$ helps the
algorithm better detect the geometry and the corresponding boundary
conditions. This makes the training algorithm agnostic to the physical
geometry of the problem, hence, keeping the implementation and training
significantly simpler. In addition, the region of interest becomes very
flexible to be chosen with its boundaries no longer required to be the
physical boundaries. Another advantage of the use of $d$ is the improvement in
the accuracy of the algorithm predictions with no additional cost.

Figure 1: _Navier-Stokes informed neural networks:_ A plain vanilla densely
connected (physics uninformed) neural network, with 10 hidden layers and 50
neurons per hidden layer per output variable (i.e., $6\times 50=300$ neurons
per hidden layer), takes the input variables $t,x,y,z$ and outputs
$c,d,u,v,w$, and $p$. As for the activation functions, we use
$\sigma(x)=\sin(x)$. For illustration purposes only, the network depicted in
this figure comprises of 2 hidden layers and 7 neurons per hidden layers. We
employ automatic differentiation to obtain the required derivatives to compute
the residual (physics informed) networks $e_{1}$, $e_{2}$, $e_{3}$, $e_{4}$,
$e_{5}$, and $e_{6}$. If a term does not appear in the blue boxes (e.g.,
$u_{xy}$ or $u_{tt}$), its coefficient is assumed to be zero. It is worth
emphasizing that unless the coefficient in front of a term is non-zero, that
term is not going to appear in the actual “compiled” computational graph and
is not going to contribute to the computational cost of a feed forward
evaluation of the resulting network. The total loss function is composed of
the regression loss of the passive scalar $c$ and its complement $d$ on the
training data, and the loss imposed by the differential equations
$e_{1}-e_{6}$. Here, $I$ denotes the identity operator and the differential
operators $\partial_{t},\partial_{x},\partial_{y}$, and $\partial_{z}$ are
computed using automatic differentiation and can be thought of as “activation
operators”. Moreover, the gradients of the loss function are backpropogated
through the entire network to train the parameters using the Adam optimizer.

Following our earlier work on physics-informed deep learning
[raissi2017physics_I, raissi2017physics_II] and _deep hidden physics models_
[raissi2018deep], we approximate the function

$(t,x,y,z)\longmapsto(c,d,u,v,w,p)$

by a deep neural network and obtain the following _Navier-Stokes informed
neural networks_ (see figure 1) corresponding to equations (2), (3), and (4);
i.e.,

$\begin{array}[]{l}e_{1}:=c_{t}+uc_{x}+vc_{y}+wc_{z}-\text{\footnotesize
Pec}^{-1}(c_{xx}+c_{yy}+c_{zz}),\\\
e_{2}:=d_{t}+ud_{x}+vd_{y}+wd_{z}-\text{\footnotesize
Pec}^{-1}(d_{xx}+d_{yy}+d_{zz}),\\\
e_{3}:=u_{t}+uu_{x}+vu_{y}+wu_{z}+p_{x}-\text{\footnotesize
Re}^{-1}(u_{xx}+u_{yy}+u_{zz}),\\\
e_{4}:=v_{t}+uv_{x}+vv_{y}+wv_{z}+p_{y}-\text{\footnotesize
Re}^{-1}(v_{xx}+v_{yy}+v_{zz}),\\\
e_{5}:=w_{t}+uw_{x}+vw_{y}+ww_{z}+p_{z}-\text{\footnotesize
Re}^{-1}(w_{xx}+w_{yy}+w_{zz}),\\\ e_{6}:=u_{x}+v_{y}+w_{z}.\end{array}$ (5)

We acquire the required derivatives to compute the residual networks $e_{1}$,
$e_{2}$, $e_{3}$, $e_{4}$, $e_{5}$, and $e_{6}$ by applying the chain rule for
differentiating compositions of functions using automatic differentiation
[baydin2015automatic]. A schematic representation of the resulting _Navier-
Stokes informed neural networks_ is given in figure 1. It is worth emphasizing
that automatic differentiation is different from, and in several respects
superior to, numerical or symbolic differentiation – two commonly encountered
techniques of computing derivatives. In its most basic description
[baydin2015automatic], automatic differentiation relies on the fact that all
numerical computations are ultimately compositions of a finite set of
elementary operations for which derivatives are known. Combining the
derivatives of the constituent operations through the chain rule gives the
derivative of the overall composition. This allows accurate evaluation of
derivatives at machine precision with ideal asymptotic efficiency and only a
small constant factor of overhead. In particular, to compute the required
derivatives we rely on Tensorflow [abadi2016tensorflow], which is a popular
and relatively well documented open-source software library for automatic
differentiation and deep learning computations. In TensorFlow, before a model
is run, its computational graph is defined statically rather than dynamically
as for instance in PyTorch [paszke2017automatic]. This is an important feature
as it allows us to create and compile the computational graph for the Navier-
Stokes informed neural networks (5) only once and keep it fixed throughout the
training procedure. This leads to significant reduction in the computational
cost of the proposed framework.

The shared parameters of the neural networks for $c$, $d$, $u$, $v$, $w$, and
$p$ can be learned by minimizing the following sum of squared errors loss
function

$\displaystyle SSE$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{N}|c(t^{n},x^{n},y^{n},z^{n})-c^{n}|^{2}$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{N}|d(t^{n},x^{n},y^{n},z^{n})-d^{n}|^{2}$
$\displaystyle+$
$\displaystyle\sum_{i=1}^{6}\sum_{n=1}^{N}|e_{i}(t^{n},x^{n},y^{n},z^{n})|^{2}.$

Here, $d^{n}=1-c^{n}$ and the first two terms correspond to the training data
on the concentration $c(t,x,y,z)$ of the passive scalar while the last term
enforces the structure imposed by equations (2), (3), and (4) at a finite set
of measurement points whose number and locations are taken to be the same as
the training data. However, it should be pointed out that the number and
locations of the points on which we enforce the set of partial differential
equations could be different from the actual training data. Although not
pursued in the current work, this could significantly reduce the required
number of training data on the concentration of the passive scalar.

##  3 Results

To generate high-resolution datasets for different benchmark problems studied
in the current work, we have employed the spectral/hp-element solver NekTar in
which the Navier-Stokes equations (3) along with the transport equation (2)
are approximated using high-order semi-orthogonal Jacobi polynomial expansions
[karniadakis2013spectral]. The numerical time integration is performed using a
third-order stiffly stable scheme until the system reaches its stationary
state. In what follows, a small portion of the resulting dataset corresponding
to this stationary solution will be used for model training, while the
remaining data will be used to validate our predictions. Our algorithm is
agnostic to the choice of initial and boundary conditions as well as the
geometry. However, we choose to provide this information for every benchmark
problem for the sake of completeness and reproducibility of the numerically
generated data. Moreover, all data and codes used in this manuscript are
publicly available on GitHub at https://github.com/maziarraissi/HFM.

To obtain the results reported in the remainder of this manuscript, we
represent each of the functions $c$, $d$, $u$, $v$, $w$, and $p$ by a 10-layer
deep neural network with 50 neurons per hidden layer (see figure 1). As for
the activation functions, we use $\sin(x)$. In general, the choice of a neural
network’s architecture (e.g., number of layers/neurons and form of activation
functions) is crucial and in many cases still remains an art that relies on
one’s ability to balance the trade off between _expressivity_ and
_trainability_ of the neural network [raghu2016expressive]. Our empirical
findings so far indicate that deeper and wider networks are usually more
expressive (i.e., they can capture a larger class of functions) but are often
more costly to train (i.e., a feed-forward evaluation of the neural network
takes more time and the optimizer requires more iterations to converge).
Moreover, the sinusoid (i.e., $\sin(x)$) activation function seems to be
numerically more stable than $\tanh(x)$, at least while computing the residual
neural networks $e_{i}$, $i=1,\ldots,6$ (see equation (5)). However, these
observations should be interpreted as conjectures rather than as firm
results111We encourage the interested reader to check out the codes
corresponding to this paper on GitHub at https://github.com/maziarraissi/HFM
and experiment with different choices for the neural networks’ architectures..
In this work, we have tried to choose the neural networks’ architectures in a
consistent fashion throughout the manuscript. Consequently, there might exist
other architectures that could possibly improve some of the results reported
in the current work.

As for the training procedure, our experience so far indicates that while
training deep neural networks, it is often useful to reduce the learning rate
as the training progresses. Specifically, the results reported in the
following are obtained after $250$, $500$, and $250$ consecutive epochs of the
Adam optimizer [kingma2014adam] with learning rates of $10^{-3}$, $10^{-4}$,
and $10^{-5}$, respectively. Each epoch corresponds to one pass through the
entire dataset. The total number of iterations of the Adam optimizer is
therefore given by $1000$ times the number of data divided by the mini-batch
size. The mini-batch size we used is $10000$ and the number of data points are
clearly specified in the following on a case by case basis. Every $10$
iterations of the optimizer takes around $1.9$ and $3.8$ seconds,
respectively, for two and three dimensional problems on a single NVIDIA Titan
X GPU card.

###  3.1 External flows

As a first example, we consider the prototypical problem of a two dimensional
flow past a circular cylinder, known to exhibit rich dynamic behavior and
transitions for different regimes of the Reynolds number
$\text{Re}=U_{\infty}D/\nu$. Assuming a non-dimensional free stream velocity
$U_{\infty}=1$, cylinder diameter $D=1$, and kinematic viscosity $\nu=0.01$,
the system exhibits a periodic steady state behavior characterized by an
asymmetrical vortex shedding pattern in the wake of the cylinder, known as the
Kármán vortex street [batchelor2000introduction]. Importantly, the passive
scalar is injected at the inlet within the interval [-2.5, 2.5] using a step
function (see figure 2). The choice of this boundary condition only depends on
the region of interest in which velocity and pressure fields are inferred
using the concentration field. Furthermore, a constant value of diffusivity
for the passive scalar is assumed to be given by $\kappa=0.01$ resulting in
$\text{Pec}=U_{\infty}\kappa/\nu=100$. Here, we use zero-slip and zero-
concentration boundary conditions on the cylinder wall. Whereas the Péclet
number is chosen to be equal to the Reynolds number, there is no restriction
on its value as shown for the internal flow benchmark problem. Physically, for
the flow of gases such as air with smoke as a passive scalar
$\text{Pec}\approx\text{Re}$. This is not, however, the case for liquid flows
with dye as the passive scalar since diffusivity of dyes are typically smaller
than most fluids, which leads to higher Péclet numbers.

Figure 2: _2D Flow past a circular cylinder:_ 2D sketch of the simulation
domain in which data on the concentration field is generated. We have assumed
a uniform free stream velocity profile at the left boundary, a zero pressure
outflow condition imposed at the right boundary located $30D$ $(D=1)$
downstream of the cylinder, and periodicity for the top and bottom boundaries
of the $[-10,30]\times[-10,10]$ domain. The passive scalar is injected at the
inlet from [-2.5, 2.5]. A completely arbitrary training domain in the shape of
a flower is depicted in the wake of the cylinder. No information on the
velocity is given for this training domain. Another training domain is also
shown by a rectangle that includes the cylinder. Information on the velocity
is only given on the left boundary of the rectangular domain shown by the
dashed line.

A representative snapshot of the input data on the concentration field is
shown in figure 2. As illustrated in this figure, the shape and extent of the
boundaries of the training domains that we choose for our analysis could be
_arbitrary_ and may vary by problem. However, it should be pointed out that
there are two important factors that need to be considered when choosing the
training domain. First, the concentration field of the passive scalar must be
present within the training domain, such that its information can be used to
infer other flow variables. Second, to avoid the need for specifying
appropriate boundary conditions for velocities, there must exist enough
concentration gradient normal to the boundaries (i.e., $\partial c/\partial
n\neq 0$) in order for our method to be able to infer a unique solution for
the velocity field. As shown in figure 3, excellent agreement can be achieved
between the predictions of our algorithm and the exact data within a
completely arbitrary training domain downstream of the cylinder, while the
inclusion of auxiliary variable $d=1-c$ can improve the accuracy of
predictions further (see figure 4). Here, no input information other than the
concentration field for the passive scalar $c$ is passed to the algorithm.

Figure 3: _2D Flow past a circular cylinder:_ A representative snapshot of the
input data on the concentration of the passive scalar is shown in the top left
panel of this figure, where on the right the same concentration field is
reconstructed based on the predictions of our algorithm. The algorithm is
capable of accurately reconstructing the velocity $u,v$ and the pressure $p$
fields shown in the third row. The exact velocity and pressure fields at the
same point in time are plotted for comparison in the second row. Note that the
pressure is off by a constant since this is an incompressible flow. Figure 4:
_2D Flow past a circular cylinder:_ Relative $L_{2}$ errors between
predictions of the model and the corresponding exact concentration, velocity,
and pressure fields for the arbitrary training domain downstream of the
cylinder. 31 million data points, scattered in space and time, are used both
to regress the concentration field and enforce the corresponding partial
differential equations. Furthermore, the training is performed with the
auxiliary variable $d$ (blue line) and without $d$ (red dashed line) for
comparison.

Next, we choose the region of interest to contain the cylinder (see the
rectangular domain in figure 2) so that the fluid forces acting on the
cylinder can be inferred. As shown in figure 5, the algorithm is capable of
accurately reconstructing the velocity and the pressure fields without having
access to sufficient observations of these fields themselves. Note that we use
a very small training domain that cannot be used in classical computational
fluid dynamics to obtain accurately the Navier-Stokes solutions. It should be
emphasized that other than the velocity on the left boundary, no other
boundary conditions are given to the algorithm. We need to impose a Dirichlet
boundary condition for the velocity at the left boundary simply because the
observations of the concentration are not providing sufficient information
(i.e., gradients) in front of the cylinder. More notably, there is no need to
impose the no-slip condition on the cylinder wall. The presence of the normal
concentration gradient naturally allows our algorithm to infer the zero
velocity condition on the walls. In regards to the predicted pressure field,
we note that due to the nature of the Navier-Stokes equations for
incompressible flows, the pressure field is only identifiable up to a
constant.

Figure 5: _2D Flow past a circular cylinder:_ A representative snapshot of the
input data on the concentration of the passive scalar is shown in the top left
panel of this figure, where on the right the same concentration field is
reconstructed based on the predictions of our algorithm. The algorithm is
capable of accurately reconstructing the velocity $u,v$ and the pressure $p$
fields shown in the third row. The exact velocity and pressure fields at the
same point in time are plotted for comparison in the second row. Note that the
pressure is off by a constant since this is an incompressible flow.

The learned pressure and velocity fields can be used to obtain the lift and
drag forces exerted on the cylinder by the fluid as shown in figure 6. The
fluid forces on the cylinder are functions of the pressure and velocity
gradients. Consequently, having trained the neural networks, we can use

$\displaystyle F_{L}=\oint\left[-pn_{y}+2\text{\footnotesize
Re}^{-1}v_{y}n_{y}+\text{\footnotesize
Re}^{-1}\left(u_{y}+v_{x}\right)n_{x}\right]ds,$ $\displaystyle
F_{D}=\oint\left[-pn_{x}+2\text{\footnotesize
Re}^{-1}u_{x}n_{x}+\text{\footnotesize
Re}^{-1}\left(u_{y}+v_{x}\right)n_{y}\right]ds,$

to obtain the lift and drag forces, respectively. Here, $\bm{n}=(n_{x},n_{y})$
is the outward normal on the cylinder and $ds$ is the arc length on the
surface of the cylinder. We use the trapezoidal rule to approximately compute
these integrals. Interestingly, the predicted lift and drag forces are in
excellent agreement with the exact ones, both in terms of the frequency of
oscillations and the amplitude. The resulting error is within 1% of the test
data. Some discrepancy, however, can be observed for the few initial and final
time instants, which can be attributed to the lack of data. This can be
further clarified by the relative $\mathcal{L}_{2}$-norm of error results
shown in figure 7. Lack of training data on $c$ for $t<0$ and $t>16$ leads to
weaker neural network predictions for the initial and final time instants.
Thus, one should take this into consideration when inference is required
within a certain time interval. Note that to compute lift and drag forces, the
gradient of the velocity has to be computed on the wall, which can be done
analytically using the parametrized surrogate velocity field (i.e., the neural
networks). Although no numerical differentiation is needed to compute the
gradients, the integration of forces on the surface of the cylinder is
approximated by a summation.

Figure 6: _Non-dimensional lift and drag forces on 2D cylinder:_ Comparison
between the exact (spectral element solution) and predicted lift (left) and
drag (right) forces per unit length of the cylinder.

In addition to the velocity and pressure fields, it is possible to discover
other unknown parameters of the flow field such as the Reynolds and Péclet
numbers. Although these parameters were prescribed in the 2D flow past the
cylinder example, we have tested a case in which both parameters are free to
be learned by the algorithm. The results are given in table 1, which shows
very good agreement with the exact values. From the practical standpoint, the
passive scalar diffusivity and the fluid viscosity (hence, their ratio i.e.,
the Prandtl number $\text{Pr}\equiv\nu/\kappa$) may be known in advance.
Therefore, discovering the Reynolds number would be sufficient whereas the
Péclet number can be computed by Pec = Re Pr.

Figure 7: _2D Flow past a circular cylinder:_ Relative $L_{2}$ errors between predictions of the model and the corresponding exact concentration, velocity, and pressure fields for the rectangular training domain that contains the cylinder. 2.8 million data points scattered within the training domain and in time are used for both regressing the concentration field and enforcing the corresponding partial differential equations. Table 1: _2D flow past a circular cylinder:_ Learned Reynolds and Péclet numbers considered as free parameters of the model. | Exact | Learned | Rel. Error
---|---|---|---
Pec | 100 | 92.39 | 7.60%
Re | 100 | 92.47 | 7.52%

The previous benchmark example was a two-dimensional (2D) flow, where we could
safely neglect the $z$-coordinate and $w$-component of the velocity from the
input and output variables, respectively. For the flow past a cylinder, if we
simply increase the Reynolds number beyond a threshold value of $\approx 185$
[ma2002low], the spanwise velocity $w$ becomes more prominent due to the
effect of the so-called “vortex stretching” [batchelor2000introduction,
karniadakis1992three]. To test the capability of the proposed algorithm in
inferring three-dimensional (3D) flow fields, we design another prototype
problem of the 3D flow past a finite-size circular cylinder confined between
two parallel plates as shown in figure 8, where similar to the previous
example, we set $\text{Re}=\text{Pec}=100$. Downstream of the cylinder, the
flow exits to an open region, which causes strong 3D effects in the wake of
the cylinder. Hence, we set the training domain in the wake of cylinder.

Figure 8: _3D flow past a circular cylinder:_ 3D sketch for the prototypical
problem of flow past a circular cylinder of diameter $D=1$. Two parallel
planes are located $10D$ apart along the $z-$axis. Flow is bounded between the
plates for $10D$ along $x-$axis. Periodic boundary conditions are prescribed
parallel to the $xz$ plane at $y=-10$ and $10$, and parallel to the $xy$ plane
at $z=0$ and $10$ starting from the location where the walls end. All the
physical boundaries are shown in black, whereas the gray box located behind
the cylinder shows the domain of interest where model training is performed.
Zero-Neumann boundary conditions are imposed for velocities and concentration
along with zero pressure at the outflow located $30D$ downstream of the
cylinder. A uniform $U_{\infty}=1$ is imposed at the inlet located $8D$
upstream of the cylinder. Passive scaler is injected at the inlet through a
finite region $[-2.5,2.5]\times[0,10]$. Zero-velocity and concentration
boundary conditions are imposed on each physical boundary. Note that no
information on the surfaces of the training (gray box) domain is given.

A representative snapshot of the input data on the concentration field in the
wake of the cylinder is plotted in the top left panel of figure 9. The iso-
surfaces for all of the fields are also plotted for comparison between the
exact data and the predictions of our algorithm. The algorithm is capable of
accurately reconstructing the velocity and the pressure fields without having
access to any observations of these fields. In particular, no information on
the velocity is given on the boundaries of the domain of interest in the
training phase. Qualitative comparison between the predictions of our
algorithm and the exact data shows good agreement for 3D flows as well,
however, the relative $\mathcal{L}_{2}$-norm errors for velocity and pressure
fields are slightly higher than the values for the 2D flow predictions as
shown in figure 10. Similar to 2D results, the neural network training lack
sufficient data for the initial and final time instants that leads to larger
errors in the predictions.

Figure 9: _3D flow past a circular cylinder:_ The iso-surfaces of exact data
and the predictions of our algorithm on the concentration of the passive
scalar are shown in top row for a representative time instant within the
selected training domain. Using the information on the concentration only, the
velocity fields $u,v,w$ are inferred, are shown in the third row, and are
compared with the exact data in the second row. In addition, the exact and
predicted pressure $p$ fields are plotted in the last row. Note the range of
contour levels are set equal between the exact and predicted iso-surfaces of
concentration fields and velocity components for better comparison. Figure 10:
_3D flow past a circular cylinder:_ Relative $L_{2}$ errors between
predictions of the model and the corresponding exact concentration, velocity,
and pressure fields. 37 million data points scattered within the training box
and in time are used for both regressing the concentration field and enforcing
the corresponding partial differential equations.

###  3.2 Internal flows

We now turn our attention to an important class of flows in confined
geometries also known as “internal flows”. While measuring average flow
velocity and pressure in ducts, pipes and even blood vessels is now a common
practice, quantifying the spatial fields, specifically the shear stresses on
the boundaries, is not a trivial task. The first benchmark problem considered
here is the transient flow over an obstacle in a 2D channel as shown in figure
11. To make the flow unsteady, we impose a velocity waveform at the inlet of
the channel. Passive scalar value on the boundaries is set to zero, where the
boundaries are assumed to be impenetrable. It is clear that the auxiliary
variable $d:=1-c$ is equal to one on the channel boundaries. In other words,
while the passive scalar concentration $c$ is convected downstream from the
inflow, the complement scalar $d$ is “virtually” infused from the boundaries
to the flow stream as shown in figure 11. The presence of the obstacle breaks
the symmetry in the flow for which analytical solutions do not exist. Thus, to
estimate the velocity field, direct measurements or forward numerical
simulations are required.

A representative snapshot of the input data on the concentration field along
with its neural network approximation are depicted in the top panel of figure
12. As shown in figure 12, predictions of the algorithm for velocity and
pressure fields are in excellent agreement with the exact data. Note that no
information for velocity is given on the boundaries of the channel or the
obstacle. It should be mentioned that while the concentration gradients normal
to the physical boundaries are sufficient to infer the flow variables without
the knowledge of the boundary conditions (i.e., no-slip velocity on the wall),
the _free_ information from the auxiliary variable $d$ in addition to $c$
helps the algorithm improve the accuracy of predictions (shown in figure 14)
at practically no additional cost.

Figure 11: _2D channel flow over an obstacle:_ Contours of the concentration
field within a 2D channel at a representative time instant. A pulsatile
velocity profile $u(t)$ ($v=0$) is imposed at the inlet (shown as the inset on
the left), and a uniform concentration $c=1$ for the passive scalar is
injected to the channel. At the outlet, a zero-Neumann boundary condition for
the velocity and concentration is considered, while the pressure is set to
zero. Furthermore, the velocity and concentration on the walls are set to
zero. The flow $\text{Re}=\bar{U}H/\nu=60$ is calculated based on the mean
inflow velocity $\bar{U}=1$ and channel hight $H=12$, whereas we choose a
smaller diffusion constant for the passive scalar leading to higher
$\text{Pec}=\bar{U}H/\kappa=180$. The training domain (white rectangle) is
considered to contain the obstacle with the upper and lower boundaries being
the physical wall boundaries. Also shown are contours of the auxiliary
variable $d$, which signifies the important role of this variable in helping
the algorithm detect the boundaries of the training domain.

One of the real-world examples of the above benchmark problem is encountered
in cardiovascular fluid mechanics, where one or multiple coronary arteries are
partially blocked by atherosclerosis plaques formed by the lipid accumulation
[hansson2005inflammation]. An accurate prognosis of these “stenotic” vessels
is an extremely important step in the treatment planning and decision-making
process for patients with cardiovascular disease [budoff2008diagnostic].
Direct measurements of pressure in the vessel is invasive and bears high risk,
whereas computational fluid modeling of blood in the coronaries requires
correct reconstruction of the geometry as well as the knowledge of all the
inflow/outflow boundary conditions. The proposed inverse approach to infer the
velocity and pressure fields with the use of a passive scalar, however, offers
a promising alternative to the conventional methods. Here, the passive scalar
could be the bolus dye that is typically injected to the blood stream for the
purpose of blood flow monitoring and medical imaging.

Figure 12: _2D channel flow over an obstacle:_ A representative snapshot of
the input data on the concentration field within the training domain is
plotted in the top left panel alongside the prediction of our algorithm. The
algorithm is capable of accurately reconstructing the velocity and the
pressure fields without having access to any observations of these fields
(shown in the second and third rows). Furthermore, no boundary conditions are
specified on the boundaries of the training domain including the physical wall
boundaries.

Using the predicted velocity fields, we are able to compute shear stresses
everywhere in the training domain. Of particular interest are wall shear
stresses, which can be computed using the following equations in 2D;

$\begin{array}[]{l}\tau_{x}=2\text{\footnotesize
Re}\left[u_{x}n_{x}+\frac{1}{2}(v_{x}+u_{y})n_{y}\right],\\\
\tau_{y}=2\text{\footnotesize
Re}\left[\frac{1}{2}(u_{y}+v_{x})n_{x}+v_{y}n_{y}\right].\end{array}$

Here, $\bm{n}=(n_{x},n_{y})$ is the outward normal on the boundary of the
domain. Note that to compute wall shear stresses
$\bm{\tau}=(\tau_{x},\tau_{y})$, the gradient of the velocity is required,
which can be computed using automatic differentiation. Wall shear stresses are
important quantities of interest in many biological processes e.g., in the
pathogenesis and progression of vascular diseases that cause aortic aneurysms.
Here, we have estimated the temporal wall shear stress magnitudes,
$\text{WSS}(x,t)=\sqrt[]{\tau_{x}^{2}+\tau_{y}^{2}}$, acting on the lower wall
using the predictions of neural networks for the 2D channel flow over the
obstacle, and plotted them against the results from the spectral/hp element
solver in figure 13. The results show excellent agreement between the
predictions and numerical estimations. Note that the numerical estimations
suffer from a slight aliasing effect (noisy oscillations in the $x$ direction)
for which dealiasing is required to retrieve a smooth distribution of wall
shear stresses. Interestingly, the neural networks prediction is smooth as the
velocity gradients of a parameterized velocity field are taken analytically.
Furthermore, similar to the previous benchmark problems, slight discrepancies
can be observed close to the initial and final times ($t\sim 0,\ 20$) due to
the lack of training data. This is further explained in figure 14, where
relative $\mathcal{L}_{2}$-norm of errors in the velocity fields show a peak
close to the initial and final time instants out of the entire prediction time
interval. Also shown in the figure is the enhancement in the model prediction
accuracy when the auxiliary variable $d$ is used as opposed to when $d$ is not
included during the training phase.

Figure 13: _2D channel flow over an obstacle:_ Wall shear stress magnitude,
$\text{WSS}(x,t)$, computed on the lower wall ($15\leq x\leq 55$) in the
training domain of figure 11 as a function of time ($0\leq t\leq 20$). Strong
wall shear stress values are observed on the surface of the cylinder ($25\leq
x\leq 35$). The plot also shows the periodicity in the flow field.

Similar to the external flow problem past the cylinder, the algorithm is able
to learn the Reynolds and Péclet numbers as free parameters. The predicted
values of these two parameters are given in table 2, which shows excellent
agreement with the exact values. Note that unlike the flow past the cylinder,
the physical Re and Pec numbers are not the same. Remarkably, the algorithm is
capable of inferring the underlying physics without knowledge of initial and
boundary conditions and physical boundaries, where the information is provided
only on the transported passive scalar.

Figure 14: _2D channel flow over an obstacle:_ Relative $L_{2}$ errors between predictions of the model and the corresponding exact concentration, velocity, and pressure fields. 10 million data points scattered within the training domain and in time are used for both regressing the concentration field and enforcing the corresponding partial differential equations. Furthermore, the training is performed with the auxiliary variable $d$ (blue line) and without $d$ (red dashed line) for comparison. Table 2: _2D channel flow over an obstacle:_ Learned Reynolds and Péclet numbers considered as free parameters of the model. | Exact | Learned | Rel. Error
---|---|---|---
Pec | 180 | 178.95 | 0.58%
Re | 60 | 59.92 | 0.12%

To further illustrate the implications of the _Navier-Stokes informed neural
networks_ in addressing real-world problems, we consider 3D physiologic blood
flow in a realistic intracranial aneurysm (ICA) shown in figure 15. The
aneurysm is located in the cavernous segment of the right internal carotid
artery at the level of the eye and beneath the brain [baek2012convergence].
Exact concentration fields are generated numerically using realistic boundary
conditions, which is a physiologic flow waveform at the inlet along with a
uniform concentration for the passive scalar. The strength of the proposed
algorithm is its ability to stay agnostic with respect to the geometry as well
as initial and boundary conditions. Hence, it is possible to focus only on the
regions where velocity and pressure fields are needed, which will
significantly reduce the size of data and the cost of training.

We first crop the aneurysm sac out from the rest of geometry, and then use
only the scalar data within the ICA sac (right panel of figure 15) for
training where no information is used for the boundary conditions. The exact
and predicted concentration, velocity and pressure fields within the ICA sac
at a sample time instant are then interpolated on two separate planes
perpendicular to $y$\- and $z$-axis, which are shown in figure 16. We observe
excellent agreement between the exact and predicted fields given the
complexity of the flow field. Furthermore, the relative $\mathcal{L}_{2}$-norm
of errors plotted in figure 17 show a significant drop in the prediction
errors away from the initial and final time instants, which justifies an
accurate analysis of model predictions in the time interval $10\leq t\leq 25$.

Figure 15: _A 3D intracranial aneurysm:_ The middle panel shows the simulation
domain and pressure field at a time instant, whereas on the right the training
domain containing only the ICA sac is shown. Two perpendicular planes have
been used to interpolate the exact data and the predicted ones for plotting 2D
contours in figure 16. A physiologic flow waveform $Q(t)$ (shown in the inset
figure) is prescribed at the inlet along with the uniform concentration for
the passive scalar. At the outlet, zero-Neumann boundary conditions are
imposed for velocities and concentration, whereas a “Windkessel” type boundary
condition is used for the pressure to represent the truncated geometry
downstream [grinberg2008outflow]. The Reynolds and Péclet numbers are
estimated based on the mean velocity and lumen diameter at the inlet. Using
the kinematic viscosity of blood and assuming a diffusivity constant $\kappa$
for the passive scalar equal to the viscosity, we obtain
$\text{Re}=\text{Pec}=98.2$. Figure 16: _A 3D intracranial aneurysm:_ Contours
of the exact fields and model predictions are plotted on two perpendicular
planes (shown in figure 15) for concentration $c$, velocity $u,v,w$, and
pressure $p$ fields in each row. The first two columns show the results
interpolated on a plane perpendicular to the $z$-axis, and the next two
columns are plotted for a plane perpendicular to the $y$-axis. Note that the
range of contour levels are set equal for all fields for better comparison.
Figure 17: _A 3D intracranial aneurysm:_ Relative $L_{2}$ errors between
predictions of the model and the corresponding exact concentration, velocity,
and pressure fields. 29 million data points scattered within the aneurysm sac
and in time are used for both regressing the concentration field and enforcing
the corresponding partial differential equations.

##  4 Discussion and Concluding Remarks

The algorithm developed here is agnostic to the geometry, initial, and
boundary conditions, which makes it highly flexible in choosing the domain of
interest for data acquisition as well as subsequent training and predictions.
Moreover, the current methodology allows us to construct computationally
efficient and fully differentiable surrogates for velocity and pressure fields
that can be further used to estimate other quantities of interest such as
shear stresses and vorticity fields. The predictions presented here are among
those cases where a pure machine learning algorithm or a mere scientific
computing approach simply cannot reproduce. A pure machine learning strategy
has no sense of the physics of the problem to begin with, and a mere
scientific computing approach relies heavily on careful specification of the
geometry as well as initial and boundary conditions.

Assuming the geometry to be known, to arrive at similar results as the ones
presented in the current work, one needs to solve significantly more expensive
optimization problems, using conventional computational methods (e.g., finite
differences, finite elements, finite volumes, spectral methods, and etc.). The
corresponding optimization problems involve some form of “parametrized”
initial and boundary conditions, appropriate loss functions, and multiple runs
of the conventional computational solvers. In this setting, one could easily
end up with very high-dimensional optimization problems that require either
backpropagating through the computational solvers [chen2018neural] or
“Bayesian” optimization techniques [shahriari2016taking] for the surrogate
models (e.g., Gaussian processes). If the geometry is further assumed to be
unknown (as is the case in this work), then its parametrization requires grid
regeneration, which makes the approach almost impractical.

It must be mentioned that we are avoiding the regimes where the Navier-Stokes
equations become chaotic and turbulent (e.g., as the Reynolds number
increases). In fact, it should not be difficult for a plain vanilla neural
network to approximate the types of complicated functions that naturally
appear in turbulence. However, as we compute the derivatives required in
equation (5), minimizing the loss function (2) might become a challenge
[raissi2018deep], where the optimizer may fail to converge to the right values
for the parameters of the neural networks. It might be the case that the
resulting optimization problem inherits the complicated nature of the
turbulent Navier-Stokes equations. Hence, inference of turbulent velocity and
pressure fields could be considered in future extensions of this line of
research.

The emphasis of this study was to demonstrate, through several prototypical
and realistic examples, the ability of the current algorithm to infer the
hidden states of the system from partial knowledge of some relevant quantities
by leveraging the known underlying dynamics of the system. We have verified
that with the use of governing physical laws, our algorithm is able to make
accurate predictions for complex 2D/3D flows. One possible limitation of the
current study is the use of synthetically generated data on the passive
scalar, which are relatively noiseless and clean compared to the realistic
measurements. Whereas 3D reconstruction of a passive scalar from a stack of 2D
projections is possible, it is not a common practice in the clinical and
industrial settings due to the technological complexities as well as
computational and processing costs. Thus, extracting information from 2D
images directly by taking into account the angle of projections seems to be a
more realistic approach. Furthermore, different imaging modalities (e.g.,
magnetic resonance vs. computed tomography angiography of bolus dye in
coronary arteries) have different spatial/temporal resolutions. To address the
lack of resolution in time, a viable strategy is to use a high-order time-
stepping scheme such as implicit Runge-Kutta for resolving the equations in
time using very few snapshots [raissi2017physics_I, raissi2017physics_II].
These are certainly important factors that have to be taken into account
moving forward using the proposed algorithm for real-world fluid mechanics
applications, and will be addressed carefully in future implementations.

In this work we have been operating under the assumption of Newtonian and
incompressible fluid flow governed by the Navier-Stokes equations. However,
the proposed algorithm can also be used when the underlying physics is non-
Newtonian, compressible, or partially known. This in fact is one of the
advantages of the algorithm in which other unknown parameters such as the
Reynolds and Péclet numbers can be inferred in addition to the velocity and
pressure fields. When the fluid is non-Newtonian (e.g., blood flow in small
vessels), one can encode the momentum equations, where the “constitutive” law
for the fluid’s stress-strain relationship is unknown, into a physics-informed
deep learning algorithm. Having data on the velocity (or even the passive
scalar), it should be possible to learn the constitutive law as well
[zhao2018active].

As demonstrated in this work, a direct implication of the current method is
quantifying the hemodynamics in the vasculature. This could potentially have
significant impact on the clinical diagnosis (especially with the non-invasive
methods) of vascular diseases associated with important pathologies such as
heart attack and stroke. Blood shear stresses acting on the vascular wall are
crucial in the prognosis of a vascular disease and their quantification is
significantly important clinically [shaaban2000wall, zarins1983carotid]. Using
the proposed method, it is possible to estimate the wall shear stresses at no
extra cost. This will simplify the complexities of the state-of-the-art
methods in which extracting the exact boundaries of the vessels from the
clinical images is required [boussel2009phase]. In addition to the shear
stresses estimation that is of interest, quantification of blood distribution
and oxygenation in cerebral arteries using modalities such as perfusion CT
[wintermark2002prognostic] and functional magnetic resonance imaging
[calamante1999measuring] has been practiced and is crucial in the assessment
of brain activity for the majority of neurodegenerative diseases such as
Alzheimer’s disease [greicius2004default] and in cognitive neuroscience
[aron2007triangulating]. Furthermore, recent advances in imaging such as the
functional photoacoustic microscopy allows 3D blood oxygenation imaging with
“capillary-level” resolution [yao2015high]. Encoding the one-dimensional
hemodynamics and oxygen transport equations into physics-informed neural
networks, and using the concentration profiles from different imaging
modalities will lead to quantifying brain hemodynamics at various length
scales.

In conclusion, we have introduced an effective algorithm based on the recently
developed _physics-informed neural networks_ framework that is capable of
encoding the underlying physical laws that govern a given dataset.
Specifically, we have applied the algorithm to an important class of physical
laws governing fluid motions. Our inverse approach in this article is unique
as we have designed a methodology to infer velocity and pressure fields merely
from the knowledge of the time-evolution of a passive scaler.

## Acknowledgements

This work received support by the DARPA EQUiPS grant N66001-15-2-4055, the
AFOSR grant FA9550-17-1-0013, the NIH grant U01HL116323, and the NSF grant
DMS-1736088. Generation of data was performed on XSEDE resources supported by
award No. TG-DMS140007. Moreover, all data and codes used in this manuscript
will be publicly available on GitHub at https://github.com/maziarraiss/HFM.


# A Parts Based Registration Loss for Detecting Knee Joint Areas

Juha Tiirola Email<EMAIL_ADDRESS>Research Unit of Mathematical
Sciences, University of Oulu, Finland

###### Abstract

In this paper, a parts based loss is considered for finetune registering knee
joint areas. Here the parts are defined as abstract feature vectors with
location and they are automatically selected from a reference image. For a
test image the detected parts are encouraged to have a similar spatial
configuration than the corresponding parts in the reference image.

Keywords: Automatic Detection, Convolutional Neural Networks, Spatial
Transformer

##  1 Introduction

Neural networks based general object detection methods need lot of supervision
in the form of bounding box annotations. Some well-known methods in this
direction are Fast R-CNN [1], Faster R-CNN [2] and Yolo [3]. In contrast to
general object detection, in the knee joint area detection problem from
bilateral PA fixed flexion X-ray images the amount of variation between
samples is small since the pose is shared between images and due to human knee
anatomy, there are no big deformations between samples.

In [4], the detection of knee joint areas was considered based on template
matching. There the registration loss was based on the use of the normalized
cross-correlation. In contrast to [4], in this paper the detection is based on
points matching. More precisely, the points here are abstract feature vectors.
By parts we refer to abstract feature vectors equipped with location and which
are extracted from a reference image. The term ’part’ was used for instance in
[5] where a dictionary of parts was learned by clustering feature map vectors.
In [6], landmark localization was considered using convolutional neural
networks where spatial configuration was integrated into heatmap regression.
In [7], human pose was estimated using a convolutional network where the
architecture could exploit structural domain constraints such as geometric
relationships between body joint locations. We use a similar idea to deal with
the parts. Since the camera pose between different bilateral fixed flexion
knee X-ray images is shared, the parts should be such that there is a common
pattern between the locations of the parts between different images. More
precisely, given a test image, the extracted patch from the test image should
be such that each part in the extracted patch is near the corresponding part
in the reference image and in addition the content of each part in the
extracted patch looks like the corresponding content in the reference image.
The term ’looks like that’ was used in [8] where explainable image
classification was considered based on selecting prototypical parts and by
making the classification based on how much each prototype is present in the
input image.

##  2 New method

(a) $a$

(b) Architecture

Figure 1: Here $a$ is the manually selected reference image, $N$ is the number
of the parts, $u(\boldsymbol{\theta})$ is the extracted patch from the input
image $u$ corresponding to the parameter $\boldsymbol{\theta}$. Heatmap
$k^{i}(\boldsymbol{\epsilon}^{i})$ is a Gaussian blob centered at
$\boldsymbol{\epsilon}^{i}+\boldsymbol{\ell}^{i}$.
$g^{i}(u(\boldsymbol{\theta}))$ is obtained by jointly denoising part
detection heatmaps
$h^{1}(u(\boldsymbol{\theta})),\dots,h^{N}(u(\boldsymbol{\theta})).$ In the
first phase, $\boldsymbol{\ell}^{i},i=1,\dotsc,N$, are learnt using the whole
dataset. In the second phase, given a test image $u$, image specific
parameters $\boldsymbol{\theta}\in\mathbb{R}^{4}$ and
$\boldsymbol{\epsilon}^{i}\in\mathbb{R}^{2},$ $i=1,\dotsc,N$, are optimized by
minimizing energy (2).

The reference image $a$ and the architecture of the new method are shown in
Figure 1. We first select a reference image $a$ manually. Given a test image
$u$, we try to extract a patch from $u$ which resembles $a$ as much as
possible. We formulate this similarity using parts. We use the VGG16 network
[9] pretrained on the ImageNet as the backbone network. We freeze the weights
of the backbone. We denote by $V\in\mathbb{R}^{28\times 28\times 512}$ the
feature tensor from the 22th layer of the VGG16 network followed by
normalization. Thus, $||V_{r,s,:}||_{2}=1$ for all $r,s$.

We first determine the parts from the reference image utilizing the whole
dataset. In the $[-1,1]^{2}$ image coordinate domain, we assume
$\boldsymbol{\ell}^{i}\in[-1,1]^{2}$, $i=1,\dotsc,N$, denote the centers of
the parts in image $a$. We use $N=9$ and initialize the
$\boldsymbol{\ell}^{i}$ as points from a regular grid
$(-0.5,-0.5)+r(0.5,0)+s(0,0.5)$ where $r,s\in\\{0,1,2\\}$. The locations
$\boldsymbol{\ell}^{i}$ are optimized over a dataset. For each $i$ we
interpolate $V(a)$ at $\boldsymbol{\ell}^{i}$. The resulting vectors are
normalized and we get $V(a)_{\boldsymbol{\ell}^{i}}\in\mathbb{R}^{512}$ which
encodes the content of the part $i$. Thus,
$||V(a)_{\boldsymbol{\ell}^{i}}||_{2}=1.$ For an image $u$ from the dataset,
we define $h^{i}(u)\in[0,1]^{28\times 28}$ by

$h^{i}(u)_{r,s}:=\langle
V(u)_{r,s},V(a)_{\boldsymbol{\ell}^{i}}\rangle_{L^{2}}$ (1)

where $\langle\cdot,\cdot\rangle_{L^{2}}$ denotes the $L^{2}$ inner product
between two vectors. Thus, $h^{i}(u)$ is the heatmap which tells how much in
$u$ there is part $i$ in each location. The locations are optimized over the
whole dataset by solving

$\min_{\\{\boldsymbol{\ell}^{i}\\}_{i=1}^{N}}\sum_{u}\sum_{i=1}^{N}||h^{i}(u)||_{1}+\lambda(1-\max_{r,s}h^{i}(u)_{r,s})$

where we use $\lambda=0.1$. Ideally, the locations should be such that in an
image $u$ the part $i$ is detected in a single location. In the above
minimization the first term encourages sparsity in $h^{i}(u)$ and the second
loss encourages that the maximum value in $h^{i}(u)$ is near $1$. After the
optimization, we freeze the locations $\\{\boldsymbol{\ell}^{i}\\}_{i=1}^{N}.$

The second step is to extract a patch from each test image $u$ such that the
extracted patch resembles $a$ visually as much as possible. This optimal
extraction is formulated as an optimization problem. There are $4+2N$
variables in the minimization problem. The first four variables are
transformation parameters $\boldsymbol{\theta}\in\mathbb{R}^{4}$ for the
spatial transformer [10] which determine the extracted patch from the input
image. The extracted patch is denoted by $u(\boldsymbol{\theta})$. The rest of
the optimization variables are disturbance vectors. Since there are slight
deformations between knee images, ideally the location of the part $i$ in a
well registered test image is the location of the part in the reference image
plus some small disturbance vector. This disturbance vector is denoted by
$\boldsymbol{\epsilon}^{i}$ for part $i$.

The heatmaps $h^{i}(u(\boldsymbol{\theta}))$ obtained using (1) denote the
detected parts $i$ in the patch $u(\boldsymbol{\theta})$. These heatmaps
$h^{i}$ are noisy and ambiguous. We use a method similar to [6] and [7] to
denoise the parts heatmaps $h^{j}$ jointly by utilizing the spatial
configuration between the locations of the parts. If image $v$ resembles the
reference image $a$ and if there is a peak in $h^{i}(v)$ near
$\boldsymbol{\ell}^{i}$, then if we move from the peak along the vector
$\boldsymbol{\ell}^{j}-\boldsymbol{\ell}^{i}$, then near the resulting
position there should be a peak in $h^{j}(v)$. The denoised version of $h^{j}$
is denoted by $g^{j}$ and is obtained by $g^{j}(v)=M_{j}(v)\odot h^{j}(v)$
where the multiplier image is

$M_{j}(v):=\frac{1}{N-1}\sum_{r=1,r\neq
j}^{N}\max_{\overline{k}}\left(G(\overline{k})\,\tau_{\boldsymbol{\ell}^{j}-\boldsymbol{\ell}^{r}+\overline{k}}(h^{r}(v))\right).$

Above $\tau$ denotes the translation operator and $G$ is the Gaussian kernel
with the standard deviation $0.08$ in $[-1,1]^{2}$ domain. In the multiplier,
$(M_{j}(v))_{a,b}$ is large when for as many $r$ as possible, $r\neq j$,
$h^{r}(v)$ is large near $(a,b)+\boldsymbol{\ell}^{j}-\boldsymbol{\ell}^{r}.$

We also model the locations of the parts directly. We assume that the location
of part $i$ in a well extracted patch from $u$ is
$\boldsymbol{\ell}^{i}+\boldsymbol{\epsilon}^{i}$ where
$||\boldsymbol{\epsilon}^{i}||$ is small. We place a Gaussian blob with the
standard deviation $0.08$ at the location
$\boldsymbol{\ell}^{i}+\boldsymbol{\epsilon}^{i}$. The heatmap obtained in
this way is denoted by $k^{i}(\boldsymbol{\epsilon}^{i})$.

The parameters $\boldsymbol{\epsilon}^{i}$, $\boldsymbol{\theta}$ are found
jointly by solving

$\displaystyle\min_{\\{\boldsymbol{\epsilon}^{i}\\}_{i=1}^{N},\boldsymbol{\theta}}-\sum_{n=1}^{N}\langle
1+k^{n}(\boldsymbol{\epsilon}^{n}),g^{n}(u(\boldsymbol{\theta}))\odot
g^{n}(u(\boldsymbol{\theta}))\rangle_{L^{2}}+\lambda\sum_{n=1}^{N}||\boldsymbol{\epsilon}^{n}||_{2}^{2}$
$\displaystyle=\min_{\\{\boldsymbol{\epsilon}^{i}\\}_{i=1}^{N},\boldsymbol{\theta}}-\sum_{n=1}^{N}\sum_{r,s}(1+k^{n}(\boldsymbol{\epsilon}^{n})_{r,s})(g^{n}(u(\boldsymbol{\theta}))_{r,s})^{2}+\lambda\sum_{n=1}^{N}||\boldsymbol{\epsilon}^{n}||_{2}^{2}.$
$\displaystyle=:\min_{\\{\boldsymbol{\epsilon}^{i}\\}_{i=1}^{N},\boldsymbol{\theta}}E(u,\boldsymbol{\epsilon}^{1},\dotsc,\boldsymbol{\epsilon}^{N},\boldsymbol{\theta}).$
(2)

The loss to be minimized incorporates the requirement that in a well extracted
patch the location of the part $i$ should be close to $\boldsymbol{\ell}^{i}$
and in the heatmap $h^{i}$ and $g^{i}$ there should be a peak near
$\boldsymbol{\ell}^{i}$.

##  3 Parametrization

We transform minimization problem (2) to an unconstrained minimization
problem. See also [4] for the interpretation of the parameters. We denote by
$(\boldsymbol{w},\boldsymbol{v})\in\mathbb{R}^{2N+4}$ the unconstrained
parameter vector. We express $\boldsymbol{\theta}$,
$\boldsymbol{\epsilon}^{1},\dotsc,\boldsymbol{\epsilon}^{N}$ as functions of
$(\boldsymbol{w},\boldsymbol{v})$. For
$\\{\boldsymbol{\epsilon}^{i}\\}_{i=1}^{N}$ we use
$(\boldsymbol{\epsilon}^{1},\dotsc,\boldsymbol{\epsilon}^{i},\dotsc,\boldsymbol{\epsilon}^{N})=\frac{1}{4}\tanh(\boldsymbol{w})$
and we initialize $\boldsymbol{w}=\boldsymbol{0}\in\mathbb{R}^{2N}$. For the
scale we use

$\boldsymbol{\theta}_{1}(\boldsymbol{v}_{1})=\left(s_{1}-\frac{1}{4}\right)+\frac{1}{2}\left(\sigma(\boldsymbol{v}_{1})\right)=\left(s_{1}-\frac{1}{4}\right)+\frac{1}{2}\left(\frac{1}{2}\left(1+\tanh(\boldsymbol{v}_{1})\right)\right)$
(3)

where $s_{1}$ is fixed. For the translations we use

$\boldsymbol{\theta}_{2}(\boldsymbol{v})=(1-\boldsymbol{\theta}_{1}(\boldsymbol{v}_{1}))(-1+2\sigma(\boldsymbol{v}_{2}))\text{
and
}\boldsymbol{\theta}_{3}(\boldsymbol{v})=(1-\boldsymbol{\theta}_{1}(\boldsymbol{v}_{1}))(-1+2\sigma(\boldsymbol{v}_{3}))$

and for the rotation we use

$\boldsymbol{\theta}_{4}(\boldsymbol{v}_{4})=\frac{1}{10}\tanh(\boldsymbol{v}_{4}).$

We initialize $\boldsymbol{v}=\boldsymbol{0}\in\mathbb{R}^{4}$.

The transformation matrix corresponding to $\boldsymbol{v}$ is given by

$A(\boldsymbol{\theta}(\boldsymbol{v})):=\begin{bmatrix}\boldsymbol{\theta}_{1}\cos(\boldsymbol{\theta}_{4})&-\boldsymbol{\theta}_{1}\sin(\boldsymbol{\theta}_{4})&\boldsymbol{\theta}_{2}\\\
\boldsymbol{\theta}_{1}\sin(\boldsymbol{\theta}_{4})&\boldsymbol{\theta}_{1}\cos(\boldsymbol{\theta}_{4})&\boldsymbol{\theta}_{3}\end{bmatrix}.$
(4)

We transform problem (2) to an unconstrained minimization problem

$\min_{\\{\boldsymbol{w}_{i}\\}_{i=1}^{N},\boldsymbol{v}}E(u,\boldsymbol{\epsilon}^{1}(\boldsymbol{w}_{1}),\dotsc,\boldsymbol{\epsilon}^{N}(\boldsymbol{w}_{N}),\boldsymbol{\theta}(\boldsymbol{v})).$
(5)

##  4 Experimental results

We test the minimization of (5) for several scales $s_{1}$ in (3). We do the
minimization for all $s_{1}\in\\{0.65,0.70,0.75,\dotsc,1.15,1.20\\}$ and we
take the patch corresponding to the smallest overall loss value. We use
$\lambda=10^{-6}$ in equation (2).

We only consider finetune registration. We use as inputs $u$ randomly slightly
enlargened outputs of the neural method [4] since our goal is to test if there
is a visual correspondence between the reference image $a$ and the minimizer
of (5).

|  |
---|---|---
|  |
|  |
|  |
|  |
|  |
|  |
|  |
|  |
|  |
|  |

Figure 2: In each column separated by vertical bar, a left image is $u$ and
the right image is $u(\boldsymbol{\theta})$.

|  |
---|---|---
|  |
|  |
|  |
|  |
|  |
|  |
|  |
|  |
|  |
|  |

Figure 3: In each column separated by vertical bar, a left image is $u$ and
the right image is $u(\boldsymbol{\theta})$.

In Figures 2 and 3 we see some numerical results. In the figures, there are
three main columns separated by vertical bars and in each row of those columns
there is a $(u,u(\boldsymbol{\theta}))$ pair where $u$ is the neural network
[4] output with some random added neighbourhood and $u(\boldsymbol{\theta})$
is a candidate for the minimizer of (2). From the results we see that there is
a closer resemblance between $u(\boldsymbol{\theta})$ and the reference image
$a$ than between $u$ and $a$. There are also less variations in sizes between
the $u(\boldsymbol{\theta})$ than between the $u$.

##  5 Conclusion

In this paper, a parts based loss was considered for finetune registering
X-ray knee joint areas using a simple template image. The parts were first
selected from a reference image. Then a patch is extracted from a test image
such that the detected parts in the patch have a similar spatial configuration
than the parts in the reference image. In the experiments, it seemed visually
that on average for input images roughly representing the knee joint areas the
distance between the minimizers of the parts based loss function and the
reference image is smaller than the distance between the input images and the
reference image.

## References

  * [1] R. Girshick.  Fast R-CNN.  In 2015 IEEE International Conference on Computer Vision (ICCV), pages 1440–1448, 2015.
  * [2] S. Ren, K. He, R. B. Girshick and J. Sun.  Faster R-CNN: Towards Real-Time Object Detection with Region Proposal Networks.  In Advances in Neural Information Processing Systems, pages 91–99, 2015.
  * [3] J. Redmon, S. K. Divvala, R. B. Girshick and A. Farhadi.  You Only Look Once: Unified, Real-Time Object Detection.  In CVPR 2016, pages 779–788, 2016.
  * [4] J. Tiirola.  A Neural Template Matching Method to Detect Knee Joint Areas.  In arXiv, 2209.11791, 2022.
  * [5] A. Kortylewski, Q. Liu, H. Wang, Z. Zhang and A. L. Yuille.  Combining Compositional Models and Deep Networks For Robust Object Classification under Occlusion.  In IEEE Winter Conference on Applications of Computer Vision, WACV 2020, pages 1322–1330, 2020.
  * [6] C. Payer, D. Štern, H. Bischof and M. Urschler.  Integrating spatial configuration into heatmap regression based CNNs for landmark localization.  In Medical Image Analysis, volume 54, pages 207–219, 2019.
  * [7] J. J. Thompson, A. Jain, Y. LeCun and C. Bregler.  Joint Training of a Convolutional Network and a Graphical Model for Human Pose Estimation.  In Advances in Neural Information Processing Systems, volume 27, pages 1799–1807, 2014.
  * [8] C. Chen, O. Li, C. Tao, A. Barnett, J. Su, C. Rudin.  This Looks Like That: Deep Learning for Interpretable Image Recognition.  In Advances in Neural Information Processing Systems, volume 32, pages 8930–8941, 2019.
  * [9] K. Simonyan and A. Zisserman.  Very Deep Convolutional Networks for Large-Scale Image Recognition.  In 3rd International Conference on Learning Representations, 2015.
  * [10] M. Jaderberg, K. Simonyan, A. Zisserman, and K. Kavukcuoglu.  Spatial transformer networks.  In Advances in Neural Information Processing Systems, volume 28, 2015.