# Confidence and Assurance of Percentiles

Sanjay M. Joshi
Independent Researcher
<EMAIL_ADDRESS>
https://www.linkedin.com/in/sanjaymjoshi/

(December 31, 2023)

###### Abstract

Confidence interval of mean is often used when quoting statistics. The same
rigor is often missing when quoting percentiles and tolerance or percentile
intervals. This article derives the expression for confidence in percentiles
of a sample population. Confidence intervals of median is compared to those of
mean for a few sample distributions. The concept of assurance from reliability
engineering is then extended to percentiles. The assurance level of sorted
samples simply matches the confidence and percentile levels. Numerical method
to compute assurance using Brent’s optimization method is provided as an open-
source python package.

_Keywords_ Confidence $\cdot$ Assurance $\cdot$ Percentiles $\cdot$ Quantiles

##  1 Introduction

Statistics are typically qualified with confidence level information.
Confidence interval of mean is commonly used to qualify mean of a random
variable. Such intervals are not commonly seen in percentiles or quantiles,
though. Likewise, the expected range of values is commonly quoted as tolerance
interval, but the related confidence in that range is often missing.

This article presents a review of techniques to compute confidence levels,
starting with derivation of equation for confidence from [1]. It then defines
a new concept of "assurance interval" that combines percentiles interval with
confidence level. When probability distribution function of a random variable
is not known, a numerical method for easy computations of assurance interval
is presented. Resources in the form of an open-source package in Python are
also provided.

##  2 Confidence interval of percentiles

Let $X_{i}$ be $i$’th sample of a random variable $X$, $1<=i<=n$. Assume that
these samples are independent and identically distributed. Their probability
distribution function need not be known. Let $\pi_{p}$ be the $p$’th
percentile (or quantile) level, $0<p<1$, for $X$.

Define a new binary random variable $Y:X<=\pi_{p}$. It represents Bernoulli
trials and has a binomial distribution. Let $j$ out of $n$ observed samples of
$X$ be all less than or equal to $\pi_{p}$. Therefore, in $n$ samples, we have
observed $j$ successes and $n-j$ failures of $Y$. The confidence that the the
probability of success, or reliability, of $Y$ is at least $p$ is [2]

$c=1-\sum_{i=0}^{n-j}\binom{n}{i}(1-p)^{i}p^{n-i}$ (1)

Substitute $i=n-k$ and rearrange summation to get

$c=1-\sum_{k=j}^{n}\binom{n}{n-k}(1-p)^{n-k}p^{k}$

Since $\binom{n}{n-k}=\binom{n}{k}$, this is the same as

$c=\sum_{k=0}^{j-1}\binom{n}{k}p^{k}(1-p)^{n-k}$ (2)

This is the same as cumulative distribution function (CDF) of variable $Y$ at
$j-1$, $\text{CDF}(j-1,n,p)$.

If the samples of $X$ were arranged in ascending order, then the above
equation gives the confidence that the $p$’th percentile is equal to or
greater than the value at $j$’th place. In other words, this is the
probability that $p$ percentile values of $X$ are less than the value at
$j$’th place.

We can write a similar equation for some $i$, $0<=i<j$. Then the confidence
that $p$ lies between $i$ and $j$ is

$c=\text{CDF}(j-1,n,p)-\text{CDF}(i-1,n,p)=\sum_{k=i}^{j-1}\binom{n}{k}p^{k}(1-p)^{n-k},$
(3)

the equation from [1].

These computations were implemented in Python and have been published as an
open-source package [4]. Example usage in the form of Jupyter notebooks is
also available. Figure 1 shows the spread of confidence interval, $i$ to $j$
for median, obtained via setting $p=0.5$. Note that the numbers on the plot
for the size of spread are for the indices, not values. For example, if 99
samples of a random variable are obtained and arranged in ascending order, the
median is 50’th sample and the 90% confidence interval is from 50-8 = 42’nd
ordered sample to 50+8 = 58’th ordered sample.

Figure 1: Spread of confidence interval of median vs. number of samples at
different levels of confidence levels

To demonstrate the range of confidence interval around the median, 1000
samples of a random variable with normal distribution (mean = 0 and variance =
1) were generated. To find the 95% confidence interval, a subset of samples
was chosen, starting with first sample. The number of samples in the subset
were increased from 15 to 1000. For each subset, the confidence interval was
computed. Figure 2 subplot a shows the 95% confidence interval around the
median. For the same subsets, the mean with 95% confidence interval is shown
in subplot b. Both the subplots show that the range narrows as the number of
samples in the subset increase. For one example full set of 1000 samples, the
95% confidence interval range is 0.14 for median and 0.12 for mean. Obviously,
the size of range may be different for a different set of 1000 samples.

(a) Median with 95% confidence interval

(b) Mean with 95% confidence interval

Figure 2: Median and mean with 95% confidence intervals shown as shaded area
for a normal variable ($\mu=0$, $\sigma=1$) at different number of samples

Does this mean that the confidence interval of median is always larger than
that for mean? To compare the confidence intervals, the experiment was
repeated several times. In each experiment, a set of 1000 random numbers was
generated using the normal distribution. The values of 95% confidence interval
of mean and median were calculated. The experiment was repeated 100 times to
generate 100 such pairs and a kernel density estimate (KDE) plot was created.
The plot is shown in left-most sub-figure of Figure 3 with a equality line as
reference. The experiment was repeated for uniform and exponential
distributions and the results are shown in other sub-figures of the same.

(a) Normal

(b) Uniform

(c) Exponential

Figure 3: KDE plots of 95% confidence intervals of median against 95%
confidence intervals of mean for different distributions

##  3 Confidence in tolerance interval

Tolerance interval is the statistical range within which a certain fraction of
population samples are expected to fall. For example, 95% confidence interval
$(a,b)$ means $a<X_{i}<b$ for 95% values of $i$, $1<=i<=n$. If the
distribution of $X$ is known, one may be tempted to use well-known formulae.
For example, if $X$ has a normal distribution with mean $\mu$ and variance
$\sigma$, then 95% tolerance interval is $\mu\pm 1.96\sigma$. However, this
assumes that $\mu$ and $\sigma$ are known precisely. If these are instead
estimated from a sample of population, then the confidence in the calculated
range is low, only 50%.

For the normal distribution, it is possible to use 95% confidence levels of
mean and variance to compute the tolerance interval with 95% confidence. The
tolerance interval for other distributions may not be easy to compute.
Further, the type of distribution is usually only an assumption. Let’s see how
we can approach this problem without making any assumption about the type of
distribution.

