# Macrophages trajectories smoothing by evolving curves

Giulia Lupi, Karol Mikula, Seol Ah Park

###### Abstract

When analyzing cell trajectories, we often have to deal with noisy data due to
the random motion of the cells and possible imperfections in cell center
detection. To smooth these trajectories, we present a mathematical model and
numerical method based on evolving open-plane curve approach in the Lagrangian
formulation. The model contains two terms: the first is the smoothing term
given by the influence of local curvature, while the other attracts the curve
to the original trajectory. We use the flowing finite volume method to
discretize the advection-diffusion partial differential equation. The PDE
includes the asymptotically uniform tangential redistribution of curve grid
points. We present results for macrophage trajectory smoothing and define a
method to compute the cell velocity for the discrete points on the smoothed
curve.

##  1 Mathematical model

In this paper, we present a mathematical model and numerical method for
macrophage trajectories smoothing based on evolving open plane curve approach
in the Lagrangian formulation, namely, we solve the equation

$\frac{\partial\textbf{x}}{\partial t}=-\delta
k\textbf{N}+\lambda[(\textbf{x}_{0}-\textbf{x})\cdot\textbf{N}]\textbf{N}+\alpha\textbf{T},$
(1)

where x is any point of the evolving curve. The terms evolving the curve in
the normal direction N and the tangential direction T are given in details
below.
The paper is structured as follows: in this section, we discuss the
mathematical model and its possible applications. In section 2, we show the
numerical discretization of the suggested model. In section 3, we expose the
numerical experiments made on cell trajectories. Moreover, the smoothed curves
allow us to compute the velocity of the cell. To do that, we propose a method
to find the new velocity for the discrete points on the evolving curve based
on the evolution of the length of the segments composing the piecewise linear
initial condition.

Let $\Gamma$ be an open plane curve with fixed endpoints

$\begin{split}\Gamma:[0&,1]\rightarrow\mathbb{R}^{2},\\\
&u\mapsto\textbf{x}(u),\end{split}$ (2)

where $\textbf{x}(u)=(x_{1}(u),x_{2}(u))$ is a point of the curve $\Gamma$. We
consider the time evolution of the curve $\Gamma$,
$\Gamma_{t}=\\{\textbf{x}(t,u),u\in[0,1]\\}$ where $t$ is the time and
$\textbf{x}(t,u)$ represents the position of the point
$\textbf{x}(u)\in\Gamma_{t}$ at time $t$. From now on, we will write x instead
of $\textbf{x}(t,u)$ for clarity.
The evolution of a point $\textbf{x}\in\Gamma$ is driven by the following
general equation

$\frac{\partial\textbf{x}}{\partial t}=\textbf{V}(\textbf{x},t),$ (3)

where $\frac{\partial\textbf{x}}{\partial t}$ is the time derivative
representing the velocity of the evolving curve. The velocity vector field V
that drives the curve motion is given by

$\textbf{V}(\textbf{x},t)=-\delta
k(\textbf{x},t)\textbf{N}(\textbf{x},t)+\lambda[(\textbf{x}_{0}-\textbf{x})\cdot\textbf{N}(\textbf{x},t)]\textbf{N}(\textbf{x},t),$
(4)

where N is the unit normal vector to the evolving curve, $k$ is its curvature,
$\textbf{x}_{0}$ is the initial condition, " $\cdot$ " is the scalar product
and $\delta,~{}\lambda$ are given positive parameters. From now on, we will
indicate by

$w(\textbf{x},t)=(\textbf{x}_{0}-\textbf{x})\cdot\textbf{N}(\textbf{x},t).$
(5)

The curvature term $-k\textbf{N}$, weighted by the parameter $\delta$,
regularizes the curve while $w$, weighted by the parameter $\lambda$, attracts
the evolving curve next to the initial condition.
The function $(\textbf{x}_{0}-\textbf{x})$ is the vector that realizes the
minimum distance between the point x and the initial curve $\textbf{x}_{0}$.
We define it in the following way: fix a point
$\textbf{x}=\textbf{x}(u)\in\Gamma_{t}$ on the evolving curve $\Gamma_{t}$,
parametrized by $u\in[0,1]$. Then define

$\begin{split}(\textbf{x}_{0}-\textbf{x})(u)=\operatorname*{arg\,min}_{\textbf{v}\in\chi_{u}}\lvert\textbf{v}\rvert,\end{split}$
(6)

where
$\chi_{u}=\\{\textbf{v}\in\mathbb{R}^{2}:\textbf{v}=\textbf{x}_{0}(q)-\textbf{x}(u),q\in[0,1]\\}$.
For every point on the evolving curve, parametrized by $u\in[0,1]$,
$(\textbf{x}_{0}-\textbf{x})$ selects the point on the original curve,
parametrized by $q\in[0,1]$, which realizes the shortest distance and
considers the vector $\textbf{v}=\textbf{x}_{0}(q)-\textbf{x}(u)$. Then, since
the tangential velocity does not influence the shape of the evolving curve, we
consider only $w(\textbf{x},t)$, the component of the vector in the normal
direction.
We rewrite equation (3) as

$\frac{\partial\textbf{x}}{\partial t}=\beta\textbf{N},$ (7)

where the normal velocity $\beta$ is given by

$\beta=-\delta k+\lambda w.$ (8)

On the other hand, from the numerical point of view, the tangential velocity
redistributes the points on the curve (see [3, 5, 8, 9]): this improves both
the stability and the robustness of numerical computations. To redistribute
the points asymptotically uniformly, we enrich the model as follows

$\frac{\partial\textbf{x}}{\partial t}=\beta\textbf{N}+\alpha\textbf{T},$ (9)

where $\beta$ is defined in (8) and $\alpha$ will be defined in section 2.1.

Equation (1) can be helpful in different applications where one needs to
smooth a curve and remove the noise while keeping the curve similar to the
original. In [2, 11, 12], models have been proposed for segmentation and
tracking of immune system cells. Once these trajectories are obtained, one may
be interested in finding the actual velocity of the cells (see for example
[4]). We applied the proposed model for smoothing macrophage trajectories. Our
goal was to derive the actual velocity of the cell starting from the curve
found connecting the cell center in every frame of a $2$D video [12].
Consequently, we needed to remove the random motion of the cell and correct
possible imperfections in cell centers detection while keeping the final
trajectory as close to the original as possible. We discuss the results in
Section 3.

##  2 Numerical Algorithm

To apply the numerical scheme, the curve is discretized to a set of points,
$\textbf{x}_{0},~{}\textbf{x}_{1},...,~{}\textbf{x}_{n+1}$ as displayed in
Fig. 1.
Let $\lvert\textbf{x}_{u}\rvert>0$, where $\textbf{x}_{u}=(\frac{\partial
x_{1}}{\partial u},\frac{\partial x_{2}}{\partial u})$ and
$g=\lvert\textbf{x}_{u}\rvert=\sqrt{(\frac{\partial x_{1}}{\partial
u})^{2}+(\frac{\partial x_{2}}{\partial u})^{2}}$. If we denote by $s$ the
unit arc-length parametrization of the curve $\Gamma$, then
$ds=\lvert\textbf{x}_{u}\rvert du=gdu$ and $du=\frac{1}{g}ds$. The unit
tangent vector T is defined as $\textbf{T}=\frac{\partial\textbf{x}}{\partial
s}=\textbf{x}_{s}$ while the unit normal vector N is
$\textbf{N}=\textbf{x}_{s}^{\perp}$ such that $\textbf{T}\wedge\textbf{N}=-1$.
If $\textbf{T}=(\frac{x_{1}}{\partial s},\frac{x_{2}}{\partial s})$ then
$\textbf{N}=(\frac{x_{2}}{\partial s},-\frac{x_{1}}{\partial s})$. From the
Frenet-Serret formulas we have $\textbf{T}_{s}=-k\textbf{N}$ and
$\textbf{N}_{s}=k\textbf{T}$, where $k$ is the curvature. It follows that
$-k\textbf{N}=\textbf{T}_{s}=(\textbf{x}_{s})_{s}=\textbf{x}_{ss}$.
We can rewrite (9) into the form of the so-called intrinsic partial
differential equation

$\textbf{x}_{t}=\delta\textbf{x}_{ss}+\alpha\textbf{x}_{s}+\lambda
w\textbf{x}_{s}^{\perp}$ (10)

which is suitable for numerical discretization. Since
$\textbf{x}=(x_{1},x_{2})$, (10) represents a system of two partial
differential equations for the two components of the position vector x.

###  2.1 Suitable choice of tangential velocity

In [3, 5, 8, 9], it has been shown how to choose the tangential velocity to
have an asymptotically uniform redistribution of points; in this section, we
recall it for completeness.
If we want to redistribute the points along the curve, we have to consider the
ratio $\frac{g}{L}$, where $g=\lvert\partial_{u}\textbf{x}\rvert$ represents
the local length of the curve while $L$ is the global length. In the discrete
form, we obtain

