A clever elimination strategy for efficient minimal solvers

2017·Arxiv

Abstract

Abstract

We present a new insight into the systematic generation of minimal solvers in computer vision, which leads to smaller and faster solvers. Many minimal problem formulations are coupled sets of linear and polynomial equations where image measurements enter the linear equations only. We show that it is useful to solve such systems by first eliminating all the unknowns that do not appear in the linear equations and then extending solutions to the rest of unknowns. This can be generalized to fully non-linear systems by linearization via lifting. We demonstrate that this approach leads to more efficient solvers in three problems of partially calibrated relative camera pose computation with unknown focal length and/or radial distortion. Our approach also generates new interesting constraints on the fundamental matrices of partially calibrated cameras, which were not known before.

1. Introduction

Computing multi-view geometry is one of the most basic and important tasks in computer vision [1]. These include minimal problem solvers [2, 3] in, e.g., structure from motion [4], visual navigation [5], large scale 3D reconstruction [6] and image based localization [7]. Fast, and efficient, minimal solvers are instrumental in RANSAC [8] based robust estimation algorithms [9].

In this paper we present a new insight into the systematic generation of minimal solvers [3], which leads to smaller and faster solvers. We explain our approach in the context of elimination theory [10] and we offer an interpretation of the theory that is useful for practice in computer vision.

Our main technical contribution is a new strategy for solving minimal problems. For many computer vision applications, that strategy allows to do more computation in an off-line stage and less computation in an on-line stage.

We exploit that many minimal problems in computer vi-

Figure 1. An illustration of the two equations (17) and (18), which define the f+E+f problem, cut by six linear equations for six image point correspondences.

sion lead to coupled sets of linear and polynomial equations where image measurements enter the linear equations only. We show how to solve such systems efficiently by first eliminating all unknowns which do not appear in the linear equations, and then extending solutions to the other unknowns. Moreover, our approach can be generalized to fully non-linear systems by linearization via monomial lifting [11].

We demonstrate that this approach leads to more efficient on-line solvers in three problems of partially calibrated relative camera pose computation with unknown focal length and/or radial distortion. Interestingly, our approach also generates new constraints on the fundamental matrices of partially calibrated cameras, which were not known before.

1.1. Related work

Historically, minimal problems [12, 13, 14, 15] addressed problems in geometry of one and two perspective cameras. Later, a more systematic approach to solving minimal problems in computer vision appeared, e.g., in [2, 16, 17, 18, 19, 20]. It developed a number of ad-hoc, as well as, systematic tools for solving polynomial systems appearing in computer vision. These were later used and improved by many researchers, e.g., [21, 22, 3, 23, 24, 25, 26, 27, 28, 29, 30, 31]. Lately, the algebraic geometry foundations for computer vision came into focus in algebraic vision [32, 33, 34, 35, 36].

One of the key elements in computer vision applications has been to design procedures for solving special polynomial systems that move the computation from the on-line stage of solving equations to an earlier off-line stage [37].

Interestingly, elimination theory [10] has not been fully exploited in such computer vision applications, although it has been used in many works [2, 20, 38, 3, 39] implicitly.

1.2. The main idea

Our main idea is to use elimination theory to do more computation in the off-line stage and less in the on-line stage.

Natural formulations of vision models often involve more unknowns than those that appear in the linear constraints depending on image measurements. For instance, the constraint det F = 0 in fundamental matrix computation does not involve any image measurements. We argue that it is advantageous to pre-process such models by computing its projection into the space of relevant unknowns. This is done by elimination. Solving the linear equations on the resulting projected variety is then fast. Subsequently, the values for the other unknowns can be determined using the Extension Theorem [10] from computer algebra.

2. Solving polynomial systems by elimination

A classical (textbook) strategy for solving systems of polynomial equations is to use elimination theory [40, 10, 41]. The strategy consists of two main steps.

1. First, the equations are “simplified” by eliminating some unknowns to get a set of equations from which the remaining unknowns can be computed. This provides a set of partial solutions.

2. Next, the partial solutions are extended to full solutions by substituting the partial solutions back into the original equations and solving for the remaining unknowns. We next explain different elimination strategies.

2.1. Elimination strategies

2.1.1 Standard textbook elimination strategy

Standard (textbook) elimination is based on the Elimination Theorem [10], which we review in Theorem 5.1 of the Appendix. It states that, for an ideal , we can read off LEX Gr¨obner bases for all elimination ideals from a LEX Gr¨obner basis G for I. Here the sequence of elimination ideals ends with . This is generated by one equation in the single unknown , which is then easy to solve numerically.

For a polynomial system with a finite number of solutions in , it is always possible [10, p. 254–255] to extend partial solutions from to . For this, we choose a single polynomial g with the lowest degree among all the univariate polynomials in after substituting the partial solution into the polynomials in .

2.1.2 Standard computer vision elimination strategy

In the existing minimal solvers, several different strategies for eliminating unknowns from the input equations were applied. These strategies were usually dependent on the spe-cific problem and were derived manually. Here we describe one strategy that was used in the vast majority of existing minimal solvers [2, 20, 38, 3, 42, 39].

Consider a system of m polynomial equations

in n unknowns . We assume that the set generates a zero dimensional ideal C[X], i.e. the system F has a finite number of solutions.

In this strategy the set F is partitioned into two subsets:

This means that contains the linear polynomials from F and contains the polynomials of higher degrees.

The linear equations can be rewritten as , where is a vector of all unknowns that appear in these equations. Then, the null space basis N of M, i.e. M N = 0, is used to parametrize the unknowns with new unknowns Y via . The parameterization is then plugged in the non-linear equations . The system is solved using, e.g., a Gr¨obner basis method and the automatic generator of efficient solvers [3]. The solutions Y are used to recover solutions .

