# Nonlinear Evolutionary PDE-Based Refinement of Optical Flow

Hirak Doshi, N. Uday Kiran
Department of Mathematics and Computer Science
Sri Sathya Sai Institute of Higher Learning, India
Email<EMAIL_ADDRESS>
###### Abstract

The goal of this paper is to propose two nonlinear variational models for
obtaining a refined motion estimation from an image sequence. Both the
proposed models can be considered as a part of a generalized framework for an
accurate estimation of physics-based flow fields such as rotational and fluid
flow. The first model is novel in the sense that it is divided into two
phases: the first phase obtains a crude estimate of the optical flow and then
the second phase refines this estimate using additional constraints. The
correctness of this model is proved using an Evolutionary PDE approach. The
second model achieves the same refinement as the first model, but in a
standard manner, using a single functional. A special feature of our models is
that they permit us to provide efficient numerical implementations through the
first-order primal-dual Chambolle-Pock scheme. Both the models are compared in
the context of accurate estimation of angle by performing an anisotropic
regularization of the divergence and curl of the flow respectively. We observe
that, although both the models obtain the same level of accuracy, the two-
phase model is more efficient. In fact, we empirically demonstrate that the
single-phase and the two-phase models have convergence rates of order
$O(1/N^{2})$ and $O(1/N)$ respectively.

######  Keywords.

Optical Flow, Evolutionary PDE, Variational Methods, Primal-Dual, Convergence

######  Mathematics Subject Classification 2020.

35A15, 35J47, 35Q68.

##  1 Introduction

Optical flow plays a key role in many advanced Computer Vision applications.
It is a rich source of information on perceptible motion in our visual world.
It’s reliable estimation is thus important and at the same time challenging.
Assuming the principle of local conservation of intensity and small temporal
variations, optical flow involves the recovery of a function
$\textbf{u}=(u,v)$ such that

$f(\textbf{x},\tau)=f(\textbf{x}+\textbf{u},\tau+\Delta\tau).$

where $f:\Omega\times[0,T]\to\mathbb{R}$ if the image sequence,
$\Omega\subset\mathbb{R}^{2}$ is open and bounded. This establishes a
correspondence between pixel motions. Using first-order approximations the
above relation can be written as

$f_{\tau}+\nabla f\cdot\textbf{u}=0,$ (1)

which is widely known as the Optical Flow Constraint (OFC). To recover the
velocity components using a variational minimization approach one writes (1)
as

$\min_{\textbf{u}}J_{1}(\textbf{u})=\int_{\Omega}(f_{t}+\nabla
f\cdot\textbf{u})^{2}.$

This is the simplest least-square minimization. This problem is ill-posed as
it leads to the aperture problem. Additional regularization terms are
necessary to ensure well-posedness. The most common regularization term is the
quadratic smoothness penalizing the gradient of the components of the flow
originally introduced by Horn and Schunck [22] in their seminal work. Cohen
[11] and Kumar et.al. [23] used the $L^{1}$ regularization which is more
robust to outliers and preserves important edge information. A new
discontinuity preserving optical flow model with $L^{1}$ norm on the OFC was
proposed and studied by Aubert et. al. [3] in the space of functions of
bounded variations $BV(\Omega)\times BV(\Omega)$. The well-posedness of the
Horn and Schunck model, as well as the Nagel model was studied by Schnörr [28]
in the space $H^{1}(\Omega)\times H^{1}(\Omega)$. Taking a step further, the
authors in [8] proposed a $L^{p}-TV/L^{p}$ ($p=1$ or $2$) model combining both
$L^{1}$ and $L^{2}$ terms. The behaviour of their regularization term is
similar to the Huber function:

$H(x;\epsilon)=\begin{cases}\frac{x^{2}}{2\epsilon},\quad
0\leq|x|\leq\epsilon\\\\[5.69046pt]
|x|-\frac{\epsilon}{2},\quad|x|>\epsilon.\end{cases}$

A detailed review and rigorous analysis of several variational optical flow
models within the framework of calculus of variations can be found in [21].

Though most of the estimation involving rigid or quasi-rigid motion can be
handled by minimizing OFC with a suitable regularization, it is insufficient
to provide an accurate estimation for fluid-based images. Traditional computer
vision techniques may not be suitable to capture these deformations of
brightness patterns because of the high spatio-temporal turbulence in these
sequences. These reasons have motivated researchers to look for an alternate
constraint that can not only preserve pixel-correspondence but also capture
certain intrinsic features of the flow. This paradigm shift hints at
constraints that are physics-dependent. A lot of work has been done involving
physics-based constraints for fluid motion estimation [12, 13, 25, 26, 27]. In
[15], we have proposed a constraint-based refinement of optical flow. Using an
image-driven evolutionary PDE model resulting from a quadratic regularization
we have shown the well-posedness of such a refinement principle. An important
characteristic of the model is the possibility of a diagonalization by the
Cauchy-Riemann operator leading to a decoupled system involving diffusion of
the curl and a multiplicative perturbation of the laplacian of the divergence
of the flow. For a specific case, it was shown that the model is close to the
physics-based model [25] using a modified augmented Lagrangian method.

A topic of major interest to the Computer Vision community is the performance
evaluation methods for optical flow estimation. Several important
contributions [4, 6, 7, 35] discuss robust methods for obtaining high accuracy
flow fields. Motivated by the harmonic-constraint based regularization [36],
we proposed a variational model for capturing rotational features and
improving angular accuracy of the flow [16].

The current work proposes a unified framework that extends our previous
methods in a nonlinear setting. The first model is a two-phase refinement
process. A crude pixel correspondence is obtained in the first phase which
constitutes a good starting solution. In the next phase, this estimate is
refined using additional constraints. The second model estimates the flow
directly in a single phase. The additional constraint is chosen motivated by
the harmonic constraint-based regularization (see [36]). This approach
replaces the oriented smoothness constraint with a weighted decomposition of
divergence and curl of flow. We aim to capture the rotational features better
by preserving edge information and improving the accuracy of the flow. Thus we
consider only the curl component in our framework with an anisotropic weight
term.

The total variation regularization leads to $\Delta_{1}$, the 1-Laplacian
operator in the Euler-Lagrange equations. Obtaining a stable convergent scheme
is a difficult task because of the singularity of the operator at the origin.
As a result most of the implementation methods often yield slower algorithms.
In this direction important contributions were made by Chambolle [9] and Zach
[34]. Chambolle and Pock [10] proposed a first-order primal-dual algorithm for
solving non-smooth convex optimization problems. This further opened up newer
directions as a large class of problems in Image Processing and Computer
Vision could be solved within this framework. In our work, we use the
Chambolle-Pock algorithm for both of our models. The numerical implementation
of the algorithm has two main steps, namely updating the primal variables by
solving a system of equations at each iteration and updating the dual
variables by computing the point-wise projection maps onto the unit ball. Both
of these steps are computationally expensive. As a result the second model
yields a slower algorithm with a convergence rate of order $O(1/N^{2})$. By
$O(1/N^{2})$ we mean if $\epsilon$ is the error threshold then the number of
iterations required to reach this threshold is $1/\epsilon^{2}$. The first
model splits the above-mentioned steps in two-phases. This leads to a faster
algorithm with a convergence rate of $O(1/N)$.