$\frac{g}{L}\approx\frac{\frac{\lvert\textbf{x}_{i}-\textbf{x}_{i-1}\rvert}{h}}{L}=\frac{\lvert\textbf{x}_{i}-\textbf{x}_{i-1}\rvert}{Lh}=\frac{\lvert\textbf{x}_{i}-\textbf{x}_{i-1}\rvert}{\frac{L}{n+1}},$

where $h=\frac{1}{n+1}$. We define the $i$-th element
$\textbf{h}_{i}=\textbf{x}_{i}-\textbf{x}_{i-1}$ with length
$h_{i}=\lvert\textbf{h}_{i}\rvert$. In the formula above, $n+1$ is the number
of elements used in the spatial discretization based on the flowing finite
volume method. The numerator represents the distance between two grid points;
on the other hand, the denominator represents the distance that would occur if
we had uniformly distributed points. It is then clear that if we want
uniformly distributed points, we will have to ask for that ratio to tend to
$1$. In continuous setting, we have to fulfill the condition
$\frac{g}{L}\rightarrow 1$.
For the time evolution of the ratio, we obtain

$(\frac{g}{L})_{t}=\frac{g_{t}L-L_{t}g}{L^{2}}.$ (11)

We observe first that the time derivative of the local length is obtained
using the Frenet-Serret formulas and is given by

$\begin{split}g_{t}&=\lvert\textbf{x}_{u}\rvert_{t}=\frac{\textbf{x}_{u}}{\lvert\textbf{x}_{u}\rvert}\cdot(\textbf{x}_{u})_{t}=\frac{g\textbf{x}_{s}}{g}\cdot(\textbf{x}_{t})_{u}=\textbf{T}\cdot(\beta\textbf{N}+\alpha\textbf{T})_{u}=\textbf{T}\cdot
g(\beta\textbf{N}+\alpha\textbf{T})_{s}=\\\ &=\textbf{T}\cdot
g(\beta_{s}\textbf{N}+\beta\textbf{N}_{s}+\alpha_{s}\textbf{T}+\alpha\textbf{T}_{s})=\textbf{T}\cdot
g(\beta_{s}\textbf{N}+k\beta\textbf{T}+\alpha_{s}\textbf{T}-k\alpha\textbf{N})=\\\
&=gk\beta+g\alpha_{s}=gk\beta+\alpha_{u}.\end{split}$ (12)

To get the formula for the global length, we have to integrate
$g=\lvert\textbf{x}_{u}\rvert$ over the curve $\Gamma$

$L=\int_{\Gamma}gdu=\int_{0}^{1}gdu,$

and then

$L_{t}=\int_{0}^{1}g_{t}du=\int_{0}^{1}gk\beta
du+\int_{0}^{1}\alpha_{u}du=\int_{0}^{1}k\beta ds+\alpha(1)-\alpha(0).$ (13)

Now, since we fixed the endpoints of the open curve, we get
$\alpha(1)=\alpha(0)=0$, and equation (13) reduces to

$L_{t}=\int_{0}^{1}k\beta ds=L\langle k\beta\rangle_{\Gamma},$ (14)

where $\langle k\beta\rangle_{\Gamma}=\frac{1}{L}\int_{\Gamma}k\beta ds$.
Substituting (12) and (14) in (11) we obtain

$(\frac{g}{L})_{t}=\frac{(gk\beta+\alpha_{u})L-Lg\langle
k\beta\rangle_{\Gamma}}{L^{2}}=\frac{g}{L}(k\beta+\alpha_{s}-\langle
k\beta\rangle_{\Gamma}).$ (15)

If we now impose that $(\frac{g}{L})_{t}=\omega(1-\frac{g}{L})$, where
$\omega$ is a parameter that determines the speed of convergence of the ratio,
we get the desired condition $\frac{g}{L}\rightarrow 1$. Finally, we obtain
the formula for the tangential velocity

$\alpha_{s}=\langle k\beta\rangle_{\Gamma}-k\beta+\omega(\frac{L}{g}-1).$ (16)

###  2.2 Numerical Discretization

Let us consider the following form of the intrinsic PDE (10)

$\textbf{x}_{t}-\alpha\textbf{x}_{s}=\delta\textbf{x}_{ss}+\lambda
w\textbf{x}_{s}^{\perp},$ (17)

where $\alpha$ and $w$ are given by (16) and (5) respectively. First, we
perform the spatial discretization based on the flowing finite volume method
[7, 9], then we discuss the semi-implicit time discretization, which is
implicit in the intrinsic diffusion term and uses the inflow-implicit/outflow-
explicit strategy for intrinsic advection term [6].

Figure 1: Open curve discretization. Finite volume
$\textbf{p}_{i}=[\textbf{x}_{i-\frac{1}{2}},\textbf{x}_{i+\frac{1}{2}}]$ is
highlighted in blue.

Consider an open curve discretized into $n+2$ grid points. We fixed the two
points $\textbf{x}_{0}$ and $\textbf{x}_{n+1}$ at the endpoints of the curve.
Consequently, for these two points, the velocity is zero.
We define $h_{i}=\lvert\textbf{x}_{i}-\textbf{x}_{i-1}\rvert$ and consider the
finite volume
$\textbf{p}_{i}=[\textbf{x}_{i-\frac{1}{2}},\textbf{x}_{i+\frac{1}{2}}]$ where
$\textbf{x}_{i+\frac{1}{2}}$ represents the middle point between
$\textbf{x}_{i}$ and $\textbf{x}_{i+1}$, i.e.,

$\textbf{x}_{i+\frac{1}{2}}=\frac{\textbf{x}_{i}+\textbf{x}_{i+1}}{2}.$

Notice that the velocities $\alpha_{i},~{}\beta_{i}$ will be referred to the
points $\textbf{x}_{i},i=0,...,n+1$ with $\alpha_{0}=0$, $\beta_{0}=0$ and
$\alpha_{n+1}=0$, $\beta_{n+1}=0$. On the other side, the curvature $k_{i}$
will be referred to the elements $\textbf{h}_{i},~{}i=1,...,n+1$.
Integrating (17) over the finite volume $\textbf{p}_{i}$, see Fig. 1, we get

$\int_{\textbf{x}_{i-\frac{1}{2}}}^{\textbf{x}_{i+\frac{1}{2}}}\textbf{x}_{t}ds-\int_{\textbf{x}_{i-\frac{1}{2}}}^{\textbf{x}_{i+\frac{1}{2}}}\alpha\textbf{x}_{s}ds=\delta\int_{\textbf{x}_{i-\frac{1}{2}}}^{\textbf{x}_{i+\frac{1}{2}}}\textbf{x}_{ss}ds+\lambda\int_{\textbf{x}_{i-\frac{1}{2}}}^{\textbf{x}_{i+\frac{1}{2}}}w\textbf{x}_{s}^{\perp}ds,$
(18)

where $\delta,~{}\lambda$ are constants and the values $\alpha,w$ are
considered to be constant over the finite volume $\textbf{p}_{i}$ and will be
indicated as $\alpha_{i},~{}w_{i}$.
The length of the interval $\textbf{p}_{i}$ is $\frac{h_{i}+h_{i+1}}{2}$
consequently, using the Newton-Leibniz formula we get the following
approximation of (18) for $i=1,...,n$

$\frac{h_{i}+h_{i+1}}{2}(\textbf{x}_{i})_{t}-\alpha_{i}[\textbf{x}]_{\textbf{x}_{i-\frac{1}{2}}}^{\textbf{x}_{i+\frac{1}{2}}}=\delta[\textbf{x}_{s}]_{\textbf{x}_{i-\frac{1}{2}}}^{\textbf{x}_{i+\frac{1}{2}}}+\lambda
w_{i}([\textbf{x}]_{\textbf{x}_{i-\frac{1}{2}}}^{\textbf{x}_{i+\frac{1}{2}}})^{\perp}.$
(19)

First of all, notice that

$\displaystyle-\alpha_{i}[\textbf{x}]_{\textbf{x}_{i-\frac{1}{2}}}^{\textbf{x}_{i+\frac{1}{2}}}$
$\displaystyle=-\alpha_{i}(\textbf{x}_{i+\frac{1}{2}}-\textbf{x}_{i-\frac{1}{2}})=-\alpha_{i}(\frac{\textbf{x}_{i}+\textbf{x}_{i+1}}{2}-\frac{\textbf{x}_{i}+\textbf{x}_{i-1}}{2})=$
(20)
$\displaystyle=-\alpha_{i}(\frac{\textbf{x}_{i+1}-\textbf{x}_{i-1}}{2})=\frac{\alpha_{i}}{2}(\textbf{x}_{i}-\textbf{x}_{i+1})-\frac{\alpha_{i}}{2}(\textbf{x}_{i}-\textbf{x}_{i-1}).$
(21)

Then, approximating $\textbf{x}_{s}$ by a finite difference, we get

$[\textbf{x}_{s}]_{\textbf{x}_{i-\frac{1}{2}}}^{\textbf{x}_{i+\frac{1}{2}}}=(\frac{\textbf{x}_{i+1}-\textbf{x}_{i}}{h_{i+1}}-\frac{\textbf{x}_{i}-\textbf{x}_{i-1}}{h_{i}}).$
(22)

