# Real-time Structure Flow

Juan David Adarve
Systems Theory and Robotics Group
Australian National University
ACT, 2601, Australia
<EMAIL_ADDRESS>
& Robert Mahony
Systems Theory and Robotics Group
Australian National University
ACT, 2601, Australia
<EMAIL_ADDRESS>

###### Abstract

This article introduces the _structure flow_ field; a flow field that can
provide high-speed robo-centric motion information for motion control of
highly dynamic robotic devices and autonomous vehicles. Structure flow is the
_angular 3D velocity_ of the scene at a given pixel. We show that structure
flow posses an elegant evolution model in the form of a Partial Differential
Equation (PDE) that enables us to create dense flow predictions forward in
time. We exploit this structure to design a predictor-update algorithm to
compute structure flow in real time using image and depth measurements. The
prediction stage takes the previous estimate of the structure flow and
propagates it forward in time using a numerical implementation of the
structure flow PDE. The predicted flow is then updated using new image and
depth data. The algorithm runs up to 600 Hz on a Desktop GPU machine for
512x512 images with flow values up to 8 pixels. We provide ground truth
validation on high-speed synthetic image sequences as well as results on real-
life video on driving scenarios.

##  1 Introduction

Reliable and robust motion control is a central component of all autonomous
robotic systems. As robotic vehicles become more mobile and begin to operate
in highly unstructured and dynamic environments, such as urban streetscapes,
the requirements on the perception subsystem of the robot are correspondingly
more demanding. Classically, most autonomous robot use a combination of
Simultaneous Localisation And Mapping (SLAM) [CCC+16] along with dedicated
sensor-control feedback for motion control. More specifically, SLAM algorithms
are used to provide large scale (sparse) environment representations
appropriate for path planning and decision processes, while real-time motion
control loops are closed using high-bandwidth direct sensor feedback. Sensor
systems capable of providing a dense, high-bandwidth, 3D spatiotemporal
(structure and motion) representation of an outdoor environment, such as
LIDAR, are expensive, and often bulky and heavy. In contrast, high frame-rate
vision systems offer a rich, dense, and dynamic sensor modality for low-weight
and low-price. Vision systems, however, require sophisticated post processing
to extract spatiotemporal information. Modern dense visual SLAM
implementations [MCR15, WSMG+16] demonstrate the potential of vision systems
to provide situational awareness to autonomous robotic vehicles, although
these algorithms still suffer from limitations such as; (relatively) low
update rates, static environment assumptions and limited map size.

We propose the following aspirational goal requirements for a modern vision
based spatiotemporal perception system capable of providing real-time
situational awareness for a robot moving in a highly unstructured and dynamic
environment such as a streetscape:

  * •

Egocentric (robo-centric) map representation.

  * •

Dense range estimation to the environment. _Benchmark: angular resolution 0.1
degrees._

  * •

Dense 3D-velocity (structure flow) estimation. _Benchmark: resolution
0.05rad.s -1, (for example, 0.05m.s-1 error in 1m.s-1 velocity at 1m
distance)._

  * •

Real-time (high-bandwdith) updates. _Benchmark: updates at 300 Hz._

Such a system would provide sufficient information for robust local path
planning and obstacle avoidance. High level path planning and decision
processes could be undertaken separately based on more traditional large scale
SLAM algorithms. The present work also continues development (see also
[AM16a]) of a new paradigm of vision processing algorithms based on real-time
filtering concepts, that we believe will provide an important tool in the
future development of autonomous robotic systems.

The technical contributions of this article are twofold. First, we introduce a
new robo-centric spatiotemporal representation of motion, the _structure
flow_. The structure flow is a 3-vector assigned to each ‘pixel’ in the image
comprised of the three-dimensional Euclidean velocity between the robot and
the environment (the scene flow [VBR+99]) scaled by the inverse range of the
scene. Analogously to optic flow, structure flow encodes information about the
visually perceived motion and ego-motion of the scene. However, structure flow
generalises optic flow by including a component in the normal direction of the
image that is closely linked to visual divergence. Structure flow will have a
similar utility as optic flow for control of robotic vehicles [HM02, RF05,
HRHM08, Sri11] with the added advantage that the normal component of structure
flow directly yields the divergent flow properties used in many of the
obstacle avoidance algorithms [MB12]. A key advantage of structure flow is
that, for reasonable assumptions on the environment and robot motion, it is
possible to derive a dynamic model of the field evolution in the form of a
Partial Differential Equation (PDE). Moreover, the structure flow PDE follows
the same pattern as the partial differential equation modeling the evolution
of the inverse depth field [BR09, ZAR12]. We exploit this to forward propagate
estimates of the structure flow and the inverse depth in image coordinates to
make predictions of future motion and structure of the environment.

The second technical contribution of the paper lies in proposing an iterative
algorithm to compute structure flow in real-time from stereo image data. The
approach taken is a highly distributed predictor-update algorithm that we
implement on a GPU. The prediction step is a numerical implementation of a PDE
integration scheme that propagates the current estimate of the structure flow
forward in time. Error in the structure flow estimate can then be estimated
across a sequence of temporal frames, and we use a simple least-squares
regularised update to iteratively correct errors in the structure flow
estimate. The resulting algorithm runs on a Nvidia GTX-780 GPU, processing
$512\times 512$ at approximately 600 Hz for flow vectors up to 8 pixels in
magnitude. The relative performance of the algorithm far outperforms (10-20
times faster) all scene flow algorithms in the literature in terms of
processing speed. Although it is less accurate than some classical algorithms,
its advantages in speed, dense estimation, and robustness make it highly
suitable for real-world mobile robotic applications.

##  2 Related Work

The estimation of motion from image sequences is a heavily studied problem in
computer vision and robotics [Bar94, BSL+11]. The importance of these motion
fields in the control of mobile robotic vehicles is equally important, and has
been studied for more than thirty years in applications such as visual control
of aerial vehicles [RF05, Sri11, HHMR12] and visual odometry [SF11].

The most common motion field considered is the _optical flow_ field, defined
as the displacement of points in the image plane [Bar94]. Algorithms for the
computation of flow can be classified as local or global based methods, from
which the seminal works of Lucas and Kanade (local) [LK81] and Horn and
Schunck [HS81] (global) marked a distinction between local and global methods.
Recent algorithms can be found in benchmark datasets such as Middlebury
[BSL+11] and Kitti [GLSU13]. In robotics, an underlying requirement for flow
algorithms is their capability to run at real-time frequencies. Real time
methods [PLBC14, AM16a] are usually local-based methods, which can be
efficiently implemented using GPU or FPGA hardware. The vast majority of
optical flow algorithms developed in by the computer-vision community in the
last ten years are focused on accuracy and not on computational efficiency
[BSL+11, GLSU13].

A more recent visual flow is _scene flow_ , defined as the three-dimensional
velocity field of points moving in the world [VBR+99] expressed in the ‘image’
frame of the robot. Stereoscopic algorithms use image sequences captured using
a stereo camera to estimate the 3D scene flow. Early work by Patras et al.
[PAT96] formulates the joint estimation of disparity and motion (scene flow)
from stereo sequences. Huguet and Devernay [HD07] proposed a variational
approach for the joint computation of disparity and optical flow. This
formulation of optical flow plus temporal disparity change is known as
_disparity flow_ [GY06]. Wedel et al. [WBV+11] decouple the variational
formulation of Huguet et. al. to compute flow and disparity in two separate
stages while preserving stereo constraints. Vogel et. al. [VSR15] compute
scene flow by means of a piece-wise planar segmentation of stereo images for
which a 3D rotation and translation (scene flow) is computed. The accuracy of
recent two-pair stereo based algorithms is compiled on the Kitti scene flow
dataset by Menze and Geiger [MG15]. Algorithms using RGB-D sensors to compute
scene flow combine monocular images and depth measurements to estimate the 3D
velocity field. Letouzey et. al. [LPB11] combine sparse features points and
dense smoothness constraints to estimate dense scene flow from RGB-D data.
Hadfield and Bowden [HB11] use a particle filter formulation to track points
in the scene. Quiroga et. al. [QDC12] use brightness and depth data in a
Lukas-Kanade type tracking framework to compute scene flow parametrized as a
rotation plus translation transform.

If only depth/range measurements are available it is still possible to compute
_range flow_. Yamamoto et. al. derived the range flow constraint [YBBR93], a
partial differential equation to model the change of depth. Differential
methods to compute range flow suffer from a 3D version of the aperture problem
similar to that found in optical flow methods [SJB02]. Herbst et. al. [HRF13]
use this differential constraint on depth and combined with color image to
derive a variational method to compute scene flow following a formulation
similar to Brox et. al. [BBPW04] for computing optical flow.

In terms of real-time capable dense 3D motion estimation algorithms, most
methods use Graphics Processing Units (GPU) for implementing per-pixel level
operations in parallel. Gong [Gon09] reports 12 Hz optical and disparity flow
estimation on QVGA image resolution ($320\times 240$). Rabe et. al. [RMWF10]
reported 25 Hz frequency on VGA ($640\times 480$) image sequences using GPU.
The algorithm of Wedel et. al. [WBV+11] reports frame rates of 20 Hz also at
QVGA resolution. RGB-D flow by Herbst et. al. [HRF13] runs between 8-30 Hz at
QVGA depending on the amount of smoothing required.

##  3 Structure Flow

We use a spherical camera model with the embedded coordinates
$S^{2}=\\{\boldsymbol{\eta}\in\mathbb{R}^{3}\;|\;\left\|\boldsymbol{\eta}\right\|=1\\}$
representing points on the sphere. The tangent space $T_{\eta}S^{2}$
associated to each point $\boldsymbol{\eta}$ is defined as

$T_{\eta}S^{2}=\\{\boldsymbol{\mu}\in\mathbb{R}^{3}\;|\;\langle\boldsymbol{\eta},\boldsymbol{\mu}\rangle=0\\}$
(1)

and the projection matrix
$\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}:\mathbb{R}^{3}\rightarrow
T_{\eta}S^{2}$ is the orthogonal projection onto $T_{\eta}S^{2}$

$\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}=(\mathbf{I}_{3}-\boldsymbol{\eta}\boldsymbol{\eta}^{\top}).$
(2)

###  3.1 Optical, scene and structure flow

Figure 1: Scene $\mathbf{v}(\boldsymbol{\eta})$, structure
$\mathbf{w}(\boldsymbol{\eta})$ and optical flow
$\boldsymbol{\Phi}(\boldsymbol{\eta})$.