The organization of the paper is as follows. In Section (2) we give the
general formulation and describe our model in detail. Next in Section (3), we
study the mathematical well-posedness of our formulation using an evolutionary
PDE approach. Subsequently, we employ the first-order primal-dual Chambolle-
Pock algorithm to our models and derive the necessary optimality conditions in
Section (4). We then discuss the implementation details, discretization of our
models and empirically demonstrate the nature of convergence in Section (5).

##  2 Our Model Description

Our general formulation is given as:

$J(\textbf{u})=\int_{\Omega}\rho(|f_{t}+\nabla
f\cdot\textbf{u}|)+\alpha\sum_{i=1}^{2}\int_{\Omega}\gamma(|\nabla
u_{i}|)+\beta\int_{\Omega}\phi(x,f,\nabla f)\psi(\textbf{u},\nabla\textbf{u})$
(2)

where $\rho:\mathbb{R}\to\mathbb{R}$ is a function of the optical flow
constraint, $\gamma:\mathbb{R}\to\mathbb{R}$ governs the regularization of the
flow. The functions $\phi$ and $\psi$ are chosen specific to applications. A
summary of some of the variational models which belong to this framework is
listed in the following table:

| $\rho(x)$ | $\gamma(x)$ | $\phi(x,f,\nabla f)$ | $\psi(\textbf{u},\nabla\textbf{u})$
---|---|---|---|---
Horn and Schunck[22] | $x^{2}$ | $x^{2}$ | 0 | 0
Cohen[11] | $x^{2}$ | $x$ | 0 | 0
Aubert[3] | $x$ | $\sqrt{1+x^{2}}$ | $c(x)$ | $\textbf{u}^{2}$
$L^{1}-TV$[34] | $x$ | $x$ | 0 | 0
Our Model (M1) | 0 | $x$ | $f^{2}$ | $(\nabla\cdot\textbf{u})^{2}$
Our Model (M2) | $x^{2}$ | $x$ | $\frac{\lambda^{2}}{\|\nabla f\|^{2}+\lambda^{2}}$ | $(\nabla_{H}\cdot\textbf{u})^{2}$
Table 1: Some choices for the functions

Based on the choice functions mentioned in the above table, the constraint-
based refinement formulation becomes

$(M1)\quad J(\textbf{u})=\alpha\sum_{i=1}^{2}\int_{\Omega}|\nabla
u_{i}|+\beta\int_{\Omega}f^{2}(\nabla\cdot\textbf{u})^{2}.$ (3)

where $|\cdot|$ is the Euclidean norm. Starting with
$\textbf{u}=\textbf{u}_{0}$, where $\textbf{u}_{0}$ is the Horn and Schunck
optical flow the above formulation obtains a refinement of $\textbf{u}_{0}$
driven by the additional constraint $f\nabla\cdot\textbf{u}$. This term is the
non-conservation term in the physics-based constraint due to non null-out of
plane components [20]. This constraint preserves the spatial characteristics
and vorticities of the flow. Thus this model can extract flow informations
from fluid-based digital imagery much better.

From the choice functions, our second formulation is given as:

$(M2)\quad J(\textbf{u})=\int_{\Omega}(f_{t}+\nabla
f\cdot\textbf{u})^{2}+\alpha\sum_{i=1}^{2}\int_{\Omega}|\nabla
u_{i}|+\beta\int_{\Omega}\frac{\lambda^{2}}{|\nabla
f|^{2}+\lambda^{2}}(\nabla_{H}\cdot\textbf{u})^{2}.$ (4)

The additional constraint in this formulation is motivated by the harmonic-
constraint based regularization discussed in [36]. By associating an
anisotropic weight term with the curl in our formulation there are two main
advantages. First, we are able to get a precise estimation of the
infinitesimal rotation within the regions. Secondly, we achieve a better
alignment of small vectors in the flow. This leads to an overall improvement
in the endpoint error.

In either case, the choice of $\phi$ decides the influence of the image term
in the regularization process. If $\phi=1$, then the additional constraint
term in the functional is flow-driven, i.e. independent of the influence of
the image data. The weight parameters $\alpha,\beta$ play an important role in
the regularization process. For rigid-body like motion which requires
important edge-information to be preserved a higher value of $\alpha$ is
preferred. For fluid-based images where there is less edge-prominence, we
choose a higher $\beta$ value.

##  3 Well-Posedness

In this section we discuss the mathematical well-posedness of the proposed
formulation. Let us denote $\textbf{u}=(u_{1},u_{2})$. The space
$W^{1,p}(\Omega),p>1$, is the reflexive Banach space

$W^{1,p}(\Omega)=\\{u\in L^{p}(\Omega):D^{\alpha}u\in
L^{p}(\Omega),|\alpha|\leq 1\\}$

with the usual norm

$\|u\|_{W^{1,p}}=\Bigg{(}\sum_{|\alpha|\leq
1}\|D^{\alpha}u\|^{p}_{L^{p}}\Bigg{)}^{1/p},1\leq p<\infty.$

To study the mathematical well-posedness of the proposed formulation we
consider the following approximation.

$J_{p,\text{R}}(\textbf{u})=\beta\int_{\Omega}\phi(x,f,\nabla
f)(\nabla\cdot\textbf{u})^{2}+\frac{\alpha}{p}\int_{\Omega}\\{|\nabla
u_{1}|^{p}+|\nabla u_{2}|^{p}\\},\quad 1<p<2$ (5)

This functional being strictly convex in $W^{1,p}(\Omega)$ admits a unique
minimizer. For this discussion we consider $\phi(x,f,\nabla f)=f^{2}$. The
first important step is to show that $J_{p,\textbf{R}}$ converges to
$J_{1,\textbf{R}}$ as $p\to 1$. For this, we refer to the discussion in
Section 3.4 in [24].

######  Lemma 1.

$\lim_{p\to 1}\frac{1}{p}\int_{\Omega}|\nabla u|^{p}=\int_{\Omega}|\nabla u|.$
(6)

######  Remark 1.

As $J_{p,\text{R}}(\textbf{u})\to J_{1,\text{R}}(\textbf{u})$, $p\to 1$, the
corresponding Euler-Lagrange equations $A_{p}=\Delta_{p}$ associated with the
regularization term also converges to $A_{1}=\Delta_{1}$.

######  Remark 2.

The case $p=2$ leads to a linear diffusion-driven refinement process. We have
previously studied and discussed this case in [15].

The associated parabolic system corresponding to the Euler-Lagrange equations
of (5) are given as

