# An Efficient and Flexible Solution for Camera Autocalibration from N≥3 Views

Authors: Rongfu Tang; Baosong Deng; Jing Li; Ye Yan
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

This paper presents an efficient and flexible solution for camera autocalibration from N≥3 views, given image correspondences and zero (or known) skew only. The knowledge is not required on camera motion, 3D information, scene, or internal constraints. Our method is essentially based only on the fundamental matrices and its main virtues are threefold. Firstly, it is shown that, in the center-orientedmetric coordinates, the focal length and aspect ratio can be estimated independent of considerable principle point shift (PPs). Thus, our method includesrecursive steps: estimating focal length and aspect ratio and then calculating the PPs. Secondly, the optimal geometric constraints are selected for calibration by using error propagation analyses. Thirdly, the Levenberg–Marquardt algorithm is adopted for the fast final refinement of four internal parameters. Our method is fast and efficient to derive a unique calibration. Besides, this method can be applied to calibrate the focal length from two views, without requiring the prior knowledge of PPs. Good performance of our method is evaluated and confirmed in both the simulation experiments and the practical tests.


---

## Body

## 1. Introduction

Camera calibration is an essential topic in photogrammetry and computer vision. Numerous works have been investigated in this subject in the last decades. Besides the image correspondences, many calibration methods need one or more additional constraints, such as camera motion [1–3] (to cite a few), scene constraints by using the vanishing points [4, 5] or plumb-lines [6, 7], 3D or 2D object information [8, 9], or partial calibration information [10, 11].Camera autocalibration is the technique to determine the internal parameters directly from multiple uncalibrated images. It requires only the image correspondences to reconstruct the metric properties of camera and scene, despite the unknown camera motion. Thus, it offers great flexibility. Autocalibration is vital when the camera information is unavailable, such as performing reconstruction by using the historical images or those downloaded from the Internet.Autocalibration was originally introduced by Faugeras et al. [12]. The methods based on Kruppa equation were later developed in [13–15] by using different numerical solutions. Stratification approach was proposed in [16] and the absolute (dual) conic was introduced as a numerical device for formulating autocalibration problem in [17]. These early works are however quite sensitive to noise and unreliable [18, 19]. The numerical solution using interval analysis is presented in [20], but this solution is insufficiently accurate and extremely slow. Recently, a novel autocalibration method was presented by using unknown isosceles right triangle in 3D scene and optimizing a nonlinear cost function in [21]. Autocalibration has also found extensive applications in intelligent transportation systems [22] and video surveillance system [23, 24]. An excellent comprehensive review on the early autocalibration work can be seen in [19]. Other relevant works were studied to estimate focal length from two views, given other internal parameters in [10, 25, 26], to calibrate the structured-light camera [27], to calibrate large-scale camera networks [28, 29]. These techniques are either inaccurate, time consuming, or hard to generalize to multiview. Generally, although it is possible to autocalibrate camera from N≥3 views and various methods have been proposed, it remains quite a difficult problem [19].In this paper, we present an efficient and flexible solution for camera autocalibration from N≥3 views (actually from N≥2 fundamental matrices), given image correspondences and zero (or known) skew only. The skew parameter is insignificant in almost all vision practice [19], and it is considered as zero here as most of the previous methods assumed. The other information is not required on camera motion, object information, scene, or internal constraints. Our method is very fast (around 0.1 second for three views) and efficient to derive a unique calibration from N≥3 views. It is highly flexible as well. It can be applied in two-view to calibrate the focal length, without the prior knowledge on principle point shifts (PPs). Any precise prior information on the internal parameters (known aspect ratio, for instance), if available, can be easily introduced into our solution.Our method is based on fundamental matrix and the two constraints of essential matrix. Its main virtues are threefold, distinctive from previous methods. Firstly, it shows that, in the center-orientedmetric coordinates, the focal length and aspect ratio can be precisely estimated independent of considerable PPs. Thus, in contrast to the conventional methods which reconstruct simultaneously all the unknown internal parameters, our method contains tworecursive steps. It first estimates focal length and aspect ratio and then calculates PPs by fixing the estimation of focal length and aspect ratio. The PPs estimation returns to contribute to refine the estimation of focal length and aspect ratio, and so on so forth. The recursion reduces the numerical complexity and improves the accuracy. Secondly, the technique of error propagation analyses is introduced to choose the optimal geometric constraints. It is shown that the selection of optimal constraints is vital for precise calibration and fast convergence. In fact, using nonoptimal geometric constraints may be one main source of the failure of the early techniques, which are similarly based on the Kruppa equation. Thirdly, Levenberg–Marquardt algorithm is adopted for the final refinement. This iteration is very fast since it is performed on only four unknown internal parameters.Another significant characteristic of the present method is using center-oriented metric coordinate system. Rather than CCD pixel coordinates, the metric coordinates are much more convenient for efficient approximation. Nevertheless, the pixel coordinates can be easily transformed to the metric coordinates and vice versa. For example, the pixel size can be simply assumed as4α×4μm, where α is the aspect ratio. The focal length and PPs can then be transferred to metric coordinates. This does not impact the metric reconstruction at all. The center-oriented coordinate is inspired by the obvious fact that the camera PPs is mostly around the image center. Although the center-oriented is not mandatory, it accelerates the convergence and is recommended.This work was partly presented in [30]. Its mathematician derivation is refined and extensive evaluations are performed in this work. The rest of the paper is organized as follows. The basic equations on projective geometry are described in the following subsection. The autocalibration method from N≥3 views is presented in Section 2. Section 3 demonstrates the performance of our method in the simulation and practical experiments. The work is concluded in the end.


### 1.1. Basic Equations

