Jae-Hyeon Park, and Dong Eui Chang

Dong Eui Chang, School of Electrical Engineering, Korea Advanced Institute of
Science and Technology, 291 Daehak-ro, Yuseong-gu, Daejeon 34141, Republic of
Korea.

# Unscented Kalman filter with stable embedding for simple, accurate and
computationally efficient state estimation of systems on manifolds in
Euclidean space

Jae-Hyeon Park  Dong Eui Chang  School of Electrical Engineering, Korea
Advanced Institute of Science and Technology, Daejeon, Korea
<EMAIL_ADDRESS>J. Park  D. E. Chang

(26 April 2016; 6 June 2016; 6 June 2016)

###### Abstract

[Summary]This paper proposes a simple, accurate and computationally efficient
method to apply the ordinary unscented Kalman filter developed in Euclidean
space to systems whose dynamics evolve on manifolds. We use the mathematical
theory called stable embedding to make a variant of unscented Kalman filter
that keeps state estimates in close proximity to the manifold while exhibiting
excellent estimation performance. We confirm the performance of our devised
filter by applying it to the satellite system model and comparing the
performance with other unscented Kalman filters devised specifically for
systems on manifolds. Our devised filter has a low estimation error, keeps the
state estimates in close proximity to the manifold as expected, and consumes a
minor amount of computation time. Also our devised filter is simple and easy
to use because our filter directly employs the off-the-shelf standard
unscented Kalman filter devised in Euclidean space without any particular
manifold-structure-preserving discretization method or coordinate
transformation.

, (2021), Unscented Kalman filter with stable embedding for simple, accurate
and computationally efficient state estimation of systems on manifolds in
Euclidean space, Int J Robust Nonlinear Control., 2021.

###### keywords:

stable embedding, unscented Kalman filter, attitude estimation, satellite,
manifold

††articletype: Article Type

##  1 Introduction

Filtering and state estimation is an important research field that is
essential in various areas of industry such as process control, tracking,
navigation, and robot control. Recently, researchers are designing filtering
and state estimation techniques specifically devised for systems on manifolds.
This is because many systems are inherently defined on manifolds, and applying
filters devised in Euclidean space directly to the systems on manifolds may
cause inaccuracy and divergence in state estimation. Examples of systems on
manifolds are quadrotors,1 winged unmanned aerial vehicles,2 satellites,3
robot arms,4 and arrival tracking in antenna array processing.5 One of the
major drawbacks of working with systems on manifolds is that it does not well
accommodate many tools available for systems in Euclidean space. There have
been efforts to generalize the unscented Kalman filter (UKF) to the filtering
and state estimation of systems on manifolds because it is one of the popular
filtering techniques for nonlinear systems.

The UKF was first devised for accurate filtering and state estimation of
nonlinear systems on Euclidean space. It became one of the major estimation
techniques because it can approximate the statistics of state estimates up to
the second order.6 The detailed explanation of UKF devised in Euclidean space
can be found in the paper by Wan and Van der Merwe7 and we will call it the
standard UKF to distinguish it from the variants of UKF later introduced. The
direct application of the standard UKF to systems on manifolds may not give
accurate estimates, because the standard UKF is derived in Euclidean space, so
state estimates deviate from the manifolds during UKF process. There are
variants of UKF that have been developed for filtering and state estimation of
systems on manifolds. Sipos 8 devised an UKF for systems on Lie group, which
is a special family of manifolds with group structure. This UKF used the
property of Lie algebra which is a tangent space at the identity of Lie group
and defined sigma points of unscented transform on the Lie algebra. Hauberg et
al. 9 devised an UKF for systems on Riemannian manifold, which is a family of
differentiable manifolds with a Riemannian metric defined. This UKF defined
sigma points of unscented transform on the tangent space of the manifold and
employed geodesics and parallel transformations. These UKFs are practical
since many systems of robotics are defined on sets that are both Lie group and
Riemannian manifold. i.e. special orthogonal group, special Euclidean group,
and n-sphere. Brossard et al. 10 used a similar framework of UKF on manifolds
as Hauberg et al. and made an open-source code that can be used for
practitioners not being familiar with manifolds and Lie groups. Cantelobre et
al. 11 applied the method of Brossard et al. on position estimation of
autonomous underwater vehicles. Sjøberg and Egeland 12 developed an UKF for
matrix Lie groups. These variants of UKF require special discretized form of
systems and manifold-structure-preserving computation steps which can be
computationally expensive.

In this article, we propose a method to directly apply the standard UKF
devised in Euclidean space to systems on manifolds. We use the mathematical
theory called stable embedding that was first founded in the paper by Chang et
al.13 Stable embedding was previously applied to extended Kalman filter (EKF)
of systems defined on special orthogonal group and the devised filter was
computationally efficient. 14 We extend systems on manifolds to an ambient
Euclidean space in such a way that the manifold becomes an invariant attractor
or an asymptotically stable invariant set for the extended system in the
Euclidean space. Then we discretize the modified system using the ordinary
Euler discretization method and apply the standard UKF derived in Euclidean
space to the resultant discrete-time system. The reason why we use the Euler
method is that it is the simplest one, but one could also use any off-the-
shelf discretization method such as the Runge-Kutta 4th-order method which
would only increase the accuracy of our proposed UKF. By inheriting the
stability property of the continuous-time extended system before the
discretization, this UKF has high estimation accuracy thanks to state
estimates staying in close proximity of the manifold. We use the satellite
attitude system model 15 to show the excellent performance of our devised UKF.
The satellite system model has its attitude quaternion state defined on
3-sphere
$\operatorname{S}^{3}=\\{\mathbf{q}\in\mathbb{H}\mid\left\|\mathbf{q}\right\|=1\\}$,
where $\mathbb{H}$ is the set of quaternions which is isomorphic to
$\mathbb{R}^{4}$ as a vector space. We have chosen the system defined on
$\operatorname{S}^{3}$, because the state space $\operatorname{S}^{3}$ is both
a Lie group and a Riemannian manifold. We compare our UKF with other UKFs
designed for systems on manifolds. There are also Kalman filters available
that apply to continuous-time systems directly.16, 17 They require integration
of differential momentum equations which can be computationally complex.
Hence, applying Kalman filters to discrete form of continuous-time systems is
frequently used in real applications due to fast computation. Since, one of
the merits of our proposed UKF is fast computation, we choose to discretize
continuous-time systems and then apply the UKFs to them.

This paper is organized as follows. Section 2 briefly explains the existing
methods of UKF for systems on manifolds. Section 3 describes our method of UKF
for systems on manifolds. Section 4 shows the performance of our method with a
simulation on satellite attitude estimation. The conclusion is presented in
Section 5. An appendix is provided at the end of the paper to give a brief
review of the standard UKF and a guideline for choosing a value of the stable
embedding constant later introduced.

##  2 Review of existing methods of UKF for systems on manifolds

Before we explain our method, we briefly introduce existing methods of
applying UKF to systems on manifolds. Every function and manifold is assumed
to be smooth in this paper unless stated otherwise. Let $M$ be a manifold with
dimension $n_{M}$ in Euclidean space $\mathbb{R}^{n_{x}}$ where $n_{M}<n_{x}$.
A general continuous-time system on $M$ is given by

$\displaystyle\dot{x}$ $\displaystyle=f(x,u),\quad x\in
M\subseteq\mathbb{R}^{n_{x}},u\in\mathbb{R}^{n_{u}}$ (1a) $\displaystyle y$
$\displaystyle=g(x,u),\quad x\in
M\subseteq\mathbb{R}^{n_{x}},u\in\mathbb{R}^{n_{u}},$ (1b)

where $f:M\times\mathbb{R}^{n_{u}}\rightarrow TM$ is the state differential
equation, $x$ is the state, $u$ is the input,
$g:M\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{R}^{n_{y}}$ is the measurement
equation, $y\in\mathbb{R}^{n_{y}}$ is the measurement. Here, $TM$ denotes the
tangent bundle of the manifold $M$.

Suppose that $f$ can be extended to an open neighborhood of the manifold $M$
in $\mathbb{R}^{n_{x}}$. If we were to naively apply the standard UKF devised
in Euclidean space, we would discretize the system (1) using the ordinary
Euler discretization without considering the manifold $M$, and apply the
standard UKF. In Appendix A, we give a brief explanation of the standard UKF
described in the paper by Wan and Van der Merwe 7 and the typical values of
parameters used in the standard UKF. Since the manifold $M$ is not closed on
Euclidean addition and subtraction, there arises the problem that $M$ may not
be stable for the discretized system, i.e, discrete-time trajectories of the
discretized system may diverge off the manifold $M$. Some researchers use the
forced projection method to compensate for the state estimates deviating from
the manifold. This method projects the state estimates $\hat{x}_{k|k}$ at the
end of each UKF update to the manifold $M$ by solving the following
optimization problem $\hat{x}_{new,k|k}=\arg\max_{x\in
M}\left\|x-\hat{x}_{k|k}\right\|^{2}$. This least-squares optimization problem
can be solved by using an appropriate gradient descent algorithm regarding the
manifold, but it can be computationally expensive. Detailed explanation of
various optimization algorithms can be found in the book by Chong and Zak.18

Some researchers transform the coordinates of state estimate points on
manifolds to vector spaces, apply the usual UKF correction steps on the vector
spaces and transform the coordinates of points in vector spaces back to the
manifolds to get the next state estimates, because many tools of Euclidean
space are available in vector spaces. The UKF from the paper by Sipos8 is
devised specifically for the case when the manifold $M$ is a Lie group. First,
discretize the continuous-time system (1) defined on the manifold $M$ into a
specific form so that the states of discretized system stay on $M$ as follows:

$x_{k+1}=f_{D}(x_{k},u_{k}),\;x_{k}\in M\;\forall k\in\mathbb{N}_{\geq 0}.$
(2)

Then apply the variant of UKF that defines sigma points on the Lie algebra of
the Lie group, where it is noted that the Lie algebra is a vector space. This
UKF needs computation of several coordinate transformations between the Lie
group and its Lie algebra at each step of estimation, which is computationally
intensive and expensive; refer to the book by Hall19 on Lie groups and Lie
algebras. The UKF from the paper by Hauberg et al.9 is devised for Riemannian
manifolds; refer to the book by Carmo20 on Riemannian manifolds. They first
discretize the continuous-time system (1) defined on the manifold $M$ into a
specific form so that the states of the discretized system stay on $M$ just
like the system (2). Then apply the variant of UKF that defines sigma points
on the tangent spaces at the state estimates on the manifold, where it is
noted that the tangent spaces of a manifold are vector spaces. This UKF also
needs computation of several coordinate transformations between the manifold
and its tangent spaces at each step of estimation, which is computationally
intensive and expensive. In summary, the above two methods proposed in the
paper by Sipos8 and the paper by Hauberg et al.9 both need coordinate changes
at each step of estimation which can be computationally expensive.

##  3 The UKF with stable embedding

Here, we devise a variant of UKF that can apply the standard UKF derived in
Euclidean space to systems on manifolds without any additional manifold-
structure-preserving discretization or coordinate transformations. Our method
is simple and computationally efficient as good as the standard UKF. We use
the mathematical theory called stable embedding in the paper by Chang et al.13
First, extend the system (1) from $M$ to $\mathbb{R}^{n_{x}}$ as follows:

$\displaystyle\dot{x}$ $\displaystyle=f_{E}(x,u),\quad
x\in\mathbb{R}^{n_{x}},u\in\mathbb{R}^{n_{u}}$ (3a) $\displaystyle y$
$\displaystyle=g_{E}(x,u),\quad
x\in\mathbb{R}^{n_{x}},u\in\mathbb{R}^{n_{u}},$ (3b)

such that

$\displaystyle f_{E}(x,u)$ $\displaystyle=f(x,u),\quad\forall\,x\in
M,u\in\mathbb{R}^{n_{u}}$ (4a) $\displaystyle g_{E}(x,u)$
$\displaystyle=g(x,u),\quad\forall x\in M,u\in\mathbb{R}^{n_{u}},$ (4b)

where
$f_{E}:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{R}^{n_{x}}$
and
$g_{E}:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{R}^{n_{y}}$.
Condition (4) implies that the extended system (3) is identical to the
original system (1) when $x\in M$. Thus, $M$ is an invariant set of the
extended system (3). Then, we modify the system (3) to make the manifold $M$
asymptotically stable, i.e. all flows of the system starting from a fixed
neighborhood of $M$ asymptotically converge to $M$ as time tends to infinity.
Suppose there is a function $V:\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}_{\geq
0}$ such that

$M=V^{-1}(0),$ (5)

and

$\langle\nabla
V(x),\,f_{E}(x,u)\rangle=0,\forall\,x\in\mathbb{R}^{n_{x}},u\in\mathbb{R}^{n_{u}},$
(6)

