11(2.5, -0.1)

This paper has been accepted for publication at the European Conference on
Computer Vision (ECCV), 2024. ©Springer

11institutetext: School of Robotics, Hunan University 22institutetext: School
of Engineering, Westlake University 33institutetext: School of Information
Science and Technology, ShanghaiTech University 44institutetext: TU Berlin,
ECDF, SCIoI Excellence Cluster and Robotics Institute Germany

# Motion and Structure from Event-based
Normal Flow

Zhongyang Ren* 11 0009-0005-5095-3226 Bangyan Liao* 22 0009-0007-7739-4879
Delei Kong 11 0000-0002-5681-587X Jinghang Li 11 0000-0001-6196-6165
Peidong Liu 22 0000-0002-9767-6220 Laurent Kneip 33 0000-0001-6727-6608
Guillermo Gallego 44 0000-0002-2672-9241 Yi Zhou🖂 11 0000-0003-3201-8873

###### Abstract

Recovering the camera motion and scene geometry from visual data is a
fundamental problem in computer vision. Its success in conventional (frame-
based) vision is attributed to the maturity of feature extraction, data
association and multi-view geometry. The emergence of asynchronous (event-
based) cameras calls for new approaches that use raw event data as input to
solve this fundamental problem. State-of-the-art solutions typically infer
data association implicitly by iteratively reversing the event data generation
process. However, the nonlinear nature of these methods limits their
applicability in real-time tasks, and the constant-motion assumption leads to
unstable results under agile motion. To this end, we reformulate the problem
in a way that aligns better with the differential working principle of event
cameras. We show that event-based normal flow can be used, via the proposed
geometric error term, as an alternative to the full (optical) flow in solving
a family of geometric problems that involve instantaneous first-order
kinematics and scene geometry. Furthermore, we develop a fast linear solver
and a continuous-time nonlinear solver on top of the proposed geometric error
term. Experiments on both synthetic and real data show the superiority of our
linear solver in terms of accuracy and efficiency, and its practicality as an
initializer for previous nonlinear solvers. Besides, our continuous-time non-
linear solver exhibits exceptional capabilities in accommodating sudden
variations in motion since it does not rely on the constant-motion assumption.
Our project page can be found at https://nail-hnu.github.io/EvLinearSolver/.

###### Keywords:

Event cameras Geometric model fitting Normal flow

00footnotetext: $\ast$ equal contribution; 🖂 corresponding author
(eeyzhou@hnu.edu.cn).

##  1 Introduction

Neuromorphic event-based cameras are bio-inspired sensors that report only
brightness changes at each pixel location asynchronously. Their high temporal
resolution enables them to sense at exactly the same rate as the scene
dynamics, and thus, event cameras are suitable for dealing with robotic
perception [36, 45, 47, 43, 42], navigation [22, 39, 38, 46], and control [7,
10] tasks that involve agile motion. A general and fundamental problem among
these tasks is to recover the camera’s ego-motion and scene geometry from
event data. One of the bottlenecks for solving this kind of problems at very
high speeds is the lack of methods for efficient and robust data association.

Unlike its standard-vision counterparts, an event camera cannot observe visual
information (_i.e_. intensity changes) when motion is parallel to edges. Such
a partial observability of event data makes optical flow estimation (as well
as feature detection and matching) an ill-posed problem. To this end, many
existing methods [15, 33, 24, 18, 20] for solving the problem of motion-model
fitting on event data propose specific metrics that assess event-data
association established in either an implicit manner or an explicitly model-
driven style. The best fit can be retrieved by iteratively optimizing the data
association, often expressed in terms of event alignment. However, we argue
two flaws of these methods: ($i$) The success of these non-linear solutions
requires a proper initialization. While a brute-force search may be
applicable, it is computationally expensive for real-time application,
specially if the dimension of the parameter is high. Efficient and
deterministic initialization methods are lacking and urgently needed. ($ii$)
The common assumption of constant velocity only holds within a short time
interval. Nevertheless, a time window that is too narrow may contain
insufficient number of events for a valid model fitting. The heuristic of
using a constant number of events may violate the constant velocity
assumption, especially in the case of sudden changes in speed.

In this paper, we rethink the problem formulation in a way that is more
consistent with the differential working principle of event cameras. Unlike
existing pipelines that use events occurring within a time interval to resolve
the relative motion, we attempt to recover the structure and instantaneous
first-order kinematic parameters involved in the motion field equation, using
as input an instantaneous observation, namely event-based normal flow. To this
end, we explore effective constraints between normal-flow observations and
model parameters for solving a family of motion-and-structure estimation
problems, such as the estimation of depth and full optical flow, instantaneous
angular velocity and linear velocity, and differential homography (see Fig.
1).

The contributions of this paper can be summarized as follows:

  * •

A normal flow constraint with a geometric connection to event data. It
overcomes the partial observability issue (_i.e_., errors caused by using
normal flow as replacement for full optical flow) in general problems of
geometric model fitting when normal flow is used as input.

  * •

Two solvers for a family of motion-and-structure estimation problems using as
input sparse event-based normal flow. Our linear solver (can be used with
RANSAC) leads to closed-form and deterministic solutions that can be used as
an initialization to existing nonlinear methods. Our nonlinear solver, under a
continuous-time formulation, can handle sudden speed variations, and thus, is
free of the constant-motion assumption.

  * •

A thorough evaluation on the two proposed solvers, comparing against state-of-
the-art methods, and also including an investigation on the performance of
existing non-linear methods initialized using our linear solver.

_Outline_ : The rest of the paper is organized as follows. We first discuss
the related work in Sec. 2. Then a compact review of preliminary knowledge is
given in Sec. 3. Our approach is presented in Sec. 4, including the normal-
flow constraint, the linear models of a family of motion-and-structure
estimation problems, and the proposed linear and nonlinear solvers.
Experimental evaluation is conducted in Sec. 5 and conclusion are drawn in
Sec. 6.

Figure 1: _Overview of our event-based motion estimation approach_. Normal
flow vectors are computed from the raw event data and are used as input to the
two solvers. The result of the linear solver consists of fitted motion models
$\\{\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{2},\cdots\\}$ at discrete
time instants, which can be used to initialize the continuous-time nonlinear
solver.

##  2 Related Work

###  2.1 Event-based Geometric Model Fitting

The major body of literature on this specific topic follows an iterative
scheme, which typically has its own special metric evaluating quantitatively
the association between events, and solve the model fitting as a nonlinear
optimization problem. Among the first works is the _Contrast Maximization_
method [15, 14], which provides a unified framework for solving a family of
motion-model fitting problems. It searches the optimal model parameters that
ultimately lead to a sharp image of warped events (IWE). The goodness of fit
is evaluated by the alignment of all involved events along point trajectories
on the image plane, and such alignment is assessed by the strength of the
contours of the IWE. The computational complexity of CMax is high because of
iteratively updating of IWE, whose cost scales with the spatial resolution and
the amount of event data.

To lessen the dependency on IWE, Nunes _et al_. proposes an entropy
minimization framework (EMin) [33] that replaces the dispersion metric of [15]
with a similarity measurement between events, enabling to evaluate the
goodness of fit in a higher dimensional space. EMin is, still, too expensive
for real-time application. Liu _et al_. [24] propose a two-step iterative
method for the problem of 3-D rotation estimation using events as input. Given
a rotation estimate, the first step determines event correspondences by
looking for the nearest neighbour in the spatio-temporal domain. The second
step employs a spatio-temporal registration which aligns corresponding events
by minimizing the sum of squared geodesic residual. The least-square nature of
the energy function in the second step leads to a computation friendly
performance. Gu _et al_. [18] present a probabilistic approach for spatio-
temporal alignment of events. The aligned events at a particular pixel are
modeled as a spatio-temporal Poisson point process and the goodness of fit can
be quantified using the likelihood of the data. The optimal model is retrieved
when the likelihood is maximized. More recently, an event-to-map alignment
scheme is proposed in [20], which adopts the time surface (TS) map as a
reference for event alignment. The proposed TS loss function drives later
events to get increasingly closer to pixel positions of corresponding earlier
events, thus, leading to a progressively sharp TS map as iterations proceed.

In summary, all the above-mentioned methods are nonlinear, and thus, a proper
initialization is necessary. Although a brute-force search is widely adopted,
sometimes it can be too expensive to locate the convergence basin in parameter
space. Moreover, [15, 33] may come across a convergence to false-positive
minima, known as event collapse [40, 41]. Hence, a quick and deterministic
solver (e.g., a closed-form solution) is desirable. Besides, most of these
methods assume a constant-motion model, which limits their applicability to
scenarios of sudden speed variations. In addition, we have noticed that a
recent work [17] can jointly estimate the line parameters and the camera’s
linear velocities. However, this method is not a general one due to the
dependence on line structures and the usage of an IMU. Thus, it is not within
the scope of our discussion.

###  2.2 Motion and Structure from Normal Flow

Due to the aperture problem, a local estimate of optical flow recovers only
the partial component of the full flow vector along the direction of image
gradient. Preliminary knowledge about normal flow can be found in the
supplementary material. Recovering motion and structure from normal flow on
intensity images is a well-studied topic. Its main body of literature
typically proposes a qualitatively deterministic solution by which a set (even
infinite) of applicable ego motions are determined via an analysis of normal
flow classification [12, 4, 21]. Only recently, researches devised ways to
formulate the Cheirality constraint as a differentiable term and use it in a
continuous optimization approach [1] or as a loss function for training an
artificial neural network [34].

In event-based perception and navigation, similar methodologies as the frame-
based ones are also witnessed, where normal flow computed from event data is
used as input to determine the Focus of Expansion (FoE) for drones’ course
estimation [11] and obstacle avoidance [6].

Most of these heuristic methods, while effective, are cumbersome for practical
application. In contrast, our method delivers a deterministic way to retrieve
motion and structure from normal flow. Additionally, it is noticed that normal
flow can replace full flow in some specific problems of geometric model
fitting, such as affine models [28, 29]. In this way, the geometric model
fitting problem can be easily solved in closed-form. Nevertheless, these
methods do not consider the observability bias in normal flows, and thus, such
naive substitution may lead to inaccurate model fitting. To this end, we
present a robust linear solver based on a carefully-designed geometric
measurement model that can overcome the partial observability of event-based
normal flow. It is worth mentioning that a recent work [26] also determines
instantaneous linear velocity using as input the normal flow. Our difference
is seen in the scope of applicability. The usage of normal flow is restricted
in the motion flow equation of [26], while our paper extends it to a wider
scope (_i.e_., a family of geometric model fitting problems) by combining
normal flow with _e.g_., differential epipolar geometry, differential
homography, _etc_.

##  3 Preliminaries

To make this paper self-contained, some preliminary knowledge is revisited.
Due to space limitations, here we only list essential concepts. A more
detailed review of preliminary knowledge can be found in the supplementary
material.

####  3.0.1 Motion Field.

Consider a perspective camera moving in a static environment, with
instantaneous angular and linear velocities in the camera frame
${\boldsymbol{\omega}}=(\omega_{x},\omega_{y},\omega_{z})^{\top}$ and
${\boldsymbol{\nu}}=(\nu_{x},\nu_{y},\nu_{z})^{\top}$, respectively. The
motion field at pixel $\mathbf{x}=(x,y)^{\top}$ produced by projecting the 3D
motion at point with depth $Z(\mathbf{x})$ is