The homogeneousmetric coordinates of image measurements are indicated by x=[x,y,1]. The camera calibration matrix K is denoted in(1)K=f0x00αfy0001where f, α, x0, and y0 are the focal length, aspect ratio, and PPs, respectively. The skew parameter is considered as zero in (1).Fundamental matrixF is the algebraic representation of epipolar geometry between two views, independent of the scene structure. It has rank of two and satisfies(2)x′TFx=0where x and x′ are the corresponding image points in two views.The epipolar geometry is described by essential matrixE when the calibration matrices are known. An important property of essential matrix is that a 3×3 matrix is an essential matrix if and only if two of its singular values are equal, and the third is zero. This property, offering two constraints, can be represented algebraically as (3)2EETE-trETEE=0where tr(∙) denotes the trace operation. The relation between F and E is given in (4).(4)E=K′TFKFor the autocalibration problem, it is assumed that the calibration matrices in different views are unknown but constant, e.g.,K′=K in (4). Since there are four unknown internal parameters in (1) and two independent constraints in (3), two fundamental matrices shall be possible to perform autocalibration.


## 1.1. Basic Equations

The homogeneousmetric coordinates of image measurements are indicated by x=[x,y,1]. The camera calibration matrix K is denoted in(1)K=f0x00αfy0001where f, α, x0, and y0 are the focal length, aspect ratio, and PPs, respectively. The skew parameter is considered as zero in (1).Fundamental matrixF is the algebraic representation of epipolar geometry between two views, independent of the scene structure. It has rank of two and satisfies(2)x′TFx=0where x and x′ are the corresponding image points in two views.The epipolar geometry is described by essential matrixE when the calibration matrices are known. An important property of essential matrix is that a 3×3 matrix is an essential matrix if and only if two of its singular values are equal, and the third is zero. This property, offering two constraints, can be represented algebraically as (3)2EETE-trETEE=0where tr(∙) denotes the trace operation. The relation between F and E is given in (4).(4)E=K′TFKFor the autocalibration problem, it is assumed that the calibration matrices in different views are unknown but constant, e.g.,K′=K in (4). Since there are four unknown internal parameters in (1) and two independent constraints in (3), two fundamental matrices shall be possible to perform autocalibration.


## 2. N-View Autocalibration

Our method on autocalibration is detailed in this section. Firstly, three basic mathematical constraints are derived from (3) and (4). The distinctive characteristics of these three constraints allow using a recursive solution, which reduces the problem complexity. Proper simplified approximation is adopted to circumvent the numerical difficulties. Selection of the optimal geometric constraints is highlighted, which is crucial for accurate calibration, fast convergence, and robustness.


### 2.1. Mathematical Derivations

LetF=USVT by using singular value decomposition and denote(5)U=uij,V=vij,S=sij,i,j=1,2,3s1≜s11,s2≜s22NoticingUUT=VVT=1 and trAB=trBA, we obtain (6)-(7) from (3)-(5).(6)2KTFKKTFTKKTFK-trKTFTKKTFKKTFK=0(7)0=2SVTKKTVSUTKKTUS-trSVTKKTVSUTKKTUS≜2SNSMS-trSNSMSwhere M=MT=UTKKTU=mij and N=NT=VTKKTV=nij, i,j=1,2,3. KKT is named as the dual image of the absolute conic (DIAC), which was important in early literatures [19]. We get following equations from (7), without giving lengthy details on mathematical deduction:(8)s12m11n11-s22m22n22=0(9)s1m12n11+s2m22n12=0(10)s1m11n12+s2m12n22=0Obviously, these three equations, of which any two are independent, correspond to the rank of two of the fundamental matrix. For instance, (10) can be derived by inserting (8) into (9) by using the transform s2m22=-s1m12n11/n12.Equations (8)-(10) are the fundamental constraints in our method. They are named asConstraints I, II,andIII, respectively. They are the functions of the three calibration parameters and thus denoted by Cif,x0,y0i=1,2,3 or Cf,x0,y0 in general.Note that these constraints are related only withmij and nij (i,j=1,2). mij (i,j=1,2) are represented in (11)-(13) as follows:(11)m11=u112+α2u212f2+u31u31+2u11x0+2u21y0+2u11u21x0y0+u112x02+u212y02=u112+α2u212f2+u31+u11x0+u21y02(12)m12=m21=u11u12+α2u21u22f2+u31u32+u12u31+u11u32x0+u22u31+u21u32y0++u21u12+u11u22x0y0+u11u12x02+u21u22y02=u11u12+α2u21u22f2+u31+u11x0+u21y0u32+u12x0+u22y0(13)m22=u122+α2u222f2+u32u32+2u12x0+2u22y0+2u12u22x0y0+u122x02+u222y02=u122+α2u222f2+u32+u12x0+u22y02Similar formulae can be derived fornij (i,j=1,2), just replacing uij by vij in (11)-(13).


### 2.2. Recursive Solution