where $V^{-1}(0)=\left\\{x\in\mathbb{R}^{n_{x}}|V(x)=0\right\\}$, and
$\langle\,,\rangle$ is the Euclidean inner product. Consider the following
modified system

$\displaystyle\dot{x}$ $\displaystyle=f_{E}(x,u)-\alpha\nabla V(x),\quad
x\in\mathbb{R}^{n_{x}},u\in\mathbb{R}^{n_{u}}$ (7a) $\displaystyle y$
$\displaystyle=g_{E}(x,u),\quad
x\in\mathbb{R}^{n_{x}},u\in\mathbb{R}^{n_{u}},$ (7b)

where $\alpha>0$ is a constant parameter called stable embedding constant.
Since the function $V$ attains its minimum value of $0$ on $M$, $\nabla
V(x)=0$ for all $x\in M$. So the system (7) coincides with the original system
(1) on $M$ and $M$ is an invariant set of (7). Since $-\nabla V$ points to the
direction of decreasing $V$, it can be shown that $M$ is asymptotically stable
for the system (7) on certain conditions. More concretely speaking, along each
trajectory of the system (7),

$\frac{dV}{dt}=\langle\nabla V,f_{E}(x,u)-\alpha\nabla
V\rangle=-\alpha\|\nabla V\|^{2}\leq 0$

by (6). Since $V$ attains its global minimum value $0$, $\nabla V(x)=0$ in a
neighborhood of $M$ if and only if $x\in M$. Hence, by Lyapunov’s stability
theory, the manifold $M$ is asymptomatically stable for the system (7). Also,
the above inequality shows that the value of parameter $\alpha$ determines the
rate of convergence to $M$; refer to Theorem 1 in the paper by Chang21 for
more on the stability conditions for $M$. An easy way to construct the
function $V$ is to set $V(x)=\|\ell(x)-c\|^{2}$ when the manifold $M$ is
defined by a level set $\ell(x)=c$ of some function $\ell(x)$.

We now add additive noises on the system (7) to produce a stochastic
continuous-time system as follows:

$\displaystyle\dot{x}$ $\displaystyle=f_{E}(x,u)-\alpha\nabla V(x)+w_{C},\quad
x\in\mathbb{R}^{n_{x}},u\in\mathbb{R}^{n_{u}},w_{C}\in\mathbb{R}^{n_{x}}$ (8a)
$\displaystyle y$ $\displaystyle=g_{E}(x,u)+v_{C},\quad
x\in\mathbb{R}^{n_{x}},u\in\mathbb{R}^{n_{u}},v_{C}\in\mathbb{R}^{n_{y}},$
(8b)

where $x(0)\sim(x_{0},P_{0})$ is the initial condition of $x$,
$w_{C}\sim(0,P_{w,C})$ is the process noise, and $v_{C}\sim(0,P_{v,C})$ is the
measurement noise. We can apply additive noises because the system (7) is
globally defined in the Euclidean space. Here, $(\mu,P)$ stands for a white
noise with mean $\mu$ and covariance $P$. It is assumed that $x(0)$, $w_{C}$,
and $v_{C}$ are mutually uncorrelated white noises.

We discretize the system (8) using the ordinary Euler discretization with time
step $h>0$ to apply the standard UKF derived in Euclidean space as follows:

$\displaystyle x_{k+1}$ $\displaystyle=x_{k}+h(f_{E}(x_{k},u_{k})-\alpha\nabla
V(x_{k}))+w_{D},$ (9a) $\displaystyle y_{k}$
$\displaystyle=g_{E}(x_{k},u_{k})+v_{D},$ (9b)

where $x_{k}=x(kh)$, $u_{k}=u(kh)$, and $y_{k}=y(kh)$ for
$k\in\mathbb{N}_{\geq 0}$. The stochastic properties are as follows:
$x_{0}\sim(x_{0},P_{0})$, $w_{D}\sim(0,P_{w,D})$ where $P_{w,D}\approx
hP_{w,C}$, and $v_{D}\sim(0,P_{v,D})$ where $P_{v,D}\approx P_{v,C}/h$. A
detailed explanation of discretization of a stochastic continuous-time system
can be found in the book by Lewis et al.22 p.84 Then we apply the standard UKF
derived in Euclidean space, described in Appendix A, to the system (9)
directly, which is valid since the system (9) is globally defined in Euclidean
space.

We call the above procedure as the UKF with stable embedding (UKF-SE). The big
merit of this method is that all trajectories of the discretized form of the
continuous-time system (7) starting from some neighborhood of $M$ stay in
close proximity of $M$.13 In the following section, we apply UKF-SE to the
satellite system model.

##  4 Simulation

In this section, we conduct numerical simulations with MATLAB to compare
performance of UKF-SE with other UKFs. The computation platform used for
conducting these numerical studies is: CPU - Intel(R) Core(TM) i7-8700 @
3.20GH, memory - 16GB RAM, and software - MATLAB R2018b. We use the satellite
system model from the paper by Kristiansen and Nicklasson15 and the paper by
Hajiyev and Soken23 for applying various UKFs. Most of the constants and
satellite parameters used in the system are from the paper by Hajiyev and
Soken.23 We organize the parameter values used in our simulation in Table 1.
We denote the Earth-centered inertial reference frame by $I$, the orbit
reference frame by $O$, and the body reference frame by $B$ as in the paper by
Kristiansen and Nicklasson.15 The stochastic continuous-time satellite system
model is given as follows:

$\displaystyle\dot{\mathbf{q}}_{BO}$
$\displaystyle=\frac{1}{2}\Lambda(\mathbf{q}_{BO})[0;\Omega_{BO}],$ (10a)
$\displaystyle\dot{\Omega}_{BI}$
$\displaystyle=\mathbb{I}^{-1}((\mathbb{I}\Omega_{BI})\times\Omega_{BI})+\mathbb{I}^{-1}\tau_{g}+w_{\Omega,C},$
(10b)

where $\mathbf{q}_{BO}\in\operatorname{S}^{3}$ is the satellite attitude
quaternion with respect to the orbit frame; $\Omega_{BI}\in\mathbb{R}^{3}$ is
the satellite angular velocity with respect to the inertial frame;
$\Lambda(\mathbf{q}_{BO})\in\mathbb{R}^{4\times 4}$ with
$\mathbf{q}_{BO}=(q_{0},q_{1},q_{2},q_{3})$ is given as

$\displaystyle\Lambda(\mathbf{q}_{BO})=\begin{bmatrix}q_{0}&-q_{1}&-q_{2}&-q_{3}\\\
q_{1}&q_{0}&-q_{3}&q_{2}\\\ q_{2}&q_{3}&q_{0}&-q_{1}\\\
q_{3}&-q_{2}&q_{1}&q_{0}\end{bmatrix};$

$\Omega_{BO}\in\mathbb{R}^{3}$ is the satellite angular velocity with respect
to the orbit frame given as

$\Omega_{BO}=\Omega_{BI}+A_{BO}\begin{bmatrix}0&-\omega_{0}&0\end{bmatrix}^{T};$

$\tau_{g}$ is the gravity gradient torque given as

$\tau_{g}=3\omega_{0}^{2}A_{BO,3}\times(\mathbb{I}A_{BO,3});$

$\mathbb{I}$ is the moment of the inertia matrix of the satellite;
$w_{\Omega,C}\in\mathbb{R}^{3}$ is the process noise which follows the
Gaussian distribution $N(0,\sigma^{2}_{\Omega,C}\cdot I_{3})$ with mean 0 and
covariance $\sigma^{2}_{\Omega,C}\cdot I_{3}$; $A_{BO}\in\mathbb{R}^{3\times
3}$ is the rotation matrix of $\mathbf{q}_{BO}=(q_{0},q_{1},q_{2},q_{3})$
given as

$\displaystyle
A_{BO}=\begin{bmatrix}(q_{0}^{2}+q_{1}^{2}-q_{2}^{2}-q_{3}^{2})&2(q_{1}q_{2}-q_{0}q_{3})&2(q_{1}q_{3}+q_{0}q_{2})\\\
2(q_{1}q_{2}+q_{0}q_{3})&(q_{0}^{2}-q_{1}^{2}+q_{2}^{2}-q_{3}^{2})&2(q_{2}q_{3}-q_{0}q_{1})\\\
2(q_{1}q_{3}-q_{0}q_{2})&2(q_{2}q_{3}+q_{0}q_{1})&(q_{0}^{2}-q_{1}^{2}-q_{2}^{2}+q_{3}^{2})\end{bmatrix};$

$\omega_{0}\approx(\mu/r^{3}_{0})^{1/2}$ is the angular velocity of the orbit
with respect to the inertial frame which is a constant; $\mu$ is the Earth’s
gravitational coefficient; $r_{0}$ is the distance between the satellite and
the Earth; and $A_{BO,3}(t)$ is the 3rd column of $A_{BO}$. The initial state
for the satellite system is $(\mathbf{q}_{BO}(0),\Omega_{BI}(0))\sim
N((\overline{\mathbf{q}}_{BO,0},\overline{\Omega}_{BI,0}),P_{0})$ where
$\overline{\mathbf{q}}_{BO,0}=(\cos(\pi/4),\sin(\pi/4),0,0)$,
$\overline{\Omega}_{BI,0}=(-0.01,0.01,0.01)[$\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}$]$,
and $P_{0}=\operatorname{diagonal}(10^{-8}\cdot I_{4},\>0.01^{2}\cdot I_{3})$.

The satellite measurement sensors are composed of a magnetometer and a rate
gyro. The magnetometer measurement model is given as follows:

$\displaystyle\begin{bmatrix}H_{x}\\\ H_{y}\\\
H_{z}\end{bmatrix}=A_{BO}\begin{bmatrix}H_{1}\\\ H_{2}\\\
H_{3}\end{bmatrix}+v_{mag,C},$ (11)

where $H_{x}$, $H_{y}$, and $H_{z}$ are the measured Earth’s magnetic field
vector components in the body frame. Here, $H_{1}$, $H_{2}$, and $H_{3}$ are
the Earth’s magnetic field vector components in the orbit frame as a function
of time given as

$\displaystyle H_{1}$
$\displaystyle=\frac{M_{e}}{r^{3}_{0}}\left\\{\cos(\omega_{0}t)[\cos(\epsilon)\sin(i)-\sin(\epsilon)\cos(i)\cos(\omega_{e}t)]-\sin(\omega_{0}t)\sin(\epsilon)\sin(\omega_{e}t)\right\\},$
$\displaystyle H_{2}$
$\displaystyle=-\frac{M_{e}}{r^{3}_{0}}[\cos(\epsilon)\cos(i)+\sin(\epsilon)\sin(i)\cos(\omega_{e}t)],$
$\displaystyle H_{3}$
$\displaystyle=\frac{2M_{e}}{r^{3}_{0}}\left\\{\sin(\omega_{0}t)[\cos(\epsilon)\sin(i)-\sin(\epsilon)\cos(i)\cos(\omega_{e}t)]-2\sin(\omega_{0}t)\sin(\epsilon)\sin(\omega_{e}t)\right\\},$

where $M_{e}$ is the magnetic dipole moment of the Earth, $\epsilon$ is the
magnetic dipole tilt, $i$ is the orbit inclination, and $\omega_{e}$ is the
spin rate of the Earth. In (11), $v_{mag,C}\in\mathbb{R}^{3}$ is the
magnetometer measurement noise which follows $N(0,\sigma^{2}_{mag,C}\cdot
I_{3})$. The rate gyros measurement model is given as follows:

$\displaystyle\Omega_{BI,meas}=\Omega_{BI}+v_{rate,C},$ (12)

where $\Omega_{BI,meas}$ is the measured angular rate of the satellite, and
$v_{rate,C}\in\mathbb{R}^{3}$ is the rate gyros measurement noise which
follows $N(0,\sigma^{2}_{rate,C}\cdot I_{3})$. All the measurements are
measured with a sampling frequency $f$.

Table 1: Parameter values used in our simulation  Notation | Parameter | Value
---|---|---
$\mathbb{I}$ | the moment of the inertia matrix of the satellite | $\operatorname{diag}(2.1,2.0,1.9)\times 10^{-3}\,[$\mathrm{k}\mathrm{g}\,\mathrm{m}^{2}$]$
$\sigma_{\Omega,C}$ | the standard deviation of the process noise | $0.01\,[$\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}$]$
$\mu$ | the Earth’s gravitational coefficient | $3.98601\times 10^{14}\,[$\mathrm{m}^{3}\mathrm{/}\mathrm{s}^{2}$]$
$r_{0}$ | the distance between the satellite and the Earth | $6928140\,[$\mathrm{m}$]$
$M_{e}$ | the magnetic dipole moment of the Earth | $7.943\times 10^{15}\,[$\mathrm{Wb}\text{\,}\mathrm{m}$]$
$\epsilon$ | the magnetic dipole tilt | $11.7$\mathrm{\SIUnitSymbolDegree}$$
$i$ | the orbit inclination | $97$\mathrm{\SIUnitSymbolDegree}$$
$\omega_{e}$ | the spin rate of the Earth | $7.29\times 10^{-5}\,[$\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}$]$
$\sigma_{mag,C}$ | the standard deviation of the magnetometer measurement noise | $200\,[$\mathrm{nT}$]$
$\sigma_{rate,C}$ | the standard deviation of the rate gyros measurement noise | $1.84\times 10^{-4}\,[$\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}$]$
$f$ | the sampling frequency of the measurement | $10\,[$\mathrm{Hz}$]$