$\mathbf{u}=\frac{1}{Z(\mathbf{x})}\mathbf{A}(\mathbf{x}){\boldsymbol{\nu}}+\mathbf{B}(\mathbf{x}){\boldsymbol{\omega}},$
(1)

where $\mathbf{A}(\mathbf{x})$ and $\mathbf{B}(\mathbf{x})$ are $2\times 3$
matrices that solely depend on $\mathbf{x}$. Equation (1) is the motion field
equation [44, Ch.8][25]. Optical flow in rigid scenes aims at estimating the
motion field from visual data, which is only possible in regions that do not
suffer from the aperture problem.

####  3.0.2 Brightness constancy

is the assumption stating that the true motion trajectories $\mathbf{x}(t)$
followed by objects on the image plane are those at which brightness is
constant: $I(\mathbf{x}(t),t)=\text{const}$. Differentiating in time (total
derivative) gives the equation at $(\mathbf{x},t)$: $\nabla
I\cdot\dot{\mathbf{x}}+\partial_{t}I=0$, where
$\dot{\mathbf{x}}=d\mathbf{x}/dt$ is the velocity of image point
$\mathbf{x}(t)$. Given $I$, this equation provides one constraint per point
$(\mathbf{x},t)$ to estimate $\dot{\mathbf{x}}$, now called optical flow.
Without additional equations, only the normal component of the optical flow
can be determined [44],

$\mathbf{n}=-\frac{\partial_{t}I}{\|\nabla I\|^{2}}\nabla I.$ (2)

####  3.0.3 Normal Flow.

Decomposing the optical flow
$\mathbf{u}=\mathbf{u}_{\perp}+\mathbf{u}_{\parallel}$ into its normal and
parallel components to the local image gradient ($\nabla I$), the normal flow
is $\mathbf{n}=\mathbf{u}_{\perp}$, and the dot-product with $\mathbf{n}$
gives

$\mathbf{n}^{\top}\mathbf{u}=\|\mathbf{n}\|^{2},\quad\text{ or equivalently
}\quad(\mathbf{u}-\mathbf{n})\cdot\mathbf{n}=0.$ (3)

####  3.0.4 Differential Homography.

(_Continuous Homography Matrix_ [27, Ch.5]). If all scene points lie on a
plane, their coordinates are related to their image projections via a
homography transformation (which encodes the information of the scene plane
$\\{\mathbf{N},d\\}$ and the camera pose). In this case, the motion field (1)
can be rewritten to adopt the form in homogeneous coordinates

$\hat{\mathbf{u}}(\mathbf{x})=\left(\mathbf{1}-\hat{\mathbf{x}}\mathbf{e}_{3}^{\top}\right)\mathbf{H}_{d}\hat{\mathbf{x}},$
(4)

where
$\mathbf{H}_{d}\doteq-\left([\boldsymbol{\omega}]_{\times}+d^{-1}\boldsymbol{\nu}\mathbf{N}^{\top}\right)$
denotes the differential homography matrix ($[\cdot]_{\times}$ refers to the
cross-product, skew-symmetric matrix), $\mathbf{1}$ is the identity matrix,
$\hat{\mathbf{x}}=(x,y,1)^{\top}$, $\hat{\mathbf{u}}=(u_{x},u_{y},0)^{\top}$
and $\mathbf{e}_{3}=(0,0,1)^{\top}$.

####  3.0.5 Differential Epipolar Geometry.

(_Continuous Epipolar Geometry_ [27, Ch.5]). The epipolar geometry examines
the distance between a potential match to the determined epipolar line, and
can be expressed as
$\hat{\mathbf{x}}^{\prime\top}\mathbf{E}\hat{\mathbf{x}}=0.$ The differential
epipolar geometry can be derived analogously, providing a relationship between
$\mathbf{x},\mathbf{u}$ and the camera velocities
$\boldsymbol{\nu},\boldsymbol{\omega}$. In homogeneous coordinates:

$\hat{\mathbf{u}}(\mathbf{x})^{\top}[\boldsymbol{\nu}]_{\times}\hat{\mathbf{x}}-\hat{\mathbf{x}}^{\top}\mathbf{s}\hat{\mathbf{x}}=0,$
(5)

where
$\mathbf{s}=\frac{1}{2}([\boldsymbol{\nu}]_{\times}[\boldsymbol{\omega}]_{\times}+[\boldsymbol{\omega}]_{\times}[\boldsymbol{\nu}]_{\times})$
is a $3\times 3$ symmetric matrix.

##  4 Methodology

Figure 1 summarizes our approach, grounded upon event-based normal flow. In
this section, we first disclose the normal-flow constraint that can overcome
the partial observability of normal flow when replacing full optical flow with
normal flow in general tasks of geometric model fitting (Sec. 4.1). By
combining the constraint and the theory in Sec. 3, we derive linear equations
for a family of motion-and-structure estimation problems (Sec. 4.2). To solve
these problems, we present a linear solver and a continuous-time nonlinear
solver (Sec. 4.3).

###  4.1 Normal Flow Constraint

Normal flow can be computed easily from the event data, for example by fitting
planes locally to the space-time “cloud” of events on image space [2, 35].
Specifically, assuming that the event data
$e_{k}\doteq\\{\mathbf{x}_{k},t_{k}\\}$ spans locally a surface in the spatio-
temporal domain, the surface can be defined by a function
$\mathrm{\Sigma}_{e}:\mathbb{R}^{2}\to\mathbb{R};\mathbf{x}_{k}\mapsto t_{k}$
that maps the pixel coordinates to the latest event’s timestamp [2, 6]. The
normal flow is parallel to the spatial derivative of $\mathrm{\Sigma}_{e}$,
namely ${\nabla}\mathrm{\Sigma}_{e}(\mathbf{x})$. More specifically, it can be
calculated by

$\mathbf{n}(\mathbf{x})=\frac{{\nabla}\mathrm{\Sigma}_{e}(\mathbf{x})}{\left\|{\nabla}\mathrm{\Sigma}_{e}(\mathbf{x})\right\|^{2}}.$
(6)

Readers may refer to the supplementary material for a more detailed
derivation. Surely there are other methods to extract event-based normal flow
and the solvers we develop in this section are independent of them, i.e., our
solvers can also be fed with the extracted normal flow data from such methods.

Event-based normal flow cannot be used directly as a replacement for the full
flow in (1) to solve for the motion parameters. However, we can leverage the
relationship between both flows (3) to define an effective geometric error
criterion that circumvents the partial observability of normal flow. If the
flow $\mathbf{u}$ depends on some unknown (“state”) variables
$\boldsymbol{\theta}$, we propose estimating them based on how well the normal
flow constraint (3) is satisfied, e.g., by minimizing some function of the
error terms

$\mathbf{n}(\mathbf{x})^{\top}\mathbf{u}(\mathbf{x};\boldsymbol{\theta})-\|\mathbf{n}(\mathbf{x})\|^{2}\doteq\epsilon_{\text{nf}}\,\in\mathbb{R}.$
(7)

(a)

(b)

(c)

(d)

Figure 2: Toy example of applying the proposed normal-flow constraint. (a): A
2D registration task using as input either optical flow or normal flow. The
groundtruth displacement is defined by a global flow of $(1.732,-1)^{\top}$.
(b)-(d): The loss landscape obtained using different geometric measurements in
the registration task. The red dot denotes the resulting displacement, and the
green cross the groundtruth displacement. The bottom-left regions in (b)-(d)
display the corresponding registration results.

To illustrate the validity of the proposed normal-flow constraint over the
naive way of simply using normal flow to replace full optical flow, we use a
2D registration task as toy example. As shown in Fig. 2, the task consists of
registering the solid-line triangle to the dashed-line one using as input
either optical flow or normal flow. We observe that using the proposed normal-
flow constraint as a geometric measurement leads to a registration result that
is almost as accurate as the one given by the full optical flow. However, due
to the partial observability of normal flow, simply using it as a replacement
for full optical flow leads to a biased registration result, as shown in Fig.
2(c). As we see in all following experiments, similar conclusions can be made
in general problems of geometric model fitting.

###  4.2 Models of Motion and Structure

Based on the derivation in Sec. 4.1, we establish linear models to a family of
motion and structure estimation problems using the event-based normal flow as
input. In particular, we first discuss how to recover optical flow and depth
by introducing the normal-flow constraint into the differential epipolar
geometry. Then we show that the combination of the normal-flow constraint with
the motion flow equation (1) results in solutions to the problem of
instantaneous angular velocity estimation, and 6-DoF motion estimation given
depth as prior, respectively. Finally, we show how to compute the differential
homography.

####  4.2.1 Optical Flow.

Given as input the accurate normal flow, and camera velocities, the full
optical flow $\mathbf{u}(\mathbf{x})$ at pixel location $\mathbf{x}$ can be
calculated by solving the following linear system obtained from the
combination of the normal-flow constraint (3) and the differential epipolar
geometry (5):

$\begin{pmatrix}{\mathbf{n}}(\mathbf{x})^{\top}\\\\[2.15277pt]
-\hat{\mathbf{x}}^{\top}[\boldsymbol{\nu}]_{\times}\end{pmatrix}{\mathbf{u}}(\mathbf{x})=\begin{pmatrix}\|\mathbf{n}(\mathbf{x})\|^{2}\\\\[2.15277pt]
\hat{\mathbf{x}}^{\top}\mathbf{s}\hat{\mathbf{x}}\end{pmatrix}.$ (8)

####  4.2.2 Depth Estimation.

Similarly, if the normal flow and camera velocities are known, the scene depth
$Z(\mathbf{x})$ can be determined by combining the normal-flow constraint (3)
with the motion field eq. (1), yielding the closed-form solution:

$Z(\mathbf{x})=\frac{\mathbf{n}(\mathbf{x})^{\top}\mathbf{A}(\mathbf{x})\boldsymbol{\nu}}{\|\mathbf{n}(\mathbf{x})\|^{2}-{\mathbf{n}}(\mathbf{x})^{\top}\mathbf{B}(\mathbf{x})\boldsymbol{\omega}}.$
(9)

####  4.2.3 Angular Velocity Estimation.

Assuming that an event camera undergoes a pure rotational motion, the optical
flow $\mathbf{u}$ in (1) can be simplified as

$\mathbf{u}(\mathbf{x})=\mathbf{B}(\mathbf{x})\boldsymbol{\omega}.$ (10)

By substituting (10) into (3), we have the linear equation in
$\boldsymbol{\omega}$:

$\mathbf{n}(\mathbf{x})^{\top}\mathbf{B}(\mathbf{x})\boldsymbol{\omega}=\|\mathbf{n}(\mathbf{x})\|^{2}.$
(11)

Since each normal flow provides only one constraint, the minimal problem
requires three normal flow vectors as input.

####  4.2.4 6-DoF Motion Estimation.

Assume the event camera undergoes a general 6-DoF motion with known depth of
the observed scene, we can solve the first-order kinematics, namely
instantaneous linear and angular velocities, by solving the following linear
system obtained from the normal-flow constraint (3) and the motion-flow
equation (1):

$\mathbf{n}(\mathbf{x})^{\top}\mathbf{D}(\mathbf{x})\begin{bmatrix}\boldsymbol{\nu}\\\
\boldsymbol{\omega}\end{bmatrix}=\|\mathbf{n}(\mathbf{x})\|^{2},$ (12)

where
$\mathbf{D}(\mathbf{x})=\begin{bmatrix}\frac{1}{Z(\mathbf{x})}\mathbf{A}(\mathbf{x})&\mathbf{B}(\mathbf{x})\end{bmatrix}$
is the $2\times 6$ feature sensitivity matrix [8, 5]. The minimal problem in
this context requires six normal flow measurements to determine the 6-DoF
motion parameters.