Consider Figure 1 where the camera is located in the centre of the body-fixed
frame. Let $\mathbf{v}_{c}$ and $\boldsymbol{a}_{c}$ denote the linear
velocity and acceleration respectively of the camera, taken with respect to an
inertial frame, and expressed in the camera frame. Let $\Omega$ denote the
angular velocity of the camera expressed in the camera frame and let
$\Omega_{\times}$ be the skew-symmetric matrix such that $\Omega\times
v=\Omega_{\times}v$ for any vector $v$. For the present derivation we will
assume that $\Omega$ and $\boldsymbol{a}_{c}$ are approximately constant. This
is a reasonable assumption given the very high frame rates of the algorithms
that we propose. With this assumption, the second order kinematics of the
camera are

$\displaystyle\frac{\mathrm{d}\mathbf{v}_{c}}{\mathrm{d}t}$
$\displaystyle=-\Omega_{\times}\mathbf{v}_{c}+\boldsymbol{a}_{c}$ (3)
$\displaystyle\frac{\mathrm{d}\Omega}{\mathrm{d}t}$ $\displaystyle=0$ (4)

Both $\Omega$ and $\boldsymbol{a}_{c}$ are readily available from a typical
inertial measurement unit and we assume they are known.

Let $\mathbf{x}=\mathbf{x}(t)\in\mathbb{R}^{3}$ denote the position of a point
in the scene at time $t$ in body fixed frame (Figure 1). Let
$\mathbf{v}_{\mathbf{x}}$ denote the velocity of the point $\mathbf{x}$ with
respect to the inertial frame and expressed in the camera body fixed frame. We
assume that the scene is moving with a constant velocity with respect to the
inertial frame. With these assumptions, the first and second order kinematics
of a point $\mathbf{x}$ in the scene expressed in the camera frame are:

$\displaystyle\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}$
$\displaystyle=-\Omega_{\times}\mathbf{x}+\mathbf{v}_{\mathbf{x}}-\mathbf{v}_{c},$
(5) $\displaystyle\frac{\mathrm{d}\mathbf{v}_{\mathbf{x}}}{\mathrm{d}t}$
$\displaystyle=-\Omega_{\times}\mathbf{v}_{\mathbf{x}}.$ (6)

The approach taken in this paper is to parameterize the environment
surrounding the robot by direction $\boldsymbol{\eta}\in S^{2}$. Thus, at time
$t$, the first visible point in the environment in direction
$\boldsymbol{\eta}$ is denoted $\mathbf{x}(\boldsymbol{\eta},t)$. The range or
depth map around the robot is a single valued depth field
$0pt:S^{2}\times\mathbb{R}\rightarrow\mathbb{R}$,

$0pt(\boldsymbol{\eta},t)=\langle\mathbf{x}(\boldsymbol{\eta},t),\mathbf{x}(\boldsymbol{\eta},t)\rangle^{\nicefrac{{1}}{{2}}}.$
(7)

Although range $0pt(\boldsymbol{\eta},t)$ appears to be the more intuitive
variable, it turns out that it is more natural in the context of the
mathematical development to use an inverse ratio of range
$\rho:S^{2}\times\mathbb{R}\rightarrow\mathbb{R}$,

$\rho(\boldsymbol{\eta},t)=\frac{0pt_{\text{ref}}}{\lambda(\boldsymbol{\eta},t)},$
(8)

where $0pt_{\text{ref}}$ is the reference range to the spherical image plane.
We will choose $0pt_{\text{ref}}=1$m for simplicity and suppress the
$0pt_{\text{ref}}$ notation in the remainder of the paper. However, one should
still think of $\rho(\boldsymbol{\eta},t)$ as a dimensionless ratio of ranges
related to the perceived visual angle of an object at range
$0pt(\boldsymbol{\eta},t)$ as observed on the sphere of radius
$0pt_{\text{ref}}$. In particular, $\rho(\boldsymbol{\eta},t)$ is
intrinsically a visual measure or angle, whereas $0pt(\boldsymbol{\eta},t)$ is
extrinsically related to the physical scene and has units of metres.

The scene flow is defined as the 3D velocity of the environment relative to
the camera [VBR+99] of a point $\mathbf{x}(\boldsymbol{\eta},t)$ in the scene.
To compute scene flow, consider a solution $\mathbf{x}(\tau)$ of (5) for
$\tau\in(t-\delta,t+\delta)$, $\delta>0$, in a small time interval around an
‘initial condition’ $\mathbf{x}(t)=\mathbf{x}(\boldsymbol{\eta},t)$ in the
desired direction $\boldsymbol{\eta}$ at time $t$. Choose
$\boldsymbol{\eta}(\tau)$ to evolve such that
$\mathbf{x}(\tau)=\mathbf{x}(\boldsymbol{\eta}(\tau),\tau)$ is the point in
the scene associated with the evolution of $\mathbf{x}(\tau)$ and note that
$\boldsymbol{\eta}(t)=\boldsymbol{\eta}$; that is $\boldsymbol{\eta}(\tau)$
evaluated at $\tau=t$ is $\boldsymbol{\eta}$. The scene flow is the vector
field $\mathbf{v}:S^{2}\times\mathbb{R}\rightarrow\mathbb{R}^{3}$ defined as

$\mathbf{v}(\boldsymbol{\eta},t)=\left.\frac{\mathrm{d}\mathbf{x}(\boldsymbol{\eta}(\tau),\tau)}{\mathrm{d}\tau}\right|_{\tau=t}=-\Omega_{\times}\mathbf{x}+\mathbf{v}_{\mathbf{x}}-\mathbf{v}_{c}.$
(9)

It is important to note that, whereas (5) tracks a particle in the
environment, the scene-flow is a field that assigns values for all points
parameterized on the sphere. In particular, $\boldsymbol{\eta}$ and $t$ are
independent variables in the expression (9) for
$\mathbf{v}(\boldsymbol{\eta},t)$.

A more common notation in the literature used for the total derivative is

$\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}(t)=\left.\frac{\mathrm{d}\mathbf{x}(\boldsymbol{\eta}(\tau),\tau)}{\mathrm{d}\tau}\right|_{\tau=t}$

where the underlying construction of the segment of trajectory
$\mathbf{x}(\tau)$ is understood. We have chosen to use a more explicit
notation in this paper to provide more clarity in the derivations.

We define the structure flow field
$\mathbf{w}:S^{2}\times\mathbb{R}\rightarrow\mathbb{R}^{3}$ as the three-
dimensional vector field consisting of the scene flow scaled by the inverse
depth

$\displaystyle\mathbf{w}(\boldsymbol{\eta},t)$
$\displaystyle=\frac{1}{0pt(\boldsymbol{\eta},t)}\mathbf{v}(\boldsymbol{\eta},t)$
$\displaystyle=-\Omega_{\times}\boldsymbol{\eta}+\frac{1}{0pt(\boldsymbol{\eta},t)}(\mathbf{v}_{\mathbf{x}}-\mathbf{v}_{c})$
(10)

Note that while scene flow has units of m.s-1, structure flow has units of
rad.s-1, where the angle is related to the ratio of distances, as would be
expected of an image based flow measure. We also introduce notation to
decompose structure flow into rotational and stabilised (linear) components,
$\mathbf{w}=\mathbf{w}_{r}+\mathbf{w}_{s}$, where

$\displaystyle\mathbf{w}_{r}(\boldsymbol{\eta},t)$
$\displaystyle=-\Omega_{\times}\boldsymbol{\eta},$ (11)
$\displaystyle\mathbf{w}_{s}(\boldsymbol{\eta},t)$
$\displaystyle=\frac{1}{0pt(\boldsymbol{\eta},t)}(\mathbf{v}_{\mathbf{x}}-\mathbf{v}_{c}).$
(12)

In order to compute evolution equations for structure flow it is necessary to
compute the total time derivative of the flow. We first compute total time
derivatives of $\boldsymbol{\eta}(\tau)$ and $0pt(\tau)$ defined along
trajectories $\mathbf{x}(\tau)=\mathbf{x}(\boldsymbol{\eta}(\tau),\tau)$
induced by the equations of motion (3), (4), (5) and (6). For depth
$0pt(\tau)$, one has

$\displaystyle\left.\frac{\mathrm{d}0pt(\tau)}{\mathrm{d}\tau}\right|_{\tau=t}$
$\displaystyle=\left.\frac{\mathrm{d}}{\mathrm{d}\tau}\langle\mathbf{x}(\tau),\mathbf{x}(\tau)\rangle^{\nicefrac{{1}}{{2}}}\right|_{\tau=t}$
$\displaystyle=\left.\frac{\mathbf{x}(\tau)^{\top}}{\langle\mathbf{x}(\tau),\mathbf{x}(\tau)\rangle^{\nicefrac{{1}}{{2}}}}\frac{\mathrm{d}\mathbf{x}(\tau)}{\mathrm{d}\tau}\right|_{\tau=t}$
$\displaystyle=\left.\langle\boldsymbol{\eta}(\tau),\frac{\mathrm{d}\mathbf{x}(\tau)}{\mathrm{d}\tau}\rangle\right|_{\tau=t}$
$\displaystyle=0pt(\boldsymbol{\eta},t)\langle\boldsymbol{\eta},\mathbf{w}(\boldsymbol{\eta},t)\rangle,$
(13)

where the last line follows by substituting for structure flow (10) at
$\tau=t$. A similar derivation for the inverse depth $\rho(t)$ yields

$\displaystyle\left.\frac{\mathrm{d}\rho(\tau)}{\mathrm{d}\tau}\right|_{\tau=t}$
$\displaystyle=-\rho(\boldsymbol{\eta},t)\langle\boldsymbol{\eta},\mathbf{w}(\boldsymbol{\eta},t)\rangle$
(14)

We note that the right hand side of both (13) and (14) depend only on the
field variables $\mathbf{w}$, $0pt$ and $\rho$ defined at
$(\boldsymbol{\eta},t)$ and not on the trajectory $\mathbf{x}(\tau)$ from
which they are derived.

Next, consider the rate of change of the direction
$\boldsymbol{\eta}(\tau)=\nicefrac{{\mathbf{x}(\tau)}}{{0pt(\tau)}}$. One has

$\displaystyle\left.\frac{\mathrm{d}\boldsymbol{\eta}(\tau)}{\mathrm{d}\tau}\right|_{\tau=t}$
$\displaystyle=\left.\frac{\mathrm{d}}{\mathrm{d}\tau}\left(\frac{\mathbf{x}(\tau)}{0pt(\tau)}\right)\right|_{\tau=t}$
(15)
$\displaystyle=\left.\frac{\frac{\mathrm{d}\mathbf{x}(\tau)}{\mathrm{d}\tau}0pt(\tau)-\mathbf{x}(\tau)\frac{\mathrm{d}0pt(\tau)}{\mathrm{d}\tau}}{0pt(\tau)^{2}}\right|_{\tau=t}$
$\displaystyle=\left.\frac{1}{0pt(\tau)}\left(\frac{\mathrm{d}\mathbf{x}(\tau)}{\mathrm{d}\tau}-\boldsymbol{\eta}(\tau)\langle\boldsymbol{\eta}(\tau),\frac{\mathrm{d}\mathbf{x}(\tau)}{\mathrm{d}\tau}\rangle\right)\right|_{\tau=t}$
$\displaystyle=\frac{1}{0pt(\boldsymbol{\eta},t)}(\mathbf{I}-\boldsymbol{\eta}\boldsymbol{\eta}^{\top})\mathbf{v}(\boldsymbol{\eta},t)$
$\displaystyle=\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\mathbf{w}(\boldsymbol{\eta},t)$
(16)