2.1.3 A clever computer vision elimination strategy

In computer vision, we often encounter polynomial systems in which only the linear equations depend on image measurements, while the nonlinear equations stay the same, regardless of image measurements. For example, the epipolar constraint [1] generates linear equations that depend on the input measurements while the singularity of the fundamental matrix [1] results in the non-linear equation det(F) = 0, which does not depend on the measurements.

Here we present a new “clever” elimination strategy, which usually allows us to do more computation in the off-line stage and less computation in the on-line stage.

Throughout this section we assume that the nonlinear equations do not depend on image measurements, i.e. for all instances these equations are the same. Later, in Sec- tion 3.3, we will show how to deal with problems when contains equations that depend on image measurements.

Let us now describe our new elimination strategy. We first divide the input equations F into the linear equations in (1) and the non-linear equations in (2). Moreover, we divide the n given unknowns X into two subsets:

The set contains the unknowns that appear in linear equations. The set contains the unknowns that appear in equations of higher degree only. We fix the following notation: , which means that and . Now, the idea of our new elimination method is to eliminate all unknowns from non-linear equations . The non-linear equations do not depend on image measurements, and are the same for all instance of a given problem. Therefore, we can perform this elimination off-line, in the pre-processing step. This elimination may be computationally demanding. However, since we do this only once in the pre-processing step, it is not an issue during the solver run time. Next, we further eliminate unknowns from using linear equations from . This is done on-line but it is fast since solving a small linear system is easy. In more detail, our method performs the following steps. Offline: 1. Let and consider the elimination ideal .

2. Compute the generators G of . These contain unknowns from only, i.e. the unknowns appearing in the linear equations . Online: 3. Rewrite the linear equations in the unknowns as , where M is a coefficient matrix and the vector contains all unknowns from .

4. Compute a null space basis N of M and re-parametrize the unknowns . If the rank of M is , i.e. the equations in are linearly independent, Y would contain new unknowns. Note that if all input equations in F were homogeneous, we could set one of the unknowns in Y to 1 (assuming it is non-zero) and then .

5. Substitute into the generators G of the elimination ideal .

6. Solve the new system of polynomial equations G(Y ) = 0 (e.g. using the Gr¨obner basis method and the precomputed elimination template for G(Y ) = 0 obtained by using the automatic generator [3]).

The main difference between our elimination strategy and the elimination strategies used before in minimal solvers (see Section 2.1.2) is that the previous strategies substitute the parametrization directly into the input nonlinear equations . This results in polynomial equations in unknowns . On the other hand, the new method eliminates unknowns from the non-linear equations and creates a system G(Y ) = 0 in k unknowns in the pre-processing step. We will show on several important problems from computer vision that solving the system G(Y ) = 0, instead of the system , is more efficient.

Before presenting our new strategy on more complicated problems from computer vision, we illustrate the key ideas of our strategy on a simpler, but still representative, example. In Appendix 5.2 we will show that the problem of estimating a 3D planar homography with unknown focal length, i.e. the projection matrix with unknown focal length, from planar points leads to a system of polynomial equations with the same structure as in this illustrative example.

2.2. Example

Let us consider the following system of nine homogeneous polynomial equations in ten unknowns X = . There are seven linear homogeneous equations in , namely

and two order equations in

Using the notation from Section 2.1.3 we have

We proceed as follows:

2. Compute the generator of the principal ideal . This is a polynomial of degree four:

The polynomial in G can be computed in the off-line pre-processing phase using the following code in the computer algebra system Macaulay2 [43]:

3. Rewrite seven linear equations form (5) as M h = 0, where and M = is 79 coefficient matrix.

4. Use a null space basis of M to reparametrize the unknowns from with two unknowns as

Since the input equations are homogeneous, we set (assuming ).

5. Substitute the new parametrization (9) into the generator (8).

6. Solve the resulting equation in one unknown .

7. Use the solutions for to recover solutions for using (9).

8. Extend the solutions for to solutions for X by substituting solutions to .

In this case, our elimination strategy generates one equation of degree four in one unknown.

On the other hand, the elimination strategy described in Section 2.1.2 generates two equations in two unknowns. More precisely, the strategy from Section 2.1.2 would substitute parametrization (9) directly into two equations from (6). This results in two equations in two unknowns and w. Solving this system of two equations in two unknowns in the on-line phase takes more time than solving a single quadratic equation.

3. Applications

3.1. f+E+f relative pose problem

The first problem that we solve using our elimination strategy is that of estimating relative pose and the common unknown focal length of two cameras from six image point correspondences. This problem is also known as the 6pt focal length problem, or the f+E+f problem. The f+E+f problem is a classical and popular problem in computer vision with many applications, e.g., in structure-from-motion [4].

The minimal f+E+f problem has 15 solutions and it was first solved by Stew`enius et al. [20] using the Gr¨obner basis method. The solver of Stew`enius consists of three G-J eliminations of three matrices of size 1233, 1633 and 1833 and the eigenvalue computation for a 1515 matrix.

More recently, two Gr¨obner basis solvers for the f+E+f problem ware proposed in [21] and [3]. The solver from [21] performs SVD decomposition of a 3450 matrix and it uses special techniques for improving the numerical stability of Gr¨obner basis solvers. The Gr¨obner basis solver generated by the automatic generator [3] performs G-J elimination of a 3146 matrix and is, to the best of our knowledge, the fastest and the most numerically stable solver for the f+E+f problem.

All the state-of-the-art (SOTA) solvers exploit that the 33 fundamental matrix satisfies

where K = diag(f, f, 1) is the diagonal calibration matrix with the unknown focal length f and E is the essential matrix [1]. The essential matrix has rank 2 and satisfies the Demazure equations [44]