####  4.2.5 Motion and Structure Estimation in Planar Scenes.

For a moving camera that observes a planar scene, the motion and structure
parameters are compatibly encoded in a homography matrix. In our case, we
attempt to recover the differential homography $\mathbf{h}_{d}$ (the vector
form of $\mathbf{H}_{d}$ in (4)) because the instantaneous linear and angular
velocities are considered. By expanding (4) and rearranging equations, we have

$\mathbf{u}(\mathbf{x})=\mathbf{C}(\mathbf{x})\mathbf{h}_{d},$ (13)

where $\mathbf{C}(\mathbf{x})$ is a $2\times 9$ matrix that solely depends on
the pixel location $\mathbf{x}$. Subsequently, substituting (3) into (13),
gives:

$\mathbf{n}(\mathbf{x})^{\top}\mathbf{C}(\mathbf{x})\mathbf{h}_{d}=\|\mathbf{n}(\mathbf{x})\|^{2}.$
(14)

As discussed in [49], retrieving the motion and structure parameters from the
decomposition of $\mathbf{h}_{d}$ is not trivial and beyond the scope of this
paper. Therefore, we only focus on the numerical solution of $\mathbf{h}_{d}$.

###  4.3 Our Solvers

According to the motion-and-structure models derived in Sec. 4.2, we have
established certain constraints between the model parameters and the normal
flow observations. This allows us to simply develop linear solvers for each
specific problem. Besides, we further propose a continuous-time nonlinear
solver that can effectively handle the presence of sudden variations in motion
and does not rely on the constant-motion assumption. To convey our solvers
concisely, we generalize the motion-and-structure models into a unified form,
as

$\mathbf{n}(\mathbf{x})^{\top}\mathbf{O}(\mathbf{x})\boldsymbol{\theta}=\|\mathbf{n}(\mathbf{x})\|^{2},$
(15)

where $\mathbf{O}(\mathbf{x})$ is the system matrix111$\mathbf{O}(\mathbf{x})$
degenerates to 1 in the problem of optical flow estimation. and
$\boldsymbol{\theta}$ are the motion/structure parameters.

####  4.3.1 Linear Solver.

Given as input a set (_e.g_., $\mathrm{K}$) of sparse normal-flow estimates
obtained at the same time222An approximation can be made based on the
assumption of constant camera motion. In this way, the normal flow estimates
can be collected within a short time interval., our linear solver works out
the following (overdetermined) linear system of equations:

$\begin{pmatrix}\mathbf{n}(\mathbf{x}_{1})^{\top}\mathbf{O}(\mathbf{x}_{1})\\\
\vdots\\\
\mathbf{n}(\mathbf{x}_{\mathrm{K}})^{\top}\mathbf{O}(\mathbf{x}_{\mathrm{K}})\end{pmatrix}\boldsymbol{\theta}=\begin{pmatrix}\|\mathbf{n}(\mathbf{x}_{1})\|^{2}\\\
\vdots\\\ \|\mathbf{n}(\mathbf{x}_{\mathrm{K}})\|^{2}\end{pmatrix}.$ (16)

For robust estimation, RANSAC [13] is adopted to deal with noise and outliers.

####  4.3.2 Continuous-Time Nonlinear Solver.

To handle the case of a sudden variation in speed that may violate the
constant-motion assumption, we further propose a continuous-time nonlinear
solver. In particular, we employ cubic B-splines to represent the continuous-
time model parameters as

${\boldsymbol{\theta}}(u(t))=\sum_{i=0}^{3}B_{i,3}(u(t)){\mathbf{c}}_{i},$
(17)

where $B_{i,k}$ denotes the basis function, $i$ the index of control points,
$k$ (=3) the order of the spline and ${\mathbf{c}}_{i}\in\mathbb{R}^{N}$ the
corresponding control point defined in the $N$-dimensional space of model
parameters. Examples of cubic B-splines applied to event camera motion
estimation are found in [31, 19]. The function $u(\cdot)$ is a normalization
operator, which transfers time $t$ to the spline’s parameter domain by means
of basis translation [9]. In this way, the model parameters at any time
instant $t$ can be obtained via an interpolation using the known control
points.

The formulation of the continuous-time model fitting problem goes as follows.
The unknown variable consists of all the involved control points,
$\boldsymbol{\mathcal{X}}\doteq[\mathbf{c}_{0},\mathbf{c}_{1},\cdots\mathbf{c}_{n}]$,
where $n$ denotes the total number of control points. By introducing the
continuous-time representation of model parameters to (15), we finally create
the following objective function:

$\boldsymbol{\mathcal{X}}^{\ast}=\arg\min_{\boldsymbol{\mathcal{X}}}\sum_{e_{k}\in\mathcal{E}}\left\|\mathbf{n}(\mathbf{x}_{k})^{\top}\mathbf{O}(\mathbf{x}_{k}){\boldsymbol{\theta}}(u(t_{k}))-\|\mathbf{n}(\mathbf{x}_{k})\|^{2}\right\|^{2},$
(18)

where $\mathcal{E}$ denotes the event set in which the norm flow observation
is available for every involved event $e_{k}$. Note that this continuous-time
formulation aligns better with the asynchronous nature of event-based normal
flow, enabling a state estimation over a larger temporal window that may lead
to a better signal-to-noise ratio. The linear solver can be used to initialize
the continuous-time nonlinear solver. An illustrative explanation of our
solvers is given in Fig. 1.

##  5 Experiments

To justify the effectiveness of the proposed algorithms, experiments are
conducted on both simulation and real data. First, we introduce the datasets
used and the applied metrics for evaluation (Sec. 5.1). Then we report the
results of our linear solver and continuous-time nonlinear solver (Sec. 5.2
and Sec. 5.3). Finally, a discussion on the computational complexity and the
limitation of our method is provided in Sec. 5.4.

###  5.1 Experimental Setup

Our evaluation utilizes four datasets, including the indoor event-camera
dataset (ECD) captured with a DAVIS240C camera and motion-capture system for
groundtruth, and a high-resolution dataset with two rotational motion
sequences using an iniVation DAVIS346 camera and Xsens IMU. The VECtor
dataset, captured with a Prophesee Gen 3 camera and an accompanying RGB-D
depth camera, is used for 6-DoF motion tracking. Additionally, a synthetic
dataset generated in Blender is employed for motion-and-structure pattern
analysis. Evaluation metrics for kinematics estimation include average error
(AE) and root mean squared error (RMSE), while differential homography
estimation is assessed using the Frobenius norm. Best results are in bold, and
second-best are underlined. The supplementary material provides detailed
descriptions of datasets and metrics, as well as implementation details of our
methods.

(a) Estimated optical flow over a map of event accumulation.

(b) Reconstruction result displayed over groundtruth structures (a v-shape
roof) and viewed from different perspectives.

Figure 3: Qualitative result of optical flow and depth on sequence
_three_wall_translation_ of our synthetic data. Note that our method returns
sparse results.

###  5.2 Evaluation of the Linear Solver

In the following evaluation, we qualitatively and quantitatively assess the
performance of our linear method (Ours) in each specific task, and compare
against different state-of-the-art (SOTA) methods (CMax[15] and Pro-STR[20])
and their customized version initialized using our linear solver.

####  5.2.1 Optical Flow and Depth Estimation.

Due to a lack of benchmark for extensive evaluation, we qualitatively evaluate
our method on the tasks of optical flow and depth estimation. As illustrated
in Fig. 3(a), our linear solver demonstrates the capability to accurately
reconstruct full flow using normal flow and known motion parameters. An
example of the depth estimation result is given in Fig. 3(b), showing that the
recovered structure aligns well with the groundtruth.

####  5.2.2 Angular Velocity Estimation.

As shown in Tab. 1, we see that our linear solver outperforms CMax and Pro-
STR. Besides, we find the combination of our linear solver and CMax leads to
the best result, indicating that the proposed linear solver is a complementary
method to existing nonlinear methods.

Table 1: Evaluation of our linear solver on the task of angular velocity
estimation.

Method _shapes_rotation_[32] _ground_rotation_ (our data) _patterns_rotation_
(our data) $e_{\omega}$ (deg/s) $\text{RMSE}_{\omega}$ (deg/s) $e_{\omega}$
(deg/s) $\text{RMSE}_{\omega}$ (deg/s) $e_{\omega}$ (deg/s)
$\text{RMSE}_{\omega}$ (deg/s) CMax[15] 15.60 26.38 9.48 15.58 5.95 9.07 Pro-
STR[20] 90.61 145.81 31.07 44.74 10.47 11.88 Ours 15.36 23.98 3.30 8.01 4.70
6.08 Ours+CMax[15] 10.12 14.29 2.85 3.74 0.35 0.73 Ours+Pro-STR[20] 19.23
40.44 9.13 17.07 4.31 5.23 Method _dynamic_rotation_[32] _boxes_rotation_ (our
data) _cubes_rotation_ (our data) $e_{\omega}$ (deg/s) $\text{RMSE}_{\omega}$
(deg/s) $e_{\omega}$ (deg/s) $\text{RMSE}_{\omega}$ (deg/s) $e_{\omega}$
(deg/s) $\text{RMSE}_{\omega}$ (deg/s) CMax[15] 16.14 26.48 8.18 13.56 6.77
10.62 Pro-STR[20] 54.64 85.63 22.72 32.05 11.07 12.79 Ours 9.41 18.06 2.60
4.78 2.31 3.02 Ours+CMax[15] 7.71 10.85 2.64 3.80 0.49 0.77 Ours+Pro-STR[20]
18.87 32.65 3.49 5.87 2.04 2.76

Table 2: Evaluation of our linear solver for differential homography
estimation.

Method _board_slow_[16] _patterns_rotation_ _patterns_6dof_ F-Norm F-Norm
F-Norm CMax[15] 1.16 1.33 1.27 Pro-STR[20] 2.07 1.22 1.37 Ours 0.24 0.20 0.43
Ours+CMax[15] 0.23 0.15 0.41 Ours+Pro-STR[20] 0.19 0.17 0.42

Table 3: Evaluation of our linear solver on the task of 6-DoF motion tracking.

Method _corner_slow_[16] _patterns_6dof_ (our data) $e_{\omega}$ (deg/s)
$\mathrm{RMSE}_{\omega}$ (deg/s) $e_{\upsilon}$ (m/s)
$\mathrm{RMSE}_{\upsilon}$ (m/s) $e_{\omega}$ (deg/s) $\mathrm{RMSE}_{\omega}$
(deg/s) $e_{\upsilon}$ (m/s) $\mathrm{RMSE}_{\upsilon}$ (m/s) CMax[15] 12.39
14.39 0.34 0.41 12.94 15.68 1.02 1.23 Pro-STR[20] 5.19 7.84 0.13 0.22 11.21
13.49 1.33 1.49 Ours 1.31 1.79 0.01 0.02 0.56 0.89 0.12 0.16 Ours+CMax[15]
2.19 2.82 0.11 0.17 0.55 0.87 0.18 0.30 Ours+Pro-STR[20] 1.06 1.34 0.04 0.05
1.07 1.56 0.12 0.15 Method _mountain_normal_[16] _cubes_6dof_ (our data)
$e_{\omega}$ (deg/s) $\mathrm{RMSE}_{\omega}$ (deg/s) $e_{\upsilon}$ (m/s)
$\mathrm{RMSE}_{\upsilon}$ (m/s) $e_{\omega}$ (deg/s) $\mathrm{RMSE}_{\omega}$
(deg/s) $e_{\upsilon}$ (m/s) $\mathrm{RMSE}_{\upsilon}$ (m/s) CMax[15] 11.01
14.17 0.44 0.52 10.74 14.33 0.78 1.05 Pro-STR[20] 5.23 10.53 0.14 0.25 9.75
11.40 1.29 1.44 Ours 1.28 1.83 0.07 0.09 0.79 1.02 0.11 0.14 Ours+CMax[15]
3.84 6.30 0.21 0.30 1.14 1.72 0.12 0.14 Ours+Pro-STR[20] 1.18 1.58 0.06 0.09
1.05 1.33 0.11 0.14