The satellite system (10) has a state $(\mathbf{q}_{BO},\Omega_{BI})$ in
$\operatorname{S}^{3}\times\mathbb{R}^{3}$ which is 6 dimensional manifold in
$\mathbb{H}\times\mathbb{R}^{3}$. We do not consider $\Omega_{BO}$ as a state
because it does not have a differential equation. We first extend the system
(10) from $\operatorname{S}^{3}\times\mathbb{R}^{3}$ to
$\mathbb{H}\times\mathbb{R}^{3}$. Define a function
$V(\mathbf{q}_{BO},\Omega_{BI}):\mathbb{H}\times\mathbb{R}^{3}\rightarrow\mathbb{R}_{\geq
0}$ as

$V(\mathbf{q}_{BO},\Omega_{BI})=\frac{1}{4}(\left\|\mathbf{q}_{BO}\right\|^{2}-1)^{2}.$

It is easy to see that $\operatorname{S}^{3}\times\mathbb{R}^{3}=V^{-1}(0)$.
The derivative of $V$ is calculated as
$\nabla_{\mathbf{q}_{BO}}{V}=(\left\|\mathbf{q}_{BO}\right\|^{2}-1)\mathbf{q}_{BO}$,
and $\nabla_{\Omega_{BI}}{V}=0$. This $V$ satistifes the conditions (5) and
(6) for $V$ in Section 3. We add $-\alpha\nabla{V}$ where $\alpha>0$ to the
extended system to get the final modified system as follows:

$\displaystyle\dot{\mathbf{q}}_{BO}$
$\displaystyle=\frac{1}{2}\Lambda(\mathbf{q}_{BO})[0;\Omega_{BO}]-\alpha(\left\|\mathbf{q}_{BO}\right\|^{2}-1)\mathbf{q}_{BO},$
(13a) $\displaystyle\dot{\Omega}_{BI}$
$\displaystyle=\mathbb{I}^{-1}((\mathbb{I}\Omega_{BI})\times\Omega_{BI})+\mathbb{I}^{-1}\tau_{g}+w_{\Omega,C},$
(13b)

where $(\mathbf{q}_{BO},\Omega_{BI})\in\mathbb{H}\times\mathbb{R}^{3}$. This
is the stably embedded version (7) of the satellite system model (10), i.e.
$\operatorname{S}^{3}\times\mathbb{R}^{3}$ is an invariant manifold of the
modified system and is exponentially stable for the modified system; refer to
Lemma 1 in Ko et al’s paper.24 Then, we may apply the standard UKF derived in
Euclidean space directly to the system model (13) and the measurement models
given in (11) and (12).

Since the satellite system in nature is a continuous-time system, we generate
a continuous-time trajectory of the satellite using (10). We integrate (10)
for $200[$\mathrm{s}$]$ and take measurements at discrete time steps of
$t=kh,k\in\mathbb{N}_{\geq 0}$ with a time step of
$h=\frac{1}{f}=0.1[$\mathrm{s}$]$. For each type of UKF, we conduct one
thousand Monte Carlo runs with each run simulated for 200 seconds. While
implementing UKF loops, we use the covariance matrix of the initial state
$P_{0}=\operatorname{diagonal}(10^{-8}\cdot I_{4},\>0.01\cdot I_{3})$ and the
covariance matrix of the process noise
$P_{w,D}=\operatorname{diagonal}(10^{-8}\cdot
I_{4},\>h\sigma^{2}_{\Omega,C}\cdot I_{3})$. We calculate the mean squared
error (MSE) of attitude estimates, and mean distance (MD) between attitude
estimates and $\operatorname{S}^{3}$ which are calculated as follows,
respectively:

$\displaystyle MSE_{\mathbf{q}}$
$\displaystyle=E[\|\mathbf{q}(kh)-\hat{\mathbf{q}}_{k}\|^{2}]\approx\frac{1}{N}\sum_{i=1}^{N}\|\mathbf{q}^{(i)}(kh)-\hat{\mathbf{q}}_{k}^{(i)}\|^{2},$
(14a) $\displaystyle MD_{\mathbf{q},\operatorname{S}^{3}}$
$\displaystyle=E[\left|\|\hat{\mathbf{q}}_{k}\|-1\right|]\approx\frac{1}{N}\sum_{i=1}^{N}\left|\|\hat{\mathbf{q}}^{(i)}_{k}\|-1\right|,$
(14b)

where $\mathbf{q}(kh)$ is the real attitude at time $t=kh$ in the continuous-
time trajectory, $\hat{\mathbf{q}}_{k}$ is the attitude estimate at time
$t=kh$ using each UKF, $E\left[\cdot\right]$ is defined as the expectation
with respect to $N=1000$ simulations, and the superscript $(i)$ stands for the
$i$-th trial of the $N$ simulations. $MSE_{\mathbf{q}}$ represents estimation
error of the attitude estimates. $MD_{\mathbf{q},\operatorname{S}^{3}}$
represents how close the attitude estimates are from $\operatorname{S}^{3}$
which is the unit quaternion set. The best $\alpha$ with a minimum estimation
error is dependent on the system and the time step, and a general guideline
for choosing a value of $\alpha$ is provided in Appendix B.

We conduct comparison simulation of our UKF-SE with other existing methods of
applying UKFs to systems on manifolds introduced in Section 2. We apply the
following types of filters on the sampled discrete-time signals of the
satellite system model:

  1. (i)

The standard UKF: Extend the satellite system model (10) from
$\operatorname{S}^{3}\times\mathbb{R}^{3}$ to $\mathbb{H}\times\mathbb{R}^{3}$
and apply the standard UKF discussed in Appendix A.

  2. (ii)

Projection by optimization: This is the forced projection method introduced in
Section 2. We project the attitude estimates $\hat{\mathbf{q}}_{k|k}$ to the
manifold $\operatorname{S}^{3}$ at the end of each estimation step of (i) by
solving
$\hat{\mathbf{q}}_{new,k|k}=\arg\max_{\mathbf{q}\in\operatorname{S}^{3}}\left\|\mathbf{q}-\hat{\mathbf{q}}_{k|k}\right\|^{2}$.
We use the Lagrangian multiplier method with the equality constraint
$\left\|\mathbf{q}\right\|^{2}-1=0,\;\mathbf{q}\in\mathbb{H}$ to solve this
optimization problem.18 We use the MATLAB command ‘fmincon’ to implement the
optimization by Lagrangian algorithm.

  3. (iii)

Projection by formula: This is for the specific case when the system is
defined on $\operatorname{S}^{3}$. After the update is over in (i), we project
the state estimate $\hat{\mathbf{q}}_{k|k}$ to the manifold
$\operatorname{S}^{3}$ by using the formula
$\hat{\mathbf{q}}_{new,k|k}=\hat{\mathbf{q}}_{k|k}/\left\|\hat{\mathbf{q}}_{k|k}\right\|$.

  4. (iv)

UKF-SE: This is our proposed UKF. From (21) in Appendix B and $h=0.1$,
$\alpha$ should satisfy $0<\alpha<\frac{20}{3}$. We choose $\alpha=2$ for our
UKF-SE.

  5. (v)

Lie algebra: This UKF is devised specifically for the systems defined on Lie
groups 8 as introduced in Section 2. The specific discretized form (2) for the
quaternion attitude system is introduced in the paper by Kraft.25 This UKF
needs several coordinate transformations between $\operatorname{S}^{3}$ and
its Lie algebra at each estimation step. The Lie algebra of
$\operatorname{S}^{3}$ is isomorphic to $\mathbb{R}^{3}$. Given a unit
quaternion $\mathbf{q}=q_{real}+\mathbf{q}_{imag}\in\operatorname{S}^{3}$
where $q_{real}$ is the real part and $\mathbf{q}_{imag}$ is the imaginary
part of the unit quaternion, the coordinate transformation of $\mathbf{q}$
into the Lie algebra is given as

$\mathbf{v}=\arctan\left(\frac{\left\|\mathbf{q}_{imag}\right\|}{q_{real}}\right)\frac{\mathbf{q}_{imag}}{\left\|\mathbf{q}_{imag}\right\|},$
(15)

where $\mathbf{v}\in\mathbb{R}^{3}$. Given an element of the Lie algebra
$\mathbf{v}\in\mathbb{R}^{3}$, the coordinate transformation of $\mathbf{v}$
into $\operatorname{S}^{3}$ is given as

$\mathbf{q}=\cos(\left\|\mathbf{v}\right\|)+\sin(\left\|\mathbf{v}\right\|)\frac{\mathbf{v}}{\left\|\mathbf{v}\right\|},$
(16)

where $\mathbf{q}\in\operatorname{S}^{3}$, and $\mathbf{v}$ is considered as a
pure imaginary quaternion. The map (15) is the inverse map of (16); see Sipos8
for more detail.

  6. (vi)

Tangent space: This UKF is devised for the systems defined on Riemannian
manifolds 9 as introduced in Section 2. It needs several coordinate
transformations between $\operatorname{S}^{3}$ and its tangent spaces at each
estimation step. We introduce the coordinate transformation formulas of the
Riemannian manifold $\operatorname{S}^{3}$, logarithm map and exponential map,
from the paper by Hauberg.26 Let us choose a base point
$\mathbf{\mu}\in\operatorname{S}^{3}$ and denote the tangent space of
$\operatorname{S}^{3}$ at $\mathbf{\mu}$ by
$T_{\mathbf{\mu}}\operatorname{S}^{3}=\\{\mathbf{u}\in\mathbb{R}^{4}\mid\langle\mathbf{u},\,\mathbf{\mu}\rangle=0\\}$.
A point $\mathbf{q}\in\operatorname{S}^{3}$ can be mapped to
$T_{\mathbf{\mu}}\operatorname{S}^{3}$ via the logarithm map,

$\displaystyle\log_{\mathbf{\mu}}(\mathbf{q})$
$\displaystyle=(\mathbf{q}-\mathbf{\mu}(\mathbf{q}^{T}\mathbf{\mu}))\frac{\theta}{\sin(\theta)},$
$\displaystyle\theta$ $\displaystyle=\arccos(\mathbf{q}^{T}\mathbf{\mu}),$
(17)

with the convention $0/\sin(0)=1$. A tangent vector $\mathbf{u}\in
T_{\mathbf{\mu}}\operatorname{S}^{3}$ can be mapped to $\operatorname{S}^{3}$
via the exponential map,

$\exp_{\mathbf{\mu}}(\mathbf{u})=\mathbf{\mu}\cos(\left\|\mathbf{u}\right\|)+\sin(\left\|\mathbf{u}\right\|)\frac{\mathbf{u}}{\left\|\mathbf{u}\right\|}.$
(18)

We use the MVIRT toolbox 27 to compute these coordinate transformations
between $\operatorname{S}^{3}$ and its tangent spaces, and weighted
interpolation of points on $\operatorname{S}^{3}$.

Figure 1: Mean squared error of $\hat{\mathbf{q}}$ versus time $t\,[$\mathrm{s}$]$ of each UKFs. Table 2: Average $MSE_{\mathbf{q}}$, average $MD_{\mathbf{q},\operatorname{S}^{3}}$, and computation time of 1000 simulations. Filters | Average $MSE_{\mathbf{q}}$ | Average $MD_{\mathbf{q},\operatorname{S}^{3}}$ | Computation time [$\mathrm{s}$]
---|---|---|---
The standard UKF | $5.513\times 10^{-7}$ | $4.026\times 10^{-4}$ | 798.5
Projection by optimization | $3.873\times 10^{-7}$ | $2.310\times 10^{-4}$ | 5882.6
Projection by formula | $2.972\times 10^{-7}$ | $3.523\times 10^{-17}$ | 842.4
UKF-SE (ours) | $2.759\times 10^{-7}$ | $1.405\times 10^{-6}$ | 869.4
Lie algebra | $1.797\times 10^{-4}$ | $5.308\times 10^{-15}$ | 3315.1
Tangent space | $3.143\times 10^{-7}$ | $3.416\times 10^{-17}$ | 5832.2