In the center-oriented metric coordinate system, the magnitude ofx0 and y0 are mostly around or smaller than 10−3 (unit: m), which corresponds to the PPs of lee than hundreds’ pixels. Thus, (14) holds for all uij≤1, by noticing U is an orthogonal matrix.(14)u31+u11x0+u21y02≈u312u32+u12x0+u22y02≈u322The benefit of approximation in (14) is that it does not change the sign of both sides. Using (14) and ignoring the terms of PPs, we have (15) and (16) as the approximations of (11) and (13), respectively.(15)m11≈u112+α2u212f2+u312(16)m22≈u122+α2u222f2+u322Similar formulae as (14)-(16) can be given for the matrix N. Note that this approximation is not applied to (12) since it may differ the sign of (12).Since Constraint I (8) contains terms mii and nii (i=1,2) only, it can be reduced to be a function of f and α from (15) and (16), when the information on PPs is unavailable. The reduced (8) is given in(17)s12u112+α2u212f2+u312v112+α2v212f2+v312-s22u122+α2u222f2+u322v122+α2v222f2+v322=0This elegant characteristic of Constraint I makes a recursive solution possible and feasible. Given N≥2 fundamental matrices, thef and α can be calculated from (17)-alike equations, without requiring any knowledge on PPs. Solving two (17)-alike equations is a plain work and the simple elimination technique is sufficient. Particularly, eliminating α gets a four-order polynomial of f2, and eliminating f obtains a six-order polynomial of α2. Subsequently, the estimates of f and α can be applied as fixed values into the constraints (8)-(10) to calculate the PPs, whose estimates can turn back to (8)-(10) to refine the calibration of f and α. This recursive process proceeds until convergence. The details on recursion are presented in the next sections.It is worth mentioning that (17) is exactly the reduced Constraint I (8) by assuming x0=y0=0. It is convenient for program implementation, with zero as the initial value of PPs. Despite this equivalence, it must be emphasized that the assumption of zero initialization is never the prerequisite for the recursive solution. Equation (17) is essentially based only on the metric coordinate system and the property of orthogonal matrix. The numerical experiments in Section 3 will also show that calibration is independent of the amount of PPs. This approximation (17), however, may not work well within the pixel coordinates, since in this case PPs are much larger than uij and vij and the approximations in (14)–(16) are less adequate.


### 2.3. PPs Calculation

Givenf and α, the constraints (8)-(10) can be represented as the function of x0 and y0, with the general form shown in (18). (18)a1x04+a2x03y0+a3x03+a4x02y02+a5x02y0+a6x02+a7x0y03+a8x0y02+a9x0y0+a10x0+a11y04+a12y03+a13y02+a14y0+a15=0In principle, any two (18)-alike equations can be used to calculate the PPs. Eliminating x0 from these two equations leads to a 44-order univariate polynomial of y0 and vice versa. However, deriving high-order polynomial is technically difficult and numerically unstable. Although the numerical methods based on the algebraic geometry techniques might be useful, they can be very slow and unstable [31]. The nonlinear iterative techniques, which require good initial values, may lead to poor convergence. In fact, a simple approximation can easily circumvent this nuisance, thanks to again the convenience of the metric coordinates. In the metric coordinates, the terms whose total degree is higher than two are numerically ignorable. Thus, (18) can be simplified as (19)a6x02++a9x0y0+a10x0+a13y02+a14y0+a15=0It becomes rather easy to solvingx0 and y0 from two (19)-alike equations. That is, eliminating x0 leads to a 4-order polynomial of y0 and vice versa.The numerical simulation shows, that (18) and (19) deliver almost the same real PPs solutions. It is stressed that this simplified approximation cannot be performed within the pixel coordinates.


### 2.4. Optimal Geometric Constraints

We know from last two subsection, that givenx0 and y0, any two constraints are sufficient to calculate f and α. As well, any two constraints are sufficient to calculate x0 and y0 by fixing f and α. Thus, multiple solutions can be derived from N≥2 fundamental matrices, since one fundamental matrix provides three constraints (8)-(10). In order to obtain a unique calibration, the proper constraints need to be selected. This can be accomplished by using the technique of error propagation analysis, which is described as follows.We perform the partial derivate for each constraint, as in(20)e=efdf2+eαdα2+exdx0+eydy0The amounts ofef, eα, ex, and ey represent the impact of erroneous f2, α2, x0, and y0 on the final accuracy. Nevertheless, the calibration parameters maintain different impacts on different constraints. Although it is hard to be proved rigorously, our simulations show the general rule as follows:(i)
Constraint I is robust to the deviation ofx0 and y0, while being sensitive to the change of f and α. Thus, Constraint I should be used to calculate the f and α.(ii)
Constraints II and III are robust to the deviation off and α, while being sensitive to the change of x0 and y0. Thus, Constraints II and III should be used to calculate the x0 and y0.This rule coincides to the analyses in Section2.2, where it is shown that the impact of PPs is insignificant on Constraint I. Therefore, it is needed only make analyses of error propagation on x0 and y0 for Constraint I and analyses on f and α for Constraints II and III.Since the amounts ofdx0 and dy0 are unknown, the error propagation e of Constraint I can be approximated by the normalized forms (21)ex0y0=exy0dx0x0+eyx0dy0y0(22)e2≈exy02+eyx02By consideringx0 and y0 are unavailable (or equal to zero) in the first recursion (see Section 2.2), (22) is reduced to be (23) in this case.(23)e2≈ex2+ey2Similarly, the error propagation off and α in Constraints II and III is approximated as (24)e2≈efα22+eαf22Then, any two constraints, whose error propagation isei2 and ej2, result in a univariate polynomial of one calibration parameter (see Sections 2.2 and 2.3). This polynomial is weighted as 1/ei2ej2. To avoid the potential dependence between Constraints II and III from single F, our strategy to select the optimal geometric constraints is as follows:(i)
Compare the weights of the Constraints II and III from singleF. The constraint of larger weight is selected. Thus, there are N constraints and N(N-1)/2 polynomials from N fundamental matrices.(ii)
Among N(N-1)/2 polynomials, the polynomials of topthree weights are selected as the optimal geometric constraints to calculate the unknown parameter. This step is applied to calculate all the calibration parameters: f, α, x0, and y0.Usually, the polynomial of the highest weight delivers the optimal solution. However, the second step of the strategy is performed to avoid the cases, where the multiview geometry is so weak that the polynomial of the highest weight might derive no rational solution.The rational solution must satisfy the following check.(25)f>0,α>0,-X≤x0≤X,-Y≤y0≤Ywhere 2X and 2Y denote the image length and width (unit: m).More practical considerations, like1/5<α<5, can be taken into account as well. Usually, more than one rational solution might be obtained after the rational check (25). Then the solution, which obtains the minimum mean square error (MSE) of all the polynomials is chosen as the unique solution.