(a) Results of the continuous-time angular velocity estimator.

(b) Comparison against two state-of-the-art nonlinear solutions using data
segments that involve sudden variations of angular velocity.

Figure 4: Angular velocity estimation using the continuous-time nonlinear
solver (Ours).

####  5.2.3 Differential Homography Estimation.

The result is consistent with that in the above task. Table 2 reveals that our
linear method outperforms CMax and Pro-STR in terms of numerical accuracy. The
combination of our linear solver and state-of-the-art nonlinear solver leads
to the best result.

####  5.2.4 6-DoF Motion Tracking.

Given known depth, we recover the linear and angular camera velocities on the
task of 6-DoF motion estimation. As shown in Tab. 3, our linear solver reports
more accurate results than both CMax and Pro-STR. In addition, we surprisingly
find that the linear solver can sometimes outperform the combination of itself
with the two nonlinear solvers. This is, according to our investigation,
attributed to the violation of constant-motion assumption in some sequences,
which is held for the two nonlinear solvers.

###  5.3 Evaluation of the Continuous-Time Solver

To justify the claimed advantage of the proposed continuous-time nonlinear
solver, we use the exemplary task of angular velocity estimation. The sequence
used (_shapes_rotation_ [32]) captures a purely rotating event camera, and
some parts of the data involves aggressive motion. As illustrated in Fig.
4(a), our nonlinear solver gives accurate estimation results with a high
stability. Moreover, we select three segments that witness a sudden variation
of angular speed within a short time interval and compare our continuous-time
nonlinear solver against CMax and Pro-STR. Note that all of them are
initialized by our linear method for a fair comparison. As Fig. 4(b) shows,
our results are more consistent with the groundtruth angular velocity than the
other two methods, which are less accurate due to the dependence on the
constant-motion assumption.

###  5.4 Discussion

####  5.4.1 Computational Complexity.

The computational complexity of the linear solver depends on the specific task
of model fitting. For optical flow and depth estimation, it is linearly
proportional to the size of the input event flow, and thus, the complexity is
$\mathcal{O}(N)$. For the other tasks, the time cost predominantly hinges on
the total number of RANSAC iterations and the computations performed per
iteration. In each iteration, we typically solve a linear overdetermined
system via QR decomposition at a complexity of $\mathcal{O}(cN^{2})$, where
$c$ is the cardinality of a minimal problem (_e.g_., $c=3$ for the problem of
angular velocity estimation).

####  5.4.2 Runtime.

All experiments are conducted on a standard computer with Intel i7-13700F CPU,
16G memory and implemented with MATLAB (version R2022a). Using the 6-DoF
tracking task as an example, we measure the runtime of the proposed linear
solver and continuous-time (CT) solver, as reported in Tab. 4.

Table 4: Runtime of 6-DoF tracking on sequence _mountain_normal_
(640$\times$480 px).

Component Runtime Normal flow extraction 20 ms (for extracting 2000 normal
flow vectors) Linear solver 25.97 ms CT solver 28.29 ms

####  5.4.3 Limitations.

To make it clear, we emphasize that our method, like other compared methods,
does not resolve the motion and structure parameters simultaneously. It should
also be noted that the success of our method relies heavily on the quality of
the input normal flow. Particularly, our normal flow extractor (including some
other open source solutions for computing normal flow from event data) may
fail in the presence of densely repetitive textures, leading to a failure of
the modeling fitting scheme. To assess the sensitivity of our method to the
quality of event-based normal flow, we propose an analysis on numerical
stability using simulation data in the supplementary material.

##  6 Conclusions

This work presents a solution to the general problem of geometric model
fitting with event data. Considering the differential working principle of
event cameras, our method establishes specific constraints between
instantaneous motion-and-structure parameters and event-based normal flow.
Experiments show that the proposed linear solver leads to closed-form and
deterministic solutions, which can be used to initialize existing state-of-
the-art nonlinear methods. The proposed continuous-time nonlinear solver,
furthermore, demonstrates the capability to handle sudden speed variations,
obviating the dependence on the constant-motion assumption. Finally, we hope
this work can inspire future research on event-based geometric model fitting,
a paramount topic to many downstream techniques, such as event-based 3-D
vision, SLAM, and scene segmentation.

## Acknowledgments

We thank Dr. Yi Yu for proof reading. This work was supported by the National
Key Research and Development Project of China under Grant 2023YFB4706600.
Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany’s Excellence Strategy – EXC 2002/1 “Science of
Intelligence” – project number 390523135.

## SUPPLEMENTARY MATERIAL

## Overview

In this supplementary material we provide the following:

  * •

A detailed review of preliminaries involved, including motion field,
differential geometry and event-based normal flow (Sec. 0.A).

  * •

Additional details about differential homography (Sec. 0.B).

  * •

Details of dataset and evaluation metrics utilized (Sec. 0.C).

  * •

Implementation details of the proposed algorithms (Sec. 0.D).

  * •

More evaluations on the proposed methods, including an analysis of numerical
stability, and the performance of the continuous-time nonlinear solver in
textured scenes. (Sec. 0.E).

##  Appendix 0.A Supplement to Preliminaries

We supplement Sec. 3 of the main paper by providing a detailed review on three
important concepts, including _motion field_ (Sec. 0.A.1), _differential
geometry_ across multiple views (Sec. 0.A.2), and _event-based normal flow_
(Sec. 0.A.3).

###  0.A.1 Motion Field

Consider a visual observer moving in a static environment, and let its
instantaneous angular velocity and linear velocity be denoted in the
observer’s body frame as
${}^{\mathtt{B}}{\boldsymbol{\omega}}=(\omega_{x},\omega_{y},\omega_{z})^{\top}$
and ${}^{\mathtt{B}}{\boldsymbol{\nu}}=(\nu_{x},\nu_{y},\nu_{z})^{\top}$,
respectively. The scene flow of a 3D point
${}^{\mathtt{B}}{\mathbf{P}}=(X,Y,Z)^{\top}$ is

$\displaystyle\dot{{}^{\mathtt{B}}{\mathbf{P}}}=-{}^{\mathtt{B}}{\boldsymbol{\omega}}\times{}^{\mathtt{B}}{\mathbf{P}}-{}^{\mathtt{B}}{\boldsymbol{\nu}}.$
(19)

Let $\mathbf{x}=(x,y)^{\top}=(X/Z,Y/Z)^{\top}$ be the image of $\mathbf{P}$
represented in calibrated camera coordinates. Projecting the scene flow vector
onto the image plane gives us the relationship between the optical flow
$\mathbf{u}$ (in calibrated coordinates) and the dynamics of the observer as

$\displaystyle\mathbf{u}=\frac{1}{Z}\mathbf{A}(\mathbf{x})\,{}^{\mathtt{B}}{\boldsymbol{\nu}}+\mathbf{B}(\mathbf{x})\,{}^{\mathtt{B}}{\boldsymbol{\omega}},$
(20) $\displaystyle\mathbf{A}(\mathbf{x})=\begin{bmatrix}-1&0&x\\\
0&-1&y\end{bmatrix},$ (21)
$\displaystyle\mathbf{B}(\mathbf{x})=\begin{bmatrix}xy&-(1+{x}^{2})&y\\\
(1+{y}^{2})&-xy&-x\end{bmatrix},$ (22)

where (20) is acknowledged as the motion field equation [25]. The depth $Z$ is
different for every pixel $\mathbf{x}$, hence we often write
$Z=Z(\mathbf{x})$.

###  0.A.2 Multi-View Differential Geometry

The theory of differential geometry across multiple views can be derived by
taking the temporal derivatives of the standard multi-view geometry [27]. It
examines the relationship between the optical flow and the instantaneous
kinematics across multiple views. We focus on two specific two-view geometric
properties: 1) the differential homography [49], and 2) the differential
epipolar geometry [48].

####  0.A.2.1 Differential Homography.

The standard homography is a projective transformation that describes the
image motion between two views of a camera observing a planar scene. Such a
mapping function, induced by the planar scene, can be expressed as

$\hat{\mathbf{x}}^{\prime}=\mathbf{H}\hat{\mathbf{x}},$ (23)

where $\hat{\mathbf{x}}$ and $\hat{\mathbf{x}}^{\prime}$ denote the
corresponding pixels in terms of homogeneous calibrated image coordinates, and
$\mathbf{H}\in\mathbb{R}^{3\times 3}$ is the homography matrix. The induced
homography is a function of the relative pose $\\{\mathtt{R},\mathbf{t}\\}$
and the structure parameters $\\{\mathbf{N},d\\}$ (i.e., parameters of the
plane in the scene), and can be expressed as
$\mathbf{H}\doteq\mathtt{R}-\frac{1}{d}\mathbf{t}\mathbf{N}^{\top}$. The
differential version replaces the potentially large displacement of the
transformation $\mathbf{x}\mapsto\mathbf{x}^{\prime}$ by its first-order
continuous-time approximation $\mathbf{x}\mapsto\mathbf{x}+\Delta
t\hat{\mathbf{u}}$, with point velocity given by:

$\hat{\mathbf{u}}(\mathbf{x})=\left(\mathbf{1}-\hat{\mathbf{x}}\mathbf{e}_{3}^{\top}\right)\mathbf{H}_{d}\hat{\mathbf{x}}.$
(24)

Here,
$\mathbf{H}_{d}\doteq-\left([\boldsymbol{\omega}]_{\times}+\frac{1}{d}\boldsymbol{\nu}\mathbf{N}^{\top}\right)$
denotes the differential homography matrix ($[\cdot]_{\times}$ denotes the
cross-product, skew-symmetric matrix), $\mathbf{e}_{3}=(0,0,1)^{\top}$,
$\hat{\mathbf{x}}=(x,y,1)^{\top}$, and
$\hat{\mathbf{u}}=(u_{x},u_{y},0)^{\top}$ the optical flow (in homogeneous
coordinates).

####  0.A.2.2 Differential Epipolar Geometry.

As a relatively relaxed constraint compared to the homography, the standard
epipolar geometry only defines a potential direction for searching the
correspondence pair
$\hat{\mathbf{x}}\leftrightarrow\hat{\mathbf{x}}^{\prime}$. Such a constraint
evaluates the distance between a potential match to the determined epipolar
line, and can be expressed as

$\hat{\mathbf{x}}^{\prime\top}\mathbf{E}\hat{\mathbf{x}}=0.$ (25)

The differential epipolar geometry can be derived similarly, yielding

$\hat{\mathbf{u}}(\mathbf{x})^{\top}[\boldsymbol{\nu}]_{\times}\hat{\mathbf{x}}-\hat{\mathbf{x}}^{\top}\mathbf{s}\hat{\mathbf{x}}=0,$
(26)

