Florian Bugarin

Adrien Bartoli

Didier Henrion

Jean-Bernard Lasserre

Jean-José Orteu

Thierry Sentenac

Rank-Constrained Fundamental Matrix Estimation by Polynomial Global Optimization Versus the Eight-Point Algorithm

Keywords: Global optimization, linear matrix inequality, fundamental matrix

The fundamental matrix can be estimated from point matches. The current gold standard is to bootstrap the eight-point algorithm and two-view projective bundle adjustment. The eight-point algorithm first computes a simple linear least squares solution by minimizing an algebraic cost and then projects the result to the closest rank-deficient matrix. We propose a single-step method that solves both steps of the eight-point algorithm. Using recent results from polynomial global optimization, our method finds the rank-deficient matrix that exactly minimizes the algebraic cost. In this special case, the optimization method is reduced to the resolution of very short sequences of convex linear problems which are computationally efficient and numerically stable. The current gold standard is known to be extremely effective but is nonetheless outperformed by our rank-constrained method for bootstrapping bundle adjustment. This is here demonstrated on simulated and standard real datasets. With our initialization, bundle adjustment consistently finds a better local minimum (achieves a lower reprojection error) and takes less iterations to converge.

Introduction

The fundamental matrix has received a great interest in the computer vision community (see for instance [START_REF] Longuet-Higgins | A computer algorithm for reconstructing a scene from two projections[END_REF][START_REF] Armangué | Overall view regarding fundamental matrix estimation[END_REF][START_REF] Tsai | Uniqueness and estimation of three-dimensional motion parameters of rigid objects with curved surfaces[END_REF][START_REF] Hartley | In Defence of the 8-point Algorithm[END_REF][START_REF] Torr | Invariant fitting of two view geometry or[END_REF][START_REF] Chojnacki | Revisiting hartley's normalized eight-point algorithm[END_REF][START_REF] Bartoli | Nonlinear estimation of the fundamental matrix with minimal parameters[END_REF]). This (3 × 3) rank-two matrix encapsulates the epipolar geometry, the projective motion between two uncalibrated perspective cameras, and serves as a basis for 3D reconstruction, motion segmentation and camera self-calibration, to name a few. Given n point matches (q i , q ′ i ), i = 1, . . . , n between two images, the fundamental matrix may be estimated in two phases. The initialization phase finds some suboptimal estimate while the refinement phase iteratively minimizes an optimal but nonlinear and nonconvex criterion. The gold standard uses the eight-point algorithm and projective bundle adjustment for these two phases, respectively. A 'good enough' initialization is necessary to avoid local minima at the refinement phase as much as possible. The main goal of this article is to improve the current state of the art regarding the initialization phase. We here focus on input point matches that do not contain mismatches (a pair of points incorrectly associated). The problem of mismatches has been specifically addressed by the use of robust methods in the literature.

The eight-point algorithm follows two steps [START_REF] Longuet-Higgins | A computer algorithm for reconstructing a scene from two projections[END_REF]. In its first step, it relaxes the rank-deficiency constraint and solves the following convex problem:

F = arg min F∈R 3×3 C(F) s.t. F 2 = 1, ( 1
)
where C is a convex, linear least squares cost, hereinafter called the algebraic cost:

C(F) = n ∑ i=1 q ′⊤ i Fq i 2 . ( 2
)
This minimization is subject to the normalization constraint F 2 = 1. This is to avoid the trivial solution F = 0. Normalization will be further discussed in section 3. The estimated matrix F is thus not a fundamental matrix yet. In its second step, the eight-point algorithm computes the closest rank-deficient matrix to F as:

F 8pt = arg min F∈R 3×3 F -F 2 s.t. det(F) = 0. (3)
Both steps can be easily solved. The first step is a simple linear least squares problem and the second step is solved by nullifying the least singular value of F. It has been shown [START_REF] Hartley | In Defence of the 8-point Algorithm[END_REF] that this simple algorithm performs extremely well in practice, provided that the image point coordinates are standardized by simply rescaling them so that they lie in [-√ 2;

√ 2] 2 . Our main contribution in this paper is an approach that solves for the fundamental matrix minimizing the algebraic cost. In other words, we find the global minimum of:

F = arg min F∈R 3×3 C(F) s.t. det(F) = 0 and F 2 = 1. ( 4
)
Perhaps more importantly, we also quantify the impact that each of F 8pt and F has when used as an initial estimate in Bundle Adjustment. Each initial estimate will lead Bundle Adjustment to its own refined estimate. The two final estimates may thus be different since, as the difference between the two initial estimates grows larger, the probability that they lie in different basins of attraction increases. Our measure quantifies:

1. how far are these two basins of attractions, 2. how many iterations will Bundle Adjustment take to converge.

The proposed algorithm uses polynomial global optimization [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF][START_REF] Henrion | Gloptipoly: Global optimization over polynomials with Matlab and SeDuMi[END_REF]. Previous attempts [START_REF] Xuelian | New fundamental matrix estimation method using global optimization[END_REF][START_REF] Chesi | Estimating the fundamental matrix via constrained least-squares: A convex approach[END_REF][START_REF] Kahl | Globally optimal estimates for geometric reconstruction problems[END_REF] in the literature differ in terms of optimization strategy and parameterization of the fundamental matrix. None solves problem (4) optimally for a general parameterization: they either do not guarantee global optimality [START_REF] Chesi | Estimating the fundamental matrix via constrained least-squares: A convex approach[END_REF][START_REF] Kahl | Globally optimal estimates for geometric reconstruction problems[END_REF] or prescribe some camera configurations [START_REF] Xuelian | New fundamental matrix estimation method using global optimization[END_REF][START_REF] Chesi | Estimating the fundamental matrix via constrained least-squares: A convex approach[END_REF][START_REF] Kahl | Globally optimal estimates for geometric reconstruction problems[END_REF] (requiring typically that the epipole in the first camera does not lie at infinity). Furthermore, the main criticism made to the optimization method we use is the resolution of a hierarchy of convex linear problems of increasing size, which is computationally ineffective and numerically unstable. The proposed solution overcomes this drawback: experiments show that, in most of cases, the proposed algorithm only requires solving the second relaxation of the sequences used in the exposed optimization method.

Our experimental evaluation on simulated and real datasets compares the difference between the eight-point algorithm and ours used as initialization to bundle adjustment. We observe that (i) bundle adjustment consistently converges within less iterations with our initialization and (ii) bundle adjustment always achieves an equal or lower reprojection error with our initialization. We provide numerous examples of real image pairs from standard datasets. They all illustrate practical cases for which our initialization method allows bundle adjustment to reach a better local minimum than the eight-point algorithm.

State of the Art

Accurately and automatically estimating the fundamental matrix from a pair of images is a major research topic. We first review a four-class categorization of existing methods, and specifically investigate the details of existing global methods. We finally state the improvements brought by our proposed global method.

Categorizing Methods

A classification of the different methods in three categories -linear, iterative and robust-was proposed in [START_REF] Armangué | Overall view regarding fundamental matrix estimation[END_REF]. Linear methods directly optimize a linear least squares cost. They include the eight-point algorithm [START_REF] Longuet-Higgins | A computer algorithm for reconstructing a scene from two projections[END_REF], SVD resolution [START_REF] Armangué | Overall view regarding fundamental matrix estimation[END_REF] and variants [START_REF] Tsai | Uniqueness and estimation of three-dimensional motion parameters of rigid objects with curved surfaces[END_REF][START_REF] Hartley | In Defence of the 8-point Algorithm[END_REF][START_REF] Torr | Invariant fitting of two view geometry or[END_REF][START_REF] Chojnacki | Revisiting hartley's normalized eight-point algorithm[END_REF]. Iterative methods iteratively optimize a nonlinear and nonconvex cost. They require, and are sensitive to the quality of, an initial estimate. The first group of iterative methods minimizes the distances between points and epipolar lines [START_REF] Salvi | An Approach to Coded Structured Light to Obtain Three Dimensional Information[END_REF][START_REF] Chen | Why not use the lm method for fundamental matrix estimation?[END_REF]. The second group minimizes some approximation of the reprojection error [START_REF] Luong | The fundamental matrix: Theory, algorithms, and stability analysis[END_REF][START_REF] Zhang | Determining the epipolar geometry and its uncertainty : A review[END_REF][START_REF] Van Den | A new constrained parameter estimator: Experiments in fundamental matrix computation[END_REF][START_REF] Chojnacki | A new approach to constrained parameter estimation applicable to some computer vision problems[END_REF]. The third group of methods minimizes the reprojection error, and are equivalent to two-view projective bundle adjustment. Iterative methods typically use a nonlinear parameterization of the fundamental matrix which guarantees that the rank-deficiency constraint is met. For instance, a minimal 7-parameter update can be used over a consistent orthogonal representation [START_REF] Bartoli | Nonlinear estimation of the fundamental matrix with minimal parameters[END_REF]. Finally, robust methods estimate the fundamental matrix while classifying each point match as inlier or outlier. Robust methods use M-Estimators [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF], LMedS (median least squares) [START_REF] Zhang | Determining the epipolar geometry and its uncertainty : A review[END_REF] or RANSAC (random sampling consensus) [START_REF] Torr | The development and comparison of robust methods for estimating the fundamental matrix[END_REF]. Both LMedS and RANSAC are stochastic.

To these three categories, we propose to add a fourth one: global methods. Global methods attempt at finding the global minimum of a nonconvex problem. Convex relaxations have been used to combine a convex cost with the rankdeficiency constraint [START_REF] Chesi | Estimating the fundamental matrix via constrained least-squares: A convex approach[END_REF]. However, these relaxations do not converge to a global minimum and the solution's optimality is not certified.

Global Methods

In theory, for a constrained optimization problem, global optimization methods do not require an initial guess and may be guaranteed to reach the global minimum, thereby certifying optimality. Such global methods can be separated in two classes. The methods of the first class describe the search space as exhaustively as possible in order to test as many candidate solutions as possible. Following this way, there are methods such as Monte-Carlo sampling, which test random elements satisfisying constraints, and reactive tabu search [START_REF] Glover | Future paths for integer programming and links to artificial intelligence[END_REF][START_REF] Hansen | Algorithms for the maximum satisfiability problem[END_REF], which continues searching even after a local minimum has been found. The major drawback of these methods is mainly in the prohibitive computation time required to have a sufficiently high probability of success. Moreover, even in case of convergence, there is no certificate of global optimality. Contrary to the methods of the first class, methods lying in the second class provide a certificate of global optimality using the mathematical theory from which they are built. Branch and Bound algorithms [START_REF] Land | An automatic method of solving discrete programming problems[END_REF] or global optimization by interval analysis [START_REF] Hansen | Global optimization using interval analysis[END_REF][START_REF] Moore | Introduction to Interval Analysis[END_REF] are some examples. However, although these methods can be faster than those of the first category, their major drawback is their lack of generality. Indeed, these methods are usually dedicated to one particular type of cost function because they use highly specific computing mechanisms to be as efficient as possible. A review of global methods may be found in [START_REF] Weise | Global Optimization Algorithms -Theory and Application[END_REF].

A good deal of research has been conducted over the last few decades on applying global optimization methods in order to solve polynomial minimization problems under polynomial constraints. The major drawback of these applications has been the difficulty to take constraints into account. But, by solving simplified problems, these approaches have mainly been used to find a starting point for local iterative methods. However, recent results in the areas of convex and polynomial optimization have facilitated the emergence of new approaches. These have attracted great interest in the computer vision community. In particular, global polynomial optimization [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF][START_REF] Henrion | Gloptipoly: Global optimization over polynomials with Matlab and SeDuMi[END_REF] has been used in combination with a finite-epipole nonlinear parameterization of the fundamental matrix [START_REF] Xuelian | New fundamental matrix estimation method using global optimization[END_REF]. This method does not consequently cover camera setups where the epipole lies at infinity. A global convex relaxation scheme [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF][START_REF] Henrion | Gloptipoly: Global optimization over polynomials with Matlab and SeDuMi[END_REF] was used to minimize the Sampson distance [START_REF] Kahl | Globally optimal estimates for geometric reconstruction problems[END_REF]. Because this implies minimizing a sum of many rational functions, the generic optimization method had to be specifically adapted and lost the property of certified global optimality.

The Proposed Method

The proposed method lies in the fourth category: it is a global method. Similarly to the eight-point algorithm, it minimizes the algebraic cost, but explicitly enforces the nonlinear rank-deficiency constraint. Contrarily to previous global methods [START_REF] Xuelian | New fundamental matrix estimation method using global optimization[END_REF][START_REF] Chesi | Estimating the fundamental matrix via constrained least-squares: A convex approach[END_REF][START_REF] Kahl | Globally optimal estimates for geometric reconstruction problems[END_REF], the proposed method handles all possible camera configurations (it does not make an assumption on the epipoles being finite or infinite) and certifies global optimality. Moreover, the presented algorithm is based on the resolution of a very short sequence of convex linear problems and is therefore computationally efficient.

A large number of attemps to introduce global optimization have been made in the literature. In [START_REF] Chesi | Estimating the fundamental matrix via constrained least-squares: A convex approach[END_REF], a dedicated hierarchy of convex relaxations is defined in order to globally solve the problem of fundamental matrix estimation. The singularity constraint is taken into account by introducing additional optimization variables, not introduced directly in the problem description. The resulting optimization algorithm is not generic and, contrary to Lasserre's hierarchy, there is no proof that the sequence of solutions of this specific hierarchy converges to the global minimum.

In [START_REF] Xuelian | New fundamental matrix estimation method using global optimization[END_REF], Lasserre's hierarchy is used jointly with the introduction of the singularity constraint in the problem description. However, the normalization constraint F 2 = 1 is replaced by fixing, a priori, one of the F coefficients to 1. Hence the other F coefficients are not bounded. Consequently, there is no guarantee that the sequence of solutions ( f t ) t∈N converges to the global minimum.

In [START_REF] Kahl | Globally optimal estimates for geometric reconstruction problems[END_REF] the authors minimize the Sampson distance by solving a specific hierarchy of convex relaxations built upon an epigraph formulation. This specific hierarchy is obtained by extending Lasserre's hierarchy, where the contraints are linear (LMIs), to a hierarchy of convex relaxations where constraints are polynomial (polynomial matrix inequality -PMI). Consequently, no asymptotic convergence of the hierarchy to the global optimum can be guaranteed. Moreover, even if the rank correction is achieved using the singularity constraint, the normalization is replaced by setting, a priori, the coefficient f 33 to 1 and so the other coefficients of F are not bounded.

Finally, in a very recent work [START_REF] Zheng | A practical rank-constrained eightpoint algorithm for fundamental matrix estimation[END_REF], the algebraic error is globally minimized thanks to the resolution of seven subproblems. Each subproblem is reduced to a polynomial equation system solved via a Gröbner basis solver. The singularity constraint is satisfied thanks to the right epipole parametrization. Although this parametrization ensures that F is singular while using the minimum number of parameters, this method is not practical since it would be necessary to solve 126 subproblems in order to cover all the 18 possible parameter sets [START_REF] Zhang | Determining the epipolar geometry and its uncertainty : A review[END_REF]. Therefore it is preferable to introduce the singularity constraint directly in the problem description rather than via some parametrization of F.

Polynomial Global Optimization

Introduction

Given a real-valued polynomial f (x) : R n → R, we are interested in solving the problem:

f ⋆ = inf x∈K f (x) (5)
where K ⊆ R n is a (not necessarily convex) compact set defined by polynomial inequalities: g j (x) ≥ 0, j = 1, . . . , m. Our optimization method is based on an idea first described in [START_REF] Lasserre | Optimisation globale et théorie des moments[END_REF]. It consists in reformulating the nonconvex global optimization problem (5) as the equivalent convex linear programming problem:

f = inf µ∈P(K) K f (x)dµ, (6)
where P(K) is the set of probability measures supported on K. Note that this reformulation is true for any continuous function (not necessarily polynomial) and any compact set

K ⊆ R n . Indeed, as f ⋆ ≤ f (x), then f ⋆ ≤ K f dµ and thus f ⋆ ≤ f . Conversely, if
x ⋆ is a global minimizer of ( 5), then the probability measure µ ⋆ △ = δ x ⋆ (the Dirac at x ⋆ ) is admissible for [START_REF] Chojnacki | Revisiting hartley's normalized eight-point algorithm[END_REF]. Moreover, because f is a solution of ( 6), the following inequality holds: K f (x)dµ ≥ f , ∀µ ∈ P(K) and thus f ⋆ = K f (x) δ x ⋆ ≥ f . Instead of optimizing over the finite-dimensional eu- clidean space K, we optimize over the infinite-dimensional set of probability measures P(K). Thus, Problem ( 6) is, in general, not easier to solve than Problem [START_REF] Torr | Invariant fitting of two view geometry or[END_REF]. However, in the special case of f being a polynomial and K being defined by polynomial inequalities, we will show how Problem ( 6) can be reduced to solving an (generically finite) sequence of convex linear matrix inequality (LMI) problems.

Notations and Definitions

First, given vectors α = (α 1 , . . . , α n ) ⊤ ∈ N n and x = (x 1 , . . . , x n ) ⊤ ∈ R n , we define the monomial x α by:

x α △ = x α 1 1 x α 2 2 . . . x α n n ( 7
)
and its degree by deg(x α )

△ = α 1 = n ∑ i=1 α i .
For t ∈ N, we define N n t the space of the n-dimensional integer vector with a norm lower than t as:

N n t △ = {α ∈ N n | α 1 ≤ t} . (8)
Then, consider the family:

{x α } α∈N n t = 1, x 1 , x 2 , . . . , x n , x 2 1 , x 1 x 2 , . . . , x 1 x n , x 2 x 3 , . . . , x 2 n , . . . , x t 1 , . . . , x t n (9)
of all the monomials x α of degree at most t, which has dimension s(t)

△ = (n + t)! t!n! .
Those monomials form the canonical basis of the vector space R t [x] of real-valued multivariate polynomials of degree at most t. Then, a polynomial p ∈ R t [x] is understood as a linear combination of monomials of degree at most t:

p(x) = ∑ α∈N n t p α x α , (10) and p
△ = (p α ) α 1 ≤t ∈ R N n t ≃ R s(t)
is the vector of its coefficients in the monomial basis {x α } α∈N n t . Its degree is equal to deg(p) Example: The polynomial

△ = max { α 1 | p α = 0}
x ∈ R 2 → p(x) = 1 + 2x 2 + 3x 2 1 + 4x 1 x 2 (11)
has a vector of coefficients p ∈ R 6 with entries p 00 = 1, p 10 = 0, p 01 = 2, p 20 = 3, p 11 = 4 and p 02 = 0.

Next, given y = (y α ) α∈N n ∈ R N n , we define the Riesz functional L y by the linear form:

L y : R [x] → R p = ∑ α∈N n p α x α → y ⊤ p = ∑ α∈N n p α y α . (12)
Thus, the Riesz functional can be seen as an operator that linearizes polynomials.

Example: For the polynomial [START_REF] Chesi | Estimating the fundamental matrix via constrained least-squares: A convex approach[END_REF], the Riesz functional reads

p(x) = 1 + 2x 2 + 3x 2 1 + 4x 1 x 2 → L y (p) = y 00 + 2y 01 + 3y 20 + 4y 11 . ( 13
)
For t ∈ N and y ∈ R N n 2t , the matrix M t (y) of size s(t) defined by:

(M t (y)) α,β = L y (x α x β ) = y α+β ∀α, β ∈ N n t ( 14
)
is called the moment matrix of order t of y. By construction, this matrix is symmetric and linear in y. Then, given q ∈ R t [x] and q ∈ R N n t the vector of its coefficients in the monomial basis, the vector:

qy △ = M t (y)q ∈ R N n t (15)
is called the shifted vector with respect to q. M t (qy), the moment matrix of order t of qy, is called the localizing matrix of degree t of q. This matrix is also symmetric and linear in y.

Example: If n = 2 then:

M 0 (y) = y 00 , M 1 (y) =  
        , (16)
and if q(x) = a + 2x Finally, recall that a symmetric matrix F ∈ S n is positive semidefinite, denoted by F 0, if and only if x ⊤ Fx ≥ 0, ∀x ∈ R n or equivalently, if and only if the minimum eigenvalue of F is non-negative. A linear matrix inequality (LMI) is a convex constraint:

F 0 + n ∑ k=1 x k F k 0, (18)
on a vector x ∈ R n , where matrices F k ∈ S m , k = 0, . . . , n are given.

Optimization Method

Let f be a real-valued multivariate polynomial, Problem (6) can be reduced to a convex linear programming problem. Indeed, if f (x) = ∑ α∈N n f α x α then:

K f dµ = K ∑ α∈N n f α x α dµ = ∑ α∈N n f α K x α dµ = L y ( f ) (19)
where each coordinate y α of the infinite sequence y ∈ R N n is equal to

K x α µ(dx),
also called the moment of order α. Consequently, if f is polynomial, then Problem ( 6) is equivalent to:

f = inf L y ( f ) s.t. y 0 = 1 y ∈ M K . (20)
with:

M K △ = y ∈ R N n | ∃µ ∈ M + (K) such that y α = K x α dµ ∀α ∈ N n , (21)
and M + (K) is the space of finite Borel measures supported on K. Remark that the constraint y 0 = 1 is added in order to impose that if y ∈ M K then y represents a measure in P(K) (and no longer in M + (K)). Although Problem ( 20) is a convex linear programming problem, it is difficult to describe the convex cone M K with simple constraints on y. But, the problem "y ∈ M K ", also called K-moment problem, is solved when K is a basic semi-algebraic set, namely:

K △ = {x ∈ R n | g 1 (x) 0, . . . , g m (x) 0} (22)
where

g j ∈ R[x], ∀ j = 1, . . . m.
Note that K is assumed to be compact. Then, without loss of generality, we assume that one of the polynomial inequalities

g j (x) 0 is of the form R 2 -x 2 2
0 where R is a sufficiently large positive constant. In this particular case, M K can be modelled using LMI conditions on the matrices M t (y) and M t (g j y), j = 1, . . . m. More precisely, thanks to the Putinar Theorem [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF][START_REF] Laurent | Sums of squares, moment matrices and optimization over polynomials[END_REF], we have:

M K = M (g 1 , . . . , g m ), (23)
where:

M (g 1 , . . . , g m ) △ = y ∈ R N n | M t (y) 0, M t (g j y) 0 ∀ j = 1, . . , m ∀t ∈ N . ( 24
) Then, Problem ( 6) is equivalent to:

f = inf y ∈ R N n L y ( f ) s.t. y 0 = 1 M t (y) 0 M t (g j y) 0 j = 1, . . , m ∀t ∈ N. (25)
To summarize, if f is polynomial and K a semi-algebraic set, then Problem ( 5) is equivalent to a convex linear programming problem with an infinite number of linear constraints on an infinite number of decision variables. Now, for t

≥ d K △ = max(d f , d g 1 , . . . , d g m ) consider the finite-dimensional truncations of Problem (25): Q t △ =            f t △ = min y ∈ R N n 2t L y ( f ) s.t. y 0 = 1 M t (y) 0, M t-d g j (g j y) 0 ∀ j ∈ {1, . . . , m} . (26)
By construction, Q t , t ∈ N generates a hierarchy of LMI relaxations of Problem (25) [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF], where each Q t , t ∈ N, is concerned with moment and localizing matrices of fixed size t. Each relaxation ( 26) can be solved by using publicdomain implementations of primal-dual interior point algorithms for semidefinite programming (SDP) [START_REF] Borchers | Csdp, a c library for semidefinite programming[END_REF][START_REF] Benson | DSDP5 user guide -software for semidefinite programming[END_REF][START_REF] Fujisawa | Sdpa (semidefinite programming algorithm) -user's manual[END_REF][START_REF] Toh | Solving semidefinite-quadraticlinear programs using sdpt3[END_REF][START_REF] Sturm | Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones[END_REF]. When the relaxation order t ∈ N tends to infinity, we obtain the following results [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF][START_REF] Nie | Optimality conditions and finite convergence of Lasserre's hierarchy[END_REF]:

f t ≤ f t+1 ≤ f and lim t→+∞ f t = f . ( 27
)
Practice reveals that this convergence is fast and very often finite, i.e. there exists a finite t 0 such that f t = f , ∀t ≥ t 0 . In fact, finite convergence is guaranteed in a number of cases (e.g. discrete optimization) and very recent results by Nie [START_REF] Nie | Optimality conditions and finite convergence of Lasserre's hierarchy[END_REF] show that finite convergence is even generic as well as an optimal solution y ⋆ t of [START_REF] Weise | Global Optimization Algorithms -Theory and Application[END_REF].

Example: Consider the polynomial optimization problem

f = min x∈R 2 -x 2 s.t. 3 -2x 2 -x 2 1 -x 2 2 ≥ 0 -x 1 -x 2 -x 1 x 2 ≥ 0 1 + x 1 x 2 ≥ 0. ( 28
)
The first LMI relaxation Q 1 is

f 1 = min y∈R 6 -
  0. ( 30
) It can be checked that f 1 = -2 ≤ f 2 = f = -1+ √ 5
2 . Note that the constraint 3 -2x 2x 2 1x 2 2 ≥ 0 certifies boundedness of the feasibility set.

However, we do not know a priori at which relaxation order t 0 the convergence occurs. Practically, to detect whether the optimal value is attained, we can use conditions on the rank of the moment and localization matrices. Indeed, let y ⋆ ∈ R N n 2t be a solution of Problem (26) at a given relaxation order t ≥ d K , if:

rank(M t (y ⋆ )) = rank(M t-d K (y ⋆ )) (31)
then f t = f . In particular, if rank(M t (y ⋆ )) = 1 then condition (31) is satisfied. Moreover, if these rank conditions are satisfied, then we can use numerical linear algebra to extract rank(M t (y ⋆ )) global optima for Problem [START_REF] Torr | Invariant fitting of two view geometry or[END_REF]. We do not describe the algorithm in this article, but the reader can refer to [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]Sections 4.3] for more advanced information. Figure 3.3 summarizes the optimization process. A Matlab interface called GloptiPoly [START_REF] Henrion | Gloptipoly 3: moments, optimization and semidefinite programming[END_REF] has been designed to construct Lasserre's LMI relaxations in a format understandable by any SDP solver interfaced via YALMIP [START_REF] Lfberg | Yalmip : A toolbox for modeling and optimization in MATLAB[END_REF]. It can be used to construct an LMI relaxation (26) of a given order corresponding to a polynomial optimization problem [START_REF] Torr | Invariant fitting of two view geometry or[END_REF] with given polynomial data entered symbolically. A numerical algorithm is implemented in GloptiPoly to detect global optimality of an LMI relaxation, using the rank