Combining (21), (22), and (19), we obtain

$\begin{split}&\frac{h_{i}+h_{i+1}}{2}(\textbf{x}_{i})_{t}+\frac{\alpha_{i}}{2}(\textbf{x}_{i}-\textbf{x}_{i+1})-\frac{\alpha_{i}}{2}(\textbf{x}_{i}-\textbf{x}_{i-1})=\\\
&=\delta_{i}(\frac{\textbf{x}_{i+1}-\textbf{x}_{i}}{h_{i+1}}-\frac{\textbf{x}_{i}-\textbf{x}_{i-1}}{h_{i}})+\lambda
w_{i}(\frac{\textbf{x}_{i+1}-\textbf{x}_{i-1}}{2})^{\perp}.\end{split}$ (23)

The basic idea of the inflow-implicit/outflow-explicit strategy is that we
treat outflow from a cell explicitly while inflow implicitly [6]. If
$\alpha_{i}<0$, i.e., $(-\alpha)$ in (17) is positive in the finite volume
$\textbf{p}_{i}$, then there is an inflow into the finite volume through its
boundary point $\textbf{x}_{i-\frac{1}{2}}$ and there is an outflow through
the other boundary point $\textbf{x}_{i+\frac{1}{2}}$. On the other hand, if
$\alpha_{i}>0$, there is outflow through $\textbf{x}_{i-\frac{1}{2}}$ and
inflow through $\textbf{x}_{i+\frac{1}{2}}$. We therefore define

$\begin{split}b_{i-\frac{1}{2}}^{in}=max(-\alpha_{i},0),b_{i-\frac{1}{2}}^{out}=min(-\alpha_{i},0),\\\
b_{i+\frac{1}{2}}^{in}=max(\alpha_{i},0),b_{i+\frac{1}{2}}^{out}=min(\alpha_{i},0),\end{split}$
(24)

and rewrite equation (23) as

$\begin{split}&\frac{h_{i}+h_{i+1}}{2}(\textbf{x}_{i})_{t}+\frac{1}{2}(b_{i+\frac{1}{2}}^{in}+b_{i+\frac{1}{2}}^{out})(\textbf{x}_{i}-\textbf{x}_{i+1})+\frac{1}{2}(b_{i-\frac{1}{2}}^{in}+b_{i-\frac{1}{2}}^{out})(\textbf{x}_{i}-\textbf{x}_{i-1})\\\
&=\delta(\frac{\textbf{x}_{i+1}-\textbf{x}_{i}}{h_{i+1}}-\frac{\textbf{x}_{i}-\textbf{x}_{i-1}}{h_{i}})+\lambda
w_{i}(\frac{\textbf{x}_{i+1}-\textbf{x}_{i-1}}{2})^{\perp}.\end{split}$ (25)

Let us approximate the time derivative by the finite difference
$\textbf{x}_{t}=\frac{\textbf{x}_{i}^{m+1}-\textbf{x}_{i}^{m}}{\tau}$: here
$m$ is the time step index and $\tau$ is the length of the discrete time step.
Remark that we take the unknowns in the inflow part of the advection term
implicitly while we take the ones in the outflow part explicitly. Moreover, we
take the diffusion term implicitly and the attracting term explicitly.
Finally, we obtain the fully discrete scheme

$\begin{split}&\textbf{x}_{i-1}^{m+1}(-\frac{\delta}{h_{i}^{m}}-\frac{b_{i-\frac{1}{2}}^{in^{m}}}{2})+\textbf{x}_{i+1}^{m+1}(-\frac{\delta}{h_{i+1}^{m}}-\frac{b_{i+\frac{1}{2}}^{in^{m}}}{2})+\\\
&\textbf{x}_{i}^{m+1}(\frac{h_{i+1}^{m}+h_{i}^{m}}{2\tau}+\frac{\delta}{h_{i}^{m}}+\frac{\delta}{h_{i+1}^{m}}+\frac{b_{i-\frac{1}{2}}^{in^{m}}}{2}+\frac{b_{i+\frac{1}{2}}^{in^{m}}}{2})=\textbf{x}_{i}^{m}\frac{h_{i+1}^{m}+h_{i}^{m}}{2\tau}\\\
&-\frac{b_{i+\frac{1}{2}}^{out^{m}}}{2}(\textbf{x}_{i}^{m}-\textbf{x}_{i+1}^{m})-\frac{b_{i-\frac{1}{2}}^{out^{m}}}{2}(\textbf{x}_{i}^{m}-\textbf{x}_{i-1}^{m})+\lambda
w_{i}^{m}(\frac{\textbf{x}_{i+1}^{m}-\textbf{x}_{i-1}^{m}}{2})^{\perp},\end{split}$
(26)

for $i=1,...,n$ where $n$ is the number of unknown grid points.
The system (26) is represented by a strictly diagonally dominant matrix, then
it is always solvable by the classical Thomas algorithm without any
restriction on the time step $\tau$.
In the numerical scheme (26), there are two parameters $\alpha_{i}^{m}$ and
$w_{i}^{m}$ given by (16) and (5) for which we have not yet defined the
discretization.
Let us first consider the discretization of $w_{i}^{m}$ referred to the curve
grid point $\textbf{x}_{i}^{m}$. Remark that the shortest distance of a point
to a curve is in the normal direction to the curve.
The algorithm is organized as follows: consider a point $\textbf{x}_{i}^{m}$
on the evolving curve. Set $D_{min}=D$ (where $D$ is a real number to be
chosen reasonably) and for $j=1,...,n+1$ repeat:

  1. 1.

consider the line $r$ passing through the points $\textbf{x}_{j-1}^{0}$ and
$\textbf{x}_{j}^{0}$ of the original curve,

  2. 2.

find the line $s$ perpendicular to $r$ and passing through
$\textbf{x}_{i}^{m}$,

  3. 3.

find the point q of intersection between the line $r$ and the line $s$,

  4. 4.

find $d=d(\textbf{q},\textbf{x}_{i}^{m})$, where $d$ is the Euclidean distance
of q and $\textbf{x}_{i}^{m}$,

  5. 5.

if $\textbf{q}\in[\textbf{x}_{j-1}^{0},\textbf{x}_{j}^{0}]$ and $d<D_{min}$
then

$\begin{split}&\textbf{x}^{0}-\textbf{x}^{m}_{i}=\textbf{q}-\textbf{x}_{i}^{m}\\\
&D_{min}=\lvert\textbf{x}^{0}-\textbf{x}^{m}_{i}\rvert\end{split}$ (27)

Due to the discretization of the curve, it may happen that the algorithm does
not find any point q. In this case, for every $j=0,...,n+1$, we check the
Euclidean distance $d(\textbf{x}_{j}^{0},\textbf{x}_{i}^{m})$ and we select
the minimum $j=j_{min}$.Then,

$\textbf{x}^{0}-\textbf{x}^{m}_{i}=\textbf{x}_{j_{min}}^{0}-\textbf{x}_{i}^{m}.$

This approach is natural since we are trying to find the minimum distance
between a point and a segment: it is either the minimum distance found in (27)
or the minimum distance between the point and one of the endpoints of the
segment.
Therefore, we get

$w_{i}^{m}=(\textbf{x}^{0}-\textbf{x}^{m}_{i})\cdot\textbf{N}_{i}^{m},$ (28)

where
$\textbf{N}_{i}^{m}=(\frac{\textbf{x}_{i+1}^{m}-\textbf{x}_{i-1}^{m}}{h_{i}^{m}+h_{i+1}^{m}})^{\perp}$.

Consider now the discretization of $\alpha_{i}^{m}$. Remark that, since we
fixed the first and the last point, we set $\alpha_{0}=0$ and
$\alpha_{n+1}=0$. We get $\alpha_{i}^{m}$ for $i=1,...,n$ by

$\alpha_{i}^{m}=\alpha_{i-1}^{m}+h_{i}^{m}\langle
k\beta\rangle_{\Gamma}^{m}-h_{i}^{m}k_{i}^{m}\beta_{i}^{m}+\omega(\frac{L^{m}}{n+1}-h_{i}^{m}).$
(29)

The curvature $k_{i}^{m}$, the normal velocity $\beta_{i}^{m}$, the mean value
$\langle k\beta\rangle_{\Gamma}^{m}$, and the total length $L^{m}$ are given
by the following formulas [1, 9, 10]

$\begin{split}&k_{i}^{m}=sgn(\textbf{h}_{i-1}^{m}\wedge\textbf{h}_{i+1}^{m})\frac{1}{2h_{i}^{m}}arccos(\frac{\textbf{h}_{i-1}^{m}\cdot\textbf{h}_{i+1}^{m}}{h_{i-1}^{m}h_{i+1}^{m}}),\\\
&\beta_{i}^{m}=-\delta k_{i}+\lambda w_{i},\\\ &\langle
k\beta\rangle_{\Gamma}^{m}=\frac{1}{L^{m}}\sum_{l=1}^{n+1}h_{l}^{m}k_{l}^{m}\beta_{i}^{m},\\\
&L^{m}=\sum_{l=1}^{n+1}h_{l}^{m},\end{split}$ (30)