$\begin{cases}\dfrac{\partial u_{1}}{\partial
t}=\Delta_{p}u_{1}+a_{0}\dfrac{\partial}{\partial
x}[f^{2}((u_{1})_{x}+(u_{2})_{y})]\text{ in
}\Omega\times(0,\infty),\\\\[14.22636pt] \dfrac{\partial u_{2}}{\partial
t}=\Delta_{p}u_{2}+a_{0}\dfrac{\partial}{\partial
y}[f^{2}((u_{1})_{x}+(u_{2})_{y})]\text{ in
}\Omega\times(0,\infty),\\\\[14.22636pt] u_{1}(x,y,0)=u_{1}^{0}\text{ in
}\Omega,\\\\[8.5359pt] u_{2}(x,y,0)=u_{2}^{0}\text{ in }\Omega,\\\\[8.5359pt]
u_{1}=0\text{ on }\partial\Omega\times(0,\infty),\\\\[8.5359pt] u_{2}=0\text{
on }\partial\Omega\times(0,\infty),\end{cases}$ (7)

where $(u_{1}^{0},u_{2}^{0})$ is the starting feasible solution obtained by
the Horn and Schunck optical flow, $a_{0}=2\beta/\alpha$. Rewriting the system
in an abstract form leads us to

$\begin{cases}\dfrac{d\textbf{u}}{dt}+\mathcal{A}_{p}\textbf{u}=0,\quad
t>0,\\\\[8.5359pt] \textbf{u}(0)=\textbf{u}_{0}\in
H^{1}(\Omega)^{2}.\end{cases}$ (8)

Here the operator $\mathcal{A}_{p}=A_{p}+F$ where

$A_{p}\textbf{u}=-\begin{bmatrix}\Delta_{p}u_{1}\\\\[8.5359pt]
\Delta_{p}u_{2}\end{bmatrix},\qquad
F\textbf{u}=-a_{0}\begin{bmatrix}\dfrac{\partial}{\partial
x}[f^{2}((u_{1})_{x}+(u_{2})_{y})]\\\\[8.5359pt] \dfrac{\partial}{\partial
y}[f^{2}((u_{1})_{x}+(u_{2})_{y})]\end{bmatrix},$

We will show that both the operators $A_{p}$ and $F$ are maximal monotone in
$W^{1,p}(\Omega)\cap L^{2}(\Omega)$ and $L^{2}(\Omega)$ respectively.

######  Lemma 2.

The operators $A_{p}$ and $F$ is maximal monotone in $W^{1,p}(\Omega)\cap
L^{2}(\Omega)$ and $L^{2}(\Omega)$ respectively.

###### Proof.

The maximal monotonicity of $A_{p}$ follows directly from the discussions in
[33]. To show monotonicity we show that $\langle
F\textbf{u},\textbf{u}\rangle\geq 0$. Indeed,

$\displaystyle\langle F\textbf{u},\textbf{u}\rangle$
$\displaystyle=-a_{0}\int_{\Omega}\Big{\\{}\frac{\partial}{\partial
x}[f^{2}((u_{1})_{x}+(u_{2})_{y})]u+\frac{\partial}{\partial
y}[f^{2}((u_{1})_{x}+(u_{2})_{y})]v\Big{\\}},$
$\displaystyle=a_{0}\int_{\Omega}\\{f^{2}((u_{1})_{x}+(u_{2})_{y})u_{x}+f^{2}((u_{1})_{x}+(u_{2})_{y})v_{y}\\},$
$\displaystyle=a_{0}\int_{\Omega}f^{2}((u_{1})_{x}+(u_{2})_{y})^{2}\geq 0,$

proving the monotonicity of $F$. To show the maximality we have to show that

$\text{Ran}(I+F)=L^{2}(\Omega)^{2},$ (9)

i.e. there exists u for all $\textbf{f}\in L^{2}(\Omega)^{2}$ such that
$\textbf{u}+F\textbf{u}=\textbf{f}$ holds. Let $\textbf{f}=(f,g)\in
L^{2}(\Omega)^{2}$ and consider the system

$u+\frac{\partial}{\partial x}[f^{2}((u_{1})_{x}+(u_{2})_{y})]=f,$
$v+\frac{\partial}{\partial y}[f^{2}((u_{1})_{x}+(u_{2})_{y})]=g,$

where $f,g\in L^{2}(\Omega)$. Applying the Cauchy-Riemann operator

$R=\begin{bmatrix}\partial_{y}&-\partial_{x}\\\\[8.5359pt]
\partial_{x}&\partial_{y}\end{bmatrix}.$

both sides we obtain the decoupled system

$\displaystyle(u_{1})_{y}-(u_{2})_{x}=f_{y}-g_{x},$ (10)
$\displaystyle(\Delta\circ k)((u_{1})_{x}+(u_{2})_{y})=f_{x}+g_{y}$ (11)

where $k$ is the image-dependent multiplicative function $k:f\mapsto
1+a_{0}f^{2}$. The first equation (10) governs the curl of the flow
$\textbf{u}=(u_{1},u_{2})$. The second equation (11) indicates a non-
homogeneous weighted diffusion process on the divergence with a weight $k$.
Let us define $h_{1}=f_{y}-g_{x}$ and $h_{2}=f_{x}+g_{y}$. Solving the second
equation gives us an expression for the divergence of the flow. Let us call
this as $h_{3}$. We thus obtain the following system

$\displaystyle(u_{1})_{y}-(u_{2})_{x}$ $\displaystyle=h_{1},$
$\displaystyle(u_{1})_{x}+(u_{2})_{y}$ $\displaystyle=h_{3}/k,$

These are the inhomogeneous Cauchy-Riemann equations. In a compact form we
rewrite them as

$R\textbf{u}=\tilde{\textbf{f}}$

where $\tilde{\textbf{f}}=(h_{1},h_{3}/k)$. The operator $R^{-1}$ is a
continuous operator of order -1 in the space $W^{1,p}(\Omega)$. Hence there
exists a unique u such that $\textbf{u}+F\textbf{u}=\textbf{f}$ holds. This
concludes the proof. ∎

Now that we have shown the maximal monotonicity let us define a function
$\Phi_{p}:L^{2}(\Omega)^{2}\to(-\infty,+\infty]$ by

$\Phi_{p}(\textbf{u})=\begin{cases}J_{p,\textbf{R}}(\textbf{u}),\qquad\textbf{u}\in[W^{1,p}(\Omega)\cap
L^{2}(\Omega)]^{2}\\\\[8.5359pt]
+\infty,\qquad\quad\textbf{u}\in[L^{2}(\Omega)\setminus
W^{1,p}(\Omega)]^{2}\end{cases}.$

Then clearly $\Phi_{p}$ is convex and lower semi-continuous. Also $\Phi_{p}$
is proper since $D(\Phi_{p})=D(A_{p})\cap D(F)\neq\emptyset$. Thus the
associated subdifferential $\partial\Phi(\textbf{u})\equiv\mathcal{A}_{p}$ is
maximal monotone. Thus there is a unique solution u of the inclusion

$0\in\textbf{u}^{\prime}(t)+\partial\Phi_{p}(\textbf{u})$

satisfying the initial conditions.

##  4 The Primal-Dual Framework