min x∈K f (x) min µ∈P(K) K f (x)dµ inf L y ( f ) △ = ∑ α f α y α s.t. y ∈ y ∈ R N n | ∃µ ∈ P(K) : y α = K x α dµ ∀α ∈ N n inf L y ( f ) s.t. y 0 = 1, M t (y) 0, M t (g j y) 0 j = 1, . . , m ∀t ∈ N. Q t            f t △ = min y ∈ R N n 2t L y ( f ) s.t. y 0 = 1 M t (y) 0, M t-d g j (g j y) 0 j = 1, . . , m.
f continuous and K compact tests [START_REF] Laurent | Sums of squares, moment matrices and optimization over polynomials[END_REF]. The algorithm also extracts numerically the global optima from the moment matrix. This approach has been successfully applied to globally solve various polynomial optimization problems (see [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF] for an overview of results and applications). In computer vision this approach was first introduced in [START_REF] Kahl | Globally optimal estimates for geometric reconstruction problems[END_REF] and used in [START_REF] Kahl | Globally optimal estimates for geometric reconstruction problems[END_REF].

f polynomial K = {x ∈ R n | g 1 (x) 0, . . . , g m (x) 0} f = ∑ α∈N n t f α K x α dµ

Application to Fundamental Matrix Estimation

An important feature of our approach is that the singularity constraint can be directly satisfied by our solution and that we do not need an initial estimate as in other methods. Nevertheless, the problem being homogeneous, an additional normalization constraint is needed to avoid the trivial solution F = 0. This is generally done by setting one of the coefficients of the F matrix to 1. However with such a normalization there is no guarantee that the feasible set is compact which is a necessary condition for the convergence of our polynomial optimization method. Moreover, this normalisation excludes a priori some geometric configurations. Thus, to guarantee compactness of the feasible set and to avoid the trivial solution, we include the additional normalization constraint F 2 = 1.