We plot $MSE_{\mathbf{q}}$ of each UKF as a semi-log plot on Figure 1. The Lie
algebra filter has $MSE_{\mathbf{q}}$ of order $10^{-4}$ which is relatively
larger than other filters, because we need to transform the points on
$\operatorname{S}^{3}$ to coordinates of $T_{1}\operatorname{S}^{3}$ and vice
versa at each UKF loop. If points on $\operatorname{S}^{3}$ are far away from
the identity, then they are far away from $T_{1}\operatorname{S}^{3}$, which
naturally incurs large errors. So, we do not plot the Lie algebra filter
result. We also do not plot MSE of the angular velocity because there is a
negligible difference on MSE of $\Omega_{BI}$ estimates among the UKF’s. We
calculate the average $MSE_{\mathbf{q}}$ along time from $t=150[$\mathrm{s}$]$
to $t=200[$\mathrm{s}$]$ and record them on the 2nd column of Table 2. The
UKF-SE has the lowest estimation error for estimating attitude among the
compared UKF’s as shown on Figure 1 and Table 2, which shows numerically that
UKF-SE is the most accurate UKF for systems on manifolds. We calculate the
average $MD_{\mathbf{q},\operatorname{S}^{3}}$ along time from
$t=150[$\mathrm{s}$]$ to $t=200[$\mathrm{s}$]$ and write on the 3rd column of
Table 2. The average $MD_{\mathbf{q},\operatorname{S}^{3}}$ of UKF-SE is
0.0035 times the average $MD_{\mathbf{q},\operatorname{S}^{3}}$ of the
standard UKF, which shows that UKF-SE keeps the attitude estimates in close
proximity to the manifold due to the stable embedding. The average
$MD_{\mathbf{q},\operatorname{S}^{3}}$’s of the projection by formula filter,
the Lie algebra filter, and the Tangent space filter are the lowest among the
UKF’s because they are composed in a specific way that constrains the attitude
estimates on the manifold. However the average
$MD_{\mathbf{q},\operatorname{S}^{3}}$’s of them are not exactly 0 because of
numerical rounding errors.

We also calculate the computation time of each UKF for 1000 Monte Carlo runs
and record them on the 4th column of Table 2. The number of sigma points is 20
for all the simulated UKFs in our paper. We use the concept of precision
operations written in the paper by Brent28 to perform a thorough computing
complexity analysis of the UKFs. Let $M(n)$ be the number of single-precision
operations required to multiply $n$-bit integers. Elementary operations, such
as addition, multiplication, and reciprocals can be evaluated in $O(M(n))$
operations, with a relative error $O(2^{-n})$.28 Trigonometric functions, such
as sin, cos, arctan can be evaluated in $O(M(n)\log(n))$ operations, with a
relative error $O(2^{-n})$.28 The UKF-SE shows 1.004 times slower performance
than the standard UKF. The reason for UKF-SE taking more time than the
standard UKF is that UKF-SE carries out the additional subtraction of the
stable embedding term
$\alpha(\left\|\mathbf{q}_{k}\right\|^{2}-1)\mathbf{q}_{k}$ for each sigma
point in each UKF loop in comparison with the standard UKF. The subtraction of
the stable embedding term has 17 elementary operations, and there are 20 sigma
points. Hence, the UKF-SE computes $O(340M(n))$ more operations than the
standard UKF in each UKF loop. The projection by optimization filter shows the
slowest computation time among the simulated filters. The reason for this is
that the projection by optimization filter solves an optimization problem in
each UKF loop. The UKF-SE shows a minor difference in computation time from
the projection by formula filter. The projection by formula filter calculates
an additional normalization formula at the end of each UKF loop in comparison
with the standard UKF, which is composed of 11 elementary operations.
Therefore, the projection by formula filter computes $O(11M(n))$ more
operations than the standard UKF in each UKF loop. It is however noted that a
general nonlinear manifold, unlike the 3-sphere $\operatorname{S}^{3}$, may
not even have a closed-form of normalization formula or, may only have
computationally expensive normalization algorithms such as polar decomposition
for special orthogonal group. The Lie algebra filter has 3.81 times slower
performance than UKF-SE because the Lie algebra filter computes additional
coordinate transformations between the Lie group and its Lie algebra, (15) and
(16), for each sigma point at each UKF loop in comparison with the standard
UKF. Computation of (15) and (16) is composed of 44 elementary operations and
3 trigonometric operations. Hence, the Lie algebra filter computes
$O((880+60\log(n))M(n))$ more operations than the standard UKF at each UKF
loop. The tangent space filter has 6.71 times slower performance than UKF-SE
because the tangent space filter computes additional coordinate
transformations between the Riemannian manifold and its tangent spaces, ((vi))
and (18), for each sigma point at each UKF loop in comparison with the
standard UKF. Computation of ((vi)) and (18) is composed of 65 elementary
operations and 4 trigonometric operations. Therefore, the tangent space filter
computes $O((1300+80\log(n))M(n))$ more operations than the standard UKF in
each UKF loop.

##  5 Conclusions

We have devised a simple, accurate, and computationally efficient unscented
Kalman filter (UKF) for the state estimation of systems on manifolds. The idea
is to extend the given system to an ambient Euclidean space, apply the special
term which makes the modified system asymptotically stable on the manifold,
and apply the standard UKF derived in Euclidean space. Our proposed UKF is
simple because it uses the off-the-shelf standard UKF devised in Euclidean
space directly while keeping the state estimates in close proximity to the
manifold. We apply this filter to the satellite system model which has
attitude state defined on the unit quaternion space. We investigate the
performance of our devised UKF by comparing with existing methods of applying
UKF to systems on manifolds. Simulation on the satellite system model shows
that our UKF is the most accurate compared to the other UKFs. Our UKF also has
the minor computation time comparable to the standard UKF devised in Euclidean
space, because our filter does not need either any manifold-structure-
preserving discretization method or coordinate transformations of other UKFs.
We plan to expand our stable embedding technique to invariant Kalman
filters.29

## Acknowledgments

This work has been supported by the NRF grant funded by the Korea government
(MSIT) (2021R1A2C2010585).

##  Appendix A The standard UKF

The UKF uses the concept of unscented transform to find stochastic properties
of state estimates. The unscented transform uses sigma points which capture
the mean and the covariance of states accurate to the second order.6 We
briefly summarize the standard UKF algorithm devised in Euclidean space
described in the paper by Wan and Van der Merwe.7 A discrete-time stochastic
system is given as follows:

$\displaystyle x_{k+1}$ $\displaystyle=f_{D}(x_{k},u_{k},w_{k}),$
$\displaystyle y_{k}$ $\displaystyle=g_{D}(x_{k},u_{k},v_{k}),$

where $k\in\mathbb{N}_{\geq 0}$, the state is $x_{k}\in\mathbb{R}^{n_{x}}$,
the input is $u_{k}\in\mathbb{R}^{n_{u}}$, the output is
$y_{k}\in\mathbb{R}^{n_{y}}$, the process noise is
$w_{k}\in\mathbb{R}^{n_{w}}$, the measurement noise is
$v_{k}\in\mathbb{R}^{n_{v}}$, $x_{0}\sim(\overline{x}_{0},P_{0})$,
$w_{k}\sim(0,P_{w,D})$, and $v_{k}\sim(0,P_{v,D})$. It is assumed that
$x_{0}$, $w_{k}$, and $v_{k}$ are mutually uncorrelated white noises. We want
state estimates $\hat{x}_{k|k}$ for the system above. First, we initialize the
state estimate as $\hat{x}_{0|0}=\overline{x}_{0}$ and the covariance estimate
as $P_{0|0}=P_{0}$. Then, we define the augmented state
$x_{k}^{a}=[x_{k}^{T},w_{k}^{T},v_{k}^{T}]^{T}\in\mathbb{R}^{L}$ for
$k\in\mathbb{N}_{\geq 0}$ where $L=n_{x}+n_{w}+n_{v}$. The augmented state
estimate is given as
$\hat{x}_{k|k}^{a}=[\hat{x}_{k|k}^{T},0,0]^{T}\in\mathbb{R}^{L}$. The
augmented covariance estimate is given as
$P_{k|k}^{a}=\operatorname{diagonal}(P_{k|k},P_{w,D},P_{v,D})$. We iteratively
generate sigma points, do time update and do measurement update at each
$k\in\mathbb{N}_{\geq 1}$. The process for generating sigma points is as
follows:

$\displaystyle\sigma_{i}^{a}=[\sigma_{i}^{x},\sigma_{i}^{w},\sigma_{i}^{v}]\text{
where
}\sigma_{i}^{x}\in\mathbb{R}^{n_{x}},\sigma_{i}^{w}\in\mathbb{R}^{n_{w}},\sigma_{i}^{v}\in\mathbb{R}^{n_{v}}\text{
and }i=0,1,2,\cdots,2L,$ $\displaystyle\sigma_{0}^{a}=\hat{x}_{k-1|k-1}^{a},$
$\displaystyle\sigma_{i}^{a}=\hat{x}_{k-1|k-1}^{a}+\sqrt{L+\lambda}\left[\sqrt{P_{k-1|k-1}^{a}}\right]_{i}\quad
i=1,\cdots,L,$
$\displaystyle\sigma_{i}^{a}=\hat{x}_{k-1|k-1}^{a}-\sqrt{L+\lambda}\left[\sqrt{P_{k-1|k-1}^{a}}\right]_{i-L}\quad
i=L+1,\cdots,2L,$ $\displaystyle W_{0}^{(m)}=\lambda/(L+\lambda),$
$\displaystyle W_{0}^{(c)}=\lambda/(L+\lambda)+(1-\gamma^{2}+\beta),$
$\displaystyle W_{i}^{(m)}=W_{i}^{(c)}=\lambda/{2(L+\lambda)}\quad
i=1,\cdots,2L,$

where $\left[\sqrt{P_{k-1|k-1}^{a}}\right]_{i}$ is the $i^{th}$ column of the
matrix, and $\lambda=\gamma^{2}(L+\kappa)-L$. The parameters are adopted from
the paper by Wan and Van der Merwe7 which are $\gamma=0.01$, $\kappa=0$, and
$\beta=2$. The process for doing time update is as follows:

$\displaystyle\sigma_{i,k|k-1}^{x}=f_{D}(\sigma_{i}^{x},u_{k-1},\sigma_{i}^{w})\quad
i=0,\cdots,2L,$
$\displaystyle\hat{x}_{k|k-1}=\sum_{i=0}^{2L}W_{i}^{(m)}\sigma_{i,k|k-1}^{x},$
$\displaystyle
P_{k|k-1}=\sum_{i=0}^{2L}W_{i}^{(c)}[\sigma_{i,k|k-1}^{x}-\hat{x}_{k|k-1}][\sigma_{i,k|k-1}^{x}-\hat{x}_{k|k-1}]^{T},$
$\displaystyle y_{i,k|k-1}=g_{D}(\sigma_{i}^{x},u_{k-1},\sigma_{i}^{v})\quad
i=0,\cdots,2L,$
$\displaystyle\hat{y}_{k|k-1}=\sum_{i=0}^{2L}W_{i}^{(m)}y_{i,k|k-1}.$

The process for doing measurement update is as follows:

$\displaystyle
P_{yy}=\sum_{i=0}^{2L}W_{i}^{(c)}[y_{i,k|k-1}-\hat{y}_{k|k-1}][y_{i,k|k-1}-\hat{y}_{k|k-1}]^{T},$
$\displaystyle
P_{xy}=\sum_{i=0}^{2L}W_{i}^{(c)}[\sigma_{i,k|k-1}^{x}-\hat{x}_{k|k-1}][y_{i,k|k-1}-\hat{y}_{k|k-1}]^{T},$
$\displaystyle K=P_{xy}P_{yy}^{-1},$
$\displaystyle\hat{x}_{k|k}=\hat{x}_{k|k-1}+K(y_{k}-\hat{y}_{k|k-1}),$
$\displaystyle P_{k|k}=P_{k|k-1}-KP_{yy}K^{T},$

where $y_{k}$ is the measurement and $\hat{x}_{k|k}$ is the state estimate at
time step $k$.

##  Appendix B A guideline for choosing a value of $\alpha$

We are given the following equation with stable embedding (13) from Section 4:

$\dot{\mathbf{q}}_{BO}=\frac{1}{2}\Lambda(\mathbf{q}_{BO})[0;\Omega_{BO}]-\alpha(\left\|\mathbf{q}_{BO}\right\|^{2}-1)\mathbf{q}_{BO}.$

We discretize the above equation using the ordinary Euler discretization with
time step $h>0$ as below