Note that $\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\mathbf{w}$
satisfies
$\langle\boldsymbol{\eta},\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\mathbf{w}\rangle=0$
and hence, is an element of the tangent space of $\boldsymbol{\eta}$. Indeed,
the total time derivative of $\boldsymbol{\eta}(\tau)$ is the optical flow
$\boldsymbol{\Phi}(\boldsymbol{\eta},t)=\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\mathbf{w}$
perceived by the spherical camera [KvD87, HM02].

Finally, consider the total derivative of the structure flow vector
$\mathbf{w}(\tau)=\mathbf{w}(\boldsymbol{\eta}(\tau),\tau)$ taken along a
trajectory induced by the equations of motion (3), (4), (5) and (6). One has

$\displaystyle\left.\frac{\mathrm{d}\mathbf{w}(\tau)}{\mathrm{d}\tau}\right|_{\tau=t}$
$\displaystyle=\left.\frac{\mathrm{d}}{\mathrm{d}\tau}\left(-\Omega_{\times}\boldsymbol{\eta}(\tau)+\frac{1}{0pt(\tau)}(\mathbf{v}_{\mathbf{x}}(\tau)-\mathbf{v}_{c}(\tau))\right)\right|_{\tau=t}$
$\displaystyle=-\left.\frac{\mathrm{d}}{\mathrm{d}\tau}\left(\Omega_{\times}\boldsymbol{\eta}(\tau)\right)\right|_{\tau=t}+\left.\frac{\mathrm{d}\rho(\tau)}{\mathrm{d}\tau}(\mathbf{v}_{\mathbf{x}}-\mathbf{v}_{c})\right|_{\tau=t}$
$\displaystyle\quad\quad\quad\quad\quad+\rho(\tau)\left.\frac{\mathrm{d}}{\mathrm{d}\tau}(\mathbf{v}_{\mathbf{x}}(\tau)-\mathbf{v}_{c}(\tau))\right|_{\tau=t}$
$\displaystyle=-\Omega_{\times}\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\mathbf{w}(\boldsymbol{\eta},t)$
$\displaystyle\quad\quad\quad-\rho(\boldsymbol{\eta})\langle\boldsymbol{\eta},\mathbf{w}(\boldsymbol{\eta},t)\rangle(\mathbf{v}_{\mathbf{x}}(t)-\mathbf{v}_{c}(t))$
$\displaystyle\quad\quad\quad\quad-\rho(\boldsymbol{\eta})\Omega_{\times}(\mathbf{v}_{\mathbf{x}}(t)-\mathbf{v}_{c}(t))+\rho(\boldsymbol{\eta})\boldsymbol{a}_{c}$
$\displaystyle=-\Omega_{\times}\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\mathbf{w}(\boldsymbol{\eta},t)-\mathbf{w}_{s}(\boldsymbol{\eta},t)\langle\boldsymbol{\eta},\mathbf{w}(\boldsymbol{\eta},t)\rangle$
$\displaystyle\quad\quad\quad-\Omega_{\times}\mathbf{w}_{s}(\boldsymbol{\eta},t)+\rho(\boldsymbol{\eta})\boldsymbol{a}_{c}$

We define a grouped acceleration term

$\displaystyle\boldsymbol{a}_{\mathbf{w}}(\boldsymbol{\eta},t)$
$\displaystyle=\rho\boldsymbol{a}_{c}-\Omega_{\times}\mathbf{w}_{s}(\boldsymbol{\eta},t)$
$\displaystyle=\rho\boldsymbol{a}_{c}-\Omega_{\times}\left(\mathbf{w}+\Omega_{\times}\boldsymbol{\eta}\right).$

This is the exogenous linear acceleration of the structure flow field at
$(\boldsymbol{\eta},t)$ due to camera motion along with the coriolis term
$\Omega_{\times}\mathbf{w}_{s}(\boldsymbol{\eta},t)$ associated with the fact
that the linear velocities are expressed in a rotating frame. Recall that
$\mathbf{w}_{s}=\mathbf{w}-\mathbf{w}_{r}$ and
$\mathbf{w}_{r}=-\Omega_{\times}\boldsymbol{\eta}$. Rewriting in terms of the
full structure flow and $\boldsymbol{a}_{\mathbf{w}}$ one has

$\displaystyle\left.\frac{\mathrm{d}\mathbf{w}(\tau)}{\mathrm{d}\tau}\right|_{\tau=t}=$
$\displaystyle-\Omega_{\times}\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\mathbf{w}(\boldsymbol{\eta},t)+\Omega_{\times}\boldsymbol{\eta}\langle\boldsymbol{\eta},\mathbf{w}(\boldsymbol{\eta},t)\rangle$
$\displaystyle\quad\quad\quad\quad-\mathbf{w}(\boldsymbol{\eta},t)\langle\boldsymbol{\eta},\mathbf{w}(\boldsymbol{\eta},t)\rangle+\boldsymbol{a}_{\mathbf{w}}(\boldsymbol{\eta},t)$
$\displaystyle=$
$\displaystyle-\mathbf{w}(\boldsymbol{\eta},t)\langle\boldsymbol{\eta},\mathbf{w}(\boldsymbol{\eta},t)\rangle-\Omega_{\times}\mathbf{w}(\boldsymbol{\eta},t)+\boldsymbol{a}_{\mathbf{w}}(\boldsymbol{\eta},t).$
(17)

Here, the term $-\Omega_{\times}\mathbf{w}$ is to be expected since
$\mathbf{w}$ is expressed with respect to the rotating camera frame. The
driving term $-\mathbf{w}\langle\boldsymbol{\eta},\mathbf{w}\rangle$ is the
autonomous growth or decrease in $\mathbf{w}$ due to its range velocity
component and the field acceleration $\boldsymbol{a}_{\mathbf{w}}$ was defined
earlier.

Figure 2 illustrates the $xyz$ components of scene, structure and optical flow
fields for a static scene and a forward moving camera (in the $z$ direction)
with no rotation. Since the scene is static, the scene flow equals
$-\mathbf{v}_{c}$ for all pixels and in particular, the $xy$ components of
scene flow are zero. The structure flow combines the scene flow and the
inverse depth at each pixel to create a velocity field that differentiates
between close and distant objects. In contrast, the optical flow generates a
divergent velocity field product of the projection onto tangent space. The
focus of expansion of the optical flow at the image center corresponds to the
direction of motion.

Figure 2: Scene, structure and optical flow fields for a forward moving camera
in a static scene.

###  3.2 Evolution equations

This section derives the fundamental partial differential equations that model
the temporal change of image brightness, structure flow and inverse depth on
the spherical camera.

####  3.2.1 Image brightness

Let $Y(\boldsymbol{\eta},t)$ denote the image brightness in the direction
$\boldsymbol{\eta}\in S^{2}$ at time $t$ associated to point $\mathbf{x}(t)$
in the scene. That is $Y:S^{2}\times\mathbb{R}\rightarrow\mathbb{R}$ is
modeled as a scalar field on the sphere. A common assumption on differential
methods for computing either optical or scene flow is that image brightness is
constant and the surface is Lambertian. That is, the brightness value
corresponding to point $\mathbf{x}(t)$ as seen by the camera does not change
with the evolution of $\mathbf{x}(t)$ in time. This assumption corresponds to
the well known constant brightness condition on the _total derivative_ of the
image brightness field,

$\frac{\mathrm{d}Y}{\mathrm{d}t}(t)=\left.\frac{\mathrm{d}Y(\boldsymbol{\eta}(\tau))}{\mathrm{d}\tau}\right|_{\tau=t}=0.$
(18)

Expressing the total derivative of $Y(\boldsymbol{\eta},t)$ in terms of its
partial derivatives one obtains

$\frac{\partial
Y}{\partial\boldsymbol{\eta}}\frac{\mathrm{d}\boldsymbol{\eta}}{\mathrm{d}t}+\frac{\partial
Y}{\partial t}=\frac{\mathrm{d}Y}{\mathrm{d}t}$ (19)

Here $\frac{\partial Y}{\partial\boldsymbol{\eta}}\in T_{\eta}S^{2}$ denotes
the image gradient and it is an element of the tangent space of
$\boldsymbol{\eta}$. Moreover,
$\frac{\mathrm{d}\boldsymbol{\eta}}{\mathrm{d}t}=\left.\frac{\mathrm{d}\boldsymbol{\eta}(\tau)}{\mathrm{d}\tau}\right|_{\tau=t}=\boldsymbol{\Phi}(\boldsymbol{\eta},t)$
is the optical flow vector at $\boldsymbol{\eta}$ defined in Equation (16).
Expressing PDE (19) in terms of structure flow, one has

$\frac{\partial Y}{\partial t}=-\frac{\partial
Y}{\partial\boldsymbol{\eta}}\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\mathbf{w}$
(20)

Equation (20) is the classical brightness conservation equation [Bar94]
expressed in spherical camera coordinates and with respect to structure flow.

####  3.2.2 Structure flow

Analogous to the brightness conservation, we consider
$\mathbf{w}(\boldsymbol{\eta},t)$ as a continuous vector field on the sphere.
The relative change with respect to $\boldsymbol{\eta}$ and $t$ is

$\frac{\partial\mathbf{w}}{\partial\boldsymbol{\eta}}\frac{\mathrm{d}\boldsymbol{\eta}}{\mathrm{d}t}+\frac{\partial\mathbf{w}}{\partial
t}=\frac{\mathrm{d}\mathbf{w}}{\mathrm{d}t}$ (21)

where $\frac{\partial\mathbf{w}}{\partial\boldsymbol{\eta}}$ is the Jacobian
matrix of $\mathbf{w}$ at $\boldsymbol{\eta}$ and
$\frac{\mathrm{d}\mathbf{w}}{\mathrm{d}t}=\left.\frac{\mathrm{d}\mathbf{w}(\boldsymbol{\eta}(\tau),\tau)}{\mathrm{d}\tau}\right|_{\tau=t}$
is the optical flow as derived in (16). In this case, the total flow is not
conserved, that is, $\frac{\mathrm{d}\mathbf{w}}{\mathrm{d}t}\neq 0$, and one
must substitute (17) to derive the associated PDE

$\frac{\partial\mathbf{w}}{\partial
t}=-\frac{\partial\mathbf{w}}{\partial\boldsymbol{\eta}}\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\mathbf{w}-\mathbf{w}\langle\boldsymbol{\eta},\mathbf{w}\rangle-\Omega_{\times}\mathbf{w}+\boldsymbol{a}_{\mathbf{w}}$
(22)