Alternatively, the singularity constraint can be inferred by parameterizing the F matrix using one or two epipoles. For instance, in [START_REF] Xuelian | New fundamental matrix estimation method using global optimization[END_REF] the authors have tried to use polynomial optimization to estimate the F matrix, parametrized by using a single epipole. However such an approach has several drawbacks. Firstly, the coefficients of F are not bounded and the convergence of the method is not guaranteed. Secondly, using the epipole explicitly increases the degree of the polynomial criterion and consequently, the size of the corresponding relaxations in the hierarchy. This results in a significant increase of the computational time. Finally, the choice of this parametrization is arbitrary and does not cover all camera configurations.

Our method is summarized in Algorithm 1 below. Its main features are:

• In contrast with [START_REF] Kahl | Globally optimal estimates for geometric reconstruction problems[END_REF][START_REF] Kahl | Globally optimal estimates for geometric reconstruction problems[END_REF] the optimization problem is formulated with an explicit Frobenius norm constraint on the decision variables. This enforces compactness of the feasibility set which is included in the Euclidean ball of radius 1. We have observed that enforcing this Frobenius norm constraint has a dramatic influence on the overall numerical behavior of the SDP solver, especially with respect to convergence and extraction of global minimizers.

• We have chosen the SDPT3 solver [START_REF] Toh | Sdpt3 -a matlab software package for semidefinite programming[END_REF][START_REF] Toh | Solving semidefinite-quadraticlinear programs using sdpt3[END_REF] since our experiments revealed that for our problem it was the most efficient and reliable solver.