$\displaystyle\mathbf{q}_{BO,k+1}$
$\displaystyle=\mathbf{q}_{BO,k}+h\left(\frac{1}{2}\Lambda(\mathbf{q}_{BO,k})[0;\Omega_{BO,k}]-\alpha(\left\|\mathbf{q}_{BO,k}\right\|^{2}-1)\mathbf{q}_{BO,k}\right),$
$\displaystyle=(1-\alpha
h(\left\|\mathbf{q}_{BO,k}\right\|^{2}-1))\mathbf{q}_{BO,k}+\frac{h}{2}\Lambda(\mathbf{q}_{BO,k})[0;\Omega_{BO,k}].$
(19)

######  Theorem B.1.

If we define the error state as below

$e_{k}=\left\|\mathbf{q}_{BO,k}\right\|^{2}-1,$

then the discrete error dynamics of (B) is given as

$e_{k+1}=-1+(1+e_{k})\left((1-\alpha
he_{k})^{2}+\frac{h^{2}\left\|\Omega_{BO,k}\right\|^{2}}{4}\right).$ (20)

######  Proof B.2.

It is easy to show the following identities:

$\displaystyle\langle\mathbf{q}_{BO,k},\,\Lambda(\mathbf{q}_{BO,k})[0;\Omega_{BO,k}]\rangle=0,$
$\displaystyle\left\|\Lambda(\mathbf{q}_{BO,k})[0;\Omega_{BO,k}]\right\|^{2}=\left\|\mathbf{q}_{BO,k}\right\|^{2}\left\|\Omega_{BO,k}\right\|^{2}.$

With the above identities and the definition of the error $e_{k}$,

$\displaystyle e_{k+1}$
$\displaystyle=\left\|\mathbf{q}_{BO,k+1}\right\|^{2}-1,$
$\displaystyle=(1-\alpha
h(\left\|\mathbf{q}_{BO,k}\right\|^{2}-1))^{2}\left\|\mathbf{q}_{BO,k}\right\|^{2}+\frac{h^{2}}{4}\left\|\mathbf{q}_{BO,k}\right\|^{2}\left\|\Omega_{BO,k}\right\|^{2}-1,$

which simplifies to (20).

The next theorem gives a general guideline for choosing a value of $\alpha$.

######  Theorem B.3.

Given $\left|e_{k}\right|\leq\epsilon<1$, and
$\left\|\Omega_{BO,k}\right\|\leq U\;\forall k$, if we choose $\alpha$, $h$,
and $\beta$ that satisfy

$0<\alpha h<\frac{2}{3},$ (21) $0<\beta<\alpha h(1-\epsilon)(2-\alpha
h\epsilon)\left(2-2\alpha h-\alpha h(2-\alpha h)\epsilon+(\alpha
h)^{2}\epsilon^{2}\right)<1,$ (22)

and

$\delta(h,\alpha h,\beta)<\epsilon,$ (23)

where

$\delta(h,\alpha h,\beta)=h\sqrt{\frac{(1-2\alpha h(1-\epsilon)+(\alpha
h)^{2}\epsilon(1+\epsilon))\epsilon(1+\epsilon)\frac{U^{2}}{2}+(1+\epsilon)^{2}\frac{h^{2}U^{4}}{16}}{\alpha
h(1-\epsilon)(2-\alpha h\epsilon)\left(2-2\alpha h-\alpha h(2-\alpha
h)\epsilon+(\alpha h)^{2}\epsilon^{2}\right)-\beta}},$

then the set

$S_{\delta}=\left\\{e_{k}\in\mathbb{R}\mid\left|e_{k}\right|<\delta(h,\alpha
h,\beta)\right\\}$ (24)

becomes a positively invariant attractor of the discrete error dynamics (20)
with a region of attraction

$S_{\epsilon}=\left\\{e_{k}\in\mathbb{R}\mid\left|e_{k}\right|\leq\epsilon\right\\}.$

######  Proof B.4.

From the condition (21), $\alpha h$ always satisfies the below inequality

$0<\alpha
h<\operatorname{min}\left(\frac{2}{2\epsilon+1},\;\frac{1+\epsilon-\sqrt{1-\epsilon^{2}}}{\epsilon(1+\epsilon)}\right),$

because the infimum of the right side is $\frac{2}{3}$ for $0<\epsilon<1$.
With $\left|e_{k}\right|\leq\epsilon<1$, the following inequalities hold

$\displaystyle 0<1-\epsilon\leq 1+e_{k}<1+\epsilon,$ $\displaystyle 0<2-\alpha
h\epsilon\leq 2-\alpha he_{k},$ $\displaystyle 0<2-2\alpha h-\alpha h(2-\alpha
h)\epsilon+(\alpha h)^{2}\epsilon^{2}\leq 2-2\alpha h-\alpha h(2-\alpha
h)e_{k}+(\alpha h)^{2}e_{k}^{2}.$ (25)

Choose a Lyapunov function candidate $V(e_{k})=e_{k}^{2}$. From (20) and
(B.4),

$\displaystyle V(e_{k+1})-V(e_{k})$ $\displaystyle=e_{k+1}^{2}-e_{k}^{2}$
$\displaystyle=-\alpha h(1+e_{k})(2-\alpha he_{k})\left(2-2\alpha h-\alpha
h(2-\alpha h)e_{k}+(\alpha h)^{2}e_{k}^{2}\right)e_{k}^{2}$
$\displaystyle\quad+(1-2\alpha h(1+e_{k})+(\alpha
h)^{2}e_{k}(1+e_{k}))e_{k}(1+e_{k})\frac{h^{2}\left\|\Omega_{BO,k}\right\|^{2}}{2}$
$\displaystyle\quad+(1+e_{k})^{2}\frac{h^{4}\left\|\Omega_{BO,k}\right\|^{4}}{16}$
$\displaystyle\leq-\alpha h(1-\epsilon)(2-\alpha h\epsilon)\left(2-2\alpha
h-\alpha h(2-\alpha h)\epsilon+(\alpha h)^{2}\epsilon^{2}\right)e_{k}^{2}$
$\displaystyle\quad+(1-2\alpha h(1-\epsilon)+(\alpha
h)^{2}\epsilon(1+\epsilon))\epsilon(1+\epsilon)\frac{h^{2}U^{2}}{2}+(1+\epsilon)^{2}\frac{h^{4}U^{4}}{16}$
(26)

since $\left\|\Omega_{BO,k}\right\|\leq U$.

If $\left|e_{k}\right|$ satisfies the inequality

$\delta(h,\alpha h,\beta)<\left|e_{k}\right|\leq\epsilon,\\\ $ (27)

which is feasible by the condition (23), then from (B.4),
$V(e_{k+1})-V(e_{k})\leq-\beta e_{k}^{2}$. Consequently,
$V(e_{k+1})\leq(1-\beta)V(e_{k})$ when $e_{k}$ satisfies (27). Therefore,
$e_{k}$ decreases exponentially in magnitude so as to get into the set
$S_{\delta}$ in finite time, since $0<\beta<1$ from the condition (22)

If $\left|e_{k}\right|$ satisfies the inequality

$\delta(h,\alpha h,0)<\left|e_{k}\right|\leq\delta(h,\alpha h,\beta),$ (28)

where

$\delta(h,\alpha h,0)=\left.\begin{matrix}\delta(h,\alpha
h,\beta)\end{matrix}\right|_{\beta=0}=h\sqrt{\frac{(1-2\alpha
h(1-\epsilon)+(\alpha
h)^{2}\epsilon(1+\epsilon))\epsilon(1+\epsilon)\frac{U^{2}}{2}+(1+\epsilon)^{2}\frac{h^{2}U^{4}}{16}}{\alpha
h(1-\epsilon)(2-\alpha h\epsilon)\left(2-2\alpha h-\alpha h(2-\alpha
h)\epsilon+(\alpha h)^{2}\epsilon^{2}\right)}},$

then $V(e_{k+1})-V(e_{k})\leq 0$ from (B.4), which means
$\left|e_{k+1}\right|\leq\left|e_{k}\right|$. Hence, $e_{k+1}\in S_{\delta}$
if $e_{k}$ satisfies (28).

If $\left|e_{k}\right|$ satisfies the inequality

$\left|e_{k}\right|\leq\delta(h,\alpha h,0),$ (29)

then by moving $V(e_{k})$ to the right side in (B.4) and using (29),

$\displaystyle V(e_{k+1})$ $\displaystyle\leq\left(1-\alpha
h(1-\epsilon)(2-\alpha h\epsilon)\left(2-2\alpha h-\alpha h(2-\alpha
h)\epsilon+(\alpha h)^{2}\epsilon^{2}\right)\right)e_{k}^{2}$
$\displaystyle\quad+(1-2\alpha h(1-\epsilon)+(\alpha
h)^{2}\epsilon(1+\epsilon))\epsilon(1+\epsilon)\frac{h^{2}U^{2}}{2}+(1+\epsilon)^{2}\frac{h^{4}U^{4}}{16}$
$\displaystyle\leq\left(1-\alpha h(1-\epsilon)(2-\alpha
h\epsilon)\left(2-2\alpha h-\alpha h(2-\alpha h)\epsilon+(\alpha
h)^{2}\epsilon^{2}\right)\right)\times\delta(h,\alpha h,0)^{2}$
$\displaystyle\quad+(1-2\alpha h(1-\epsilon)+(\alpha
h)^{2}\epsilon(1+\epsilon))\epsilon(1+\epsilon)\frac{h^{2}U^{2}}{2}+(1+\epsilon)^{2}\frac{h^{4}U^{4}}{16}$
$\displaystyle=\delta(h,\alpha h,0)^{2},$ (30)

which implies

$\left|e_{k+1}\right|\leq\delta(h,\alpha h,0),$

if $e_{k}$ satisfies (29). Therefore, $e_{k+1}\in S_{\delta}$ if $e_{k}$
satisfies (29). To sum up, every trajectory of the discrete error dynamics
(20) starting in the set $S_{\delta}$ remains in $S_{\delta}$ for all forward
time, and every trajectory starting in $S_{\epsilon}$ enters the set
$S_{\delta}$ in finite time and remains in $S_{\delta}$ from then on. In other
words, $S_{\delta}$ is a positively invariant attractor of the error dynamics
with $S_{\epsilon}$ a region of attraction.

######  Remark B.5.

For any $0<\epsilon<1$ and $U>0$, we can find $\alpha$, $h$, and $\beta$ that
satisfy the conditions (21), (22), and (23). First, we can always find a value
for $\alpha h$ that satisfies the third inequality in (22) by choosing a
sufficiently small $\alpha h$, for example. Choose a value of $\beta$ that
satisfies the first and second inequalities of (22). Once $\alpha h$ and
$\beta$ are fixed, we can always choose a value of $h$ that satisfies (23).
Then, find a value of $\alpha$ that satisfies (21).