where $\textbf{h}_{i}^{m}=\textbf{x}_{i}^{m}-\textbf{x}_{i-1}^{m}$,
$\lvert\textbf{h}_{i}^{m}\rvert=h_{i}^{m}$. For the first and last elements,
since we don’t have the values of $h_{i-1}^{m}$, respectively $h_{i+1}^{m}$,
we set $k_{1}^{m}=k_{2}^{m}$ and $k_{n+1}^{m}=k_{n}^{m}$.

##  3 Numerical experiments

###  3.1 Smoothed trajectories

First, we considered the simple curve shown in Figs. 2-4 to see the influence
of the term which attracts the curve to the initial condition. In Fig. 2, we
considered the evolution of the curve driven by the PDE (1) with $\lambda=0$
and $\lambda=1$, $\delta=0.001$, $\omega=1$, and $\tau=0.001$ for $10000$ time
steps. The red line represents the initial condition while the blue lines are
the results of the smoothing plotted every $1000$ time step. The curve evolved
only by the curvature influence goes to the straight line, while if we
consider $\lambda=1$, the curve is attracted to the initial condition and
stays closer to it. In Fig. 3 we considered the same initial curve as in Fig.
2 and we plotted the results of the smoothing after $400$ time steps for
$\lambda=0$, $\lambda=1$, $\lambda=5$ and $\lambda=10$, $\delta=0.001$,
$\omega=1$, and $\tau=0.001$: we observe that considering bigger $\lambda$,
the smoothed curve stays closer to the original one. To better understand the
influence of the attracting term in Fig. 4 we plotted the vectors driving one
of the discrete points of the evolving curve: the blue arrow is the vector
$(\textbf{x}_{0}-\textbf{x})$, the red arrow is its projection in the normal
direction to the evolving curve, and the green one is the diffusion term. We
observe that the attracting vector always points towards the original curve:
it is clear that, if we consider bigger $\lambda$, the result will stay closer
to the initial condition.

Figure 2: Evolution of the initial curve (red) with $\lambda=1$ (top) and
$\lambda=0$ (bottom), $\delta=0.001$, $\omega=1$, and $\tau=0.001$. The blue
curves are the evolved curves plotted every $1000$ time step, from time step
$1000$ to $10000$ Figure 3: Comparison of the evolution of the initial curve
(red) for different values of $\lambda$ after $400$ time steps. Results are
shown for $\lambda=0$ (blue), $\lambda=1$ (light blue), $\lambda=5$ (pink) and
$\lambda=10$ (green), $\delta=0.001$, $\omega=1$, and $\tau=0.001$. Figure 4:
Plotting of vectors driving one of the curve grid points during the evolution.
$-\delta k\textbf{N}$ (green), $(\textbf{x}_{0}-\textbf{x})$(blue), and its
projection in normal direction to the curve
$[(\textbf{x}_{0}-\textbf{x})\cdot\textbf{N}]\textbf{N}$(red).

Then, we focused on real data of macrophage motion. In the following
experiments, we applied the proposed model to macrophage trajectories. The
piecewise linear initial condition was found connecting the cell center of the
macrophage for every frame of the $2$D$+$time data [12]. First, we added grid
points in the original trajectory so that we had almost uniformly distributed
points. For every trajectory, the parameters were chosen as follows:
$\lambda=1$, $\delta=0.005$, $\omega=1$, and $\tau=0.0001$. Figs. 5-9 show the
result of the smoothing for two different trajectories. The red lines in Figs.
5 and 8 are the original trajectories for two different macrophages. As one
can observe, due to the random motion of the cells, these trajectories usually
have noisy parts: Figs 6, 7, and 9 show the details of these parts. The blue
lines are the smoothed trajectories after different numbers of time steps: the
evolving curve is kept close to the original trajectory by the attracting term
but we also observe the smoothing of the noisy parts due to the diffusion
term. In Figs. 6 and 7, we show the result of the smoothing until the random
part of the trajectory is completely smoothed: the black dot represents the
last point of the trajectory.

As a stopping criterion, we considered the mean Hausdorff distance between two
discrete curves. Consider the discrete sets of points
$\mathcal{A}=\\{a_{0},...,a_{n+1}\\}$, $\mathcal{B}=\\{b_{0},...,b_{n+1}\\}$:
we will indicate the elements of the discrete curves by
$\textbf{a}_{i}=a_{i}-a_{i-1}$ and $\textbf{b}_{i}=b_{i}-b_{i-1}$,
$i=1,...,n+1$ respectively. Let us indicate by
$A=\\{\textbf{a}_{1},...,\textbf{a}_{n+1}\\}$ and
$B=\\{\textbf{b}_{1},...,\textbf{b}_{n+1}\\}$ the sets of elements. The mean
Hausdorff distance $\bar{d}_{H}(\mathcal{A},\mathcal{B})$ is then defined as

$\begin{split}\bar{\delta}_{H}(\mathcal{A},\mathcal{B})&=\frac{1}{n}\sum_{i=1}^{n}\min_{\textbf{b}\in
B}d(a_{i},\textbf{b}),\\\
\bar{\delta}_{H}(\mathcal{B},\mathcal{A})&=\frac{1}{n}\sum_{i=1}^{n}\min_{\textbf{a}\in
A}d(b_{i},\textbf{a}),\\\
\bar{d}_{H}(\mathcal{A},\mathcal{B})&=\frac{\bar{\delta}_{H}(\mathcal{A},\mathcal{B})+\bar{\delta}_{H}(\mathcal{B},\mathcal{A})}{2},\\\
\end{split}$ (31)

where $d(a_{i},\textbf{b})$ is the distance between a point
$a_{i}\in\mathcal{A}$ and an element $\textbf{b}\in B$ defined in section 2.2.
Since in our model the endpoints $a_{0},~{}a_{n+1},~{}b_{0},~{}b_{n+1}$ are
fixed, we do not consider them in the computation of the mean Hausdorff
distance. In order not to slow down the computation too much this distance was
calculated only every $p$ number of time steps; we chose the tolerance
$\epsilon=0.000065$ and stopped iterating when
$d_{H}(\textbf{x}_{n},\textbf{x}_{n+p})<\epsilon$. For the trajectory in Fig.
5, the number of time steps needed to stop the evolution was $580$, while for
the one in Fig. 8, it was $340$.

Figure 5: Original trajectory (red line) and smoothed trajectory (blue line).
Results (from top to bottom) after $50$, $200$, $400$ and $580$ time steps.

Figure 6: Detail of noisy part of trajectory in Fig. 5. Original trajectory
(red line) and smoothed trajectory (blue line). Results (from top to bottom)
after $50$, $200$, $400$ and $580$ time steps.

Figure 7: Detail of noisy part of trajectory in Fig. 5. Original trajectory
(red line) and smoothed trajectory (blue line). Results (from top to bottom)
after $1000$, $1500$, $2500$, and $5000$ time steps.

Figure 8: Original trajectory (red line) and smoothed trajectory (blue line).
Results (from top to bottom) after $50$, $100$, $200$ and $340$ time steps.

Figure 9: Detail of noisy part of trajectory in Fig. 8. Original trajectory
(red line) and smoothed trajectory (blue line). Results (from top to bottom)
after $50$, $100$, $200$ and $340$ time steps.

###  3.2 Velocities on smoothed curves

Our final goal was to compute the velocity for every grid point on the evolved
curve. In the $2$D$+$time data, the real-time between each frame was $1$ time
unit between the frames, which is equal to 2 minutes. That means, in the
original piecewise linear trajectory, the time between two endpoints of a
segment was $1$ time unit (see Fig. 10). First, when we added the grid points
inside the initial segments, we did it in such a way that the new added points
were uniformly distributed inside the segment. Consequently, we calculated the
original velocity for the grid points as the original length of the segment
divided by the time interval $\Delta_{t}$, that in our case was $1$ (time
unit), and consider it constant inside the segment. Therefore, in the
beginning, all the added grid points inside a segment have the same velocity.
We assume this property holds also for the new velocities of the grid points
on the evolved curve. Our goal was to derive the velocity for every grid point
on the evolved curve but we couldn’t do it directly because, due to the
stability of the numerical computations, we added the tangential velocity in
our model and this caused an asymptotically uniform redistribution of the
points. If the trajectory is not too complicated, the tangential velocity does
not move the points too much. But, in our case, some of the trajectories had
noisy parts with many points, which were moved a lot during the evolution, so
we were losing almost completely the spatial information.
Then, to get the new velocity estimation, first we studied the evolution of
the length of the segments. For the total length, we have already shown in
section 2.1 that it holds

$L_{t}=\int_{0}^{1}k\beta ds+\int_{0}^{1}\alpha_{s}ds=\int_{0}^{1}k\beta
ds+\alpha(1)-\alpha(0).$ (32)