### 2.5. Nonlinear Optimization

The calibration solution of the recursive procedure is usually suboptimal. In order to take advantage of the whole multiview geometry and further refine the calibration, the nonlinear optimization is recommended for final refinement, by using the recursive solution as the initial values. The Levenberg–Marquardt algorithm can adopted and can be run very fast since there are only four unknown internal parameters. As usual, the initial value of the damping factorλ can be set as 10−3. It is fairly expected that, with the good values obtained by our recursive method, the nonlinear optimization is able to obtain the global optimal solution.


### 2.6. Two-View Calibration

There are two independent constraints from two views (single fundamental matrix), Constraint I and Constraint II (or III). As mentioned previously in Sections2.2 and 2.4, Constraint I should be chosen to calculate the focal length by using (17), by given the aspect ratio only.Compared to the other techniques [10, 25, 26], which require known PPs and aspect ratio, our method can calculate the focal length from two views without any knowledge on the PPs. Thus, our method is most flexible in two-view application.


### 2.7. Summarizations

Our method is based on the fundamental matrix of epipolar geometry. It is beyond this work on how to derive the fundamental matrix from the image correspondences. Many methods are described in [19]. Although it is not optimal, the fast normalized 8-point algorithm [32] is employed in this study. Our method is summarized as follows.Inputs: image correspondences from N≥3 images, the image dimension2X and 2Y, maximum recursive number L (and zero skew).Outputs: Calibration matrixK.Algorithm:(1)
If the image correspondences,2X and 2Y, are in pixel coordinates, they are transformed to be those in metric coordinates (unit: m, by assuming the pixel size is 4α×4μm for instance).(2)
Compute the fundamental matrices of the image pairs by using the normalized 8-point algorithm or others [19].(3)
Setl=1,f0=α0=x00=y00=0.(4)
Use all the Constraints I to derive the uniquefl and αl by fixing x0l-1and y0l-1 (see Sections 2.2 and 2.4).(5)
Use all the Constraints II and III to derive the uniquex0l and y0l by fixing fl and αl (see Sections 2.3 and 2.4).(6)
If either (26) or (27) is satisfied, go to Step (6).(26)fl-fl-1fl<10-3,αl-αl-1αl<10-3,x0l-x0l-1x0l<10-3,y0l-y0l-1y0l<10-3(27)l>LElse,l=l+1 and go back to Step (4).(7)
Perform Levenberg–Marquardt algorithm on the all constraints for the final refinement.(8)
The calibration matrix in metric coordinates is transformed to be that in pixel coordinates, if necessary.Any prior information, if available, can be introduced into this method. For example, if aspect ratio is known, thenα is fixed in the above procedure; if the good initial values are available, they shall be used in the first recursion.


## 2.1. Mathematical Derivations

LetF=USVT by using singular value decomposition and denote(5)U=uij,V=vij,S=sij,i,j=1,2,3s1≜s11,s2≜s22NoticingUUT=VVT=1 and trAB=trBA, we obtain (6)-(7) from (3)-(5).(6)2KTFKKTFTKKTFK-trKTFTKKTFKKTFK=0(7)0=2SVTKKTVSUTKKTUS-trSVTKKTVSUTKKTUS≜2SNSMS-trSNSMSwhere M=MT=UTKKTU=mij and N=NT=VTKKTV=nij, i,j=1,2,3. KKT is named as the dual image of the absolute conic (DIAC), which was important in early literatures [19]. We get following equations from (7), without giving lengthy details on mathematical deduction:(8)s12m11n11-s22m22n22=0(9)s1m12n11+s2m22n12=0(10)s1m11n12+s2m12n22=0Obviously, these three equations, of which any two are independent, correspond to the rank of two of the fundamental matrix. For instance, (10) can be derived by inserting (8) into (9) by using the transform s2m22=-s1m12n11/n12.Equations (8)-(10) are the fundamental constraints in our method. They are named asConstraints I, II,andIII, respectively. They are the functions of the three calibration parameters and thus denoted by Cif,x0,y0i=1,2,3 or Cf,x0,y0 in general.Note that these constraints are related only withmij and nij (i,j=1,2). mij (i,j=1,2) are represented in (11)-(13) as follows:(11)m11=u112+α2u212f2+u31u31+2u11x0+2u21y0+2u11u21x0y0+u112x02+u212y02=u112+α2u212f2+u31+u11x0+u21y02(12)m12=m21=u11u12+α2u21u22f2+u31u32+u12u31+u11u32x0+u22u31+u21u32y0++u21u12+u11u22x0y0+u11u12x02+u21u22y02=u11u12+α2u21u22f2+u31+u11x0+u21y0u32+u12x0+u22y0(13)m22=u122+α2u222f2+u32u32+2u12x0+2u22y0+2u12u22x0y0+u122x02+u222y02=u122+α2u222f2+u32+u12x0+u22y02Similar formulae can be derived fornij (i,j=1,2), just replacing uij by vij in (11)-(13).


## 2.2. Recursive Solution

