Triangulation: Why Optimize?
Seong Hun Lee Javier Civera
cs.CV

Computational geometric methods for preferential clustering of particle suspensions.
    [
    December 30, 2023
====================================================================================






For decades, it has been widely accepted that the gold standard for two-view triangulation is to minimize the cost based on reprojection errors.
In this work, we challenge this idea.
We propose a novel alternative to the classic midpoint method that leads to significantly lower 2D errors and parallax errors.
It provides a numerically stable closed-form solution based solely on a pair of backprojected rays. Since our solution is rotationally invariant, it can also be applied for fisheye and omnidirectional cameras.
We show that for small parallax angles, our method outperforms the state-of-the-art in terms of combined 2D, 3D and parallax accuracy, while achieving comparable speed.






§ INTRODUCTION




Locating the 3D point given its projections in multiple views is called triangulation.
This classic yet fundamental problem in computer vision has immediate relevance to many applications, including visual odometry <cit.>, simultaneous localization and mapping (SLAM) <cit.> and structure-from-motion (SfM) <cit.>.
As such, achieving fast and accurate triangulation has been a goal of many research endeavors in the past decades.

For two views of known calibration and pose, the problem 
could be solved ideally if one finds the intersection of two backprojected rays corresponding to the same point.
However, the two rays are most likely skew due to noisy measurements and inaccurate camera model.
Since it is not obvious how to estimate the 3D position of the point from two skew rays, different methods have been proposed.
Mainly, they can be classified into three types:
(1) midpoint methods <cit.> that find the (weighted) midpoint of the common perpendicular between the two rays, (2) linear least squares methods <cit.>, and (3) optimal methods that “minimally” correct the two rays to make them intersect <cit.>.
Note that all these three types of methods produce solutions that minimize some cost function;
the (weighted) midpoint minimizes the (weighted) sum of squared distances to each ray, linear least squares methods minimize the algebraic errors, and optimal methods minimize a cost function based on either image reprojection errors <cit.> or angular reprojection errors <cit.>.
The most common cost functions are the L_1 norm (sum of magnitude), L_2 norm (sum of squares) and L_∞ norm (maximum) of the reprojection errors.



In this work, we suggest a different approach.
Instead of minimizing geometric or algebraic errors, we find a midpoint between a certain pair of points on each ray.
Like the classic midpoint method, our method takes the two rays as input.
Therefore, it is invariant to changes of camera rotation and applicable for perspective, fisheye and omnidirectional cameras.
Unlike the classic method, however, the two points on each ray are not necessarily on the common perpendicular.
We will see that our midpoint method bears a striking similarity to the classic method in the formulation, and yet it offers a significant performance gain in 2D and parallax accuracy. 
Compared to the optimal methods, our method yields much lower 3D errors at low parallax and similar 2D errors to those of L_2 and L_∞ optimal methods.
This motivates the question: 
In two-view triangulation, why optimize if there is a better way?

The main contributions of this paper are the following:


0em
    
  * We propose a novel method belonging to a group called the generalized weighted midpoint (GWM) method.
    We show that our method outperforms existing ones (including the classic midpoint method and the state-of-the-art optimal methods) in terms of combined 2D, 3D and parallax accuracy. 
    

    
  * Additionally, we propose a test of the adequacy (similar to the cheirality check) that identifies unreliable results and a weighting scheme that enhances 2D accuracy. 
    

    
  * We perform an extensive evaluation and analysis of various methods, revealing an intricate link between 3D accuracy and parallax estimation.
    This will provide an intuitive explanation of why our midpoint method performs better than the others.














§ RELATED WORK




One of the earliest works that addressed the two-view triangulation problem is <cit.> where the depth of a 3D point is estimated using simple algebra.
For robustness, later works mostly adopted geometric approaches, such as the midpoint method <cit.> and the minimization of the epipolar distance <cit.> or reprojection error <cit.>.
Among those, the last approach has become the de facto standard in computer vision <cit.>.

Optimal methods refer to those triangulation methods that minimize the cost based on reprojection errors.  
Assuming that the image measurements are independently perturbed by the noise in the same distribution of certain types, the optimal methods find the maximum likelihood (ML) solution.
For Gaussian and Laplacian distribution, the ML solution is to minimize the L_2 norm or L_1 norm of the reprojection errors, respectively <cit.>.
This can be found in closed form by solving a polynomial of degree six or eight <cit.>.
Alternatively, the L_2 solution can be obtained using iterative correction methods <cit.>.
While these iterative methods do not guarantee global optimality, they were shown to be faster and more stable.
For a uniform distribution, minimizing the L_∞ norm leads to the ML estimate for the lower bound of the noise <cit.>, and the solution is obtained by solving a quartic polynomial <cit.>.
Unlike the L_1 and L_2 cost, the L_∞ cost has a simple shape with a single minimum, but it is relatively more sensitive to noise and outliers <cit.>.


These optimal methods assume that the image measurement errors follow certain distributions. 
However, this assumption is neither justified nor likely <cit.>.
An equally (if not more) justified alternative is to assume that the noise model applies to the bearing measurements instead of the image.
For fisheye or omnidirectional cameras, the angular reprojection errors are more suitable than the image reprojection errors <cit.>.
Also, formulating the triangulation problem in terms of angular errors leads to much simpler ML solutions <cit.>.

Although the existing optimal methods can provide relatively good 3D results in many cases <cit.>, none of them are theoretically optimal in terms of 3D errors.
In fact, the discrepancy between 2D optimality and 3D accuracy has already been reported by Hartley and Sturm <cit.>.
They found that in Euclidean reconstruction, the midpoint and the linear least squares method achieve higher 3D accuracy than the L_1 or L_2 optimal methods, despite consistently (and sometimes significantly) larger 2D errors. 
In this work, we provide additional insights on this matter.
Furthermore, we show that a simple modification to the midpoint method can substantially reduce the 2D errors while maintaining 3D accuracy. 







§ PRELIMINARIES







Throughout the paper, we use bold letters for vectors and matrices, and light letters for scalars.
The Euclidean norm of a vector 𝐯 is denoted by ‖𝐯‖, and the unit vector by 𝐯=𝐯/‖𝐯‖.
The angle between two lines 𝐋_0 and 𝐋_1 is denoted by ∠(𝐋_0,  𝐋_1)∈[0, π/2].

Consider a 3D point observed by two cameras C_0 and C_1.
We define 𝐱_0 = [x_0, y_0, z_0]^⊺ and 𝐱_1 = [x_1, y_1, z_1]^⊺ as the unknown 3D coordinates of the point in the camera reference frame C_0 and C_1, respectively. 
Let 𝐑 and 𝐭 be the known rotation and translation between the two cameras, such that 𝐱_1=𝐑𝐱_0+𝐭.
Assuming that the camera calibration matrix 𝐊 is known, the normalized image coordinates 𝐟_0=[x_0/z_0, y_0/z_0, 1]^⊺ and 𝐟_1=[x_1/z_1, y_1/z_1, 1]^⊺ can be obtained by 𝐟_0 = 𝐊^-1𝐮_0 and 𝐟_1 = 𝐊^-1𝐮_1, where 𝐮_0 = (u_0, v_0,1)^⊺ and 𝐮_1=(u_1, v_1,1)^⊺ are the homogeneous pixel coordinates of the point observation in each frame.

In the ideal situation (Fig. <ref>a), the two backprojected rays intersect, satisfying the epipolar constraint <cit.>, i.e.,
𝐟_1·(𝐭×𝐑𝐟_0)=0.
Then, the intersection is given by 𝐱_1 = λ_0𝐑𝐟_0+𝐭 or 𝐱_1 = λ_1𝐟_1 for some scalar depth λ_0 and λ_1.
However, this rarely happens due to inaccuracies in the image measurements and the camera model.
Inferring a 3D point from two skew rays requires a nontrivial method.

Once the estimate of the 3D point (𝐱'_1) is obtained using some triangulation method, its accuracy can be evaluated in several ways.
One way is to compute the 3D error, i.e., e_3D = ‖𝐱'_1-𝐱_true‖.
Another way is to compute the 2D error (aka the reprojection error), i.e.,



    d_i 
        = ‖𝐊(𝐟_i-𝐟_i')‖
        =‖𝐊(𝐟_i-([0  0  1]  𝐱'_i)^-1𝐱'_i)‖  for   i=0,1,

where 𝐱'_0 = 𝐑^⊺(𝐱'_1-𝐭).
These two errors are illustrated in Fig. <ref>b.
Note that the 2D error represents the deviation from the measurement, whereas the 3D error represents the deviation from the ground truth.
Also, unlike the 3D error, the 2D error of a 3D point can be evaluated in different norms, e.g., L_1 norm (d_0+d_1), L_2 norm (√( d_0^2+d_1^2)) and L_∞ norm (max(d_0, d_1)).
Besides 2D and 3D accuracy, we can also evaluate the accuracy of the resulting parallax angle (see Fig. <ref>c). 
The parallax error is defined as follows:



    e_β 
        = |β_true-β'|
        = |∠(𝐱_true, 𝐱_true-𝐭)-∠(𝐱'_1,𝐱'_1-𝐭)|.

We define the “raw parallax” as the angle between the original backprojected rays:



    β_raw=∠(𝐑𝐟_0, 𝐟_1).

This gives a rough estimate of the parallax angle independently of the translation and the triangulation method. 




§ PROPOSED METHOD




 §.§ Generalized Weighted Midpoint (GWM) Method



R0.5
 

 
    < g r a p h i c s >

The classic midpoint and another example of the generalized weighted midpoint.





A GWM method consists of three steps:
(1) Given two backprojected rays corresponding to the same point, estimate the depth along each ray (λ_0, λ_1) using some method.
(2) Compute the 3D point on each ray at depth λ_0 and λ_1, i.e., 𝐭+λ_0𝐑𝐟_0 and λ_1𝐟_1 in C_1.
(3) Obtain the final estimate of the  3D point by computing their weighted average. 

The classic midpoint method <cit.> is one such type of method where the two points on each ray are the closest pair of points with the equal weight.
Fig. <ref> shows another possible example of the generalized weighted midpoint.




 §.§ Alternative Midpoint Method



We propose an alternative midpoint method that belongs to the generalized weighted midpoint method.
First, consider the case where the two backprojected rays happen to intersect (see Fig. <ref>a).
In this case, the most sensible solution is the point of intersection, and the corresponding depths along the rays can be obtained using the sine rule:   

    λ_0 = 
        sin(∠( 𝐟_1, 𝐭))/sin(∠(𝐑𝐟_0, 𝐟_1,))‖𝐭‖=‖𝐟_1×𝐭‖/‖𝐑𝐟_0×𝐟_1‖,   λ_1 = 
        sin(∠( 𝐑𝐟_0, 𝐭))/sin(∠(𝐑𝐟_0, 𝐟_1,))‖𝐭‖=‖𝐑𝐟_0×𝐭‖/‖𝐑𝐟_0×𝐟_1‖.

We use this formula to estimate the depths even when the two rays are skew.
Computing the 3D points on each ray at depth λ_0 and λ_1, respectively, we get

    𝐭+λ_0𝐑𝐟_0 = 𝐭+ ‖𝐟_1×𝐭‖/‖𝐑𝐟_0×𝐟_1‖𝐑𝐟_0   and  λ_1𝐟_1 = ‖𝐑𝐟_0×𝐭‖/‖𝐑𝐟_0×𝐟_1‖𝐟_1,

Taking the midpoint between these two points leads to

    𝐱'_1 = 1/2(𝐭+ ‖𝐟_1×𝐭‖/‖𝐑𝐟_0×𝐟_1‖𝐑𝐟_0+‖𝐑𝐟_0×𝐭‖/‖𝐑𝐟_0×𝐟_1‖𝐟_1).

Note that letting 𝐩 = 𝐑𝐟_0×𝐟_1, 𝐪 = 𝐑𝐟_0×𝐭 and 𝐫 =  𝐟_1×𝐭 allows us to write (<ref>) as

    λ_0 = ‖𝐫‖/‖𝐩‖,
          λ_1 = ‖𝐪‖/‖𝐩‖.

Interestingly, these are in a similar form to the depths given by the classic midpoint method[We provide the derivation in the Appendix.]:

    λ_mid_0 = 𝐩·𝐫/‖𝐩‖,
          λ_mid_1 = 𝐩·𝐪/‖𝐩‖.

The difference is in the numerator; 
(<ref>) has the magnitude of 𝐫 and 𝐪, whereas (<ref>) has their projection onto 𝐩.
As a result, we always get λ_0 ≥λ_mid_0 and λ_1 ≥λ_mid_1.
In most cases, this means that our midpoint will be located farther than the classic midpoint.
Fig. <ref> depicts one such example.
When we estimate that a point is located farther away from the cameras, it usually results in a lower estimate of the parallax angle, as will be shown in Section <ref>.





 §.§ Cheirality and Test of Adequacy




We say that the cheirality constraint <cit.> is violated when a triangulated point has negative depth(s).
This can happen for many reasons, such as spurious data association or the noise in the image point near the epipole.
Normally, it does not pose a serious problem because we can easily check the cheirality for each point and discard the bad ones.
For the classic midpoint method, this can be done by checking the signs of the depths given by (<ref>).

For our midpoint method, however, this is not possible because the depths given by (<ref>) are always positive. 
Fig <ref>a and  <ref>b illustrate the difference between the two methods. 
In our method, the depths alone cannot tell us whether or not the triangulation result is reliable.

Therefore, we use a different method to test the adequacy;
we discard the point correspondence if changing the sign of at least one depth to negative leads to a smaller distance between the two points on each ray, i.e.,



    =0mu=0mu=0mu
        ‖𝐭+λ_0𝐑𝐟_0-λ_1𝐟_1‖^2 
         ≥ min(
        ‖𝐭+λ_0𝐑𝐟_0+λ_1𝐟_1‖^2, 
        ‖𝐭-λ_0𝐑𝐟_0-λ_1𝐟_1‖^2,
        ‖𝐭-λ_0𝐑𝐟_0+λ_1𝐟_1‖^2
        )

For the classic midpoint method, letting λ_0 = |λ_mid_0| and λ_1 = |λ_mid_1|  gives effectively the same result as the cheirality check.
For example, (<ref>) holds in Fig <ref>a because the two points are closest when λ_0 = -|λ_mid_0| and λ_1 = -|λ_mid_1|.





 §.§ Inverse Depth Weighted Midpoint



The unweighted midpoint given by (<ref>) often entails disproportionate reprojection errors in the two images.
Fig. <ref>c shows an example.
Notice that the ray with a smaller depth tends to yield a larger reprojection error.
To compensate this imbalance, we propose to use the inverse depth <cit.> as a weight:



    𝐱'_1 
        = λ_0^-1(𝐭+λ_0𝐑𝐟_0)+λ_1^-1(λ_1𝐟_1)/λ_0^-1+λ_1^-1(<ref>)=‖𝐪‖/‖𝐪‖+‖𝐫‖(𝐭+‖𝐫‖/‖𝐩‖(𝐑𝐟_0+𝐟_1)).








§ EVALUATION RESULTS




We evaluate the following methods:
Lee and Civera's L_1, L_2 and L_∞ optimal angular methods (`L_1 ang', `L_2 ang', `L_∞ ang') <cit.>, Hartley and Sturm's L_1 and L_2 optimal methods (`L_1 img', `L_2 img') and linear methods (`DLT', `LinLS') <cit.>, Lindstrom's L_2 method with five iterations (`L_2 img (5 it.)') <cit.>, Níster's L_∞ method (`L_∞ img') <cit.>, the classic midpoint method (`Mid') <cit.>, and our method without and with the weighting (`Mid2', `wMid2').
The evaluation was performed on synthetic datasets generated as follows:
A set of 8×8 point clouds of 5,000 points each are generated with a Gaussian radial distribution 𝒩(0,(d/4)^2) where d is the distance from the world origin.
Each point cloud is centered at [0,0,d]^⊺ for d=2^n with n=-1,0,...,+6, and their image projections are perturbed by Gaussian noise 𝒩(0,σ^2) for σ=1,2...,8.
The size and the focal length of the images are 1,024^2 pixels and 512 pixel, respectively.
We have four configurations for the camera poses: (1) `orbital' - two cameras at [±0.5,0,0]^⊺ pointing at the point cloud center, (2) `lateral' - two cameras at [±0.5,0,0]^⊺ pointing at [0,0,∞]^⊺, (3) `forward' - two cameras at [0,0,±0.5]^⊺ pointing at the point cloud center, and (4) `diagonal' - two cameras at ±[√(3)/6,√(3)/6,√(3)/6]^⊺ pointing at [0,0,∞]^⊺.
The poses are slightly perturbed with uniform noise 𝒰(0, 0.01).
In total, the datasets provide over a million unique triangulation problems.




We aggregate the results in Fig. <ref> and <ref>.
Our observations agree with <cit.> in that:


0em
    
  * Generally, greater noise and lower parallax lead to larger 3D errors. 
    All methods yield almost equally low 3D errors for high-parallax points (> 4 deg). 

    
  * 2D and 3D errors are not well correlated.
    For example, LinLS and Mid perform best in 3D, but worst in 2D.




Additionally, we report the following findings:


0em
    
  * It is difficult to tell which method is the best in terms of 2D accuracy. 
    For example, L_1 methods yield the lowest 2D errors in L_1 norm, but relatively larger errors in L_2 and L_∞ norm.
    It is not obvious which norm is more important.
    That said, some methods can still perform consistently better than others;
    wMid2, L_2 and L_∞ methods consistently outperform LinLS, Mid, Mid2 and DLT in all 2D error criteria.
    

    
  * As shown in Fig. <ref>, LinLS and Mid are clearly more biased to overestimate the small parallax angles (< 4 deg). 
    This explains their relatively low 3D errors at low parallax. 
    Fig. <ref> illustrates a simplified example of this effect.
    

    
  * Our methods (Mid2 and wMid2) achieve the best overall accuracy in 3D + parallax. 



R0.45
 

 
    < g r a p h i c s >


Top: In this 2D example, the distance is proportional to (β/2).
Bottom: Relative impact is defined as |f(β)/g(β)| where f(β) and g(β) are the distance errors when the parallax is under- and overestimated by 0.5 deg, respectively.
Underestimating the parallax angle yields a larger distance error than overestimating it by the same degree, and especially more so for smaller parallax.






The last finding can be solely attributed to the low-parallax points, for which our methods show similar 3D accuracy to that of Mid, but much better parallax accuracy.
The latter can be explained by the fact that Mid2 always yields larger depths than Mid (as discussed in Section <ref>), which in effect lowers the estimated parallax angle on average (as shown in Fig. <ref>).
Since Mid tends to overestimate small parallax angles, this works to our advantage at low parallax.


 
The first plot in Fig. <ref> shows that Mid, LinLS and our two methods underestimate small parallax angles slightly less than the rest.
It is no coincidence that precisely these four methods achieve the best 3D accuracy; 
Fig. <ref> suggests that underestimating a small parallax angle has a severe impact on 3D accuracy.
At low parallax, it seems that our methods hit the sweet spot by (1) underestimating less than the optimal methods and DLT (thus achieving lower 3D errors) and (2) overestimating less than Mid and LinLS (thus achieving lower parallax errors). 




We suspect that some of the large 2D errors of Mid and LinLS at low parallax are related to their large parallax overestimation (see the second row of Fig. <ref>).
Large reprojection errors mean that the rays were corrected by a large amount, and pivoting two almost parallel rays (i.e., rays that correspond to a low-parallax point) will most likely increase the angle between them.
This link between the 2D and the parallax error could partially explain why our methods yield smaller 2D errors than Mid and LinLS for low-parallax points.

We also found that the 2D errors of wMid2 resemble those of L_∞ methods (see the bottom row of Fig. <ref>).
Note that the L_∞ optimal solutions yield the equal reprojection errors in the two views <cit.>.  
This suggests that our inverse depth weighting (Section <ref>) not only reduces the overall 2D errors, but also balances the reprojection errors in the two images.

Fig. <ref> compares the speed of each method. 
We included the test of the adequacy (Section <ref>) in both our methods.
All methods were implemented in C++ using the Eigen library <cit.>, compiled using GCC <cit.> with -O3 level optimization, and run on a laptop CPU (Intel i7-4810MQ, 2.8 GHz).
Although wMid2 is almost two times slower than Mid2 and three times than Mid, it is still at least ten times faster than the state-of-the-art method <cit.>.
Notice that even though (<ref>) and (<ref>) are similar, Mid is still almost twice faster than Mid2.
This happens for two reasons: 
First, we avoid computing the square root in (<ref>) by multiplying the numerator and denominator by ‖𝐩‖ and get λ_mid_0 = 𝐩^⊺𝐫/𝐩^⊺𝐩 and λ_mid_1 = 𝐩^⊺𝐪/𝐩^⊺𝐩.
Second, the test of adequacy described in Section <ref> takes longer than the standard cheirality check.
Nevertheless, given that the computational cost of triangulation is relatively small compared to other operations (e.g., point matching, pose estimation and structure refinement) <cit.>, our methods provide an excellent trade-off between speed and accuracy.





§ DISCUSSIONS






 §.§ On Optimality

A reasonable doubt is that the proposed (weighted) midpoint method may be in fact optimal in some error criterion. 
If we consider algebraic errors, it is obvious that our midpoint minimizes ‖𝐱'_1-𝐠(𝐑,𝐭,𝐟_0, 𝐟_1)‖ where 𝐠(·) is the vector-valued function given by the right-hand side of (<ref>) or (<ref>).
However, there are infinitely many other error functions for which the global minimum is 𝐠(𝐑,𝐭,𝐟_0, 𝐟_1), and it is hard to tell which one is geometrically meaningful (at least, we could not find it or disprove its existence).

Another reasonable doubt is that, among the GWM methods, there could be a better one than ours, since we have not proved the optimality of the proposed weighted midpoint.
A more fundamental question is then: 
What accuracy measure should we choose to define the optimality, knowing that there is a discrepancy between different types of errors, e.g., image/angular reprojection errors in different norms, 3D and parallax errors?
At least for our method, it seems that the weighting affects mostly the 2D errors, so there might be a certain weighting scheme that guarantees 2D optimality without compromising 3D and parallax accuracy.
A more elaborate error analysis and the comparison of different weightings within the GWM framework remain for future work.
For more discussion of the optimality in geometric vision, we refer to <cit.> and Appendix 3 of <cit.>.





 §.§ On Practical Implications

In many SfM pipelines, two-view triangulation is used to initialize the 3D map points prior to the bundle adjustment <cit.>.
Since the points with small parallax angles are associated with large 3D uncertainty, they are usually discarded.
This strategy is viable if there are enough correspondences with large parallax angles.
However, it is not ideal, as (1) low-parallax points can be useful for camera orientation estimation <cit.> and (2) if the goal itself is to reconstruct the scene with large depths compared to the baseline.
For problems such as reconstruction from small-baseline (or accidental) motions <cit.>, small parallax angles are quite common, so our method could be relevant.
Extending our method to multiple views for reconstructing low-parallax scenes would be an interesting future direction.





§ CONCLUSIONS


Triangulation from two views with known calibration and pose is an age-old problem in computer vision.
Existing methods formulate the problem as the minimization of some cost function, most commonly reprojection errors.
In this paper, we asked ourselves if this is really the best approach.
To this end, we proposed a novel variant of the classic midpoint method that does not minimize geometric or algebraic errors.

We found that all the existing methods we evaluated perform poorly at low parallax, producing large errors in either 2D, 3D or parallax. 
On the other hand, our midpoint method achieves very good overall accuracy.
We also showed that incorporating the inverse depth weighting can further reduce the 2D errors.
Although our method is not theoretically optimal, it provides, with speed and simplicity, a superior balance of 2D, 3D and parallax accuracy in practice.



§ ACKNOWLEDGEMENT

We thank Christopher Mei at Microsoft for valuable discussions.
This work was partially supported by the Spanish government (project PGC2018-096367-B-I00) and the Aragón regional government (Grupo DGA-T45_17R/FSE).



§ APPENDIX

In the following, we derive (<ref>).
In literature, the formula has been used many times without derivation <cit.>.
For the sake of completeness, we present the full derivations here.
In doing so, we will use the following properties of the dot product and the cross product operations:

    𝐚× (𝐚×𝐛)=𝐚 (𝐚·𝐛)-𝐛.  
    
        (𝐚×𝐛)·(𝐚×𝐜)=𝐛·𝐜-(𝐚·𝐛)(𝐚·𝐜). 
    
        (𝐚×𝐛)×(𝐚×𝐜)=(𝐚·(𝐛×𝐜))𝐚.

Next, we introduce the following lemma:


Consider two skew lines 𝐋_0(s_0) = 𝐜_0+s_0𝐦_0 and 𝐋_1(s_1) = 𝐜_1+s_1𝐦_1 in 3D space. 
Let 𝐭=𝐜_0-𝐜_1 and (𝐫_0, 𝐫_1) be the two points on each line that form the closest pair.
Then, 

    𝐫_0 = 𝐜_0+(𝐦_0×𝐦_1)·(𝐦_1×𝐭)/‖𝐦_0×𝐦_1‖^2𝐦_0

and 

    𝐫_1 = 𝐜_1+(𝐦_0×𝐦_1)·(𝐦_0×𝐭)/‖𝐦_0×𝐦_1‖^2𝐦_1.


Proof.   
In geometry, it is a well-known fact that the closest pair of points on two skew lines lie on the common perpendicular to both lines.
In other words, the vector 𝐫_0-𝐫_1 is perpendicular to both 𝐋_0 and 𝐋_1.
Therefore, for some scalar τ, 

    𝐫_0-𝐫_1 = τ(𝐦_0×𝐦_1).

Since point 𝐫_0 and 𝐫_1 are respectively located along 𝐋_0 and 𝐋_1, we can write 

    𝐫_0=𝐜_0+λ_0𝐦_0   and  𝐫_1=𝐜_1+λ_1𝐦_1.

for some scalar λ_0 and λ_1. 
Then, (<ref>) becomes

    𝐭+λ_0𝐦_0-λ_1𝐦_1=τ𝐧,

where 𝐧=𝐦_0×𝐦_1.
This makes a system of three equations (in each coordinate x, y and z) with three unknowns  λ_0, λ_1, and τ.
Removing τ from the equations leads to

    t_x+λ_0m_0x-λ_1m_1x/n_x
        =
        t_y+λ_0m_0y-λ_1m_1y/n_y

and

    t_y+λ_0m_0y-λ_1m_1y/n_y
        =
        t_z+λ_0m_0z-λ_1m_1z/n_z.

Note that 𝐭=[t_x, t_y, t_z]^⊺, 𝐧=[n_x, n_y, n_z]^⊺, 𝐦_0 = [m_0x, m_0y, m_0z]^⊺ and 𝐦_1 = [m_1x, m_1y, m_1z]^⊺.
From (<ref>) and (<ref>), we get

    λ_0 = λ_1(m_1xn_y-m_1yn_x)+t_yn_x-t_xn_y/m_0xn_y-m_0yn_x

and

    λ_0 = λ_1(m_1yn_z-m_1zn_y)+t_zn_y-t_yn_z/m_0yn_z-m_0zn_y.

Equating the right-hand sides of (<ref>) and (<ref>) leads to

    λ_1 = A-B/C-D,

where

    A   = (t_zn_y-t_yn_z)(m_0xn_y-m_0yn_x)
    
        B   = (t_yn_x-t_xn_y)(m_0yn_z-m_0zn_y)
    
        C   = (m_1xn_y-m_1yn_x)(m_0yn_z-m_0zn_y) 
       =(𝐦_1×𝐧)_z(𝐦_0×𝐧)_x,
    
        D   = (m_1yn_z-m_1zn_y)(m_0xn_y-m_0yn_x)
       = (𝐦_1×𝐧)_x(𝐦_0×𝐧)_z.

We can rearrange A-B into

    A-B = n_y𝐭·[ m_0yn_z-m_0zn_y; m_0zn_x-m_0xn_z; m_0xn_y-m_0yn_x ].

The latter term in the dot product of (<ref>) is equal to 𝐦_0×𝐧. 
Thus,

    A-B    = n_y𝐭·(𝐦_0×(𝐦_0×𝐦_1))
       (<ref>)=n_y𝐭·(𝐦_0(𝐦_0·𝐦_1)-𝐦_1)
       =n_y((𝐦_0·𝐦_1)(𝐦_0·𝐭)-𝐦_1·𝐭)
       (<ref>)=-n_y(𝐦_0×𝐦_1)·(𝐦_0×𝐭).

We can rearrange C-D into

    C-D   =(𝐦_1×𝐧)_z(𝐦_0×𝐧)_x-(𝐦_1×𝐧)_x(𝐦_0×𝐧)_z 
       =((𝐦_1×𝐧)×(𝐦_0×𝐧))_y 
       (<ref>)=((𝐧·(𝐦_1×𝐦_0))𝐧)_y 
       =(-‖𝐦_0×𝐦_1‖^2𝐧)_y 
       =-‖𝐦_0×𝐦_1‖^2 n_y.

Substituting (<ref>) and (<ref>) into (<ref>) gives

    λ_1 = (𝐦_0×𝐦_1)·(𝐦_0×𝐭)/‖𝐦_0×𝐦_1‖^2.

Finally, substituting (<ref>) into (<ref>) leads to (<ref>).
Equation (<ref>) is derived analogously. 




By substituting 𝐑𝐟_0 and 𝐟_1 into 𝐦_0 and 𝐦_1, we can use lemma <ref> to obtain λ_mid_0 and λ_mid_1 in (<ref>).


    λ_mid_0 = (𝐑𝐟_0×𝐟_1)·(𝐟_1×𝐭)/‖𝐑𝐟_0×𝐟_1‖^2,
          λ_mid_1 = (𝐑𝐟_0×𝐟_1)·(𝐑𝐟_0×𝐭)/‖𝐑𝐟_0×𝐟_1‖^2.

Letting 𝐩 = 𝐑𝐟_0×𝐟_1, 𝐪 = 𝐑𝐟_0×𝐭 and 𝐫 =  𝐟_1×𝐭 leads to (<ref>).

30


urlstyle

[eig()]eigen_library
Eigen Library v3.
<http://eigen.tuxfamily.org>.

[gcc()]gcc
GNU Compiler Collection (gcc).
<https://www.gnu.org/software/gcc>.

[Beardsley et al.(1994)Beardsley, Zisserman, and
  Murray]beardsley_midpoint1
P. A. Beardsley, A. Zisserman, and D. W. Murray.
Navigation using affine structure from motion.
In Eur. Conf. on Computer Vision, pages 85–96, 1994.

[Beardsley et al.(1997)Beardsley, Zisserman, and
  Murray]beardsley_midpoint2
P. A. Beardsley, A. Zisserman, and D. W. Murray.
Sequential updating of projective and affine structure from motion.
Intl. Journal of Computer Vision, 230 (3):0
  235–259, 1997.

[Civera et al.(2008)Civera, Davison, and
  Montiel]civera2008interacting
J. Civera, A. J. Davison, and J. M. M. Montiel.
Interacting multiple model monocular slam.
In IEEE Intl. Conf. on Robotics and Automation, pages
  3704–3709, 2008.

[Civera et al.(2008)Civera, Davison, and Montiel]inverse_depth
J. Civera, A. J. Davison, and J. M. M. Montiel.
Inverse depth parametrization for monocular SLAM.
IEEE Trans. Robot., 240 (5):0 932–945, 2008.
ISSN 1552-3098.
10.1109/TRO.2008.2003276.

[Delon and Rougé(2007)]delon2007small
Julie Delon and Bernard Rougé.
Small baseline stereovision.
Journal of Mathematical Imaging and Vision, 280
  (3):0 209–223, 2007.

[Forster et al.(2014)Forster, Pizzoli, and Scaramuzza]svo
C. Forster, M. Pizzoli, and D. Scaramuzza.
SVO: fast semi-direct monocular visual odometry.
In IEEE Intl. Conf. on Robotics and Automation, pages 15–22,
  2014.

[Harris(1987)]harris1
C. G. Harris.
Determination of ego-motion from matched points.
In Third Alvey Vision Conference, pages 189–192, 1987.

[Harris and Pike(1988)]harris2
C.G. Harris and J.M. Pike.
3D positional integration from image sequences.
Image and Vision Computing, 60 (2):0 87 – 90,
  1988.

[Hartley and Kahl(2007)]optimal_algorithms_hartley
R. Hartley and F. Kahl.
Optimal algorithms in multiview geometry.
In Asian Conf. on Computer Vision, pages 13–34, 2007.

[Hartley and Schaffalitzky(2004)]L_inf_hartley
R. Hartley and F. Schaffalitzky.
L_∞ minimization in geometric reconstruction problems.
In IEEE Conf. on Computer Vision and Pattern Recognition,
  pages 504–509, 2004.

[Hartley and Zisserman(2003)]hartley_book
R. Hartley and A. Zisserman.
Multiple View Geometry in Computer Vision.
Cambridge University Press, New York, NY, USA, 2 edition, 2003.

[Hartley(1994)]hartley_euclidean
R. I. Hartley.
Euclidean reconstruction from uncalibrated views.
In Proceedings of the Second Joint European - US Workshop on
  Applications of Invariance in Computer Vision, pages 237–256, 1994.

[Hartley and Sturm(1997)]hartley_triangulation
R. I. Hartley and P. Sturm.
Triangulation.
Comput. Vis. Image Underst., 680 (2):0
  146–157, 1997.

[Kanatani et al.(2008)Kanatani, Sugaya, and Niitsuma]kanatani
K. Kanatani, Y. Sugaya, and H. Niitsuma.
Triangulation from two views revisited: Hartley-Sturm vs. optimal
  correction.
In British Machine Vision Conference, pages 173–182, 2008.

[Kanazawa and Kanatani(1995)]kanazawa
Y. Kanazawa and K. Kanatani.
Reliability of 3-D reconstruction by stereo vision.
IECE Trans. Inf. & Syst., E78-D0 (10):0
  1301–1306, 1995.

[Ke and Kanade(2003)]L1_mle
Qifa Ke and Takeo Kanade.
Robust subspace computation using L1 norm.
School of Computer Science, Carnegie Mellon University Pittsburgh,
  PA, 2003.

[Lee and
  Civera(2019)]closed_form_optimal_triangulation_based_angular_errors
S. Lee and J. Civera.
Closed-form optimal two-view triangulation based on angular errors.
In IEEE Int. Conf. on Computer Vision, 2019.

[Lindstrom(2010)]lindstrom
P. Lindstrom.
Triangulation made easy.
In IEEE Conf. on Computer Vision and Pattern Recognition,
  pages 1554–1561, 2010.

[Longuet-Higgins(1981)]eight_point
H. C. Longuet-Higgins.
A computer algorithm for reconstructing a scene from two
  projections.
Nature, 2930 (5828):0 133–135, 1981.

[Mičušík and Pajdla(2006)]sfm_wideangle
B. Mičušík and T. Pajdla.
Structure from motion with wide circular field of view cameras.
IEEE Trans. Pattern Anal. Mach. Intell., 280
  (7):0 1135–1149, 2006.

[Mur-Artal et al.(2015)Mur-Artal, Montiel, and Tardós]orb-slam
R. Mur-Artal, J. M. M. Montiel, and J. D. Tardós.
ORB-SLAM: A versatile and accurate monocular SLAM system.
IEEE Trans. Robot., 310 (5):0 1147–1163,
  2015.
10.1109/tro.2015.2463671.

[Níster(2001)]nister_phd_thesis
D. Níster.
Automatic Dense Reconstruction from Uncalibrated Video
  Sequences.
PhD thesis, KTH, 2001.

[Nousias et al.(2019)Nousias, Lourakis, and Bergeles]nousias2019large
S. Nousias, M. Lourakis, and C. Bergeles.
Large-scale, metric structure from motion for unordered light fields.
In IEEE Conf. on Computer Vision and Pattern Recognition,
  2019.

[Oliensis(2002)]closed_form_oliensis
J. Oliensis.
Exact two-image structure from motion.
IEEE Trans. Pattern Anal. Mach. Intell., 240
  (12):0 1618–1633, 2002.

[Pirchheim et al.(2013)Pirchheim, Schmalstieg, and
  Reitmayr]pirchheim2013handling
C. Pirchheim, D. Schmalstieg, and G. Reitmayr.
Handling pure camera rotation in keyframe-based SLAM.
In IEEE International Symposium on Mixed and Augmented
  Reality, pages 229–238, 2013.

[Schönberger and Frahm(2016)]sfm_revisited
J. L. Schönberger and J. Frahm.
Structure-from-Motion revisited.
In IEEE Conf. on Computer Vision and Pattern Recognition,
  pages 4104–4113, 2016.

[Snavely et al.(2006)Snavely, Seitz, and Szeliski]snavely2006photo
N. Snavely, S. M. Seitz, and R. Szeliski.
Photo Tourism: Exploring photo collections in 3D.
ACM Trans. Graph., 250 (3):0 835–846, 2006.

[Yang et al.(2019)Yang, Fang, Zhao, and
  Deng]iteratively_reweighted_midpoint
K. Yang, W. Fang, Y. Zhao, and N. Deng.
Iteratively reweighted midpoint method for fast multiple view
  triangulation.
IEEE Robot. Autom. Lett., 40 (2):0 708–715,
  2019.

[Yu and Gallup(2014)]accidental_motion
F. Yu and D. Gallup.
3D reconstruction from accidental motion.
In IEEE Conf. on Computer Vision and Pattern Recognition,
  pages 3986–3993, 2014.