Equation (22) describes a non-linear transport process where the structure
flow is propagated on the sphere at a velocity equal to the induced optical
flow $\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\mathbf{w}$.
Additionally, the source terms inject or remove energy from the structure flow
field according to the kinematics of the robot and the scene relative to it.

####  3.2.3 Depth and inverse depth

Similar partial differential equations can be derived for the conservation of
the depth and inverse depth fields, $0pt(\boldsymbol{\eta},t)$ and
$\rho(\boldsymbol{\eta},t)$, respectively. Considering the total derivatives
of depth and inverse depth developed in Equations (13) and (14), we obtain:

$\displaystyle\frac{\partial 0pt}{\partial t}$ $\displaystyle=-\frac{\partial
0pt}{\partial\boldsymbol{\eta}}\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\mathbf{w}+0pt\langle\boldsymbol{\eta},\mathbf{w}\rangle$
(23) $\displaystyle\frac{\partial\rho}{\partial t}$
$\displaystyle=-\frac{\partial\rho}{\partial\boldsymbol{\eta}}\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\mathbf{w}-\rho\langle\boldsymbol{\eta},\mathbf{w}\rangle$
(24)

The conservation of inverse depth has previously been reported in the
literature [BR09, ZAR12].

In practice, Equation (24) is preferred since $\rho$ is a dimensionless ratio
(8) that has an intrinsic interpretation as visual measure. Note also that
(24) has the same structure as the first two terms of (22) which leads to
desirable computational properties of the combined system of equations
considered in Section 4.

##  4 Filter Architecture

This section introduces a filtering approach for the computation of structure
flow in real time using brightness and depth measurements. We use a pyramidal
structure of filter blocks to support large pixel displacements per frame. Let
$H$ denote the number of levels in the pyramid structure indexed as
$h=1,\dots,H$, where $h=1$ denotes the original resolution level.

The filter state at discrete time index $k$ is denoted by the set

$X^{k}=\left\\{\begin{Bmatrix}\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{1}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{k}}}\\\
\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{1}}]{\rho}{{}^{\scalebox{0.7}{k}}}\end{Bmatrix},\begin{Bmatrix}\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{2}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{k}}}\\\
\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{2}}]{\rho}{{}^{\scalebox{0.7}{k}}}\end{Bmatrix},\dots,\begin{Bmatrix}\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k}}}\\\
\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\rho}{{}^{\scalebox{0.7}{k}}}\end{Bmatrix}\right\\}$
(25)

Each level $h$ contains an estimate of the inverse depth field
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{k}}}$
using data of the corresponding pyramid level. State
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k}}}$
at top level $H$ represents a coarse estimate of the structure flow based on
low resolution down-sampled data. The flow at each level, $h$, of the pyramid
is denoted
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k}}}$,
however, the flow itself is not used as the dynamic state of the filter for
the lower levels of the pyramid. Instead, the lower level filter states
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{k}}}$
define the increment to the flow
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h+1}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k}}}$
given higher resolution data at level $h$ and the the structure flow is
reconstructed by

$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k}}}=\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h+1:h}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k}}}+\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{k}}}$
(26)

applied recursively for $h=H-1,\dots,1$. Here
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h+1:h}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k}}}$
is the flow at level $h+1$ up-sampled to level $h$.

(a) Top level filter loop. The structure flow
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k}}}$
is used in the prediction stage $[P]$ to compute predictions
$\\{\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k+}}},\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\rho}{{}^{\scalebox{0.7}{k+}}}\\}$.
New raw image and inverse depth measurements are processed at $[I_{Y}]$ and
$[I_{\rho}]$ blocks to extract model parameters. These parameters are used in
the update stage $[U]$ to compute new state estimate
$\\{\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k+1}}},\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\rho}{{}^{\scalebox{0.7}{k+1}}}\\}$.

(b) Lower levels filter loop. The reconstruction block $[R]$ reconstructs the
flow at this level using the flow from level $h+1$ and state
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{k+1}}}$.
The reconstructed flow is used in the prediction stage $[P_{[]}]$ to create
predictions
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{Y}{{}^{\scalebox{0.7}{k+}}},\\{\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{k+}}},\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{k+}}}\\}$
and
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k+}}}$.
These predictions are combined with image and inverse depth model parameters
at the update stage $[\Delta U]$ to create a new estimate
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{k+1}}}$
and
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{k+1}}}$.

Figure 3: Pyramidal filter architecture.

Figure 3(a) illustrates the filter structure for the topmost level $H$. Blocks
$[I_{Y}]$ and $[I_{\rho}]$ extract linear model parameters from image and
inverse depth raw measurement data. Details of these processes are provided in
Sections 4.2 and 4.3. In the prediction block $[P]$ (Section 4.4), the current
estimate
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k}}}$
is used to create a prediction
$\\{\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k+}}},\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\rho}{{}^{\scalebox{0.7}{k+}}}\\}$
for the flow and inverse depth state at time $k+1$. We use notation $k+$ to
refer to the state before new measurements at time $k+1$ are incorporated into
it. The update block $[U]$ (Section 4.5) takes the predicted state
$\\{\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k+}}},\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\rho}{{}^{\scalebox{0.7}{k+}}}\\}$
and brightness and inverse depth parameters at $k+1$ to create a new state
estimate
$\\{\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k+1}}},\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\rho}{{}^{\scalebox{0.7}{k+1}}}\\}$.

Figure 3(b) illustrates the filter structure for levels $h=1,\dots,H-1$. The
structure flow output is computed at the reconstruction block $[R]$ using the
output of the update block $[\Delta U]$ and the flow from level $h+1$
(Equation (26)). In the prediction block $[P_{[]}]$, the reconstructed flow is
used to create predictions of flow
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k+}}}$,
image brightness
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{Y}{{}^{\scalebox{0.7}{k+}}}$
and state
$\\{\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{k+}}},\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{k+}}}\\}$.
These predictions are plugged into the update block $[\Delta U]$ and combined
with new measurement parameters from $[I_{Y}]$ and $[I_{\rho}]$ blocks to
compute a new estimate
$\\{\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{k+1}}},\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{k+1}}}\\}$.

###  4.1 Spherical image data structure

Numerical implementation of the different stages of the algorithm are
performed directly on the Spherepix data structure [AM16b], a data structure
based on spherical image geometry. This allows us to preserve the original
geometry used to derive the conservation equations in Section 3.2, ensures
robustness and consistency of the algorithm over wide field of view, and
provides data in an ideal format for incorporation into robotic control
algorithms.

Figure 4: Spherepix grid. Spherical coordinate $\boldsymbol{\eta}_{ij}$ and
its associated tangent plane $T_{\eta_{ij}}S^{2}$. Vector $\boldsymbol{\mu}\in
T_{\eta_{ij}}S^{2}$ can be represented as a 2-vector $\boldsymbol{\beta}$ in
the (row, column) direction.

Figure 4 illustrates the relevant properties of a Spherepix image. At a given
pixel $\boldsymbol{\eta}_{ij}\in S^{2}$, neighboring points
$\boldsymbol{\eta}$ are projected onto tangent space coordinates
$\boldsymbol{\mu}$ using the orthographic projection

$\boldsymbol{\mu}=\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}_{ij}$}}\boldsymbol{\eta}.$
(27)

As points $\boldsymbol{\mu}\in T_{\eta_{ij}}S^{2}$ only have two degrees of
freedom, it is convenient to express them as 2D plane coordinates. Let
$\boldsymbol{\beta}\in\mathbb{R}^{2}$ be the 2D coordinate representation of
$\boldsymbol{\mu}$. One has

$\boldsymbol{\beta}=\mathbf{B}_{\scalebox{0.6}{$\boldsymbol{\eta}_{ij}$}}^{\top}\boldsymbol{\mu}$
(28)

where $\mathbf{B}_{\scalebox{0.6}{$\boldsymbol{\eta}_{ij}$}}$ is an
orthonormal basis matrix at $\boldsymbol{\eta}_{ij}$. The row entries of
$\mathbf{B}_{\scalebox{0.6}{$\boldsymbol{\eta}_{ij}$}}$ correspond to the
normalized projections of pixels $\boldsymbol{\eta}_{i,j+1}$ and
$\boldsymbol{\eta}_{i+1,j}$ on the tangent space of $\boldsymbol{\eta}_{ij}$.
The transformation from 2D to 3D tangent space coordinates is

$\boldsymbol{\mu}=\mathbf{B}_{\scalebox{0.6}{$\boldsymbol{\eta}_{ij}$}}\boldsymbol{\beta}$
(29)

since $\mathbf{B}_{\scalebox{0.6}{$\boldsymbol{\eta}_{ij}$}}$ is orthonormal
and $\boldsymbol{\mu}$ lies in the tangent space. Finally, a scalar
$\Delta\mu_{ij}=\left\|\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}_{ij}$}}\boldsymbol{\eta}_{i,j+1}\right\|$
measures the pixel separation between two neighbors in tangent space.

###  4.2 Brightness parameters extraction

Raw image measurements $Y$ from the camera are used to fit a linear brightness
model for each pixel $p=(i,j)$. The model parameters correspond to the image
gradient $\partial_{\boldsymbol{\eta}}\hat{Y}_{p}\in T_{\eta_{p}}S^{2}$ and
the constant term $\hat{Y}_{p}\in\mathbb{R}$ at $p$.

Computation of these parameters is performed in two stages. First, the cost
function (30) is minimized for parameters
$\\{\partial_{\boldsymbol{\beta}}\hat{Y}_{p},\hat{Y}_{p}\\}$ expressed in
local 2D plane coordinates.

$\varepsilon_{p}=\sum_{q\in\Lambda_{p}}G(p,q)\|Y_{q}-\partial_{\boldsymbol{\beta}}\hat{Y}_{p}(p-q)-\hat{Y}_{p}\|^{2},$
(30)

where $\Lambda_{p}$ is a $5\times 5$ support window of pixels centered on $p$,
and $G(p,q)$ is a Gaussian weight mask $G=g^{\top}g$ with $g=[1,4,6,4,1]/16$.
Thanks to the spatial symmetry of the support window and weight function, it
is possible to efficiently solve for
$\\{\partial_{\boldsymbol{\beta}}\hat{Y}_{p},\hat{Y}_{p}\\}$ as a series of 1D
convolutions in the row and column axis. The two-dimensional image gradient
$\partial_{\boldsymbol{\beta}}\hat{Y}_{p}$ is then expressed as a 3D tangent
vector $\partial_{\boldsymbol{\eta}}\hat{Y}_{p}$ following Equation (29)

$\partial_{\boldsymbol{\eta}}\hat{Y}_{p}=\Delta\mu_{p}\mathbf{B}_{\scalebox{0.6}{$p$}}^{\top}\partial_{\boldsymbol{\beta}}\hat{Y}_{p}.$
(31)

