Triangulation

COMPUTER VISION AND IMAGE UNDERSTANDING Vol. 68, No. 2, November, pp. 146–157, 1997 ARTICLE NO. IV970547 Triangulation Richard I. Hartley GE-CRD, Room K1-5C39, P.O. Box 8, Schenectady, New York 12301 and Peter Sturm GRAVIR-IMAG & INRIA Rhoˆne-Alpes, 655 Avenue de l’Europe, 38330 Montbonnot, France Received April 19, 1995; accepted August 24, 1996 point of the common perpendicular to the two rays (the midpoint method). Perhaps a better choice would be toIn this paper, we consider the problem of finding the position of a point in space given its position in two images taken with divide the common perpendicular in proportion to the cameras with known calibration and pose. This process requires distance from the two camera centers, since this would the intersection of two known rays in space and is commonly more closely equalize the angular error. Nevertheless, this known as triangulation. In the absence of noise, this problem method will not give optimal results, because of various is trivial. When noise is present, the two rays will not generally approximations (for instance, the angles will not be pre- meet, in which case it is necessary to find the best point of cisely equal in the two cases). In the case of projective intersection. This problem is especially critical in affine and reconstruction, or affine reconstruction, however, the cam-projective reconstruction in which there is no meaningful metric era matrices will be known in a projective frame of refer-information about the object space. It is desirable to find a ence, in which concepts such as common perpendicular ortriangulation method that is invariant to projective transforma- midpoint (in the projective case) have no sense. In thistions of space. This paper solves that problem by assuming a case, the simple midpoint method here will not work.Gaussian noise model for perturbation of the image coordinates. The importance of a good method for triangulation isThe triangulation problem may then be formulated as a least- squares minimization problem. In this paper a noniterative clearly shown by Beardsley et al. who demonstrate that solution is given that finds the global minimum. It is shown the midpoint method gives bad results. In [2, 3] they suggest that in certain configurations, local minima occur, which are an alternative method based on ‘‘quasi-Euclidean’’ recon- avoided by the new method. Extensive comparisons of the new struction. In this method, an approximation to the correct method with several other methods show that it consistently Euclidean frame is selected and the midpoint method is gives superior results.  1997 Academic Press carried out in this frame. The disadvantage of this method is that an approximate calibration of the camera is needed. It is also clearly suboptimal.1. THE TRIANGULATION PROBLEM This paper is an extended version of [9] which describes a new algorithm that gives an optimal global solution toWe suppose that a point x in R3 is visible in two images. the triangulation problem, equally valid in both the affineThe two camera matrices P and P9 corresponding to the and the projective reconstruction cases. The solution reliestwo images are supposed known. Let u and u9 be projec- on the concepts of epipolar correspondence and the funda-tions of the point x in the two images. From these data, mental matrix [4]. The algorithm is noniterative and simplethe two rays in space corresponding to the two image points in concept, relying on techniques of elementary calculusmay easily be computed. The triangulation problem is to to minimimze the chosen cost function. It is also moderatefind the intersection of the two lines in space. At first sight in computation requirements. In a series of experiments,this is a trivial problem, since intersecting two lines in space the algorithm is extensively tested against many otherdoes not present significant difficulties. Unfortunately, in methods of triangulation and found to give consistent supe-the presence of noise these rays cannot be guaranteed to rior performance. No knowledge of camera calibration iscross, and we need to find the best solution under some needed.assumed noise model. A commonly suggested method [2] is to choose the mid- The triangulation problem is a small cog in the machin- 146 1077-3142/97 $25.00 Copyright  1997 by Academic Press All rights of reproduction in any form reserved. TRIANGULATION 147 ery of computer vision, but in many applications of scene 3. THE MINIMIZATION CRITERION reconstruction it is a critical one, on which ultimate accu- We assume that the camera matrices, and hence theracy depends [2]. fundamental matrix, are known exactly, or at least with great accuracy compared with a pair of matching points2. TRANSFORMATIONAL INVARIANCE in the two images. A formula is given in [6] for computing the fundamental matrix given a pair of camera matrices.In the past few years, there has been considerable inter- The two rays corresponding to a matching pair of pointsest in the subject of affine or projective reconstruction [4, u ←→ u9 will meet in space if and only if the points satisfy8, 10, 12, 13, 15, 16]. In such reconstruction methods, a 3D the familiar [11] relationshipscene is to be reconstructed up to an unknown transforma- tion from the given class. Normally, in such a situation, u9ÁFu 5 0. (1)instead of knowing the correct pair of camera matrices P and P9, one has a pair PH21 and P9H21 where H is an It is clear, particularly for projective reconstruction, thatunknown transformation of the considered class. it is inappropriate to minimize errors in the 3D projectiveFor instance, in the method of projective reconstruction space, P 3. For instance, the method that finds the midpointgiven in [8] one starts with a set of image point correspon- of the common perpendicular to the two rays in space isdences ui ←→ u9i . From these correspondences, one can com- not suitable for projective reconstruction, since conceptspute the fundamental matrix F, and hence a pair of camera such as distance and perpendicularity are not valid in thematrices Pˆ and Pˆ9. In the method of [8], the pair of camera context of projective geometry. In fact, in projective recon-matrices differ from the true ones by an unknown transfor- struction, this method will give different results dependingmation H, and Pˆ is normalized so that Pˆ 5 (I u 0). Finally, on which particular projective reconstruction is consid-the 3D space points can be computed by triangulation. If ered—the method is not projective invariant.desired, the true Euclidean reconstruction of the scene Normally, errors occur not in placement of a feature inmay then be accomplished by the use of ground control space, but in its location in the two images, due to digitiza-points to determine the unknown transformation, H, and tion errors, or the exact identification of a feature in thehence the true camera matrices, P and P9. Similarly, in [7] image. It is common to assume that features in the imagesone of the steps of a projective reconstruction algorithm are subject to Gaussian noise which displaces the featureis the reconstruction of points from three views, normalized from its correct location in the image. We assume thatso that the first camera matrix has the form (I u 0). Given noise model in this paper.three or more views, an initial projective reconstruction A typical observation consists of a noisy point correspon-may be transformed to a Euclidean reconstruction under dence u ←→ u9 which does not in general satisfy the epipolarthe assumption that the images are taken all with the same constraint (1). In reality, the correct values of the corre-camera [5]. sponding image points should be points uˆ ←→ uˆ9 lying closeA desirable feature of the method of triangulation used to the measured points u ←→ u9 and satisfying the equationis that it should be invariant under transformations of the uˆ9ÁFuˆ exactly. We seek the points uˆ and uˆ9 that minimizeappropriate class. Thus, denote by t a triangulation method the functionused to compute a 3D space point x from a point correspon- dence u ←→ u9 and a pair of camera matrices P and P9. d(u, uˆ)2 1 d(u9, uˆ9)2, (2) We write where d(p, p) represents Euclidean distance, subject to the x 5 t(u, u9, P, P9). epipolar constraint uˆ9ÁFuˆ 5 0.The triangulation is said to be invariant under a transfor- mation H if Assuming a Gaussian error distribution, the points uˆ9 and uˆ are the most likely values for true image point correspon- t(u, u9, P, P9) 5 H21t(u, u9, PH21, P9H21). dences. Once uˆ9 and uˆ are found, the point x may be found by any triangulation method, since the corresponding rays This means that triangulation using the transformed cam- will meet precisely in space. eras results in the transformed point. If the camera matrices are known only up to an affine (or projective) transforma- 4. AN OPTIMAL METHOD OF TRIANGULATION tion, then it is clearly desirable to use an affine- (resp. projective)-invariant triangulation method to compute the In this section, we describe a method of triangulation that finds the global minimum of the cost function (2) using3D space points. 148 HARTLEY AND STURM a noniterative algorithm. If the Gaussian noise model can 4.2. Details of Minimization be assumed to be correct, this triangulation method is then If both of the image points correspond with the epipoles, provably optimal. This new method will be referred to as then the point in space lies on the line joining the camera the polynomial method, since it requires the solution of a centers. In this case it is impossible to determine the posi- sixth-order polynomial. tion of the point in space. If only one of the corresponding point lies at an epipole, then we conclude that the point 4.1. Reformulation of the Minimization Problem in space must coincide with the other camera center. Con- sequently, we assume in the following that neither of theGiven a measured correspondence u ←→ u9, we seek a two image points u and u9 corresponds with an epipole.pair of points uˆ and uˆ9 that minimize the sum of squared In this case, we may simplify the analysis by applying adistances (2) subject to the epipolar constraint uˆ9ÁFuˆ 5 0. rigid transformation to each image in order to place bothAny pair of points satisfying the epipolar constraint must points u and u9 at the origin, (0, 0, 1)Á in homogeneouslie on a pair of corresponding epipolar lines in the two coordinates. Furthermore, the epipoles may be placed onimages. Thus, in particular, the optimum point uˆ lies on the x-axis at points (1, 0, f )Á and (1, 0, f 9)Á, respectively.an epipolar line lˆ and uˆ9 lies on the corresponding epipolar A value f equal to 0 means that the epipole is at infinity.line lˆ9. The details on how to determine these rigid transforma-Now, we consider a pair of corresponding epipolar lines tions are given in Section 4.3. In the following, we assumel and l9. Of all pairs of points on l and l9 it is, of course, that in homogeneous coordinates, u 5 u9 5 (0, 0, 1)Á andthe pair of orthogonal projections of u on l, respectively, that the two epipoles are at points (1, 0, f )Á and (1, 0, f 9)Á.u9 on l9 which minimizes the sum of squared distances Applying these rigid transformations has no effect on(2). Let (u, u9) be the pair of these orthogonal projections. the sum-of-squares distance function (2) and hence doesWe may write d(u, u) 5 d(u, l), where d(u, l) represents not change the minimization problem. However, the funda-the perpendicular distance from the point u to the line l. mental matrix must be adapted according to these transfor-A similar expression holds for d(u9, u9). mations. Since F(1, 0, f )Á 5 (1, 0, f 9)F 5 0, the fundamentalIn view of the previous paragraph, we may reformulate matrix has a special form (how to compute this matrix fromthe minimization problem as follows. We seek to minimize the original fundamental matrix is described in Section 4.3): d(u, l)2 1 d(u9, l9)2, (3) F 5 1 ff 9d 2f 9c 2f 9d 2fb a b 2fd c d 2 . (4)where l and l9 range over all choices of corresponding epipolar lines. Suppose we have determined the pair of corresponding Consider an epipolar line in the first image passing epipolar lines lˆ and lˆ9 which minimize (3). The searched through the point (0, t, 1)Á (still in homogeneous coordi- points uˆ and uˆ9 are then just the orthogonal projections of nates) and the epipole (1, 0, f )Á. We denote this epipolar u on lˆ, respectively, u9 on lˆ9. line by l(t). The vector representing this line is given by Our strategy for minimizing (3) is as follows: the cross product (0, t, 1)T 3 (1, 0, f )T 5 (tf, 1, 2t)T, so the squared distance from the line to the origin is1. Parameterize the pencil of epipolar lines in the first image by a parameter t. Thus an epipolar line in the first d(u, l(t))2 5 t2 1 1 (tf )2 .image may be written as l(t). 2. Using the fundamental matrix F, compute the corre- sponding epipolar line l9(t) in the second image. Using the fundamental matrix to find the corresponding epipolar line in the other image, we see that3. Express the distance function d(u, l(t))2 1 d(u9, l9(t))2 explicitly as a function of t. l9(t) 5 F(0, t, 1)Á 5 (2f 9(ct 1 d), at 1 b, ct 1 d)Á. 4. Find the value of t that minimizes this function. This is the representation of the line l9(t) as a homoge- In this way, the problem is reduced to that of finding neous vector. The squared distance of this line from the the minimum of a function of a single variable, t. It will origin is equal to be seen that for a suitable parameterization of the pencil of epipolar lines the distance function is a rational polynomial d(u9, l9(t))2 5 (ct 1 d)2 (at 1 b)2 1 f 92(ct 1 d)2 .function of t. Using techniques of elementary calculus, the minimization problem reduces to finding the real roots of a polynomial of degree 6. The total squared distance is therefore given by TRIANGULATION 149 s(t) 5 t2 1 1 (tf )2 1 (ct 1 d)2 (at 1 b)2 1 f 92(ct 1 d)2 . (5) Our task is to find the minimum of this function. We may find the minimum using techniques of elemen- tary calculus, as follows. We compute the derivative s9(t) 5 2t2 (1 1 (tf )2)2 2 2(ad 2 bc)(at 1 b)(ct 1 d) ((at 1 b)2 1 f 92(ct 1 d)2)2 . (6) FIG. 1. Example of a cost function with three minima. Maxima and minima of s(t) will occur when s9(t) 5 0. Collecting the two terms in s9(t) over a common denomina- tor and equating the numerator to 0 gives a condition The (by L translated) epipole e 5 (e1 , e2 , e3)T is rotated on the x-axis, ifr(t) 5 t((at 1 b)2 1 f 92(ct 1 d)2)2 2 (ad 2 bc)(1 1 (tf )2)2(at 1 b)(ct 1 d) (7) RLe P (1, 0, f )T 5 0. for some f. Developing the left-hand side, we obtain the The minima and maxima of s(t) will occur at the roots of equation this polynomial. This is a polynomial of degree 6, which may have up to six real roots, corresponding to three min- sin Q(e1 2 e3u1) 1 cos Q(e2 2 e3u2) 5 0, ima and three maxima of the function s(t). The absolute minimum of the function s(t) may be found by finding the which allows us to determine the rotation angle Q. The roots of r(t) and evaluating the function s(t) given by (5) complete rigid transformation in the first image is given at each of the real roots. More simply, one checks the by T 5 RL. A transformation T9 for the second image is value of s(t) at the real part of each root (complex or real) determined analogously. of r(t), which saves the trouble of determining if a root is The fundamental matrix for the transformed images (the real or complex. One should also check the asymptotic same as in (4)) is then given by value of s(t) as t R y to see if the minimum distance occurs when t 5 y, corresponding to an epipolar line u 5 1/ f in F 5 T9F0T21, the first image. where F0 here denotes the fundamental matrix before car-4.3. Determining the Rigid Transformations rying out the transformations T and T9. We first carry out a translation which takes the point u 4.4. Local Minimato the origin. If u is given by u 5 (u1 , u2 , 1)T, the translation is represented by The fact that r(t) in (7) has degree 6 means that s(t) may have as many as three minima. In fact, this is indeed possible, as the following case shows. Setting f 5 f 9 5 1 and a 5 2, b 5 3, c 5 3, d 5 4 givesL 5 1 1 0 2u1 0 1 2u2 0 0 1 2. F 5 1 4 23 24 23 2 3 24 3 4 2Now, in order to place the epipole e on the x-axis, we rotate around the origin by an angle Q which is determined as shown in the following. A rotation around the origin can be represented by a matrix and a function s(t) 5 t2 1 1 t2 1 (3t 1 4)2 (2t 1 3)2 1 (3t 1 4)2R 5 1 cos Q 2sin Q 0 sin Q cos Q 0 0 0 1 2. with graph as shown in Fig. 1 (in this graph and also Fig. 150 HARTLEY AND STURM reconstruction; here we only use it to obtain the parameter- ization in Section 4.2. Of course, in practice the projection matrices or fundamental matrix are not exactly known. Correcting these improves the accuracy of the reconstruc- tion, but this requires iterative methods (cf. Section 5.5) and usually a good initialization. 5. OTHER TRIANGULATION METHODS In this section, we discuss several other triangulation FIG. 2. This is the cost function for a perfect point match, which methods that will be compared with the Polynomial nevertheless has two minima. method. 5.1. Linear Triangulation 2 we make the substitution t 5 tan(u) and plot for u in the The linear triangulation method is the most common range 2f/2 # u # f/2, so as to show the whole infinite one, described, for instance, in [8]. Suppose that u 5 Px. range of t). The three minima are clearly shown. We write in homogeneous coordinates u 5 w(u, v, 1)Á, As a second example, we consider the case where f 5 where (u, v) are the observed point coordinates and w is f 9 5 1, and a 5 2, b 5 21, c 5 1, d 5 0, i.e., an unknown scale factor. Now, denoting by pÁi the ith row of the matrix P, the equation u 5 Px may be written as wu 5 pÁ1 x, wv 5 p Á 2 x, w 5 p Á 3 x.F 5 1 0 21 0 1 2 21 0 1 0 2. Eliminating w using the third equation, we arrive at In this case, the function s(t) is given by upÁ3 x 5 p Á 1 x (8) vpÁ3 x 5 p Á 2 x. s(t) 5 t2 t2 1 1 1 t2 t2 1 (2t 2 1)2 . From two views, we obtain a total of four linear equations in the coordinates of x, which may be written in the formBoth terms of the cost function vanish for a value of t 5 Ax 5 0 for a suitable 4 3 4 matrix, A. These equations0, which means that the corresponding points u and u9 define x only up to an indeterminant scale factor, and weexactly satisfy the epipolar constraint. This can be verified seek a nonzero solution for x. Of course, with noisy data,by observing that u9ÁFu 5 0. Thus the two points are the equations will not be satisfied precisely, and we seekexactly matched. A graph of the cost function s(t) is shown a best solution.in Fig. 2. One sees apart from the absolute minimum at t 5 0 also a local minimum at t 5 1 (u 5 P/4). Thus, even The Linear–Eigen Method. There are many ways to in the case of perfect matches local minima may occur. solve for x to satisfy Ax 5 0. In one popular method, one This example shows that an algorithm that attempts to finds x to minimize iAxi subject to the condition ixi 5 1. minimize the cost function (2) or equivalently (3) by an The solution is the unit eigenvector corresponding to the iterative search beginning from an arbitrary initial point smallest eigenvalue of the matrix AÁA. This problem may is in danger of finding a local minimum, even in the case be solved using the singular value decomposition or Ja- of perfect point matches. cobi’s method for finding eigenvalues of symmetric matri- ces [1, 14]. 4.5. Optimality The Linear–LS Method. By setting x 5 (x, y, z, 1)Á Under the assumption of an unbiased Gaussian noise one reduces the set of homogeneous equations, Ax 5 0, model, the most probable reconstruction is the one that to a set of four nonhomoge