• We force the interior-point algorithm to increase the accuracy as much as possible, overruling the default parameter set in SDPT3. Then the solver runs as long as it can make progress.

• Generically, polynomial optimization problems on a compact set have a unique global minimizer; this is confirmed in our numerical experiments as the moment matrix has almost always rank-one (which certifies global optimality) at the second SDP relaxation of the hierarchy. In some (very few) cases, due to the numerical extraction, the global minimum is not fully accurate but yet largely satisfactory.

Algorithm 1 Polynomial optimization for fundamental matrix estimation Require: Matched points (q i , q ′ i ), i = 1, . . . , n 1: Create the cost function Crit = ∑ n i=1 q ′⊤ i Fq i 2 :

mpol('F',3,3); for k = 1:size(q1) n(k) = (q2'*F*q1)^2; end; Crit = sum(n); 2: Create the constraints det(F) = 0 and F 2 = 1: K_det = det(F) == 0; K_fro = trace(F*F') == 1; 3: Fix the accuracy of the solver to 0, then the solver runs as long as it can make progress: pars.eps = 0; mset(pars); 4: Change the default SDP solver to SDPT3: mset('yalmip',true); mset(sdpsettings('solver','sdpt3')); 5: Form the second LMI relaxation of the problem: P = msdp(min(crit),K_det,K_fro,2); 6: Solve the second LMI relaxation: msol(P);