###  4.3 Inverse depth parameters extraction

In this paper we assume that depth measurements are available and are inverted
for the data fitting process. We fit a linear model comprised of parameters
$\partial_{\boldsymbol{\eta}}\hat{\rho}_{p}\in T_{\eta_{p}}S^{2}$ for the
inverse depth gradient vector and $\hat{\rho}_{p}\in\mathbb{R}$ for the
constant term. In particular for the computation of
$\partial_{\boldsymbol{\eta}}\hat{\rho}_{p}$, one must properly handle object
occlusions to remove undesired large gradient vectors. At each pixel location,
we compute the gradient vector on both sides of a discontinuity and select the
appropriate one. That is, the two-dimensional gradient vector
$\partial_{\boldsymbol{\beta}}\hat{\rho}_{p}=(\partial_{\boldsymbol{\beta}_{1}}\hat{\rho}_{p},\partial_{\boldsymbol{\beta}_{2}}\hat{\rho}_{p})^{\top}$
is computed as follows

$\displaystyle\partial_{\boldsymbol{\beta}_{1}}\hat{\rho}_{p}$
$\displaystyle=\left\\{\begin{matrix}\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{1}^{+}$}}\rho_{p}&\text{if}~{}~{}|\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{1}^{+}$}}\rho_{p}|\leq|\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{1}^{-}$}}\rho_{p}|\\\
\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{1}^{-}$}}\rho_{p}&\text{otherwise}\end{matrix}\right.$
(32) $\displaystyle\partial_{\boldsymbol{\beta}_{2}}\hat{\rho}_{p}$
$\displaystyle=\left\\{\begin{matrix}\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{2}^{+}$}}\rho_{p}&\text{if}~{}~{}|\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{2}^{+}$}}\rho_{p}|\leq|\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{2}^{-}$}}\rho_{p}|\\\
\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{2}^{-}$}}\rho_{p}&\text{otherwise}\end{matrix}\right.$
(33)

where $\delta_{\scalebox{0.5}{${}^{+}$}}$, $\delta_{\scalebox{0.5}{${}^{-}$}}$
denote forward and backward difference operators in the
$\boldsymbol{\beta}_{1}$ (column) and $\boldsymbol{\beta}_{2}$ (row) axis.
Table 1 defines these as well as the central difference operator.

| $\delta_{\boldsymbol{\beta}_{1}}u_{ij}$ | $\delta_{\boldsymbol{\beta}_{2}}u_{ij}$
---|---|---
backward | $\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{1}^{-}$}}=u_{j}-u_{j-1}$ | $\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{2}^{-}$}}=u_{i}-u_{i-1}$
central | $\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{1}^{0}$}}=u_{j+1}-u_{j-1}$ | $\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{2}^{0}$}}=u_{i+1}-u_{i-1}$
forward | $\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{1}^{+}$}}=u_{j+1}-u_{j}$ | $\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{2}^{+}$}}=u_{i+1}-u_{i}$
Table 1: Difference operators in beta coordinates.

The 3D vector representation $\partial_{\boldsymbol{\eta}}\hat{\rho}_{p}\in
T_{\eta_{p}}S^{2}$ is

$\partial_{\boldsymbol{\eta}}\hat{\rho}_{p}=\Delta\mu_{p}\mathbf{B}_{\scalebox{0.6}{$p$}}^{\top}\partial_{\boldsymbol{\eta}}\hat{\rho}_{p}$
(34)

The constant parameter term is given by the raw inverse depth data; that is,
$\hat{\rho}_{p}=\rho_{p}$.

###  4.4 State prediction ($k\rightarrow k+$)

The prediction stage of the filter uses the current state estimate to create
predictions of the structure flow forward in time. In the present paper, we
choose to ignore the source terms derived from inertial measurements in the
PDE (36). That is, we assume

$\Omega_{\times}\mathbf{w}-\boldsymbol{a}_{\mathbf{w}}\approx 0.$ (35)

In practice, for high frame rates and the scenarios that we consider, these
terms have sub-pixel magnitudes and are at least an order of magnitude less
significant than the other terms in the propagation algorithm. In particular,
the contribution of this term is proportional to the magnitude of $\Omega$ and
$\boldsymbol{a}_{c}$ and inversely proportional to frame rate. Thus, for high
frame rates and moderate motion these terms become insignificant. Figure 5
plots the $xyz$ components of Equation (35) for a structure flow field
computed from realistic velocities and acceleration values to provide an
indication of the validity of this assumption.

For the top level $H$, the prediction block uses the current coarse estimate
of structure flow and the inverse depth field as initial conditions
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}(0):=\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k}}}$
and
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\rho}{{}^{\scalebox{0.7}{}}}(0):=\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\rho}{{}^{\scalebox{0.7}{k}}}$
for the PDEs (22) and (24) with assumption (35). We numerically solve these
PDEs to create the state predictions
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k+}}}:=\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}(1)$
and
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\rho}{{}^{\scalebox{0.7}{k+}}}:=\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\rho}{{}^{\scalebox{0.7}{}}}(1)$
for the next frame $k+1$. The explicit PDEs that we solve at level $H$ are

$\displaystyle\frac{\partial\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}}{\partial
t}$
$\displaystyle=-\frac{\partial\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}}{\partial\boldsymbol{\eta}}\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}-\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}\langle\boldsymbol{\eta},\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}\rangle$
(36)
$\displaystyle\frac{\partial\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\rho}{{}^{\scalebox{0.7}{}}}}{\partial
t}$
$\displaystyle=-\frac{\partial\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\rho}{{}^{\scalebox{0.7}{}}}}{\partial\boldsymbol{\eta}}\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}-\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\rho}{{}^{\scalebox{0.7}{}}}\langle\boldsymbol{\eta},\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}\rangle$
(37)

For lower levels $h=1,\dots,H-1$, the propagation scheme is defined by a PDE
system modeling the transport of the structure flow, image brightness and
state
$\\{\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{}}},\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{}}}\\}$
by the reconstructed structure flow
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}$
at each level.

$\displaystyle\frac{\partial\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}}{\partial
t}$
$\displaystyle=-\frac{\partial\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}}{\partial\boldsymbol{\eta}}\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}-\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}\langle\boldsymbol{\eta},\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}\rangle$
(38)
$\displaystyle\frac{\partial\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{}}}}{\partial
t}$
$\displaystyle=-\frac{\partial\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{}}}}{\partial\boldsymbol{\eta}}\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}-\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{}}}\langle\boldsymbol{\eta},\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}\rangle$
(39)
$\displaystyle\frac{\partial\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{}}}}{\partial
t}$
$\displaystyle=-\frac{\partial\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{}}}}{\partial\boldsymbol{\eta}}\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}-\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{}}}\langle\boldsymbol{\eta},\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}\rangle$
(40)
$\displaystyle\frac{\partial\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{Y}{{}^{\scalebox{0.7}{}}}}{\partial
t}$
$\displaystyle=-\frac{\partial\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{Y}{{}^{\scalebox{0.7}{}}}}{\partial\boldsymbol{\eta}}\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}$
(41)

Initial conditions are set to
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\mathbf{w}}{{}^{\scalebox{0.7}{}}}(0):=\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k}}}$,
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{}}}(0):=\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{k}}}$,
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{}}}(0):=\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{k}}}$
and
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{Y}{{}^{\scalebox{0.7}{}}}(0):=\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{Y}{{}^{\scalebox{0.7}{k}}}$
and the system is run for one time step unit. Details on the numerical
solution to these equations are provided in Section 4.6.

Figure 5: Contribution of the source terms in Equation (35). The simulated
camera runs at 300 Hz and moves with a linear velocity of
$5\nicefrac{{m}}{{s}}$, acceleration of $1\nicefrac{{m}}{{s^{2}}}$ both in the
$z$ direction and angular velocity of $180\nicefrac{{\deg}}{{s}}$ in $y$ axis.

###  4.5 State update ($k+\rightarrow k+1$)

In the update stage of the filter, predictions made for time $k+1$ are
corrected using new image and inverse depth measurements. There are two update
schemes used; one for top level $H$ to create the updated state
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k+1}}}$
and one applied to the remaining lower levels to update
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{K+1}}}$.
The update of the inverse depth state
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{k+1}}}$
follows the same scheme for all pyramid levels.

For level $H$, the update stage creates a new estimate
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{K+1}}}$
using low resolution data available at level $H$. Cost function (42) uses low
resolution image and inverse depth measurements, as well as the predicted flow
from the propagation stage to create a new structure flow estimate
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k+1}}}$.
The cost function

$^{H}\varepsilon_{\mathbf{w}}=\gamma_{1}\|^{H}E_{Y}\|^{2}+\gamma_{2}\|^{H}E_{\rho}\|^{2}+\gamma_{3}\|^{H}E_{t}\|^{2}$
(42)

consists of three data terms: ${}^{H}E_{Y}$ (image brightness conservation),
${}^{H}E_{\rho}$ (depth measurements) and ${}^{H}E_{t}$ (regularisation based
on previous estimate). Scalar gains $\gamma_{1}$, $\gamma_{2}$, $\gamma_{3}$
control the relative weight of each term.

The image data term $E_{Y}$ is based on the image conservation PDE (20). It is
defined as

$^{H}E_{Y}=\partial_{\boldsymbol{\eta}}\hat{Y}^{k+1}\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\mathbf{w}^{k+1}+\Delta\mu^{2}(\hat{Y}^{k+1}-\hat{Y}^{k})$
(43)

where $\partial_{\boldsymbol{\eta}}\hat{Y}^{k+1}$, $\hat{Y}^{k+1}$ and
$\hat{Y}^{k}$ are the image gradient and constant terms computed at the image
model block. The temporal difference term is multiplied by $\Delta\mu^{2}$ to
compensate for the transformation of gradient vector
$\partial_{\boldsymbol{\eta}}\hat{Y}^{k+1}$ to tangent space coordinates. This
cost term only constrains one degree of freedom of $\mathbf{w}$, corresponding
to the tangential component of structure flow in the direction of the image
gradient vector $\partial_{\boldsymbol{\eta}}\hat{Y}$. This is the well known
aperture problem in optical flow estimation.

The inverse depth data term $E_{\rho}$ uses the inverse depth conservation PDE
(24) to recover the structure flow from state $\rho^{k}$ and measurements and
$\hat{\rho}^{k+1}$. The data term is

${}^{H}E_{\rho}=~{}$
$\displaystyle\partial_{\boldsymbol{\eta}}\hat{\rho}^{k+1}\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\mathbf{w}^{k+1}$
$\displaystyle+\Delta\mu^{2}\left((\hat{\rho}^{k+1}-\rho^{k})+\hat{\rho}^{k+1}\langle\boldsymbol{\eta},\mathbf{w}^{k+1}\rangle\right)$
(44)

