In many 3D vision problems that range from camera calibration to robot navigation, it is required to represent 3D straight lines. Plücker coordinates are one of the most used formulations article:pottman:2001
. Vectors of Plücker coordinates are built-up by stacking both direction and moment of the respective lines (both 3D vectors), giving a 6D vector. Direction and moment vectors must satisfy the so-called Klein quadric constraint, that is, they need to be orthogonal to each other. In this paper we address the problem of estimating Pücker coordinates from general (unconstrained)vectors. This is called the Plücker correction problem.
For many reasons, especially when considering data with noise, it is hard to include the orthogonal constraint in the estimator (frequently, it requires non-linear procedures). Some authors do not consider this constraint on their methods or propose optimization techniques involving this constraint but, to avoid unnecessary computational effort, non-linear procedures are stopped before fulfilling the respective constraint. Several examples can be found in the literature, e.g. camera models (mapping between pixels and 3D straight lines) article:miraldo:2013 ; triangulation of 3D lines article:josephson:2008 ; structure-from-motion using lines Bartoli05 ; article:lemaire:2007 ; and 3D reconstruction of lines using a single image of a non-central catadioptric camera article:lemaire:2007 . Thus, to get Plücker coordinates on these cases, a Plücker correction needs to be applied.
The state-of-art method for solving the Plücker correction problem is due to Bartoli and Sturm at Bartoli05 , hereafter called BS method. They find the closest Plücker coordinates (in Euclidean sense) from an unconstrained six-dimensional vector by orthogonally projecting the input vector onto the Klein quadric (and so it verifies the orthogonal constraint). Their approach involves SVD decompositions and can be found in (Bartoli05, , p. 425). It should be noticed that in the description of the BS Plücker correction algorithm (Table 2 in Bartoli05 ) there is a typo in the entries of the matrix . The correct matrix is given by
We performed a detailed analysis of the proof of the BS method ((Bartoli05, , p. 425)) and propose some clarifications. For the proof of BS method the following identity is used:
where and are matrices of sizes and , respectively, with , and is an matrix with orthonormal columns (i.e . , but ). This identity is related with the invariance of the Frobenius norm for left multiplication by matrices with orthonormal columns.
However, (2) does not hold in general. It is valid if, and only if, the space of columns of coincides with the space of columns of , in which case , for some matrix of size . We notice that, in the proof of the BS method, this requirement is not explicitly mentioned.
Paper Wu15 also addresses the same problem. However they do not address the problem of the computational efficiency and the proposed approach is quite different from the method described in this paper. Moreover, even if the solution proposed in Wu15 is indeed a global minimum, there is no formal proof of such a fact. We shall recall that Lagrange multipliers yield only local minima, unless the existence of a global minimum is guaranteed.
1.1 Notations and Problem Definition
Column vectors are represented by bold small letters (e.g. for an -dimensional vector). Bold capital letters denote matrices (e.g. for an matrix). Regular small letters denoted zero dimensional elements (e.g. ). denotes the Frobenius norm for matrices or the -norm for vectors. Recall that, for any matrix , the Frobenius norm is given by and, for any vector , the -norm satisfies .
The Plücker coordinates of a 3D straight line can be represented by a six dimensional vector:
where are, respectively, the direction and moment of the line, verifying the Klein quadric constraint
see article:pottman:2001 .
Let and be given vectors in (not necessarily satisfying the orthogonality constraint), and assume that and denote vectors in . Mathematically, the Plücker correction problem corresponds to solving the nonlinear constrained optimization problem formulated as
While the objective function
is convex, the Klein quadric constraint is not. The optimization problem (5) belongs to a class of non-convex problems known in the literature as quadratically constrained quadratic programs. This means, in particular, that the existence of a global minimum may be a non trivial problem.
The main goals of this paper are to prove the existence of a global minimum for (5) and to give an explicit formula for computing such a minimum. Our approach is essentially based on the application of the classical Lagrange multipliers to the constrained problem (5). It does not involve singular value decompositions which turns it faster than the BS method. Our results are supported by mathematical proofs and numerical experiments. In addition our method is designed to deal with general -dimensional vectors and , and not exclusively with -dimensional vectors.
1.2 The Explicit Formula
Let us consider two general 3D vectors and , such that and . As it will be shown in the following sections, the proposed solution for the global minimum of (5) is given by
2 Plücker Correction using Lagrange Multipliers
Suppose again that and satisfy and , and let denote the Lagrangian function associated to (5). Some calculation yields
where the real number is the Lagrange multiplier, and are the aimed solutions. The partial derivatives of the Lagrangian with respect to , and are (for formulae of derivatives with respect to vectors see Lutkepohl96 )
where and are
Hence, the solutions of (17) are
Denoting, for ,
we know that and are the candidates to be a local minimum of the Lagrangian . However, since the gradient of the constraint involved in (5) annihilates at , one also needs to consider this non regular and non stationary point. Thus a local minimum of (5) must be attained at one of the following three points: , or .
Now we shall note the following facts.
The assumption of and being such that , with , guarantees, in particular, that cannot be , which ensures that and in (16) are well defined.
In (19), is always a real number because . Indeed,
Moreover, if , then and are orthogonal, which means that the global minimum of the objective function is attained at .
We end this section by showing that the objective function (6) satisfies
which proves that given in (20) is the local minimum where the objective function attains the smallest value. In Sec. 2.1 is is shown that is also the global minimum. Substituting and given in (16) in the objective function (6), some calculation yields the following real function depending on the real variable :
Now (21) can be rewritten as
Let us write the function in the form
with and .
In fact, because () satisfies the quadratic equation (17), one has . Therefore,
Finally, considering the equality and the fact that the quadratic function is non negative, for all real numbers , the second inequality of (23) arises as a consequence of the equivalences
2.1 Existence of a Global Minimum
Given a unit vector in the Euclidean space , consider the line through the origin with the direction of . Given another vector , it is well-known that the unique vector of that is closest to (in Euclidean sense) is the orthogonal projection of onto , which is the vector (see, for instance, (Meyer00, , p.435)).
Hence, the constrained problem (5
) corresponds to find orthonormal vectorsand minimizing the sum of the distance , between and the line through the origin with the direction of , with the distance , between and the line through the origin with the direction of (see Fig. 1). This means that (5) can be reformulated as
This formulation of the Plücker correction problem is apparently more complicated and less practical than (5), but it is helpful to show that (5) has in fact a global minimum. To see this, one just needs to observe that the constraints in (25) define a closed and bounded (and, consequently, compact) set in endowed with the Euclidean metric. Thus, by the classical Weierstrass theorem (see, for instance, (Luenberger08, , Appendix A.6)), there exists at least a global minimum. This proves that the analysis carried out in Section 2, using the Lagrange multipliers, guarantees that the objective function (6) attains a global minimum at
where is defined as in (19).
As far as we know, the cases when
rarely occur in practical problems of computer vision. However, it is worth to make some comments on this particular case.
If with , the solution of (5) is obviously and . If but then replacing by in the Lagrangian (9), we get a simpler expression. Finding the first order optimality conditions and solving them, one easily concludes that any pair of vectors of the form , with and arbitrary, gives a local minimum. The value of the Lagrangian at all these local minima is and does not depend on . Choosing, for instance, , it follows that the pair is a local minimum of (5). Using a similar argument to that of Sec. 2.1, this pair is also a global minimum.
Similarly, if , with , it can be concluded that provides also a global minimum for (5).
In this section, we compare the method based on the explicit formula (7), with the method of Bartoli and Sturm Bartoli05 , in terms of computational effort. Both methods were implemented using MATLAB. The codes will be available in the author’s website. We consider three different algorithms:
LMPC which denotes to the method derived in this paper;
BS which corresponds to the Bartoli and Sturm’s approach; and
BS-LSVD which denotes to the method proposed by Bartoli and Sturm, where the SVD is computed using closed-form techniques.
A detailed description of each algorithm is shown in A
To compare the methods we use the following procedure: we randomly generated unconstrained vectors (Klein quadric is not enforced). For each trial, we apply both Plücker correction algorithms to the respective six-dimensional vectors, storing the values of the corresponding objective functions (6). The results in terms of computational time required for each algorithm is shown in Tab. 1.
|Algorithm||For all trials||For each trial (median)|
Both LMPC and BS-LSVD methods can be implemented using only closed-form steps. However, while LMPC can be computed with a few steps (it only requires six lines of code), BS-LSVD requires more algebraic operations and takes a longer time to run, see Tab. 1. Indeed, BS-LSVD is about two times slower than LMPC. On the other hand, if we consider the original Bartoli and Sturm approach BS, which requires the computation of two singular value decompositions, it requires iterative steps and, as a result, this method is significantly slower. From our experiments, one can see that this method is more than six times slower than our approach.
3.1 Discussion of the Experimental Results
In many applications there are estimates of 3D line coordinates that are not obtained from points. 3D lines can be estimated from intersections of planes, for example, or they can be obtained from non-conventional sensors, such as non-central generic cameras. In non-central generic cameras calibration, it consists in estimating a 3D line for each pixel article:sturm:2011 ; article:miraldo:2013 . If we represent 3D lines using Plücker coordinates, the Klein quadric constraint has to be enforced — in this case this correction has to be applied to all and each pixel, which corresponds to call the Plücker correction algorithm a large number of times. Let us consider a camera system, with a standard image size of , which contains a total of pixels. For this case and from the experimental results, one can conclude that the Plücker correction step using the method proposed in this paper will be, at least, two times faster than Bartoli and Sturm method, which consists in saving more than eight seconds.
Other example is, using perspective cameras, the estimation of 3D lines from the intersections of the back-projecting planes from two or more images.
In this paper we have addressed the Plücker correction problem, by minimizing the Frobenius norm between the estimated and input vectors. By solving the corresponding optimization problem using Lagrange multipliers, a simple solution, that can be computed in closed-form with a very small number of operations, was proposed. Contrarily to the state-of-the-art method, we have proved theoretically that our method computes always the global minimum. In addition, the special cases (where a solution cannot be computed) are analysed. As the experimental results show, the proposed method is faster.
- (1) H. Lütkepohl, Handbook of Matrices, Jonh Wiley and Sons, 1996.
- (2) J. Nocedal and S. Wright, Numerical Optimization, Springer-Verlag New York, 1999.
- (3) C. D. Meyer, Matrix Analysis and Applied Linear Algebra, SIAM, Philadelphia, 2000.
- (4) H. Pottmann and J. Wallner, “Computational Line Geometry”, Springer-Verlag Berlin Heidelberg, 2001.
- (5) R. Hartley and A. Zisserman, Multiple View Geometry in Computer Vision, 2nd Edition, Cambridge University Press, 2004.
- (6) A. Bartoli and P. Sturm, “Structure-from-Motion Using Lines: Representation, Triangulation and Bundle Adjustment”, Computer Vision and Image Understanding, 2005.
- (7) D. Lanman, M. Wachs G. Taubin, and F. Cukierman, “Reconstructing a 3D Line from a Single Catadioptric Image”, IEEE Int’l Symposium on 3D Data Processing, Visualization, and Transmission (3DPVT), 2006
- (8) T. Lemaire and S. Lacroix, “Monocular-vision based SLAM using Line Segments”, IEEE Proc. Int’l Conf. Robotics and Automation (ICRA), 2007.
- (9) D. Luenberger and Y. Ye, Linear and Nonlinear Programming, Third Ed., Springer, 2008.
- (10) K. Josephson and F. Kah, “Triangulation of Points, Lines and Conics”, J Math Imaging Vis, 2008.
- (11) P. Sturm, S. Ramalingam, “A Generic Concept for Camera Calibration”, European Conf. Computer Vision (ECCV), 2011
- (12) P. Miraldo and H. Araujo, “Calibration of Smooth Camera Models”, IEEE Trans. Pattern Analisys and Machine Intelligence, 2013.
- (13) P. Miraldo and H. Araujo, “Planar Pose Estimation for General Cameras using Known 3D Lines”, IEEE/RSJ Proc. Int’l Conf. Intelligent Robots and Systems (IROS), 2014.
- (14) P. Miraldo, H. Araujo, and N. Gonçalves, “Pose Estimation for General Cameras using Lines”, IEEE Trans. Cybernetics, 2015.
- (15) F. Wu, M. Zhang, G. Wang, and Z. Hu, “Algebraic Error Based Triangulation and Metric of Lines”, PLoS ONE 10(7), e0132354, 2015.
Appendix A Algorithms
In this section we show the code used in the experiments. Considering two general vectors and that do not verify the Klein constraint, using our method, the closest orthogonal vectors are given by and such that:
[gobble=2,numbers=left,frame=lines,label=LMPC] p = a1*b1 + a2*b2 + a3*b3; q = a1*a1 + a2*a2 + a3*a3 + b1*b1 + b2*b2 + b3*b3; mu = 2*p/(q+sqrt(q*q-4*p*p)); u_= 1/(1-mu*mu); x1 = (a1-mu*b1)*u_; x2 = (a2-mu*b2)*u_; x3 = (a3-mu*b3)*u_; y1 = (b1-mu*a1)*u_; y2 = (b2-mu*a2)*u_; y3 = (b3-mu*a3)*u_;
The algorithm proposed by Bartoli and Sturm is based on the singular value decomposition and can be implemented as: [gobble=2,numbers=left,frame=lines,label=BS] [U,S,V] = svd(A,0); Z = S*V’; z11 = Z(1,1); z21 = Z(2,1); z22 = Z(2,2); z12 = Z(1,2); T = [z12, z22; z21, -z11]; [ ,St,V_] = svd(T); hv = V_(:,2); hV = [hv(1),-hv(2);hv(2),hv(1)]; R = U*hV*diag(diag(hV’*S*V’)); x = R(:,1); y = R(:,2);
Usually, the singular value decomposition requires iterative techniques. However, since in this case we are dealing with matrices, is it possible to derive a closed form solution for this decompositions. Thus, in our experiments we also implemented a closed form solution for the Bartoli and Sturm method: [gobble=2,numbers=left,frame=lines,label=BS-LSVD] s1 = sqrt((a11*a11*a11*a11 + 2*a11*a11 *a12*a12 + 2*a11*a11 *a21*a21 - . . . ) s2 = sqrt(a11*a11 /2 - (a11*a11*a11*a11 + 2*a11*a11 *a12*a12 + . . . ); s2 = sqrt(s2); v11 = -((a11*a12 + a21*a22 + a31*a32))/(a11^2 + a21^2 + a31^2 - s1*s1); v21 = 1; nv = (v11*v11 + v21*v21)^(1/2); v11 = v11/nv; v21 = v21/nv; v12 = v21; v22 = -v11; u11 = (a12*v21 + a11*v11)/s1; u12 = (a12*v22 + a11*v12)/s2; u21 = (a21*v11 + a22*v21)/s1; u22 = (a21*v12 + a22*v22)/s2; u31 = (a31*v11 + a32*v21)/s1; u32 = (a31*v12 + a32*v22)/s2; z11 = s1*v11; z12 = s1*v21; z21 = s2*v12; z22 = s2*v22; t11 = z12; t12 = z22; t21 = z21; t22 = -z11; st1 = t11/2 + t22/2 - (t11*t11 - 2*t11*t22 + t22*t22 + 4*t12*t21)^(1/2)/2; st2 = t11/2 + t22/2 + (t11*t11 - 2*t11*t22 + t22*t22 + 4*t12*t21)^(1/2)/2; if st1 ¡ st2 v1 = (t12*1)/(st1 - t11); v2 = 1; nv = (v1*v1 + v2*v2)^(1/2); v1 = v1/nv; v2 = v2/nv; else v1 = (t12*1)/(st2 - t11); v2 = 1; nv = (v1*v1 + v2*v2)^(1/2); v1 = v1/nv; v2 = v2/nv; end h11 = v1; h12 = -v2; h21 = v2; h22 = v1; x1 = (u11*h11 + u12*h21)*(h11*s1*v11 + h21*s2*v12); y1 = (u11*h12 + u12*h22)*(h12*s1*v21 + h22*s2*v22); x2 = (u21*h11 + u22*h21)*(h11*s1*v11 + h21*s2*v12); y2 = (u21*h12 + u22*h22)*(h12*s1*v21 + h22*s2*v22); x3 = (u31*h11 + u32*h21)*(h11*s1*v11 + h21*s2*v12); y3 = (u31*h12 + u32*h22)*(h12*s1*v21 + h22*s2*v22);
As it can be easily seen in the previous codes, our method LMPC requires much less operations when compared with the closed form of Bartoli and Sturm algorithm BS-LSVD. On the other hand, it does not require iterative techniques, such as the BS method, which usually makes the algorithm slower.