张正友相机标定--A Flexible New Technique for Camera Calibration

【实例简介】
张正友相机标定方法，经典论文，英文原文----A Flexible New Technique for Camera Calibration
A Flexible New Technique for Camera Calibration Abstract We propose a flexible new technique to easily calibrate a camera. It is well suited for use without specialized knowledge of 3D geometry or computer vision. The technique only requires the camera to observe a planar pattern shown at a few(at least two)different orientations. Either the camera or the planar pattern can be freely moved. The motion need not be known. Radial lens distortion is modeled. The proposed procedure consists of a closed-form solution, followed by a non linear refinement based on the maximum likelihood criterion both computer simulation and real data have been used to test the proposed technique and very good results have been obtained Compared with classical techniques which use expensive equipment such as two or three orthog- onal planes the proposed technique is easy to use and flexible. It advances 3d computer vision one step from laboratory environments to real world use Index Terms--Camera calibration, calibration from planes, 2D pattern, absolute conic, projective mapping, lens distortion, closed-form solution, maximum likelihood estimation, flexible setup 1 Motivations Camera calibration is a necessary step in 3D computer vision in order to extract metric information from 2D images. Much work has been done, starting in the photogrammetry community(see [2 4] to cite a few), and more recently in computer vision( [9, 8, 23, 7, 26, 24, 17, 6] to cite a few We can classify those techniques roughly into two categories: photogrammetric calibration and self calibration Photogrammetric calibration. Camera calibration is performed by observing a calibration object whose geometry in 3-D space is known with very good precision. Calibration can be done very efficiently [5]. The calibration object usually consists of two or three planes orthogonal to each other. Sometimes, a plane undergoing a precisely known translation is also used [23]. These approaches require an expensive calibration apparatus, and an elaborate setup Selfi-calibration. Techniques in this category do not use any calibration object. Just by moving a camera in a static scene, the rigidity of the scene provides in general two constraints [17, 15 on the cameras'internal paramcters from one camera displacement by using image informa- tion alone. Therefore, if images are taken by the same camera with fixed internal parameters correspondences between three images are sufficient to recover both the internal and external parameters which allow us to reconstruct 3-D structure up to a similarity [16, 13]. while this ap proach is very flexible, it is not yet mature [1]. Because there are many parameters to estimate we cannot always obtain reliable results Other techniques exist: vanishing points for orthogonal directions [3, 14], and calibration from pure rotation [11, 21] Our current research is focused on a desktop vision system (Dvs)since the potential for using DVSs is large. Cameras are becoming cheap and ubiquitous. A DVs aims at the general public who are not experts in computer vision. a typical computer user will perform vision tasks only from time to time, so will not be willing to invest money for expensive equipment. Therefore, flexibility, robustness and low cost are important. The camera calibration technique described in this paper was developed with these considerations in mind The proposed technique only requires the camera to observe a planar pattern shown at a few(at least two)different orientations. The pattern can be printed on a laser printer and attached to a"rea- sonable planar surface(e.g, a hard book cover). Either the camera or the planar pattern can be moved by hand. The motion need not be known. The proposed approach lies between the photogrammet- ric calibration and self-calibration, because we use 2D metric information rather than 3D or purely Imp plicit one. Both computer simulation and real data have been used to test the proposed technique and very good results have been obtained. Compared with classical techniques, the proposed tech- nique is considerably more flexible. Compared with self-calibration, it gains considerable degree of robustness. We believe the new technique advances 3D computer vision one step from laboratory environments to the real world Note that Bill Triggs [22 recently developed a self-calibration technique from at least 5 views of a planar scene. His technique is more flexible than ours, but has difficulty to initialize. Liebowitz and Zisserman [14] described a technique of metric rectification for perspective images of planes using metric information such as a known angle, two equal though unknown angles, and a known length ratio. They also mentioned that calibration of the internal camera parameters is possible provided at least three such rectified planes, although no experimental results were shown The paper is organized as follows. Section 2 describes the basic constraints from observing a single planc. Section 3 describes the calibration procedure. We start with a closed-form solution followed by nonlinear optimization. Radial lens distortion is also modeled. Section 4 studies con figurations in which the proposed calibration technique fails. It is very easy to avoid such situations in practice. Section 5 provides the experimental results. Both computer simulation and real data are used to validate the proposed technique. In the Appendix, we provides a number of details, including the techniques for estimating the homography between the model plane and its image. 2 Basic Equations We examine the constraints on the cameras intrinsic parameters provided by observing a single plane We start with the notation used in this paper 2.1 Notation A 2D point is denoted by m-u, U. A 3D point is denoted by M-X,Y, Z. We use x to denote the augmented vector by adding I as the last element: m-u, u, 1]"andM-[X,Y, Z, 1] .A camera is modeled by the usual pinhole: the relationship between a 3d point m and its image projection m is sm=AR tM, (1) where s is an arbitrary scale factor,(R, t), called the extrinsic parameters, is the rotation and trans lation which relates the world coordinate system to the camera coordinate system, and A, called the camera intrinsic matrix, is given b A=0B0 with(uo, vo) the coordinates of the principal point, a and B the scale factors in image u and v axes and y the parameter describing the skewness of the two image axes We use the abbreviation A-for(A-l)or(A)-1 2.2 Homography between the model plane and its image Without loss of generality, we assume the model plane is on Z=0 of the world coordinate system Let s denote the i column of the rotation matrix R by ri. From(1), we have r2 rat A r By abuse of notation, we still use M to denote a point on the model plane, but M=X, Y since Z is always equal to O. In turn, M=[X, Y, 1. Therefore, a model point M and its image m is related by a homography H sm=HM with H=A rI r2 t As is clear, the 3 x 3 matrix H is defined up to a scale factor 2.3 Constraints on the intrinsic parameters Given an image of the model plane, an homography can be estimated(see Appendix a). Let's denote it by H=h1 h2 h3. From(2), we have hI h2 h=XA rl r2t where a is an arbitrary scalar. Using the knowledge that ri and r2 are orthonormal, we have hiRAh=0 hiaAh=haha These are the two basic constraints on the intrinsic parameters, given one homography. Because a homography has 8 degrees of freedom and there are 6 extrinsic parameters (3 for rotation and 3 for translation), we can only obtain 2 constraints on the intrinsic parameters. Note that a-TA-l actually describes the image of the absolute conic [16]. In the next subsection, we will give an geometric interpretation 2.4 Geometric Interpretation We are now relating (3)and (4 )to the absolute conic It is not difficult to verify that the model plane, under our convention, is described in the camera coordinate system by the following equation: where w=o for points at infinity and w= l otherwise. This plane intersects the plane at infinity at a line, and we can casily scc that r1 and/2 "e two particular points on that line. Any point on it is a linear combination of these two points, i.e ar1+br 0 0 pbi: Now, let,'s compute the intersection of the above line with the absolute conic. By definition, the Int xoo, known as the circularpoint, satisfies,r> =0,1 r1+br2)2(ar1+br2)=0,ora2+b2=0 The solution is b-+ ai, where i2--1. That is, the two intersection points are r1土2 C Their projection in the image plane is then given, up to a scale factor, by A(r1±ir2)=h1士 Point moo is on the image of the absolute conic, described by A-TA-[16]. This gives (h±ih2)2A-A-1(h1±ih2)=0 Requiring that both real and imaginary parts be zero yields (3)and(4) 3 Solving Camera Calibration This section provides the details how to effectively solve the camera calibration problem. We start with an analytical solution, followed by a nonlinear optimization technique based on the maximum likelihood criterion. Finally, we take into account lens distortion, giving both analytical and nonlinear solutions 3.1 Closed-form solution B11B1231 B=AA B12 B22 B B13 B23 Bs r(U U0 +2+1 Note that B is symmetric, defined by a 6D vector b=[B1,B12,B2,B13,B23,B3 (6) Let the itn column vector of H be hi=[hal, hi2, hi3].Then, we have h Bh Vi=[hihg,hi1hj2+hi2hgl, hi2hi 2 Pishi1+ hi1hj 3, hi3hj2+hi] h 3, h:3h;31/ Therefore, the two fundamental constraints (3)and (4), from a given homography, can be rewritten as 2 homogeneous equations in b N12 T (v11-V22 If n images of the model plane are observed, by stacking n such equations as( 8 )we have Vb=0, (9 where V is a 2n x 6 matrix. If n> 3, we will have in general a unique solution b defined up to a scale factor. If n =2, we can impose the skewless constraint y =0,i.e,0,1,0,0.0,0b=0, which is added as an additional equation to(9).(If n=l, we can only solve two camera intrinsic parameters e.g., a and B, assuming uo and vo are known(e. g, at the image center) and y=0, and that is indee what we did in [19] for head pose determination bascd on the fact that eyes and mouth are reasonably coplanar. The solution to(9) is well known as the eigenvector of vv associated with the smallest eigenvalue(equivalently, the right singular vector of V associated with the smallest singular value) Once b is estimated, we can compute all camera intrinsic matrix A. Sce Appendix b for the details Once a is known, the extrinsic parameters for each image is readily computed. From(2), we have 1=XA hi r2=入A-h 3=r1×r: t=AA h3 with A-1/A h1ll-1/A h2l. Of course, because of noise in data, the so-computed matrix R.=r1, r2, ra does not in gencral satisfy the properties of a rotation matrix. Appendix C describes a method to estimate the best rotation matrix from a general 3 x 3 matrix 3.2 Maximum likelihood estimation The above solution is obtained through minimizing an algebraic distance which is not physicall meaningful. We can refine it through maximum likelihood inference We are given n images of a model plane and there are m points on the model plane. Assume that the image points are corrupted by independent and identically distributed noise. The maximum likelihood estimate can be obtained by minimizing the following functional ∑∑m-m(AR,t,N川2 (10) where In(a, Ri, ti, Mi)is the projection of point M; in image i, according to equation(2). A rotation R is parameterized by a vector of 3 parameters, denoted by r, which is parallel to the rotation axis and whose magnitude is equal to the rotation angle. R and r are related by the rodrigues formula [5] Minimizing (10)is a nonlinear minimization problem, which is solved with the Levenberg- Marquardt Algorithm as implemented in Minpack [18]. It requires an initial guess of A, Ri, tili=1.ny which can be obtained using the technique described in the previous subsection 3.3 Dealing with radial distortion Up to now, we have not considered lens distortion of a camera. However, a desktop camera usually exhibits significant lens distortion, especially radial distortion. In this section, we only consider the first two terms of radial distortion. The reader is referred to [20, 2, 4, 26] for more elaborated models Based on the reports in the literature [2, 23, 25], it is likely that the distortion function is totall dominated by the radial components, and especially dominated by the first term. It has also been found that any more elaborated modeling not only would not help(negligible when compared with sensor quantization), but also would cause numerical instability [23, 25 t(u, u) be the ideal(nonobservable distortion-free) pixel image coordinates, and(u, i) the corresponding real observed image coordinates. The ideal points are the projection of the model points according to the pinhole model. Similarly, (., g) and(i, y)are the ideal(distortion-free)and real(distorted) normalized image coordinates. We have [2, 25 =x+xk1(x2+y2)k2(x2十 y=y+y(k1(x2+y2)+k2(x2+y2)2], where k1 and kg are the coefficients of the radial distortion. The center of the radial distortion is the same as the principal point Fromu=uo + a i t ry and v=vo I By and assuming y =0, we have +(-0)[k1(x2+y2)+k2(x2-y2)2 =0+(v-)k1(x2+y2)+k2(x2+y2)2] (12) Estimating Radial Distortion by Alternation. As the radial distortion is expected to be small,one would expect to estimate the other five intrinsic parameters, using the technique described in Sect. 3.2 reasonable well by simply ignoring distortion. One strategy is then to estimate hi and k2 after having and(12), we have two equations for each point in each imagc el coordinates(u, u). Then, from(11) ,-0 CHAPI kz Given m points in n images, we can stack all equations together to obtain in total 2mn equations, or in matrix form as Dk=d, where k=k1, k2. The linear least-squares solution is given by k=(DD)Dd (13) Once hi and k2 are estimated, one can refine the estimate of the other parameters by solving(10)with m(A,Ri, ti, Mi)replaced by(11)and(12). We can alternate these two procedures until convergence Complete Maximum Likelihood Estimation. Experimentally, we found the convergence of the above alternation technique is slow. A natural extension to (10)is then to estimate the complete set of parameters by minimizing the following functional ∑∑‖m;-血(A,k,k,R,t,M川2 (14) a typo was reported by Johannes Koester johannes koester(@ uni-dortmund de] via email on Aug. 13, 2008 where Im(A,k:1, k 2. Ri, ti, Mi)is the projection of point M, in image i according to equation(2), followed by distortion according to(ll) and( 12). This is a nonlinear minimization problem, which is solved with the Levenberg-Marquardt Algorithm as implemented in Mirpack [18]. A rotation is again parameterized by a 3-vector r, as in Sect. 3. 2. An initial guess of A and Ri, tii=1.n/can be obtained using the technique described in Sect. 3. 1 or in Sect. 3.2. An initial guess of h1 and k2 can be obtained with the technique described in the last paragraph, or simply by setting them to O 3.4 Summary The recommended calibration procedure is as follows 1. Print a pattern and attach it to a planar surface 2. Take a few images of the model plane under different orientations by moving either the plane or the camera 3. Detect the feature points in the images 4. Estimate the five intrinsic parameters and all the extrinsic parameters using the closed-form solution as described in sect. 3.1 5. Estimate the coefficients of the radial distortion by solving the linear least-squares(13); 6. Refine all parameters by minimizing(14) 4 Degenerate Configurations We study in this section configurations in which additional images do not provide more constraints on the camera intrinsic parameters. Because(3)and(4) are derived from the properties of the rotation matrix,if R2 is not independent of R1, then image 2 does not provide additional constraints. In particular, if a plane undergoes a pure translation, then R2= Ri and image 2 is not helpful for camera calibration. In the following, we consider a more complex configuration Proposition 1. I/the model plane ul the second position is parallel to its first position, then che second homography does not provide additional constraints Proof. Under our convention, R2 and Ri are related by a rotation around z-axis. That is 0- sing 0 RI g cos 9 0 0 0 where a is the angle of the relative rotation. We will use superscript (1)nd(2) to denote vectors related to image I and 2, respectively. It is clear that we have h(2=X(2)(Ar()cos01 Ar(2)sin 0)= 入 A()(ni cos 6+(1) sin e (2)_x(2)(Ar sin g I Ar-co A4( hi sin/6(cos 8) Then, the first constraint(3) from image 2 becomes hiAAh22=1[(cos20-sin20)() cos 0 0(hdA-TA-h(D)-hoDA-TA-Iha) which is a linear combination of the two constraints provided by H1. Similarly, we can show that the second constraint from image 2 is also a linear combination of the two constraints provided by H1 Therefore, we do not gain any constraint from H2 The result is self-evident because parallel planes intersect with the plane at infinity at the same circular points, and thus according to Sect. 2. 4 they provide the same constraints In practice, it is very easy to avoid the degenerate configuration: we only need to change the orientation of the model plane from one snapshot to another Although the proposed technique will not work if the model plane undergoes pure translation camera calibration is still possible if the translation is known. Please refer to Appendix d 5 Experimental Results The proposed algorithm has been tested on both computer simulated data and real data. The closed form solution involves finding a singular value decomposition of a small 21 x 6 matrix, where r is the number of images. The nonlinear refinement within the Levenberg-Marquardt algorithm takes 3 to 5 iterations to converge 5.1 Computer Simulations The simulated camera has the following property: a= 1250, B=900, y=1.09083(equivalent to 89.950), 20= 255, U0= 255. The image resolution is 512x512. The model plane is a checker pattern containing 10x 11=110 corner points (so we usually have more data in the v direction than in the u direction). The size of pattern is 18cmx25cm. The orientation of the plane is represented by a 3d vector r, which is parallel to the rotation axis and whose magnitude is equal to the rotation angle. Its position is represented by a 3D vector t (unit in centimeters) Performance w.r.t. the noise level. In this experiment, we use three planes with rI=[20, 0, 0- -9.-1255002,r2=0.20,0r,t2=(9,-125.510,r3=55-30,-30,-1 tg=[-105, -12.5, 525. Gaussian noise with 0 mean and o standard deviation is added to the projected image points. The estimated camera parameters are then compared with the ground truth We measure the relative error for a and B, and absolute error for uo and vo. We vary the noise level rom OI pixels to 1.5 pixels. For each noise level, we perform 100 independent trials, and the results shown are the average. As we can see from Fig. l, errors increase linearly with the noise level. (The error for y is not shown, but has the same property. For o=0.5(which is larger than the normal noise in practical calibration), the errors in a and b are less than 0.3%, and the errors in uo and vo are around l pixel. The error in uo is larger than that in o. The main reason is that there are less data in the a direction than in the o direction as we said before Performance w.r.t. the number of planes. This experiment investigates the performance with re- spect to the number of planes(more precisely, the number of images of the model plane). The orien- tation and position of the model plane for the first three images are the same as in the last subsection From the fourth image, we first randomly choose a rotation axis in a uniform sphere, then apply a rotation angle of 30, we vary the number of images from 2 to 16. For each number. 100 trials of independent plane orientations(except for the first three) and independent noise with mean o and 【实例截图】
【核心代码】