Analogously to the image cost term, the scalar component of Equation (44) is
multiplied by $\Delta\mu^{2}$ to compensate for the transformation of vector
quantities to 3D tangent space coordinates. As PDE (24) includes the effect of
$\mathbf{w}$ along normal direction $\boldsymbol{\eta}$, it is possible to
recover both tangent and normal components of the flow field. For the tangent
component, it still suffers from a form of aperture problem where only flow in
the direction of the inverse depth gradient can be recovered [SJB02].

Finally, the temporal smoothing or regularisation term, $E_{t}$, includes the
predicted flow
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k+}}}$
as prior solution to the cost function

${}^{H}E_{t}=\mathbf{w}^{k+1}+\mathbf{w}^{k+}$ (45)

This regularization term guarantees there is a solution for the structure flow
at each pixel. In regions where it is not possible to recover flow from the
data, the predicted flow
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k+}}}$
will act as best estimate for time $k+1$.

Cost function (46) describes a linear least squares problem with respect to
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k+1}}}$.
Solving for the minimum of (42) can be expressed as a $3\times 3$ linear
system

$A\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k+1}}}=b$
(46)

After solving Equation (46), we apply an average smoothing filter of size
$5\times 5$ to the resulting flow field
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{H}}]{\mathbf{w}}{{}^{\scalebox{0.7}{k+1}}}$
in order to accelerate the diffusion of flow estimates to textureless regions
of the image.

For levels $h=0,\dots,H-1$, the update stage computes new state estimates
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{k+1}}}$
provided with the predicted state and new measurement data. Cost function (47)
is defined as

$^{h}\varepsilon_{\mathbf{w}}=\gamma_{1}\|^{h}E_{Y}\|^{2}+\gamma_{2}\|^{h}E_{\rho}\|^{2}+\gamma_{3}\|^{h}E_{t}\|^{2}$
(47)

where the data terms for image, inverse depth and temporal smoothing defined
are

${}^{h}E_{Y}=~{}$
$\displaystyle\partial_{\boldsymbol{\eta}}\hat{Y}^{k+1}\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{k+1}}}+\Delta\mu^{2}(\hat{Y}^{k+1}-\hat{Y}^{k+})$
(48) ${}^{h}E_{\rho}=~{}$
$\displaystyle\partial_{\boldsymbol{\eta}}\hat{\rho}^{k+1}\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{k+1}}}$
$\displaystyle+\Delta\mu^{2}[(\hat{\rho}^{k+1}-\rho^{k+})+\hat{\rho}^{k+1}\langle\boldsymbol{\eta},\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{k+1}}}\rangle]$
(49) ${}^{h}E_{t}=~{}$
$\displaystyle\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{k+1}}}-\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{k+}}}$
(50)

Equation (47) describes a linear least squares problem with respect to unknown
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{k+1}}}$,
and can be solved following the same procedure used for top level. Once
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\scalebox{0.6}{$\Delta$}\mathbf{w}}{{}^{\scalebox{0.7}{k+1}}}$
has been computed, the structure flow at level $h$ is reconstructed using
Equation (26). Once completed, the resulting flow is cascaded to the next
level below until level $h=1$ is reached and the flow at the original
resolution is computed.

The formulation of the inverse depth state update is equal for all pyramid
levels. The proposed cost function

$^{h}\varepsilon_{\rho}=\gamma_{4}\|\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{k+1}}}-\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\hat{\rho}}{{}^{\scalebox{0.7}{k+1}}}\|^{2}+\gamma_{5}\|\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{k+1}}}-\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{k+}}}\|^{2}$
(51)

considers both current measurements and the predicted inverse depth from
previous frame for estimating new state
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{k+1}}}$.
Here $\gamma_{4}$ and $\gamma_{5}$ are the relative weights between
measurements and prediction. The least squares solution of
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{k+1}}}$
is

$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{k+1}}}=\frac{\gamma_{4}\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\hat{\rho}}{{}^{\scalebox{0.7}{k+1}}}+\gamma_{5}\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{k+}}}}{\gamma_{4}+\gamma_{5}}$
(52)

In places where no measurements are available, $\gamma_{4}=0$ and the updated
inverse depth is
$\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{k+}}}=\tensor*[_{\scalebox{0.6}{}}^{\scalebox{0.6}{h}}]{\rho}{{}^{\scalebox{0.7}{k+1}}}$.

###  4.6 Numerical prediction

Considering the real-time nature of the structure flow filtering algorithm, we
desire a numerical scheme that is fast and robust to noisy input for solving
the propagation equations in Section 4.4. We propose a numerical scheme based
on upwind finite differences [Tho95] following a similar approach to our
previous work on optical flow algorithms [AM16a].

Consider Equation (53) as a discrete version of prediction Equations (36) and
(38).

$\frac{\mathbf{w}^{n+1}-\mathbf{w}^{n}}{\Delta
t}=-\frac{\partial\mathbf{w}^{n}}{\partial\boldsymbol{\eta}}\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\mathbf{w}^{n}-\mathbf{w}^{n}\langle\boldsymbol{\eta},\mathbf{w}^{n}\rangle$
(53)

where $\mathbf{w}^{n+1}$ is the unknown flow at time index $n+1$. For
convenience in the notation, we omit pyramid level index $h$ and pixel
coordinate $(i,j)$, as the scheme runs identically at every level and for each
pixel. Index $n$ refers to an internal time index for the numerical scheme
running for $n=0,\dots,N-1$ with $\Delta t=\nicefrac{{1}}{{N}}$ and $N$ being
a parameter describing the number of iterations. Initial conditions of the
problem are set to $\mathbf{w}^{0}:=\mathbf{w}^{k}$ and the output of the
scheme is $\mathbf{w}^{N-1}=:\mathbf{w}^{k+}$.

Equation (53) is a non-linear PDE on $\mathbf{w}$ for which direct application
of finite difference methods are not reliable [Tho99, p. 140]. To break this
non-linearity we use the idea of _dominant optical flow_. Optical flow in 2D
tangent plane coordinates and pixel units is computed as

$\Phi_{ij}:=\begin{pmatrix}u_{ij}\\\
v_{ij}\end{pmatrix}=\frac{1}{\Delta\mu}\mathbf{B}_{\scalebox{0.6}{$\boldsymbol{\eta}_{ij}$}}\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}_{ij}$}}\mathbf{w}_{ij}$
(54)

The dominant optical flow
$\hat{\Phi}_{ij}=(\hat{u}^{n}_{ij},\hat{v}^{n}_{ij})$ is calculated as the
largest flow component in the row and column axis

$\displaystyle\hat{u}_{ij}$
$\displaystyle=\left\\{\begin{matrix}u_{i,j-1}&\text{if}~{}~{}\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{1}^{0}$}}|u_{ij}|>0\\\
u_{i,j+1}&\text{otherwise}\end{matrix}\right.$ (55)
$\displaystyle\hat{v}_{ij}$
$\displaystyle=\left\\{\begin{matrix}v_{i-1,j}&\text{if}~{}~{}\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{2}^{0}$}}|v_{ij}|>0\\\
v_{i+1,j}&\text{otherwise}\end{matrix}\right.$ (56)

with $\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{1}^{0}$}}$ and
$\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{2}^{0}$}}$ the central difference
operators defined in Table 1.

Functions $\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{1}$}}(\cdot;\hat{u})$
and $\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{2}$}}(\cdot;\hat{v})$ define
the upwind finite difference operators in column and row axis. They are
computed as

$\displaystyle\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{1}$}}(f;\hat{u})$
$\displaystyle=\left\\{\begin{matrix}\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{1}^{-}$}}f&\text{if}~{}~{}\hat{u}>0\\\
\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{1}^{+}$}}f&\text{otherwise}\end{matrix}\right.$
(57)
$\displaystyle\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{2}$}}(f;\hat{v})$
$\displaystyle=\left\\{\begin{matrix}\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{2}^{-}$}}f&\text{if}~{}~{}\hat{v}>0\\\
\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{2}^{+}$}}f&\text{otherwise}\end{matrix}\right.$
(58)

These operators evaluate the direction of the dominant flow in the column and
row axis and select the adequate upwind difference operator. By using the
dominant optic flow the remaining PDE evolution is linear and has a stable
solution.

The iterations of the numerical scheme are as follows. First, the dominant
optical flow $\hat{u}^{n}$ is calculated from $\mathbf{w}^{n}$. Then, the
propagation takes place in the $\boldsymbol{\beta}_{1}$ (column) axis for all
points in the spherical grid as follows

$\displaystyle\mathbf{w}^{*}$ $\displaystyle=\mathbf{w}^{n}-\Delta
t\left[\hat{u}^{n}\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{1}$}}(\mathbf{w}^{n};\hat{u}^{n})+\mathbf{w}^{n}\langle\boldsymbol{\eta},\mathbf{w}^{n}\rangle\right]$
(59) $\displaystyle\scalebox{0.6}{$\Delta$}\mathbf{w}^{*}$
$\displaystyle=\scalebox{0.6}{$\Delta$}\mathbf{w}^{n}-\Delta
t\left[\hat{u}^{n}\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{1}$}}(\scalebox{0.6}{$\Delta$}\mathbf{w}^{n};\hat{u}^{n})+\scalebox{0.6}{$\Delta$}\mathbf{w}^{n}\langle\boldsymbol{\eta},\mathbf{w}^{n}\rangle\right]$
(60) $\displaystyle\rho^{*}$ $\displaystyle=\rho^{n}-\Delta
t\left[\hat{u}^{n}\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{1}$}}(\rho^{n};\hat{u}^{n})+\rho^{n}\langle\boldsymbol{\eta},\mathbf{w}^{n}\rangle\right]$
(61) $\displaystyle Y^{*}$ $\displaystyle=Y^{n}-\Delta
t\hat{u}^{n}\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{1}$}}(Y^{n};\hat{u}^{n})$
(62)

producing the propagated fields $\mathbf{w}^{*}$,
$\scalebox{0.6}{$\Delta$}\mathbf{w}^{*}$, $\rho^{*}$ and $Y^{*}$. Next, the
dominant flow $\hat{v}^{*}$ is calculated from $\mathbf{w}^{*}$ and the
propagation takes place in the $\boldsymbol{\beta}_{2}$ (row) axis