Let $X_{i}$ be $i$’th sample of a random variable $X$, $1<=i<=n$. Assume that
these samples are independent and identically distributed. Their probability
distribution function need not be known. After the samples are sorted, then
the confidence or probability that the the $p$’th percentile (or quantile)
level, $0<p<1$, of these samples lies between $i$’th and $j$’th sample is [1]

$c=\sum_{k=i}^{j-1}\binom{n}{k}p^{k}(1-p)^{n-k}$ (4)

Note that $i$ and $j$ do not need to include the median index. The convention
for tolerance interval, though, is that it is the smallest central region of
the population within which the relevant percentile of the population samples
lie.

Define the $p$’th percentile ($0<p<1$) tolerance interval $(\pi_{g},\pi_{h})$
of a population sample of a random variable $X$ with a confidence $c$,
$0<c<1$, as the values of sorted $n$ samples of $X$ at indices $g$ and $h$
($1<=g<h<=n$), respectively, where $(g,h)$ is the smallest interval around the
median index such that

$c/2<=\sum_{k=1}^{h-1}\binom{n}{k}p^{k}(1-p)^{n-k},$ (5)

for the upper limit $h$ and

$1-c/2<=\sum_{k=1}^{g-1}\binom{n}{k}p^{k}(1-p)^{n-k}$ (6)

for the lower limit $g$.

In other words, no more than $1-c$ fraction of the samples outside of $(g,h)$
are spread outside the interval. Before we see illustrations of these
intervals, let’s introduce the concept of assurance interval.

##  4 Assurance interval for sorted samples

Quoting two numbers, percentiles and confidences, indeed gives more complete
view of the statistics. However, it may not make much sense when there is a
large mismatch between these two. For example, having a really low confidence,
say less than 50%, for a 95% tolerance interval is usually not sufficient.
Likewise, having high confidence, say 90%, in a smaller tolerance interval,
say 20%, may not be always worth it.

The concept of assurance was introduced in reliability engineering by Fulton
[3] by simply setting reliability equal to confidence. For example, assurance
of 90% means 90% reliability with 90% confidence. Building on this concept,
define assurance interval as an interval $(g,h)$ spanning central $a$
percentile values of $X$, with confidence level $a$. Mathematically,
substitute $c=a$ and $p=a$ in Equations 5 and 6 to get

$a/2<=\sum_{k=1}^{h-1}\binom{n}{k}a^{k}(1-a)^{n-k},$

for the upper limit $h$ and

$1-a/2<=\sum_{k=1}^{g-1}\binom{n}{k}a^{k}(1-a)^{n-k}$

for the lower limit $g$.

Given $n$, we can solve this equation for $a$ using numerical approach.
Brent’s optimization [5] method (function brentq in SciPy library) is used to
generate assurance intervals.

The above logic was implemented in Python package [4]. Figure 4 shows how
tolerance intervals expand with increasing percentiles and confidence levels.
It also compares assurance intervals with tolerance intervals.

Figure 4: Tolerance Intervals with confidence and Assurance Intervals for 200
sorted samples

##  5 Conclusions

In this article, the expression for confidence in percentiles for a population
was derived. Like the confidence interval of mean of a population, the
confidence interval of a median (or any percentile) decreases as the number of
samples increase.

The confidence interval of median was observed to be larger than that of mean
for normal and uniform distributions. Both these distributions are
symmetrical, so the mean and median are expected to be the same. For
exponential distribution, which is asymmetrical, the median is expected to be
smaller than mean. In the example shown, the confidence intervals of median
and mean were about the same. However, this relationship depends on the
parameters of the distributions.

We reviewed confidence in tolerance interval and defined an expression to
compute the central tolerance interval given desired percentiles, minimum
confidence, and number of samples. We observed that tolerance interval expands
as the percentile level increases at a constant confidence level. Tolerance
interval also expands as the confidence level increases at a constant
percentile level.

The concept of assurance from reliability engineering was extended to
tolerance interval and we defined the assurance interval with percentiles
matching minimum confidence. Numerical method using Brent’s optimization to
compute assurance interval was demonstrated. The tolerance intervals expand
only slightly as confidence increases, as opposed to increase in percentiles.
Assurance intervals match percentiles with confidence and are easier for
communication.

A Python open-source library to implement the expressions discussed in this
article was presented. A Jupyter notebook version to demonstrate sample usage
is included with the library.

## References

  * [1] Pennsylvania State University, Department of Statistics.  STAT 415: Introduction to Mathematical Statistics Lesson 19.1
  * [2] International Organization for Standardization (ISO)  _ISO 16269-6:2014 Statistical interpretation of data — Part 6: Determination of statistical tolerance intervals_ ,  Ed. 2, 2014
  * [3] S.N. Luko.  "Attribute Reliability and the Success Run: A Review",  DOI:10.4271/972753.
  * [4] S.M. Joshi.  _Reliable Statistics Package_ ,  https://github.com/sanjaymjoshi/relistats,  accessed Jan 1, 2024.
  * [5] R.P. Brent.  "Chapter 4: An Algorithm with Guaranteed Convergence for Finding a Zero of a Function",  _Algorithms for Minimization without Derivatives_ , Englewood Cliffs, NJ: Prentice-Hall, ISBN 0-13-022335-2, 1973


# Robust 360-8PA: Redesigning The Normalized 8-point Algorithm for 360-FoV
Images

Bolivar Solarte1
<EMAIL_ADDRESS>Chin-Hsuan Wu1
<EMAIL_ADDRESS>Kuan-Wei Lu1
<EMAIL_ADDRESS>Min Sun1
<EMAIL_ADDRESS>Wei-Chen Chiu2
<EMAIL_ADDRESS>Yi-Hsuan Tsai3
<EMAIL_ADDRESS>
###### Abstract

In this paper, we present a novel preconditioning strategy for the classic
8-point algorithm (8-PA) for estimating an essential matrix from 360-FoV
images (i.e., equirectangular images) in spherical projection. To alleviate
the effect of uneven key-feature distributions and outlier correspondences,
which can potentially decrease the accuracy of an essential matrix, our method
optimizes a non-rigid transformation to deform a spherical camera into a new
spatial domain, defining a new constraint and a more robust and accurate
solution for an essential matrix. Through several experiments using random
synthetic points, 360-FoV, and fish-eye images, we demonstrate that our
normalization can increase the camera pose accuracy about 20% without
significantly overhead the computation time. In addition, we present further
benefits of our method through both a constant weighted least-square
optimization that improves further the well known Gold Standard Method (GSM)
(i.e., the non-linear optimization by using epipolar errors); and a relaxation
of the number of RANSAC iterations, both showing that our normalization
outcomes a more reliable, robust, and accurate solution.