(also known as the trace constraint). In all SOTA solvers [20, 21, 3], the linear equations from the epipolar constraints

for six image point correspondences in two views are first rewritten in a matrix form

where M is a 69 coefficient matrix and f is a vector of 9 elements of the fundamental matrix F. For six (generic) image correspondences in two views, the coefficient matrix M has a three-dimensional null space. Therefore, the fundamental matrix can be parametrized by two unknowns as

where are matrices created from the three-dimensional null space of M and x and y are new unknowns.

We use the parametrization (14), the rank constraint for the fundamental matrix

and the trace constraint (11) for the essential matrix, together with (10) in the following form:

This results in ten third- and fifth-order polynomial equations in three unknowns x, y and . In (16) we set . We note that the trace constraint (16) can be simplified by multiplying it with .

The ten equations (15) and (16) in three unknowns x, y and w were used as the input equations in all SOTA Gr¨obner basis solvers to the f+E+f problem [20, 21, 3].

Note, that all SOTA solvers followed the elimination method described in Section 2.1.2 and they differ only in the method used for solving the final non-linear system .

Next, we present a new solver for the f+E+f problem created using our elimination strategy in Section 2.1.3. This strategy not only generates a more efficient solver, but it also reveals new interesting constraints on the fundamental matrices of cameras with unknown focal length.

Elimination ideal formulation For the f+E+f problem we start with the ideal generated by ten equations from the rank constraint (15) and the trace constraint (11) with the essential matrix (10).

Since the epipolar constraint (12) gives us linear equations in , we have . Hence the strategy presented in Section 2.1.3 will first eliminate the unknown focal length f.

To compute the generators of the elimination ideal , i.e. the elements that do not contain the focal length f, we use the following Macaulay2 [43] code:

The output tells us that the variety of G has dimension 6 and degree 15, and that G is the complete intersection of two hypersurfaces in , cut out by the cubic

and the quintic

The vanishing of (17) and (18), together with the equation for extracting the unknown focal length from the fundamental matrix [45] (see also Section 5.3.2) completely describe the f+E+f problem. Therefore we can formulate the following result.

Result 3.1 The zero set of (17) and (18) equals the space of all fundamental matrices F, i.e. the singular 33 matrices, that can be decomposed into , where K = diag(f, f, 1) for some non-zero and E is an essential matrix. By intersecting this variety with six hyperplanes given by the epipolar constraints (12) for six image point correspondences, we obtain up to 15 real solutions for the fundamental matrix (see Figure 1).

In our new efficient on-line solver for the f+E+f problem, we first use the linear equations from the epipolar constraint (12) for six image point correspondences to

Figure 2. Numerical stability: Logof the relative error of the focal length for the (a) f+E+f problem; EI-fEf solver (blue), Kukelova08 [3] (red). (b) E+f problem; EI-Ef solver (blue), Buj-nak09 [39] (red).

parametrize the fundamental matrix F with two new unknowns x and y (14). After substituting this parametrization into the two generators (17) and (18) of , we get two equations (of degree 3 and 5) in two unknowns x and y. By solving these two equations we get up to 15 real solutions for the fundamental matrix F.

These two equations in two unknowns can be solved either using a Sylvester resultant [10] or using the Gr¨obner basis method, which was used in all SOTA solvers for the f+E+f problem. The Gr¨obner basis solver for these two equations, generated using the automatic generator [3], performs G-J elimination of a 2136 matrix. This matrix contains almost less nonzero elements than the matrix from the smallest 3146 SOTA solver [3] that was also generated using the automatic generator, however for the original formulation with ten equations in three unknowns. The sparsity patterns of both these solvers are shown in Section 5.5.

Experiments Since the new solver for the f+E+f problem is algebraically equivalent to the SOTA solvers, we have evaluated the new f+E+f solver on synthetic noise free data only.

We studied the behavior of the new f+E+f solver (EIfEf) based on the new elimination strategy presented in Section 2.1.3 on noise-free data to check its numerical stability. We compared it to the results of the SOTA Gr¨obner basis solver Kukelova08 [3]. In this experiment, we generated 10000 synthetic scenes with 3D points distributed at random in a cube. Each 3D point was projected by two cameras with random but feasible orientation and position and with random focal length . Figure 2(a) shows of the relative error of the focal length f obtained by selecting the real root closest to the ground truth value . In the case of the new EI-fEf solver, the focal length f was extracted from the computed F using the formula presented in Section 5.3.2. However, any method for extracting the focal length from F, e.g. the SVD-based method [45], can be used here.

The new EI-fFf solver (blue) is slightly less stable than Kukelova08 (red). However, both solvers provide very stable results without larger errors and in the presence of noise and in real applications their performance is in fact equivalent. Moreover, the new solver is smaller and more efficient.

3.2. E+f 6pt relative pose problem

The second problem that we solve using the new elimination strategy is the problem of estimating relative pose of one calibrated and one up to focal length calibrated camera from six image point correspondences, i.e. the E+f problem.

The minimal E+f problem was first solved by Bujnak et al. [39] using Gr¨obner bases and the polynomial eigenvalue method. Their solver performs G-J eliminations on a 2130 matrix and the eigenvalue computation for a 99 matrix.

For the E+f problem the first camera is calibrated up to an unknown focal length and the second camera is fully calibrated. Therefore, the relationship between the essential and the fundamental matrix has the form

where K = diag(f, f, 1) is a diagonal calibration matrix of the first camera, containing the unknown focal length f. By substituting this relationship into the trace constraint for the essential matrix (11), and setting Q = K K, we obtain

The SOTA solver [39] uses the elimination strategy from Section 2.1.2 and starts by rewriting the epipolar constraint (12) as M f = 0. Then, the fundamental matrix is parametrized by two unknowns x and y as in the f+E+f case (Eq. (14)). With this formulation, the rank constraint (15) and the trace constraint (20) result in ten third and fourth order polynomial equations in three unknowns x, y and . These ten equations are solved in [39] by using the automatic generator of Gr¨obner basis solvers [3].