where
$\mathbf{s}=\frac{1}{2}([\boldsymbol{\nu}]_{\times}[\boldsymbol{\omega}]_{\times}+[\boldsymbol{\omega}]_{\times}[\boldsymbol{\nu}]_{\times})$.
Note that (26) can be derived directly from (20) by eliminating the depth
information.

###  0.A.3 Event-Based Normal Flow

####  0.A.3.1 Normal Flow.

While the motion field is the true projection of the 3D scene velocities on
the image plane, its determination from visual quantities measured at the
image plane is called “optical flow”. The optical flow is the apparent motion
of patterns on the image plane. As is well known, the optical flow can only be
a good approximation for the true motion (i.e., the motion field) at image
points that are unambiguous to track, i.e., that are surrounded by edge
patterns that accurately determine their location in both image plane
directions ($x$ and $y$). A common assumption that is leveraged to compute
optical flow is “brightness constancy”, that is, the approximation that image
brightness remains constant as the image point $\mathbf{x}$ moves. This is an
approximation for points corresponding to matt and Lambertian objects; it is
not a sensible assumptions for points on shiny and specular objects.

Letting $I(\mathbf{x},t)$ be the brightness function on the image plane as a
function of time, the brightness constancy assumption can be compactly written
as:

$I(\mathbf{x}(t),t)=\text{const}.$ (27)

A first-order Taylor expansion gives

$I(\mathbf{x}+\Delta\mathbf{x},t+\Delta t)\approx I(\mathbf{x},t)+(\nabla
I(\mathbf{x},t))^{\top}\dot{\mathbf{x}}\Delta
t+\partial_{t}I(\mathbf{x},t)\Delta t,$ (28)

where $\dot{\mathbf{x}}\equiv\dot{\mathbf{x}}(t)=(dx/dt,dy/dt)^{\top}$ is the
velocity of image point $\mathbf{x}(t)$,
$\nabla=(\partial_{x},\partial_{y})^{\top}$ are the spatial derivatives and
$\partial t$ is the temporal derivative. Moving terms around,

$I(\mathbf{x}+\Delta\mathbf{x},t+\Delta t)-I(\mathbf{x},t)\approx(\nabla
I(\mathbf{x},t))^{\top}\dot{\mathbf{x}}\Delta
t+\partial_{t}I(\mathbf{x},t)\Delta t.$ (29)

Brightness constancy states that for the true motion curves $\mathbf{x}(t)$,
the intensity remains the same, $I(\mathbf{x}+\Delta\mathbf{x},t+\Delta
t)\approx I(\mathbf{x},t)$, and therefore (since $\Delta t>0$)

$(\nabla
I(\mathbf{x},t))^{\top}\dot{\mathbf{x}}+\partial_{t}I(\mathbf{x},t)\approx 0.$
(30)

Assuming that the brightness $I$ is known (i.e., measured by a camera), (30)
is one equation in two unknowns ($\dot{\mathbf{x}}\in\mathbb{R}^{2}$). More
equations are needed to be able to determine the velocity $\dot{\mathbf{x}}$
at point $(\mathbf{x},t)$ in space-time. Without additional equations or
information, all that can be determined is the velocity component that is
parallel to $\nabla I(\mathbf{x},t)$ (i.e., perpendicular to the edge) due to
the dot product operation. That is, decomposing
$\dot{\mathbf{x}}=\dot{\mathbf{x}}_{\parallel}+\dot{\mathbf{x}}_{\perp}$ into
its parallel and perpendicular components to the local edge $\nabla I$,
respectively, we have $\nabla
I(\mathbf{x},t)^{\top}\dot{\mathbf{x}}_{\parallel}=0$. Substituting in (30)
gives:

$\displaystyle 0$ $\displaystyle\approx$ $\displaystyle(\nabla
I(\mathbf{x},t))^{\top}\left(\dot{\mathbf{x}}_{\parallel}+\dot{\mathbf{x}}_{\perp}\right)+\partial_{t}I(\mathbf{x},t)$
(31) $\displaystyle=$ $\displaystyle\underbrace{(\nabla
I(\mathbf{x},t))^{\top}\dot{\mathbf{x}}_{\parallel}}_{0}+(\nabla
I(\mathbf{x},t))^{\top}\dot{\mathbf{x}}_{\perp}+\partial_{t}I(\mathbf{x},t)$
(32) $\displaystyle=$ $\displaystyle(\nabla
I(\mathbf{x},t))^{\top}\dot{\mathbf{x}}_{\perp}+\partial_{t}I(\mathbf{x},t)$
(33)

From here, we can work out
$\mathbf{n}\equiv\dot{\mathbf{x}}_{\perp}\propto\nabla I$. We know that

$\mathbf{n}\equiv\dot{\mathbf{x}}_{\perp}=\|\dot{\mathbf{x}}_{\perp}\|\frac{\nabla
I}{\|\nabla I\|}.$ (34)

Substituting in (33) (and omitting the evaluation point $(\mathbf{x},t)$, for
clarity),

$(\nabla I)^{\top}\|\dot{\mathbf{x}}_{\perp}\|\frac{\nabla I}{\|\nabla
I\|}=-\partial_{t}I,$ (35)

which gives

$\|\dot{\mathbf{x}}_{\perp}\|=-\partial_{t}I\frac{\|\nabla I\|}{(\nabla
I)^{\top}\nabla I}=-\partial_{t}I\frac{1}{\|\nabla I\|}.$ (36)

Substituting in (34), it finally gives

$\mathbf{n}(\mathbf{x},t)=-\frac{\partial_{t}I(\mathbf{x},t)}{\|\nabla
I(\mathbf{x},t)\|^{2}}\nabla I(\mathbf{x},t).$ (37)

This component of the optical flow $\dot{\mathbf{x}}(\mathbf{x},t)$, which is
perpendicular to the edge $\nabla I(\mathbf{x},t)$, is called _normal flow_.
While (37) is well-known in conventional (frame-based) computer vision, it is
not directly applicable to event cameras because the latter do not directly
provide the means to compute all terms required (the spatial and temporal
derivatives of the brightness). Next, we present how to compute the normal
flow in the case of event cameras.

####  0.A.3.2 Computation of Event-based Normal Flow.

The output of an event camera is a stream of events, where each event
$\boldsymbol{e}_{k}\doteq(\mathbf{x}_{k},t_{k},p_{k})$ consists of the space-
time coordinates $(\mathbf{x}_{k},t_{k})$ at which the intensity change of
predefined size happened and the sign of the change (_i.e_., polarity
$p_{k}\in\\{+1,-1\\}$). Following the definition in [2], we utilize the
differential mapping function
$\mathrm{\Sigma}_{e}:\mathbb{R}^{2}\to\mathbb{R}$ that maps the pixel
coordinate to the latest event’s timestamp, i.e., $\mathbf{x}\mapsto
t_{\text{last}}(\mathbf{x})$. $\mathrm{\Sigma}_{e}$ is a 2D image, also called
a time map; and it is often referred to as _time surface_ (TS) because when
interpreting the image as an elevation map, a moving edge produces events and
such points
$(\mathbf{x}_{k},t_{\text{last}}(\mathbf{x}_{k}))\subset\mathbb{R}^{3}$
approximately form a surface (ignoring the discrete pixel lattice) [2, 23]. It
also goes by the name of _Surface of Active Events_ (SAE).

_Lemma._ _Given the mapping function $\mathrm{\Sigma}_{e}$ at pixel coordinate
$\mathbf{x}$ and its finite neighbourhood, the direction of the first-order
partial derivative
${\nabla}\mathrm{\Sigma}_{e}(\mathbf{x})=\frac{\partial\mathrm{\Sigma}_{e}}{\partial\mathbf{x}}(\mathbf{x})$
is identical to that of the normal flow at $\mathbf{x}$_.

_Proof._ The directional derivative
${\nabla}_{\mathbf{d}}\mathrm{\Sigma}_{e}(\mathbf{x})$ can be defined by the
dot product of direction $\mathbf{d}$ and
${\nabla}\mathrm{\Sigma}_{e}(\mathbf{x})$, as
${\nabla}_{\mathbf{d}}\mathrm{\Sigma}_{e}(\mathbf{x})\doteq\mathbf{d}\cdot{\nabla}\mathrm{\Sigma}_{e}(\mathbf{x})$,
whose maximum is achieved when
$\mathbf{d}_{\max}=\frac{{\nabla}\mathrm{\Sigma}_{e}(\mathbf{x})}{\left\|{\nabla}\mathrm{\Sigma}_{e}(\mathbf{x})\right\|}$.
The direction $\mathbf{d}_{\max}$ is proved to be perpendicular to the level
set $\mathcal{S}\doteq\\{\mathbf{x}|\mathrm{\Sigma}_{e}(\mathbf{x})=t\\}$
which corresponds to the latent edge pixels at time $t$. Hence, direction
$\mathbf{d}_{\max}$ is in parallel to the direction of local image gradient,
namely the normal flow’s direction. $\hfill\square$

The lemma demonstrates how to determine the direction of normal flow, but
there is still a unknown scale factor between the magnitude of
$\mathbf{d}_{\max}$ and that of the normal flow.

_Proposition._ _Given the mapping function $\mathrm{\Sigma}_{e}$ at pixel
coordinate $\mathbf{x}$ and its finite neighbourhood, and the first-order
partial derivative ${\nabla}\mathrm{\Sigma}_{e}(\mathbf{x})$, normal flow
$\mathbf{n}(\mathbf{x})$ at pixel $\mathbf{x}$ can be calculated by_

$\mathbf{n}(\mathbf{x})=\frac{{\nabla}\mathrm{\Sigma}_{e}(\mathbf{x})}{\left\|{\nabla}\mathrm{\Sigma}_{e}(\mathbf{x})\right\|^{2}}.$
(38)

_Proof._ Direction $\mathbf{d}_{\max}$ owns the algebra meaning: the maximum
time increase magnitude
$\mathrm{\Delta}{t}_{\max}=\mathbf{d}_{\text{max}}\cdot\nabla\mathrm{\Sigma}_{e}(\mathbf{x})=\left\|{\nabla}\mathrm{\Sigma}_{e}(\mathbf{x})\right\|$
given a unit pixel displacement ($\|\mathbf{d}_{\max}\|=1$). Since the
magnitude of the normal flow is equal to the pixel displacement per unit time,
dividing $\mathbf{d}_{\max}$ by $\mathrm{\Delta}{t}_{\max}$ is the normal
flow, i.e.,
$\mathbf{n}(\mathbf{x})=\frac{\mathbf{d}_{\max}}{\mathrm{\Delta}{t}_{\max}}=\frac{\nabla\mathrm{\Sigma}_{e}(\mathbf{x})}{\|\nabla\mathrm{\Sigma}_{e}(\mathbf{x})\|^{2}}$.
$\hfill\square$

The above normal flow calculation requires the smooth and differentiable
mapping function, we can approximate its partial derivative operator with the
modern image processing techniques such as _Sobel filter_ or _local plane
fitting_. All above-mentioned geometric elements are illustrated in Fig. 5.
Based on this, we present the normal flow constraint
$\mathbf{n}(\mathbf{x})^{\top}\mathbf{u}(\mathbf{x})=\|\mathbf{n}(\mathbf{x})\|^{2}$
for robust estimation of motion and structure parameters in Sec. 4 of our
paper.