00footnotetext: 1 _National Tsing Hua University_ 00footnotetext: 2 _National
Chiao Tung University_ 00footnotetext: 3 _NEC Labs America_

##  1 Introduction

Estimating the relative pose between different views of the same scene has
been studied for decades in Computer Vision and Robotics. For instances,
Visual Odometry (VO), Simultaneous Localization and Mapping (SLAM), Structure
from Motion (SFM), among others, generally leverage this primary estimation
for initializing the first camera poses, triangulate landmarks in 3D, re-
localize the camera pose, and prune out outliers correspondences from the
system.

In practice, the goal of estimating a camera pose between two images relies on
finding a geometry constraint for their pixels. This is known as the essential
matrix or epipolar constraint. In general, the procedure of calculating an
essential matrix includes two main steps: 1) An abundant number of salient
pixels, i.e., key-features, are extracted from each image, followed by
matching them across different views; 2) Based on that correspondences, the
essential matrix is calculated satisfying the epipolar constraints. Finally,
after the essential matrix is derived, the relative camera pose can be
recovered by singular value decomposition [24, 8].

Figure 1: Illustration for the distribution of key-features in a 360-FoV
image. (a) and (b) present a 360-FoV image and its top view (point cloud)
respectively, with its key-features denoted as landmarks in green color. In
(c), the same landmarks presented in (a) and (b) are projected into a unit
sphere. In (d), our proposed normalization scheme (cf. Section 4.1) has been
applied to the spherical features in (c). Note that the key-features in the
ovoid surface (d) are geometrically shifted compared to (c) due to our
normalization. As a result, (d) expands the relative angles between every key-
feature, which in turn leads to a more stable solution when using the 8-PA
method [14].

Several algorithms have been proposed to find an essential matrix, however the
most widely used solutions are the five-point (5-PA) [20] and eight-point
(8-PA) [14] algorithms. Despite the former uses the minimum number of
correspondence points needed for calibrated cameras, its implementation
usually relies on a polynomial approximation with multiple solutions. In
contrast, the 8-PA is a linear method without the ambiguity in its outcome. In
general, the 8-PA is mostly used for 360-FoV images (e.g., [7, 18, 26, 4, 6,
27]), due to its simplicity and proved stability under large field of views
[2]. On the other hand, unlike the 5-PA, the 8-PA requires more iterations for
outlier removal using a RANSAC evaluation, hence increasing its computation
time for a large ratio of outliers [5]. This defines a clear disadvantage of
the 8-PA.

Although studies in [2] show that a wider FoV may increase the stability of
the 8-PA for spherical cameras, these assume that matched key-points in the
image are well-distributed in the whole FoV. However, in practice, that
distribution mainly depends on external factors, which sometimes yields to
clustered or uneven distribution of key-points (see Fig. 1) (e.g. [27]).

In this work, we improve the 8-PA [14] for spherical cameras by re-defining
the pre-conditioning strategy proposed by [9] to be applied to spherical
projection, which effectively deals with outliers and uneven key-feature
distributions for 360-FoV images. Additionally, we extend the usage of our
novel pre-conditioning by proposing a constant weighted least-square solution,
which improves the Gold Standard Method (GSM) [21, 8] that is usually used to
refine the camera pose. Lastly, we also present results comparing our solution
under a RANSAC evaluation, showing that our preconditioning is capable to
effectively deal with outliers, hence potentially reducing the number of
required iterations.

We evaluate our methods under both sequences of real 360-FoV and fish-eye
scenarios, where the former is our own dataset, collected from Matterport3D
[1] and rendered using MINOS [23], while the latter is the TUM-VI dataset
[25]. We show that our method significantly outperforms the state-of-the-art
8-PA for spherical cameras without the overhead in computation time,
demonstrating the robustness of our method against noise and outliers. In
favor of the research community, the source code is available at
https://github.com/EnriqueSolarte/robust_360_8PA.

##  2 Related Work

Several approaches have been developed to estimate an essential matrix based
on matched key-features. Nevertheless, the most well-known approaches are
still the 5-PA [20] and 8-PA [14], where the former has been improved largely
in past years, e.g., [10, 13, 12, 3]. However, all of them rely on a
polynomial approximation, which inevitably leads to multiple essential matrix
solutions. In contrast, the 8-PA [14] is a simpler and linear method that
always outputs a unique result.

In [9], a normalization strategy is introduced upon 8-PA [14], improving
homogeneity in the input data, and robustness against noise. However, this
preconditioning was mainly designed for uncalibrated pinhole cameras. Later,
in [19], the normalization [9] is further explained in terms of a generalized
total least square problem, which opens the idea of a general normalization.
However, in this work, the authors focus only on perspective cameras, proving
that the normalization [9] (i.e., an isotropic and non-isotropic
normalization), indeed, reduces the effects of noise in the solution of the
8-PA for key-point described in a homogeneous plane. In the literature there
is no evidence for normalization for spherical projection.

For spherical projection, the 8-PA has been widely used as a standard solution
(e.g., [17, 4, 6, 21, 29, 7]), that estimates a quick initial guess, where
other methods can efficiently refine it further [24, 8, 17]. Note that a
reliable and fast estimation of an initial guess is always desired.

Figure 2: Normalized 8-PA for spherical cameras. Our method takes two raw
360-FoV images as input, from where we extract several key-features and
matching their correspondences across the images as shown in (a). Next, in
(b), the same landmarks are projected into a unit sphere by spherical
projection (cf. 3-A). Afterward, the sphere is deformed by our normalization
method (cf. 4) in order to make a better distribution of key-features (c).
Lastly, we estimate the essential matrix by using DLT over the normalized
domain (cf. 4.1).

Several works have been proposed to explain the stability of the 8-PA as a
linear least-square problem, e.g., [32, 19]. However, most of them assume a
known noise distribution in the data. Instead, [2] recently demonstrates that
without any noise assumption, the 8-PA stability is highly related to the FoV
of the used image. However, the statement assumes a uniform distribution of
key-features along the FoV.

In general, a key-feature distribution is highly defined by external factors,
such as poor scene illumination, lack of texture, among others [8, 17, 22].
Indeed, for spherical cameras, due to the high distortion, uneven feature
distributions are a critical issue for matching correspondences as studied in
[27]. Therefore, using key-features directly from 360-FoV images may define an
uneven distribution.

##  3 Preliminaries

###  3.1 Spherical Projection and Bearing Vectors

Unlike perspective images where a homogeneous plane is used, 360-FoV images
are generally represented by a centralized spherical projection, which can be
described as:

$\begin{bmatrix}\theta\\\
\phi\end{bmatrix}=\begin{bmatrix}{2\pi}/{W}&&0&&-\pi\\\
0&&{-\pi}/{H}&&{\pi}/{2}\end{bmatrix}\begin{bmatrix}u\\\ v\\\ 1\end{bmatrix}$
(1) $\mathbf{q_{n}}=\begin{bmatrix}x_{n}\\\ y_{n}\\\
z_{n}\end{bmatrix}=\begin{bmatrix}cos(\phi)sin(\theta)\\\ -sin(\phi)\\\
cos(\phi)cos(\theta)\\\ \end{bmatrix},$ (2)

where (1) projects a pixel $(u,v)$ into the spherical coordinate
$(\theta,\phi)$, and then (2) transforms that coordinate into a unit vector
$\mathbf{q_{n}}$. Hereinafter, we name $\mathbf{q_{n}}$ as a bearing vector.
Note that $W$ and $H$ are the width and height of the 360-FoV image.

###  3.2 Epipolar Constraint and the Eight-Point Algorithm

Without loss of generality, the same epipolar constraint defined for
perspective images can be applied to spherical projection [24]. Therefore,
given a pair of 360-FoV images, two unit sphere cameras can be projected
(i.e., $C_{1}$ and $C_{2}$ in Fig. 2-(d)), from which the bearing vectors
$\mathbf{q_{1}}$ and $\mathbf{q_{2}}$ can be also defined. Then, the epipolar
constraint, which defines the coplanarity of $\mathbf{q_{1}}$ and
$\mathbf{q_{2}}$ onto an epipolar plane $\pi$, can be described as follows:

$\mathbf{q_{2}^{\top}Eq_{1}=0},\text{ with }\mathbf{E=[t~{}]_{\times}R}~{},$
(3)

where $\mathbf{E}$ represents the essential matrix $\in\mathbb{R}^{3\times 3}$
with rank of 2; $[\mathbf{t}~{}]_{\times}$ stands for a skew-symmetric matrix
coming from $\mathbf{t}\in\mathbb{R}^{3}$; and $\mathbf{R}$ defines a relative
camera rotation $\in SO(3)$.

We can then reformulate (3) into a Total Least Square Problem as follows:

$\mathbf{A}[\mathbf{E}]_{\textsc{v}}=0~{},$ (4)

where $[\mathbf{E}]_{\textsc{v}}$ is a vector of the coefficients in matrix
$\mathbf{E}$ by a row-wise concatenation; $\mathbf{A}$ is an $n\times 9$
matrix which is built upon stacking at least $n\geq 8$ correspondences of
bearing vectors (usually called observation matrix [8]), where an $i^{th}$ row
of this matrix can be defined by the Kronecker product of two correspondence
bearing vectors as
$\mathbf{A_{i}}=\mathbf{q^{i}_{1}}\otimes\mathbf{q^{i}_{2}}$.

Afterward, Singular Value Decomposition (SVD) is applied to $\mathbf{A}$ for
finding the solution of (4), in which the last column of the orthogonal
subspace of $\mathbf{A}$ defines $[\mathbf{E}]_{\textsc{v}}$, and can then be
reshaped into $\mathbb{R}^{3\times 3}$ and forced into rank(2). This procedure
is called Direct Linear Transformation (DLT) [8].

##  4 Robust Spherical Normalization

As the key-features are primarily located at areas with high texture, a wider
FoV image cannot always provide an uniformly distributed locations of key-
features. Thus, we propose our normalization strategy for the classic 8-PA to
reduce the effect when the key-features are not uniformly distributed.

First, we derive the mechanism introduced in [9] but with our normalization
applied, and then we define a new vector space for stabler DLT solution of an
essential matrix. We also discuss the reasons why our normalization increases
the stability of the 8-PA. Lastly, we present two non-linear optimizations
aiming to improve both the 8-PA and GSM methods without overhead computation
time.

###  4.1 Normalization

Unlike [9], we define our normalization as a transformation that deforms
bearing vectors from an unit sphere camera into a ovoid surfaces (see Fig.
3-(c)). This matrix transformation is defined as follows:

$\displaystyle\mathbf{N}=\begin{bmatrix}S&&0&&0\\\ 0&&S&&0\\\
0&&0&&K\end{bmatrix}~{},$ (5)
$\displaystyle\mathbf{\hat{q}}_{i}=\mathbf{N}~{}\mathbf{{q}}_{i}~{},$ (6)

where $S,K\in\mathbb{R}$ and $|\mathbf{N}|\neq 0$. For convenience, our
normalization is designed as a diagonal matrix $\mathbb{R}^{3\times 3}$,
controlled by two parameters, $S$ and $K$, as presented in (5), which allows
us to deform bearing vectors along XY and Z directions, respectively. Thus, we
represent the normalization (6) with the epipolar constraint presented in (3):

$\displaystyle\mathbf{\hat{q}}_{2}^{\top}~{}\mathbf{N}_{2}^{-\top}~{}\mathbf{E}~{}\mathbf{N}_{1}^{-1}~{}\mathbf{\hat{q}}_{1}=0~{}.$
(7)

By further arranging (7), we can build our normalized constraint as follows:

$\displaystyle\mathbf{\hat{q}}_{2}^{\top}\hat{\mathbf{E}}~{}\mathbf{\hat{q}}_{1}=0~{},$
(8)
$\displaystyle\mathbf{E}=\mathbf{N}^{\top}\hat{\mathbf{E}}~{}\mathbf{N}~{}.$
(9)

Note that (8) is embodied by $\mathbf{\hat{E}}$, which stands for our
normalized essential matrix in the normalized domain; thus, we can find
$\mathbf{\hat{E}}$ by using the standard DLT procedure described in Sec. 3. In
this way, we can recover the original essential matrix $\mathbf{E}$ by left-
right multiplying the matrix $\mathbf{N}$ as described in (9), which we call
denormalization. Finally, the relative camera pose $\mathbf{T}\in SE(3)$ can
be recovered from $\mathbf{E}$ by SVD. To be clear, we present the procedure
in Fig. 2.

###  4.2 Deformation of Spherical Projection

Intuitively, normalizing bearing vectors through a non-rigid transformation
(5) deforms the unit sphere into a different spatial domain which in turns
defines a new observation matrix $\mathbf{A}$. To explain why this
normalization leads to more stable DLT solution for (8), we present two
properties that define stability for 8-PA in the context of spherical cameras.