Next we present a new solution to the E+f problem that uses the new elimination strategy from Section 2.1.3.

Elimination ideal formulation We start with the ideal generated by ten equations from the rank constraint (15) and the trace constraint (11), with the essential matrix (19). As for the f+E+f problem, the epipolar constraint (12) gives linear equations in . Therefore we again will eliminate only the unknown focal length f.

To compute the generators of the elimination ideal , i.e. the generators that do not contain f, we can use a similar Macaulay2 code as for the f+E+f problem, just by replacing line E = K*F*K with line E = F*K.

For the E+f problem, the variety of G has dimension 6 and degree 9 in and is defined by one cubic and three quartics (see Appendix 5.4).

In the online solver, the epipolar constraint (12) for six image point correspondences is used to parametrize the fundamental matrix F with two new unknowns x and y (14). This parametrization, applied to the four generators of the elimination ideal , gives four equations of degree three and four in two unknowns. We solve these four equations in two unknowns using the Gr¨obner basis method [3]. The Gr¨obner basis solver, generated using the automatic generator [3], performs G-J elimination of a 615 matrix. This matrix is much smaller than the elimination template matrix from the SOTA solver [39], which has the size 2130.

Experiments We studied the behavior of the new E+f elimination ideal based solver (EI-Ef) on noise-free data and compared it to the results of the SOTA Gr¨obner basis solver Bujnak09 [39].

We generated 10000 synthetic scenes with 3D points distributed at random in a cube. Each 3D point was projected by two cameras with random but feasible orientation and position. The focal length of the first camera was randomly drawn from the interval and the focal length of the second camera was set to 1, i.e. the second camera was considered as calibrated. Figure 2(b) shows of the relative error of the focal length f obtained by selecting the real root closest to the ground truth value . For the new EI-Ef solver, the focal length f was extracted from the computed F using the formula presented in Appendix 5.3.1.

The new EI-Ef solver (blue) is not only smaller but also slightly more stable than Bujnak09 [39] (red). Both solvers provide very stable results without larger errors.

3.3. E+f+k 7pt relative pose problem

The last problem that we will formulate and solve using the new elimination strategy presented in Section 2.1.3 is the problem of estimating the epipolar geometry of one calibrated camera and one camera with unknown focal length and unknown radial distortion, i.e. uncalibrated camera with radial distortion. We denote this problem by E+f+k.

A popular model for radial distortion is the one-parameter division model [46]. This is an undistortion model that can handle even quite pronounced radial distortions:

In this model are the homogeneous coordinates of the measured (and radially distorted) image points and is the distortion parameter. This model was used in the first 7pt minimal solution to the E+f+k problem presented in [47].

For the E+f+k problem, the epipolar constraint reads as

where are the homogeneous coordinates of corresponding ideally projected image points, i.e., points not corrupted by radial distortion [1]. Note that for the right camera we do not know the camera calibration parameters and we measure distorted image points. Therefore, to use these distorted image points in the epipolar constraint, we first need to undistort them using the model (21).

The epipolar constraint (22) together with the trace (20) and the rank constraint (15) form a quite complicated system of polynomial equations. Note that all equations in this system are non-linar and therefore the method from Section 2.1.2 cannot be directly applied.

In the SOTA solver [47], this system is first simplified by manually eliminating some unknowns. First, the authors set , which implies that their solver does not work for motions where . Next, they use the epipolar constraint (22) for six image point correspondences to eliminate six unknowns , which appear linearly in the epipolar constraint, from the equations. Then, the remaining equation from the epipolar constraint for the seventh image point correspondence, together with the trace (20) and the rank constraint (15), form a system of 11 (one quadratic, four and six degree) equations in four unknowns . Then, the equations are again manually simplified. They generate the elimination template by multiplying 11 input equations by a set of monomials such that the maximum degree of the monomials in the resulting equations is 8. The resulting elimination template has size 200231. The authors of this solver observed that by using automatic strategies from [3] or by further reducing the size of the elimination template, the numerical stability of their solver deteriorates. To improve the numerical stability of the final Gr¨obner basis solver, the authors further choose 40 monomials instead of necessary 19 for basis selection.

It can be seen that the E+f+k problem requires a very careful manual manipulation of the input equations to get a numerically stable solver. However, still, the final solver is quite large and not really useful in real applications.

Here we will show that using the new elimination strategy presented in Section 2.1.3 we can solve this problem ef-ficiently without the need for any special manipulation and treatment of input equations. Moreover, the final solver obtained using this new method is much more efficient and numerically stable than the SOTA solver [47].

Elimination ideal formulation Unfortunately, for the E+f+k problem the epipolar constraint (22) does not give us linear equations. Therefore, we can’t directly apply the method presented in Section 2.1.3. However, in this case we can easily linearize the equations from the epipolar constraint (22).

The epipolar constraint (22) contains monomials .

To linearize the equations (22) we set

Now the equations from (22) can be seen as linear homogeneous equations in the 12 unknowns .

Another view on this linearization is that the distorted image points are lifted to 4D space and the fundamental matrix F is enriched by one column to

where is the column of F. The 34 fundamental matrix (26) was introduced in [48] and is known as the one-sided radial distortion matrix. With this matrix, the epipolar constraint (22) can be written as

For the E+f+k problem, our method starts with the ideal generated by 13 equations, i.e. three equations from the constraints (23)-(25), the rank constraint (15) and the nine equations from the trace constraint (11), with the essential matrix of the form (19). These 13 equations form the set from (1) in our elimination strategy.