The primal-dual method is a numerical tool for solving optimization problems.
The main idea is to replace a primal problem with an equivalent saddle point
problem by introducing dual variables and employ efficient algorithms to
obtain the desired convergence. In the recent past several saddle point
frameworks have been proposed for variational problems in image processing and
computer vision [10, 31, 35]. As our formulation involves non-smooth convex
functionals, the most suitable framework is the one proposed by Chambolle and
Pock [10]. Let $\Omega\subset\mathbb{R}^{2}$ be an open, bounded set,
$\mathcal{X},\mathcal{Y}$ be two finite-dimensional vector spaces with the
scalar product $(\cdot,\cdot)$ and the norm $\|\cdot\|$. Denote the primal
variable $\textbf{u}=(u_{1},u_{2})$ and the dual variable
$\textbf{d}=(d_{1},d_{2},d_{3})$. We first consider the variational problem in
the following form

$\operatorname*{arg\,min}_{\textbf{u}}G(\textbf{u})+F(K\textbf{u}).$ (12)

where $F,G:\mathcal{X}\to[0,\infty]$ are convex, proper and lower-
semicontinuous functionals, $K:\mathcal{X}\to\mathcal{Y}$ is a continuous,
linear operator. The equivalent primal-dual formulation is given as

$\operatorname*{arg\,min}_{\textbf{u}}\operatorname*{arg\,max}_{\textbf{d}}\>G(\textbf{u})+(K\textbf{u},\textbf{d})-F^{*}(\textbf{d}).$
(13)

where $F^{*}$ is the convex conjugate of $F$. Table (2) gives a summary of
each term of our model using the above notations. Given a $\tau,\sigma>0$, an
initial $(\textbf{u}^{0},\textbf{d}^{0})\in\mathcal{X}\times\mathcal{Y}$, the
Chambolle-Pock algorithm solves the saddle point problem (13) by the following
algorithm:

$\displaystyle\textbf{d}^{k+1}$
$\displaystyle=\operatorname{\textbf{prox}}_{\sigma
F^{*}}(\textbf{d}^{k}+\sigma K\bar{\textbf{u}}^{k}),$
$\displaystyle\textbf{u}^{k+1}$
$\displaystyle=\operatorname{\textbf{prox}}_{\tau G}(\textbf{u}^{k}-\tau
K^{*}\textbf{d}^{k+1}),$ $\displaystyle\bar{\textbf{u}}_{k+1}$
$\displaystyle=2\textbf{u}_{k+1}-\textbf{u}_{k}\qquad\text{(over-
relaxation)},$

where

$\operatorname{\textbf{prox}}_{\tau G}(\textbf{d})=(I+\tau\partial
G)^{-1}(\textbf{d})=\operatorname*{arg\,min}_{\textbf{u}}\Bigg{\\{}\frac{1}{2}\|\textbf{u}-\textbf{d}\|^{2}+\tau
G(\textbf{u})\Bigg{\\}}$

is the proximal or the resolvent operator. This can be thought of as a trade-
off between minimizing $G$ and being close to d. We now employ the above
algorithm for our problem and derive the necessary optimality conditions.

###  4.1 Optimality Condition for Our Model $(M1)$

In this case we have

$G(\textbf{u})=0,\quad
F(K\textbf{u})=\frac{1}{2}\int_{\Omega}f^{2}(\nabla\cdot\textbf{u})^{2}+\sum_{i=1}^{2}\int_{\Omega}|\nabla
u_{i}|.$

The Operator $K$ is given as

$K\textbf{u}=\begin{bmatrix}\nabla&0\\\ 0&\nabla\\\
f^{2}\partial_{x}&f^{2}\partial_{y}\end{bmatrix}\begin{bmatrix}u_{1}\\\\[11.38092pt]
u_{2}\end{bmatrix}.$

Therefore,

$K^{*}\textbf{d}=-\begin{bmatrix}\nabla\cdot&0&\partial_{x}(f^{2}\cdot)\\\
0&\nabla\cdot&\partial_{y}(f^{2}\cdot)\end{bmatrix}\begin{bmatrix}d_{1}\\\
d_{2}\\\ d_{3}\end{bmatrix}.$

Model | $\phi(x,f,\nabla f)$ | $G(\textbf{u})$ | $F(K\textbf{u})$ | $K$
---|---|---|---|---
Model (M1) | $f^{2}$ | 0 | $\begin{aligned} \frac{1}{2}\int_{\Omega}\phi(x,f,\nabla f)(\nabla\cdot\textbf{u})^{2}\\\ &+\sum_{i=1}^{2}\int_{\Omega}|\nabla u_{i}|\end{aligned}$ | $\begin{pmatrix}\nabla&0\\\ 0&\nabla\\\ \phi\partial_{x}&\phi\partial_{y}\end{pmatrix}$
Model (M2) | $\dfrac{\lambda^{2}}{\lambda^{2}+\|\nabla f\|^{2}}$ | $\displaystyle\frac{1}{2}\int_{\Omega}(f_{t}+\nabla f\cdot\textbf{u})^{2}$ | $\begin{aligned} \frac{1}{2}\int_{\Omega}\phi(x,f,\nabla f)(\nabla_{H}\cdot\textbf{u})^{2}\\\ &+\sum_{i=1}^{2}\int_{\Omega}|\nabla u_{i}|\end{aligned}$ | $\begin{pmatrix}\nabla&0\\\ 0&\nabla\\\ \phi\partial_{y}&-\phi\partial_{x}\end{pmatrix}$
Table 2: Summary of the terms in the Primal-Dual Formulation

Using standard dual identities the convex conjugate $F^{*}(\textbf{d})$ is
computed as

$F^{*}(\textbf{d})=\frac{1}{2}\|d_{3}\|^{2}_{2}+\alpha\sum_{i=1}^{2}\delta_{B(L^{\infty})}(d_{i}/\alpha),$

where $B(L^{\infty})$ denotes the unit ball in $L^{\infty}(\Omega)$ and

$\delta_{B(L^{\infty})}(x^{*})=\begin{cases}0,\text{ if }x^{*}\in
B(L^{\infty})\\\ +\infty,\text{ otherwise}\end{cases}.$

Thus the primal-dual formulation is given as

$\operatorname*{arg\,min}_{\textbf{u}}\operatorname*{arg\,max}_{\textbf{d}}\>\>(\textbf{u},K^{*}\textbf{d})-\frac{1}{2\beta}\|d_{3}\|^{2}_{2}-\alpha\sum_{i=1}^{2}\delta_{B(L^{\infty})}(d_{i}/\alpha).$

Accordingly, the Chambolle-Pock algorithm for this primal-dual problem is
given as:

$\displaystyle\tilde{\textbf{d}}^{k+1}$ $\displaystyle=\textbf{d}^{k}+\sigma
K\bar{\textbf{u}},$ $\displaystyle\textbf{d}_{1,2}^{k+1}$
$\displaystyle=\operatorname*{arg\,min}_{\textbf{d}}\Bigg{\\{}\frac{1}{2}\|\textbf{d}-\tilde{\textbf{d}}_{1,2}^{k+1}\|_{2}^{2}+\alpha\sigma\delta_{B(L^{\infty})}(\textbf{d}/\alpha)\Bigg{\\}},$
$\displaystyle d_{3}^{k+1}$
$\displaystyle=\operatorname*{arg\,min}_{d}\Bigg{\\{}\frac{1}{2}\|d-\tilde{d}_{3}^{k+1}\|_{2}^{2}+\frac{\sigma}{2\beta}\|d\|_{2}^{2}\Bigg{\\}},$
$\displaystyle\tilde{\textbf{u}}^{k+1}$ $\displaystyle=\textbf{u}^{k}-\tau
K^{*}\textbf{d}^{k+1},$ $\displaystyle\bar{\textbf{u}}_{k+1}$
$\displaystyle=2\textbf{u}_{k+1}-\textbf{u}_{k}.$

To derive the optimality condition for the dual variables $d_{3}$, consider
the functional

$J(d_{3})=\frac{1}{2}\int_{\Omega}(d_{3}-\tilde{d}_{3})^{2}+\frac{\sigma}{2\beta}\int_{\Omega}d_{3}^{2}.$

Therefore setting $d_{\theta}J=0$ we get

$d_{3}-\tilde{d}_{3}+\frac{\sigma}{\beta}d_{3}=0.$

Rearranging we get

$d_{3}^{k+1}=\frac{\beta}{\beta+\sigma}\tilde{d}_{3}^{k+1},$

The solution for the indicator function $\delta$ is given by the point-wise
projections of
$\tilde{\textbf{d}}^{k+1},\operatorname*{\textbf{proj}}\nolimits_{\alpha}(\tilde{\textbf{d}}^{k+1})$
onto the unit ball, see [8, 17]. Thus the iterative scheme for the Chambolle-
Pock is given as

$\displaystyle\tilde{\textbf{d}}^{k+1}$ $\displaystyle=\textbf{d}^{k}+\sigma
K\bar{\textbf{u}},$ $\displaystyle\textbf{d}_{1,2}^{k+1}$
$\displaystyle=\operatorname*{\textbf{proj}}\nolimits_{\alpha}\Big{(}\tilde{\textbf{d}}_{1,2}^{k+1}\Big{)},$
$\displaystyle d_{3}^{k+1}$
$\displaystyle=\frac{\beta}{\beta+\sigma}\tilde{d}_{3}^{k+1},$
$\displaystyle\tilde{\textbf{u}}^{k+1}$ $\displaystyle=\textbf{u}^{k}-\tau
K^{*}\textbf{d}^{k+1},$ $\displaystyle\bar{\textbf{u}}_{k+1}$
$\displaystyle=2\textbf{u}_{k+1}-\textbf{u}_{k}.$

###  4.2 Optimality Condition for Our Model $(M2)$

In this case we have

$G(\textbf{u})=\frac{1}{2}\int_{\Omega}(f_{t}+\nabla
f\cdot\textbf{u})^{2},\quad
F(K\textbf{u})=\frac{1}{2}\int_{\Omega}\phi(f,\nabla
f)(\nabla_{H}\cdot\textbf{u})^{2}+\sum_{i=1}^{2}\int_{\Omega}|\nabla u_{i}|.$

The Operator $K$ is given as

$K\textbf{u}=\begin{bmatrix}\nabla&0\\\ 0&\nabla\\\
\phi\partial_{y}&-\phi\partial_{x}\end{bmatrix}\begin{bmatrix}u_{1}\\\\[11.38092pt]
u_{2}\end{bmatrix}.$

Therefore,

$K^{*}\textbf{d}=-\begin{bmatrix}\nabla\cdot&0&\partial_{y}(\phi\cdot)\\\
0&\nabla\cdot&-\partial_{x}(\phi\cdot)\end{bmatrix}\begin{bmatrix}d_{1}\\\
d_{2}\\\ d_{3}\end{bmatrix}.$

As before, the convex conjugate $F^{*}(\textbf{d})$ is computed as

$F^{*}(\textbf{d})=\frac{1}{2}\|d_{3}\|^{2}_{2}+\alpha\sum_{i=1}^{2}\delta_{B(L^{\infty})}(d_{i}/\alpha),$

Thus the primal-dual formulation is given as

$\operatorname*{arg\,min}_{\textbf{u}}\operatorname*{arg\,max}_{\textbf{d}}\>\>\frac{1}{2}\int_{\Omega}(f_{t}+\nabla
f\cdot\textbf{u})^{2}+(\textbf{u},K^{*}\textbf{d})-\frac{1}{2\beta}\|d_{3}\|^{2}_{2}-\alpha\sum_{i=1}^{2}\delta_{B(L^{\infty})}(d_{i}/\alpha).$

Accordingly, the Chambolle-Pock algorithm for this primal-dual problem is
given as:

$\displaystyle\tilde{\textbf{d}}^{k+1}$ $\displaystyle=\textbf{d}^{k}+\sigma
K\bar{\textbf{u}},$ $\displaystyle\textbf{d}_{1,2}^{k+1}$
$\displaystyle=\operatorname*{arg\,min}_{\textbf{d}}\Bigg{\\{}\frac{1}{2}\|\textbf{d}-\tilde{\textbf{d}}_{1,2}^{k+1}\|_{2}^{2}+\alpha\sigma\delta_{B(L^{\infty})}(\textbf{d}/\alpha)\Bigg{\\}},$
$\displaystyle d_{3}^{k+1}$
$\displaystyle=\operatorname*{arg\,min}_{d}\Bigg{\\{}\frac{1}{2}\|d-\tilde{d}_{3}^{k+1}\|_{2}^{2}+\frac{\sigma}{2\beta}\|d\|_{2}^{2}\Bigg{\\}},$
$\displaystyle\tilde{\textbf{u}}^{k+1}$ $\displaystyle=\textbf{u}^{k}-\tau
K^{*}\textbf{d}^{k+1},$ $\displaystyle\textbf{u}^{k+1}$
$\displaystyle=\operatorname*{arg\,min}_{\textbf{u}}\Bigg{\\{}\frac{1}{2}\|\textbf{u}-\tilde{\textbf{u}}^{k+1}\|_{2}^{2}+\frac{\tau}{2}\int_{\Omega}(f_{t}+\nabla
f\cdot\textbf{u})^{2}\Bigg{\\}},$ $\displaystyle\bar{\textbf{u}}_{k+1}$
$\displaystyle=2\textbf{u}_{k+1}-\textbf{u}_{k}.$

The optimality conditions for the dual variables follow directly from above.
For the primal variable u the optimality condition can be obtained directly by
a quadratic minimization, see [17] for more details. The equations are given
as

$\displaystyle(1+\tau f_{x}^{2})u_{1}+\tau f_{x}f_{y}u_{2}$
$\displaystyle=\tilde{u}^{k+1}_{1}-\tau f_{x}f_{t},$ $\displaystyle\tau
f_{x}f_{y}u_{1}+(1+\tau f_{y}^{2})u_{2}$
$\displaystyle=\tilde{u}^{k+1}_{2}-\tau f_{y}f_{t}.$

