Cayley Transformation and Numerical Stability of ... - Semantic Scholar
Int J Comput Vis DOI 10.1007/s11263-008-0193-x
Cayley Transformation and Numerical Stability of Calibration Equation F.C. Wu · Z.H. Wang · Z.Y. Hu
Received: 25 May 2008 / Accepted: 22 October 2008 © US Government 2008
Abstract The application of Cayley transformation to enhance the numerical stability of camera calibration is investigated. First, a new calibration equation, called the standard calibration equation, is introduced using the Cayley transformation and its analytical solution is obtained. The standard calibration equation is equivalent to the classical calibration equation, but it exhibits remarkable better numerical stability. Second, a one-parameter calibration family, called the Cayley calibration family which is equivalent to the standard calibration equation, is obtained using also the Cayley transformation and it is found that this family is composed of those infinite homographies whose rotation has the same axis with the rotation between the two given views. The condition number of equations in the Cayley calibration family varies with the parameter value, and an algorithm to determine the best parameter value is provided. Third, the generalized Cayley calibration families equivalent to the standard calibration equation are also introduced via generalized Cayley transformations. An example of the generalized Cayley transformations is illustrated, called the S-Cayley calibration family. As in the Cayley calibration family, the numerical stability of equations in a generalized Cayley calibration family also depends on the parameter value. In addition, a more generic calibration family is also proposed and F.C. Wu () · Z.H. Wang · Z.Y. Hu National Laboratory of Pattern Recognition, Institute of Automation, Chinese Academy of Sciences, PO Box 2728, Beijing 100190, People’s Republic of China e-mail: [email protected] Z.H. Wang e-mail: [email protected] Z.Y. Hu e-mail: [email protected]
it is proved that the standard calibration equation, the Cayley calibration family and the S-Cayley calibration family are all some special cases of this generic calibration family. Keywords Camera calibration · The absolute conic · Calibration equation · Cayley transformation · Numerical stability
1 Introduction Camera calibration is a necessary step to recover 3D metric information from 2D images. Many methods and techniques for camera calibration have been proposed in the last decades. The existing methods can be roughly classified into two categories: methods based on a reference calibration object and self-calibration methods. In the first kind, there are mainly the methods based on 3D object (Brown 1971; Faig 1975; Tsai 1986; Faugeras 1993), 2D pattern (Zhang 1999, 2000; Sturm and Maybank 1999) and 1D object (Zhang 2004; Wu et al. 2005; Wang et al. 2007). In these methods, camera calibration is performed by observing a 3D calibration object, a 2D calibration pattern or a 1D calibration object. Since the Euclidean geometries of calibration objects are used, the methods can provide a calibration of high-precision. In the second kind, the methods widely used are the calibration based on the Kruppa equation (Faugeras et al. 1992; Maybank and Faugeras 1992; Hartley 1997a; Luong and Faugeras 1997), the absolute conic (Pollefeys et al. 1996; Hartley 1997b; Pollefeys and Gool 1999) and the dual absolute quadric (Triggs 1997; Ponce et al. 2000). These methods do not need any reference calibration objects and use only the correspondences of image entities, such as points, lines. It is remarkable and desirable that the self-calibration methods require no reference
Int J Comput Vis
objects. Regrettably, in addition to their low calibration precision, the numerical stabilities of self-calibration methods are usually poor. Hartley and Zisserman (2000) pointed out that the algorithms absolutely depended on self-calibration should be restrainedly used in practice. Though the people have paid great attention to improve the precision and the numerical stability of self-calibration methods, no satisfactory results are obtained up to now. In this paper, we mainly investigate the numerical stabilities of the self-calibration methods based on the absolute conic. The calibration equation induced by the abT ωH = ω or, equivalently H T ω = ωH −1 , solute conic, H∞ ∞ ∞ ∞ is called the classical calibration equation in this paper. Here H∞ is the infinite homography and ω the image of the absolute conic. The previous studies are mainly focused on the conditions and computational means to solve the infinite homography (such as, Pollefeys et al. 1996; Pollefeys and Gool 1999; Hartley 1997b, etc.). For a long time the numerical stability of the calibration equation itself has been neglected in research community. In computation term, generally speaking different forms of the same equation can usually have different numerical stabilities. Thus we have the following problem: For the self-calibration based on the absolute conic, do there exist more stable forms than the classical calibration equation in computation? This problem is answered affirmatively in this paper. More specifically, we use the Cayley transformation to enhance the numerical stability of camera calibration and here are our main contributions: • Using the Cayley transformation of the infinite homography, a new calibration equation, called the standard calibration equation, is introduced and its analytical solution is obtained. The Cayley transformation of the infinite homography transforms image point to its conjugate point on the image of absolute conic. The standard calibration equation is equivalent to the classical calibration equation, but it exhibits remarkable better numerical stability. • A one-parameter family equivalent to the standard calibration equation, called the Cayley calibration family, is introduced using also the Cayley transformation and it is shown that the Cayley calibration family is composed of those infinite homographies whose rotation has the same axis with the rotation between the two given views. The condition number of equations in the Cayley calibration family varies with the parameter value, and an algorithm is also provided to determine the best parameter which gives the best numerical stability among the Cayley family. • Generalized Cayley calibration family equivalent to the standard calibration equation is introduced via generalized Cayley transformations. Given a generalized Cayley transformation, a generalized Cayley calibration family can be obtained. Similarly to the Cayley calibration
family, the numerical stability of equations in this general family varies also with the parameter value. An example of the generalized Cayley transformations is illustrated, called the S-Cayley calibration family. • A more generic calibration family is proposed and it is proved that the standard calibration equation, the Cayley calibration family and the S-Cayley calibration family are all some special cases of this generic calibration family. This paper is organized as follows. Section 2 reviews briefly the classical calibration algorithm based on the absolute conic. In the Sect. 3, the Cayley transformation is introduced. The standard calibration equation, the Cayley calibration family, and the generalized Cayley calibration family are elaborated in Sects. 4, 5 and 6, respectively. Section 7 gives more generic equivalent equations. Section 8 reports the experiment results. Section 9 concludes this paper.
2 Classical Calibration Equation The infinite homography between the two views is denoted as H∞ in this paper. If the motion between the two views is (R, t), R is the rotation and t the translation, and the intrinsic matrix of camera is K, then the infinite homography can be expressed as H∞ = KRK −1 .
(1)
Let ω be the image of the absolute conic (IAC) i.e., ω = K −T K −1 . Then there is the classical calibration equation (CCE): T ωH∞ = ω. H∞
(2)
Evidently, this equation can be changed to the following equivalent form: T −1 ω = ωH∞ . H∞
(3)
For the convenience of statement, (2) is called CCE of the first type (CCE-I), and (3) is called CCE of the second type (CCE-II). Hartley (1997b) used the CCE-II to compute ω and pointed out that the CCE can only provide 4 independent constraints on ω. Thus at least 3 views are necessary to linearly solve ω. After ω is determined, the intrinsic matrix K can be determined by the Cholesky factorization of ω, then the rotation R can in turn be determined. The classical calibration algorithm can be sketched as: Classical Algorithm Goal Given N views of a scene and assume the infinite homographies Hi∞ between the 0th and ith views are known. Compute the intrinsic matrix K and the rotation Ri between the 0th and ith views.
Int J Comput Vis
Algorithm The algorithm consists of the following three steps: 1. Let H˜ i∞ = Det−1/3 (Hi∞ )Hi∞ . Compute the least squares solution of the following set of linear equations on ω: T H˜ i∞ ωH˜ i∞ = ω
−1 T (or H˜ i∞ ω = ωH˜ i∞ ),
i = 1, 2, . . . , N − 1. 2. Determine the intrinsic matrix K by taking the Cholesky factorization of ω. 3. Determine the rotation Ri by finding the best rotaK = K −1 H˜ tion approximation of H˜ i∞ i∞ K, i.e., compute firstly the singular value decomposition of the K , H˜ K = U diag(t , t , t )V T , then set up matrix H˜ i∞ 1 2 3 i∞ T Ri = U V . Remark 1 Since the infinite homography computed in practice is equal to matrix KRK −1 up to a scale factor, it is necessary that the determinant of the computed infinite homography be normalized to 1 in the step 1, otherwise (2) and (3) do not hold. Computing Infinite Homography Here, we sketch three different methods for the infinite homography computation. Method 1 Computation with infinite correspondences (i.e., image point correspondences from infinite plane). Given two views of a scene and a set of infinite correspondences: (1) m(0) j ↔ mj ,
j = 1, 2, . . . , M.
We can linearly compute the infinite homography H∞ from the following set of linear equations: (1)
Here e is the epipole on the 2nd view, F T e = 0, a is an unknown 3-vector. Given an infinite correspondence m(0) ↔ m(1) , we have: [m(1) ]× ([e ]F + e a T )m(0) = 0. In the above three equations on a, there are only one independent constraint. And thus, at least 3 infinite correspondences are necessary to determine linearly the vector a. Method 3 Computation using the modulus constraint. For more than 2 views, using the modulus constraint shown in Pollefeys et al. (1996), i.e., the three eigenvalues of infinite homography have the same modulus, the infinite homographies between views can also be derived without resort to the infinite correspondences. But it needs to solve a set of quartic equations in three variables, and thus is highly non-linear. The computation steps are as follows. First, by the trifocal tensor method in Hartley and Zisserman (2000) or the matrix factorization methods in Sturm and Triggs (1996), Oliensis and Hartley (2007), a projective reconstruction is computed with a set of image correspondences from 3D space points: P0 = (I |0),
Pi = (Hi |ei )
(j = 1, 2, . . . , N − 1).
Then, the infinite homography between the 0th view and ith view is expressed as: Hi∞ = Hi + ei a T
(i = 1, 2, . . . , N − 1).
Here a is an unknown 3-vector, the direction of infinite plane in the reconstructed projective space. And thus, the vector a can be non-linearly computed by using the modulus constraints.
(0)
[mj ]× H∞ mj = 0 (j = 1, 2, . . . , M). This method needs at least 4 infinite correspondences because each infinite correspondence can only give rise to two independent constraints on H∞ . Hartley (1997b) pointed out that if the motion between two views is a rotation around the camera center, then infinite correspondences in the above can be replaced by image correspondences from 3D space points. Method 2 Computation with the epipolar geometry. Given the foundational matrix between two views, F , which can be computed with image correspondences from 3D space points, three infinite correspondences are needed in order to linearly compute the infinite homography H∞ . This is because the infinite homography H∞ can be expressed as: H∞ = [e ]F + e a T .
3 Cayley Transformation The following function of complex variable is called the Cayley transformation (Bell 1992): ξ = ϕ(λ) =
1−λ . 1+λ
(4)
It is a conformal mapping on the extended complex plane. It is ready to verify that λ = ϕ −1 (ξ ) =
1−ξ = ϕ(ξ ). 1+ξ
So, the inverse transformation is also the Cayley transformation. For a matrix A ∈ C n×n , let σ (A) denote the set of all eigenvalues of A. Assume −1 ∈ / σ (A) i.e., det(I + A) = 0,
Int J Comput Vis
then the Cayley transformation of A is defined by
By right-multiplying the two sides of this equation by C, and notice that B T C + CB = 0, have
B = ϕ(A) = (I − A)(I + A)−1 −1
= (I + A)
(I − A).
(5)
AT C(I − B) = C(I + B). So,
From this equation, we have (I + A)B = (I − A) or, equivalently, A(I + B) = (I − B).
(6)
If Bx = −x, then by (5) there must be x = 0, and thus, −1 ∈ / σ (B). Therefore, from (6) we have
= (I + B)−1 (I − B).
(7)
This indicates that the matrix A is also the Cayley transformation of B. According to the above discussions, the Cayley transformation of matrices is a one-to-one mapping on the matrix set {A ∈ C n×n | − 1 ∈ / σ (A)} and we can prove that: ⇔
β=
= C(I + B)(I + B)−1 = C. Hence, the necessary condition holds too.
4 Standard Calibration Equation
A = ϕ(B) = (I − B)(I + B)−1
λ ∈ σ (A)
AT CA = AT C(I − B)(I + B)−1
1−λ ∈ σ (ϕ(A)). 1+λ
(8)
For the Cayley transformation of matrices, we have the following proposition, which is the fundamental basis of our current work. Proposition 1 Assume matrices A, B ∈ satisfy (5) and (7) and let C ∈ C n×n be an arbitrary matrix, then, C n×n
4.1 Cayley Transformation of H∞ The infinite homography H∞ satisfies det(I + H∞ ) = det(K(I + R)K −1 ) = det(I + R). If the rotation angle of R is not equal to iπ (i ∈ Z), det(I + R) = 0. So, in this case, there is the Cayley transformation of H∞ : H = ϕ(H∞ ) = (I − H∞ )(I + H∞ )−1 = (I + H∞ )−1 (I − H∞ ).
(10)
And H∞ is also the Cayley transformation of H : H∞ = ϕ(H ) = (I − H )(I + H )−1 = (I + H )−1 (I − H ).
(11)
(9)
In this paper, we always assume that the rotation angle θ of R is not equal to iπ (i ∈ Z).
Proof By left-multiplying the two sides of (5) by (I + A), and transposed, we obtain
Proposition 2 H is a conjugate of an anti-symmetric matrix.
B T (I + A)T = (I − A)T .
Proof Substituting H∞ = KRK −1 into (10) yields
AT CA = C
⇔
B T C + CB = 0.
By right-multiplying the two sides of this equation by CA, and notice that AT CA = C, we have B T C(I + A) = −C(I − A).
H = (I − KRK −1 )(I + KRK −1 )−1 = K (I − R)(I + R)−1 K −1 = Kϕ(R)K −1 .
(12)
ϕ(R)
So,
Where ϕ(R) is the Cayley transformation of rotation R. Since
B T C = −C(I − A)(I + A)−1 = −CB.
(ϕ(R))T = (I + R T )−1 (I − R T )
Therefore the sufficient condition holds. Now, we prove the necessary condition next. By left-multiplying the two sides of (7) by (I + B), and transposed, we obtain AT (I + B)T = (I − B)T .
= (I + R −1 )−1 (I − R −1 ) = (R + I )−1 RR −1 (R − I ) = (R + I )−1 (R − I ) = −ϕ(R),
Int J Comput Vis
ϕ(R) is an anti-symmetric matrix, and thus by (12) the proposition holds. {1, eiθ , e−iθ }.
Let the eigenvalues of rotation R be Then it is not difficult to see that ϕ(R) has the eigenvalues {0, i tan(θ/2), −i tan(θ/2)} and has the same eigenvectors with R, i.e., ⎧ ⎧ ⎨ϕ(R)x = 0 ⎨ Rx = x ϕ(R)y = i tan(θ/2)y ⇔ Ry = eiθ y ⎩ ⎩ ϕ(R)y = i tan(θ/2)y ϕ(R)y = e−iθ y. According to Proposition 2, H has also the eigenvalues {0, i tan(θ/2), −i tan(θ/2)} and the corresponding eigenvectors are {Kx, Ky, Ky}. Corollary 1 (a) rank(H ) = 2. (b) Let H K = K −1 H K. Then the rotation R is the Cayley transformation of H K : R = ϕ(H K ) = (I − H K )(I + H K )−1 .
(13)
Proposition 3 Let the eigenvalues of H be {0, is, −is}, then, (a) H can be factorized as H = P [a]× P −1 ,
(14)
where P is a real matrix, [a]× is the anti-symmetric matrix defined by the vector a = (s, 0, 0)T , i.e. ⎛ ⎞ 0 0 0 [a]× = ⎝ 0 0 s ⎠ . 0 −s 0 Equation (14) is called the standard anti-symmetric factorization of H in this paper. (b) If H has another standard anti-symmetric factorization H = P [a]× P −1 , then P must be of the form P = P W (α, β, γ ), where ⎛
α ⎝ W= 0 0
0 β −γ
(15) ⎞
0 γ⎠ β
is an invertible matrix. ¯ be the unit eigenvectors of H correProof Let {p1 , q, q} sponding to the eigenvalues {0, is, −is}, i.e. ⎧ ⎨ H p1 = 0 H q = isq ⎩ ¯ H q¯ = −is q.
Here q¯ is the conjugate of q. Then, we obtain ¯ = (p 1 , q, q)diag(0, ¯ is, −is). H (p 1 , q, q) By right-multiplying the two sides of this equation by ⎞ ⎛ 1 0 0 √ √ ⎟ ⎜ U = ⎝ 0 1/ 2 1/i 2 ⎠ , √ √ 0 1/ 2 −/i 2 since
q + q¯ q − q¯ ¯ (p 1 , q, q)U ∈ R 3×3 = p1 , √ , √ 2 i 2
and
¯ ¯ [a]× is, −is)U = (p 1 , q, q)U (p 1 , q, q)diag(0, q + q¯ q − q¯ = p1 , √ , √ [a]× , 2 i 2 we have q + q¯ q − q¯ q + q¯ q − q¯ H p1 , √ , √ = p1 , √ , √ [a]× 2 i 2 2 i 2 P
P
and thus (14) holds. We now prove (b). Assume H = P [a]× P −1 is another standard anti-symmetric factorization, then P [a]× P −1 = P [a]× P −1 , and thus, −1 −1 P P[a]× − [a]× P P = 0. W
W
According to Lemma 1 below, solutions of W are ⎞ ⎛ x 0 0 ⎟ ⎜ W = xI + y[a]× + z[a]2× = ⎝ 0 x − zs 2 ys ⎠ , 0 −ys x − zs 2 x(y 2 s 2 + (x − zs 2 )2 ) = 0. Let α = x, β = x − zs 2 , γ = ys, then, W = W (α, β, γ ) = P −1 P . Hence, P = P W .
Lemma 1 Assume c is a nonzero 3-vector, then solutions of the equation [c]× X − X[c]× = 0 are: X = xI + y[c]× + z[c]2× ,
x, y, z ∈ R.
(16)
Proof Suppose A, B ∈ R n×n and let vec(B) = (bT1 , bT2 , . . . , bTn )T , where bTi is the ith row of B, then we can obtain vec(AB) = (A ⊗ In )vec(B),
(17)
Int J Comput Vis
vec(BA) = (In ⊗ AT )vec(B).
(18)
Here ⊗ denotes the tensor product of matrices. By applying (17) and (18), [c]× X − X[c]× = 0 can be changed as ([c]× ⊗ I3 + I3 ⊗ [c]× )vec(X) = 0.
(19)
Since the three eigenvalues of anti-symmetric matrix [c]× are {0, i c 2 , −i c 2 }, the coefficient matrix of (19) has eigenvalue zero of multiplicity 3. It is easy to verify that {vec(I ), vec([c]× ), vec([c]2× )} are three linear independent solutions of (19). And thus its solutions are
(b) Symmetrical positive-definite solutions of the SCE are: ω = xP −T diag(1, w, w)P −1 ,
x > 0, w > 0.
(22)
Proof Substituting H = P [a]× P −1 into (20) and by some simple manipulations, we have −[a]× P T ωP + P T ωP [a]× = 0.
(23)
By denoting X = P T ωP in (23), we have [a]× X − X[a]× = 0.
(24)
According to Lemma 1, its solutions are vec(X) = xvec(I ) + yvec([c]× ) + zvec([c]2× ),
X = xI + z[a]× + y[a]2× ,
x, y, z ∈ R.
x, y, z ∈ R.
Since It yields that the solutions of [c]× X − X[c]× = 0 are X = xI + y[c]× + z[c]2× ,
x, y, z ∈ R.
X − XT = 0
⇔
− (xI + z[a]× + y[a]2× )T = 0
4.2 Standard Calibration Equation ⇔ By Proposition 1, the following proposition holds.
(20)
This equation is called the standard calibration equation (SCE) in this paper, which has the similar form with the CCE-II. Though the CCE-I and the SCE are equivalent to each other, their forms have remarkable differences. First, the CCE-I is of a symmetrical form, while the SCE is of an anti-symmetric form because it implies that the matrix ωH is anti-symmetric. Second, the CCE-I is homogeneous on only the matrix ω, while the SCE is homogeneous on not only the matrix ω but also the matrix H . Third, the coefficients of ω are quadratic function of the H∞ ’s elements in the CCE-I, and are linear function of the H ’s elements in the SCE. Besides these, the SCE has a nice property as stated in the following proposition. Proposition 5 Let H = P [a]× P −1 is a standard antisymmetric factorization of H (see Proposition 3), then, (a) Symmetrical solutions of the SCE are: ω = P −T diag(x, x − ys 2 , x − ys 2 )P −1 , x, y ∈ R.
2z[a]× = 0
⇔
z = 0,
symmetrical solutions of (24) are
Proposition 4 The CCE-I (or CCE-II) is equivalent to the equation below: H T ω + ωH = 0.
(xI + z[a]× + y[a]2× )
(21)
X = xI + y[a]2× ,
x, y ∈ R.
Notice that [a]2× = diag(0, −s 2 , −s 2 ), we have X = diag(x, x − ys 2 , x − ys 2 ),
x, y ∈ R.
Hence symmetrical solutions of the SCE are ω = P −T diag(x, x − ys 2 , x − ys 2 )P −1 ,
x, y ∈ R.
Evidently, ω is a symmetrical positive-definite matrix if and only if X is a symmetrical positive-definite matrix, and thus symmetrical positive-definite solutions of the SCE are ω = P −T diag(x, x − ys 2 , x − ys 2 )P −1 ,
x > 0, y < xs −2 .
Let w = 1 − x −1 ys 2 , we have ω = xP −T diag(1, w, w)P −1 ,
x > 0, w > 0.
Remark 2 The symmetrical solutions (21) and symmetrical positive-definite solutions (22) of the SCE do not depend on the specific selection of matrix P in the standard antisymmetric factorization of H . This is because: Assume H = P [a]× P −1 is another standard anti-symmetric factorization, then according to Proposition 3(b), P = P W (α, β, γ ). It is not difficult to verify that W diag(a, b, b) = diag(a, b, b)W,
∀a, b
Int J Comput Vis
or, equivalently, W diag(a, b, b)W −1 = diag(a, b, b),
Proof By Proposition 5, the symmetrical positive-definite solutions of HiT ω + ωHi = 0 are
∀a, b.
And thus, P −T diag(a, b, b)P −1 = P −T diag(a, b, b)P −1 . Hence, the symmetrical solutions (21) and symmetrical positive-definite solutions (22) do not depend on the selection of matrix P in the standard anti-symmetric factorization. Proposition 5 tells us that from the SCE of a pair of views, the IAC ω can be determined up to a one-parameter family. Thus under the 4-parameter camera model with zero skew, the IAC ω can be uniquely determined from the SCE of a pair of views and its computation is very easy with (22): q11 q12 , ω = P −T diag 1, − q21 q22 + q31 q32 q11 q12 P −1 , − q21 q22 + q31 q32 where qij is the (i, j )-element of P −1 . Given two SCEs from three views, the following proposition provides an analytical solution of ω. Proposition 6 Assume Ri is the rotation between the 0th and ith views and the corresponding SCE is: HiT ω + ωHi
= 0,
i = 1, 2.
(25)
Then (25) has a unique solution if and only if that the rotation axes of the two rotations are not parallel. And if we let the standard anti-symmetric factorization of Hi be Hi = Pi [ai ]× Pi−1 , then its unique solution is: ω = P1−T diag(1, w1 , w1 )P1−1 = x2 P2−T diag(1, w2 , w2 )P2−1 , where 2 (Q21 + Q231 ) 3j =2 (Q22j + Q23j ) , w1 = 3 2 + q2 ) 2 2 (q21 j =2 (q2j + q3j ) 31 Q211 + w12 (Q212 + Q213 ) x2 = , 2 + q2 + q2 ) det2 (Q)(q11 12 13 Q221 + Q231 w2 = , 2 + q2 ) det2 (Q)(q21 31 (qij ) = (q 1 , q 2 , q 3 )
T
= P2−1 P1 = Q,
(Qi1 , Qi2 , Qi3 )T = q j × q k , i = 1, 2, 3; j, k = i and j < k.
(26)
ω = xi Pi−T diag(1, wi , wi )Pi−1 (xi > 0, wi > 0),
i = 1, 2. (27)
Assume x1 = 1. Then from (27), we have diag(1, w1 , w1 ) = x2 QT diag(1, w2 , w2 )Q = x2 (q 1 q T1 + w2 (q 2 q T2 + q 3 q T3 )). It is easy to see that this equation is equivalent to the following equation: ⎧ ⎪ ⎨diag(1, w1 , w1 )(q 2 × q 3 ) = x2 det(Q)q 1 diag(1, w1 , w1 )(q 1 × q 3 ) = −w2 det(Q)q 2 (28) ⎪ ⎩ diag(1, w1 , w1 )(q 1 × q 2 ) = w2 det(Q)q 3 . (a) If the rotation axes of the two rotations are parallel to each other, then from the proof of Propositions 2 and 3, the first rows of P1 and P2 are equal up to a nonzero scale. And thus Q can be expressed as ⎛ ⎞ α ∗ ∗ Q = P2−1 P1 = ⎝ 0 ∗ ∗ ⎠ , α = 0. 0 ∗ ∗ Hence, (28) is equivalent to the equation below: ⎧ ⎨det(Q)x2 = α Q23 w1 = −q23 w2 Q22 w1 = −q22 w2 , ⎩ Q32 w1 = −q32 w2 , Q33 w1 = −q33 w2 . Clearly, its solution is not unique, and thus the solution of (25) is not unique either. (b) If the rotation axes of two rotations are not parallel, then the two elements q21 and q31 of Q are not all zero. Rewriting the 2nd and 3rd equations in (28) as (Q21 , Q22 w1 , Q23 w1 ) = − det(Q)(q21 w2 , q22 w2 , q23 w2 ) (Q31 , Q32 w1 , Q33 w1 ) = det(Q)(q31 w2 , q32 w2 , q33 w2 ). 2 + q 2 = 0, we have Since q21 31 Q221 + Q231 w2 = , 2 2 + q2 ) det (Q)(q21 31 2 (Q21 + Q231 ) 3j =2 (Q22j + Q23j ) . w1 = 3 2 + q2 ) 2 2 (q21 j =2 (q2j + q3j ) 31
Then, from the first equation in (28) we have Q211 + w12 (Q212 + Q213 ) x2 = . 2 + q2 + q2 ) det2 (Q)(q11 12 13 Hence, (25) has the unique solution (26).
Int J Comput Vis
It is well-known that having the parallel rotation axis is a degenerate configuration of the set of CCEs (Hartley and Zisserman 2000). However, Proposition 6 more explicitly shows this is the unique degenerate configuration of the set of CCEs (or, SCEs). Since an exact factorization, H = P [a]× P −1 , is not possible in practice due to the data error, we need a numerical algorithm to determine the IAC ω.
Second, though the rotation R can be computed directly with the infinite homography H∞ once the intrinsic matrix K is known, we recommend that the rotation R be computed with the formula R = ϕ(H K ). Since the intrinsic matrix K here is obtained from the SCE, R should be computed by R = ϕ(H K ). The combination of such an R and K can best match the image data, and thus can give a better computation of the rotation.
4.3 Standard Algorithm
4.4 Geometrical Interpretation of H
Goal Given N views of a scene and assume the infinite homographies Hi∞ between the 0th and ith views are known. Compute the intrinsic matrix K and the rotation Ri between the 0th and ith views. Algorithm The algorithm consists of the following three steps: 1. Let H˜ i∞ = Det−1/3 (Hi∞ )Hi∞ and compute its Cayley transformation Hi = ϕ(H˜ i∞ ). Then find the least squares solution of the set of linear equations on ω: HiT ω + ωHi = 0
(i = 1, 2, . . . , N − 1).
2. Compute the intrinsic matrix K by the Cholesky factorization of ω. 3. Let HiK = K −1 Hi K and compute its Cayley transformation, R˜ i = ϕ(H K ). Then, determine the rotation Ri by finding the best rotation approximation of the matrix R˜ i , i.e., firstly compute the singular value decomposition of the matrix R˜ i , R˜ i = U diag(t1 , t2 , t3 )V T , then set up Ri = U V T . This algorithm is called the standard algorithm in this paper. In the absence of data error and with the infinite precision of computation, the CCE and the SCE should give the same ω. However, in the presence of data error and with numerical computation, their numerical performances manifest a great difference (see Sect. 5.3). Remark 3 First, the matrix H has an important constraint of rank 2. In practice due to the data error, the computed H does not satisfy this constraint. If in the step 1 this constraint is enforced, both success rate and calibration precision of the algorithm can be improved. The following shows an Rank-enforced method: Firstly, compute the eigenvalues factorization of the matrix H , H = U diag(s, t, t)U −1 , then let its real eigenvalue s = 0, and replace H with the matrix = Re(U diag(0, t, t)U −1 ). H Here Re(A) denotes the real part of A. In theory U diag(0, t, t)U −1 is a real matrix, but a very small imaginary part is possible owing to the numerical computation.
We will give a geometrical interpretation of H in this subsection. Let m ↔ m be an image correspondence, both m and m being projections of an infinite point X∞ = (XT , 0)T . Then, αm = KX,
βm = KRX,
and thus, αm − βm = K(I − R)X, αm + βm = K(I + R)X. By the above two equations, we obtain αm − βm = K(I − R)(I + R)−1 K −1 (αm + βm ). Therefore, m − γ m = H (m + γ m )
(γ = β/α).
(29)
From (29), we can see that the matrix H transforms the linear combination m + γ m of an infinite correspondence m ↔ m into another linear combination m − γ m . This indicates that the matrix H is not a transformation between the two views. Equation (29) cannot be used for a linear computation of the matrix H from infinite correspondences due to the unknown scale γ . But a nonlinear solution is possible. By eliminating the unknown scale γ , (29) becomes [m + H m ]× (m − H m) = 0.
(30)
It is a nonlinear constraint on H . If we regard (I − H∞ ) and (I + H∞ ) as two transformation from the first view into the second view, then H = (I + H∞ )−1 (I − H∞ ) is a transformation on the first view, and is also a transformation on the second view because (I + H∞ )−1 is commutative with (I − H∞ ) i.e., H = (I + H∞ )−1 (I − H∞ ) = (I − H∞ )(I + H∞ )−1 .
Int J Comput Vis
(a)
(b)
Fig. 1 The intrinsic geometrical meaning of the transformation H : the vanishing point v ∞ of the rotation axis of R is a 1-dimensional zero space of H and the vanishing line l ∞ of plane orthogonal to the rotation axis of R is 2-dimensional invariant subspace of H . (a) For
an arbitrary x ∈ / l ∞ ∪ {v ∞ }, y = H x ∈ l ∞ is a conjugate point of x on the IAC and H is a transposition on {y, H y}. (b) Except for v ∞ , H transforms all points on l to its conjugate point c(l) on the IAC and H transforms c(l) to the intersection of l and l ∞
By considering H as a transformation on the same view, we can give its intrinsic geometrical meaning as stated in the following proposition.
By the proof of Proposition 2, H has the same eigenvectors with H∞ . Since l ∞ is a 2-dimensional invariant subspace of H∞ , it is also a 2-dimensional invariant subspace of H , and thus {H x, H 2 x, H 3 x} ⊂ l ∞ . By (a), both H 2 x and H 3 x are the conjugate points of H x and H 2 x = H x, therefore there must be H 3 x = H x. Hence, H is a transposition on {H x, H 2 x}. (c) Let x be a fixed point on the line l but not equal to v ∞ , then ∀y ∈ l\{v ∞ } there is
Proposition 7 Let v ∞ be the vanishing point of the rotation axis of R and l ∞ be the vanishing line of plane orthogonal with the rotation axis of R. Assume x is an arbitrary image point, then, (a) If x is not the vanishing point v ∞ , then H transforms x to its conjugate point y on the IAC, i.e., if y = H x, then (x, y) satisfies x T ωy = 0. (b) ∀x ∈ / l ∞ ∪ {v ∞ }, there must be H x ∈ l ∞ and H is a transposition on {H x, H 2 x}, i.e., H (H x) = H 2 x and H (H 2 x) = H x, as shown in Fig. 1(a). (c) Let l be an arbitrary line passing through v ∞ , and c(l) be its conjugate point on the IAC. Then, except for the vanishing point v ∞ , H transforms all points on the line l to the same point c(l) and H transforms c(l) to the intersection of the line l and the line l ∞ , q = l ∞ ∩ l, as shown in Fig. 1(b). Proof (a) From the SCE (20), the matrix ωH is antisymmetric. And thus for any image point x, there must be x T ωH x = 0. Note that m ∈ N (H ) if and only if m = H∞ m, namely, m is the vanishing point v ∞ . Therefore, H x = 0, and thus y = H x is a conjugate point of x on the IAC. (b) ∀x ∈ / l ∞ ∪ {v ∞ }, there must be x = αv ∞ + βx 1 + γ x 2 , where x 1 , x 2 ∈ l ∞ and α, β, γ = 0. Then, we have H x = αH v ∞ + βH x 1 + γ H x 2 = βH x 1 + γ H x 2 ∈ l ∞ .
y = αv ∞ + βx,
β = 0.
And thus, H y = βH x, i.e., H y and H x are the same point. By (a), H x is the conjugate point of the line l on the IAC, c(l). By (b), {c(l), H c(l)} ⊂ l ∞ and H c(l) is a conjugate point of c(l) on the IAC, and thus H c(l) must be the intersection of l and l ∞ . Remark 4 This is an interesting proposition. First, this proposition tells us that using H we can obtain the constraints on the IAC from the vanishing points of two orthogonal directions to each other. Obviously, by the matrix H we cannot obtain the constraints on the IAC from all orthogonal directions, for example, the constraints from two orthogonal directions of the rotation axis of R and its orthogonal line. But, from (b) and (c), it is not difficult to see that H can provide the four independent constraints. Second, this proposition tells also us that H 2 is a 2-dimensional projection with the center v ∞ and the image line l ∞ as shown Fig. 2(a), and H 3 is H self, i.e., H 3 = H . Observing the transformation property of H on two lines passing through v ∞ and a pair of conjugate points on l ∞ , we can have the Euclidean properties of H , as shown in Fig. 2(b). In essence, the constraints on the IAC induced by H are from the orthogonality of 3D structure, while H∞ are from the orthogonality of 3D motion R. Hence, the Cayley transformation, H = ϕ(H∞ ), is a conversion between the
Int J Comput Vis
(a)
(b)
Fig. 2 (a) H 2 is a 2-dimensional projection with the center v ∞ and the image line l ∞ . It maps point x to the intersection of line l ∞ and the line passing through v ∞ and x. (b) Let L(m) and (l) be the directions of the line and the plane in 3D space corresponding to m and l and assume p and q are a pair of conjugate points on the IAC, then
L(v ∞ )⊥ (l ∞ ) and L(p)⊥L(q). Since x is on line v ∞ ×p, L(x) is on the plane L(v ∞ ) × L(p), similarly L(y) is on the plane L(v ∞ ) × L(q). Proposition 7 shows that H transforms L(x) to its orthogonal direction L(q), L(y) to its orthogonal direction L(p), and transposes the two orthogonal directions L(p) and L(q) to each other
3D motion orthogonality and the 3D structure orthogonality. Third, this proposition gives a linear constraint on H as follows. If a line L and a plane P are orthogonal to each other in 3D space, then their vanishing point x and vanishing line l satisfy l T H x = 0. The only exception is the case of rotation axis of R and its orthogonal plane. A proof is ready: Assume the line L is not parallel to the rotation axis of R, then the relation of its vanishing point x ∈ / N (H ) and the vanishing line of its orthogonal plane is l = ωx. By (a), x T ω(H x) = 0, i.e., l T H x = 0 which is a linear constraint on H .
Equation (31) is called the Cayley transformation family of H . Next, we show some properties of the Cayley transformation family. Proposition 8 (a) ϕ(1 · H ) = H∞ ,
ϕ(0 · H ) = I.
(32)
(b) For τ ∈ (−∞, +∞), ϕ −1 (τ H ) = ϕ(−τ H )
or,
ϕ(τ H ) · ϕ(−τ H ) = I, τ H = ϕ 2 (τ H ).
(33) (34)
5 Cayley Calibration Family
Proof We only prove (33) here. Both (I − τ H ) and (I + τ H ) are invertible for τ ∈ (−∞, +∞), therefore
5.1 Cayley Transformation Family of H
ϕ −1 (τ H ) = (I + τ H )(I − τ H )−1
Since H = (I − H∞ )(I + H∞ )−1 has the eigenvalues:
= (I − (−τ H ))(I + (−τ H ))−1 = ϕ(−τ H ).
{0, i tan(θ/2), −i tan(θ/2)}, for τ ∈ (−∞, +∞) the eigenvalues of matrix I ± τ H must be {1, 1 ± iτ tan(θ/2), 1 ∓ iτ tan(θ/2)}.
The matrix ϕ(τ H ) can be computed from the infinite homography H∞ by (35) via the following proposition. Proposition 9 ϕ(τ H ) = ((1 + τ )I + (1 − τ )H∞ )−1 × ((1 − τ )I + (1 + τ )H∞ ).
Then, det(1 ± τ H ) = 1 + τ 2 tan2 (θ/2) = 0.
Proof
Hence (τ H ) has the Cayley transformation:
ϕ(τ H ) = (I + τ H )−1 (I − τ H ) = (I + τ (I + H∞ )−1 (I − H∞ ))−1
ϕ(τ H ) = (I − τ H )(I + τ H )−1 for τ ∈ (−∞, +∞).
(31)
× (I − τ (I + H∞ )−1 (I − H∞ ))
(35)
Int J Comput Vis
= ((I + H∞ ) + τ (I − H∞ ))−1 (I + H∞ )
equations are equivalent to the CCE-I. Evidently, the CCF-I has the following equivalent form:
× (I + H∞ )−1 ((I + H∞ ) − τ (I − H∞ ))
ϕ T (τ H )ω = ωϕ −1 (τ H )
= ((1 + τ )I + (1 − τ )H∞ )−1 × ((1 − τ )I + (1 + τ )H∞ ).
Let ϕK (τ H ) = K −1 ϕ(τ H )K. If the intrinsic matrix K is known, then we can compute the rotation R by (36) via the following proposition. Proposition 10 R = ((τ + 1)I + (τ − 1)ϕK (τ H ))−1 × ((τ − 1)I + (τ + 1)ϕK (τ H )).
(36)
Proof From (35), we have ϕK (τ H ) = K −1 ϕ(τ H )K
× ((1 − τ )I + (1 + τ )R) = ((1 − τ )I + (1 + τ )R) × ((1 + τ )I + (1 − τ )R)
By Proposition 8, ϕ −1 (τ H ) = ϕ(−τ H ), and thus the above equation can be rewritten as ϕ T (τ H )ω = ωϕ(−τ H )
(τ = 0).
(38)
The advantage of this form is to avoid the matrix inverse operations, hence improve numerical stability. Equation (38) is called Cayley calibration family of the second type, denoted as CCF-II(τ ). Next, we give an explanation of the CCF in term of the infinite homography. In other words, we will prove that for τ = 0, the Cayley transformation ϕ(τ H ) is an infinite homography. Proposition 12 For τ = 0, the Cayley transformation ϕ(τ H ) of τ H is an infinite homography and the corresponding rotation has the same rotation axis with the rotation between the two views.
= ((1 + τ )I + (1 − τ )R)−1
−1
(τ = 0).
Proof Let R be the rotation between the two given views. Then according to Proposition 9, we have
.
ϕ(τ H )
Hence,
= K ((1 + τ )I + (1 − τ )R)−1 ((1 − τ )I + (1 + τ )R) K −1 .
ϕK (τ H )((1 + τ )I + (1 − τ )R)
U (τ )
= ((1 − τ )I + (1 + τ )R). By some simple manipulations, we can obtain (36).
(39)
5.2 Cayley Calibration Family Since for τ ∈ (−∞, +∞), ϕ(τ H ) is the Cayley transformation of τ H , according to Proposition 1 the equation ϕ T (τ H )ωϕ(τ H ) = ω is equivalent to the equation (τ H )T ω + ω(τ H ) = 0. When τ = 0, both the two equations become the identity equation, and thus we do not consider this trivial case. Evidently, for τ = 0 the SCE is equivalent to the equation (τ H )T ω + ω(τ H ) = 0. Thus, we have the following proposition.
In order to prove Proposition 12, we need only to prove that U (τ ) is a rotation and it has the same rotation axis with R. For this purpose, let the eigenvalues factorization of R be: R = V diag(1, eiθ , e−iθ )V −1 .
(40)
Proposition 11 The SCE is equivalent to the following equation:
Then θ is the rotation angle of R, V = (v 1 , v 2 , v 3 ) is a unitary matrix and v 1 is the rotation axis of R. Substituting (40) into U (τ ) in (39) and by some simple manipulations, we can obtain 1 − τ + (1 + τ )eiθ U (τ ) = V diag 1, , 1 + τ + (1 − τ )eiθ 1 − τ + (1 + τ )e−iθ V −1 . 1 + τ + (1 − τ )e−iθ
ϕ T (τ H )ωϕ(τ H ) = ω
It is not difficult to prove that
(τ = 0).
(37)
The above equation gives a one-parameter family equivalent to the SCE. Since it has the same form with the CCE-I, we call it Cayley calibration family of the first type, denoted as CCF-I(τ ). In the CCF-I(τ ), the equation corresponding to τ = 1is just the CCE-I, i.e. CCF-I(1) = CCE-I, and all the
1 − τ + (1 + τ )eiθ 1 + τ + (1 − τ )eiθ 2 cos(θ/2) + iτ sin(θ/2) = = z(τ ), cos2 (θ/2) + τ 2 sin2 (θ/2)
Int J Comput Vis
1 − τ + (1 + τ )e−iθ 1 + τ + (1 − τ )e−iθ 2 cos(θ/2) − iτ sin(θ/2) = = z(τ ). cos2 (θ/2) + τ 2 sin2 (θ/2) Thus, U (τ ) can be rewritten as U (τ ) = V diag(1, z(τ ), z(τ ))V T .
(41)
Since z(τ )is a unit complex number, U (τ ) is a rotation and v 1 is its rotation axis. Remark 5 From the proof of Proposition 12, we can see that the rotation angle of rotation U (τ ) corresponding to infinite homography ϕ(τ H ) is θ (τ ) = 2 Arctan(τ tan(θ/2)).
(42)
Thus if τ varies in (−∞, +∞), we can obtain all the rotations with the fixed rotation axis v 1 . In other words, given an infinite homography H∞ with the rotation R (rotation axis v and rotation angle θ ), by the Cayley transformation of τ ϕ(H∞ ) or, τ H we can obtain all the infinite homographies with the fixed rotation axis v but arbitrary rotation angle, and thus Proposition 6 tells us again that the configuration of cameras with the same rotation axis is a degenerate configuration. Of course, by changing the rotation angle from the eigenvalues factorization of H∞ , we can also obtain these infinite homographies, but it needs to compute the eigensystem of H∞ . In the contrast, we here need merely taking two times of the Cayley transformation, ϕ(τ ϕ(H∞ )), in order to change the rotation angle. 5.3 Improving Calibration Conditioning Given the infinite homography H∞ between two views, we can obtain the CCF-I(τ ) equivalent to the CCE-I: CCF-I(τ ) :
ϕ T (τ H )ωϕ(τ H ) − ω = 0,
τ = 0,
and CCE-I = CCF-I(1) (ϕ(1 · H ) = H∞ ). In the CCF-I(τ ), if we define CCF-I(0) by H T ω + ωH = 0 (SCE), then by Proposition 10 we obtain an equivalent family with τ in (−∞, +∞), denoted as CCF-I also. Similarly, for the CCF-II(τ ), we define CCF-II(0) also by the SCE, and thus obtain an equivalent family of the second type with τ in (−∞, +∞). The numerical stability of a homogeneous linear system, such as Ax = 0, can be assessed by examining the condition number of its coefficient matrix A, it is also called the condition number of this homogeneous linear system. The condition number is a typical measure of the sensitivity of homogeneous linear system to noise. A large condition number
usually means less numerical stability, and a small condition number indicates good numerical stability. If the rank of the matrix A is k, then its condition number is defined as the ratio of the largest singular value to the kth smallest singular value (Hartley 1997c). And thus the condition number of the matrix A is computed by taking the singular value decomposition of this matrix. Figure 3 gives the condition numbers of equations in 5 CCF-Is and 5 CCF-IIs. From this figure, we can see that different equation in the CCF-I (or CCF-II) has different condition number, and thus has different numerical stability. Hence, we can improve the numerical stability of calibration by selecting a parameter τ such that corresponding equation has a smaller condition number. We next give an algorithm for selecting the best parameter value. Let the set of m CCF-Is (or CCF-IIs) be ϕ T (τi Hi )ωϕ(τi Hi ) − ω = 0,
i = 1, 2, . . . , m
(or, ϕ T (τi Hi )ω − ωϕ(−τi Hi ) = 0, i = 1, 2, . . . , m). The m-parameter family of the above homogeneous linear systems can be rewritten as the following matrix form: τ = (τ1 , τ2 , . . . , τm ) ∈ (−∞, +∞)×m ,
A(τ )c = 0,
(43)
here c = (ω11 , ω12 , ω13 , ω22 , ω23 , ω33 )T . In theory, the rank of matrix A is 5, therefore our problem is to solve the following minimization: σ1 ×m min J = (τ )|τ ∈ (−∞, +∞) , (44) σ5 where σ1 ≥ σ2 ≥ · · · ≥ σ5 ≥ σ6 are singular values of A(τ ). For solving the above problem, we firstly compute the gradient of J . By taking the singular value decomposition of A, A = U DV T , then σj = uTj Av j where uj , v j are the j th row of the orthogonal matrices U and V , respectively. And thus we have ∂(uTj Av j ) ∂σj = ∂τi ∂τi =
∂uTj ∂τi
= σj Notice have
Av j + uTj
∂uTj ∂τi
∂v Tj ∂τi
∂A ∂v n v j + uTj A ∂τi ∂τi
uj + uTj
v j = uTj
∂uj ∂τi
∂v j ∂A v j + σj v Tj . ∂τi ∂τi = 0 from uTj uj = v Tj v j = 1, we
∂σj ∂A = uTj vj . ∂τi ∂τi Hence the gradient of J on τi is: ∂A ∂A ∂J 1 = 2 σ5 uT1 v 1 − σ1 uT5 v5 . ∂τi ∂τi ∂τi σ5
(45)
Int J Comput Vis
(a) Condition numbers of CCF-Is
(b) Condition numbers of CCF-IIs
Fig. 3 Condition numbers of the CCF-Is and CCF-IIs corresponding to 5 infinite homographies H∞ = KRK −1 where K : fu = fv = 1000, s = 0.0, pu = pv = 512 and the rotation Ri = Rx (θi1 ) · Ry (θi2 ) · Rz (π/θi3 ) (i = 1, 2, 3, 4, 5), θij is randomly generated by the uniform distribution on the interval [−π/4, π/4]. The abscissa and the ordinate are parameter τ and condition number respectively. We can see that: For the CCF-I, the condition numbers of equations vary with the parameter τ and the condition number reaches its minimum at τ = 0 i.e.,
the SCE has the smallest condition number among the CCF-I. For the CCF-II, the condition numbers of equations vary with the parameter τ also, but the condition numbers reach their minimums at different τ . The condition numbers of the CCF-I and the CCF-II for the same parameter are different and the CCF-II has always a smaller condition number, thus the numerical performance of CCF-II is better than that of CCF-I in general
With the computed gradient ∇J (τ ) of J , we can use some numerical methods to solve the minimization problem (44), such as steepest descent method or conjugate gradient method (Press et al. 1988). In the experiments of this paper, we use the steepest descent method as outlined below.
1. Let s0 = 0, s1 = 1, n = 1 and given a termination δ. 2. Compute e = J (τ k + sn−1 p k ) − J (τ k + sn p k ). 3. In the case e > 0, set up sn+1 = sn + 10λ and let n ← n + 1, return to Step 2. In the case e < 0, if |sn+1 − sn | < δ, terminate iteration and put out sn , or else do Step 4. 4. (a) t = λ/10. (b) Set up sn+1 = sn + t and do Step 2. (c) In the case e > 0, let n ← n + 1 and return to (b). In the case e < 0, if |sn+1 − sn | < δ, terminate iteration and put out sn , or else let t ← t/10 and return to (b).
Determining the Parameters τ The algorithm consists of the following four steps: 1. Given an initialization τ 0 and a termination ε. Let k = 0. 2. Compute the gradgent ∇J (τ k ) of J with (45), ∂J ∂J ∂J , ,..., . ∇J (τ k ) = ∂τ1 ∂τ2 ∂τm τ =τ k If ∇J (τ k ) < ε, terminate iteration and put out τ k , or else do Step 3. 3. Let p k = −∇J (τ k ). Solve sk such that J (τ k + sk pk ) = min J (τ k + sp k )|s ≥ 0 . 4. Let τ k+1 = τ k + sk p k and k ← k + 1. Return to Step 2. Remark 6 First, in this paper the initialization is set to τ 0 = (1, 1, . . . , 1), corresponding to the set of CCEs. Second, in Step 3, the exact solution s cannot be obtained since J (τ k + spk ) has no analytical expression. We can use an algorithm similar to Levenberg-Marquardt algorithm (Press et al. 1988) to find an approximate solution of sk as follows:
In our work, δ = 10−5 and λ = 0.1 are used. Third, in practice finding an exact solution of the minimization (44) is unnecessary. This is because the condition number is a qualitative analysis of numerical stability; the condition numbers with the same order of magnitude are unable to lead to a big difference in numerical stability. And thus, we can terminate iteration once the condition number is less than a given threshold. In this work, 50 is used as the threshold. 5.4 Cayley Algorithm By the discussions in Sects. 5.1, 5.2 and 5.3, we can also give a calibration algorithm similar to the standard one, called the Cayley algorithm in this paper. Goal Given N views of a scene and assume the infinite homographies Hi∞ between the 0th and ith views are known. Compute the intrinsic matrix K and the rotation Ri between the 0th and ith views.
Int J Comput Vis
Algorithm The algorithm consists of the following four steps: 1. Let H˜ i∞ = det−1/3 (Hi∞ )Hi∞ and compute its Cayley transformation Hi = ϕ(H˜ i∞ ) and CCF-I (or CCF-II): ϕ T (τi Hi )ωϕ(τi Hi ) = ω,
i = 1, 2, . . . , N
(or ϕ T (τi Hi )ω = ωϕ(−τi Hi ), i = 1, 2, . . . , N ). 2. Determine the parameters τ = (τ1∗ , τ2∗ , . . . , τN∗ ) by solving the minimization (44), then find the least squares solution of the set of linear equations on ω: ϕ T (τi∗ Hi )ωϕ(τi∗ Hi ) = ω,
i = 1, 2, . . . , N
(or ϕ T (τi∗ Hi )ω = ωϕ(−τi∗ Hi ), i = 1, 2, . . . , N ). 3. Determine the intrinsic matrix K by taking the Cholesky factorization of ω. 4. Let ϕK (τi∗ H ) = K −1 ϕ(τi∗ Hi )K and compute
× ((τi∗ − 1)I + (τi∗ + 1)ϕK (τi∗ Hi )).
6.1 Generalized Cayley Transformation If a function of complex variable, ξ = φ(λ), satisfies the following three conditions, then it is called a generalized Cayley transformation: (a) φ(λ) is an analytical function with real coefficients in |λ| < δ; (b) φ(λ)φ(−λ) = 1 for |λ| < δ; (c) φ(0) > 0 and φ (0) = 0. It is easy to see that the Cayley transformation is a special case of the generalized Cayley transformation. Below are two examples of generalized Cayley transformations:
(49)
Then for −∞ < τ < ∞ the transformation (48) is a generalized Cayley transformation of matrix A. If r(A) ≤ 1 (here r(A) = max{|λ| : λ ∈ σ (A)} is called the spectrum radius of matrix A), then for −π/2 < τ < π/2, the transformation (49) is also a generalized Cayley transformation of matrix A. For generalized Cayley transformations, we have the following proposition. Proposition 13 Assume φ(τ A) (−ε < τ < ε) is a generalized Cayley transformation of matrix A ∈ C n×n and C ∈ C n×n is an arbitrary matrix. Then φ(τ A)T Cφ(τ A) = C
for 0 < |τ | < ε
AT C + CA = 0.
(50)
∂ ∂ g(τ A) = φ(τ A)T Cφ(τ A) ∂τ ∂τ ∂ φ(τ AT )Cφ(τ A) = ∂τ
Assume equation φ(τ A)T Cφ(τ A) = C holds for 0 < |τ | < ε, then for 0 < |τ | < ε there must exist AT φτ (τ AT )Cφ(τ A) + φ(τ AT )CAφτ (τ A) = 0.
(52)
Because φ(λ) is analytical in −ε < τ < ε, the left side in (52) is a continuous mapping on the variable τ . Thus, letting τ → 0 in (52) leads to φ (0)φ(0)(AT C + CA) = 0. By the conditions (b) and (c), φ (0)φ(0) = 0, and thus there must be AT C + CA = 0.
1 − sin(λ) φs (λ) = . 1 + sin(λ)
For a matrix A ∈ C n×n , a generalized Cayley transformation of A can be defined as (46)
Now assume equation AT C + CA = 0 holds, then we can verify that φ(τ AT )C = Cφ(−τ A),
φτ (τ AT )C = Cφτ (−τ A)
holds for any analytical function φ(λ). Then we obtain AT φτ (τ AT )Cφ(τ A) = AT Cφτ (−τ A)φ(τ A),
where ε is a positive number such that σ (εA) ⊂ {λ : |λ| < δ}.
φs (τ A) = (I − sin(τ A))(I + sin(τ A))−1 .
= AT φτ (τ AT )Cφ(τ A) + φ(τ AT )CAφτ (τ A). (51)
6 Generalized Cayley Calibration Families
(−ε < τ < ε),
(48)
Proof Let g(τ A) = φ(τ A)T Cφ(τ A), then,
Then, determine the rotation Ri by finding the best rotation approximation of the matrix R˜ i , as shown in the standard algorithm in Sect. 4.3.
Bτ = φ(τ A)
φe (τ A) = (I + e−τ A )(I + eτ A )−1 ,
⇔
Ri = ((τi∗ + 1)I + (τi∗ − 1)ϕK (τi∗ Hi ))−1
1 + e−λ , φe (λ) = 1 + eλ
For example, let
(47)
φ(τ AT )CAφτ (τ A) = CAφ(−τ A)φτ (τ A).
Int J Comput Vis
So, by (51), we have ∂ g(τ A) = AT Cφτ (−τ A)φ(τ A) + CAφ(−τ A)φτ (τ A). ∂τ
Proof Let R be the rotation between the two given views and its eigenvalues factorization be:
By the condition (b), there is φτ (τ A)φ(−τ A) = φ(τ A)φτ (−τ A).
R = V diag(1, eiθ , e−iθ )V −1 ,
And thus, for 0 < |τ | < ε,
where θ is its rotation angle, V = (v 1 , v 2 , v 3 ) is an unitary matrix and v 1 is the rotation axis of R. Then according to Proposition 2,
∂ g(τ A) = (AT C + CA)φ(τ A)φτ (−τ A) = 0. ∂τ Therefore, for 0 < |τ | < ε,
τ H = Kτ ϕ(R)K −1
g(τ A) = φ(τ A)T Cφ(τ A) = constant matrix.
= KVdiag(0, iτ tan(θ/2), −iτ tan(θ/2))V −1 K −1 .
By the continuity of φ(τ A) at τ = 0, we have
Then,
φ(τ A)T Cφ(τ A) = g(τ A) = φ 2 (0)C
φ(τ H ) = K V diag(φ(0), φ(is(τ )), φ(−is(τ )))V −1 K −1 ,
for 0 < |τ | < ε.
Q(τ )
By the condition (b), φ 2 (0) = 1, and thus there must be φ(τ A)T Cφ(τ A) = C.
6.2 Generalized Cayley Calibration Family Applying a generalized Cayley transformation to the Cayley transformation of the infinite homography, we can obtain the following proposition according to the Proposition 13. Proposition 14 The SCE is equivalent to the following oneparameter family φ(τ H )T ωφ(τ H ) = ω
Proposition 15 The generalized Cayley transformation φ(τ H ) is an infinite homography and its rotation has the same rotation axis with the rotation between the two views.
for 0 < |τ | < ε.
(53a)
Given a generalized Cayley transformation, by (53a) we obtain a one-parameter family equivalent to the SCE. Since it has the same form with the CCE-I and CCF-I, it is called the generalized Cayley calibration family of the first type (GCCF-I). Because φ −1 (τ H ) = φ(−τ H ) for −ε < τ < ε, the SCCF-I has the following equivalent form, it is called the generalized Cayley calibration family of the second type (GCCF-II): φ(τ H )T ω = ωφ(−τ H ) for 0 < |τ | < ε.
(53b)
As in the CCFs, the numerical stability of equations in a GCCF varies also with the parameter τ . For a given generalized Cayley calibration family, we can use the method in Sect. 5.3 to reduce the condition number of calibration equations by selecting a proper parameter value and can give a calibration algorithm similar to the algorithm in Sect. 5.3, called as generalized Cayley algorithm. Next, we give an explanation of (53a) also in terms of the infinite homography.
here s(τ ) = τ tan(θ/2). In order to prove Proposition 15, we need only to prove that Q(τ ) is a rotation and it has the same rotation axis with R. From the conditions (b) and (c) of the function φ, there is φ(0) = 1. By the conditions (b) of the function φ, we have φ(is(τ )) · φ(−is(τ )) = 1. Then from the conditions (a) of the function φ, we have φ(−is(τ )) = φ(is(τ )). Therefore {φ(is(τ )), φ(−is(τ ))} is a conjugate pair of unit complex numbers, and Q(τ ) can be written as: Q(τ ) = V diag(1, eiσ (τ ) , e−iσ (τ ) )V −1 , where σ (τ ) =
1 1 log(φ(is(τ ))) = log(φ(iτ tan(θ/2))). i i
Hence, Q(τ ) is a rotation and v 1 is its rotation axis.
Remark 7 The GCCF-I contains also the CCE-I. Let θ be the rotation angle between the two given views, then in the GCCF-I the parameter τ corresponding to the CCE-I is the solution of the following equation: log(φ(iτ tan(θ/2))) = iθ. 6.3 S-Cayley Calibration Family As an example, we examine the generalized Cayley calibration family induced by the transformation (49): B T (τ H )ωB(τ H ) = ω
for 0 < |τ |
Cayley Transformation and Numerical Stability of Calibration Equation F.C. Wu · Z.H. Wang · Z.Y. Hu
Received: 25 May 2008 / Accepted: 22 October 2008 © US Government 2008
Abstract The application of Cayley transformation to enhance the numerical stability of camera calibration is investigated. First, a new calibration equation, called the standard calibration equation, is introduced using the Cayley transformation and its analytical solution is obtained. The standard calibration equation is equivalent to the classical calibration equation, but it exhibits remarkable better numerical stability. Second, a one-parameter calibration family, called the Cayley calibration family which is equivalent to the standard calibration equation, is obtained using also the Cayley transformation and it is found that this family is composed of those infinite homographies whose rotation has the same axis with the rotation between the two given views. The condition number of equations in the Cayley calibration family varies with the parameter value, and an algorithm to determine the best parameter value is provided. Third, the generalized Cayley calibration families equivalent to the standard calibration equation are also introduced via generalized Cayley transformations. An example of the generalized Cayley transformations is illustrated, called the S-Cayley calibration family. As in the Cayley calibration family, the numerical stability of equations in a generalized Cayley calibration family also depends on the parameter value. In addition, a more generic calibration family is also proposed and F.C. Wu () · Z.H. Wang · Z.Y. Hu National Laboratory of Pattern Recognition, Institute of Automation, Chinese Academy of Sciences, PO Box 2728, Beijing 100190, People’s Republic of China e-mail: [email protected] Z.H. Wang e-mail: [email protected] Z.Y. Hu e-mail: [email protected]
it is proved that the standard calibration equation, the Cayley calibration family and the S-Cayley calibration family are all some special cases of this generic calibration family. Keywords Camera calibration · The absolute conic · Calibration equation · Cayley transformation · Numerical stability
1 Introduction Camera calibration is a necessary step to recover 3D metric information from 2D images. Many methods and techniques for camera calibration have been proposed in the last decades. The existing methods can be roughly classified into two categories: methods based on a reference calibration object and self-calibration methods. In the first kind, there are mainly the methods based on 3D object (Brown 1971; Faig 1975; Tsai 1986; Faugeras 1993), 2D pattern (Zhang 1999, 2000; Sturm and Maybank 1999) and 1D object (Zhang 2004; Wu et al. 2005; Wang et al. 2007). In these methods, camera calibration is performed by observing a 3D calibration object, a 2D calibration pattern or a 1D calibration object. Since the Euclidean geometries of calibration objects are used, the methods can provide a calibration of high-precision. In the second kind, the methods widely used are the calibration based on the Kruppa equation (Faugeras et al. 1992; Maybank and Faugeras 1992; Hartley 1997a; Luong and Faugeras 1997), the absolute conic (Pollefeys et al. 1996; Hartley 1997b; Pollefeys and Gool 1999) and the dual absolute quadric (Triggs 1997; Ponce et al. 2000). These methods do not need any reference calibration objects and use only the correspondences of image entities, such as points, lines. It is remarkable and desirable that the self-calibration methods require no reference
Int J Comput Vis
objects. Regrettably, in addition to their low calibration precision, the numerical stabilities of self-calibration methods are usually poor. Hartley and Zisserman (2000) pointed out that the algorithms absolutely depended on self-calibration should be restrainedly used in practice. Though the people have paid great attention to improve the precision and the numerical stability of self-calibration methods, no satisfactory results are obtained up to now. In this paper, we mainly investigate the numerical stabilities of the self-calibration methods based on the absolute conic. The calibration equation induced by the abT ωH = ω or, equivalently H T ω = ωH −1 , solute conic, H∞ ∞ ∞ ∞ is called the classical calibration equation in this paper. Here H∞ is the infinite homography and ω the image of the absolute conic. The previous studies are mainly focused on the conditions and computational means to solve the infinite homography (such as, Pollefeys et al. 1996; Pollefeys and Gool 1999; Hartley 1997b, etc.). For a long time the numerical stability of the calibration equation itself has been neglected in research community. In computation term, generally speaking different forms of the same equation can usually have different numerical stabilities. Thus we have the following problem: For the self-calibration based on the absolute conic, do there exist more stable forms than the classical calibration equation in computation? This problem is answered affirmatively in this paper. More specifically, we use the Cayley transformation to enhance the numerical stability of camera calibration and here are our main contributions: • Using the Cayley transformation of the infinite homography, a new calibration equation, called the standard calibration equation, is introduced and its analytical solution is obtained. The Cayley transformation of the infinite homography transforms image point to its conjugate point on the image of absolute conic. The standard calibration equation is equivalent to the classical calibration equation, but it exhibits remarkable better numerical stability. • A one-parameter family equivalent to the standard calibration equation, called the Cayley calibration family, is introduced using also the Cayley transformation and it is shown that the Cayley calibration family is composed of those infinite homographies whose rotation has the same axis with the rotation between the two given views. The condition number of equations in the Cayley calibration family varies with the parameter value, and an algorithm is also provided to determine the best parameter which gives the best numerical stability among the Cayley family. • Generalized Cayley calibration family equivalent to the standard calibration equation is introduced via generalized Cayley transformations. Given a generalized Cayley transformation, a generalized Cayley calibration family can be obtained. Similarly to the Cayley calibration
family, the numerical stability of equations in this general family varies also with the parameter value. An example of the generalized Cayley transformations is illustrated, called the S-Cayley calibration family. • A more generic calibration family is proposed and it is proved that the standard calibration equation, the Cayley calibration family and the S-Cayley calibration family are all some special cases of this generic calibration family. This paper is organized as follows. Section 2 reviews briefly the classical calibration algorithm based on the absolute conic. In the Sect. 3, the Cayley transformation is introduced. The standard calibration equation, the Cayley calibration family, and the generalized Cayley calibration family are elaborated in Sects. 4, 5 and 6, respectively. Section 7 gives more generic equivalent equations. Section 8 reports the experiment results. Section 9 concludes this paper.
2 Classical Calibration Equation The infinite homography between the two views is denoted as H∞ in this paper. If the motion between the two views is (R, t), R is the rotation and t the translation, and the intrinsic matrix of camera is K, then the infinite homography can be expressed as H∞ = KRK −1 .
(1)
Let ω be the image of the absolute conic (IAC) i.e., ω = K −T K −1 . Then there is the classical calibration equation (CCE): T ωH∞ = ω. H∞
(2)
Evidently, this equation can be changed to the following equivalent form: T −1 ω = ωH∞ . H∞
(3)
For the convenience of statement, (2) is called CCE of the first type (CCE-I), and (3) is called CCE of the second type (CCE-II). Hartley (1997b) used the CCE-II to compute ω and pointed out that the CCE can only provide 4 independent constraints on ω. Thus at least 3 views are necessary to linearly solve ω. After ω is determined, the intrinsic matrix K can be determined by the Cholesky factorization of ω, then the rotation R can in turn be determined. The classical calibration algorithm can be sketched as: Classical Algorithm Goal Given N views of a scene and assume the infinite homographies Hi∞ between the 0th and ith views are known. Compute the intrinsic matrix K and the rotation Ri between the 0th and ith views.
Int J Comput Vis
Algorithm The algorithm consists of the following three steps: 1. Let H˜ i∞ = Det−1/3 (Hi∞ )Hi∞ . Compute the least squares solution of the following set of linear equations on ω: T H˜ i∞ ωH˜ i∞ = ω
−1 T (or H˜ i∞ ω = ωH˜ i∞ ),
i = 1, 2, . . . , N − 1. 2. Determine the intrinsic matrix K by taking the Cholesky factorization of ω. 3. Determine the rotation Ri by finding the best rotaK = K −1 H˜ tion approximation of H˜ i∞ i∞ K, i.e., compute firstly the singular value decomposition of the K , H˜ K = U diag(t , t , t )V T , then set up matrix H˜ i∞ 1 2 3 i∞ T Ri = U V . Remark 1 Since the infinite homography computed in practice is equal to matrix KRK −1 up to a scale factor, it is necessary that the determinant of the computed infinite homography be normalized to 1 in the step 1, otherwise (2) and (3) do not hold. Computing Infinite Homography Here, we sketch three different methods for the infinite homography computation. Method 1 Computation with infinite correspondences (i.e., image point correspondences from infinite plane). Given two views of a scene and a set of infinite correspondences: (1) m(0) j ↔ mj ,
j = 1, 2, . . . , M.
We can linearly compute the infinite homography H∞ from the following set of linear equations: (1)
Here e is the epipole on the 2nd view, F T e = 0, a is an unknown 3-vector. Given an infinite correspondence m(0) ↔ m(1) , we have: [m(1) ]× ([e ]F + e a T )m(0) = 0. In the above three equations on a, there are only one independent constraint. And thus, at least 3 infinite correspondences are necessary to determine linearly the vector a. Method 3 Computation using the modulus constraint. For more than 2 views, using the modulus constraint shown in Pollefeys et al. (1996), i.e., the three eigenvalues of infinite homography have the same modulus, the infinite homographies between views can also be derived without resort to the infinite correspondences. But it needs to solve a set of quartic equations in three variables, and thus is highly non-linear. The computation steps are as follows. First, by the trifocal tensor method in Hartley and Zisserman (2000) or the matrix factorization methods in Sturm and Triggs (1996), Oliensis and Hartley (2007), a projective reconstruction is computed with a set of image correspondences from 3D space points: P0 = (I |0),
Pi = (Hi |ei )
(j = 1, 2, . . . , N − 1).
Then, the infinite homography between the 0th view and ith view is expressed as: Hi∞ = Hi + ei a T
(i = 1, 2, . . . , N − 1).
Here a is an unknown 3-vector, the direction of infinite plane in the reconstructed projective space. And thus, the vector a can be non-linearly computed by using the modulus constraints.
(0)
[mj ]× H∞ mj = 0 (j = 1, 2, . . . , M). This method needs at least 4 infinite correspondences because each infinite correspondence can only give rise to two independent constraints on H∞ . Hartley (1997b) pointed out that if the motion between two views is a rotation around the camera center, then infinite correspondences in the above can be replaced by image correspondences from 3D space points. Method 2 Computation with the epipolar geometry. Given the foundational matrix between two views, F , which can be computed with image correspondences from 3D space points, three infinite correspondences are needed in order to linearly compute the infinite homography H∞ . This is because the infinite homography H∞ can be expressed as: H∞ = [e ]F + e a T .
3 Cayley Transformation The following function of complex variable is called the Cayley transformation (Bell 1992): ξ = ϕ(λ) =
1−λ . 1+λ
(4)
It is a conformal mapping on the extended complex plane. It is ready to verify that λ = ϕ −1 (ξ ) =
1−ξ = ϕ(ξ ). 1+ξ
So, the inverse transformation is also the Cayley transformation. For a matrix A ∈ C n×n , let σ (A) denote the set of all eigenvalues of A. Assume −1 ∈ / σ (A) i.e., det(I + A) = 0,
Int J Comput Vis
then the Cayley transformation of A is defined by
By right-multiplying the two sides of this equation by C, and notice that B T C + CB = 0, have
B = ϕ(A) = (I − A)(I + A)−1 −1
= (I + A)
(I − A).
(5)
AT C(I − B) = C(I + B). So,
From this equation, we have (I + A)B = (I − A) or, equivalently, A(I + B) = (I − B).
(6)
If Bx = −x, then by (5) there must be x = 0, and thus, −1 ∈ / σ (B). Therefore, from (6) we have
= (I + B)−1 (I − B).
(7)
This indicates that the matrix A is also the Cayley transformation of B. According to the above discussions, the Cayley transformation of matrices is a one-to-one mapping on the matrix set {A ∈ C n×n | − 1 ∈ / σ (A)} and we can prove that: ⇔
β=
= C(I + B)(I + B)−1 = C. Hence, the necessary condition holds too.
4 Standard Calibration Equation
A = ϕ(B) = (I − B)(I + B)−1
λ ∈ σ (A)
AT CA = AT C(I − B)(I + B)−1
1−λ ∈ σ (ϕ(A)). 1+λ
(8)
For the Cayley transformation of matrices, we have the following proposition, which is the fundamental basis of our current work. Proposition 1 Assume matrices A, B ∈ satisfy (5) and (7) and let C ∈ C n×n be an arbitrary matrix, then, C n×n
4.1 Cayley Transformation of H∞ The infinite homography H∞ satisfies det(I + H∞ ) = det(K(I + R)K −1 ) = det(I + R). If the rotation angle of R is not equal to iπ (i ∈ Z), det(I + R) = 0. So, in this case, there is the Cayley transformation of H∞ : H = ϕ(H∞ ) = (I − H∞ )(I + H∞ )−1 = (I + H∞ )−1 (I − H∞ ).
(10)
And H∞ is also the Cayley transformation of H : H∞ = ϕ(H ) = (I − H )(I + H )−1 = (I + H )−1 (I − H ).
(11)
(9)
In this paper, we always assume that the rotation angle θ of R is not equal to iπ (i ∈ Z).
Proof By left-multiplying the two sides of (5) by (I + A), and transposed, we obtain
Proposition 2 H is a conjugate of an anti-symmetric matrix.
B T (I + A)T = (I − A)T .
Proof Substituting H∞ = KRK −1 into (10) yields
AT CA = C
⇔
B T C + CB = 0.
By right-multiplying the two sides of this equation by CA, and notice that AT CA = C, we have B T C(I + A) = −C(I − A).
H = (I − KRK −1 )(I + KRK −1 )−1 = K (I − R)(I + R)−1 K −1 = Kϕ(R)K −1 .
(12)
ϕ(R)
So,
Where ϕ(R) is the Cayley transformation of rotation R. Since
B T C = −C(I − A)(I + A)−1 = −CB.
(ϕ(R))T = (I + R T )−1 (I − R T )
Therefore the sufficient condition holds. Now, we prove the necessary condition next. By left-multiplying the two sides of (7) by (I + B), and transposed, we obtain AT (I + B)T = (I − B)T .
= (I + R −1 )−1 (I − R −1 ) = (R + I )−1 RR −1 (R − I ) = (R + I )−1 (R − I ) = −ϕ(R),
Int J Comput Vis
ϕ(R) is an anti-symmetric matrix, and thus by (12) the proposition holds. {1, eiθ , e−iθ }.
Let the eigenvalues of rotation R be Then it is not difficult to see that ϕ(R) has the eigenvalues {0, i tan(θ/2), −i tan(θ/2)} and has the same eigenvectors with R, i.e., ⎧ ⎧ ⎨ϕ(R)x = 0 ⎨ Rx = x ϕ(R)y = i tan(θ/2)y ⇔ Ry = eiθ y ⎩ ⎩ ϕ(R)y = i tan(θ/2)y ϕ(R)y = e−iθ y. According to Proposition 2, H has also the eigenvalues {0, i tan(θ/2), −i tan(θ/2)} and the corresponding eigenvectors are {Kx, Ky, Ky}. Corollary 1 (a) rank(H ) = 2. (b) Let H K = K −1 H K. Then the rotation R is the Cayley transformation of H K : R = ϕ(H K ) = (I − H K )(I + H K )−1 .
(13)
Proposition 3 Let the eigenvalues of H be {0, is, −is}, then, (a) H can be factorized as H = P [a]× P −1 ,
(14)
where P is a real matrix, [a]× is the anti-symmetric matrix defined by the vector a = (s, 0, 0)T , i.e. ⎛ ⎞ 0 0 0 [a]× = ⎝ 0 0 s ⎠ . 0 −s 0 Equation (14) is called the standard anti-symmetric factorization of H in this paper. (b) If H has another standard anti-symmetric factorization H = P [a]× P −1 , then P must be of the form P = P W (α, β, γ ), where ⎛
α ⎝ W= 0 0
0 β −γ
(15) ⎞
0 γ⎠ β
is an invertible matrix. ¯ be the unit eigenvectors of H correProof Let {p1 , q, q} sponding to the eigenvalues {0, is, −is}, i.e. ⎧ ⎨ H p1 = 0 H q = isq ⎩ ¯ H q¯ = −is q.
Here q¯ is the conjugate of q. Then, we obtain ¯ = (p 1 , q, q)diag(0, ¯ is, −is). H (p 1 , q, q) By right-multiplying the two sides of this equation by ⎞ ⎛ 1 0 0 √ √ ⎟ ⎜ U = ⎝ 0 1/ 2 1/i 2 ⎠ , √ √ 0 1/ 2 −/i 2 since
q + q¯ q − q¯ ¯ (p 1 , q, q)U ∈ R 3×3 = p1 , √ , √ 2 i 2
and
¯ ¯ [a]× is, −is)U = (p 1 , q, q)U (p 1 , q, q)diag(0, q + q¯ q − q¯ = p1 , √ , √ [a]× , 2 i 2 we have q + q¯ q − q¯ q + q¯ q − q¯ H p1 , √ , √ = p1 , √ , √ [a]× 2 i 2 2 i 2 P
P
and thus (14) holds. We now prove (b). Assume H = P [a]× P −1 is another standard anti-symmetric factorization, then P [a]× P −1 = P [a]× P −1 , and thus, −1 −1 P P[a]× − [a]× P P = 0. W
W
According to Lemma 1 below, solutions of W are ⎞ ⎛ x 0 0 ⎟ ⎜ W = xI + y[a]× + z[a]2× = ⎝ 0 x − zs 2 ys ⎠ , 0 −ys x − zs 2 x(y 2 s 2 + (x − zs 2 )2 ) = 0. Let α = x, β = x − zs 2 , γ = ys, then, W = W (α, β, γ ) = P −1 P . Hence, P = P W .
Lemma 1 Assume c is a nonzero 3-vector, then solutions of the equation [c]× X − X[c]× = 0 are: X = xI + y[c]× + z[c]2× ,
x, y, z ∈ R.
(16)
Proof Suppose A, B ∈ R n×n and let vec(B) = (bT1 , bT2 , . . . , bTn )T , where bTi is the ith row of B, then we can obtain vec(AB) = (A ⊗ In )vec(B),
(17)
Int J Comput Vis
vec(BA) = (In ⊗ AT )vec(B).
(18)
Here ⊗ denotes the tensor product of matrices. By applying (17) and (18), [c]× X − X[c]× = 0 can be changed as ([c]× ⊗ I3 + I3 ⊗ [c]× )vec(X) = 0.
(19)
Since the three eigenvalues of anti-symmetric matrix [c]× are {0, i c 2 , −i c 2 }, the coefficient matrix of (19) has eigenvalue zero of multiplicity 3. It is easy to verify that {vec(I ), vec([c]× ), vec([c]2× )} are three linear independent solutions of (19). And thus its solutions are
(b) Symmetrical positive-definite solutions of the SCE are: ω = xP −T diag(1, w, w)P −1 ,
x > 0, w > 0.
(22)
Proof Substituting H = P [a]× P −1 into (20) and by some simple manipulations, we have −[a]× P T ωP + P T ωP [a]× = 0.
(23)
By denoting X = P T ωP in (23), we have [a]× X − X[a]× = 0.
(24)
According to Lemma 1, its solutions are vec(X) = xvec(I ) + yvec([c]× ) + zvec([c]2× ),
X = xI + z[a]× + y[a]2× ,
x, y, z ∈ R.
x, y, z ∈ R.
Since It yields that the solutions of [c]× X − X[c]× = 0 are X = xI + y[c]× + z[c]2× ,
x, y, z ∈ R.
X − XT = 0
⇔
− (xI + z[a]× + y[a]2× )T = 0
4.2 Standard Calibration Equation ⇔ By Proposition 1, the following proposition holds.
(20)
This equation is called the standard calibration equation (SCE) in this paper, which has the similar form with the CCE-II. Though the CCE-I and the SCE are equivalent to each other, their forms have remarkable differences. First, the CCE-I is of a symmetrical form, while the SCE is of an anti-symmetric form because it implies that the matrix ωH is anti-symmetric. Second, the CCE-I is homogeneous on only the matrix ω, while the SCE is homogeneous on not only the matrix ω but also the matrix H . Third, the coefficients of ω are quadratic function of the H∞ ’s elements in the CCE-I, and are linear function of the H ’s elements in the SCE. Besides these, the SCE has a nice property as stated in the following proposition. Proposition 5 Let H = P [a]× P −1 is a standard antisymmetric factorization of H (see Proposition 3), then, (a) Symmetrical solutions of the SCE are: ω = P −T diag(x, x − ys 2 , x − ys 2 )P −1 , x, y ∈ R.
2z[a]× = 0
⇔
z = 0,
symmetrical solutions of (24) are
Proposition 4 The CCE-I (or CCE-II) is equivalent to the equation below: H T ω + ωH = 0.
(xI + z[a]× + y[a]2× )
(21)
X = xI + y[a]2× ,
x, y ∈ R.
Notice that [a]2× = diag(0, −s 2 , −s 2 ), we have X = diag(x, x − ys 2 , x − ys 2 ),
x, y ∈ R.
Hence symmetrical solutions of the SCE are ω = P −T diag(x, x − ys 2 , x − ys 2 )P −1 ,
x, y ∈ R.
Evidently, ω is a symmetrical positive-definite matrix if and only if X is a symmetrical positive-definite matrix, and thus symmetrical positive-definite solutions of the SCE are ω = P −T diag(x, x − ys 2 , x − ys 2 )P −1 ,
x > 0, y < xs −2 .
Let w = 1 − x −1 ys 2 , we have ω = xP −T diag(1, w, w)P −1 ,
x > 0, w > 0.
Remark 2 The symmetrical solutions (21) and symmetrical positive-definite solutions (22) of the SCE do not depend on the specific selection of matrix P in the standard antisymmetric factorization of H . This is because: Assume H = P [a]× P −1 is another standard anti-symmetric factorization, then according to Proposition 3(b), P = P W (α, β, γ ). It is not difficult to verify that W diag(a, b, b) = diag(a, b, b)W,
∀a, b
Int J Comput Vis
or, equivalently, W diag(a, b, b)W −1 = diag(a, b, b),
Proof By Proposition 5, the symmetrical positive-definite solutions of HiT ω + ωHi = 0 are
∀a, b.
And thus, P −T diag(a, b, b)P −1 = P −T diag(a, b, b)P −1 . Hence, the symmetrical solutions (21) and symmetrical positive-definite solutions (22) do not depend on the selection of matrix P in the standard anti-symmetric factorization. Proposition 5 tells us that from the SCE of a pair of views, the IAC ω can be determined up to a one-parameter family. Thus under the 4-parameter camera model with zero skew, the IAC ω can be uniquely determined from the SCE of a pair of views and its computation is very easy with (22): q11 q12 , ω = P −T diag 1, − q21 q22 + q31 q32 q11 q12 P −1 , − q21 q22 + q31 q32 where qij is the (i, j )-element of P −1 . Given two SCEs from three views, the following proposition provides an analytical solution of ω. Proposition 6 Assume Ri is the rotation between the 0th and ith views and the corresponding SCE is: HiT ω + ωHi
= 0,
i = 1, 2.
(25)
Then (25) has a unique solution if and only if that the rotation axes of the two rotations are not parallel. And if we let the standard anti-symmetric factorization of Hi be Hi = Pi [ai ]× Pi−1 , then its unique solution is: ω = P1−T diag(1, w1 , w1 )P1−1 = x2 P2−T diag(1, w2 , w2 )P2−1 , where 2 (Q21 + Q231 ) 3j =2 (Q22j + Q23j ) , w1 = 3 2 + q2 ) 2 2 (q21 j =2 (q2j + q3j ) 31 Q211 + w12 (Q212 + Q213 ) x2 = , 2 + q2 + q2 ) det2 (Q)(q11 12 13 Q221 + Q231 w2 = , 2 + q2 ) det2 (Q)(q21 31 (qij ) = (q 1 , q 2 , q 3 )
T
= P2−1 P1 = Q,
(Qi1 , Qi2 , Qi3 )T = q j × q k , i = 1, 2, 3; j, k = i and j < k.
(26)
ω = xi Pi−T diag(1, wi , wi )Pi−1 (xi > 0, wi > 0),
i = 1, 2. (27)
Assume x1 = 1. Then from (27), we have diag(1, w1 , w1 ) = x2 QT diag(1, w2 , w2 )Q = x2 (q 1 q T1 + w2 (q 2 q T2 + q 3 q T3 )). It is easy to see that this equation is equivalent to the following equation: ⎧ ⎪ ⎨diag(1, w1 , w1 )(q 2 × q 3 ) = x2 det(Q)q 1 diag(1, w1 , w1 )(q 1 × q 3 ) = −w2 det(Q)q 2 (28) ⎪ ⎩ diag(1, w1 , w1 )(q 1 × q 2 ) = w2 det(Q)q 3 . (a) If the rotation axes of the two rotations are parallel to each other, then from the proof of Propositions 2 and 3, the first rows of P1 and P2 are equal up to a nonzero scale. And thus Q can be expressed as ⎛ ⎞ α ∗ ∗ Q = P2−1 P1 = ⎝ 0 ∗ ∗ ⎠ , α = 0. 0 ∗ ∗ Hence, (28) is equivalent to the equation below: ⎧ ⎨det(Q)x2 = α Q23 w1 = −q23 w2 Q22 w1 = −q22 w2 , ⎩ Q32 w1 = −q32 w2 , Q33 w1 = −q33 w2 . Clearly, its solution is not unique, and thus the solution of (25) is not unique either. (b) If the rotation axes of two rotations are not parallel, then the two elements q21 and q31 of Q are not all zero. Rewriting the 2nd and 3rd equations in (28) as (Q21 , Q22 w1 , Q23 w1 ) = − det(Q)(q21 w2 , q22 w2 , q23 w2 ) (Q31 , Q32 w1 , Q33 w1 ) = det(Q)(q31 w2 , q32 w2 , q33 w2 ). 2 + q 2 = 0, we have Since q21 31 Q221 + Q231 w2 = , 2 2 + q2 ) det (Q)(q21 31 2 (Q21 + Q231 ) 3j =2 (Q22j + Q23j ) . w1 = 3 2 + q2 ) 2 2 (q21 j =2 (q2j + q3j ) 31
Then, from the first equation in (28) we have Q211 + w12 (Q212 + Q213 ) x2 = . 2 + q2 + q2 ) det2 (Q)(q11 12 13 Hence, (25) has the unique solution (26).
Int J Comput Vis
It is well-known that having the parallel rotation axis is a degenerate configuration of the set of CCEs (Hartley and Zisserman 2000). However, Proposition 6 more explicitly shows this is the unique degenerate configuration of the set of CCEs (or, SCEs). Since an exact factorization, H = P [a]× P −1 , is not possible in practice due to the data error, we need a numerical algorithm to determine the IAC ω.
Second, though the rotation R can be computed directly with the infinite homography H∞ once the intrinsic matrix K is known, we recommend that the rotation R be computed with the formula R = ϕ(H K ). Since the intrinsic matrix K here is obtained from the SCE, R should be computed by R = ϕ(H K ). The combination of such an R and K can best match the image data, and thus can give a better computation of the rotation.
4.3 Standard Algorithm
4.4 Geometrical Interpretation of H
Goal Given N views of a scene and assume the infinite homographies Hi∞ between the 0th and ith views are known. Compute the intrinsic matrix K and the rotation Ri between the 0th and ith views. Algorithm The algorithm consists of the following three steps: 1. Let H˜ i∞ = Det−1/3 (Hi∞ )Hi∞ and compute its Cayley transformation Hi = ϕ(H˜ i∞ ). Then find the least squares solution of the set of linear equations on ω: HiT ω + ωHi = 0
(i = 1, 2, . . . , N − 1).
2. Compute the intrinsic matrix K by the Cholesky factorization of ω. 3. Let HiK = K −1 Hi K and compute its Cayley transformation, R˜ i = ϕ(H K ). Then, determine the rotation Ri by finding the best rotation approximation of the matrix R˜ i , i.e., firstly compute the singular value decomposition of the matrix R˜ i , R˜ i = U diag(t1 , t2 , t3 )V T , then set up Ri = U V T . This algorithm is called the standard algorithm in this paper. In the absence of data error and with the infinite precision of computation, the CCE and the SCE should give the same ω. However, in the presence of data error and with numerical computation, their numerical performances manifest a great difference (see Sect. 5.3). Remark 3 First, the matrix H has an important constraint of rank 2. In practice due to the data error, the computed H does not satisfy this constraint. If in the step 1 this constraint is enforced, both success rate and calibration precision of the algorithm can be improved. The following shows an Rank-enforced method: Firstly, compute the eigenvalues factorization of the matrix H , H = U diag(s, t, t)U −1 , then let its real eigenvalue s = 0, and replace H with the matrix = Re(U diag(0, t, t)U −1 ). H Here Re(A) denotes the real part of A. In theory U diag(0, t, t)U −1 is a real matrix, but a very small imaginary part is possible owing to the numerical computation.
We will give a geometrical interpretation of H in this subsection. Let m ↔ m be an image correspondence, both m and m being projections of an infinite point X∞ = (XT , 0)T . Then, αm = KX,
βm = KRX,
and thus, αm − βm = K(I − R)X, αm + βm = K(I + R)X. By the above two equations, we obtain αm − βm = K(I − R)(I + R)−1 K −1 (αm + βm ). Therefore, m − γ m = H (m + γ m )
(γ = β/α).
(29)
From (29), we can see that the matrix H transforms the linear combination m + γ m of an infinite correspondence m ↔ m into another linear combination m − γ m . This indicates that the matrix H is not a transformation between the two views. Equation (29) cannot be used for a linear computation of the matrix H from infinite correspondences due to the unknown scale γ . But a nonlinear solution is possible. By eliminating the unknown scale γ , (29) becomes [m + H m ]× (m − H m) = 0.
(30)
It is a nonlinear constraint on H . If we regard (I − H∞ ) and (I + H∞ ) as two transformation from the first view into the second view, then H = (I + H∞ )−1 (I − H∞ ) is a transformation on the first view, and is also a transformation on the second view because (I + H∞ )−1 is commutative with (I − H∞ ) i.e., H = (I + H∞ )−1 (I − H∞ ) = (I − H∞ )(I + H∞ )−1 .
Int J Comput Vis
(a)
(b)
Fig. 1 The intrinsic geometrical meaning of the transformation H : the vanishing point v ∞ of the rotation axis of R is a 1-dimensional zero space of H and the vanishing line l ∞ of plane orthogonal to the rotation axis of R is 2-dimensional invariant subspace of H . (a) For
an arbitrary x ∈ / l ∞ ∪ {v ∞ }, y = H x ∈ l ∞ is a conjugate point of x on the IAC and H is a transposition on {y, H y}. (b) Except for v ∞ , H transforms all points on l to its conjugate point c(l) on the IAC and H transforms c(l) to the intersection of l and l ∞
By considering H as a transformation on the same view, we can give its intrinsic geometrical meaning as stated in the following proposition.
By the proof of Proposition 2, H has the same eigenvectors with H∞ . Since l ∞ is a 2-dimensional invariant subspace of H∞ , it is also a 2-dimensional invariant subspace of H , and thus {H x, H 2 x, H 3 x} ⊂ l ∞ . By (a), both H 2 x and H 3 x are the conjugate points of H x and H 2 x = H x, therefore there must be H 3 x = H x. Hence, H is a transposition on {H x, H 2 x}. (c) Let x be a fixed point on the line l but not equal to v ∞ , then ∀y ∈ l\{v ∞ } there is
Proposition 7 Let v ∞ be the vanishing point of the rotation axis of R and l ∞ be the vanishing line of plane orthogonal with the rotation axis of R. Assume x is an arbitrary image point, then, (a) If x is not the vanishing point v ∞ , then H transforms x to its conjugate point y on the IAC, i.e., if y = H x, then (x, y) satisfies x T ωy = 0. (b) ∀x ∈ / l ∞ ∪ {v ∞ }, there must be H x ∈ l ∞ and H is a transposition on {H x, H 2 x}, i.e., H (H x) = H 2 x and H (H 2 x) = H x, as shown in Fig. 1(a). (c) Let l be an arbitrary line passing through v ∞ , and c(l) be its conjugate point on the IAC. Then, except for the vanishing point v ∞ , H transforms all points on the line l to the same point c(l) and H transforms c(l) to the intersection of the line l and the line l ∞ , q = l ∞ ∩ l, as shown in Fig. 1(b). Proof (a) From the SCE (20), the matrix ωH is antisymmetric. And thus for any image point x, there must be x T ωH x = 0. Note that m ∈ N (H ) if and only if m = H∞ m, namely, m is the vanishing point v ∞ . Therefore, H x = 0, and thus y = H x is a conjugate point of x on the IAC. (b) ∀x ∈ / l ∞ ∪ {v ∞ }, there must be x = αv ∞ + βx 1 + γ x 2 , where x 1 , x 2 ∈ l ∞ and α, β, γ = 0. Then, we have H x = αH v ∞ + βH x 1 + γ H x 2 = βH x 1 + γ H x 2 ∈ l ∞ .
y = αv ∞ + βx,
β = 0.
And thus, H y = βH x, i.e., H y and H x are the same point. By (a), H x is the conjugate point of the line l on the IAC, c(l). By (b), {c(l), H c(l)} ⊂ l ∞ and H c(l) is a conjugate point of c(l) on the IAC, and thus H c(l) must be the intersection of l and l ∞ . Remark 4 This is an interesting proposition. First, this proposition tells us that using H we can obtain the constraints on the IAC from the vanishing points of two orthogonal directions to each other. Obviously, by the matrix H we cannot obtain the constraints on the IAC from all orthogonal directions, for example, the constraints from two orthogonal directions of the rotation axis of R and its orthogonal line. But, from (b) and (c), it is not difficult to see that H can provide the four independent constraints. Second, this proposition tells also us that H 2 is a 2-dimensional projection with the center v ∞ and the image line l ∞ as shown Fig. 2(a), and H 3 is H self, i.e., H 3 = H . Observing the transformation property of H on two lines passing through v ∞ and a pair of conjugate points on l ∞ , we can have the Euclidean properties of H , as shown in Fig. 2(b). In essence, the constraints on the IAC induced by H are from the orthogonality of 3D structure, while H∞ are from the orthogonality of 3D motion R. Hence, the Cayley transformation, H = ϕ(H∞ ), is a conversion between the
Int J Comput Vis
(a)
(b)
Fig. 2 (a) H 2 is a 2-dimensional projection with the center v ∞ and the image line l ∞ . It maps point x to the intersection of line l ∞ and the line passing through v ∞ and x. (b) Let L(m) and (l) be the directions of the line and the plane in 3D space corresponding to m and l and assume p and q are a pair of conjugate points on the IAC, then
L(v ∞ )⊥ (l ∞ ) and L(p)⊥L(q). Since x is on line v ∞ ×p, L(x) is on the plane L(v ∞ ) × L(p), similarly L(y) is on the plane L(v ∞ ) × L(q). Proposition 7 shows that H transforms L(x) to its orthogonal direction L(q), L(y) to its orthogonal direction L(p), and transposes the two orthogonal directions L(p) and L(q) to each other
3D motion orthogonality and the 3D structure orthogonality. Third, this proposition gives a linear constraint on H as follows. If a line L and a plane P are orthogonal to each other in 3D space, then their vanishing point x and vanishing line l satisfy l T H x = 0. The only exception is the case of rotation axis of R and its orthogonal plane. A proof is ready: Assume the line L is not parallel to the rotation axis of R, then the relation of its vanishing point x ∈ / N (H ) and the vanishing line of its orthogonal plane is l = ωx. By (a), x T ω(H x) = 0, i.e., l T H x = 0 which is a linear constraint on H .
Equation (31) is called the Cayley transformation family of H . Next, we show some properties of the Cayley transformation family. Proposition 8 (a) ϕ(1 · H ) = H∞ ,
ϕ(0 · H ) = I.
(32)
(b) For τ ∈ (−∞, +∞), ϕ −1 (τ H ) = ϕ(−τ H )
or,
ϕ(τ H ) · ϕ(−τ H ) = I, τ H = ϕ 2 (τ H ).
(33) (34)
5 Cayley Calibration Family
Proof We only prove (33) here. Both (I − τ H ) and (I + τ H ) are invertible for τ ∈ (−∞, +∞), therefore
5.1 Cayley Transformation Family of H
ϕ −1 (τ H ) = (I + τ H )(I − τ H )−1
Since H = (I − H∞ )(I + H∞ )−1 has the eigenvalues:
= (I − (−τ H ))(I + (−τ H ))−1 = ϕ(−τ H ).
{0, i tan(θ/2), −i tan(θ/2)}, for τ ∈ (−∞, +∞) the eigenvalues of matrix I ± τ H must be {1, 1 ± iτ tan(θ/2), 1 ∓ iτ tan(θ/2)}.
The matrix ϕ(τ H ) can be computed from the infinite homography H∞ by (35) via the following proposition. Proposition 9 ϕ(τ H ) = ((1 + τ )I + (1 − τ )H∞ )−1 × ((1 − τ )I + (1 + τ )H∞ ).
Then, det(1 ± τ H ) = 1 + τ 2 tan2 (θ/2) = 0.
Proof
Hence (τ H ) has the Cayley transformation:
ϕ(τ H ) = (I + τ H )−1 (I − τ H ) = (I + τ (I + H∞ )−1 (I − H∞ ))−1
ϕ(τ H ) = (I − τ H )(I + τ H )−1 for τ ∈ (−∞, +∞).
(31)
× (I − τ (I + H∞ )−1 (I − H∞ ))
(35)
Int J Comput Vis
= ((I + H∞ ) + τ (I − H∞ ))−1 (I + H∞ )
equations are equivalent to the CCE-I. Evidently, the CCF-I has the following equivalent form:
× (I + H∞ )−1 ((I + H∞ ) − τ (I − H∞ ))
ϕ T (τ H )ω = ωϕ −1 (τ H )
= ((1 + τ )I + (1 − τ )H∞ )−1 × ((1 − τ )I + (1 + τ )H∞ ).
Let ϕK (τ H ) = K −1 ϕ(τ H )K. If the intrinsic matrix K is known, then we can compute the rotation R by (36) via the following proposition. Proposition 10 R = ((τ + 1)I + (τ − 1)ϕK (τ H ))−1 × ((τ − 1)I + (τ + 1)ϕK (τ H )).
(36)
Proof From (35), we have ϕK (τ H ) = K −1 ϕ(τ H )K
× ((1 − τ )I + (1 + τ )R) = ((1 − τ )I + (1 + τ )R) × ((1 + τ )I + (1 − τ )R)
By Proposition 8, ϕ −1 (τ H ) = ϕ(−τ H ), and thus the above equation can be rewritten as ϕ T (τ H )ω = ωϕ(−τ H )
(τ = 0).
(38)
The advantage of this form is to avoid the matrix inverse operations, hence improve numerical stability. Equation (38) is called Cayley calibration family of the second type, denoted as CCF-II(τ ). Next, we give an explanation of the CCF in term of the infinite homography. In other words, we will prove that for τ = 0, the Cayley transformation ϕ(τ H ) is an infinite homography. Proposition 12 For τ = 0, the Cayley transformation ϕ(τ H ) of τ H is an infinite homography and the corresponding rotation has the same rotation axis with the rotation between the two views.
= ((1 + τ )I + (1 − τ )R)−1
−1
(τ = 0).
Proof Let R be the rotation between the two given views. Then according to Proposition 9, we have
.
ϕ(τ H )
Hence,
= K ((1 + τ )I + (1 − τ )R)−1 ((1 − τ )I + (1 + τ )R) K −1 .
ϕK (τ H )((1 + τ )I + (1 − τ )R)
U (τ )
= ((1 − τ )I + (1 + τ )R). By some simple manipulations, we can obtain (36).
(39)
5.2 Cayley Calibration Family Since for τ ∈ (−∞, +∞), ϕ(τ H ) is the Cayley transformation of τ H , according to Proposition 1 the equation ϕ T (τ H )ωϕ(τ H ) = ω is equivalent to the equation (τ H )T ω + ω(τ H ) = 0. When τ = 0, both the two equations become the identity equation, and thus we do not consider this trivial case. Evidently, for τ = 0 the SCE is equivalent to the equation (τ H )T ω + ω(τ H ) = 0. Thus, we have the following proposition.
In order to prove Proposition 12, we need only to prove that U (τ ) is a rotation and it has the same rotation axis with R. For this purpose, let the eigenvalues factorization of R be: R = V diag(1, eiθ , e−iθ )V −1 .
(40)
Proposition 11 The SCE is equivalent to the following equation:
Then θ is the rotation angle of R, V = (v 1 , v 2 , v 3 ) is a unitary matrix and v 1 is the rotation axis of R. Substituting (40) into U (τ ) in (39) and by some simple manipulations, we can obtain 1 − τ + (1 + τ )eiθ U (τ ) = V diag 1, , 1 + τ + (1 − τ )eiθ 1 − τ + (1 + τ )e−iθ V −1 . 1 + τ + (1 − τ )e−iθ
ϕ T (τ H )ωϕ(τ H ) = ω
It is not difficult to prove that
(τ = 0).
(37)
The above equation gives a one-parameter family equivalent to the SCE. Since it has the same form with the CCE-I, we call it Cayley calibration family of the first type, denoted as CCF-I(τ ). In the CCF-I(τ ), the equation corresponding to τ = 1is just the CCE-I, i.e. CCF-I(1) = CCE-I, and all the
1 − τ + (1 + τ )eiθ 1 + τ + (1 − τ )eiθ 2 cos(θ/2) + iτ sin(θ/2) = = z(τ ), cos2 (θ/2) + τ 2 sin2 (θ/2)
Int J Comput Vis
1 − τ + (1 + τ )e−iθ 1 + τ + (1 − τ )e−iθ 2 cos(θ/2) − iτ sin(θ/2) = = z(τ ). cos2 (θ/2) + τ 2 sin2 (θ/2) Thus, U (τ ) can be rewritten as U (τ ) = V diag(1, z(τ ), z(τ ))V T .
(41)
Since z(τ )is a unit complex number, U (τ ) is a rotation and v 1 is its rotation axis. Remark 5 From the proof of Proposition 12, we can see that the rotation angle of rotation U (τ ) corresponding to infinite homography ϕ(τ H ) is θ (τ ) = 2 Arctan(τ tan(θ/2)).
(42)
Thus if τ varies in (−∞, +∞), we can obtain all the rotations with the fixed rotation axis v 1 . In other words, given an infinite homography H∞ with the rotation R (rotation axis v and rotation angle θ ), by the Cayley transformation of τ ϕ(H∞ ) or, τ H we can obtain all the infinite homographies with the fixed rotation axis v but arbitrary rotation angle, and thus Proposition 6 tells us again that the configuration of cameras with the same rotation axis is a degenerate configuration. Of course, by changing the rotation angle from the eigenvalues factorization of H∞ , we can also obtain these infinite homographies, but it needs to compute the eigensystem of H∞ . In the contrast, we here need merely taking two times of the Cayley transformation, ϕ(τ ϕ(H∞ )), in order to change the rotation angle. 5.3 Improving Calibration Conditioning Given the infinite homography H∞ between two views, we can obtain the CCF-I(τ ) equivalent to the CCE-I: CCF-I(τ ) :
ϕ T (τ H )ωϕ(τ H ) − ω = 0,
τ = 0,
and CCE-I = CCF-I(1) (ϕ(1 · H ) = H∞ ). In the CCF-I(τ ), if we define CCF-I(0) by H T ω + ωH = 0 (SCE), then by Proposition 10 we obtain an equivalent family with τ in (−∞, +∞), denoted as CCF-I also. Similarly, for the CCF-II(τ ), we define CCF-II(0) also by the SCE, and thus obtain an equivalent family of the second type with τ in (−∞, +∞). The numerical stability of a homogeneous linear system, such as Ax = 0, can be assessed by examining the condition number of its coefficient matrix A, it is also called the condition number of this homogeneous linear system. The condition number is a typical measure of the sensitivity of homogeneous linear system to noise. A large condition number
usually means less numerical stability, and a small condition number indicates good numerical stability. If the rank of the matrix A is k, then its condition number is defined as the ratio of the largest singular value to the kth smallest singular value (Hartley 1997c). And thus the condition number of the matrix A is computed by taking the singular value decomposition of this matrix. Figure 3 gives the condition numbers of equations in 5 CCF-Is and 5 CCF-IIs. From this figure, we can see that different equation in the CCF-I (or CCF-II) has different condition number, and thus has different numerical stability. Hence, we can improve the numerical stability of calibration by selecting a parameter τ such that corresponding equation has a smaller condition number. We next give an algorithm for selecting the best parameter value. Let the set of m CCF-Is (or CCF-IIs) be ϕ T (τi Hi )ωϕ(τi Hi ) − ω = 0,
i = 1, 2, . . . , m
(or, ϕ T (τi Hi )ω − ωϕ(−τi Hi ) = 0, i = 1, 2, . . . , m). The m-parameter family of the above homogeneous linear systems can be rewritten as the following matrix form: τ = (τ1 , τ2 , . . . , τm ) ∈ (−∞, +∞)×m ,
A(τ )c = 0,
(43)
here c = (ω11 , ω12 , ω13 , ω22 , ω23 , ω33 )T . In theory, the rank of matrix A is 5, therefore our problem is to solve the following minimization: σ1 ×m min J = (τ )|τ ∈ (−∞, +∞) , (44) σ5 where σ1 ≥ σ2 ≥ · · · ≥ σ5 ≥ σ6 are singular values of A(τ ). For solving the above problem, we firstly compute the gradient of J . By taking the singular value decomposition of A, A = U DV T , then σj = uTj Av j where uj , v j are the j th row of the orthogonal matrices U and V , respectively. And thus we have ∂(uTj Av j ) ∂σj = ∂τi ∂τi =
∂uTj ∂τi
= σj Notice have
Av j + uTj
∂uTj ∂τi
∂v Tj ∂τi
∂A ∂v n v j + uTj A ∂τi ∂τi
uj + uTj
v j = uTj
∂uj ∂τi
∂v j ∂A v j + σj v Tj . ∂τi ∂τi = 0 from uTj uj = v Tj v j = 1, we
∂σj ∂A = uTj vj . ∂τi ∂τi Hence the gradient of J on τi is: ∂A ∂A ∂J 1 = 2 σ5 uT1 v 1 − σ1 uT5 v5 . ∂τi ∂τi ∂τi σ5
(45)
Int J Comput Vis
(a) Condition numbers of CCF-Is
(b) Condition numbers of CCF-IIs
Fig. 3 Condition numbers of the CCF-Is and CCF-IIs corresponding to 5 infinite homographies H∞ = KRK −1 where K : fu = fv = 1000, s = 0.0, pu = pv = 512 and the rotation Ri = Rx (θi1 ) · Ry (θi2 ) · Rz (π/θi3 ) (i = 1, 2, 3, 4, 5), θij is randomly generated by the uniform distribution on the interval [−π/4, π/4]. The abscissa and the ordinate are parameter τ and condition number respectively. We can see that: For the CCF-I, the condition numbers of equations vary with the parameter τ and the condition number reaches its minimum at τ = 0 i.e.,
the SCE has the smallest condition number among the CCF-I. For the CCF-II, the condition numbers of equations vary with the parameter τ also, but the condition numbers reach their minimums at different τ . The condition numbers of the CCF-I and the CCF-II for the same parameter are different and the CCF-II has always a smaller condition number, thus the numerical performance of CCF-II is better than that of CCF-I in general
With the computed gradient ∇J (τ ) of J , we can use some numerical methods to solve the minimization problem (44), such as steepest descent method or conjugate gradient method (Press et al. 1988). In the experiments of this paper, we use the steepest descent method as outlined below.
1. Let s0 = 0, s1 = 1, n = 1 and given a termination δ. 2. Compute e = J (τ k + sn−1 p k ) − J (τ k + sn p k ). 3. In the case e > 0, set up sn+1 = sn + 10λ and let n ← n + 1, return to Step 2. In the case e < 0, if |sn+1 − sn | < δ, terminate iteration and put out sn , or else do Step 4. 4. (a) t = λ/10. (b) Set up sn+1 = sn + t and do Step 2. (c) In the case e > 0, let n ← n + 1 and return to (b). In the case e < 0, if |sn+1 − sn | < δ, terminate iteration and put out sn , or else let t ← t/10 and return to (b).
Determining the Parameters τ The algorithm consists of the following four steps: 1. Given an initialization τ 0 and a termination ε. Let k = 0. 2. Compute the gradgent ∇J (τ k ) of J with (45), ∂J ∂J ∂J , ,..., . ∇J (τ k ) = ∂τ1 ∂τ2 ∂τm τ =τ k If ∇J (τ k ) < ε, terminate iteration and put out τ k , or else do Step 3. 3. Let p k = −∇J (τ k ). Solve sk such that J (τ k + sk pk ) = min J (τ k + sp k )|s ≥ 0 . 4. Let τ k+1 = τ k + sk p k and k ← k + 1. Return to Step 2. Remark 6 First, in this paper the initialization is set to τ 0 = (1, 1, . . . , 1), corresponding to the set of CCEs. Second, in Step 3, the exact solution s cannot be obtained since J (τ k + spk ) has no analytical expression. We can use an algorithm similar to Levenberg-Marquardt algorithm (Press et al. 1988) to find an approximate solution of sk as follows:
In our work, δ = 10−5 and λ = 0.1 are used. Third, in practice finding an exact solution of the minimization (44) is unnecessary. This is because the condition number is a qualitative analysis of numerical stability; the condition numbers with the same order of magnitude are unable to lead to a big difference in numerical stability. And thus, we can terminate iteration once the condition number is less than a given threshold. In this work, 50 is used as the threshold. 5.4 Cayley Algorithm By the discussions in Sects. 5.1, 5.2 and 5.3, we can also give a calibration algorithm similar to the standard one, called the Cayley algorithm in this paper. Goal Given N views of a scene and assume the infinite homographies Hi∞ between the 0th and ith views are known. Compute the intrinsic matrix K and the rotation Ri between the 0th and ith views.
Int J Comput Vis
Algorithm The algorithm consists of the following four steps: 1. Let H˜ i∞ = det−1/3 (Hi∞ )Hi∞ and compute its Cayley transformation Hi = ϕ(H˜ i∞ ) and CCF-I (or CCF-II): ϕ T (τi Hi )ωϕ(τi Hi ) = ω,
i = 1, 2, . . . , N
(or ϕ T (τi Hi )ω = ωϕ(−τi Hi ), i = 1, 2, . . . , N ). 2. Determine the parameters τ = (τ1∗ , τ2∗ , . . . , τN∗ ) by solving the minimization (44), then find the least squares solution of the set of linear equations on ω: ϕ T (τi∗ Hi )ωϕ(τi∗ Hi ) = ω,
i = 1, 2, . . . , N
(or ϕ T (τi∗ Hi )ω = ωϕ(−τi∗ Hi ), i = 1, 2, . . . , N ). 3. Determine the intrinsic matrix K by taking the Cholesky factorization of ω. 4. Let ϕK (τi∗ H ) = K −1 ϕ(τi∗ Hi )K and compute
× ((τi∗ − 1)I + (τi∗ + 1)ϕK (τi∗ Hi )).
6.1 Generalized Cayley Transformation If a function of complex variable, ξ = φ(λ), satisfies the following three conditions, then it is called a generalized Cayley transformation: (a) φ(λ) is an analytical function with real coefficients in |λ| < δ; (b) φ(λ)φ(−λ) = 1 for |λ| < δ; (c) φ(0) > 0 and φ (0) = 0. It is easy to see that the Cayley transformation is a special case of the generalized Cayley transformation. Below are two examples of generalized Cayley transformations:
(49)
Then for −∞ < τ < ∞ the transformation (48) is a generalized Cayley transformation of matrix A. If r(A) ≤ 1 (here r(A) = max{|λ| : λ ∈ σ (A)} is called the spectrum radius of matrix A), then for −π/2 < τ < π/2, the transformation (49) is also a generalized Cayley transformation of matrix A. For generalized Cayley transformations, we have the following proposition. Proposition 13 Assume φ(τ A) (−ε < τ < ε) is a generalized Cayley transformation of matrix A ∈ C n×n and C ∈ C n×n is an arbitrary matrix. Then φ(τ A)T Cφ(τ A) = C
for 0 < |τ | < ε
AT C + CA = 0.
(50)
∂ ∂ g(τ A) = φ(τ A)T Cφ(τ A) ∂τ ∂τ ∂ φ(τ AT )Cφ(τ A) = ∂τ
Assume equation φ(τ A)T Cφ(τ A) = C holds for 0 < |τ | < ε, then for 0 < |τ | < ε there must exist AT φτ (τ AT )Cφ(τ A) + φ(τ AT )CAφτ (τ A) = 0.
(52)
Because φ(λ) is analytical in −ε < τ < ε, the left side in (52) is a continuous mapping on the variable τ . Thus, letting τ → 0 in (52) leads to φ (0)φ(0)(AT C + CA) = 0. By the conditions (b) and (c), φ (0)φ(0) = 0, and thus there must be AT C + CA = 0.
1 − sin(λ) φs (λ) = . 1 + sin(λ)
For a matrix A ∈ C n×n , a generalized Cayley transformation of A can be defined as (46)
Now assume equation AT C + CA = 0 holds, then we can verify that φ(τ AT )C = Cφ(−τ A),
φτ (τ AT )C = Cφτ (−τ A)
holds for any analytical function φ(λ). Then we obtain AT φτ (τ AT )Cφ(τ A) = AT Cφτ (−τ A)φ(τ A),
where ε is a positive number such that σ (εA) ⊂ {λ : |λ| < δ}.
φs (τ A) = (I − sin(τ A))(I + sin(τ A))−1 .
= AT φτ (τ AT )Cφ(τ A) + φ(τ AT )CAφτ (τ A). (51)
6 Generalized Cayley Calibration Families
(−ε < τ < ε),
(48)
Proof Let g(τ A) = φ(τ A)T Cφ(τ A), then,
Then, determine the rotation Ri by finding the best rotation approximation of the matrix R˜ i , as shown in the standard algorithm in Sect. 4.3.
Bτ = φ(τ A)
φe (τ A) = (I + e−τ A )(I + eτ A )−1 ,
⇔
Ri = ((τi∗ + 1)I + (τi∗ − 1)ϕK (τi∗ Hi ))−1
1 + e−λ , φe (λ) = 1 + eλ
For example, let
(47)
φ(τ AT )CAφτ (τ A) = CAφ(−τ A)φτ (τ A).
Int J Comput Vis
So, by (51), we have ∂ g(τ A) = AT Cφτ (−τ A)φ(τ A) + CAφ(−τ A)φτ (τ A). ∂τ
Proof Let R be the rotation between the two given views and its eigenvalues factorization be:
By the condition (b), there is φτ (τ A)φ(−τ A) = φ(τ A)φτ (−τ A).
R = V diag(1, eiθ , e−iθ )V −1 ,
And thus, for 0 < |τ | < ε,
where θ is its rotation angle, V = (v 1 , v 2 , v 3 ) is an unitary matrix and v 1 is the rotation axis of R. Then according to Proposition 2,
∂ g(τ A) = (AT C + CA)φ(τ A)φτ (−τ A) = 0. ∂τ Therefore, for 0 < |τ | < ε,
τ H = Kτ ϕ(R)K −1
g(τ A) = φ(τ A)T Cφ(τ A) = constant matrix.
= KVdiag(0, iτ tan(θ/2), −iτ tan(θ/2))V −1 K −1 .
By the continuity of φ(τ A) at τ = 0, we have
Then,
φ(τ A)T Cφ(τ A) = g(τ A) = φ 2 (0)C
φ(τ H ) = K V diag(φ(0), φ(is(τ )), φ(−is(τ )))V −1 K −1 ,
for 0 < |τ | < ε.
Q(τ )
By the condition (b), φ 2 (0) = 1, and thus there must be φ(τ A)T Cφ(τ A) = C.
6.2 Generalized Cayley Calibration Family Applying a generalized Cayley transformation to the Cayley transformation of the infinite homography, we can obtain the following proposition according to the Proposition 13. Proposition 14 The SCE is equivalent to the following oneparameter family φ(τ H )T ωφ(τ H ) = ω
Proposition 15 The generalized Cayley transformation φ(τ H ) is an infinite homography and its rotation has the same rotation axis with the rotation between the two views.
for 0 < |τ | < ε.
(53a)
Given a generalized Cayley transformation, by (53a) we obtain a one-parameter family equivalent to the SCE. Since it has the same form with the CCE-I and CCF-I, it is called the generalized Cayley calibration family of the first type (GCCF-I). Because φ −1 (τ H ) = φ(−τ H ) for −ε < τ < ε, the SCCF-I has the following equivalent form, it is called the generalized Cayley calibration family of the second type (GCCF-II): φ(τ H )T ω = ωφ(−τ H ) for 0 < |τ | < ε.
(53b)
As in the CCFs, the numerical stability of equations in a GCCF varies also with the parameter τ . For a given generalized Cayley calibration family, we can use the method in Sect. 5.3 to reduce the condition number of calibration equations by selecting a proper parameter value and can give a calibration algorithm similar to the algorithm in Sect. 5.3, called as generalized Cayley algorithm. Next, we give an explanation of (53a) also in terms of the infinite homography.
here s(τ ) = τ tan(θ/2). In order to prove Proposition 15, we need only to prove that Q(τ ) is a rotation and it has the same rotation axis with R. From the conditions (b) and (c) of the function φ, there is φ(0) = 1. By the conditions (b) of the function φ, we have φ(is(τ )) · φ(−is(τ )) = 1. Then from the conditions (a) of the function φ, we have φ(−is(τ )) = φ(is(τ )). Therefore {φ(is(τ )), φ(−is(τ ))} is a conjugate pair of unit complex numbers, and Q(τ ) can be written as: Q(τ ) = V diag(1, eiσ (τ ) , e−iσ (τ ) )V −1 , where σ (τ ) =
1 1 log(φ(is(τ ))) = log(φ(iτ tan(θ/2))). i i
Hence, Q(τ ) is a rotation and v 1 is its rotation axis.
Remark 7 The GCCF-I contains also the CCE-I. Let θ be the rotation angle between the two given views, then in the GCCF-I the parameter τ corresponding to the CCE-I is the solution of the following equation: log(φ(iτ tan(θ/2))) = iθ. 6.3 S-Cayley Calibration Family As an example, we examine the generalized Cayley calibration family induced by the transformation (49): B T (τ H )ωB(τ H ) = ω
for 0 < |τ |