In this case, the “lifted” epipolar constraint (27) gives us linear equations in 12 elements of the radial distortion fundamental matrix (26), i.e. linear equations in . Therefore, for the E+f+k problem, we use the new elimination strategy to eliminate two unknowns, the focal length f and the radial distortion parameter , i.e. .

, i.e. the generators that do not contain f and , we can use a similar Macaulay2 code as for the E+f problem. We only need to replace the first line with

R = QQ[f,k,f11,f12,f13,f21,f22,f23,f31,f32, f33,y13,y23,y33];

and add one additional line at the end

For the E+f+k problem the variety of Gu has dimension 7 and degree 19 in . In addition to the three quadrics of the form , the ideal generators for Gu are

Figure 3. Numerical stability E+f+k problem : (a) Logrelative error of the focal length (b) Logof the relative error of the radial distortion; EI-Efk solver (blue), Kuang14 [47](red).

two cubics and nine quartics, i.e. altogether 14 polynomials. Although this system of 14 polynomial equations looks quite complex it is much easier to solve than the original system with and f that was used in the SOTA [47].

The “lifted” epipolar constraint (27) for seven general image point correspondences can be rewritten as , where M is 712 coefficient matrix and is a 121 vector containing the elements of the one-sided distortion fundamental matrix . This means that the one-sided distortion fundamental matrix can be parametrized by four new unknowns and , using the 5-dimensional null space of M, as

Substituting (28) into the 14 generators of the elimination ideal gives 14 equations in four unknowns. We solve these equations using the Gr¨obner basis method [3]. The Gr¨obner basis solver, generated using the automatic generator [3], performs G-J elimination of a 5170 matrix. This matrix is much smaller than the elimination template matrix from the SOTA solver [47], which has the size 200231. Moreover, the new solver doesn’t require an application of the methods for improving numerical stability of the Gr¨obner basis solver that were used in [47].

After solving 14 equations in four unknowns , we reconstruct solutions for F using (28), solutions for using (23), and solutions for f using the formula presented in Appendix 5.3.1.

Experiments We first studied the numerical stability of the new E+f+k solver (EI-Efk) on noise-free data and compared it to the results of the SOTA Gr¨obner basis solver Kuang14 [47]. In this experiment, we generated 10000 synthetic scenes in the same way as in the E+f experiment and image points in the first camera were corrupted by radial distortion following the one-parameter division model. The radial distortion parameter was drawn at random from the interval . Figure 3(a) shows of the relative error of the focal length f obtained by selecting the real root closest to the ground truth value . Figure 3(b) shows of the

Figure 4. Comparison of the new EI-Efk solver (blue) with the SOTA Kuang14 [47] solver (red). Boxplots of estimated for different noise levels and

relative error of the radial distortion obtained by selecting the real root closest to the ground truth value . In this case the new EI-Efk solver (blue) is not only significantly smaller but also significantly more stable than the SOTA Gr¨obner basis solver Kuang14 [47] (red). What is really important for real applications is that the new EI-Efk solver provides very stable results without larger errors, while for the SOTA solver Kuang14 [47] we observe many failures.

Next, Figure 4 shows the results of experiments with noise simulation for the E+f+k problem. We show the estimated radial distortion parameters for the ground truth radial distortion and 200 runs for each noise level. We compared our new E+f+k solver with the SOTA Kuang solver [47]. Figure 4 shows resuts by MATLAB boxplot. In the presence of noise, our new EI-Efk solver (blue) gives similar or even better estimates than the SOTA solver Kuang14 [47] for which we observed more failures (crosses).

3.4. Computational complexity

Here, we show a comparison of the computational effi-ciency of the new elimination-based solvers (EI-fEf, EI-Ef, EI-Efk) and the SOTA solvers [3, 39, 47]. Since we do not have comparable implementations of the SOTA solvers, we compare the sizes of the G-J eliminations (QR decompositions) performed by these solver. G-J elimination is one of the most time consuming steps for all solvers. The comparison of sizes is reported in the Tab. 1. The last row of this table displays the ratio of the number of non-zero elements of the template matrices of SOTA solvers () and the number of non-zero elements of the template matrices of our new elimination-based solvers ().