Thus the iterative scheme for the Chambolle-Pock is given as

$\displaystyle\tilde{\textbf{d}}^{k+1}$ $\displaystyle=\textbf{d}^{k}+\sigma
K\bar{\textbf{u}},$ $\displaystyle\textbf{d}_{1,2}^{k+1}$
$\displaystyle=\operatorname*{\textbf{proj}}\nolimits_{\alpha}\Big{(}\tilde{\textbf{d}}_{1,2}^{k+1}\Big{)},$
$\displaystyle d_{3}^{k+1}$
$\displaystyle=\frac{\beta}{\beta+\sigma}\tilde{d}_{3}^{k+1},$
$\displaystyle\tilde{\textbf{u}}^{k+1}$ $\displaystyle=\textbf{u}^{k}-\tau
K^{*}\textbf{d}^{k+1},$ $\displaystyle\textbf{u}^{k+1}$
$\displaystyle=\Bigg{(}\frac{b_{1}c_{3}-c_{2}b_{2}}{c_{1}c_{3}-c_{2}^{2}},\frac{b_{2}c_{1}-c_{2}b_{1}}{c_{1}c_{3}-c_{2}^{2}}\Bigg{)},$
$\displaystyle\bar{\textbf{u}}_{k+1}$
$\displaystyle=2\textbf{u}_{k+1}-\textbf{u}_{k}.$

where $c_{1},c_{2},c_{3}$ are the elements of the coefficient matrix given by
$c_{1}=1+\tau f_{x}^{2},c_{2}=\tau f_{x}f_{y},c_{3}=1+\tau f_{y}^{2}$,
$b_{1},b_{2}$ are the right hand side values given by
$b_{1}=\tilde{u}^{k+1}_{1}-\tau f_{x}f_{t},b_{2}=\tilde{u}^{k+1}_{2}-\tau
f_{y}f_{t}$. In the next section we will look at the numerical discretization
and other implementation details.

##  5 Results

Having obtained the Chambolle-Pock algorithm for solving the saddle-point
problem, we now look at the implementation details. Algorithm 1 shows the
Chambolle-Pock algorithm for our nonlinear constraint-based refinement model.

Algorithm 1

1:Initialize $\tau,\sigma\leftarrow 1/\sqrt{8},1/\sqrt{8}$

2:Initialize $\textbf{u}^{0}\leftarrow$ HS($f_{1},f_{2}$),
$\textbf{d}^{0}\leftarrow$ 0

3:Initialize operator $K$

4:repeat

5: $\textbf{u}_{\text{old}}\leftarrow\textbf{u}$

6: Build Matrix $K$

7: $\tilde{\textbf{d}}\leftarrow\textbf{d}+\sigma K\bar{\textbf{u}}$

8:
$d_{1,2}\leftarrow\operatorname*{\textbf{proj}}_{\sigma/\alpha}(\tilde{d}_{1,2})$

9: $d_{3}\leftarrow\frac{\beta}{\beta+\sigma}\tilde{d}_{3}$

10: Build Matrix $K^{*}$

11: $\tilde{\textbf{u}}\leftarrow\textbf{u}-\tau K^{*}\textbf{d}$

12: $\textbf{u}\leftarrow\tilde{\textbf{u}}$

13: $\bar{\textbf{u}}\leftarrow 2\textbf{u}-\textbf{u}_{\text{old}}$

14:until convergence

As mentioned previously this model works in two phases wherein the first phase
we obtain a crude-pixel correspondence and subsequently refine this estimate
in the next phase driven by additional constraints. The initial Horn and
Schunck flow was computed using the Chambolle-Pock algorithm, see [10, 17].
Here we observed that using a forward difference scheme for both spatial and
temporal image derivatives $f_{x},f_{y}$ and $f_{t}$ respectively does not
yield a stable discretization. Instead, a forward difference scheme for
$f_{t}$ and a central difference scheme for $f_{x},f_{y}$ does yield a stable
numerical scheme. In the next step, the operator matrix $K$ is constructed for
updating the dual variables $d_{1},d_{2},d_{3}$ shown in steps 9 and 10. This
requires solving two sub-problems, one for $d_{1},d_{2}$ and the other for
$d_{3}$.

Now $\nabla u_{i}=(u_{{i}_{x}},u_{{i}_{y}}),i=1,2$. The associated dual
variable is $d_{i}=(d_{i,1},d_{i,2})$. The primal formulation had a total-
variation regularization. Accordingly,

$|\nabla u_{i}|_{L^{1}}=|u_{{i}_{x}}|+|u_{{i}_{y}}|.$

Thus the associated dual norm for the variable $d_{i}$ gives

$\|d_{i}\|_{L^{\infty}}=\max\\{|d_{i,1}|,|d_{i,2}|\\}.$

The solution for this minimization is the point-wise projection onto the unit
ball corresponding to the dual norm. As shown previously, the convex conjugate
of the total variation term is the indicator function
$\delta_{L^{\infty}}(d_{i}/\alpha)$. The associated convex set is defined by

$\\{d_{i}:\|d_{i}/\alpha\|\leq 1\\}=\\{d_{i}:\|d_{i}\|\leq\alpha\\},i=1,2.$

Thus the dual update for $d_{1,2}$ can be obtained by the point-wise
projection of $\tilde{d}_{1,2}$ onto $[-\alpha,\alpha]$ given as

$d_{1,2}=\operatorname*{\textbf{proj}}\nolimits_{\sigma/\alpha}(\tilde{d}_{1,2})=\min(\alpha,\max(-\alpha,\tilde{d}_{1,2})).$

The sub-problem for $d_{3}$ is a linear quadratic minimization problem as
discussed in the previous section. In the next step the adjoint operator
$K^{*}$ is constructed to update the primal variable u. The subsequent over-
relaxation step $\bar{\textbf{u}}\leftarrow
2\textbf{u}-\textbf{u}_{\text{old}}$ is a particular case for $\theta=1$ in
[10] for easier estimates of the convergence. The algorithm is further
simplified if the regularization term is linear. In this case the only
difference that occur is the updation of the dual variable $d_{1,2}$ leading
to the optimality condition

$d_{1,2}=\frac{\alpha}{\alpha+\sigma}\tilde{d}_{1,2}.$

The stopping criterion is determined by computing the normalized error from
the primal and dual residues. This error metric was introduced by the authors
in [18] and is numerically less expensive. Let
$\textbf{u}^{(k)},\textbf{d}^{(k)}$ be the primal and the dual updates after
$k$ iterations respectively. Then the primal and dual residues at the $k^{th}$
iteration are computed by the formula:

$\displaystyle p_{\text{res}}^{(k)}:$
$\displaystyle=\Bigg{|}\frac{\textbf{u}^{(k)}-\textbf{u}^{(k+1)}}{\tau}-K^{*}(\textbf{d}^{(k)}-\textbf{d}^{(k+1)})\Bigg{|},$
$\displaystyle d_{\text{res}}^{(k)}:$
$\displaystyle=\Bigg{|}\frac{\textbf{d}^{(k)}-\textbf{d}^{(k+1)}}{\sigma}-K(\textbf{u}^{(k)}-\textbf{u}^{(k+1)})\Bigg{|}.$