## References

  * 1 Fresk E, Nikolakopoulos G. Full quaternion based attitude control for a quadrotor. In: European Control Conference. IEEE. ; 2013; Zurich, Switzerland: 3864–3869.
  * 2 Guo A, Zhou Z, Zhu X, Bai F. Low-Cost Sensors State Estimation Algorithm for a Small Hand-Launched Solar-Powered UAV. Sensors 2019; 19(21): 4627.
  * 3 Soken HE, Hajiyev C, Sakai Si. Robust Kalman filtering for small satellite attitude estimation in the presence of measurement faults. European Journal of Control 2014; 20(2): 64–72.
  * 4 Jeon S, Tomizuka M, Katou T. Kinematic Kalman filter (KKF) for robot end-effector sensing. Journal of Dynamic Systems, Measurement, and Control 2009; 131(2).
  * 5 Rentmeesters Q, Absil PA, Van Dooren P. A filtering technique on the Grassmann manifold. In: Proceedings of the 19th International Symposium on Mathematical Theory of Networks and Systems. ; 2010; Budapest, Hungary.
  * 6 Julier SJ, Uhlmann JK. New extension of the Kalman filter to nonlinear systems. In: Signal Processing, Sensor Fusion, and Target Recognition VI. International Society for Optics and Photonics. ; 1997; Orlando, FL, United States: 182–193.
  * 7 Wan EA, Van Der Merwe R. The unscented Kalman filter for nonlinear estimation. In: Proceedings of the IEEE 2000 Adaptive Systems for Signal Processing, Communications, and Control Symposium. IEEE. ; 2000; Lake Louise, Alberta, Canada: 153–158.
  * 8 Sipos BJ. Application of the manifold-constrained unscented Kalman filter. In: IEEE/ION Position, Location and Navigation Symposium. IEEE. ; 2008; Monterey, CA, USA: 30–43.
  * 9 Hauberg S, Lauze F, Pedersen KS. Unscented Kalman filtering on Riemannian manifolds. Journal of mathematical imaging and vision 2013; 46(1): 103–120.
  * 10 Brossard M, Barrau A, Bonnabel S. A code for unscented kalman filtering on manifolds (UKF-M). In: International Conference on Robotics and Automation. IEEE. ; 2020: 5701–5708.
  * 11 Cantelobre T, Chahbazian C, Croux A, Bonnabel S. A real-time unscented Kalman filter on manifolds for challenging AUV navigation. In: International Conference on Intelligent Robots and Systems. IEEE. ; 2020: 2309–2316.
  * 12 Sjøberg AM, Egeland O. Lie algebraic unscented Kalman filter for pose estimation. IEEE Transactions on Automatic Control 2021: 1–1.
  * 13 Chang DE, Jiménez F, Perlmutter M. Feedback integrators. Journal of Nonlinear Science 2016; 26(6): 1693–1721.
  * 14 Park JH, Phogat KS, Kim W, Chang DE. Transversely stable extended Kalman filters for systems on manifolds in Euclidean spaces. Journal of Dynamic Systems, Measurement, and Control 2021; 143(6): 064501\.
  * 15 Kristiansen R, Nicklasson PJ. Satellite attitude control by quaternion-based backstepping. In: Proceedings of the 2005, American Control Conference. IEEE. ; 2005; Portland, OR, USA: 907–912.
  * 16 Frogerais P, Bellanger JJ, Senhadji L. Various ways to compute the continuous-discrete extended Kalman filter. IEEE Transactions on Automatic Control 2011; 57(4): 1000–1004.
  * 17 Sarkka S. On unscented Kalman filtering for state estimation of continuous-time nonlinear systems. IEEE Transactions on Automatic Control 2007; 52(9): 1631–1641.
  * 18 Chong EK, Zak SH. An introduction to optimization.  Canada: John Wiley & Sons .  2013\.  ISBN 978-1-118-27901-4.
  * 19 Hall B. Lie groups, Lie algebras, and representations: an elementary introduction.  Switzerland: Springer .  2015\.  ISBN 978-3-319-13466-6.
  * 20 Do Carmo MP. Riemannian geometry.  Boston: Birkhäuser .  1992\.  ISBN 978-0-8176-3490-2.
  * 21 Chang DE. On controller design for systems on manifolds in Euclidean space. International Journal of Robust and Nonlinear Control 2018; 28(16): 4981–4998.
  * 22 Lewis FL, Xie L, Popa D. Optimal and robust estimation: with an introduction to stochastic control theory.  Boca Raton, Florida, USA: CRC press .  2017\.  ISBN 978-0849390081.
  * 23 Hajiyev C, Soken HE. Robust adaptive unscented Kalman filter for attitude estimation of pico satellites. International Journal of Adaptive Control and Signal Processing 2014; 28(2): 107–120.
  * 24 Ko W, Phogat KS, Petit N, Chang DE. Tracking controller design for satellite attitude under unknown constant disturbance using stable embedding. Journal of Electrical Engineering & Technology 2021; 16(2): 1089–1097.
  * 25 Kraft E. A quaternion-based unscented Kalman filter for orientation tracking. In: Proceedings of the Sixth International Conference of Information Fusion. IEEE. ; 2003: 47–54.
  * 26 Hauberg S. Directional statistics with the spherical normal distribution. In: 21st International Conference on Information Fusion (FUSION). IEEE. ; 2018: 704–711.
  * 27 Bergmann R. MVIRT, a toolbox for manifold-valued image restoration. In: IEEE International Conference on Image Processing. IEEE. ; 2017; Beijing, China: 201–205.
  * 28 Brent RP. Fast multiple-precision evaluation of elementary functions. Journal of the ACM (JACM) 1976; 23(2): 242–251.
  * 29 Trawny N, Roumeliotis SI. Indirect Kalman filter for 3D attitude estimation. University of Minnesota, Dept. of Comp. Sci. & Eng., Tech. Rep 2005; 2: 2005.


# Fast Event-based Double Integral for Real-time Robotics

Shijie Lin1,5, Yingqiang Zhang1, Dongyue Huang2,5, Bin Zhou3,5, Xiaowei Luo4
and Jia Pan1,† 1S. Lin, Y. Zhang, and J. Pan are with the Department of
Computer Science, The University of Hong Kong, Hong Kong SAR, China.
<EMAIL_ADDRESS>Huang is with The Department of
Mechanical and Automation Engineering, The Chinese University of Hong Kong,
Hong Kong SAR, China<EMAIL_ADDRESS>Zhou is with the School of
Computer Science and Engineering, Beihang University, Beijing, China.
<EMAIL_ADDRESS>Luo is with the Department of Architecture and Civil
Engineering, City University of Hong Kong, Hong Kong SAR, China.
<EMAIL_ADDRESS>Peng Cheng Laboratory, Shenzhen, Guangdong, China.†
Corresponding author.

###### Abstract

Motion deblurring is a critical ill-posed problem that is important in many
vision-based robotics applications. The recently proposed event-based double
integral (EDI) provides a theoretical framework for solving the deblurring
problem with the event camera and generating clear images at high frame-rate.
However, the original EDI is mainly designed for offline computation and does
not support real-time requirement in many robotics applications. In this
paper, we propose the fast EDI, an efficient implementation of EDI that can
achieve real-time online computation on single-core CPU devices, which is
common for physical robotic platforms used in practice. In experiments, our
method can handle event rates at as high as 13 million event per second in a
wide variety of challenging lighting conditions. We demonstrate the benefit on
multiple downstream real-time applications, including localization, visual tag
detection, and feature matching.

##  I Introduction

Clear image acquisition is the prerequisite of various vision-based algorithms
to operate in robotics systems. However, due to the inherent blurry effect of
the active pixel sensor, motion artifacts can heavily degrade the image when
relative motion exists between the camera and scenes. This effect is further
intensified in low light conditions due to the demand for longer exposure
time. Condisering only frames, motion deblurring is an ill-posed problem [1]
and poorly addressed by existing methods [2, 3, 4]. Luckily, a novel type of
bio-inspired neuromorphic vision sensor, _i.e_., the dynamic vision sensor
(DVS) [5, 6], can encode light variation in high temporal resolution events,
allowing the motion deblurring being effectively addressed [7, 8, 9, 10, 11].

Pan _et al_. [7] proposed the event-based double integral (EDI) model that
bridges the formation of events and images, showing the effective capability
in deblurring and high-rate video reconstruction. EDI was later adopted by
various works [12, 13, 14] for offline processing. However, for real-time
robotics, existing EDI implementations are too slow ($1.5\text{\,}\mathrm{s}$
per image [7]). To employ the promising deblurring capability brought by EDI
for various robotics applications, we need to run it fast, making millions of
events to be processed online without jamming the whole system. However, this
is still quite challenging for the following reasons.

  * •

High Computational Complexity. Vanilla EDI [7, 15] is formulated using latent
images. Without modification, the computation complexity is coupled with the
image size and requires multiple redundant frame-wise operations,
significantly limiting its efficiency.

  * •

Infeasible to Real-time. When running EDI for real-time robotic systems,
millions of events must be processed within frame intervals to ensure
continuous captured data won’t jam the whole system. However, the original
implementation of EDI was targeted for offline processing and its code is
difficult to be adapted for online and onboard robotic system.

  * •

Unknown Contrast Parameters. To deal with the contrast parameter, related
works using EDI for image deblurring usually involve time-consuming
optimization [7, 16] or learning-based inference [17, 8], both of which are
highly time-consuming and infeasible for real-time processing with CPU-only
devices.

In this work, we propose fast EDI, an efficient reformulated version of EDI
that achieves real-time deblurring with an event rate up to 13 Million Ev/s in
a single core CPU. To achieve this performance, we propose to evenly
distribute the workload for online processing and attain the goal by
reformulating the EDI model and implementing it using a novel list-based
container. We further introduce a method to estimate the contrast from the
hardware parameters allowing robust and efficient parameter determination. In
experiments, we compare our methods with existing approaches in terms of
perceptual quality and runtime, showing state-of-the-art performance. Our fast
EDI allows real-time deblurring and benefits various downstream applications,
including visual tag detection in highly dynamic scenes, SLAM and feature
matching in low light conditions. We fully release the real-time processing
code and dataset with detailed hardware parameters.

Contribution: We propose the fast EDI, an efficient way of implementing the
event-based double integral (EDI), unlocking the high-speed sensing capability
of event cameras for clear image acquisition. We demonstrate that the EDI can
be implemented in real-time using a single core CPU, enabling the use of
various applications like localization, mapping, tracking, _etc_. We release
an efficient C++ implementation and dataset with hardware settings. Webpage:
github.com/eleboss/fast_EDI

##  II Related Work

Motion deblurring and reconstruction. The problems of motion deblurring and
image reconstruction have been widely researched over decades. Early
approaches are mostly based on the assumption of the static scene and
leverages gradient [18, 19, 20] or non-gradient priors [21, 22, 23]. Recent
works improved the performance through learning-based methods [2, 3, 4]. But
their performance is generally not guaranteed under all conditions.
Challenging motions or non-informative scenes could lead to degraded
performance. With events, the motion deblurring becomes much easier to
achieve. Early works fuse events and images using the asynchronous
complementary filter [24], manifold regularisation [25], or direct integration
[26]. But as the image formation model is not considered, their performance is
generally worse than EDI [7], which takes into account the event-frame
generation model and shows promising performance in image deblurring and high-
rate video reconstruction. Based on EDI, later works improved their
performance using learning frameworks. Songnan _et al_. [27] leveraged the
discrete EDI as a physical insight when designing the learning framework,
which was later improved by Xu _et al_. [11] by employing the optical flow of
the event and frame for learning and training on the real data. Recently,
Xiang and Yu [8] leveraged the EDI model to design a learnable double integral
network, which can provide a generalizable model for deblurring and frame
interpolation. However, previous works that rely on EDI require time-consuming
optimization or learning-based inference, making real-time processing
infeasible for robots with only CPU mounted. Even worse, according to our
knowledge, there is no solution existing for efficient online processing of
EDI. These issues significantly limit the usage of EDI in real-time robotics
like localization, feature tracking, visual tag detection, _etc_.

EDI related applications. As an important model for event-frame relation, EDI
has been adopted in various applications, including image denoising [16], high
dynamic range imaging [28], event-based super-resolution [14], human-pose
estimation [12], event-based optical flow estimation [29], and depth
estimation [30]. However, existing works primarily work in an offline manner
that does not require real-time processing. Whenever the number of events for
processing within each second easily exceeds millions, EDI computation could
consume tremendous time and bottleneck the whole system. Our fast EDI can
overcome this limitation and allow the use of the EDI in a wide variety of
applications.

Figure 1: Workload distribution. Our fast event-based double integral (EDI)
can distribute the overall workload during acquisition time, achieving real-
time processing. Vanilla EDI processes all events after the acquisition, which
cannot complete before the next image is ready, leading to system jamming.

##  III Preliminaries

We first briefly review the formation model of events and images, and then
formulate the event-based double integral (EDI) model.

#### Event generation

The event camera works asynchronously in responding to changes of log
intensity and triggers timestamped events whenever the log-scale intensity
change exceeds the contrast parameter $c>0$, _i.e_.,
$\log\left(L(\mathbf{x},t)\right)-\log\left(L(\mathbf{x},\tau)\right)=p\cdot
c,$ where $L(\mathbf{x},t)$ and $L(\mathbf{x},\tau)$ indicate the
instantaneous latent image at time $t$ and $\tau$ at pixel position
$\mathbf{x}$, and the polarity $p\in\\{+1,-1\\}$ denotes the direction of
intensity changes. Thus the event is a 4-dimensional tuple
$\mathbf{e}_{k}\triangleq(\mathbf{x}_{k},t_{k},p_{k})$, where
$\mathbf{x}_{k}=(x_{k},y_{k})^{T}$ is the pixel position, $t_{k}$ is the
triggering time and $p_{k}$ is the polarity of $k$-th events, respectively.
Finally, the brightness motion can be encode in the events set
$\mathcal{E}=\left\\{\mathbf{e}_{k}\mid k=1,\ldots,N_{ev}\right\\}$, where
$N_{ev}$ is the number of events. Leveraging events, we can give the
relationship between latent images $L(f)$ [7] (omitting the pixel positions
$\mathbf{x}$ for readability):

$\displaystyle L(t)=L(f)\exp\left(c\int_{f}^{t}e(s)ds\right),$ (1)

where $e(t)\doteq p\cdot\delta(t-\tau)$ is the continuous sampling function of
events with the Dirac function $\delta(\cdot)$.

#### EDI model

On the other hand, the blurry images are formulated using latent images $L(t)$
within the exposure time $T$:

$\displaystyle B=\frac{1}{T}\int_{t\in\mathcal{T}}L(t)dt,$ (2)

where the blurry image is averaged results of latent images with the exposure
duration $\mathcal{T}$. Then the latent images is reformulated as [7]:
$L(f)=\frac{B}{E(f,\mathcal{T})}$, with