Figure 3:  Stability of the spherical projection. In panel (a) given the
locations of the landmarks $\textbf{P}_{i}$ and $\textbf{P}_{j}$ the internal
angles $\beta^{1}_{ij}$ and $\beta^{2}_{ij}$ as well as the motion parallax
$\alpha_{i}$ and $\alpha_{j}$ are defined. In panels (b)-(d), geometric
visualizations of our normalization are presented (cf. Sec.4.1). Panel (b)
shows a spherical camera without normalization as a reference. In (c), the
effects of setting $S=2$ are illustrated, while in (d) the effects of $K=2$.
Note that normalizing a spherical camera using $(S,~{}K)$ can increase the
motion parallax in the data. Lastly, in panel (e) we present a surface loss
for translation error build upon several combinations of $S$ and $K$ values
from a grid. We can verify that there exists a combination of $(S,~{}K)$ that
reduces error.

First, based on [2], we assert that the FoV of an image has a strong
correlation with the stability of an essential matrix estimation when using
the 8-PA. The main reason is because large FoV images are prone to define
large internal angles between its bearing vectors, i.e. the angles
$\beta^{1}_{ij}$ and $\beta^{2}_{ij}$ in Fig. 3(a). This can be mathematically
justified by representing these internal angles in terms of the observation
matrix $\mathbf{A}$ as shown in (10), which in turn evaluates a bound for the
second least singular value of $\mathbf{A}$ (i.e. $\sigma_{8}$) as presented
in (11). For more details, we refer to Sec. 3.4. in [2]:

$\displaystyle\left\|\mathbf{A}\mathbf{A}^{\top}\right\|^{2}_{F}=\sum_{i}\sum_{j}{(cos^{2}_{i}\beta^{1}_{j})(cos^{2}_{i}\beta^{2}_{j})}~{},$
(10)
$\displaystyle\sigma_{8}\leq\sqrt{\frac{n}{8}-\frac{1}{8}\sqrt{\frac{8\left\|\mathbf{A}\mathbf{A}^{\top}\right\|^{2}_{F}-n^{2}}{7}}}~{}.$
(11)

Based on [2, 32], the error in an essential matrix $\mathbf{E}$ can be defined
as a function of $\sigma_{8}$ as follows:

$\displaystyle\Delta\mathbf{E}\leq\frac{\left\|\mathbf{Q}\right\|_{2}}{\sigma_{8}}~{},$
(12)

where $\Delta\mathbf{E}$ is the error in the essential matrix, and
$\mathbf{Q}$ represents the perturbation matrix which embodies the noise and
outliers in bearing vectors. In practice, based on (12), we assert that larger
internal angles, $\beta^{1}_{ij}$ and $\beta^{2}_{ij}$, can evaluate larger
$\sigma_{8}$ values, which in turn result in a lower $\Delta\mathbf{E}$ error.
Therefore, if we deform the spherical projection by increasing the internal
angles between bearing vectors, we will have larger $\sigma_{8}$ values and
thus obtain a more stable solution. In addition, based on the translation
vector $\mathbf{t}$ between cameras and the distance to some landmarks in the
scene, we can define the motion parallax for the spherical projection as the
angles $\alpha_{j}$ and $\alpha_{i}$ in Fig. 3(a). Therefore, analogously to
the motion parallax defined for perspective projection (e.g. [8, 17]), we can
assert that when the angles $\alpha_{j}$ and $\alpha_{i}$ are close to zero,
the DLT solution of the 8-PA is unable to recover a camera pose. Thus, the
lack of motion parallax is one of the degenerative conditions for the 8-PA
(e.g. [14, 9, 17]). Here, if we properly dislocate every bearing vector in a
particular direction, we can modify the motion parallax and further define a
more reliable estimation of a camera pose.

As illustrated in Fig. 3(c), we compute a grid of $S$ and $K$ values defining
different normalizing matrices (5). By using random synthetic landmarks in 3D,
we apply our proposed normalization to two known spherical cameras as
described in Sec. 4.1. Thus, we can visualize that there exists a set of $S$
and $K$ values which improves a camera pose estimation by normalizing bearing
vectors.

In the same synthetic environment, we evaluate the motion parallax for every
matched bearing vector and plot it in a scale of color in Fig. 3(b)-(d). By
deforming the spherical projection along a specific direction, we can magnify
the motion parallax in a set of bearing vectors.

###  4.3 Non-linear Optimization over S and K

Although the normalized 8-PA proposed by [9] gives us a mechanism to apply our
normalization (5) to spherical features, the method does not provide a valid
procedure to evaluate a normalized matrix under spherical projection.

Therefore, we propose a non-linear optimization to reduce errors in the
epipolar constraint by locally finding an optimal matrix $\mathbf{N}$,
parameterized by $S$ and $K$, which effectively normalizes bearing vectors in
spherical projection.

In practice, based on the projected distance in [21], we define our residual
error in the epipolar constraint as follows:

$\displaystyle\epsilon(\Theta)=\frac{|\mathbf{q}_{2}^{\top}~{}\mathbf{E}_{\Theta}~{}\mathbf{q}_{1}|}{\left\|\mathbf{q}_{2}\right\|\left\|\mathbf{E}_{\Theta}~{}\mathbf{q}_{1}\right\|}~{},$
(13)

where ${\mathbf{E}}_{\Theta}$ represents the essential matrix evaluated by a
particular set of $\Theta$ parameters. Thus, for our proposed method, $\Theta$
represents $S$ and $K$ while as for the optimization of GSM, it is defined as
$\xi\in\mathbb{R}^{6}$, i.e. $\mathbf{T}\in SE(3)$ mapped into the Lie group
$\xi\in\mathfrak{se}(3)$ [17]. Thus, we can define our non-linear optimization
to find a set of $S^{*}$ and $K^{*}$ parameters as follows:

$\displaystyle
S^{*},K^{*}=\operatorname*{argmin}_{S,K\in\mathbb{R}^{2}}\left\|\epsilon(S,K)\right\|^{2}_{1}.$
(14)

Unlike GSM, which is defined over 6-DoF while minimizing the least-square
errors of (13), our proposed optimization (14) is defined over 2-DoF only
(i.e. $S$ and $K$), which do not evaluate a camera pose directly; thus an
initial evaluation of the 8-PA is not needed for an initial guess, which in
turns does not add the overhead in time. To deploy our optimization (14), we
leverage the LM algorithm [15] with both $S$ and $K$ parameters set to 1 as a
trivial initial point.

Figure 4:  Synthetic Points Evaluation. Unless noted otherwise, all of the
experiments described in this figure uses: 200 matched/correspondences bearing
vectors in spherical projection as input; von Misses-Fisher noise of
$\kappa=500$ (3.21$\degree$ of error). In panel (a), a constant outliers ratio
of 20% has been added, while in (b), the input data is evaluated as free-
outliers data.