Experimental Results

This section presents results obtained by the test procedure presented below with the 8-point method and our global method. First, criteria to evaluate the performance of a fundamental matrix estimate are described. Next, the evaluation methodology is detailed. Experiments were then carried out on synthetic data to test the sensitivity to noise and to the number of point matches. Finally, experiments on real data were performed to confirm previous results and to study the influence of the type ofmotion between the two images.

Evaluation Criteria

Various evaluation criteria were proposed in the literature [START_REF] Zhang | Determining the epipolar geometry and its uncertainty : A review[END_REF] to evaluate the quality of a fundamental matrix estimate. Driven by practice, a fundamental matrix estimate F is evaluated with respect to the behavior of its subsequent refinement by projective bundle adjustment. Bundle Adjustment from two uncalibrated views is described as a minimization problem. The cost function is the RMS reprojection errors. The unknowns are the 3D points Q i , i = 1, . . . , n and the projection matrices P and P ′ . The criteria we use are:

1. The initial reprojection error written e Init (F).

2. The final reprojection error e BA (F).

The number of iterations taken by Bundle Adjustment to converge, Iter(F).

These three criteria assess whether the estimates provided by the two methods, denoted by F 8pt and F Gp , are in a 'good' basin of attraction. Indeed, the number of iterations gives an indication of the distance between the estimate and the optimum while e BA (F) gives an indication on the quality of the optimum.