In the center-oriented metric coordinate system, the magnitude ofx0 and y0 are mostly around or smaller than 10−3 (unit: m), which corresponds to the PPs of lee than hundreds’ pixels. Thus, (14) holds for all uij≤1, by noticing U is an orthogonal matrix.(14)u31+u11x0+u21y02≈u312u32+u12x0+u22y02≈u322The benefit of approximation in (14) is that it does not change the sign of both sides. Using (14) and ignoring the terms of PPs, we have (15) and (16) as the approximations of (11) and (13), respectively.(15)m11≈u112+α2u212f2+u312(16)m22≈u122+α2u222f2+u322Similar formulae as (14)-(16) can be given for the matrix N. Note that this approximation is not applied to (12) since it may differ the sign of (12).Since Constraint I (8) contains terms mii and nii (i=1,2) only, it can be reduced to be a function of f and α from (15) and (16), when the information on PPs is unavailable. The reduced (8) is given in(17)s12u112+α2u212f2+u312v112+α2v212f2+v312-s22u122+α2u222f2+u322v122+α2v222f2+v322=0This elegant characteristic of Constraint I makes a recursive solution possible and feasible. Given N≥2 fundamental matrices, thef and α can be calculated from (17)-alike equations, without requiring any knowledge on PPs. Solving two (17)-alike equations is a plain work and the simple elimination technique is sufficient. Particularly, eliminating α gets a four-order polynomial of f2, and eliminating f obtains a six-order polynomial of α2. Subsequently, the estimates of f and α can be applied as fixed values into the constraints (8)-(10) to calculate the PPs, whose estimates can turn back to (8)-(10) to refine the calibration of f and α. This recursive process proceeds until convergence. The details on recursion are presented in the next sections.It is worth mentioning that (17) is exactly the reduced Constraint I (8) by assuming x0=y0=0. It is convenient for program implementation, with zero as the initial value of PPs. Despite this equivalence, it must be emphasized that the assumption of zero initialization is never the prerequisite for the recursive solution. Equation (17) is essentially based only on the metric coordinate system and the property of orthogonal matrix. The numerical experiments in Section 3 will also show that calibration is independent of the amount of PPs. This approximation (17), however, may not work well within the pixel coordinates, since in this case PPs are much larger than uij and vij and the approximations in (14)–(16) are less adequate.


## 2.3. PPs Calculation

Givenf and α, the constraints (8)-(10) can be represented as the function of x0 and y0, with the general form shown in (18). (18)a1x04+a2x03y0+a3x03+a4x02y02+a5x02y0+a6x02+a7x0y03+a8x0y02+a9x0y0+a10x0+a11y04+a12y03+a13y02+a14y0+a15=0In principle, any two (18)-alike equations can be used to calculate the PPs. Eliminating x0 from these two equations leads to a 44-order univariate polynomial of y0 and vice versa. However, deriving high-order polynomial is technically difficult and numerically unstable. Although the numerical methods based on the algebraic geometry techniques might be useful, they can be very slow and unstable [31]. The nonlinear iterative techniques, which require good initial values, may lead to poor convergence. In fact, a simple approximation can easily circumvent this nuisance, thanks to again the convenience of the metric coordinates. In the metric coordinates, the terms whose total degree is higher than two are numerically ignorable. Thus, (18) can be simplified as (19)a6x02++a9x0y0+a10x0+a13y02+a14y0+a15=0It becomes rather easy to solvingx0 and y0 from two (19)-alike equations. That is, eliminating x0 leads to a 4-order polynomial of y0 and vice versa.The numerical simulation shows, that (18) and (19) deliver almost the same real PPs solutions. It is stressed that this simplified approximation cannot be performed within the pixel coordinates.


## 2.4. Optimal Geometric Constraints

We know from last two subsection, that givenx0 and y0, any two constraints are sufficient to calculate f and α. As well, any two constraints are sufficient to calculate x0 and y0 by fixing f and α. Thus, multiple solutions can be derived from N≥2 fundamental matrices, since one fundamental matrix provides three constraints (8)-(10). In order to obtain a unique calibration, the proper constraints need to be selected. This can be accomplished by using the technique of error propagation analysis, which is described as follows.We perform the partial derivate for each constraint, as in(20)e=efdf2+eαdα2+exdx0+eydy0The amounts ofef, eα, ex, and ey represent the impact of erroneous f2, α2, x0, and y0 on the final accuracy. Nevertheless, the calibration parameters maintain different impacts on different constraints. Although it is hard to be proved rigorously, our simulations show the general rule as follows:(i)
Constraint I is robust to the deviation ofx0 and y0, while being sensitive to the change of f and α. Thus, Constraint I should be used to calculate the f and α.(ii)
Constraints II and III are robust to the deviation off and α, while being sensitive to the change of x0 and y0. Thus, Constraints II and III should be used to calculate the x0 and y0.This rule coincides to the analyses in Section2.2, where it is shown that the impact of PPs is insignificant on Constraint I. Therefore, it is needed only make analyses of error propagation on x0 and y0 for Constraint I and analyses on f and α for Constraints II and III.Since the amounts ofdx0 and dy0 are unknown, the error propagation e of Constraint I can be approximated by the normalized forms (21)ex0y0=exy0dx0x0+eyx0dy0y0(22)e2≈exy02+eyx02By consideringx0 and y0 are unavailable (or equal to zero) in the first recursion (see Section 2.2), (22) is reduced to be (23) in this case.(23)e2≈ex2+ey2Similarly, the error propagation off and α in Constraints II and III is approximated as (24)e2≈efα22+eαf22Then, any two constraints, whose error propagation isei2 and ej2, result in a univariate polynomial of one calibration parameter (see Sections 2.2 and 2.3). This polynomial is weighted as 1/ei2ej2. To avoid the potential dependence between Constraints II and III from single F, our strategy to select the optimal geometric constraints is as follows:(i)
Compare the weights of the Constraints II and III from singleF. The constraint of larger weight is selected. Thus, there are N constraints and N(N-1)/2 polynomials from N fundamental matrices.(ii)
Among N(N-1)/2 polynomials, the polynomials of topthree weights are selected as the optimal geometric constraints to calculate the unknown parameter. This step is applied to calculate all the calibration parameters: f, α, x0, and y0.Usually, the polynomial of the highest weight delivers the optimal solution. However, the second step of the strategy is performed to avoid the cases, where the multiview geometry is so weak that the polynomial of the highest weight might derive no rational solution.The rational solution must satisfy the following check.(25)f>0,α>0,-X≤x0≤X,-Y≤y0≤Ywhere 2X and 2Y denote the image length and width (unit: m).More practical considerations, like1/5<α<5, can be taken into account as well. Usually, more than one rational solution might be obtained after the rational check (25). Then the solution, which obtains the minimum mean square error (MSE) of all the polynomials is chosen as the unique solution.