To further show the versatility of our normalization upon the 8-PA solution
for 360∘ images, we also propose a robust constant weighting function to
improve the GSM accuracy without increasing its computation time. This
optimization is described as follows:

$\displaystyle\xi^{*}=\operatorname*{argmin}_{\xi\in\mathbb{R}^{6}}\sum_{i}\omega_{i{SK}}~{}\epsilon_{i}(\xi)^{2}~{},$
(15)

where $\omega_{SK}$ is a constant vector build upon the residuals
$\epsilon(S^{*},K^{*})$, which defines the confidence for every matched
bearing vector as a probabilistic function $P(\epsilon(S^{*},K^{*}))$.
Combined with our proposed optimization (14), we compute a robust essential
matrix $\mathbf{E}_{SK}$ and evaluate a residuals vector
$\epsilon(S^{*},K^{*})$, from where $\omega_{SK}$ is computed as a normal
distribution $\mathcal{N}(\epsilon(S^{*},K^{*})|\mu,\sigma)$, evaluating $\mu$
and $\sigma$ as the mean and standard deviation of the residuals vector
$\epsilon(S^{*},K^{*})$, respectively. To deploy the optimization (15), the LM
algorithm [15] is also used.

##  5 Experimental Results

In this section, we conduct experiments to demonstrate that our method
achieves a more accurate camera pose estimation for 360-FoV images. We define
two environments: synthetic random points (Sec. 5.1) and sequence of real-
images (Sec. 5.2 and 5.3).

For experiments in Sec. 5.3, we use both 360-FoV images in spherical
projection and fish-eye images in the unified camera model [31]. For the
former, we collect our own dataset (named 360-MP3D-VO) by rendering the
Matterport3D data [1] on Minos [23] to project a continues sequence of 360-FoV
frames; for the latter, we use the TUM-VI dataset [25].

The evaluated baselines are the classic 8-PA [14] in spherical projection and
the GSM solution [21]. For the former, we compare with our normalized approach
using the optimization (14), which we refer to as $Opt~{}\epsilon(S,K)$. To
compare with GSM, we use our proposed optimization (15), where we refer to
$\omega_{\xi}$ and $\omega_{SK}$ as a weighting function evaluated by
residuals from a camera pose $\xi$ (computed by using the 8-PA[9]) and
residuals obtained from our normalization $Opt~{}\epsilon(S,K)$, respectively.

Similar to [3, 2, 7], we use $\epsilon_{R}$ as the rotation error and
$\epsilon_{t}$ as the translation error for camera pose estimation:

$\epsilon_{R}=\frac{1}{\pi}~{}cos^{-1}\left(\frac{tr(\mathbf{R}^{\top}\mathbf{\tilde{R}})-1}{2}\right),\quad\epsilon_{\mathbf{t}}=\frac{1}{\pi}~{}cos^{-1}(\mathbf{t}^{\top}\mathbf{\tilde{t}})~{}.$
(16)

###  5.1 Noise and Outlier Evaluation

In this experiment, we project several landmarks in 3D by using the ground
truth depth data from our 360-MP3D-VO dataset. To generate ground truth camera
poses, we randomly sample translation vectors $\mathbf{t}$ in a uniform range
of $[-1,1]\in\mathbb{R}^{3}$. Moreover, for each axis, we generate random
rotation matrices $\mathbf{R}\in SO(3)$ by sampling Euler-angles in a uniform
range of $[-\pi/4,\pi/4]$. Then, we construct our camera poses as homogeneous
transformations $\mathbf{T}=[\mathbf{R}|\mathbf{t}]\in SE(3)$.

To show the effect of using different numbers of point, we uniformly sample
between 8 to 200 3D-landmarks for each evaluation. Then, similar to [7, 2], we
add a constant von Misses-Fisher noise (vMF) of $\kappa=500$ (i.e.,
3.21$\degree$ of error) as well as a constant outlier ratio of 20%. The
results of 15K evaluations are presented in Fig. 4-(a). In addition, to
evaluate camera pose under different levels of noise, we sample 200
3D-landmarks of 15K different scenes, and then we incrementally add a vMF
noise defined by a $\kappa$ equals to 100, 200, 500, 1000, and 10000,
representing 10.21$\degree$, 5.22$\degree$, 3.21$\degree$, 2.27$\degree$,
1.60$\degree$ and 0.72$\degree$ of error, respectively. This experiment is
presented in Fig. 4-(b). To further evaluate our method against outliers, we
use a constant number of 200 3D-landmarks and a vMF noise of $\kappa=500$,
then we increase outliers ratio from 0% to 70%. This is shown in Fig. 4-(c).

In Fig. 4, we show favorable results against others in terms of noise level,
outlier level, and number of points. For instance, our normalization is
capable of reducing the effect of outliers around 10% compared with the
baseline, i.e., using $\omega_{\xi}$ without our normalization to build
$\omega$.

###  5.2 Ablation Study

Table 1: Comparisons with 8-PA methods.

$\sigma_{8}$$\shortuparrow$ $\alpha$$\shortuparrow$ Rot-e $\times 10^{-3}$
Tran-e $\times 10^{-3}$ Time (ms) $Opt~{}\epsilon(S,K)$ 20.025 0.858 11.429
25.416 45.52 8-PA [14] 0.650 0.731 13.387 32.367 40.33 Isotropic (a) [9] 1.232
0.723 15.937 49.414 40.33 Non-Isotropic (b) [9] 0.012 0.723 13.439 34.354
40.33 (a) + (b) [19] 0.107 0.879 14.173 38.808 40.88

Table 2: Ablation study.

Rot-e $\times 10^{-3}$ Tran-e $\times 10^{-3}$ Time (ms) GSM [21] (a) 9.1029
19.7902 45.318 Gaussian (a) + not const $\omega_{\xi}$ 2.0951 5.2065 218.598
(a) + not const $\omega_{SK}$ 2.0943 5.1894 212.924 (a) + $\omega_{\xi}$
3.9048 9.2592 50.279 (a) + $\omega_{SK}$ 3.6802 7.5717 58.651 t-distribution
(a) + not const $\omega_{\xi}$ 9.1810 19.8160 48.126 (a) + not const
$\omega_{SK}$ 8.9496 18.9284 58.836 (a) + $\omega_{\xi}$ 9.1810 19.8160 50.279
(a) + $\omega_{SK}$ 9.0581 19.5573 53.795

Table 3: Experimental results in real scenes on MP3D-VO and TUM-VI.