Figure 5:  Illustration of event-based normal flow, explained on the spatio-
temporal profile induced by a moving edge. Given a pixel coordinate
$\mathbf{x}$, when the increment $\mathrm{\Delta}\mathbf{x}$’s direction
$\mathbf{d}$ equals $\mathbf{d}_{\max}$ ($\mathbf{d}_{\max}$ is the one that
maximizes the directional derivative
$\nabla_{\mathbf{d}}\mathrm{\Sigma}_{e}(\mathbf{x})=\mathbf{d}\cdot\nabla\mathrm{\Sigma}_{e}(\mathbf{x})$,
i.e.,
$\mathbf{d}_{\max}=\frac{\nabla\mathrm{\Sigma}_{e}(\mathbf{x})}{\|\nabla\mathrm{\Sigma}_{e}(\mathbf{x})\|}$),
it aligns with the direction of the time surface gradient
$\nabla\mathrm{\Sigma}_{e}(\mathbf{x})$, thereby determining the direction of
the normal flow $\mathbf{n}(\mathbf{x})$. Since the time increment (lifetime)
corresponding to $\mathbf{d}_{\max}$ is
$\mathrm{\Delta}{t}_{\max}=\mathbf{d}_{\max}\cdot\nabla\mathrm{\Sigma}_{e}(\mathbf{x})=\|\nabla\mathrm{\Sigma}_{e}(\mathbf{x})\|$,
the normal flow
$\mathbf{n}(\mathbf{x})=\frac{\mathbf{d}_{\max}}{\mathrm{\Delta}{t}_{\max}}=\frac{\nabla\mathrm{\Sigma}_{e}(\mathbf{x})}{\|\nabla\mathrm{\Sigma}_{e}(\mathbf{x})\|^{2}}$.

##  Appendix 0.B More About Differential Homography

In this section, a detailed explanation to Sec. 4 of our paper is added. In
particular, we first discuss the way to recover the true differential
homography from our linear solution (Sec. 0.B.1). Then, we disclose how to
retrieve the motion and scene parameters by decomposing differential
homography (Sec. 0.B.2).

###  0.B.1 Recover the True Differential Homography from Our Linear Solution

As is well-known [27, 49], Eq. (24) can only recover $\mathbf{H}_{L}$ =
$\mathbf{H}_{d}$ +$\epsilon\mathbf{I}$ with an unknown scale $\epsilon$ since
$\mathbf{H}_{d}$ has a one-dimensional null space. So our linear solver faces
this problem too, here we could recover the true differential homography
matrix as follows: after solving for $\mathbf{H}_{L}$, let
$\mathbf{M}_{L}=\mathbf{H}_{L}+\mathbf{H}_{L}^{\top}$, then the second largest
eigenvalue of $\mathbf{M}_{L}$ equals $2\epsilon$, which enabling recovery of
differential homography matrix by
$\mathbf{H}_{d}=\mathbf{H}_{L}-\epsilon\mathbf{I}$.

_Proof._ Let

$\displaystyle\mathbf{M}_{L}$
$\displaystyle=-\left(\frac{\boldsymbol{\nu}}{d}\mathbf{n}^{\top}+\mathbf{n}{\frac{\boldsymbol{\nu}}{d}}^{\top}\right)+2\epsilon\mathbf{I}$
(39) $\displaystyle=-\mathbf{M}+2\epsilon\mathbf{I}.$

We define the eigenvalue of $\mathbf{M}$ as
$\lambda_{\max},\lambda_{\text{mid}},\lambda_{\min}$ with the decreasing
order. Obviously $\mathbf{M}$ is a symmetric matrix, so according to the
Rayleigh quotient theorem, for any vector
$\\{\mathbf{x}|\mathbf{x}\in\mathbb{R}^{3},\left\|\mathbf{x}\right\|=1\\}$ we
have

$\lambda_{\min}\mathbf{x}^{\top}\mathbf{x}\leq\mathbf{x}^{\top}\mathbf{M}\mathbf{x}\leq\lambda_{\max}\mathbf{x}^{\top}\mathbf{x}.$
(40)

For simplicity, we omit $d$ and define a new variable
$\mathbf{v}={\boldsymbol{\nu}}/{d}$. Taking a vector $\mathbf{x}$
perpendicular to ($\mathbf{v}+\mathbf{n}$), and not perpendicular to
$\mathbf{v}$ and $\mathbf{n}$ at the same time, we can prove that

$\displaystyle\lambda_{\min}\mathbf{x}^{\top}\mathbf{x}$
$\displaystyle\leq\mathbf{x}^{\top}\left(\mathbf{n}\mathbf{v}^{\top}+\mathbf{v}\mathbf{n}^{\top}\right)\mathbf{x}$
(41)
$\displaystyle=\mathbf{x}^{\top}((\mathbf{v}+\mathbf{n})(\mathbf{v}+\mathbf{n})^{\top}-(\mathbf{v}\mathbf{v}^{\top}+\mathbf{n}\mathbf{n}^{\top}))\mathbf{x}$
$\displaystyle=-\mathbf{x}^{\top}(\mathbf{v}\mathbf{v}^{\top}+\mathbf{n}\mathbf{n}^{\top})\mathbf{x}$
$\displaystyle=-[(\mathbf{x}^{\top}\mathbf{v})^{2}+(\mathbf{x}^{\top}\mathbf{n})^{2}]<0.$

And we know that the inner product of the vector itself
$\mathbf{x}^{\top}\mathbf{x}$ is always positive. We can then conclude that
$\lambda_{\min}<0$.

Similarly, for another $\mathbf{x}$ perpendicular to
($\mathbf{v}-\mathbf{n}$), and not perpendicular to $\mathbf{v}$ and
$\mathbf{n}$ at the same time we can prove that

$\displaystyle\lambda_{\max}\mathbf{x}^{\top}\mathbf{x}$
$\displaystyle\geq\mathbf{x}^{\top}(\mathbf{n}\mathbf{v}^{\top}+\mathbf{v}\mathbf{n}^{\top})\mathbf{x}$
(42)
$\displaystyle=\mathbf{x}^{\top}((\mathbf{v}\mathbf{v}^{\top}+\mathbf{n}\mathbf{n}^{\top})-(\mathbf{v}-\mathbf{n})(\mathbf{v}-\mathbf{n})^{\top})\mathbf{x}$
$\displaystyle=\mathbf{x}^{\top}(\mathbf{v}\mathbf{v}^{\top}+\mathbf{n}\mathbf{n}^{\top})\mathbf{x}$
$\displaystyle=[(\mathbf{x}^{\top}\mathbf{v})^{2}+(\mathbf{x}^{\top}\mathbf{n})^{2}]>0.$

Similarly, we can prove that $\lambda_{\max}>0$.

For any rank deficient square matrix, its determinant must equal to zero and
there are at least one eigenvalue equals to zero. The rank of the $\mathbf{M}$
is less than two since it contains two outer product of vectors

$\text{rank}(\mathbf{M})=\text{rank}\left(\frac{\boldsymbol{\nu}}{d}\mathbf{n}^{\top}+\mathbf{n}\frac{\boldsymbol{\nu}}{d}^{\top}\right)\leq
2.$ (43)

Besides, $\lambda_{\max}>0$ and $\lambda_{\min}<0$, we can come to the
conclusion that

$\lambda_{\text{mid}}=0.$ (44)

Consequently the second largest value of
$\mathbf{M}_{L}=-\mathbf{M}+2\epsilon\mathbf{I}$ equals to $2\epsilon$.
$\square$

###  0.B.2 Retrieving Motion and Structure Parameters by Decomposing
Differential Homography

After recovering the real differential homography matrix $\mathbf{H}_{d}$,
this section introduces the decomposition method to get the motion and
structure parameters.

Firstly, the eigen-decomposition of $\mathbf{M}$ can be written as:

$\mathbf{M}=\mathbf{Q}\boldsymbol{\Lambda}\mathbf{Q}^{\top}.$ (45)

where
$\mathbf{Q}=\\{\mathbf{q}_{\max},\mathbf{q}_{\text{mid}},\mathbf{q}_{\min}\\}$
is the orthogonal matrix which contains all eigenvectors of $\mathbf{M}$ and
$\text{diag}(\boldsymbol{\Lambda})=\\{\lambda_{\max},\lambda_{\text{mid}},\lambda_{\min}\\}$
which collects all eigenvalues in the diagonal elements.

We then can define two auxiliary vectors $\mathbf{j}$ and $\mathbf{k}$

$\displaystyle\mathbf{j}$
$\displaystyle=\sqrt{\frac{\lambda_{\max}}{2}}\mathbf{q}_{\max}+\sqrt{\frac{-\lambda_{\min}}{2}}\mathbf{q}_{\min},$
(46) $\displaystyle\mathbf{k}$
$\displaystyle=\sqrt{\frac{\lambda_{\max}}{2}}\mathbf{q}_{\max}-\sqrt{\frac{-\lambda_{\min}}{2}}\mathbf{q}_{\min},$

such that

$\mathbf{M}=\mathbf{j}\mathbf{k}^{\top}+\mathbf{k}\mathbf{j}^{\top}.$ (47)

_Proof._

In Sec. 0.B.1, we already concluded that the egienvalues of $\mathbf{M}$ have
following properties:

$\lambda_{\max}>\lambda_{\text{mid}}=0>\lambda_{\min}.$ (48)

Then, we can substitute this result into Eq. (45):

$\displaystyle\mathbf{M}$
$\displaystyle=\begin{bmatrix}\mathbf{q}_{\max}&\mathbf{q}_{\text{mid}}&\mathbf{q}_{\min}\end{bmatrix}\begin{bmatrix}\lambda_{\max}&&\\\
&0&\\\
&&\lambda_{\min}\end{bmatrix}\begin{bmatrix}{\mathbf{q}_{\max}}^{\top}\\\
{\mathbf{q}_{\text{mid}}}^{\top}\\\ {\mathbf{q}_{\min}}^{\top}\end{bmatrix}$
(49)
$\displaystyle=\lambda_{\max}\mathbf{q}_{\max}{\mathbf{q}_{\max}}^{\top}+\lambda_{\min}\mathbf{q}_{\min}{\mathbf{q}_{\min}}^{\top}$
$\displaystyle=\frac{\lambda_{\max}}{2}\mathbf{q}_{\max}{\mathbf{q}_{\max}}^{\top}+\frac{\lambda_{\max}}{2}\mathbf{q}_{\max}{\mathbf{q}_{\max}}^{\top}+$
$\displaystyle~{}~{}~{}~{}\frac{\lambda_{\min}}{2}\mathbf{q}_{\min}{\mathbf{q}_{\min}}^{\top}+\frac{\lambda_{\min}}{2}\mathbf{q}_{\min}{\mathbf{q}_{\min}}^{\top}$
$\displaystyle=(\frac{\lambda_{\max}}{2}\mathbf{q}_{\max}{\mathbf{q}_{\max}}^{\top}+\sqrt{\frac{-\lambda_{\min}\lambda_{\max}}{4}}\mathbf{q}_{\min}{\mathbf{q}_{\max}}^{\top}$
$\displaystyle~{}~{}~{}~{}-\sqrt{\frac{-\lambda_{\min}\lambda_{\max}}{4}}\mathbf{q}_{\max}{\mathbf{q}_{\min}}^{\top}-\frac{-\lambda_{\min}}{2}\mathbf{q}_{\min}{\mathbf{q}_{\min}}^{\top})$
$\displaystyle~{}~{}~{}~{}+(\frac{\lambda_{\max}}{2}\mathbf{q}_{\max}{\mathbf{q}_{\max}}^{\top}-\sqrt{\frac{-\lambda_{\min}\lambda_{\max}}{4}}\mathbf{q}_{\min}{\mathbf{q}_{\max}}^{\top}$
$\displaystyle~{}~{}~{}~{}+\sqrt{\frac{-\lambda_{\min}\lambda_{\max}}{4}}\mathbf{q}_{\max}{\mathbf{q}_{\min}}^{\top}-\frac{-\lambda_{\min}}{2}\mathbf{q}_{\min}{\mathbf{q}_{\min}}^{\top}).$