In the case of the total length, we had fixed the endpoints of the trajectory
so we had $\alpha(1)=\alpha(0)=0$. Similarly, for the evolution of every
segment, we can consider the same formula. We will denote by $j$ the index for
the segments and by $i$ the index for the grid points. Consider the $j$-th
segment: we will indicate by $u_{j-1}$ and $u_{j}$ the first and last endpoint
of the segment parametrization, respectively. Therefore, for the evolution of
the length $L_{j}$ of the $j$-th segment, it holds

$(L_{j})_{t}=\int_{u_{j-1}}^{u_{j}}k\beta
ds+\int_{u_{j-1}}^{u_{j}}\alpha_{s}ds=\int_{u_{j-1}}^{u_{j}}k\beta
ds+\alpha(u_{j})-\alpha(u_{j-1}).$ (33)

Notice that in this case, in general, we have $\alpha(u_{j})\neq 0$ and
$\alpha(u_{j-1})\neq 0$. But, as we have already underlined, the tangential
velocity is used only for the stability of the numerical computations and it
does not influence the shape of the evolving curve. So, to not lose the
spatial information, we will calculate the new length of the segment as if we
have applied to the trajectory the model with $\alpha=0$, i.e.,
$\alpha(u_{j})=\alpha(u_{j-1})=0$. Consequently, we consider the following
model for the evolution of the $j$-th segment

$(L_{j})_{t}=\int_{u_{j-1}}^{u_{j}}k\beta ds.$ (34)

Let’s now consider the discretization: let $m$ be the time step index and
$\tau$ the length of the discrete time step. We approximate the time
derivative by the finite difference and consider the formulas for $k,~{}\beta$
defined in (30). To relate the index $j$ of the segments to the index $i$ of
the grid points we will indicate by $\mathcal{I}(u_{j-1})$,
$\mathcal{I}(u_{j})$ the corresponding $i$ index of the endpoint
$\textbf{x}(u_{j-1})$, respectively $\textbf{x}(u_{j})$ (see Fig. 10). Then,
we obtain

$\frac{L_{j}^{m+1}-L_{j}^{m}}{\tau}=\sum_{i=\mathcal{I}(u_{j-1})+1}^{\mathcal{I}(u_{j})}h_{i}^{m}k_{i}^{m}\beta_{i}^{m}.$
(35)

Figure 10: Visualization of curve discretization. The green dots represent the
grid points and the red dots represent the endpoints of the segments.

Notice that we consider the $i$ index starting from $\mathcal{I}(u_{j-1})+1$:
that is because $h_{\mathcal{I}(u_{j-1})+1}$ is defined as
$h_{\mathcal{I}(u_{j-1})+1}=\lvert\textbf{x}_{\mathcal{I}(u_{j-1})+1}-\textbf{x}_{\mathcal{I}(u_{j-1})}\rvert$,
which is indeed the first element $h_{i}$ of the segment with endpoints
$\textbf{x}(u_{j-1})$ and $\textbf{x}(u_{j})$. We obtain the fully discrete
approximation of equation (34) for the new length of the $j$-th segment

$L_{j}^{m+1}=L_{j}^{m}+\tau\sum_{i=\mathcal{I}(u_{j-1})+1}^{\mathcal{I}(u_{j})}h_{i}^{m}k_{i}^{m}\beta_{i}^{m}.$
(36)

For the first and last element of each segment, we use a different computation
of $k$ and $\beta$ that does not consider informations from the previous and
next segment.
If we consider the evolution of the curve driven only by the curvature
influence, then $\beta_{i}=-\delta k_{i}$ and we obtain

$\sum_{i=\mathcal{I}(u_{j-1})+1}^{\mathcal{I}(u_{j})}h_{i}^{m}k_{i}^{m}\beta_{i}^{m}=-\sum_{i=\mathcal{I}(u_{j-1})+1}^{\mathcal{I}(u_{j})}h_{i}^{m}\delta(k_{i}^{m})^{2}<0.$

Therefore $L_{j}^{m+1}<L_{j}^{m}$ for every segment $j$ and for every time
step index $m$. Experimentally, this property holds also if $\beta_{i}=-\delta
k_{i}+\lambda w_{i}$. If the segment is in a part of the trajectory which has
high curvature (for example in the noisy parts where the cell shows random
motion), the length of that segment will decrease faster than the length of a
segment in a part of the trajectory with small curvature (for example in the
parts where the cell shows directional motion). It is then clear that some of
the segments may eventually disappear and we may obtain $L_{j}^{m+1}=0$.

Let us first focus on the simple case where the length of one segment is
decreasing but the segment is not disappearing. In the first time step,
$\mathcal{I}(u_{j-1})$ and $\mathcal{I}(u_{j})$ are loaded from the original
trajectory and are the indexes of the endpoints of the original segments.
Later, the length of every segment changes so, after applying the formula in
(36), we need to find the new $i$ indexes for $\mathcal{I}(u_{j-1})$ and
$\mathcal{I}(u_{j})$. If we indicate by

$\begin{split}L_{s}^{m+1}&=\sum_{j=1}^{M-1}L_{j}^{m+1}\\\ \end{split}$ (37)

where $M-1$ is the number of original segments. Ideally,
$L_{s}^{m+1}=L^{m+1}$, where $L^{m+1}$ is defined in (30). However, in the
numerical computations, in later time steps, it holds $L_{s}^{m+1}\neq
L^{m+1}$. We define the ratio of the $j$-th segment as

$r_{j}^{m+1}=\frac{L_{j}^{m+1}}{L_{s}^{m+1}},$ (38)

and the discrete length of the segment is defined as

$L^{d}_{j}=r_{j}^{m+1}L^{m+1},$ (39)

so that $\sum_{j=1}^{M-1}L^{d}_{j}=L^{m+1}$. Once we obtain the new lengths
$L^{d}_{j}$, to get the new $i$ indexes $\mathcal{I}(u_{j-1})$ and
$\mathcal{I}(u_{j})$ for the endpoints of the segments we apply the following
algorithm. Set $i=1$, $j=1$, $S=L^{d}_{j}$, $S_{e}=0$ and repeat until $j=M-1$
and $i=n+1$:

  * •

$S_{e}=S_{e}+h_{i}$

  * •

if $S_{e}\geq S$

    * –

if $S_{e}>S$ then $\mathcal{I}(u_{j})=i-1$

    * –

if $S_{e}=S$ then $\mathcal{I}(u_{j})=i$

    * –

update $S=S+L^{d}_{j+1}$.

    * –

$j=j+1$

  * •

$i=i+1$

Notice that we focused only on the index of the last endpoint of each segment
$\mathcal{I}(u_{j})$: that is because the last endpoint of the segment $j$ is
the first endpoint of the segment $j+1$.

Let us now focus on the disappearing segments. We add the following condition:
if $L_{j}^{m+1}<h_{min}$, then $L_{j}^{m+1}=0$, where

$h_{min}=\min_{1\leq i\leq n+1}h_{i}.$

If we obtain $L_{j}^{m+1}=0$ that means the first and the last endpoints of
that segment at time step index $m+1$ became the same point, then
$\mathcal{I}(u_{j})=\mathcal{I}(u_{j-1})$. So we modify the algorithm for
finding the new indexes $i$ for $\mathcal{I}(u_{j})$ adding the following
condition

  * •

if $r_{j}^{m+1}=\frac{L_{j}^{m+1}}{L_{s}}=0$ then
$\mathcal{I}(u_{j})=\mathcal{I}(u_{j-1})$.

This condition will cause that, when we apply formula (36) for every time step
index $m+k$, with $k>0$, we will obtain

$L_{j}^{m+k}=L_{j}^{m+1}=0.$

To get the new velocity estimation consider $(\Delta_{t})_{j}$ the time
between the two endpoints $\textbf{x}(u_{j-1})$ and $\textbf{x}(u_{j})$ of the
segment $j$. In the beginning, every segment has the same value of
$(\Delta_{t})_{j}$; later, if a segment does not disappear, we consider the
same $(\Delta_{t})_{j}$ as in the beginning. On the other hand, if a segment
$j$ disappears, we add $\frac{(\Delta_{t})_{j}}{2}$ to the first previous not
disappeared segment,

$(\Delta_{t})_{j-1}=(\Delta_{t})_{j-1}+\frac{(\Delta_{t})_{j}}{2},$

and $\frac{(\Delta_{t})_{j}}{2}$ to the first following not disappeared
segment

$(\Delta_{t})_{j+1}=(\Delta_{t})_{j+1}+\frac{(\Delta_{t})_{j}}{2}.$

To find the new velocity estimation for the grid points on the evolved curve,
we apply the following algorithm. Set $j=1$, $i=1$ and repeat until $j=M-1$:

  * •

If $r_{j}^{m+1}>0$

$\lvert\textbf{v}\rvert=\frac{L^{d}_{j}}{(\Delta_{t})_{j}}.$ (40)

  * •

for $i=\mathcal{I}(u_{j-1}),...,\mathcal{I}(u_{j})$

$\textbf{v}_{i}^{m+1}=\lvert\textbf{v}\rvert\frac{(\textbf{x}_{i}^{m+1}-\textbf{x}_{i-1}^{m+1})}{\lvert\textbf{x}_{i}^{m+1}-\textbf{x}_{i-1}^{m+1}\rvert}.$
(41)