## 2.5. Nonlinear Optimization

The calibration solution of the recursive procedure is usually suboptimal. In order to take advantage of the whole multiview geometry and further refine the calibration, the nonlinear optimization is recommended for final refinement, by using the recursive solution as the initial values. The Levenberg–Marquardt algorithm can adopted and can be run very fast since there are only four unknown internal parameters. As usual, the initial value of the damping factorλ can be set as 10−3. It is fairly expected that, with the good values obtained by our recursive method, the nonlinear optimization is able to obtain the global optimal solution.


## 2.6. Two-View Calibration

There are two independent constraints from two views (single fundamental matrix), Constraint I and Constraint II (or III). As mentioned previously in Sections2.2 and 2.4, Constraint I should be chosen to calculate the focal length by using (17), by given the aspect ratio only.Compared to the other techniques [10, 25, 26], which require known PPs and aspect ratio, our method can calculate the focal length from two views without any knowledge on the PPs. Thus, our method is most flexible in two-view application.


## 2.7. Summarizations

Our method is based on the fundamental matrix of epipolar geometry. It is beyond this work on how to derive the fundamental matrix from the image correspondences. Many methods are described in [19]. Although it is not optimal, the fast normalized 8-point algorithm [32] is employed in this study. Our method is summarized as follows.Inputs: image correspondences from N≥3 images, the image dimension2X and 2Y, maximum recursive number L (and zero skew).Outputs: Calibration matrixK.Algorithm:(1)
If the image correspondences,2X and 2Y, are in pixel coordinates, they are transformed to be those in metric coordinates (unit: m, by assuming the pixel size is 4α×4μm for instance).(2)
Compute the fundamental matrices of the image pairs by using the normalized 8-point algorithm or others [19].(3)
Setl=1,f0=α0=x00=y00=0.(4)
Use all the Constraints I to derive the uniquefl and αl by fixing x0l-1and y0l-1 (see Sections 2.2 and 2.4).(5)
Use all the Constraints II and III to derive the uniquex0l and y0l by fixing fl and αl (see Sections 2.3 and 2.4).(6)
If either (26) or (27) is satisfied, go to Step (6).(26)fl-fl-1fl<10-3,αl-αl-1αl<10-3,x0l-x0l-1x0l<10-3,y0l-y0l-1y0l<10-3(27)l>LElse,l=l+1 and go back to Step (4).(7)
Perform Levenberg–Marquardt algorithm on the all constraints for the final refinement.(8)
The calibration matrix in metric coordinates is transformed to be that in pixel coordinates, if necessary.Any prior information, if available, can be introduced into this method. For example, if aspect ratio is known, thenα is fixed in the above procedure; if the good initial values are available, they shall be used in the first recursion.


## 3. Experiments

Both simulated and practical experiments are carried out to evaluate our method. All the experiments are implemented in MATLAB, in one PC with 32-bit Windows Seven system, CPU 2.5GHz, 3GB RAM.


### 3.1. Simulation Experiments

The presented autocalibration method is evaluated in simulation studies, with respect to different noise levels, PPs, and number of image correspondences. The image resolution of the virtual camera is 2000×1600, focal length 2000 pixels (angle of view 52°×43°), and aspect ratio 1.2. PPs vary in the experiments. The pixel size is thus assumed as4α×4μm and all the image correspondences transformed to be the metric coordinates. Three typically posed images are simulated.Each simulation experiment is run 100 times for statistical purpose. The accuracy of focal length and aspect ratio is quantified by relative error (%), and PPs are evaluated by absolute error (unit: pixel). The fundamental matrices are computed by using the fast normalized 8-point algorithm in both simulation and practical experiments.Given Perfect Fundamental Matrices. Since the normalized 8-point algorithm is not a mandatory part, we first evaluate the performance of our calibration method when all fundamental matrices are perfectly given. The results are shown in Figure 1, with respect to PPs varying from 0 pixel to 150 pixels. Both the relative errors of f and α are impressively smaller than 10−7%, and the absolute error of PPs is less than 2×10−6 pixel. The running time is around 0.1 second in all cases. Do note that the performance is independent of the PPs amount.Figure 1
Performance evaluation with respect to PPs variations, by given perfect fundamental matrices (left:f and α (%); middle: PPs (pixels); right: running time (sec)).Noise Effect. The computation of fundamental matrix is impacted by the noise, whose effect on our calibration method is evaluated here. There are 100 correspondences for each image pair and the real PPs are 50 pixels. The results are illustrated in Figure 2 with the noise increasing from 0.0 to 1.0 pixel.Figure 2
Performance evaluation with respect to different noise levels (left:f and α (%); middle: PPs (pixels); right: running time (sec)).It is clear from Figure2 that the relative errors of f and α and the absolute error of PPs increase linearly with the noise level. For the noise smaller than 0.2 pixel, the relative errors of f and α are smaller than 0.1% and the PPs error is smaller than 5 pixels. The running time is around 0.1 second, independent of the noise level.Number of Image Correspondences. The computation of fundamental matrix is impacted as well by the number of the image correspondences, whose effect on calibration is illustrated in Figure 3. The noise level is 0.1 pixel and the real PPs are 50 pixels.Figure 3
Performance evaluation with respect to the number of image correspondences (left:f and α (%); middle: PPs (pixels); right: running time (sec)).From Figure3, the calibration accuracy increase substantially with the number of correspondences. When there are around 100 image correspondences, the relative errors of f and α are around 0.05% and the absolute errors of PPs are 1.5 pixels. The running time is about 0.05 second in all cases.PPs Effect. Here we again demonstrate in Figure 4 the calibration results with respect to PPs, under the noise of 0.1 pixel and 100 correspondences. The real PPs vary from 10 to 150 pixels. From Figure 4, the relative errors of f and α are around 0.05% and the absolute error of PPs is smaller than 2 pixels. It is obvious that the calibration performance is fully independent of the PPs deviation, coincident to the analyses in Section 2.2 and the results in Section 3.1.Figure 4
Performance evaluation with respect to PPs (left:f and α (%); middle: PPs (pixels); right: running time (sec)).