It is easy to prove

$\displaystyle\mathbf{j}\mathbf{k}^{\top}=$
$\displaystyle\left(\frac{\lambda_{\max}}{2}\mathbf{q}_{\max}{\mathbf{q}_{\max}}^{\top}+\sqrt{\frac{-\lambda_{\min}\lambda_{\max}}{4}}\mathbf{q}_{\min}{\mathbf{q}_{\max}}^{\top}\right.$
(50)
$\displaystyle\quad\left.-\sqrt{\frac{-\lambda_{\min}\lambda_{\max}}{4}}\mathbf{q}_{\max}{\mathbf{q}_{\min}}^{\top}-\frac{-\lambda_{\min}}{2}\mathbf{q}_{\text{min}}{\mathbf{q}_{\min}}^{\top}\right),$
$\displaystyle\mathbf{k}\mathbf{j}^{\top}=$
$\displaystyle\left(\frac{\lambda_{\max}}{2}\mathbf{q}_{\max}{\mathbf{q}_{\max}}^{\top}-\sqrt{\frac{-\lambda_{\min}\lambda_{\max}}{4}}\mathbf{q}_{\min}{\mathbf{q}_{\max}}^{\top}\right.$
$\displaystyle\quad\left.+\sqrt{\frac{-\lambda_{\min}\lambda_{\max}}{4}}\mathbf{q}_{\max}{\mathbf{q}_{\min}}^{\top}-\frac{-\lambda_{\min}}{2}\mathbf{q}_{\min}{\mathbf{q}_{\min}}^{\top}\right).$

We then can conclude that
$\mathbf{M}=\mathbf{j}\mathbf{k}^{\top}+\mathbf{k}\mathbf{j}^{\top}$.
$\hfill\square$

Since $\mathbf{M}$ is a symmetric matrix, we cannot directly distinguish
${\boldsymbol{\nu}}/{d}$ and $\mathbf{n}$ from each other. Note that
$\mathbf{n}$ is a normal vector which has a unique norm. We can thus decompose
all candidates as

$\left\\{\begin{array}[]{ll}\frac{1}{d}\boldsymbol{\nu}=\mathbf{j}\|\mathbf{k}\|,\quad\mathbf{n}=\mathbf{k}/\|\mathbf{k}\|,&\quad\boldsymbol{\omega}=-(\mathbf{H}_{d}+\mathbf{j}\mathbf{k}^{\top})^{\vee}\\\\[4.30554pt]
\frac{1}{d}\boldsymbol{\nu}=\mathbf{k}\|\mathbf{j}\|,\quad\mathbf{n}=\mathbf{j}/\|\mathbf{j}\|,&\quad\boldsymbol{\omega}=-(\mathbf{H}_{d}+\mathbf{k}\mathbf{j}^{\top})^{\vee},\end{array}\right.$
(51)

where the $(.)^{\vee}$ denotes the inverse operator of $(.)_{\times}$, which
converts a skew-symmetric matrix into a vector.

Above all, the decomposition of the differential homography matrix results in
two sets of possible solutions that are indistinguishable in the absence of
prior scene or motion information. However, when integrated with data from
other sensors, this decomposition remains valuable, enabling the estimation of
otherwise unknown motion and scene parameters.

##  Appendix 0.C Dataset and Evaluation Metrics

To supplement Sec. 5 of our paper, we detail our dataset as well as the
evaluation metrics used in the experiments.

###  0.C.1 Our Dataset

####  0.C.1.1 Real Data.

For our study on rotational motion estimation, we initially performed
experiments using the ECD dataset[32], which was captured using an event
camera[3] with a relatively low spatial resolution (a DAVIS240C with
240$\times$180 px). To assess our proposed algorithm with more modern sensors,
we recorded two similar sequences, _ground_rotation_ and _boxed_rotation_ ,
using an event camera (iniVation DAVIS346) with a higher spatial resolution
(346$\times$260 px). Additionally, we employ an Xsens Mti-630 IMU to provide
ground-truth angular velocity. Fig. 6 offers an insight into the recorded
data. Meanwhile, Fig. 7 depicts the ground-truth angular velocity information,
featuring high dynamics of the event camera in each sequence.

####  0.C.1.2 Synthetic Data.

In the paper, we conduct a comprehensive set of experiments using a synthetic
dataset created with a simulator [37]. This dataset comprises sequences,
including _two_wall_translation_ , _patterns_rotation_ , _cubes_rotation_ ,
_patterns_6dof_ , and _cubes_6dof_ , synthesized using an event camera with a
spatial resolution of 640$\times$480 px. Additionally, the simulator provides
us with ground truth data of the optical flow, depth and camera trajectory for
each sequence. All sequences are recorded over a time span of 0.5 seconds
under different motion patterns. Specifically, sequences _patterns_6dof_ and
_cubes_6dof_ feature a combined motion pattern of a linear translation plus a
rotation around a certain axis. In contrast, sequences _patterns_rotation_ and
_cubes_rotation_ involve pure rotation at a constant angular velocity, while
sequence _two_wall_translation_ exhibits a pure translation. Figure 8 shows
the simulated scenes and the pose of the camera (represented by an orange
wireframe pyramid).

(a) _ground_rotation_

(b) _boxes_rotation_

Figure 6: Event data (represented with a naive accumulation) and corresponding
frames from our data collected using a DAVIS-346 event camera.

(a) _ground_rotation_

(b) _boxes_rotation_

Figure 7: Angular velocity measurements from an IMU, used as ground truth.

(a) _cubes_rotation_

(b) _patterns_rotation_

Figure 8: An event camera (orange) capturing different sequences in simulated
scenes.

###  0.C.2 Evaluation Metrics

In the quantitative evaluation, we utilize the average angular velocity error
($e_{w}$) and root mean square error (RMSE) to measure the angular velocity
estimation error. The latter is defined as

${\rm
RMSE}_{\omega}=\sqrt{\frac{1}{3m}\sum_{i=1}^{m}\left(e_{\omega_{x}i}^{2}+e_{\omega_{y}i}^{2}+e_{\omega_{z}i}^{2}\right)}.$
(52)

For linear velocity evaluation, we also evaluate error by means of the average
($e_{v}$) and the RMSE, with

${\rm
RMSE}_{v}=\sqrt{\frac{1}{3m}\sum_{i=1}^{m}\left(e_{v_{x}i}^{2}+e_{v_{y}i}^{2}+e_{v_{z}i}^{2}\right)}.$
(53)

As for the evaluation metric for the differential homography, we assess the
errors in the estimated differential homography matrix by means of the
Frobenius norm of the difference with respect to the ground truth, rather than
directly comparing motion and structure parameters. This is because it is non-
trivial to determine the correct one from the resulting two-set motion and
structure parameters (see 0.B.2), and it is out of the scope of this paper.

##  Appendix 0.D Implementation Details

This section provides implementation details to supplement Sec. 5 of the
paper.

####  0.D.0.1 Normal Flow Extraction.

Since our minimal solver is built on normal flow constraint, the pre-stage
normal flow calculation becomes crucial. Here, we directly calculate gradients
from raw event data by local plane fitting [2, 30]. To further refine these
calculations, RANSAC is applied. Ultimately, we employ (38) to accurately
determine the partial normal flow.

####  0.D.0.2 Linear Solver.

From the event batch, we extracted numerous partial normal flow constraints,
which enabled us to construct an overdetermined system of equations. However,
the estimated normal flow is still subject to noise due to the limitations of
existing methods. To address this, we employed RANSAC for robust sampling.

####  0.D.0.3 Hyper-Parameter Settings.

Parameters including window size and threshold can have a crucial influence on
the performance of each solver, which therefore needs clarification. Here, we
have configured the spatial and temporal window sizes for normal flow
extraction at $7\text{px}\times 7\text{px}\times 0.04\text{s}$, respectively.
Additionally, the RANSAC threshold of plane fitting is set to $10^{-5}$, while
the RANSAC threshold of linear solver is established at $10^{-4}$.

##  Appendix 0.E Additional Results

To further investigate and discuss the performance of our solver, we conduct
some extensive experiments to supplement Sec. 5 of the paper.

(a) Optical Flow

(b) Depth

(c) Rotational Motion

(d) Differential Homography

(e) 6-DoF Motion

Figure 9: Numerical stability analysis on the tasks of optical flow, depth
estimation, rotational motion estimation, differential homography, and 6-DoF
motion tracking, respectively.  Figure 10: Results of the continuous-time
angular velocity estimator on sequence _dynamic_rotation_.

###  0.E.1 Numerical Stability Analysis

We conduct the numerical stability analysis of the linear solver in the five
problems tackled: 1) Optical flow estimation, 2) Depth estimation, 3)
Rotational estimation, 4) Differential homography estimation, and 5) 6-DoF
motion tracking. We gradually increase the noise level from 0.01 px to 100 px
on the normal flow observations, and assess how sensitive the linear solver is
in each problem.

As shown in Fig. 9, for optical flow estimation, the error increases with the
noise level and the linear solver can estimate optical flow under 1-pixel
error in normal flow observation. Also, the linear solver can estimate depth
with high accuracy under the same noise level. For tasks involving rotation,
differential homography, and 6-DoF motion estimation, the linear solver
demonstrates greater resilience, effectively functioning under more severe
noise conditions. In general, accurate results are witnessed in all five
problems under a noise level of 1 pixel, thereby substantiating the robustness
of our linear solver.

###  0.E.2 Continuous-Time Nonlinear Solver

The experiments on sequence _shapes_rotation_ from the dataset [32] have
already demonstrated the superiority of the continuous-time nonlinear solver
in handling aggressive motion. However, the relatively simple texture of this
sequence raises questions about the solver’s performance in more complex
visual environments. To address this, we conduct an additional experiment on
sequence _dynamic_rotation_ , also from the dataset [32], which features a
more natural and higher textured scene. The results, illustrated in Fig. 10,
confirm the solver’s capability to perform accurately in everyday-textured
scenes, indicating its robustness across diverse scene complexities.