Evaluation Method

The following algorithm summarizes our evaluation method:

Evaluate(F)
1. Inputs: fundamental matrix estimate F, n point matches (q i , q ′ i ), i = 1, . . . , n

Form initial projective cameras [41]:

Find the second epipole from F ⊤ e ′ ∼ 0 (3×1)

Find the canonical plane homography

H * ∼ [e] × F Set P ∼ I (3×3) 0 (3×1) and P ′ ∼ [H * e ′ ]
3. Form initial 3D points [42]:

Triangulate each point independently by minimizing the reprojection error For each simulation series and for each parameter of interest (noise, number of points and number of motions), the same methodology is applied with the following four steps:

1. For two given motions between two successives images ([R k t k ]) and for a given matrix K of internal parameters, a set of 3D points (Q i ) i , i = 1, . . . , n is generated and two projection matrices P and P ′ are defined. In practice, the rotations matrices, R 1 and R 2 , of two motions are defined by:

R k △ =   cos(θ k ) 0 sin(θ k ) 0 1 0 -sin(θ k ) 0 cos(θ k )   with        θ 1 = π 3 and θ 2 = π 6 ( 32
)
and their translation vectors by t 1 = (20, 0, 5) ⊤ and t 2 = (6, 0, 0) ⊤ . These matrices are chosen such that [R 1 ,t 1 ] is a large movement and [R 2 ,t 2 ] is a small movement (see Figure 2). We simulated points lying in a cube with 10 meter side length. The first camera looks at the center of the cube and it is located 15 meters from the center of the cube. The focal length of the camera is 700 pixels and the resolution is 640 × 480 pixels.

Thanks to projection matrices

P = K [R 1 , t 1 ] and P ′ = K [R 2 , t 2 ]
, the set of 3D points (Q i ) i is projected into the two images as (q i , q ′ i ) i . At each of their pixel coordinates, a centered Gaussian noise with a variance σ 2 is added. In order to have statistical evidence, the results are averaged over 100 trials.

3. The resulting noisy points ( q i , q ′ i ) i are used to estimate F by our method F Gp and the reference 8-point method F 8pt . 4. Finally, via our evaluation procedure we evaluate the estimation error with respect to the noise standard deviation σ and the number of points n.

Sensitivity to Noise

We tested in two simulation series the influence of σ ranging from 0 to 2 pixels. The number of simulated points is 50. The first (resp. second) simulation series is based on the first motion [R 1 t 1 ] (resp. the second motion [R 2 t 2 ]). For the two motions, re-projection errors, e Init (F) or e BA (F), increase with the same slope when the noise level increases. Notice that for both movements, the Bundle-Adjustment step does not improve the results. Indeed, the noise gaussian noise is added to the projections (q i , q ′ i ) i . So this is noise which in practice would be produced by the extraction points process. Thus the solution produced by the resolution of the linear system is very close to the optimum and does not need to be refined. The initial solution provided by the triangulation step is then very close to a local minimum of the Bundle Adjustment problem. Moreover, the variation of the errors of initial re-projection before (8pt -Init and Gp -Init) and after (8pt -BA and Gp -BA) Bundle Adjustment versus the noise standard deviation is linear. However, the number of iterations needed for convergence is different in the two methods. The initial estimate of the triangulation computed from F Gp is closer to the local minimum than that obtained from F 8pt . For the first motion (large displacement between camera 1 and 2), the number of iterations of the global method (in green) remains smaller than for the 8-point method (in blue) even though their difference seems to decrease when the noise level is high (σ > 1). For a significant displacement the quality of the estimate F by the global method remains better even though the difference in quality diminishes with the noise level. Conversely, for the second motion (small displacement between the camera 1 and 2) both methods are equivalent since the difference in quality is only significant for a high level of noise (σ > 1). This is logical as the movement is less important. As a conclusion, the 8-point method provides a solution equivalent to that obtained with the global method when the displacement is not too important. For more significant movements the provided solution is not so close even though still in the same basin of attraction of a local minimum.

Influence of the Number of Points

In this experiment, we kept the noise level constant with a standard deviation σ 2 = 0.5 pixels. We tested the influence of the number of matches (q i , q ′ i ) i on the quality of the resulting estimate of F. The number of points N varied from 10 to 100. Two simulation series are also carried out with the two motions. Figure 4 brings with the same organization the evaluation criteria. It displays the influence of the number of matches for estimating F on the re-projection errors and on the number of iterations. For both motions and for a sufficiently high number of matches (N > 50), re-projection errors, before and after refinement with Bundle Adjustment, or the number of iterations versus the number of matches converge to the same asymptote. From a high number of matches, the initial estimate from triangulation computed with F 8pt and with F Gp are both in the same basin of attraction for the Bundle Adjustment problem. However, for a number of matches smaller than 50, the number of iterations to converge is smaller for given re-projection errors. The quality of the estimation by the global method seems better. The initial estimate from triangulation computed with F 8pt goes away from the basin of convergence whereas the one computed with F Gp remains in the basin.

Influence of the Number of Points with Wide Baseline