$\displaystyle
E(f,\mathcal{T})=\frac{1}{T}\int_{t\in\mathcal{T}}\exp{\left(c\int_{f}^{t}e(s)ds\right)dt}$
(3)

indicating the physical relation between latent images and blurry images from
concurrent events, which is also refer to the event-based double integral
(EDI) model.

Figure 2: EDI implemented with conventional array-like container requires
frame-wise additions for the outer integral and thus is inefficient.

(a) Step 1: inner integral (during exposure)

(b) Step 2: outer integral (after exposure).

Figure 3: EDI implemented with the proposed list-based container can
significantly reduce the number of addition operations and decouple it from
the image size, thus achieving efficient real-time processing. In step 1, the
inner integral operates pixel-wise and pushes back the latest latent pixel
value, _i.e_., $\text{Exp}(\cdot)$ and the counter value in the latent pixels
lists. In step 2, the outer integral could be computed by summing the
multiplication between latent pixel value (green exp value) and the difference
of two consecutive counter numbers.

#### Current issues

Although the continuous EDI model is given in previous literature [7, 15],
efficient online implementation was never specified. To achieve this goal, we
need to overcome the following difficulties:

  * •

Workload imbalance. The EDI computation requires an accurate selection to
separate events within the exposure interval. However, for online processing,
the camera generates data sequentially (Fig. 1), meaning the selection cannot
be made before the image’s exposure is completed. Thus, vanilla EDI starts
processing after the exposure is completed. However, the exposure time of an
image could exceed $30\text{\,}\mathrm{m}\mathrm{s}$, and millions of events
could be triggered during this time. In a real-time system, all data need to
be processed before subsequent data is ready. If the computation of vanilla
EDI cannot be done within the frame interval, the subsequent data will cause
jamming, slowing down the system and leading to an out-of-memory crash.

  * •

Redundant operations. Unfortunately, the vanilla EDI usually cannot complete
all computations within the frame interval due to redundant operations.
Specifically, computing Eq. 3 requires repeatedly frame-like reconstruction of
the latent image (inner integral) over the whole exposure interval (outer
integral), making the real-time processing highly challenging.

  * •

Unknown contrast parameter. Previous works [7, 15] use time-consuming
optimization or inaccurate manual setting to handle the contrast parameters
$c$. Both of them cannot satisfy the need for real-time robotics systems.

Therefore, we need to develop the whole computation pipeline from scratch to
implement the EDI for efficient online processing.

##  IV Methods

To improve the EDI for real-time processing, we developed the fast EDI, whose
is implementation is based on a novel list-based container to achieve a
balanced workload distribution. We further introduce the method of contrast
estimation using hardware parameter [31], allowing a robust and efficient
determination.

###  IV-A Fast event-based double integral (F-EDI)

####  IV-A1 Workload distribution

As the events are generated asynchronously from the sensor, one key to
boosting the processing speed is to make each event processed right after it
is triggered. However, for robotics systems, various applications are needed
to consume the processing power of the CPU, and cores in the low-end processor
are limited. Thus it is best if we can complete all computation using a single
thread. To do that, we first compute the inner integral and exponentiation for
each event _i.e_. $\exp\left(c\int_{f}^{t}e(s)ds\right)$, and record the
integral result for quick retrieval. Once the exposure is completed, the outer
integral only needs to retrieve all containers and calculate the outer
integral. This workload distribution allows us to fully leverage the CPU’s
computing resource to process millions of events in real-time (Fig. 1),
thereby alleviating the jamming problem.

####  IV-A2 Efficient list-based EDI

However, even with the desired workload, the EDI computation is still
inefficient for real-time processing. The vanilla EDI is formulated as a
latent image, and the inner integral is responsible for rendering the latent
image. Then the outer integral needs to sum up all latent images repetitively.
However, since exponentiation exists between the inner and outer integral,
frame-wise additions are required for summing up all latent images. As shown
in Fig. 2, given the latest triggered events marked in red, the inner integral
renders the event images, then exponentiation transforms the event image into
the latent image, and the outer integral need to sum up all latent images. The
first two operations, _i.e_. inner integral and exponentiation, are efficient
as we can use the pixel-wise operation. But after the exponentiation, the
latent image is nonlinear. Thus the outer integral needs to repetitively sums
up all latent images in a frame-wise manner. Meaning for each event, equal to
frame size $N_{x}$ additions are required for the outer integral leading to
$\mathcal{O}(N_{ev}\times N_{x})$ complexity.

We find this process could be simplified by redesigning the container. Since
each time an event is generated, only one pixel alters its value in the latent
image, most pixels for the additions in outer integral remain unchanged. Thus,
we can temporally record each pixel’s inner integral results and the number of
additions during the exposure time. Then, we can convert the repetitive
additions to summing several multiplications in a single retrieval after the
exposure time. Therefore, we propose the list-based container, as shown in
Fig. 3. With this container, the fast EDI can be done in two steps. Step 1:
During the exposure time, it conducts the inner integral using the latest
activated events, taking its exponentiation, and increasing the global
addition counter by one. Then two values (_i.e_., $<$exp value, counter$>$)
will be pushed back in a list marked by specific pixels positions (the red
cube in Fig. 3). Step 2: After exposure, our fast EDI loops over all lists at
specific pixels to retrieve the integral values and compute the outer integral
in a each list by summing multiplications between the latent pixel value
(green exp value in Fig. 3)) and the difference of two consecutive counter
numbers (blue value in Fig. 3).

In this way, we can turn millions of additions into simply hundreds of
multiplications and additions, thereby reducing the processing time. Finally,
the time complexity is purely linear to the number of events, _i.e_.,
$\mathcal{O}(N_{ev})$.

###  IV-B Event Contrast from Hardware Parameters

Vanilla EDI [7] uses an energy function based on image sharpness for
optimization to find the optimal contrast. But such an optimization algorithm
has two drawbacks. First, if the highly blurred image lacks enough high-
contrast edges, the optimization is likely to converge to a local optimum.
Second, the optimization is inefficient, making it difficult to deploy on
real-time robotic systems. To address these issues, we introduce an estimation
method based on hardware parameters. From hardware settings, we can directly
calculate its contrast:

$\displaystyle
C_{\text{OFF}}=\frac{\kappa_{n}\mathcal{C}_{2}}{\kappa_{p}^{2}\mathcal{C}_{1}}\ln(\frac{\mathcal{I}_{\text{OFF}}}{\mathcal{I}_{d}}),C_{\text{ON}}=\frac{\kappa_{n}\mathcal{C}_{2}}{\kappa_{p}^{2}\mathcal{C}_{1}}\ln(\frac{\mathcal{I}_{\text{ON}}}{\mathcal{I}_{d}}),$
(4)

where $\kappa_{n}=\kappa_{p}=0.7$ are the back gate coefficients of n and p
FET transistors. $\mathcal{C}_{1}/\mathcal{C}_{2}$ is the capacitor ratio of
DVS ($130/6$ for our DAVIS346). $\mathcal{I}_{\text{d}}$,
$\mathcal{I}_{\text{ON}}$, and $\mathcal{I}_{\text{OFF}}$ are the bias current
set by the user in a coarse-fine bias generator [32], which can be computed
with jAER toolbox [33]. In this way, we can greatly simplify the optimization
and increase the efficiency.

##  V Experiments

In experiments, we first evaluate the primary performance of the fast EDI,
_i.e_., runtime and deblurring capability. Then we evaluate the fast EDI in
three real-time downstream applications, _i.e_., feature tracking, visual tag
detection, and SLAM. We encourage the reader to view the supplementary video
for more information.

TABLE I: Average scores of deblurring (Best: $\downarrow$)

Metric | Method
---|---
| [34] | [35] | [7] | Ours
PIQE [36] | 62.92 | 54.71 | 42.55 | 38.12
SSEQ [37] | 31.64 | 37.58 | 32.04 | 24.01
BRISQUE [38] | 32.53 | 37.39 | 35.04 | 28.75

TABLE II: Average matched inliers of feature tracking (Better: $\uparrow$)

Feature | Trajectory 1 | Trajectory 2 | Trajectory 3
---|---|---|---
| Blur | Ours | Blur | Ours | Blur | Ours
SURF [39] | 83.5 | 92.1 | 90.2 | 100.4 | 95.5 | 106.8
FAST [40] | 39.3 | 41.6 | 42.7 | 46.1 | 71.2 | 73.6
MSER [41] | 39.2 | 42.3 | 48.8 | 51.6 | 82.0 | 84.6
BRISK [42] | 68.3 | 67.8 | 73.2 | 73.3 | 112.2 | 112.4

(a) Events

(b) Images

(c) [34]

(d) [35]

(e) [7]

(f) Ours

Figure 4: Image deblurring results in two sequences, line (first row) and
scissor (second row). (c) Conventional method [34] and (d) learning-based
method [35] fail to recover (a) the heavily blur images. (e) [7] works well
when the background is clean but is unable to tackle heavily blurred images.
In contrast, (f) our method can stably generate clean images.

(a) Blurred image

(b) Our deblurred image

(c) Matching results enhanced by our fast EDI

Figure 5: Results of SURF [39] feature matching. Our method can effectively
deblur (a) the blurry image and generate clear images for (b) accurate feature
matching.

(a) Blur

(b) Result 1

(c) Result 2

Figure 6: Results of Apriltag [43] detection. Our method can deblur (a) the
blurry image cause by rapid shaking motions and generate (b) (c) clear images,
making consecutive successful Apriltag detection possible.

(a) Tracking results using original blur images

(b) Tracking results using our deblurred images

Figure 7: Tracking result using VINS-mono SLAM [44]. It could be clearly
observed that (a) the tracking using blurry images cannot extract effective
features. Most of them are in low confidence (blue dots). In contrast,
tracking using images deblurred by our methods has much more effective
features with high confidence (red dots).

###  V-A Experimental Setups

####  V-A1 Hardware

We conducted all the experiments using Intel NUC mini PC with Intel
<EMAIL_ADDRESS>GHz CPU, and DAVIS346 [6] with contrast parameter set as as
$C_{\text{ON}}=0.26$ and $C_{\text{OFF}}=-0.26$ to fit the single contrast
model used in [7]. The NUC mini PC is a common choice for various robotics
systems and has multi-CPU cores for running other robotics algorithms like
recognition, planning, and localization. We use C++ to implement our
algorithms on the DV Platform with a ROS node to publish the real-time
results.

####  V-A2 Dataset

Since existing datasets do not include details of hardware settings for
contrast estimation, during the real-time experiments, we recorded the online
data as an evaluation dataset, which contains different motions (_e.g_.,
shaking, random move) and scene illuminations (_e.g_., low lighting
conditions, sunlight, artificial light). Details are in the appendix.

###  V-B Runtime Performance

The vanilla EDI is implemented using Matlab with C++ warper and cannot support
further modification. Its performance also can not support real-time
processing (1.5 seconds per image, according to [7]). Thus, we compare our
fast EDI with EDI implemented with the conventional array-like container using
C++ and the same setup. During experiments, we record the average events rate
of both methods. Our fast EDI could exceed 13 million Ev/s event rates using a
single core CPU. In comparison, vanilla EDI can only achieve 49 thousand Ev/s.
This result shows our fast EDI is $260\times$ faster than vanilla EDI, which
demonstrates the efficiency of our method.

###  V-C Deblurring

We compare our approach with other state-of-the-art image deblurring methods,
including conventional method [34], learning-based method [35], and event-
frame reconstruction method [7]. Qualitative comparisons are shown in Fig. 4.
Multiple no-reference metrics, _i.e_., SSEQ [37], PIQE [36], and BRISQUE [38]
are adopted to provide quantitative comparisons, with the average scores on
our dataset shown in Table I, where the lower is the score, the better is the
performance. According to Table I, our method achieves the best score
regarding all metrics. Other methods get suboptimal scores because the heavily
blurred image can hardly provide any useful information for reconstruction.
From Fig. 4, we can see that both our method and [7] handle the partially
blurred images well, _e.g_., the blur of a moving scissor in Fig. 4(b).
However, for completely blurred images, only our method can reconstruct clean
images (Fig. 4(f)). [7] converges to sub-optimal contrast parameters, and the
other two methods [35, 34] fail to converge. It is worth noting that, with a
carefully manually adjusted contrast parameter, the deblurring performance of
vanilla EDI is identical to our fast EDI. Thus our method also shares the
defects of EDI, like easily affected by noise. Some of these defects have been
addressed in recently proposed learning-based methods [8], but such methods
rely on GPU for fast inference. As our work aims at introducing an efficient
implementation of EDI, comparison with them is out of the scope of this work.

###  V-D Real-time Applications

####  V-D1 Feature matching and tracking