Therefore the normalized error at $k^{th}$ step is obtained as:

$e^{(k)}=\frac{p_{\text{res}}^{(k)}+d_{\text{res}}^{(k)}}{\mu(\Omega)},$

where $\mu(\Omega)$ refers to the measure of the domain $\Omega$. Chambolle
and Pock [10] also showed that the convergence criterion is fulfilled when
$\tau\sigma\|K\|^{2}<1,\theta=1$. Thus $\tau$ and $\sigma$ needs to be chosen
accordingly. An optimal numerical upper-bound was obtained by Chambolle [9]
which satisfies the above criterion. Accordingly we set
$\tau=\sigma=1/\sqrt{8}$.

The Chmabolle-Pock algorithm for the angular accuracy model follows in a
similar manner. The main difference lies in the primal update step 13 because
of the explicit presence of the data term in the functional. The subroutine
_solvePrimal_ updates the primal variable u by solving a quadratic
minimization problem discussed previously leading to the following update
step,

$\textbf{u}^{k+1}=\Bigg{(}\frac{b_{1}c_{3}-c_{2}b_{2}}{c_{1}c_{3}-c_{2}^{2}},\frac{b_{2}c_{1}-c_{2}b_{1}}{c_{1}c_{3}-c_{2}^{2}}\Bigg{)},$

where $c_{1}=1+\tau f_{x}^{2},c_{2}=\tau f_{x}f_{y},c_{3}=1+\tau f_{y}^{2}$,
$b_{1}=\tilde{u}^{k+1}_{1}-\tau f_{x}f_{t},b_{2}=\tilde{u}^{k+1}_{2}-\tau
f_{y}f_{t}$. We now show the results for obtained by implementing the
algorithm. The first sequence is the Oseen vortex pair. For more details on
the sequence we refer to [25].

(a)

(b)

Figure 1: Oseen Vortex Pair [25]. Figure 2: Velocity magnitude plot for the
Oseen vortex pair with $\alpha=0.1,\beta=0.01,\text{iter}=50$.

Figure (2) shows the velocity magnitude plot obtained for the oseen vortex
pair. The algorithm produces dense flow fields while correctly estimating the
vortex cores. We now empirically demonstrate the rate of convergence of both
models.

| $\epsilon=0.1$ | $\epsilon=0.01$
---|---|---
Model (M1) | Model (M2) | Model (M1) | Model (M2)
Oseen Vortex Pair | 78 | 627 | 755 | 53753
Cloud Sequence | 20 | 99 | 450 | 16068
Sphere Sequence | 10 | 316 | 561 | 18009
Hydrangea | 103 | 346 | 937 | 12574
Rubberwhale | 42 | 351 | 617 | 13590
Table 3: Total number of iterations required by the algorithms to reach the
threshold of $\epsilon$.

Table (3) shows the number of iterations required by the algorithm to reach
the error threshold of $\epsilon$. By an order $O(1/N)$ convergence we mean
that the number of iterations required to reach a tolerance $\epsilon$ is
$O(1/\epsilon)$. This is validated from the above table. For Model (M1), for
$\epsilon=0.1$, the number of iterations required to reach the threshold of
$0.1$ is a multiple of 10. For Model (M2) it requires a multiple of $10^{2}$
iterations. The table also shows that roughly it requires a multiple of $100$
iterations for the Model (M1) algorithm to reach the threshold of
$\epsilon=0.01$ and a multiple of $100^{2}$ iterations for Model (M2). The
reason for this efficiency can be explained from the fact that in phase 1, the
Horn and Schunck initialization brings the solution within a close error
range. As a result in the second phase we observe a $O(1/N)$ convergence as
mentioned in [10].

###  5.1 Modern Implementation Principles

Recent developments in optical flow computation reveal that the flow estimates
can be significantly improved by incorporating certain principles. To
accommodate these principles in our framework Algorithm 1 is suitably modified
to make it a part of a larger implementation procedure.

The computation of image derivatives follows a weighted averaging principle
[29, 7]. The current flow estimates are used to warp the second image towards
the first using bi-cubic interpolation. The time derivative is the difference
between the first image and the warped image. The spatial derivatives are
obtained as a weighted average of the first image and the warped image. The
weight coefficient is called the blending ratio was chosen between 0 and 1.

Figure 3: Estimated flow field for the _rubberwhale_ sequence with
$\alpha=10,\beta=1$.

To account for large displacements of pixel motions, a coarse-to-fine
pyramidal scheme is employed [32, 7, 17]. The flow field is first computed at
the coarsest level. This estimate is upsampled to the next finer level via
interpolation and is used to warp the second image towards the first image.
The flow increments at this finer level are then computed between the first
image and the warped image. This process continues till the finest resolution
level is reached. At each pyramid level, 10 warping steps are performed. after
each warping iteration, a $5\times 5$ median filter is applied on the flow
estimates to remove outliers. Figure (3) shows the obtained flow field for
nonlinear refinement using Algorithm 1 for the rubberwhale sequence.

| Model 1† | Model 2†
---|---|---
AAE | EPE | AAE | EPE
Linear Refinement | 3.412 | 0.105 | 3.410 | 0.105
Nonlinear Refinement | 3.397 | 0.104 | 3.355 | 0.103
Table 4: Comparison of the Average Angular Error (AAE) and End Point Error
(EPE) for _rubberwhale_ sequence. († refers to the algorithm incorporated with
the modern implementation principles)

An improvement is seen in the average angular error (AAE) and the end-point
error (EPE) for the nonlinear refinement compared to the linear refinement for
the _rubberwhale_ sequence as shown in Table (4). The improvement however is
not significant which indicates that incorporating the modern implementation
principles like coarse-to-fine warping, median filtering and so on also
improves the accuracy of the linear refinement process. There are however
image sequences for which edge-information are not well-preserved by the
linear refinement process. We demonstrate this with the sphere sequence [2].

(a)

(b)

Figure 4: Sphere sequence [2]

(a) Linear Refinement

(b) Nonlinear Refinement

Figure 5: Estimated flow fields of the _sphere_ sequence from [2] using
Algorithm 1† with $\alpha=1,\beta=0.1$ († refers to the algorithm incorporated
with the modern implementation principles).

Figure (5) shows the color-coded flow estimate for the _sphere_ sequence using
the Middlebury color coding. The isotropic behaviour is seen in the linear
case because of which the edges are not well preserved.

## Conclusion

In this paper we have proposed two nonlinear variational models for obtaining
an accurate estimation of physics-based flow fields such as rotational and
fluid flow. The first model is a novel two-phase refinement process where in
the first phase a crude estimate is obtained and subsequently refined using
additional constraints in the second phase. We have studied the well-posedness
of this model using an Evolutionary PDE approach. The second model performs
the same refinement using a single functional. We used the first-order primal-
dual Chambolle-Pock algorithm for the numerical implementation of the above
models. We further empirically demonstrated that the two-phase model leads to
a faster convergence rate of the order $O(1/N)$ compared to the second model
which has a convergence rate of the order $O(1/N^{2})$.