$\displaystyle\mathbf{w}^{n+1}$ $\displaystyle=\mathbf{w}^{*}-\Delta
t\left[\hat{v}^{*}\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{2}$}}(\mathbf{w}^{*};\hat{v}^{*})+\mathbf{w}^{*}\langle\boldsymbol{\eta},\mathbf{w}^{*}\rangle\right]$
(63) $\displaystyle\scalebox{0.6}{$\Delta$}\mathbf{w}^{n+1}$
$\displaystyle=\scalebox{0.6}{$\Delta$}\mathbf{w}^{*}-\Delta
t\left[\hat{v}^{*}\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{2}$}}(\scalebox{0.6}{$\Delta$}\mathbf{w}^{*};\hat{v}^{*})+\scalebox{0.6}{$\Delta$}\mathbf{w}^{*}\langle\boldsymbol{\eta},\mathbf{w}^{*}\rangle\right]$
(64) $\displaystyle\rho^{n+1}$ $\displaystyle=\rho^{*}-\Delta
t\left[\hat{v}^{*}\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{1}$}}(\rho^{*};\hat{u}^{*})+\rho^{*}\langle\boldsymbol{\eta},\mathbf{w}^{*}\rangle\right]$
(65) $\displaystyle Y^{n+1}$ $\displaystyle=Y^{*}-\Delta
t\hat{v}^{*}\delta_{\scalebox{0.5}{$\boldsymbol{\beta}_{2}$}}(Y^{*};\hat{v}^{*})$
(66)

The numerical scheme is stable if inequality in Equation (67) is satisfied for
all points in the grid for all $n$.

$\Delta
t\max_{n,i,j}\left\\{|\hat{u}^{n}_{ij}|,|\hat{v}^{n}_{ij}|\right\\}\leq 1$
(67)

Given this stability condition where $\Delta t=\nicefrac{{1}}{{N}}$, one needs
to choose an suitable number of iterations $N$ of the numerical scheme
depending on the maximum flow expected in the application. In particular, as
the camera frame rate increases, the magnitude of the estimated flow decreases
and hence, less numerical iterations are required, leading to excellent
numerical properties for real-time high-speed vision applications.

##  5 Experimental Validation

###  5.1 GPU implementation

We developed a GPU implementation of the filtering algorithm for its
deployment in real-time applications. Computations within each stage of the
algorithm are pixel independent and can easily be parallelized on GPU
hardware. The different stages of the algorithm are programmed and connected
using a pipeline pattern. We run the implementation on a Desktop computer with
an Intel i7-4790K CPU and a Nvidia GTX-780 graphics card with 2304 CUDA cores
and 3GB of memory.

###  5.2 Ground truth evaluation

To validate the accuracy of the proposed algorithm, we require video sequences
recorded at high frame rates that sample the dynamics of the environment
smoothly. We used the Blender source files of the Urban Canyon Dataset from
Zhang et. al. [ZRFS16] to render a photo-realistic city environment in which a
300Hz camera is attached to the front part of a moving vehicle. For the
evaluation, we consider a subset of the trajectory (Figure 6(b)) composed of
10000 frames. Images are rendered using a 360 degrees panoramic camera and
then mapped to a front facing spherepix grid of $512\times 512$ pixels.

The structure flow ground truth $\mathbf{w}_{\text{gt}}$ is computed for each
pixel from the rendered depth map and the camera linear and angular velocity
using Equation (10) with $\mathbf{v}_{\mathbf{x}}\equiv 0$.

(a) Panoramic rendering.

(b) Vehicle trajectory.

Figure 6: Urban Canyon dataset.

We use _Root Mean Square_ error (RMSE) and _Absolute Angular Error_ (AAE) to
measure the accuracy of our algorithm relative to ground truth. The RMSE is
calculated by

$\text{RMSE}=\left\|\frac{\mathbf{w}_{\text{gt}}-\mathbf{w}}{\Delta\mu}\right\|$
(68)

and the AAE is computed as

$\text{AAE}=\arccos\left(\frac{1+\langle\mathbf{w}_{\text{gt}},\mathbf{w}\rangle}{\sqrt{1+\mathbf{w}_{\text{gt}}}\sqrt{1+\mathbf{w}}}\right)$
(69)

For displaying purposes, we split the 3D structure flow $\mathbf{w}$ into
tangent and normal components. Tangent flow
$\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\mathbf{w}(\boldsymbol{\eta},t)$,
that is, the structure flow projected onto the tangent space, is transformed
to 2D plane coordinates and is expressed in pixel units

$\mathbf{w}_{\perp}=\frac{\mathbf{B}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}^{\top}\mathbf{P}_{\scalebox{0.6}{$\boldsymbol{\eta}$}}\mathbf{w}}{\Delta\mu}.$
(70)

It is possible to use the standard optical flow color wheel encoding [BSL+11]
to represent the tangent flow vector field as a color image. The normal flow,
Equation (71), is the component of the structure flow along direction
$\boldsymbol{\eta}$

$\mathbf{w}_{\|}=\frac{\langle\boldsymbol{\eta},\mathbf{w}\rangle}{\Delta\mu}$
(71)

It is also measured in pixel units and plotted as a scalar field.

Figure 7 shows a selection of the first 300 images corresponding to 1 second
of real-life time. The initial condition of the structure flow is set to zero.
The flow is first identified at regions with brightness or depth
discontinuities and it is diffused to textureless regions. Over time, a dense
flow field is estimated and maintained using new measurements. Convergence of
the algorithm can be visualized through the RMSE field plot (bottom row).
Regions with higher error are localized at occlusion boundaries, where the
algorithm needs to re-identify flow.

Figure 7: Results on the Urban Canyon dataset. The sequence plots the
estimated structure flow and ground truth for the first second (300 frames) of
video.

Figure 8(a) plots the average RMSE error for the evaluated sequence. After
approximately 150 frames, the filter algorithm has converged to a dense flow
estimate. After this, the RMSE remains stable around 0.3 pixels and increases
in parts of the sequence with high camera rotation, as the algorithm needs to
adapt to the changes in velocity. Figure 8(b) displays the average AAE, which
stays approximately constant at 20 degrees.

(a) Average Root Mean Square Error.

(b) Average Absolute Angular Error.

Figure 8: Error metrics for the Urban Canyon dataset.

In terms of runtime performance, Tables 2 and 3 shows the estimated runtime
and frame rates of the algorithm configured with one and two pyramid levels
respectively. For each configuration, we vary the maximum flow allowed in the
algorithm, which affects the number of numerical iterations required in the
propagation stage. Runtime is considered as the sum of time required to map
image and depth measurements onto the spherepix image plus the runtime of the
flow filter itself. We exclude memory transfer operations between CPU and GPU
memory spaces, as they can be overlaid with the execution of GPU kernel code.

Max flow | Runtime | Frequency | Throughput
---|---|---|---
(pixels) | (ms) | (Hz) | (Mpix/s)
1 | 0.929 $\nicefrac{{+}}{{-}}$ 0.3 | 1133.9 $\nicefrac{{+}}{{-}}$ 176.9 | 283.483 $\nicefrac{{+}}{{-}}$ 44.2
2 | 1.029 $\nicefrac{{+}}{{-}}$ 0.2 | 998.1 $\nicefrac{{+}}{{-}}$ 111.5 | 249.521 $\nicefrac{{+}}{{-}}$ 27.9
4 | 1.341 $\nicefrac{{+}}{{-}}$ 0.2 | 755.6 $\nicefrac{{+}}{{-}}$ 67.1 | 188.896 $\nicefrac{{+}}{{-}}$ 16.8
Table 2: Runtime performance with one pyramid level at $512\times 512$ resolution. Smooth iterations set to 2. Max flow | Runtime | Frequency | Throughput
---|---|---|---
(pixels) | (ms) | (Hz) | (Mpix/s)
1 | 1.355 $\nicefrac{{+}}{{-}}$ 0.2 | 747.4 $\nicefrac{{+}}{{-}}$ 65.5 | 186.838 $\nicefrac{{+}}{{-}}$ 16.4
2 | 1.405 $\nicefrac{{+}}{{-}}$ 0.2 | 720.5 $\nicefrac{{+}}{{-}}$ 63.4 | 180.126 $\nicefrac{{+}}{{-}}$ 15.9
4 | 1.501 $\nicefrac{{+}}{{-}}$ 0.2 | 674.3 $\nicefrac{{+}}{{-}}$ 59.0 | 168.584 $\nicefrac{{+}}{{-}}$ 14.7
8 | 1.700 $\nicefrac{{+}}{{-}}$ 0.2 | 594.7 $\nicefrac{{+}}{{-}}$ 51.2 | 148.680 $\nicefrac{{+}}{{-}}$ 12.8
Table 3: Runtime performance with two pyramid levels at $512\times 512$
resolution. Smooth iterations set to [2, 4] for bottom and top level
respectively.

On average, the algorithm runs approximately at 600 Hz with 8 pixels maximum
flow. Although the number of operations per pixel is deterministic given
parameters such as smooth iterations and maximum flow, there is some
variability in the reported runtime, which can be attributed to GPU task
swapping and CPU process scheduling.

###  5.3 Evaluation on real-life data

We performed a set of experiments using real life data captured using a ZED
Stereo Camera. Stereo video is captured at a resolution of $1440\times 720$
per camera at 60 Hz. We used the provided software to extract depth from the
stereo video and use it in our structure flow algorithm. In this case, we use
a spherepix patch of $1024\times 1024$ resolution which approximately matches
the resolution and field of view of the input camera. The algorithm is
configured with two pyramid levels and maximum flow of 8 pixels. Figure 9
displays the results for some selected scenes in driving scenarios with moving
vehicles.

An important difference of the estimated flow compared to conventional optical
flow is the ability to directly distinguish objects getting closer or farther
from the camera, as can be observed in the normal flow plots (column 4). Table
4 shows the average runtime for this configuration varying the maximum flow.
On average, our method achieves frame rates of 172 Hz, which is easily
sufficient for use in control loops for even highly dynamic mobile robotic
vehicles.

Max flow | Runtime | Frequency | Throughput
---|---|---|---
(pixels) | (ms) | (Hz) | (Mpix/s)
1 | 4.683 $\nicefrac{{+}}{{-}}$ 0.2 | 213.889 $\nicefrac{{+}}{{-}}$ 8.2 | 213.889 $\nicefrac{{+}}{{-}}$ 8.2
2 | 4.826 $\nicefrac{{+}}{{-}}$ 0.2 | 207.480 $\nicefrac{{+}}{{-}}$ 6.9 | 207.480 $\nicefrac{{+}}{{-}}$ 6.9
4 | 5.154 $\nicefrac{{+}}{{-}}$ 0.2 | 194.258 $\nicefrac{{+}}{{-}}$ 6.3 | 194.258 $\nicefrac{{+}}{{-}}$ 6.3
8 | 5.815 $\nicefrac{{+}}{{-}}$ 0.2 | 172.149 $\nicefrac{{+}}{{-}}$ 5.1 | 172.149 $\nicefrac{{+}}{{-}}$ 5.1
Table 4: Runtime performance with two pyramid levels at $1024\times 1024$
resolution. Smooth iterations set to [2, 4] for bottom and top level
respectively. Figure 9: Results on real-life stereo video sequences. The
camera captures $1440\times 720$ images at 60 Hz and computes depth. The
brightness and estimated depth are used in the filter algorithm to compute the
underlying structure flow. The normal flow component of the structure flow
offers immediate information about objects getting close to or far from the
camera (rows 1 and 3).