If the segment does not disappear, we consider the grid points inside that
segment to have constant velocity equal to the new length of the segment
divided by the time interval $(\Delta_{t})_{j}$. On the other hand, if the
$j$-th segment disappears, the time $(\Delta_{t})_{j}$ is distributed between
the previous and following segments that have not disappeared.
For the numerical experiments, to first check our method, we set as initial
condition a semi-ellipse and we considered the evolution of the curve with
$\lambda=0$, $\delta=0.05$, $\omega=1$, and $\tau=0.001$. We calculated the
velocities using the algorithm described above: the result is shown in Fig. 11
after $1000$ time steps. Once we got the velocities for every grid point, we
normalized the length of the vectors to obtain better visualization: the color
represents the norm. Then, if the arrow is yellow it means in that grid point
the velocity is high. On the other side, if the arrow is blue or light blue,
at that point the velocity is slower. As we expected, the length of the
segment in the center of the discretized semi-ellipse decreases faster than
the others. Indeed, that is the part of the ellipse with the highest
curvature.

Figure 11: Top: initial trajectory and plotting of the velocity vectors.
Bottom: initial trajectory (red line), evolved curve, and plotting of the
velocity vectors. The vectors are normalized for better visualization, color
represents the norm.

Then, we focused on the real trajectories of the macrophages. We considered
the evolution of the trajectories with $\lambda=1$, $\delta=0.005$,
$\omega=1$, and $\tau=0.0001$. Results for three different trajectories are
shown in Figs. 12-14. We considered the evolved curve obtained using the mean
Hausdorff distance as a stopping criterion as described in the last paragraph
of section 3.1. We observe that the velocity is slowed down in the random
parts of the original trajectory that are smoothed during the evolution.

Figure 12: Top: initial trajectory and plotting of the velocity vectors.
Bottom: initial trajectory (red line), evolved curve, and plotting of the
velocity vectors. The vectors are normalized for better visualization, color
represents the norm.

Figure 13: Top: initial trajectory and plotting of the velocity vectors.
Bottom: initial trajectory (red line), evolved curve, and plotting of the
velocity vectors. The vectors are normalized for better visualization, color
represents the norm.

Figure 14: Top: initial trajectory and plotting of the velocity vectors.
Bottom: initial trajectory (red line), evolved curve, and plotting of the
velocity vectors. The vectors are normalized for better visualization, color
represents the norm.

##  4 Acknowledgement

This work has received funding from the European Union’s Horizon 2020 research
and innovation programme under the Marie Sklodowska-Curie grant agreement No
955576 and by the grants APVV-19-0460, VEGA 1/0436/20.

## References

  * [1] M.Ambroz, M.Kollár, K.Mikula, Semi-implicit scheme for semiautomatic segmentation in NaturaSat software, ALGORITMY 2020- 21st Conference on Scientific Computing, Vysoké Tatry - Podbanské, Slovakia, September 10 - 15, 2020, Proceedings of contributed papers, Editors: P.Frolkovič, K.Mikula and D.Ševčovič, Vydavateĺstvo SPEKTRUM STU, 2020, ISBN: 978-80-227-5032-5, pp. 171-180.
  * [2] K.M.Henry, L.Pase, C.F.Ramos-Lopez, G.J.Lieschke, S.A.Renshaw, C.C.Reyes-Aldasoro, PhagoSight: an open-source MATLAB® package for the analysis of fluorescent neutrophil and macrophage migration in a zebrafish model, 2013, PLoS ONE, 8(8): e72636.
  * [3] T.Y.Hou, J.S.Lowengrub, M.J.Shelley, Removing the stiffness from interfacial flows with surface tension, Journal of Computational Physics, 1994, Elsevier, 114(2): pp. 312-338.
  * [4] V.Kadirkamanathan,S.R.Anderson,S.A.Billings,X.Zhang,G.R.Holmes, C.C.Reyes-Aldasoro, P.M.Elks,S.A.Renshaw, The neutrophil’s eye-view: Inference and visualisation of the chemoattractant field driving cell chemotaxis in vivo, 2012, PloS ONE, 7(4): e35182.
  * [5] M.Kimura, Numerical analysis of moving boundary problems using the boundary tracking method, Japan journal of industrial and applied mathematics, 1997, Springer, 14(3): pp. 373-398.
  * [6] K.Mikula, M.Ohlberger, Inflow-Implicit/Outflow-Explicit Scheme for Solving Advection Equations, Finite Volumes in Complex Applications VI, Problems & Perspectives, Eds. J.Fořt et al. (Proceedings of the Sixth International Conference on Finite Volumes in Complex Applications, Prague, June 6-10, 2011), 2011, Springer Verlag, pp. 683-692.
  * [7] K.Mikula, D.Ševčovič, Evolution of plane curves driven by a nonlinear function of curvature and anisotropy, SIAM Journal on Applied Mathematics, 61(5), 2001, pp. 1473-1501.
  * [8] K.Mikula, D.Ševčovič, A direct method for solving an anisotropic mean curvature flow of plane curves with an external force, Mathematical Methods in the Applied Sciences, 27(13), 2004: pp. 1545-1565.
  * [9] K.Mikula, D.Ševčovič, M.Balazovjech, A simple, fast and stabilized flowing finite volume method for solving general curve evolution equations, Communications in Computational Physics, 7(1), 2010: pp. 195-211.
  * [10] K.Mikula, J.Urbán, M.Kollár, M.Ambroz, I. Jarolímek, J.Šibík, M.Šibíková, An automated segmentation of NATURA 2000 habitats from Sentinel-2 optical data, Discrete and Continuous Dynamical Systems - Series S, 2021, 14(3): 1017-1032.
  * [11] Seol Ah Park, T.Sipka, Z.Krivá, M.Ambroz, M.Kollár, B.Kósa, M.Nguyen-Chi, G.Lutfalla, K.Mikula, Macrophage image segmentation by thresholding and subjective surface method, Tatra Mountains Mathematical Publications, Vol. 75, 2020, pp. 103-120.
  * [12] Seol Ah Park, T.Sipka, Z.Krivá, M.Nguyen-Chi, G.Lutfalla, K.Mikula, Segmentation-based tracking of macrophages in 2d+time microscopy movies inside a living animal, Computers in Biology and Medicine, Submitted, 2022.


# HyperSound: Generating Implicit Neural Representations of Audio Signals with
Hypernetworks

Filip Szatkowski1,2∗ Karol J. Piczak4∗ Przemysław Spurek4
Jacek Tabor4 Tomasz Trzciński1,2,3,4
1Warsaw University of Technology 2IDEAS NCBR 3Tooploox
4Faculty of Mathematics and Computer Science, Jagiellonian University
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>

###### Abstract

Implicit neural representations (INRs) are a rapidly growing research field,
which provides alternative ways to represent multimedia signals. Recent
applications of INRs include image super-resolution, compression of high-
dimensional signals, or 3D rendering. However, these solutions usually focus
on visual data, and adapting them to the audio domain is not trivial.
Moreover, it requires a separately trained model for every data sample. To
address this limitation, we propose HyperSound, a meta-learning method
leveraging hypernetworks to produce INRs for audio signals unseen at training
time. We show that our approach can reconstruct sound waves with quality
comparable to other state-of-the-art models.

##  1 Introduction

Implicit neural representations (INRs) are coordinate-based representations of
multimedia signals, where the signal is modeled with a neural network. Such
representations are decoupled from the spatial resolution, so the signal can
be resampled at any arbitrary frequency. At the same time, the memory
requirements for its storage remain constant. The field of INRs is rapidly
growing, and their applications include super-resolution [1, 2], compression
[1, 2] or 3D rendering [3]. However, in the audio domain, the evaluation of
these approaches was so far limited to learning individual INRs of particular
input recordings [1], which is highly ineffective.

Meta-learning methods, such as hypernetworks [4], can generate INRs for any
arbitrary signal with a single model. Hypernetworks learn to generate weights
for smaller target networks, which can serve as INRs. Hypernetworks were
successfully applied in obtaining INRs for images [2, 5], point clouds [2, 6]
and videos [2]. However, creating INRs for high-dimensional and high-variance
data such as audio is difficult, and training hypernetworks for this task can
be unstable.

In this work, we propose a meta-learning approach based on hypernetworks,
where we learn a general recipe for creating INRs of arbitrary audio samples
from outside of the training dataset. To our knowledge, our model is the first
application of hypernetwork-based INRs to the audio domain.

Figure 1: Overview of the HyperSound framework. We use a single hypernetwork
model to produce distinct INRs based on arbitrary audio signals provided as
input.

##  2 Related Works

Implicit neural representations (INRs) are applied across many multimedia
domains. NeRF [3]-based INRs are currently state-of-the-art for image-based
rendering and view synthesis. SIREN [1] shows how to produce good quality INRs
for various signals, such as image, video, and sound signals, using neural
networks with periodic activation functions and dedicated weight
initialization scheme. Modulated Periodic Activations [2] further improve the
quality of SIREN-based INRs for high-resolution signals by introducing
modulation and synthesis sub-networks. However, in the audio domain, the
evaluation of these approaches was so far limited to learning individual INRs
of particular input recordings.