Rotation error $\times 10^{-3}$ Translation error $\times 10^{-3}$ 8-PA [14]
$Opt~{}\epsilon(S,K)$ GSM [21] $\omega_{\xi}$ $\omega_{SK}$ 8-PA [14]
$Opt~{}\epsilon(S,K)$ GSM [21] $\omega_{\xi}$ $\omega_{SK}$ MP3D-VO Q75 23.577
20.504 15.286 9.909 7.964 91.442 72.693 43.707 29.945 20.653 Q50 10.827 9.515
7.389 3.523 3.072 36.814 29.645 19.453 9.985 7.966 Q25 4.624 4.075 3.277 1.638
1.524 14.745 12.276 8.501 4.611 4.009 TUM-VI [25] Q75 44.479 40.514 49.291
39.938 38.256 221.131 199.042 166.594 169.122 140.846 Q50 27.211 24.969 27.226
21.075 19.241 117.701 104.742 82.457 74.661 62.743 Q25 15.512 14.462 14.595
10.507 9.708 58.791 53.283 40.899 33.724 28.800

Table 4: Experimental results under different thresholds with RANSAC.

Threshold: 2.30E-04 Iterations: 590 Outliers in data: 1% Threshold: 1.10E-03
Iterations: 66 Outliers in data: 20% 8-PA [14] $Opt~{}\epsilon(S,K)$ GSM [21]
$\omega_{\xi}$ $\omega_{SK}$ 8-PA [14] $Opt~{}\epsilon(S,K)$ GSM [21]
$\omega_{\xi}$ $\omega_{SK}$ Rot-e $\times 10^{-3}$ 3.957 3.950 3.691 2.994
2.919 10.658 9.963 9.978 4.292 4.065 Tran-e $\times 10^{-3}$ 8.757 8.689 8.239
7.660 7.477 22.656 19.806 20.394 10.540 9.452 Residual 2.72E-06 2.72E-06
2.88E-06 2.94E-06 2.94E-06 2.72E-05 2.74E-05 2.77E-05 2.91E-05 2.92E-05 Time
(sec) 0.149 0.153 0.181 0.183 0.189 0.015 0.021 0.050 0.065 0.055

In this experiment, we show results of our solution by using tracked key-
features from a sequence of 360-FoV images in a 6-DoF camera motion. For each
frame in this sequence, we extract around 200 key-features using Shi-Tomasi
key-points [28] and the KLT tracker [16]. For every evaluation, we ensure that
there exits at least a minimum distance of 0.5m between frames. We evaluate
over 2K different environments and compute the median values of the estimated
errors.

In Table 1, we compare the proposed normalized solution with the 8-PA
algorithms [14, 9, 19], all of them in spherical projection. Note that
isotropic and non-isotropic pre-conditioning are the normalization strategies
proposed by [9], which are particularly used for uncalibrated cameras in
perspective view, but can still be used for spherical projection. In the
results, our normalization is capable of reducing errors in camera pose
without significantly adding the overhead time. Note that, our normalization
obtains larger values of $\sigma_{8}$ (second least singular value of an
observation matrix $\mathbf{A}$) as well as $\alpha$ (motion parallax in the
normalizing domain), showing that our normalization truly increases the
stability in the DLT solution of an essential matrix.

In Table 2, we compare the effect of our normalization in the weighted non-
linear optimization (15) upon the GSM solution [21]. In rows 2-5, we find that
using a Gaussian distribution constantly performs better than t-distributions.
In rows 2,3 and 6,7, we evaluate a non-constant weighing function $\omega$
which is updated at every iteration in the LM optimization [30, 11, 7]; this
approach is also known as Iterative Re-weighted Least-Square method (IRLS).
Although IRLS achieves the lowest error, it is still the most time-consuming
method among all evaluations.

###  5.3 Experiments on Real Scenes

We evaluate a sequence of 360-FoV and fish-eye images with our own 360-MP3D-VO
and the TUM-VI [25] datasets, respectively. Similar to Sec. 5.2, we track 200
key-features between frames, evaluating around 15K pairs of frames using
360-MP3D-VO, and 16K paired images using TUM-VI.

For the evaluations on the TUM-VI [25] dataset, we only use the scenes that
contain a complete ground truth camera pose, i.e., room-1 to room-6. Moreover,
since this dataset is mainly used for visual-inertial tasks, some frames in
our evaluation has been skipped due to drastic camera movements and severe
changes in illumination.

Results of camera pose errors on both datasets are presented in Table 3, using
the quantiles Q25, Q50 and Q75 of error evaluations. From the results, we can
verify that our approaches constantly outperform the baselines, with the
lowest errors in every evaluation. Moreover, for the averaged condition (i.e.,
Q50 results), our strategies are capable of reducing errors up to 10% by
normalizing bearing vectors, using $Opt~{}\epsilon(S,K)$; and up to 50% by
using our weighted optimization $\omega_{SK}$.

Additionally, we evaluate our proposed approaches $Opt~{}\epsilon(S,K)$ and
$\omega_{SK}$ in the context of a RANSAC evaluation. We design two settings
using the same amount of correspondence bearing vectors, noise, and outliers
(e.g., 400 3D-landmark, vMF noise of $\kappa=500$, and 50% of outliers), but
with two different thresholds for RANSAC. In addition, upon the RANSAC
results, we evaluate the final essential matrix using only the detected
inliers by using our proposed methods as well as the baselines. In Table 4,
results over 2K evaluations are presented.

On the left side of Table 4, we set a threshold for RANSAC as $2.3\times
10^{-4}$, which successfully detects all the inliers after 590 iterations;
whereas on the right side, the previous threshold is relaxed until $1.1\times
10^{-3}$, speeding up the evaluation by detecting 70% of inliers in 66 RANSAC
iterations. However, since we know in advance the outlier ratio of our data,
we can assert that there is 20% of the remaining outliers. Therefore,
comparing columns 2 and 11, we verify that our proposed $\omega_{SK}$ performs
similar to the one using RANSAC with a tight threshold, but 3 times faster,
showing the benefit of our robust estimation.

##  6 Conclusions

In this paper, we propose a novel pre-conditioning strategy to the classic
8-point algorithm [14] for estimating an essential matrix in spherical
projection. Our solution redesigns the normalizing algorithm [9] to alleviate
the effect of poor/uneven distribution of key-features, increasing stability
and robustness against outliers. We also extend our approach, towards
improving the well-known Gold Standard Method [21] for spherical projection by
proposing a constant weighted non-linear optimization, built upon our
normalization strategy. Based on extensive experiments under different scene
conditions, our proposed methods outperform the baselines, increasing the
accuracy in camera pose up to 30% without a significant impact on the
computation time.

Acknowledgement.

This project is supported by MOST Joint Research Center for AI Technology and
All Vista Healthcare with program MOST 110-2634-F-007-016 and MOST
110-2636-E-009-001.