######  Acknowledgements.

The authors dedicate this paper to Bhagawan Sri Sathya Sai Baba, Revered
Founder Chancellor, SSSIHL.

## References

  * [1] https://vision.middlebury.edu/flow/data/
  * [2] http://of-eval.sourceforge.net/
  * [3] Aubert G., Deriche R., Kornprobst P.: Computing Optical Flow via Variational Techniques, SIAM Journal of Applied Mathematics, Vol. 60, 156-182 (1999).
  * [4] Barron J.L., Fleet J., Beauchemin S., Performance of Optical Flow Techniques, International Journal of Computer Vision, 12(1), 43-77 (1994).
  * [5] Brezis H., Functional Analysis, Sobolev Spaces and Partial Differential Equations, Springer, New York (2011).
  * [6] Brox T., Bruhn A., Papenberg N., Weickert J., High Accuracy Optical Flow Estimation Based on a Theory of Warping, European Conference of Computer Vision, 25-36 (2004).
  * [7] Bruhn A., Weickert J., Towards Ultimate Motion Estimation: Combining Highest Accuracy with Real-time Performance, Tenth IEEE International Conference on Computer Vision, Vol. 1, 749-755 (2005).
  * [8] Burger M., Dirks H., Frerking L., On Optical Flow Models for Variational Motion Estimation In: M. Bergounioux, G. Peyré, C. Schnörr, J. Caillau & T. Haberkorn (Ed.), Variational Methods: In Imaging and Geometric Control, Berlin, Boston: De Gruyter, 225-251 (2017).
  * [9] Chambolle A., An Algorithm for Total Variation Minimization and Applications, Journal of Mathematical Imaging and Vision, 20, 89-97 (2004).
  * [10] Chambolle A., Pock T., A First-Order Primal-Dual Algorithm for Convex Problems with Applications in Imaging, Journal of Mathematical Imaging and Vision, 40, 120-145 (2011).
  * [11] Cohen I., Nonlinear Variational Method for Optical Flow Computation, Proceedings of the 8th SCIA, 523-530 (1993).
  * [12] Corpetti T., Mémin E., Pérez P.: Estimating Fluid Optical Flow, Proceedings of the 15th International Conference on Pattern Recognition, ICPR-2000, Vol. 3, 1033-1036 (2000).
  * [13] Corpetti T., Heitz D., Arroyo G., Mémin E., Santa-Cruz A.: Fluid Experimental Flow Estimation based on an Optical Flow Scheme, Experiments in Fluids, (2006).
  * [14] Dirks H., Variational Methods for Joint Motion Estimation and Image Reconstruction, Ph.D. Thesis, Wilhems-Universität (2015).
  * [15] Doshi H, N. Uday Kiran: Constraint-Based Refinement of Optical Flow, https://arxiv.org/abs/2011.12267.
  * [16] Doshi H., N. Uday Kiran, A Variational Optical Flow Model for Improving Angular Accuracy, https://arxiv.org/abs/2102.00487.
  * [17] Dirks H., Variational Methods for Joint Motion Estimation and Image Reconstruction, Ph.D. Thesis, Wilhems-Universität (2015).
  * [18] Goldstein T., Esser E., Baraniuk R., Adaptive Primal-Dual Hybrid Gradient Methods for Saddle-Point Problems, arXiv:1305.0546 (2013).
  * [19] Hauer D., Mazon J.: Regularizing effects of homogeneous evolution equations: the case of homogenity of order zero, Journal of Evolution Equations, 19, 965-996 (2019).
  * [20] Heitz D., Mémin E., Schnörr C.: Variational Fluid Flow measurements from Image Sequences: Synopsis and Perspectives, Experiments in Fluids, 48, 369–393 (2010).
  * [21] Hinterberger W., Scherzer O., Schnörr C., Weickert J.: Analysis of Optical Flow Models in the Framework of Calculus of Variations, Numer. Funct. Anal. and Optimiz., 23(1&2), 69-89 (2002).
  * [22] Horn B.K.P., Schunck B.G.: Determining Optical Flow, Artificial Intelligence, Vol. 17, 185-203 (1981).
  * [23] Kumar A., Tannenbaum A., Balas G., Optical flow: A curve evolution approach, IEEE Transactions of Image Processing 5, 598-611 (1996).
  * [24] Martin A., Schiavi E., de Léon S.S.: On 1-Laplacian elliptic equations modeling magnetic resonance imaging Rician denoising, Journal of Mathematical Imaging and Vision, 57:202-224 (2017).
  * [25] Liu T.: OpenOpticalFlow: An Open Source Program for Extraction of Velocity Fields from Flow Visualization Images, Journal of Open Research Software, 5:29 (2017).
  * [26] Liu T., Shen L.: Fluid Flow and Optical Flow, Journal of Fluid Mechanics, 614, 253-291 (2008).
  * [27] Luttman A., Bollt E.M., Basnayake R., Kramer S., Tufillaro N.B.: A Framework for Estimating Potential Fluid Flow from Digital Imagery, Chaos: An Interdisciplinary Journal of Nonlinear Science, Vol. 23, 3 (2013).
  * [28] Schnörr C.: Determining optical flow for irregular domains by minimizing quadratic functionals of a certain class, International Journal of Computer Vision, Vol. 6, 25-38 (1991).
  * [29] Sun D., Roth R., Black M.J., Secrets of Optical Flow Estimation and their Principles, IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2432-2439 (2010).
  * [30] Sun, D., Roth, S., Black, M. J., A Quantitative Analysis of Current Practices in Optical Flow Estimation and the Principles Behind them, International Journal of Computer Vision, 106(2), 115-137 (2014).
  * [31] Zhang X., Burger M., Osher S., A Unified Primal-Dual Algorithm Framework Based on Bregman Iteration, Journal of Scientific Computing, 46, 20-46 (2011).
  * [32] Wedel A., Pock T., Zach C., Cremers D., Bischof H., An Improved Algorithm for TV-$L^{1}$ Optical Flow, Dagstuhl Motion Workshop (2008).
  * [33] Wei L., Agarwal R. P., Wong P.J.Y.: Existence and iterative construction of solutions to non-linear Dirichlet boundary value problems with $p$-Laplacian operator, Complex Variables and Elliptic Equations: An International Journal, 55:5-6, 601-608 (2010).
  * [34] Zach C., Pock T., Bischof H., A Duality Based Approach for Realtime TV-$L^{1}$ Optical Flow, LNCS 4713, 214-223, DAGM (2007).
  * [35] Zhang B., Zhu Z., A Primal-Dual Algorithm Framework for Convex Saddle-Point Optimization, Journal of Inequalities and Applications, 267 (2017).
  * [36] Zhao J., Wang Y., Wang H., Optical flow with Harmonic Constraint and Oriented Smoothness, Sixth International Conference on Image and Graphics, 94-99 (2011).