## References

  * [1] Barranco, F., Fermüller, C., Aloimonos, Y., Ros, E.: Joint direct estimation of 3D geometry and 3D motion using spatio temporal gradients. Pattern Recog. 113, 107759 (2021). https://doi.org/10.1016/j.patcog.2020.107759
  * [2] Benosman, R., Clercq, C., Lagorce, X., Ieng, S.H., Bartolozzi, C.: Event-based visual flow. IEEE Trans. Neural Netw. Learn. Syst. 25(2), 407–417 (2014). https://doi.org/10.1109/TNNLS.2013.2273537
  * [3] Brandli, C., Berner, R., Yang, M., Liu, S.C., Delbruck, T.: A 240x180 130dB 3$\mu$s latency global shutter spatiotemporal vision sensor. IEEE J. Solid-State Circuits 49(10), 2333–2341 (2014). https://doi.org/10.1109/JSSC.2014.2342715
  * [4] Brodskỳ, T., Fermüller, C., Aloimonos, Y.: Structure from motion: Beyond the epipolar constraint. Int. J. Comput. Vis. 37, 231–258 (2000)
  * [5] Bryner, S., Gallego, G., Rebecq, H., Scaramuzza, D.: Event-based, direct camera tracking from a photometric 3D map using nonlinear optimization. In: IEEE Int. Conf. Robot. Autom. (ICRA) (2019). https://doi.org/10.1109/ICRA.2019.8794255
  * [6] Clady, X., Clercq, C., Ieng, S.H., Houseini, F., Randazzo, M., Natale, L., Bartolozzi, C., Benosman, R.: Asynchronous visual event-based time-to-contact. Front. Neurosci. 8(9) (2014). https://doi.org/10.3389/fnins.2014.00009
  * [7] Conradt, J., Cook, M., Berner, R., Lichtsteiner, P., Douglas, R.J., Delbruck, T.: A pencil balancing robot using a pair of AER dynamic vision sensors. In: IEEE Int. Symp. Circuits Syst. (ISCAS). pp. 781–784 (2009). https://doi.org/10.1109/ISCAS.2009.5117867
  * [8] Corke, P.: Robotics, Vision and Control: Fundamental Algorithms in MATLAB. Springer Tracts in Advanced Robotics, Springer (2017). https://doi.org/10.1007/978-3-319-54413-7
  * [9] De Boor, C.: On calculating with b-splines. Journal of Approximation theory 6(1), 50–62 (1972)
  * [10] Delbruck, T., Lang, M.: Robotic goalie with 3ms reaction time at 4% CPU load using event-based dynamic vision sensor. Front. Neurosci. 7, 223 (2013). https://doi.org/10.3389/fnins.2013.00223
  * [11] Dinaux, R., Wessendorp, N., Dupeyroux, J., De Croon, G.C.: Faith: Fast iterative half-plane focus of expansion estimation using optic flow. IEEE Robotics and Automation Letters 6(4), 7627–7634 (2021)
  * [12] Fermüller, C.: Passive navigation as a pattern recognition problem. Int. J. Comput. Vis. 14(2), 147–158 (1995)
  * [13] Fischler, M.A., Bolles, R.C.: Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography. Commun. ACM 24(6), 381–395 (1981). https://doi.org/10.1145/358669.358692
  * [14] Gallego, G., Gehrig, M., Scaramuzza, D.: Focus is all you need: Loss functions for event-based vision. In: IEEE Conf. Comput. Vis. Pattern Recog. (CVPR). pp. 12272–12281 (2019). https://doi.org/10.1109/CVPR.2019.01256
  * [15] Gallego, G., Rebecq, H., Scaramuzza, D.: A unifying contrast maximization framework for event cameras, with applications to motion, depth, and optical flow estimation. In: IEEE Conf. Comput. Vis. Pattern Recog. (CVPR). pp. 3867–3876 (2018). https://doi.org/10.1109/CVPR.2018.00407
  * [16] Gao, L., Liang, Y., Yang, J., Wu, S., Wang, C., Chen, J., Kneip, L.: Vector: A versatile event-centric benchmark for multi-sensor slam. IEEE Robotics and Automation Letters 7(3), 8217–8224 (2022). https://doi.org/10.1109/LRA.2022.3186770
  * [17] Gao, L., Su, H., Gehrig, D., Cannici, M., Scaramuzza, D., Kneip, L.: A 5-point minimal solver for event camera relative motion estimation. In: Int. Conf. Comput. Vis. (ICCV). pp. 8015–8025 (2023). https://doi.org/10.1109/ICCV51070.2023.00739
  * [18] Gu, C., Learned-Miller, E., Sheldon, D., Gallego, G., Bideau, P.: The spatio-temporal Poisson point process: A simple model for the alignment of event camera data. In: Int. Conf. Comput. Vis. (ICCV). pp. 13495–13504 (2021). https://doi.org/10.1109/ICCV48922.2021.01324
  * [19] Guo, S., Gallego, G.: CMax-SLAM: Event-based rotational-motion bundle adjustment and SLAM system using contrast maximization. IEEE Trans. Robot. pp. 1–20 (2024). https://doi.org/10.1109/TRO.2024.3378443
  * [20] Huang, X., Zhang, Y., Xiong, Z.: Progressive spatio-temporal alignment for efficient event-based motion estimation. In: IEEE Conf. Comput. Vis. Pattern Recog. (CVPR). pp. 1537–1546 (2023)
  * [21] Ji, H., Fermuller, C.: A 3D shape constraint on video. IEEE Trans. Pattern Anal. Mach. Intell. 28(6), 1018–1023 (2006)
  * [22] Kim, H., Leutenegger, S., Davison, A.J.: Real-time 3D reconstruction and 6-DoF tracking with an event camera. In: Eur. Conf. Comput. Vis. (ECCV). pp. 349–364 (2016). https://doi.org/10.1007/978-3-319-46466-4_21
  * [23] Lagorce, X., Orchard, G., Gallupi, F., Shi, B.E., Benosman, R.: HOTS: A hierarchy of event-based time-surfaces for pattern recognition. IEEE Trans. Pattern Anal. Mach. Intell. 39(7), 1346–1359 (Jul 2017). https://doi.org/10.1109/TPAMI.2016.2574707
  * [24] Liu, D., Parra, A., Chin, T.J.: Spatiotemporal registration for event-based visual odometry. In: IEEE Conf. Comput. Vis. Pattern Recog. (CVPR). pp. 4935–4944 (2021). https://doi.org/10.1109/CVPR46437.2021.00490
  * [25] Longuet-Higgins, H.C., Prazdny, K.: The interpretation of a moving retinal image. Proceedings of the Royal Society of London. Series B. Biological Sciences 208(1173), 385–397 (1980). https://doi.org/10.1098/rspb.1980.0057
  * [26] Lu, X., Zhou, Y., Niu, J., Zhong, S., Shen, S.: Event-based visual inertial velometer (2024)
  * [27] Ma, Y., Soatto, S., Košecká, J., Sastry, S.: An Invitation to 3-D Vision: From Images to Geometric Models. Springer (2004)
  * [28] Meyer, F., Bouthemy, P.: Estimation of time-to-collision maps from first order motion models and normal flows. In: 1992 11th IAPR International Conference on Pattern Recognition. vol. 1, pp. 78–82. IEEE Computer Society (1992)
  * [29] Meyer, F.G.: Time-to-collision from first-order models of the motion field. IEEE transactions on robotics and automation 10(6), 792–798 (1994)
  * [30] Mueggler, E., Forster, C., Baumli, N., Gallego, G., Scaramuzza, D.: Lifetime estimation of events from dynamic vision sensors. In: IEEE Int. Conf. Robot. Autom. (ICRA). pp. 4874–4881 (2015). https://doi.org/10.1109/ICRA.2015.7139876
  * [31] Mueggler, E., Gallego, G., Rebecq, H., Scaramuzza, D.: Continuous-time visual-inertial odometry for event cameras. IEEE Trans. Robot. 34(6), 1425–1440 (Dec 2018). https://doi.org/10.1109/tro.2018.2858287
  * [32] Mueggler, E., Rebecq, H., Gallego, G., Delbruck, T., Scaramuzza, D.: The event-camera dataset and simulator: Event-based data for pose estimation, visual odometry, and SLAM. Int. J. Robot. Research 36(2), 142–149 (2017). https://doi.org/10.1177/0278364917691115
  * [33] Nunes, U.M., Demiris, Y.: Entropy minimisation framework for event-based vision model estimation. In: Eur. Conf. Comput. Vis. (ECCV). pp. 161–176 (2020). https://doi.org/10.1007/978-3-030-58558-7_10
  * [34] Parameshwara, C.M., Hari, G., Fermüller, C., Sanket, N.J., Aloimonos, Y.: Diffposenet: direct differentiable camera pose estimation. In: IEEE Conf. Comput. Vis. Pattern Recog. (CVPR). pp. 6845–6854 (2022)
  * [35] Pijnacker Hordijk, B.J., Scheper, K.Y., De Croon, G.C.: Vertical landing for micro air vehicles using event-based optical flow. J. Field Robot. 35(1), 69–90 (2018)
  * [36] Rebecq, H., Gallego, G., Scaramuzza, D.: EMVS: Event-based multi-view stereo. In: British Mach. Vis. Conf. (BMVC) (2016). https://doi.org/10.5244/C.30.63
  * [37] Rebecq, H., Gehrig, D., Scaramuzza, D.: ESIM: an open event camera simulator. In: Conf. on Robotics Learning (CoRL) (2018)
  * [38] Rebecq, H., Horstschaefer, T., Scaramuzza, D.: Real-time visual-inertial odometry for event cameras using keyframe-based nonlinear optimization. In: British Mach. Vis. Conf. (BMVC) (2017)
  * [39] Rebecq, H., Horstschäfer, T., Gallego, G., Scaramuzza, D.: EVO: A geometric approach to event-based 6-DOF parallel tracking and mapping in real-time. IEEE Robot. Autom. Lett. 2(2), 593–600 (2017). https://doi.org/10.1109/LRA.2016.2645143
  * [40] Shiba, S., Aoki, Y., Gallego, G.: Event collapse in contrast maximization frameworks. Sensors 22(14), 1–20 (2022). https://doi.org/10.3390/s22145190
  * [41] Shiba, S., Aoki, Y., Gallego, G.: A fast geometric regularizer to mitigate event collapse in the contrast maximization framework. Adv. Intell. Syst. p. 2200251 (2022). https://doi.org/10.1002/aisy.202200251
  * [42] Shiba, S., Aoki, Y., Gallego, G.: Secrets of event-based optical flow. In: Eur. Conf. Comput. Vis. (ECCV). pp. 628–645 (2022). https://doi.org/10.1007/978-3-031-19797-0_36
  * [43] Stoffregen, T., Gallego, G., Drummond, T., Kleeman, L., Scaramuzza, D.: Event-based motion segmentation by motion compensation. In: Int. Conf. Comput. Vis. (ICCV). pp. 7243–7252 (2019). https://doi.org/10.1109/ICCV.2019.00734
  * [44] Trucco, E., Verri, A.: Introductory Techniques for 3-D Computer Vision. Prentice Hall PTR, Upper Saddle River, NJ, USA (1998)
  * [45] Zhou, Y., Gallego, G., Rebecq, H., Kneip, L., Li, H., Scaramuzza, D.: Semi-dense 3D reconstruction with a stereo event camera. In: Eur. Conf. Comput. Vis. (ECCV). pp. 242–258 (2018). https://doi.org/10.1007/978-3-030-01246-5_15
  * [46] Zhou, Y., Gallego, G., Shen, S.: Event-based stereo visual odometry. IEEE Trans. Robot. 37(5), 1433–1450 (2021). https://doi.org/10.1109/TRO.2021.3062252
  * [47] Zhu, A.Z., Chen, Y., Daniilidis, K.: Realtime time synchronized event-based stereo. In: Eur. Conf. Comput. Vis. (ECCV). pp. 438–452 (2018). https://doi.org/10.1007/978-3-030-01231-1_27
  * [48] Zhuang, B., Cheong, L.F., Lee, G.H.: Rolling-shutter-aware differential sfm and image rectification. In: Int. Conf. Comput. Vis. (ICCV). pp. 948–956 (2017). https://doi.org/10.1109/ICCV.2017.108
  * [49] Zhuang, B., Tran, Q.H.: Image stitching and rectification for hand-held cameras. In: Eur. Conf. Comput. Vis. (ECCV) (2020). https://doi.org/10.1007/978-3-030-58571-6_15