### 3.2. Practical Experiments

We carry out practical experiments for autocalibration from two datasets. One dataset is taken in our lab environment and another is from the EPFL public vision dataset [8].In the first experiment, three images are taken from our lab as shown in Figure5. These images are almost (not fully) distortion free. It must be stressed that NO target information is used in this experiment, though these images look a bit like the checkerboard in Zhang’s method [9]. We implement such experiment since very accurate and robust coded correspondences can be obtained by using the commercial software Australis (http://www.solve3d.net/ENG/Aus_eproducts.htm) which gets the coded targets with noise level less than 0.1 pixels. The autocalibration results are compared to those derived by the standard bundle adjustment. It must be noticed that the good initial values of calibration parameters and external orientation are indispensable for bundle adjustment method. The true aspect ratio is fixed in this adjustment implementation.Figure 5
Three practical images of the first dataset for autocalibration tests.The results in the first dataset are shown in Table1. The relative differences of f and α between autocalibration and bundle adjustment are 0.15 % and 0.067%, respectively. The absolute differences of x0 and y0 are 1.24 and 0.63 pixels. The running time of autocalibration is 0.15 second while it is 3.20 seconds for bundle adjustment. The autocalibration accuracy, referenced to the adjustment results, coincides with that in the simulation study.Table 1
Calibration results of bundle adjustment and our autocalibration method in the first dataset.
Parameter α f x 0 y 0 Time (sec) Bundle adjustment 1.00 (fixed) 2302.82 1238.55 1027.04 3.20  Autocalibration 0.99933 2306.26 1237.66 1027.67 0.15In the second experiment, from EPFL dataset we select three images of the castle, fountain, and herzjesu collections, respectively, as illustrated in Figure6. These three sets of images are taken using a single camera and the image resolution is 3072×2048. Image correspondences are obtained by detecting the SIFT features and matching features using the RASAC strategy with constraint of the fundamental matrix.Figure 6
Practical images of the second dataset from EPFL (from top: castle, fountain, and herzjesu).The calibration results are demonstrated in Table2. It is found that our method obtains very accurate calibration results, and the average errors of f, x0, and y0 reduce to 11.8, 15.0, and 35.2 pixels, respectively. The good performance of our method is thus confirmed in these experiments.Table 2
Calibration results of our autocalibration method in the EPFL dataset.
data/Parameter α f x 0 y 0 Time (sec) Ground truth 1 2362.12 1520.69 1006.81  Castle 0.993 2375.7 1535.5 966.2 0.18  Fountain 0.996 2351.5 1509.7 980.3 0.15  Herzjesu 0.991 2352.8 1501.2 1045.4 0.21It is also observed that the errors ofy0 calibration (vertical direction) are much larger than those of x0 (horizontal direction). This is mainly due to the geometry configurations that the rotation variations in y0 direction are much less than those in x0 direction in all the three image sets. In other words, inferior calibration of y0 is mainly caused by the fact that the calibration geometry in y0 direction is much weaker than x0 direction.


## 3.1. Simulation Experiments

The presented autocalibration method is evaluated in simulation studies, with respect to different noise levels, PPs, and number of image correspondences. The image resolution of the virtual camera is 2000×1600, focal length 2000 pixels (angle of view 52°×43°), and aspect ratio 1.2. PPs vary in the experiments. The pixel size is thus assumed as4α×4μm and all the image correspondences transformed to be the metric coordinates. Three typically posed images are simulated.Each simulation experiment is run 100 times for statistical purpose. The accuracy of focal length and aspect ratio is quantified by relative error (%), and PPs are evaluated by absolute error (unit: pixel). The fundamental matrices are computed by using the fast normalized 8-point algorithm in both simulation and practical experiments.Given Perfect Fundamental Matrices. Since the normalized 8-point algorithm is not a mandatory part, we first evaluate the performance of our calibration method when all fundamental matrices are perfectly given. The results are shown in Figure 1, with respect to PPs varying from 0 pixel to 150 pixels. Both the relative errors of f and α are impressively smaller than 10−7%, and the absolute error of PPs is less than 2×10−6 pixel. The running time is around 0.1 second in all cases. Do note that the performance is independent of the PPs amount.Figure 1
Performance evaluation with respect to PPs variations, by given perfect fundamental matrices (left:f and α (%); middle: PPs (pixels); right: running time (sec)).Noise Effect. The computation of fundamental matrix is impacted by the noise, whose effect on our calibration method is evaluated here. There are 100 correspondences for each image pair and the real PPs are 50 pixels. The results are illustrated in Figure 2 with the noise increasing from 0.0 to 1.0 pixel.Figure 2
Performance evaluation with respect to different noise levels (left:f and α (%); middle: PPs (pixels); right: running time (sec)).It is clear from Figure2 that the relative errors of f and α and the absolute error of PPs increase linearly with the noise level. For the noise smaller than 0.2 pixel, the relative errors of f and α are smaller than 0.1% and the PPs error is smaller than 5 pixels. The running time is around 0.1 second, independent of the noise level.Number of Image Correspondences. The computation of fundamental matrix is impacted as well by the number of the image correspondences, whose effect on calibration is illustrated in Figure 3. The noise level is 0.1 pixel and the real PPs are 50 pixels.Figure 3
Performance evaluation with respect to the number of image correspondences (left:f and α (%); middle: PPs (pixels); right: running time (sec)).From Figure3, the calibration accuracy increase substantially with the number of correspondences. When there are around 100 image correspondences, the relative errors of f and α are around 0.05% and the absolute errors of PPs are 1.5 pixels. The running time is about 0.05 second in all cases.PPs Effect. Here we again demonstrate in Figure 4 the calibration results with respect to PPs, under the noise of 0.1 pixel and 100 correspondences. The real PPs vary from 10 to 150 pixels. From Figure 4, the relative errors of f and α are around 0.05% and the absolute error of PPs is smaller than 2 pixels. It is obvious that the calibration performance is fully independent of the PPs deviation, coincident to the analyses in Section 2.2 and the results in Section 3.1.Figure 4
Performance evaluation with respect to PPs (left:f and α (%); middle: PPs (pixels); right: running time (sec)).


## 3.2. Practical Experiments

We carry out practical experiments for autocalibration from two datasets. One dataset is taken in our lab environment and another is from the EPFL public vision dataset [8].In the first experiment, three images are taken from our lab as shown in Figure5. These images are almost (not fully) distortion free. It must be stressed that NO target information is used in this experiment, though these images look a bit like the checkerboard in Zhang’s method [9]. We implement such experiment since very accurate and robust coded correspondences can be obtained by using the commercial software Australis (http://www.solve3d.net/ENG/Aus_eproducts.htm) which gets the coded targets with noise level less than 0.1 pixels. The autocalibration results are compared to those derived by the standard bundle adjustment. It must be noticed that the good initial values of calibration parameters and external orientation are indispensable for bundle adjustment method. The true aspect ratio is fixed in this adjustment implementation.Figure 5
Three practical images of the first dataset for autocalibration tests.The results in the first dataset are shown in Table1. The relative differences of f and α between autocalibration and bundle adjustment are 0.15 % and 0.067%, respectively. The absolute differences of x0 and y0 are 1.24 and 0.63 pixels. The running time of autocalibration is 0.15 second while it is 3.20 seconds for bundle adjustment. The autocalibration accuracy, referenced to the adjustment results, coincides with that in the simulation study.Table 1
Calibration results of bundle adjustment and our autocalibration method in the first dataset.
Parameter α f x 0 y 0 Time (sec) Bundle adjustment 1.00 (fixed) 2302.82 1238.55 1027.04 3.20  Autocalibration 0.99933 2306.26 1237.66 1027.67 0.15In the second experiment, from EPFL dataset we select three images of the castle, fountain, and herzjesu collections, respectively, as illustrated in Figure6. These three sets of images are taken using a single camera and the image resolution is 3072×2048. Image correspondences are obtained by detecting the SIFT features and matching features using the RASAC strategy with constraint of the fundamental matrix.Figure 6
Practical images of the second dataset from EPFL (from top: castle, fountain, and herzjesu).The calibration results are demonstrated in Table2. It is found that our method obtains very accurate calibration results, and the average errors of f, x0, and y0 reduce to 11.8, 15.0, and 35.2 pixels, respectively. The good performance of our method is thus confirmed in these experiments.Table 2
Calibration results of our autocalibration method in the EPFL dataset.
data/Parameter α f x 0 y 0 Time (sec) Ground truth 1 2362.12 1520.69 1006.81  Castle 0.993 2375.7 1535.5 966.2 0.18  Fountain 0.996 2351.5 1509.7 980.3 0.15  Herzjesu 0.991 2352.8 1501.2 1045.4 0.21It is also observed that the errors ofy0 calibration (vertical direction) are much larger than those of x0 (horizontal direction). This is mainly due to the geometry configurations that the rotation variations in y0 direction are much less than those in x0 direction in all the three image sets. In other words, inferior calibration of y0 is mainly caused by the fact that the calibration geometry in y0 direction is much weaker than x0 direction.


## 4. Conclusions

We presented an efficient method for camera autocalibration from N≥3 views, by given zero skew and image correspondences only. This method is based on fundamental matrix only. The focal length, aspect ratio, and PPs can be computed from N≥2 fundamental matrices. Our method contains two main parts: analytical recursive solution and final nonlinear optimization. It is fast (around 0.1 second from 3 views in all our tests) and flexible. Any prior information on internal parameters can be easily introduced in our method. As a byproduct, an analytical variation of our method can be applied to calculate the focal length from two views, without requiring any knowledge on principal points.Both the simulation and practical studies show our autocalibration method is quite efficient. Generally speaking, the performance of autocalibration is dependent on the quality of fundamental matrices and multiview geometry, but independent of the PPs deviation. On the one hand, the computation of fundamental matrix is influenced by the noise level and the number of the image correspondences, whose impact on autocalibration was demonstrated in this paper. The fundamental matrix can also be influenced by the lens distortion. On the other hand, multiview geometry has significant influences on the calibration. Since autocalibration exploits only the image correspondences only, it may require relatively stronger geometry than the other calibration methods. For example, three images at least are needed for autocalibration; pure camera translation (without rotation) and planar scene can definitely corrupt the autocalibration. The impact of the view geometry on autocalibration will be studied in our next work.Calibration is an indispensable prerequisite for orientation and reconstruction. This work makes it practically possible and feasible to perform calibration, orientation, and reconstruction, by using only the multiview image correspondences. One future work is further improving the robustness of the present method.



---

*Source: 8348340-2018-11-25.xml*