##  6 Conclusions and Future Work

This article introduced the structure flow field as the three-dimensional
vector field encoding the velocity of the surrounding environment relative to
the camera body fixed frame and scaled by the inverse depth. We presented
partial differential equations expressed in spherical geometry modeling the
temporal evolution of the structure flow field, image brightness, depth and
inverse depth fields. These PDEs can be used both to derive predictors to
propagate these quantities forward in time, or as innovation terms to estimate
the structure flow from measured data.

A real-time algorithm for computing structure flow is proposed. The algorithm
is based on a filtering architecture that incrementally computes the flow from
temporal sequences of brightness and depth data. The algorithm achieves frame
rates of the order of 600 Hz at $512\times 512$ pixel resolution and 172 Hz at
$1024\times 1024$ for flow vectors up to 8 pixels. Experimental validation of
the algorithm using simulated high-speed video with ground truth is provided,
as well as results on real-life video sequences using a stereo camera.

Assumption (35) is a practical assumption that allows implementation of the
algorithm in the absence of a calibrated and time-synchronised inertial
measurement unit. Such a system was unavailable for the field tests that we
undertook and would have involved considerable complexity in systems
integration. It is clear that inclusion of the exogenous acceleration terms
would improve results, however, we do not believe that the improvement will be
particularly significant due to the relatively small value of these terms
(Figure 5). More significantly, the current formulation of the filter
algorithm relies on the availability of depth measurements. A natural
extension is to express the cost functions relying on depth measurements in
terms of image disparity for a coupled stereo pair. Such approach will allow
us to formulate the state update equations purely on image data and is the
topic of ongoing investigation.

##  7 Acknowledgments

Thanks to Henri Rebecq from the Robotic and Perception Group at University of
Zurich for sharing the Blender source files of the Urban Canyon dataset. This
research was supported by the Australian Research Council through the
“Australian Centre of Excellence for Robotic Vision” CE140100016.

## References

  * [AM16a] J. D. Adarve and R. Mahony.  A filter formulation for computing real time optical flow.  IEEE Robotics and Automation Letters, 1(2):1192–1199, 2016.
  * [AM16b] J. D. Adarve and R. Mahony.  Spherepix: a data structure for spherical image processing.  IEEE Robotics and Automation Letters (submitted), 2016.
  * [Bar94] J.L. Barron.  System and experiment. performance of optical flow techniques.  Int. J. Comput. Vision, 12:1:43–77, 1994.
  * [BBPW04] Thomas Brox, Andres Brhun, Nils Papenberg, and Joachim Weickert.  High accuracy optical flow estimation based on a theory for warping.  In Proc. European Conf. Comput. Vision, volume 4, pages 25–36, 2004\.
  * [BR09] Silvère Bonnabel and Pierre Rouchon.  Fusion of inertial and visual: a geometrical observer-based approach.  In Proc. Mediterranean Conf. on Intelligent Systems and Automation, volume 1107, page 54, 2009.
  * [BSL+11] Simon Baker, Daniel Scharstein, J. Lewis, Stefan Roth, Michael Black, and Richard Szeliski.  A database and evaluation methodology for optical flow.  Int. J. Comput. Vision, 92:1–31, 2011.
  * [CCC+16] Cesar Cadena, Luca Carlone, Henry Carrillo, Yasir Latif, Davide Scaramuzza, José Neira, Ian D Reid, and John J Leonard.  Past, present, and future of simultaneous localization and mapping: Towards the robust-perception age.  2016\.
  * [GLSU13] A Geiger, P Lenz, C Stiller, and R Urtasun.  Vision meets robotics: The kitti dataset.  Int. J. Robotics Reseach, 32(11):1231–1237, Sep. 2013.
  * [Gon09] Minglun Gong.  Real-time joint disparity and disparity flow estimation on programmable graphics hardware.  Comput. Vision and Image Understanding, 113(1):90 – 100, 2009.
  * [GY06] Minglun Gong and Yee-Hong Yang.  Disparity flow estimation using orthogonal reliability-based dynamic programming.  In Proc. IEEE Int. Conf. Pattern Recognition, volume 2, pages 70–73, 2006.
  * [HB11] S. Hadfield and R. Bowden.  Kinecting the dots: Particle based scene flow from depth sensors.  In Proc. IEEE Conf. Comput. Vision, pages 2290–2295, 2011.
  * [HD07] F. Huguet and F. Devernay.  A variational method for scene flow estimation from stereo sequences.  In Proc. IEEE Conf. Comput. Vision, pages 1–7, Oct 2007.
  * [HHMR12] B. Herisse, T. Hamel, R. Mahony, and F.-X. Russotto.  Landing a vtol unmanned aerial vehicle on a moving platform using optical flow.  IEEE Trans. Robotics, 28(1):77–89, 2012.
  * [HM02] T. Hamel and R. Mahony.  Visual servoing of an under-actuated dynamic rigid-body system: an image-based approach.  IEEE Trans. Robot. Autom, 18(2):187–198, Apr 2002.
  * [HRF13] E. Herbst, X. Ren, and D. Fox.  Rgb-d flow: Dense 3-d motion estimation using color and depth.  In Proc. IEEE Int. Conf. Robotics and Automation, pages 2276–2282, 2013.
  * [HRHM08] B. Herisse, F.-X. Russotto, T. Hamel, and R. Mahony.  Hovering flight and vertical landing control of a vtol unmanned aerial vehicle using optical flow.  In Proc. IEEE/RSJ Int. Conf. Intelligent Robots and Systems, pages 801–806, 2008.
  * [HS81] Berthold K.P. Horn and Brian G. Schunck.  Determining optical flow.  Artificial Intelligence, 17(1):185 – 203, 1981.
  * [KvD87] J. J. Koenderink and Andrea J. van Doorn.  Facts on optic flow.  Biological Cybernetics, 56(4):247–254, 1987.
  * [LK81] Bruce D. Lucas and Takeo Kanade.  An iterative image registration technique with an application to stereo vision.  In Proc. 7th Int. Joint Conf. Artif. Intell., pages 674–679, 1981\.
  * [LPB11] Antoine Letouzey, Benjamin Petit, and Edmond Boyer.  Surface flow from depth and color images.  In British Machine Vision Conference, pages 46–1, 2011.
  * [MB12] C. McCarthy and N. Barnes.  A unified strategy for landing and docking using spherical flow divergence.  IEEE Trans. Pattern Anal. Mach. Intell., 34(5):1024–1031, May 2012\.
  * [MCR15] Maxime Meilland, Andrew I. Comport, and Patrick Rives.  Dense omnidirectional rgb-d mapping of large-scale outdoor environments for real-time localization and autonomous navigation.  J. of Field Robotics, 32(4):474–503, 2015.
  * [MG15] Moritz Menze and Andreas Geiger.  Object scene flow for autonomous vehicles.  In Proc. IEEE Conf. Comput Vision and Pattern Recognition, June 2015\.
  * [PAT96] I. Patras, N. Alvertos, and G. Tziritas.  Joint disparity and motion field estimation in stereoscopic image sequences.  In Proc. IEEE Int. Conf. on Pattern Recognition, volume 1, pages 359–363 vol.1, Aug 1996.
  * [PLBC14] Aurélien Plyer, Guy Le Besnerais, and Frédéric Champagnat.  Massively parallel lucas kanade optical flow for real-time video processing applications.  J. Real-Time Img. Processing, pages 1–18, Apr. 2014.  special issue paper.
  * [QDC12] J. Quiroga, F. Devernay, and J. Crowley.  Scene flow by tracking in intensity and depth data.  In Proc. IEEE Conf. Comput. Vision and Pattern Recognition Workshops, pages 50–57, June 2012.
  * [RF05] Franck Ruffier and Nicolas Franceschini.  Optic flow regulation: the key to aircraft automatic guidance.  Robotics and Autonomous Systems, 50(4):177 – 194, 2005.
  * [RMWF10] Clemens Rabe, Thomas Müller, Andreas Wedel, and Uwe Franke.  Dense, robust, and accurate motion field estimation from stereo image sequences in real-time.  In Proc. European Conf. Comput. Vision, pages 582–595, 2010.
  * [SF11] D. Scaramuzza and F. Fraundorfer.  Visual odometry [tutorial].  IEEE Robot. Autom. Mag., 18(4):80–92, Dec 2011.
  * [SJB02] Hagen Spies, Bernd Jähne, and John L. Barron.  Range flow estimation.  Comput. Vision and Image Understanding, 85(3):209 – 231, 2002.
  * [Sri11] Mandyam V Srinivasan.  Visual control of navigation in insects and its relevance for robotics.  Current Opinion in Neurobiology, 21(4):535 – 543, 2011.
  * [Tho95] James William Thomas.  Numerical partial differential equations: finite difference methods, volume 22.  Springer, 1995.
  * [Tho99] James William Thomas.  Numerical partial differential equations: conservation laws and elliptic equations, volume 3.  Springer Berlin, 1999.
  * [VBR+99] S. Vedula, S. Baker, P. Rander, R. Collins, and T. Kanade.  Three-dimensional scene flow.  In Proc. IEEE Conf. Comput. Vision, volume 2, pages 722–729 vol.2, 1999.
  * [VSR15] Christoph Vogel, Konrad Schindler, and Stefan Roth.  3d scene flow estimation with a piecewise rigid scene model.  Int. J. Comput. Vision, 115(1):1–28, 2015.
  * [WBV+11] Andreas Wedel, Thomas Brox, Tobi Vaudrey, Clemens Rabe, Uwe Franke, and Daniel Cremers.  Stereoscopic scene flow computation for 3d motion understanding.  Int. J. Comput. Vision, 95(1):29–51, 2011.
  * [WSMG+16] Thomas Whelan, Renato F Salas-Moreno, Ben Glocker, Andrew J Davison, and Stefan Leutenegger.  Elasticfusion: Real-time dense slam and light source estimation.  Int. J. Robotics Reseach, 2016.
  * [YBBR93] M. Yamamoto, P. Boulanger, J. A. Beraldin, and M. Rioux.  Direct estimation of range flow on deformable shape from a video rate range camera.  IEEE Trans. Pattern Anal. Mach. Intell., 15(1):82–89, 1993.
  * [ZAR12] N. Zarrouati, E. Aldea, and P. Rouchon.  So(3)-invariant asymptotic observers for dense depth field estimation based on visual data and known camera motion.  In American Control Conference, pages 4116–4123, 2012.
  * [ZRFS16] Zichao Zhang, H. Rebecq, C. Forster, and D. Scaramuzza.  Benefit of large field-of-view cameras for visual odometry.  In Proc. IEEE Int. Conf. Robotics and Automation., pages 801–808, 2016.