Our method can also benefit the feature matching and tracking by improving
image quality. As shown in Fig. 5(a), images taken in low lighting suffer from
motion blur, which can barely provide any information for feature matching. In
contrast, our fast EDI recovers a clean image in Fig. 5(b), allowing
successful SURF [39] feature matching in Fig. 5(c). Quantitative comparison
results are summarized in Table II, where our fast EDI improves the matching
for all three types of features, which again verifies the effectiveness of our
method.

####  V-D2 Visual tag detection

Our fast EDI also benefits visual tag detection due to its effectiveness in
real-time image deblurring. As shown in Fig. 6(a), large motions will blur the
visual tag, while our reconstructed clean image can properly preserve
important information for successful Apriltag [43] detection in a sequence of
206 Apriltag images, as shown in Fig. 6(b) and Fig. 6(c). In particular, the
number of detected Apriltag increases from 73 to 198 after using our method.

####  V-D3 SLAM

Simultaneous Localization and Mapping (SLAM) is an important robotics
application that relies on real-time image processing. We use the VINS-Mono
([44], no loop closure) in a dark underground garage with a handheld DAVIS346
camera for evaluation. In experiments, we started and ended from the same
position marked by visual tags. As the environment is quite dark ($<$
$50\text{\,}\mathrm{L}\mathrm{u}\mathrm{x}$), the frame-based camera requires
a longer exposure time for imaging, introducing motion blur and destablizing
the tracking front-end (Fig. 7(a)). Our fast EDI can provide real-time
deblurred results to facilitate tracking (Fig. 7(b)), gaining a much higher
localization accuracy compared with trajectory using original blur image, as
shown in Fig. 8(b). Using blurry images, the translation error between
starting and ending point is
$44.16\text{\,}\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{e}\mathrm{r}$. This error
is reduced to
$5.22\text{\,}\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{e}\mathrm{r}$ with the
assistance of fast EDI.

(a) Trajectories of VINS-Mono (XY-axis)

(b) Trajectories of VINS-Mono (Z-axis) Figure 8: Result of SLAM using VINS-
mono [44]. VINS runs with images deblurred by our fast EDI (green trajectory)
has much lower drifting compared with results using original blurred image
(red trajectory). The starting and ending position are the same and marked by
a calibration board. Best view in color.

###  V-E Limitations

The memory consumption of this work increases linear to the number of
triggered events. We therefore suggest that a minimum of
$8\text{\,}\mathrm{G}\mathrm{B}$ memory be allocated for the system. When
using event cameras from manufacturers other than Inivation, it is noteworthy
that the jAER toolbox may not yield accurate contrast parameter values. Hence,
it is advised that users either contact the manufacturers directly for inquiry
or employ calibration methods [45].

##  VI Conclusion

In this paper, we propose the fast EDI, a novel way to efficiently implement
the EDI to achieve real-time image deblurring using event cameras. In
comparison with the conventional EDI model, we boost the computational speed
260 times on single core CPU and successfully bring the deblurring capability
of EDI in improving the performance of various real-time applications,
including SLAM, visual tag detection, and feature matching.

## References

  * [1] K. Purohit, A. Shah, and A. Rajagopalan, “Bringing alive blurred moments,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2019, pp. 6830–6839.
  * [2] J. Sun, W. Cao, Z. Xu, and J. Ponce, “Learning a convolutional neural network for non-uniform motion blur removal,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2015, pp. 769–777.
  * [3] X. Tao, H. Gao, X. Shen, J. Wang, and J. Jia, “Scale-recurrent network for deep image deblurring,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2018, pp. 8174–8182.
  * [4] L. Pan, Y. Dai, M. Liu, and F. Porikli, “Depth map completion by jointly exploiting blurry color images and sparse depth maps,” in _2018 IEEE Winter Conference on Applications of Computer Vision (WACV)_. IEEE, 2018, pp. 1377–1386.
  * [5] P. Lichtsteiner, C. Posch, and T. Delbruck, “A 128$\times$128 120 db 15$\mu$s latency asynchronous temporal contrast vision sensor,” _IEEE Journal of Solid-State Circuits_ , vol. 43, no. 2, pp. 566–576, 2008.
  * [6] C. Brandli, R. Berner, M. Yang, S.-C. Liu, and T. Delbruck, “A 240$\times$180 130 db 3 $\mu$s latency global shutter spatiotemporal vision sensor,” _IEEE Journal of Solid-State Circuits_ , vol. 49, no. 10, pp. 2333–2341, 2014\.
  * [7] L. Pan, C. Scheerlinck, X. Yu, R. Hartley, M. Liu, and Y. Dai, “Bringing a blurry frame alive at high frame-rate with an event camera,” in _2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2019, pp. 6813–6822.
  * [8] X. Zhang and L. Yu, “Unifying motion deblurring and frame interpolation with events,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2022, pp. 17 765–17 774.
  * [9] C. Scheerlinck, H. Rebecq, D. Gehrig, N. Barnes, R. Mahony, and D. Scaramuzza, “Fast image reconstruction with an event camera,” in _Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision_ , 2020, pp. 156–163.
  * [10] B. Wang, J. He, L. Yu, G.-S. Xia, and W. Yang, “Event enhanced high-quality image recovery,” in _European Conference on Computer Vision_. Springer, 2020, pp. 155–171.
  * [11] F. Xu, L. Yu, B. Wang, W. Yang, G.-S. Xia, X. Jia, Z. Qiao, and J. Liu, “Motion deblurring with real events,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV)_ , 2021, pp. 2583–2592.
  * [12] L. Xu, W. Xu, V. Golyanik, M. Habermann, L. Fang, and C. Theobalt, “Eventcap: Monocular 3d capture of high-speed human motions using an event camera,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 4968–4978.
  * [13] D. Gu, J. Li, Y. Zhang, and Y. Tian, “How to learn a domain-adaptive event simulator?” in _Proceedings of the 29th ACM International Conference on Multimedia_ , 2021, pp. 1275–1283.
  * [14] L. Wang, T.-K. Kim, and K.-J. Yoon, “Eventsr: From asynchronous events to image reconstruction, restoration, and super-resolution via end-to-end adversarial learning,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2020, pp. 8315–8325.
  * [15] L. Pan, R. Hartley, C. Scheerlinck, M. Liu, X. Yu, and Y. Dai, “High frame rate video reconstruction based on an event camera,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2020.
  * [16] Z. W. Wang, P. Duan, O. Cossairt, A. Katsaggelos, T. Huang, and B. Shi, “Joint filtering of intensity images and neuromorphic events for high-resolution noise-robust imaging,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2020, pp. 1609–1619.
  * [17] S. Lin, J. Zhang, J. Pan, Z. Jiang, D. Zou, Y. Wang, J. Chen, and J. Ren, “Learning event-driven video deblurring and interpolation,” in _European Conference on Computer Vision_. Springer, 2020, pp. 695–710.
  * [18] R. Fergus, B. Singh, A. Hertzmann, S. T. Roweis, and W. T. Freeman, “Removing camera shake from a single photograph,” in _ACM SIGGRAPH 2006 Papers_ , 2006, pp. 787–794.
  * [19] D. Krishnan, T. Tay, and R. Fergus, “Blind deconvolution using a normalized sparsity measure,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2011, pp. 233–240.
  * [20] L. Sun, S. Cho, J. Wang, and J. Hays, “Edge-based blur kernel estimation using patch priors,” in _IEEE international conference on computational photography (ICCP)_. IEEE, 2013, pp. 1–8.
  * [21] W.-S. Lai, J.-J. Ding, Y.-Y. Lin, and Y.-Y. Chuang, “Blur kernel estimation using normalized color-line prior,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2015, pp. 64–72.
  * [22] J. Pan, D. Sun, H. Pfister, and M.-H. Yang, “Deblurring images via dark channel prior,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 40, no. 10, pp. 2315–2328, 2017.
  * [23] Y. Yan, W. Ren, Y. Guo, R. Wang, and X. Cao, “Image deblurring via extreme channels prior,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2017, pp. 4003–4011.
  * [24] C. Scheerlinck, N. Barnes, and R. Mahony, “Continuous-time intensity estimation using event cameras,” in _Asian Conference on Computer Vision (ACCV)_. Springer, 2018, pp. 308–324.
  * [25] G. Munda, C. Reinbacher, and T. Pock, “Real-time intensity-image reconstruction for event cameras using manifold regularisation,” _International Journal of Computer Vision_ , vol. 126, no. 12, pp. 1381–1393, 2018.
  * [26] C. Brandli, L. Muller, and T. Delbruck, “Real-time, high-speed video decompression using a frame- and event-based davis sensor,” in _2014 IEEE International Symposium on Circuits and Systems (ISCAS)_ , 2014, pp. 686–689.
  * [27] S. Lin, J. Zhang, J. Pan, Z. Jiang, D. Zou, Y. Wang, J. Chen, and J. Ren, “Learning event-driven video deblurring and interpolation,” in _Proceedings of the European Conference on Computer Vision (ECCV)_. Springer, 2020, pp. 695–710.
  * [28] Z. Wang, Y. Ng, C. Scheerlinck, and R. Mahony, “An asynchronous kalman filter for hybrid event cameras,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV)_ , 2021, pp. 448–457.
  * [29] L. Pan, M. Liu, and R. Hartley, “Single image optical flow estimation with an event camera,” in _2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_. IEEE, 2020, pp. 1669–1678.
  * [30] D. Gehrig, M. Rüegg, M. Gehrig, J. Hidalgo-Carrió, and D. Scaramuzza, “Combining events and frames using recurrent asynchronous multimodal networks for monocular depth prediction,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 2822–2829, 2021.
  * [31] Y. Nozaki and T. Delbruck, “Temperature and parasitic photocurrent effects in dynamic vision sensors,” _IEEE Transactions on Electron Devices_ , vol. 64, no. 8, pp. 3239–3245, 2017.
  * [32] T. Delbruck, R. Berner, P. Lichtsteiner, and C. Dualibe, “32-bit configurable bias current generator with sub-off-current capability,” in _2010 IEEE International Symposium on Circuits and Systems (ISCAS)_ , 2010, pp. 1647–1650.
  * [33] “jaer project,” http://jaerproject.org.
  * [34] J. Pan, D. Sun, H. Pfister, and M.-H. Yang, “Deblurring images via dark channel prior,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 40, no. 10, pp. 2315–2328, 2018.
  * [35] M. Jin, G. Meishvili, and P. Favaro, “Learning to extract a video sequence from a single motion-blurred image,” in _2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2018, pp. 6334–6342.
  * [36] N. Venkatanath, D. Praneeth, M. C. Bh, S. S. Channappayya, and S. S. Medasani, “Blind image quality evaluation using perception based features,” in _2015 Twenty First National Conference on Communications (NCC)_ , 2015, pp. 1–6.
  * [37] L. Liu, B. Liu, H. Huang, and A. C. Bovik, “No-reference image quality assessment based on spatial and spectral entropies,” _Signal Processing: Image Communication_ , vol. 29, no. 8, pp. 856–863, 2014.
  * [38] A. Mittal, A. K. Moorthy, and A. C. Bovik, “No-reference image quality assessment in the spatial domain,” _IEEE Transactions on Image Processing_ , vol. 21, no. 12, pp. 4695–4708, 2012.
  * [39] H. Bay, A. Ess, T. Tuytelaars, and L. Van Gool, “Speeded-up robust features (surf),” _Computer vision and image understanding_ , vol. 110, no. 3, pp. 346–359, 2008.
  * [40] E. Rosten and T. Drummond, “Fusing points and lines for high performance tracking,” in _Tenth IEEE International Conference on Computer Vision (ICCV) Volume 1_ , vol. 2, 2005, pp. 1508–1515.
  * [41] J. Matas, O. Chum, M. Urban, and T. Pajdla, “Robust wide-baseline stereo from maximally stable extremal regions,” _Image and Vision Computing_ , vol. 22, no. 10, pp. 761–767, 2004.
  * [42] S. Leutenegger, M. Chli, and R. Y. Siegwart, “Brisk: Binary robust invariant scalable keypoints,” in _2011 International Conference on Computer Vision (ICCV)_ , 2011, pp. 2548–2555.
  * [43] E. Olson, “Apriltag: A robust and flexible visual fiducial system,” in _2011 IEEE International Conference on Robotics and Automation (ICRA)_ , 2011, pp. 3400–3407.
  * [44] T. Qin, P. Li, and S. Shen, “Vins-mono: A robust and versatile monocular visual-inertial state estimator,” _IEEE Transactions on Robotics_ , vol. 34, no. 4, pp. 1004–1020, 2018.
  * [45] Z. Wang, Y. Ng, P. van Goor, and R. Mahony, “Event camera calibration of per-pixel biased contrast threshold,” _arXiv preprint arXiv:2012.09378_ , 2020.