Hypernetworks [4] are a meta-learning framework, where one network
(hypernetwork) generates the weights for another network (target network).
Hypernetworks can be used for a variety of tasks such as model compression
[7], continual learning [8] or generating INRs. In particular, hypernetworks
were used to generate INRs for images [5], shapes [6] and videos [2]. To our
knowledge, our work is the first application of hypernetworks for audio INRs
generation.

First successful attempts at raw waveform processing with deep neural networks
were models such as WaveNet [9] and SampleRNN [10], but their autoregressive
nature makes them slow and prone to accumulation of errors. Later
architectures such as ParallelWaveNet [11], NSynth [12], MelGAN [13] or SING
[14] proposed non-autoregressive architectures for audio generation. Recent
autoencoder-based models such as RAVE [15] or SoundStream [16] are able to
process high-resolution signals in an end-to-end fashion, producing audio of
very good perceptual quality.

##  3 Model overview

Sound waves are traditionally represented digitally as a collection of
amplitude values sampled at regular intervals, which approximates a continuous
real function $x(t)$. Our goal is to obtain a meta-recipe for generating audio
INRs replicating such functions. While creating INRs for particular audio
samples can be quite easily done with gradient descent, finding a general
solution is much harder due to the inherent complexity of audio time series.
Therefore, similar to [5], we model such functions with neural networks $T$
(target networks), parameterized by weights generated by another neural
network $H$ (hypernetwork). Our framework, shown in Fig. 1, can be described
as

$\bm{\theta}_{\bm{x}}=H(\bm{x}),\vspace*{-.3em}$ (1)
$\hat{x}(t)=T(t,{\bm{\theta}}_{\bm{x}}).$ (2)

###  3.1 Hypernetwork architecture

Typical audio recordings contain several thousands of samples, so the
hypernetwork is composed of a convolutional encoder that produces a latent
representation of a lower dimensionality, and fully connected layers that
transform this representation to weights $\bm{\theta}$ of the target network.

We use an encoder based on SoundStream [16], and the fully connected part of
the hypernetwork is composed of six fully-connected layers with biases and ELU
[17] activation, where the last layer produces the flattened weights of the
target network.

###  3.2 Approximating sound waves with neural networks

A target network should be as small as possible to avoid overfitting and
explosion of weights in the hypernetwork, but must be expressive enough to
represent a wide variety of audio recordings. Our target network has one input
and one output, and consists of a positional embedding layer followed by four
fully-connected layers of $256$ neurons with biases and ReLU activation.
Inspired by NeRF [3], we define embedding vectors $\bm{\gamma}$ as

$\bm{\gamma}(t)=\left[\sin(2^{0}\pi t),\,\cos(2^{0}\pi t),\,\sin(2^{1}\pi
t),\,\cos(2^{1}\pi t),\,\ldots,\,\sin({2}^{L-1}\pi t),\,\cos({2}^{L-1}\pi
t)\right],$ (3)

where $t$ is the time coordinate and $L$ denotes the embedding size. We
rescale input coordinates to a range of $[0,1]$ and use $L=16$.

###  3.3 Optimization

We train the hypernetwork in a supervised fashion using backpropagation. To
obtain audio results that are more perceptually pleasant, we use a loss
function that penalizes the reconstruction error both in the time and
frequency domains. Given an original recording $\bm{x}$ and its reconstruction
$\hat{\bm{x}}$ generated with a target network, we compute the loss function
as

$L(\bm{x},\hat{\bm{x}})={\lambda}_{SL1}*L_{SL1}(\bm{x},\hat{\bm{x}})+{\lambda}_{STFT}*L_{STFT}(\bm{x},\hat{\bm{x}}),$
(4)

where $L_{SL1}$ is a smooth L1 loss [18] with $\beta=0.1$, $L_{STFT}$ is a
multi-resolution mel-scale STFT loss introduced in ParallelWaveGAN [19] and
${\lambda}_{SL1}$, ${\lambda}_{STFT}$ are the weights of these two losses. We
use ${\lambda}_{SL1}=1$ and ${\lambda}_{STFT}=1$. For STFT loss, we use $128$
mel bins and FFT sizes of $[512,1024,2048]$ with matching window sizes and an
overlap of $87.5$%.

##  4 Experiments

Ground truthReconstructionGround truthReconstruction Figure 2: Examples of
VCTK validation samples reconstructed with HyperSound.

We test the reconstruction quality of our model on the VCTK dataset
downsampled to $f=$22\,050\text{\,}\mathrm{H}\mathrm{z}$$, with the recordings
of the last $10$ speakers retained as a validation set. We set the recording
length to $32\,768$ samples, and use data augmentations such as random crop,
phase mangle, or dequantization proposed in RAVE [15]. We train the models for
1.25M steps using the AdamW optimizer [20] with a learning rate of $510-5$ and
a batch size of $16$. In Fig. 2, we show waveforms and spectrograms of the
reconstructions obtained with our model on samples selected from the
validation set.

Since there is no consensus on a single approach for quantitative evaluation
of audio quality, we assess the reconstruction results with multiple metrics
such as MSE, Log-Spectral Distance (LSD) [21], SI-SNR [22], PESQ [23], STOI
[24] and CDPAM [25]. We also compare the reconstruction quality of our model
with the RAVE baseline. Moreover, we test the quality of resampling performed
by our model. Finally, we investigate the impact of the target network size
and the employed loss function on the reconstruction quality. The results of
our experiments are shown in Tab. 1. Unless explicitly mentioned, model
hyperparameters in experiments are as described in Sec. 3.

Comparision with RAVE. We compare the reconstruction error of HyperSound with
RAVE [15] trained for 3M steps on the same dataset, with an identical sampling
rate of $22\,050\text{\,}\mathrm{H}\mathrm{z}$. We find that our model
produces reconstructions closer to the originals in the spectral domain and
obtains better perceptual scores of PESQ and STOI. However, RAVE
reconstructions achieve better MSE, SI-SNR, and perceptual CDPAM scores, with
slightly less perceptible noise and robotic artifacts.

Reconstruction quality while resampling. We also test the reconstruction error
of the model trained with a sampling rate of
$22\,050\text{\,}\mathrm{H}\mathrm{z}$ when used for downsampling and
upsampling $22\,050\text{\,}\mathrm{H}\mathrm{z}$ recordings to other commonly
used frequencies: $8000\text{\,}\mathrm{H}\mathrm{z}$,
$16\,000\text{\,}\mathrm{H}\mathrm{z}$ and
$44\,100\text{\,}\mathrm{H}\mathrm{z}$. We obtain the ground truth by
downsampling the original VCTK $48\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$
recordings to the desired sampling frequency using soxr high-quality setting.
As shown by the LSD levels, the hypernetwork approach works best when
reconstructing signals close to the sampling rate used for training.
Nevertheless, it does not collapse when queried with time coordinates from
outside the training domain. However, we hypothesize that for proper super-
resolution capabilities, our learning regime would require an introduction of
an additional loss term directed specifically at the quality of representation
in higher frequency bands.

Impact of target network architecture. We compare our model which uses the
baseline target network (4 layers of $256$ neurons each, 206K parameters) with
target networks that use $4$ layers of $64$ neurons (14K parameters) and $6$
layers of $384$ neurons (752K parameters). We find that the base version
presents an optimal trade-off between the reconstruction quality and
computational requirements, as the size of the last hypernetwork layer scales
linearly with the number of parameters in the target network. The smaller
variant of the target network achieves comparable results to the base using
fewer parameters than the number of samples in the original signal.

Impact of the loss function. We compare HyperSound trained with the loss
function as described in Sec. 3.3 with models trained using the STFT loss
where we do not apply the mel-scale, but use a wider array of FFT sizes of
$[128,256,512,1024,2048]$. We also try training without Smooth L1 loss
(${\lambda}_{SL1}=0,{\lambda}_{STFT}=1$). Moreover, we find that STFT is vital
for stable training and obtaining perceptually plausible reconstructions, as
our training runs with ${\lambda}_{SL1}=1.$ and ${\lambda}_{STFT}=[0.1,0.01]$
collapsed. As evidenced by the results, the L1 part of the loss function leads
to slightly better results and faster training, as it enables the model to
correctly learn the DC offset.