## References

  * [1] Angel Chang, Angela Dai, Thomas Funkhouser, Maciej Halber, Matthias Niessner, Manolis Savva, Shuran Song, Andy Zeng, and Yinda Zhang.  Matterport3D: Learning from RGB-D data in indoor environments.  In International Conference on 3D Vision (3DV), 2017.
  * [2] Thiago LT da Silveira and Claudio R Jung.  Perturbation analysis of the 8-point algorithm: a case study for wide fov cameras.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 11757–11766, 2019.
  * [3] Kaveh Fathian, J Pablo Ramirez-paredes, Emily A Doucette, J Willard Curtis, and Nicholas R Gans.  QuEst: A Quaternion-Based Approach for Camera.  IEEE Robotics and Automation Letters (RA-L), 2018.
  * [4] Cornelia Fermüller and Yiannis Aloimonos.  Ambiguity in structure from motion: Sphere versus plane.  International Journal of Computer Vision (IJCV), 1998.
  * [5] Friedrich Fraundorfer and Davide Scaramuzza.  Visual odometry: Part ii: Matching, robustness, optimization, and applications.  IEEE Robotics and Automation Magazine (RA-M), 2012.
  * [6] Jun Fujiki, Akihiko Torii, and Shotaro Akaho.  Epipolar geometry via rectification of spherical images.  In International Conference on Computer Vision/Computer Graphics Collaboration Techniques and Applications, 2007.
  * [7] H. Guan and W. A. P. Smith.  Structure-from-motion in spherical video using the von mises-fisher distribution.  IEEE Transactions on Image Processing, 26(2):711–723, 2017.
  * [8] Richard Hartley and Andrew Zisserman.  Multiple view geometry in computer vision.  Cambridge university press, 2003.
  * [9] Richard I Hartley.  In defense of the eight-point algorithm.  IEEE Transactions on Pattern Analysis and Machine Intelligence (TPAMI), 1997.
  * [10] Hongdong Li and R. Hartley.  Five-point motion estimation made easy.  In 18th International Conference on Pattern Recognition (ICPR’06), volume 1, pages 630–633, 2006.
  * [11] Christian Kerl, Jürgen Sturm, and Daniel Cremers.  Robust odometry estimation for rgb-d cameras.  In 2013 IEEE International Conference on Robotics and Automation, pages 3748–3754. IEEE, 2013.
  * [12] Bo Li, Lionel Heng, Gim Hee Lee, and Marc Pollefeys.  A 4-point algorithm for relative pose estimation of a calibrated camera with a known relative rotation angle.  In IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2013.
  * [13] Hongdong Li and Richard Hartley.  Five-point motion estimation made easy.  In International Conference on Pattern Recognition (ICPR), 2006\.
  * [14] H Christopher Longuet-Higgins.  A computer algorithm for reconstructing a scene from two projections.  Nature, 1981.
  * [15] M.I.A. Lourakis.  levmar: Levenberg-marquardt nonlinear least squares algorithms in C/C++.  [web page] `http://www.ics.forth.gr/~lourakis/levmar/`, Jul. 2004\.  [Accessed on 31 Jan. 2005.].
  * [16] Bruce D Lucas, Takeo Kanade, et al.  An iterative image registration technique with an application to stereo vision.  International Joint Conferences on Artificial Intelligence, 1981\.
  * [17] Yi Ma, Stefano Soatto, Jana Kosecka, and S Shankar Sastry.  An invitation to 3-d vision: from images to geometric models, volume 26.  Springer Science & Business Media, 2012.
  * [18] Pierre Moulon, Pascal Monasse, Romuald Perrot, and Renaud Marlet.  Openmvg: Open multiple view geometry.  In International Workshop on Reproducible Research in Pattern Recognition, 2016.
  * [19] Matthias Mühlich and Rudolf Mester.  The role of total least squares in motion analysis.  In European Conference on Computer Vision, pages 305–321. Springer, 1998.
  * [20] David Nistér.  An efficient solution to the five-point relative pose problem.  IEEE Transactions on Pattern Analysis and Machine Intelligence (TPAMI), 2004.
  * [21] A. Pagani and D. Stricker.  Structure from motion using full spherical panoramic cameras.  In 2011 IEEE International Conference on Computer Vision Workshops (ICCV Workshops), pages 375–382, 2011.
  * [22] Ethan Rublee, Vincent Rabaud, Kurt Konolige, and Gary Bradski.  Orb: An efficient alternative to sift or surf.  In IEEE International Conference on Computer Vision (ICCV), 2011\.
  * [23] Manolis Savva, Angel X. Chang, Alexey Dosovitskiy, Thomas Funkhouser, and Vladlen Koltun.  MINOS: Multimodal indoor simulator for navigation in complex environments.  ArXiv:1712.03931, 2017.
  * [24] Davide Scaramuzza and Friedrich Fraundorfer.  Tutorial: visual odometry.  IEEE Robotics and Automation Magazine (RA-M), 2011.
  * [25] D. Schubert, T. Goll, N. Demmel, V. Usenko, J. Stueckler, and D. Cremers.  The tum vi benchmark for evaluating visual-inertial odometry.  In International Conference on Intelligent Robots and Systems (IROS), October 2018.
  * [26] Shinya Sumikura, Mikiya Shibuya, and Ken Sakurada.  OpenVSLAM: A Versatile Visual SLAM Framework.  In ACM Conference on Multimedia (MM), 2019.
  * [27] Hajime Taira, Yuki Inoue, Akihiko Torii, and Masatoshi Okutomi.  Robust feature matching for distorted projection by spherical cameras.  IPSJ Transactions on Computer Vision and Applications, 2015.
  * [28] Tiziano Tommasini, Andrea Fusiello, Emanuele Trucco, and Vito Roberto.  Making good features track better.  In Proceedings. 1998 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (Cat. No. 98CB36231), pages 178–183. IEEE, 1998.
  * [29] Akihiko Torii, Atsushi Imiya, and Naoya Ohnishi.  Two-and three-view geometry for spherical cameras.  In Proceedings of the sixth workshop on omnidirectional vision, camera networks and non-classical cameras, 2005.
  * [30] Philip HS Torr and David William Murray.  The development and comparison of robust methods for estimating the fundamental matrix.  International journal of computer vision, 24(3):271–300, 1997.
  * [31] Vladyslav Usenko, Nikolaus Demmel, and Daniel Cremers.  The double sphere camera model.  In 2018 International Conference on 3D Vision (3DV), pages 552–560. IEEE, 2018.
  * [32] Per-Åke Wedin.  Perturbation bounds in connection with singular value decomposition.  BIT Numerical Mathematics, 12(1):99–111, 1972.