Table 1. New EI solvers are much smaller than SOTA solvers. Lines SOTA vs EI (new) show that clever solvers eliminate much smaller matrices. They also manipulate much fewer numbers. See ratios of non-zero numbers in the SOTA ((

4. Conclusion

We have presented a new insight into minimal solver construction based on elimination theory. By eliminating separately linear and non-linear equations and combining that later, we were able to generate much smaller solvers than before, see Tab. 1. We also generated an interesting new constraint, Eq. (18), on partially calibrated camera pairs. Our method was first motivated by the idea of exploiting linear equations (1) of our systems but we also demonstrated that it can produce efficient solvers (Sec. 3.3) by linearizing fully non-linear situations.

5. Appendix

This appendix includes (1) additional details on the Elimination Theorem (in Sec. 5.1), (2) derivation of the constraints on the projection for planar scenes by cameras with unknown focal length (in Sec. 5.2), (3) details of focal length extraction (in Sec. 5.3), (4) detailed presentation of the generators of E+f problem (in Sec. 5.4), and (5) results on the solvers’ sparsity (in Sec. 5.5).

5.1. Details on the elimination theorem

Here we provide additional details for

Theorem 5.1 (Elimination theorem [10]) Let

be an ideal and let G be a Gr¨obner basis of I with respect to the lexicographic monomial order where . Then, for every , the set is a Gr¨obner basis of the l-th elimination ideal .

See [10] for a full account of the theory.

The ring stands for all polynomials in n unknowns with complex coefficients. In computer vision applications, however, coefficients of polynomial systems are always real (in fact, rational) numbers and our systems consist of a finite number s of polynomial equations .

The ideal generated by s polynomials (generators) is the set of all polynomial linear combinations of the polynomials . Here the multipliers are polynomials. All elements in the ideal I evaluate to zero (are satisfied) at the solutions to the equations .

The Gr¨obner basis of an ideal I is a particularly convenient set of the generators of I, which can be used to find solutions to the original system in an easy way. For instance, for linear (polynomial) equations, a Gr¨obner basis of the ideal generated by the linear polynomials is obtained by Gaussian elimination. After Gaussian elimination, equations appear in a triangular form allowing one to solve for one unknown after another. This pattern carries on in a similar way to (some) Gr¨obner bases of general polynomial systems and thus it makes Gr¨obner bases a convenient tool for solving general polynomial systems.

Algorithmic construction of Gr¨obner bases relies on an ordering of monomials to specify in which order to deal with monomials of a polynomial. Lexicographic monomial order (LEX) is a particularly convenient order, which can be used to produce Gr¨obner bases that are in the triangular form. LEX orders monomials as words in a dictionary. An important parameter of a LEX order (i.e. ordering of words) is the order of the unknowns (i.e. ordering of letters). For instance, monomial when x > y > z (i.e. xyyz is before xyzz in a standard dictionary). However, when x < y < z, then . We see that there are n! possible LEX orders when dealing with n unknowns.

The set contains all the polynomials in Gr¨obner basis G that contain only unknowns . For instance, if G is a Gr¨obner basis in the triangular form, then contains polynomials in one, two, ..., l unknowns.

The polynomials generate the elimination ideal , containing all polynomials from I that use the unknowns only. Hence, for each of n! orderings, we get n elimination ideals .

5.2. 3D planar homograpy with unknown focal length

We assume that a planar object (say, simply a plane) is observed by an unknown camera with the projection matrix [1]

where K = diag(f, f, 1) is the calibration matrix with the unknown focal length is the unknown rotation, and the unknown translation.

Without loss of generality, we assume that the plane is defined by z = 0, i.e. all 3D points with homogeneous coordinates have the coordinate . Then, the image points and the corresponding 3D points are related by

where are unknown scalars, , and H = is a homography matrix that has the form

where is the column of the projection matrix P (29).

Next, from the projection equation (30), we eliminate the scalar values . This can be done by multiplying (30) by the skew symmetric matrix [1] to get

The matrix equation (32) contains three polynomial equations, two of which are linearly independent. This means that we need at least 3.5 2D 3D point correspondences to estimate the unknown homography H, because H has 7 degrees of freedom: three parameters for the rotation, three parameters for the translation and also the focal length.

For the 3.5 point correspondences, matrix equation (32) results in seven linearly independent linear homogeneous equations in nine elements of the homography matrix H.

Moreover, we have here two additional polynomial constraints on elements of H. For the first two columns of the rotation matrix R, there holds

This means that the elements of the first two columns of the homography matrix (31) satisfy

where w = 1/f.

Hence, estimating 3D planar homography with unknown focal length results in seven linear homogeneous equations and two non-linear homogeneous equations in X = . This system of nine homogeneous equations has the same form as that presented in Section 2.2. Therefore this system can be efficiently solved using the new elimination strategy presented in Section 2.1.3. This strategy results in solving one fourth-degree equation in one unknown (see Section 2.2).

5.3. Extraction of the focal length

In this section we present formulas for extracting the focal length from a given fundamental matrix F for two cases

1. E = F K

2. E = K F K

where K = diag(f, f, 1) is a diagonal calibration matrix. Unlike most of the existing formulas and methods for extracting the focal length from the fundamental matrix F, the presented formulas contain directly elements of the fundamental matrix. They don’t require an SVD decomposition of the fundamental matrix or computation of the epipoles.

5.3.1 E+f problem

Here we will assume that the principal points [1] are at the origin (which can be always achieved by shifting the known principal points) and use the recent result [49, Lemma 5.1] which we restate in our notation:

Lemma 5.2 Let F be a fundamental matrix of the form that satisfies E = F K. Then there are exactly two pairs of essential matrix and focal length (X = E, f) and (X = diag. The positive f is recovered from by the following formula

5.3.2 f+E+f problem

To derive formulas for the extraction of f from F computed from images with the same unknown focal length, we follow methods developed in [49]. In this case, the result is the following formula for , namely: quantity divided by

, which can be obtained by the following Macaulay2 code

5.4. The elimination ideal for the E+f problem

We consider the E+f problem from Section 3.2, i.e. the problem of estimating epipolar geometry of one calibrated and one up to focal length calibrated camera. Here, in this case

where K = diag(f, f, 1) is a diagonal calibration matrix for the first camera, containing the unknown focal length f. Here, F is the 33 fundamental matrix and E is the 33 essential matrix [1]

For the E+f problem, we have the ideal generated by ten equations, one cubic from the rank constraint

and nine polynomials from the trace constraint

where Q = K K.

For this problem, the new elimination strategy from Section 2.1.3 leads to computing the generators of the elimination ideal , i.e. the generators that do not contain f. To compute these generators we can use the following Macaulay2 [43] code:

degree 9 in and is defined by one cubic and three quartics. It can be verified that these four polynomials correspond to the four maximal minors of the 34 matrix:

5.5. Sparsity patterns of solvers

Here, we show a comparison of the sparsity patterns of our new elimination-based solvers (EI-fEf, EI-Ef, EI-Efk) and of the SOTA solvers [3, 39, 47].

Figure 5 shows the sparsity patterns of the (a) state-of-the-art (SOTA) 3146 Kukelova08 [3] solver for the f+E+f problem and (b) the new 2136 EI-fEf solver for this problem. In this case the new EI-fEf solver is not only smaller but also sparser. The ratio of the number of non-zero elements of the 3146 template matrix of the SOTA solver Kukelova08 [3] () and the number of non-zero elements of the 2136 matrix of the EI-fEf solver () is 3.

Figure 6 shows the sparsity patterns of the (a) SOTA 2130 Bujnak09 [39] solver and (b) the new 615 EI-Ef solver for the E+f problem. Here the ratio of the number of non-zero elements of the template matrix of the SOTA solver [39] and the number of non-zero elements of the template matrix of our new EI-Ef solver is 5.2.

Finally, Figure 7 shows the sparsity patterns of the (a) SOTA 200231 Kuang14 [47] solver and (b) the new 5170 EI-Efk solver for the E+f+k problem. Here, the ratio is approximately 2.8.

Acknowledgement

Z. Kukelova was supported by The Czech Science Foundation Project GACR P103/12/G084. T. Pajdla was supported

Figure 5. Sparsity patterns for the solvers to the f+E+f prob- lem: (a) state-of-the-art 3146 Kukelova08 [3] solver (b) the new 2136 EI-fEf solver.

Figure 6. Sparsity patterns for the solvers to the E+f problem: (a) state-of-the-art 2130 Bujnak09 [39] solver (b) the new 6EI-Ef solver.

Figure 7. Sparsity patterns for the solvers to the E+f+k prob- lem: (a) state-of-the-art 200231 Kuang14 [47] solver (b) the new 5170 EI-Efk solver.

by the EU-H2020 project LADIO No. 731970. J. Kileel and B. Sturmfels were supported by the US National Science Foundation (DMS-1419018) and the Max-Planck Institute for Mathematics in the Sciences, Leipzig, Germany.

References

[1] R. Hartley and A. Zisserman. Multiple View Geometry in Computer Vision. Cambridge, 2nd edition, 2003. 1, 2, 4, 7, 9, 10

[2] D. Nist´er. An efficient solution to the five-point rel- ative pose problem. IEEE Transactions on Pattern Analysis and Machine Intelligence, 26(6):756–770, June 2004. 1, 2

[3] Z. Kukelova, M. Bujnak, and T. Pajdla. Automatic generator of minimal problem solvers. In ECCV 2008

– European Conference on Computer Vision, volume 5304 of Lecture Notes in Computer Science. Springer, 2008. 1, 2, 3, 4, 5, 6, 7, 8, 11

[4] N. Snavely, S. M. Seitz, and R. Szeliski. Modeling the world from internet photo collections. Int. J. Computer Vision, 80(2):189–210, 2008. 1, 4

[5] D. Scaramuzza and F. Fraundorfer. Visual odometry. IEEE Robot. Automat. Mag., 18(4):80–92, 2011. 1

[6] J. Heinly, J. L. Sch¨onberger, E. Dunn, and J.-M. Frahm. Reconstructing the world in six days. In CVPR – IEEE Conference on Computer Vision and Pattern Recognition, pages 3287–3295, 2015. 1

[7] T. Sattler, B. Leibe, and L. Kobbelt. Efficient and effective prioritized matching for large-scale imagebased laocalization. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2016. 1

[8] M. A. Fischler and R. C. Bolles. Random sample con- sensus: A paradigm for model fitting with applications to image analysis and automated cartography. Commun. ACM, 24(6):381–395, 1981. 1

[9] R. Raguram, O. Chum, M. Pollefeys, J. Matas, and J.-M. Frahm. USAC: A universal framework for random sample consensus. IEEE Trans. Pattern Analysis Machine Intelligence, 35(8):2022–2038, 2013. 1

[10] D. A. Cox, J. Little, and D. O’Shea. Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra. Springer, 2015. 1, 2, 5, 9

[11] N. Courtois, A. Klimov, J. Patarin, and A. Shamir. Efficient algorithms for solving overdefined systems of multivariate polynomial equations. In Advances in Cryptology - EUROCRYPT 2000, International Conference on the Theory and Application of Cryptographic Techniques, pages 392–407, 2000. 1

[12] J. A. Grunert. Das Pothenotische Problem in erweiterter Gestalt nebst ¨uber seine Anwendungen in der Geod¨asie. Grunerts Archiv f¨ur Mathematik und Physik, 1:238–248, 1841. 1

[13] R. Sturm. Das Problem der Projektivit¨at und seine Anwendung auf die Fl¨achen zweiten Grades. Math. Annalen, 1:533–573, 1869. 1

[14] H. Longuet-Higgins. A computer algorithm for re- constructing a scene from two projections. Nature, 293(5828), 1981. 1

[15] H. Longuet-Higgins. A method of obtaining the rel- ative positions of four points from three perspective projections. Image Vision Computing, 10:266–270, 1992. 1

[16] D. Nist´er and F. Schaffalitzky. Four points in two or three calibrated views: Theory and practice. Int. J. Computer Vision, 67(2):211–231, 2006. 1

[17] L. Quan, B. Triggs, and B. Mourrain. Some results on minimal euclidean reconstruction from four points. Journal of Mathematical Imaging and Vision, 24(3):341–348, 2006. 1

[18] Z. Kukelova and T. Pajdla. A minimal solution to the autocalibration of radial distortion. In CVPR – IEEE Conference on Computer Vision and Pattern Recognition, 2007. 1

[19] D. Nist´er and H. Stew´enius. A minimal solution to the generalised 3-point pose problem. Journal of Mathematical Imaging and Vision, 27(1):67–79, 2007. 1

[20] H. Stew´enius, D. Nist´er, F. Kahl, and F. Schaffalitzky. A minimal solution for relative pose with unknown focal length. Image and Vision Computing, 26(7):871– 877, 2008. 1, 2, 4

[21] M. Byr¨od, K. Josephson, and K. ˚Astr¨om. Improving numerical accuracy of Gr¨obner basis polynomial equation solver. In ICCV – IEEE International Conference on Computer Vision, pages 1–8, 2007. 2, 4

[22] M. Byr¨od, K. Josephson, and K. ˚Astr¨om. A columnpivoting based strategy for monomial ordering in numerical Gr¨obner basis calculations. In ECCV – European Conference on Computer Vision, volume 5305, pages 130–143. Springer, 2008. 2

[23] S. Ramalingam and P. F. Sturm. Minimal solutions for generic imaging models. In CVPR – IEEE Conference on Computer Vision and Pattern Recognition, 2008. 2

[24] L. Kneip, D. Scaramuzza, and R. Siegwart. A novel parametrization of the perspective-three-point problem for a direct computation of absolute camera position and orientation. In CVPR – IEEE Conference on Computer Vision and Pattern Recognition, pages 2969–2976, 2011. 2

[25] O. Naroditsky and K. Daniilidis. Optimizing polyno- mial solvers for minimal geometry problems. In ICCV –IEEE International Conference on Computer Vision, pages 975–982, 2011. 2

[26] R. Hartley and H. Li. An efficient hidden variable approach to minimal-case camera motion estimation.

IEEE Transactions on Pattern Analysis and Machine Intelligence, 34(12):2303–2314, 2012. 2

[27] F. Camposeco, T. Sattler, and M. Pollefeys. Minimal solvers for generalized pose and scale estimation from two rays and one point. In ECCV – European Conference on Computer Vision, pages 202–218, 2016. 2

[28] L. Kneip, R. Siegwart, and M. Pollefeys. Finding the exact rotation between two images independently of the translation. In ECCV – European Conference on Computer Vision, pages 696–709, 2012. 2

[29] Z. Kukelova, J. Heller, M. Bujnak, A. W. Fitzgibbon, and T. Pajdla. Efficient solution to the epipolar geometry for radially distorted cameras. In ICCV – IEEE International Conference on Computer Vision, pages 2309–2317, 2015. 2

[30] J. A. Hesch and S. I. Roumeliotis. A direct leastsquares (DLS) method for PnP. In ICCV – IEEE International Conference on Computer Vision, pages 383– 390, 2011. 2

[31] V. Larsson and K. ˚Astr¨om. Uncovering symmetries in polynomial systems. In ECCV – European Conference on Computer Vision, pages 252–267, 2016. 2

[32] Ch. Aholt, B. Sturmfels, and R. R. Thomas. A Hilbert scheme in computer vision. Canadian Journal of Mathematics, 65:961–988, 2013. 2

[33] Ch. Aholt and L. Oeding. The ideal of the trifocal variety. Math. of Computation, 83:2553–2574, 2014. 2

[34] S. Agarwal, H. Lee, B. Sturmfels, and R. R. Thomas. On the existence of epipolar matrices. Int. J. Computer Vision, 121:403–415, 2017. 2

[35] M. Joswig, J. Kileel, B. Sturmfels, and A. Wagner. Rigid multiview varieties. Internat. J. Algebra Comput., 26(4):775–788, 2016. 2

[36] M. Trager, J. Ponce, and M. Hebert. Trinocular geom- etry revisited. Int. J. Computer Vision, 120(2):134– 152, 2016. 2

[37] Z. Kukelova. Algebraic Methods in Computer Vision. PhD thesis, Czech Technical Univ. Prague, 2013. 2

[38] H. Stewenius, C. Engels, and D. Nist´er. Recent de- velopments on direct relative orientation. ISPRS J. of Photogrammetry and Remote Sensing, 60, 2006. 2

[39] M. Bujnak, Z. Kukelova, and T. Pajdla. 3D reconstruc- tion from image collections with a single known focal length. In ICCV – IEEE International Conference on Computer Vision, 2009. 2, 5, 6, 8, 11

[40] B. Sturmfels. Solving Systems of Polynomial Equations. American Mathematical Society, CBMS Regional Conferences Series, No. 97, 2002. 2

[41] T. Becker and V. Weispfenning. Gr¨obner Bases: A Computational Approach to Commutative Algebra. Springer, 2003. 2

[42] M. Bujnak. Algebraic Solutions to Absolute Pose Problems. PhD thesis, Czech Technical University in Prague, 2012. 2

[43] D. Grayson and M. Stillman. Macaulay2, a software system for research in algebraic geometry. available at www.math.uiuc.edu/Macaulay2/. 3, 5, 11

[44] M. D´emazure. Sur deux probl`emes de reconstruction. Technical Report 882, INRIA, Rocquencourt, 1988. 4

[45] G. N. Newsam, D. Q. Huynh, M. J. Brooks, and H. P. Pan. Recovering unknown focal lengths in selfcalibration: An essentially linear algorithm and degenerate configurations. In Int. Arch. Photogrammetry and Remote Sensing, pages 575–580, 1996. 5

[46] A. W. Fitzgibbon. Simultaneous linear estimation of multiple view geometry and lens distortion. In CVPR – IEEE Conference on Computer Vision and Pattern Recognition, volume 1. IEEE, 2001. 6

[47] Y. Kuang, J. E. Solem, F. Kahl, and K. ˚Astr¨om. Minimal solvers for relative pose with a single unknown radial distortion. In CVPR – IEEE Conference on Computer Vision and Pattern Recognition, 2014. 6, 7, 8, 11

[48] J. H. Brito, C. Zach, K. K¨oser, M. J. Ferreira, and M. Pollefeys. One-sided radial fundamental matrix estimation. In BMVC – British Machine Vision Conference, 2012. 7

[49] J. Kileel, Z. Kukelova, T. Pajdla, and B. Sturmfels. Distortion varieties. http://arxiv.org/abs/1610.01860, 2016. 10

designed for accessibility and to further open science