In order to sustain the previous behavior, for a noise standard deviation of σ 2 = 1 pixel and for the significant displacement [R 1 t 1 ], the influence of the number of matching points on the re-projection errors and on the number of iterations was tested. In this difficult context, Figure 5 demonstrates that the initial estimate computed from F Gp is always closer to the local minimum than that computed from F 8pt . No matter what is the number of matching points, the number of iterations needed to converge is always smaller.

As a conclusion, the quality of solutions obtained by both methods is almost identical when the movement is not too important, the number of matching points is sufficiently large, and the noise level is not too high. However, when one of these three parameters varies then the 8-point method lacks precision whereas the global method still allows Bundle Adjustment to convergence to the global minimum. The 8-point method computes the projection of an unconstrained local minimizer on the feasible set whereas the global method provides a global minimizer of the constrained optimization problem. It is already surprising that even for good values of the three parameters the resulting solutions are not too far apart. But for worst values of the parameters it would be even more surprising.

Experiments on Real Data

The evaluation criteria remain the same, e Init (F), e BA (F) and Iter(F) and the computation time is added. Two experiments were carried out with two sets of images that illustrate different motions between two successive images.

Experiment 1

The first set of four images (see Table 6) shows all possible epipolar configurations (right or left epipole at infinity . . . every pair of images, 60 matches are available to compute an estimate of F. The values of the evaluation criteria are summarized in table 6. No matter what pair of images is used, the re-projection errors and the number of iterations are almost always better when F Gp is used as initial guess. In addition, for three motions (A -C, A -D and C -D), in contrast with the initial guess F Gp , the initial guess from the 8-point method is not in a better basin of attraction. This may explain why the initial re-projection errors e Init (F) are sometimes larger for F Gp as the initial guess may be in a good basin of attraction but with a larger re-projection error. For the four motions A -B, B -C, B -D and B -D, both initializations are in the same basin of attraction but the number of iterations demonstrates that the initial guess from the global method is always closer to the local minimizer. Finally, even though the computation time of the latter is significantly larger than for the 8-point method, it still remains compatible with a practical use.

Experiment 2

The second experiment compares the two methods on large motions. It is based on many series of images. First we test our algorithm with the classic series Library, Merton, dinosaur and house that are available at www.robots.ox.ac. uk/ ~vgg/data/data-mview.html. For the set of three images of Library and Merton serie, Table 7, 8, 9 and Table 10 demonstrate that the quality of the solution achieved by the global method is always better than with the 8-point method (in some cases both solutions are very close).

We also conducted the same tests on other pairs of images. For the first pair, we used images from a standard cylinder graciously provided by the company NOOMEO (http://www.noomeo.eu/). This cylinder is use to evaluate the accuracy of 3D reconstructions. Matched points are calculated with digital image correlation method. They are located in a window inside the cylinder. Thus, we have 6609 pairs (q i , q i ) i matched to sub-pixel precision. Results are presented in the Table 11. We observe that the computation time of the 8-point method exceeds one second. This is due to the large number of matched points which leads to the resolution of a large linear system. However, as the points are precisely matched, this system is well conditioned. But the quality of the fundamental matrix estimated with the 8-point method is not sufficient to properly initialize the Bundle-Adjustment because the final re-projection error is 1.47 pixels. At the same time, even if the number of iterations is larger, our global method supplies a good estimation because the final re-projection error is 0.25 pixels. Furthermore, the calculation time remains constant in approximately 2 seconds. For the second pair, we use images taken by an endoscope. Table 12 shows the results obtained on this difficult case. As for the previous example, we observe that the fundamental matrix estimated by our global method is good quality because the final error is 0.93 pixels. At the same time, Bundle Adjustement puts more iterations to converge on a less precise solution when we use F 8pt to initialize it.

For the set of 36 images of the Dinosaur series and 9 images of the house series, we tested the influence of motion amplitude between a pair of image on the quality of the resulting estimates obtained by both methods. For this purpose, we had both estimates with all possible motions ((0, 1), (1, 2), (2, 3), . . .) with 1-image distance, then all possible motions ((0, 2), (1, 3), (2, 4), . . .) with 2-image distance, and so on. With this process, we can measure the influence of the average angle on the quality of the fundamental matrix estimated by both methods. Figures from 13 to 18 shows the average of re-projection errors and the average of number of iterations with respect to average angle for the two series. The re-projection error after Bundle Adjustment is always smaller with the global method and with always a smaller number of iterations. Next, the larger the movement the more the solution by both methods deteriorates. But the deterioration is larger for the 8-point method than for the global method. One may also observe that in some cases the re-projection error before Bundle Adjustment is in favor of the 8-point method. In analogy with the real cases studied before, this may be due to the fact that for these cases the initial guess F Gp is in a basin of attraction with a better local minimum than in the basin of attraction associated with F 8Pt , but the 'distance' between the initial guess and the corresponding local minimizer is larger for F Gp than for F 8Pt . Indeed in such cases the number of iterations is larger for F Gp than for F 8Pt .

Conclusion

We have studied the problem of estimating globally the fundamental matrix over nine parameters and under rank and normalisation constraints. We have proposed a polynomial-based approach which enables one to estimate the fundamental matrix with good precision. More generally, we have shown how to modify the constraints on the numerical certificate of optimality to obtain fast and robust convergence. The method converges in a reasonable amount of time compared to other global optimization methods.

From computational experiments conducted on both simulated and real data we conclude that the global method always provides an accurate initial estimation for the subsequent bundle adjustment step. Moreover, we have shown that if the eight-point method has a lower computational cost, its resulting estimate is frequently far from the global optimum obtained by the global method.

- Views nb Points e Init (F) e BA (F) Iter(F) Time (s) Views nb Points e Init (F) e BA(F) Iter(F) Time (s)

A - -B - -C - -D - Views Epipoles e Init (F) e BA (F) Iter(F) Time (s) e e ′ F 8Pt F Gp F 8Pt F Gp F 8Pt F Gp F 8Pt F Gp A B ∞ ∞ 0.
F 8Pt F Gp F 8Pt F Gp F 8Pt F Gp F 8Pt F
F 8Pt F Gp F 8Pt F Gp F 8Pt F Gp F 8Pt F

Figure 1 :

 1 Figure 1: Polynomial optimization process; see the main text for details.

4 .

 4 Compute e Init (F) 5. Run two-view uncalibrated bundle adjustment [19] 6. Compute e BA (F) and Iter(F) 7. Outputs: e Init (F), e BA (F) and Iter(F) 4.3 Experiments on Simulated Data 4.3.1 Simulation Procedure

  Figure 3 gathers the influence of noise on the evaluation criteria. The first line shows the reproduction errors before, e Init (F), and after e BA (F) refinement through Bundle Adjustment with respect to the noise standard deviation. The second line shows the number of iterations Iter(F) of the Bundle Adjustment versus the noise standard deviation. The first (resp. second) row concerns the first (resp. second) motion.

Figure 2 :

 2 Figure 2: Projection of the cube in the camera on initial position (I) and in the camera after applying the rigid transformation [R 1 t 1 ] (II) and the rigid transformation [R 2 t 2 ] (III).

Figure 3 :

 3 Figure 3: For two movements, [R 1 t 1 ] (left column) and [R 2 t 2 ] (right column), reprojection errors and number or iterations measured against image noise.

Figure 4 :Figure 5 :

 45 Figure 4: For two movements, [R 1 t 1 ] (left column) and [R 2 t 2 ] (right column), reprojection errors and number or iterations measured against number of points for a gaussian noise with a variance fixed to 0.5.

Figure 7 :

 7 Figure 7: Reprojection Error before (e Init (F)) and after Bundle Adjustment (e BA (F)), Number of Iterations (Iter(F)), and CPU time to compute F (Time), obtained when combining pairs of images of the Library series

Figure 9 :

 9 Figure 9: Reprojection Error before (e Init (F)) and after Bundle Adjustment (e BA (F)), Number of Iterations (Iter(F)), and CPU time to compute F (Time), obtained when combining pairs of images of the Merton2 series

Figure 11 :Figure 12 :Figure 13 :Figure 14 :Figure 15 :Figure 16 :Figure 17 :Figure 18 :

 1112131415161718 Figure 11: Reprojection Error before (e Init (F)) and after Bundle Adjustment (e BA (F)), Number of Iterations (Iter(F)), and CPU time to compute F (Time), obtained when combining pairs of images of the Cylinder series. The matched points are located in blue bounding boxes.

  01 y 20 y 11 y 02 y 10 y 20 y 11 y 30 y 21 y 12 y 01 y 11 y 02 y 21 y 12 y 03 y 20 y 30 y 21 y 40 y 31 y 22 y 11 y 21 y 12 y 31 y 22 y 13 y 02 y 12 y 03 y 22 y 13 y 04

			y 00 y 10	y 10 y 01 y 20 y 11	  ,
	M 2 (y) =	       	y 01 y 00	y 11 y 02 y 10 y

  2y 20 + 3y 02 ay 10 + 2y 30 + 3y 12 ay 01 + 2y 21 + 3y 03 ay 10 + 2y 30 + 3y 12 ay 20 + 2y 40 + 3y 22 ay 11 + 2y 31 + 3y 13 ay 01 + 2y 21 + 3y 03 ay 11 + 2y 31 + 3y 13 ay 02 + 2y 22 + 3y 04

	M 1 (qy) =	 	2 1 + 3x 2 2 then: ay 00 +



 .

[START_REF] Van Den | A new constrained parameter estimator: Experiments in fundamental matrix computation[END_REF]

  10 y 01 y 20 y 11 y 02 y 10 y 20 y 11 y 30 y 21 y 12 y 01 y 11 y 02 y 21 y 12 y 03 y 20 y 30 y 21 y 40 y 31 y 22 y 11 y 21 y 12 y 31 y 22 y 13 y 02 y 12 y 03 y 22 y 13 y 04 2y 01y 20y 02 3y 10 -2y 11y 30y 12 3y 01 -2y 02y 21y 3y 10 -2y 11y 30y 12 3y 20 -2y 21y 40y 22 3y 11 -2y 12y 31y 3y 01 -2y 02y 21y 03 3y 11 -2y 12y 31y 13 3y 02 -2y 03y 22y 10y 01y 11 -y 20y 11y 21 -y 11y 02y 12 -y 20y 11y 21 -y 30y 21y 31 -y 21y 31y 21 -y 11y 02y 12 -y 21y 12y 22 -y 12y 03y 13 00 + y 11 y 10 + y 21 y 01 + y 12 y 10 + y 21 y 20 + y 31 y 11 + y 22 y 01 + y 12 y 11 + y 22 y 02 + y 13

	and the second LMI relaxation Q 2 is
	f 2 = min y∈R 15	-y 01
	s.t.	y 00 = 1         y 00 y        	0,
		 	3y 00 -  0
		  -y   0
			y
		
			y 01
			s.t. y 00 = 1   y 00 y 10 y 01 y 10 y 20 y 11	  0	(29)
			y 01 y 11 y 02
			3y 00 -2y 01 -y 20 -y 02 ≥ 0 -y 10 -y 01 -y 11 ≥ 0 y 00 + y 11 ≥ 0,

  Reprojection Error before (e Init (F)) and after Bundle Adjustment (e BA (F)), Number of Iterations (Iter(F)), and CPU time to compute F (Time), obtained when combining pairs of images of the Merton1 series Gp F 8Pt F Gp F 8Pt F Gp 1 2 345 24.7158 29.6530 14.0481 2.7599 69 42 0.037 2.91 2 3 197 143.431 154.880 73.8776 72.8977 14 19 0.026 2.52 1 3 270 59.8109 77.9734 30.3824 15.6367 38 32 0.031 2.37

	1 2 2 3 1 3 Views nb Points 485 439 384 Figure 8: 1	15.1771 61.84817 63.79435 41.81854 40.2224 29 15.3126 9.8859 9.7360 77 51.3128 2.3690 7.5245 0.56180 5 2 3 e Init (F) e BA (F) Iter(F) F 8Pt F Gp F 8Pt F	Gp 55 0.048 2.02 28 0.020 2.22 52 0.040 2.71 Time (s)

  Reprojection Error before (e Init (F)) and after Bundle Adjustment (e BA (F)), Number of Iterations (Iter(F)), and CPU time to compute F (Time), obtained when combining pairs of images of the Merton3 series Gp F 8Pt F Gp F 8Pt F Gp 1 2 6609 6.3142 6.4021 1.4701 0.2531 401 237 1.23 2.33

	1 2 2 3 1 3 Figure 10: Views nb Points 400 197 264	62.7290 135.950 116.7659 118.1754 24.7184 13.5496 68.6030 52.6557 23.1740 140.801 83.3365 76.8614 13 5 9 e Init (F) e BA (F) Iter(F) F 8Pt F Gp F 8Pt F	Gp 32 0.041 3.03 14 0.025 2.52 11 0.032 2.39 Time (s)