Table 1: HyperSound evaluation based on reconstructions of the VCTK validation set |  | MSE | LSD | SI-SNR | PESQ | STOI | CDPAM
---|---|---|---|---|---|---|---
| Ideal metric behavior | $\to 0$ | $\to 0$ | $\to 100$ | $\to 4.5$ | $\to 1$ | $\to 0$
Comparison with RAVE | Model |
HyperSound∗ | 0.049 | 0.99 | -23.50 | 1.75 | 0.87 | 0.35
RAVE | 0.040 | 1.19 | -22.42 | 1.22 | 0.74 | 0.19
Reconstruction quality while resampling | Target SR [Hz] |
$8000$ | 0.048 | 1.24 | -23.76 | 1.67 | 0.83 | 0.37
$16\,000$ | 0.047 | 1.19 | -23.80 | 1.66 | 0.86 | 0.35
22 050∗ | 0.049 | 0.99 | -23.50 | 1.75 | 0.87 | 0.35
$44\,100$ | 0.047 | 1.38 | -24.16 | 1.72 | 0.87 | 0.38
Impact of target network size | Network size |
Small [4 x 64] | 0.048 | 1.07 | -24.56 | 1.56 | 0.84 | 0.35
Base [4 x 256]∗ | 0.049 | 0.99 | -23.50 | 1.75 | 0.87 | 0.35
Large [6 x 384] | 0.047 | 0.96 | -23.70 | 1.86 | 0.89 | 0.36
Impact of the loss function | Loss function |
L1$\,+\,$STFT | 0.052 | 0.94 | -25.05 | 1.35 | 0.81 | 0.25
Only STFT | 0.055 | 0.94 | -25.63 | 1.34 | 0.81 | 0.23
L1$\,+\,$MelSTFT∗ | 0.049 | 0.99 | -23.50 | 1.75 | 0.87 | 0.35
Only MelSTFT | 0.050 | 1.02 | -23.31 | 1.70 | 0.87 | 0.29
∗default HyperSound hyperparameters

##  5 Conclusion

We demonstrate the possibility of applying hypernetworks to the generation of
implicit neural representations for audio signals. Reconstructions generated
with our model are quantitatively comparable to the state-of-the-art model
RAVE. However, we find that the perceptual quality of our reconstructions is
still slightly lacking. We hope that our work can be further improved by
optimizing the hypernetwork architecture and designing target networks better
suited to the audio domain. Initial results for signal compression are also
promising but require further investigation.

## Acknowledgments and Disclosure of Funding

This work was supported by Foundation for Polish Science (grant no
POIR.04.04.00-00-14DE/18-00) carried out within the Team-Net program co-
financed by the European Union under the European Regional Development Fund,
as well as the National Centre of Science (Poland) Grant No.
2020/39/B/ST6/01511.

## References

  * Sitzmann et al. [2020] Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein.  Implicit Neural Representations with Periodic Activation Functions.  _Advances in Neural Information Processing Systems_ , 33:7462–7473, 2020.
  * Mehta et al. [2021] Ishit Mehta, Michaël Gharbi, Connelly Barnes, Eli Shechtman, Ravi Ramamoorthi, and Manmohan Chandraker.  Modulated Periodic Activations for Generalizable Local Functional Representations.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 14214–14223, 2021.
  * Mildenhall et al. [2021] Ben Mildenhall, Pratul P Srinivasan, Matthew Tancik, Jonathan T Barron, Ravi Ramamoorthi, and Ren Ng.  NeRF: Representing Scenes as Neural Radiance Fields for View Synthesis.  _Communications of the ACM_ , 65(1):99–106, 2021\.
  * Ha et al. [2016] David Ha, Andrew Dai, and Quoc V Le.  HyperNetworks.  _arXiv preprint arXiv:1609.09106_ , 2016.
  * Klocek et al. [2019] Sylwester Klocek, Łukasz Maziarka, Maciej Wołczyk, Jacek Tabor, Jakub Nowak, and Marek Śmieja.  Hypernetwork functional image representation.  In _International Conference on Artificial Neural Networks_ , pages 496–510. Springer, 2019.
  * Spurek et al. [2020] Przemysław Spurek, Sebastian Winczowski, Jacek Tabor, Maciej Zamorski, Maciej Zięba, and Tomasz Trzciński.  Hypernetwork approach to generating point clouds.  _arXiv preprint arXiv:2003.00802_ , 2020.
  * Zhao et al. [2020] Dominic Zhao, Johannes von Oswald, Seijin Kobayashi, João Sacramento, and Benjamin F Grewe.  Meta-Learning via Hypernetworks.  2020\.
  * Von Oswald et al. [2019] Johannes Von Oswald, Christian Henning, João Sacramento, and Benjamin F Grewe.  Continual learning with hypernetworks.  _arXiv preprint arXiv:1906.00695_ , 2019.
  * Oord et al. [2016] Aaron van den Oord, Sander Dieleman, Heiga Zen, Karen Simonyan, Oriol Vinyals, Alex Graves, Nal Kalchbrenner, Andrew Senior, and Koray Kavukcuoglu.  WaveNet: A Generative Model for Raw Audio.  _arXiv preprint arXiv:1609.03499_ , 2016.
  * Mehri et al. [2016] Soroush Mehri, Kundan Kumar, Ishaan Gulrajani, Rithesh Kumar, Shubham Jain, Jose Sotelo, Aaron Courville, and Yoshua Bengio.  SampleRNN: An Unconditional End-to-End Neural Audio Generation Model.  _arXiv preprint arXiv:1612.07837_ , 2016.
  * Oord et al. [2018] Aaron Oord, Yazhe Li, Igor Babuschkin, Karen Simonyan, Oriol Vinyals, Koray Kavukcuoglu, George Driessche, Edward Lockhart, Luis Cobo, Florian Stimberg, et al.  Parallel WaveNet: Fast High-Fidelity Speech Synthesis.  In _International conference on machine learning_ , pages 3918–3926. PMLR, 2018.
  * Engel et al. [2017] Jesse Engel, Cinjon Resnick, Adam Roberts, Sander Dieleman, Mohammad Norouzi, Douglas Eck, and Karen Simonyan.  Neural Audio Synthesis of Musical Notes with WaveNet Autoencoders.  In _International Conference on Machine Learning_ , pages 1068–1077. PMLR, 2017.
  * Kumar et al. [2019] Kundan Kumar, Rithesh Kumar, Thibault de Boissiere, Lucas Gestin, Wei Zhen Teoh, Jose Sotelo, Alexandre de Brébisson, Yoshua Bengio, and Aaron C Courville.  MelGAN: Generative Adversarial Networks for Conditional Waveform Synthesis.  _Advances in neural information processing systems_ , 32, 2019.
  * Défossez et al. [2018] Alexandre Défossez, Neil Zeghidour, Nicolas Usunier, Léon Bottou, and Francis Bach.  SING: Symbol-to-Instrument Neural Generator.  _Advances in neural information processing systems_ , 31, 2018.
  * Caillon and Esling [2021] Antoine Caillon and Philippe Esling.  RAVE: A variational autoencoder for fast and high-quality neural audio synthesis.  _arXiv preprint arXiv:2111.05011_ , 2021.
  * Zeghidour et al. [2021] Neil Zeghidour, Alejandro Luebs, Ahmed Omran, Jan Skoglund, and Marco Tagliasacchi.  SoundStream: An End-to-End Neural Audio Codec.  _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , 30:495–507, 2021.
  * Clevert et al. [2015] Djork-Arné Clevert, Thomas Unterthiner, and Sepp Hochreiter.  Fast and Accurate Deep Network Learning by Exponential Linear Units (ELUs).  _arXiv preprint arXiv:1511.07289_ , 2015.
  * Girshick [2015] Ross Girshick.  Fast R-CNN.  In _Proceedings of the IEEE international conference on computer vision_ , pages 1440–1448, 2015.
  * Yamamoto et al. [2020] Ryuichi Yamamoto, Eunwoo Song, and Jae-Min Kim.  Parallel WaveGAN: A fast waveform generation model based on generative adversarial networks with multi-resolution spectrogram.  In _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 6199–6203. IEEE, 2020.
  * Loshchilov and Hutter [2017] Ilya Loshchilov and Frank Hutter.  Decoupled Weight Decay Regularization.  _arXiv preprint arXiv:1711.05101_ , 2017.
  * Liu et al. [2022] Haohe Liu, Woosung Choi, Xubo Liu, Qiuqiang Kong, Qiao Tian, and DeLiang Wang.  Neural Vocoder is All You Need for Speech Super-resolution.  _arXiv preprint arXiv:2203.14941_ , 2022.
  * Luo and Mesgarani [2018] Yi Luo and Nima Mesgarani.  TasNet: time-domain audio separation network for real-time, single-channel speech separation.  In _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 696–700. IEEE, 2018.
  * Rix et al. [2001] Antony W Rix, John G Beerends, Michael P Hollier, and Andries P Hekstra.  Perceptual evaluation of speech quality (PESQ)-a new method for speech quality assessment of telephone networks and codecs.  In _2001 IEEE international conference on acoustics, speech, and signal processing. Proceedings (Cat. No. 01CH37221)_ , volume 2, pages 749–752. IEEE, 2001.
  * Taal et al. [2010] Cees H Taal, Richard C Hendriks, Richard Heusdens, and Jesper Jensen.  A short-time objective intelligibility measure for time-frequency weighted noisy speech.  In _2010 IEEE international conference on acoustics, speech and signal processing_ , pages 4214–4217. IEEE, 2010.
  * Manocha et al. [2021] Pranay Manocha, Zeyu Jin, Richard Zhang, and Adam Finkelstein.  CDPAM: Contrastive learning for perceptual audio similarity.  In _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 196–200. IEEE, 2021.