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
本文档为【Triangulation】,请使用软件OFFICE或WPS软件打开。作品中的文字与图均可以修改和编辑,
图片更改请在作品中右键图片并更换,文字修改请直接点击文字进行修改,也可以新增和删除文档中的内容。
该文档来自用户分享,如有侵权行为请发邮件ishare@vip.sina.com联系网站客服,我们会及时删除。
[版权声明] 本站所有资料为用户分享产生,若发现您的权利被侵害,请联系客服邮件isharekefu@iask.cn,我们尽快处理。
本作品所展示的图片、画像、字体、音乐的版权可能需版权方额外授权,请谨慎使用。
网站提供的党政主题相关内容(国旗、国徽、党徽..)目的在于配合国家政策宣传,仅限个人学习分享使用,禁止用于任何广告和商用目的。