# A Decoupled Calibration Method for Camera Intrinsic Parameters and Distortion Coefficients

Authors: Kun Yan; Hong Tian; Enhai Liu; Rujin Zhao; Yuzhen Hong; Dan Zuo
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Camera calibration is a necessary process in the field of vision measurement. In this paper, we propose a flexible and high-accuracy method to calibrate a camera. Firstly, we compute the center of radial distortion, which is important to obtain optimal results. Then, based on the radial distortion of the division model, the camera intrinsic parameters and distortion coefficients are solved in a linear way independently. Finally, the intrinsic parameters of the camera are optimized via the Levenberg-Marquardt algorithm. In the proposed method, the distortion coefficients and intrinsic parameters are successfully decoupled; calibration accuracy is further improved through the subsequent optimization process. Moreover, whether it is for relatively small image distortion or distortion larger image, utilizing our method can get a good result. Both simulation and real data experiment demonstrate the robustness and accuracy of the proposed method. Experimental results show that the proposed method can be obtaining a higher accuracy than the classical methods.


---

## Body

## 1. Introduction

Camera calibration is an important part of the application of photogrammetry that aims to compute the camera model parameters from two-dimensional images [1, 2]. Camera model parameters include intrinsic and extrinsic parameters. The intrinsic parameters describe the geometry of imaging process, while the extrinsic parameters indicate camera position and attitude in the world coordinate system. Camera calibration precision directly affects the measuring accuracy of vision measurement system. Therefore, the study of a flexible and high-precision camera calibration method has very important significance.Nowadays, the techniques of camera calibration can be divided into two categories: the traditional camera calibration and the self-calibration of the camera. The traditional camera calibration methods use scene information, including the points or lines with precise coordinates, to solve camera parameters, while self-calibration methods only use the relationship between the sequences of images to solve the camera parameters. In general, to satisfy the high-accuracy requirement, we always make use of the traditional camera calibration method. Besides, the traditional method consists of three types: the linear method, nonlinear optimization method, and the two-step method. Hall et al. [3] introduced the first linear method by computing the 3×4 transformation matrix based on the pinhole model. Later, the nonlinear calibration method is developed by introducing a variety of lens distortion [4–10]. Zhang [11] proposed a flexible calibration technique for desktop vision system (DVS) by using a printed planar calibration pattern, which is a typical representative of the two-step method.Recently, more research is dedicated to improving the performance of the camera calibration. Wang et al. [12] proposed a new calibration model of camera lens distortion, which is according to a transform from ideal image plane to real sensor array plane. Ahmed and Farag [13] proposed a robust approach to distortion calibration by using a least-median-of-squares estimator based on the analysis of distorted straight lines in the images. Later, Ricolfe-Viala and Sánchez-Salmerón [14] presented a robust metric calibration method with nonlinear camera lens distortion, which computes the camera lens distortion isolated from the camera calibration process under stable conditions, independently of the computed lens distortion model or the number of parameters. Many other researchers [15–19] also proposed the nonlinear objective function in three-dimensional space or distortion free space to minimize the calibration error to improve the accuracy of the calibration.However, in the above literatures, most of the calibration methods employ the same nonlinear optimization methods. And almost completely, the lens distortion coefficients and other intrinsic and extrinsic camera parameters are estimated in an optimization framework at the same time. According to the report of Hartley [20], such nonlinear iteration can be troublesome and may converge to the local minima without selecting a good initial value. Moreover, due to the distortion often coupling parameters in the camera’s internal parameters and external parameters, methods which extend the calibration of the pinhole model to obtain the camera distortion parameters lead to high internal parameters miscalculation [7]. In addition, the coupling between the different parameters can make the estimation results quite unreliable [7, 13]. Therefore, it is important to use a different method to estimate the camera distortion coefficients apart from the pinhole model. In [14], Ricolfe-Viala and Sánchez-Salmerón computed lens distortion model or the number of distortion parameters independently. But they solved the distortion center by a nonlinear minimization method. Although Ahmed and Farag [13] successfully decoupled intrinsic parameters and distortion coefficient, they relied on the fact that straight lines in the scene must always perspectively project to straight lines in the image. Moreover, they assumed that the distortion center is known. In [21], Fitzgibbon proposed a noniterative method to estimate radial distortion by introducing the division model. He also successfully decoupled intrinsic parameters and distortion coefficient. However, the radial distortion includes only one parameter in his work.In this paper, a high-accuracy calibration method is proposed. Different from [13, 14, 21], this method computes the center of radial distortion firstly, which is important to obtain optimal results. Afterwards, based on the division model, the interior parameter and distortion coefficients of the camera are estimated by linear method, which provide a good initial value for the subsequent optimization. Moreover, our method can calculate camera intrinsic parameters and any number of distortion coefficients. In addition, the intrinsic parameters of the camera are optimized via the Levenberg-Marquardt algorithm, and then we obtain more accurate results. Experimental results show that the proposed method can be obtaining a higher accuracy than the classical methods.This paper is organized as follows. Section2 gives a brief description of the camera model and the lens distortion model. Section 3 describes the detail procedure of the proposed camera calibration method. Section 4 verifies the proposed calibration method by simulation experiments and real data experiments. At last, this paper ends with several important conclusions in Section 5.


## 2. Camera Model

The camera model is shown in Figure1. Given one homogeneous coordinates point Pw=(Xw,Yw,Zw,1)T in 3D space. Pu=(Xu,Yu,1)T is the perspective projection of Pw in the image plane based on pinhole model without lens distortion. Pd=(Xd,Yd,1)T is the real projection of Pw in the image plane of camera considering the lens distortion.Figure 1
The model of camera.Without consideration of lens distortion, the mapping between the 3D pointPw and 2D image point Pu is given by(1)sPu=KRtPw,K=fusku00fvv0001,where s represents the nonzero scale factors; K is the camera intrinsic matrix, with (u0,v0) being the coordinates of the principal point, fu and fv being the effective focal length in pixels, and sk being the parameter describing the skew of the two image axes; Rt represent the rotation and translation vector from the world coordinate system to the camera coordinate system, respectively.Because of the plane calibration model used in our calibration procedure, we assume that the model planeZw=0 with no loss of generality. From (1), we have(2)sXuYu1=Kr1r2r3tXwYw01=Kr1r2tXwYw1,where ri is the ith column of the rotation matrix R. By abuse of notation, we still denote the point on the model plane by Pw=(Xw,Yw,1)T. Therefore, 3D point Pw and its image point Pu are related by a homography H,(3)sPu=HPw,H=Kr1r2t.With the influence of lens distortion, as shown in Figure1, the actual image point Pd is not the point Pu, which is the image plane intersection of the connection of the 3D point Pw and the optical center Oc, but it has some deviation. Theoretical calculations will be carried out through the actual image coordinates after correction for the ideal image coordinates; the calibration process can be achieved through distortion compensation model. Camera lens distortion was first introduced by Conrady in 1919 with the decentring lens distortion. Afterwards, Brown [22] proposed the radial, decentring, and prism distortion model which has been widely used. In general, the radial distortion is sufficient for a high-accuracy measurement. Brown proposed the Polynomial Model (PM), which is the most popular model to describe radial distortion:(4)Pu-e=Pd-e·Lrd,k,withLrd,k=1+k1rd2+k2rd4+⋯,where k1,k2,… are the distortion coefficients, e=(du0,dv0,1)T is the homogeneous coordinate of the center of distortion (COD), and rd is the pixel radius to e.PM works best for lens with small distortions. For the wide-angle lens or fisheye lens, there is considerable distortion and they often need too many terms compared to the actual. So, Fitzgibbon [21] suggested the division model (DM) to describe the radial distortion model,(5)Pu-e=Pd-eLrd,k,where L(rd,k) is the same as in (4). Combining (4) and (5) we can get a more generic Rational Model (RM):(6)Pu-e=Pd-e·L1rd,k1L2rd,k2.However, the calculation speed of the imaging model based on the RM is relatively slow, and the convergence results are more dependent on the accuracy of the initial value. On the other hand, the most remarkable advantage of the DM over the PM is that it is able to express high distortion at much lower order. So, our work for camera calibration is based on the DM.


## 3. Camera Calibration Procedure

In this section, the process of the proposed calibration method is described in detail. The whole process is divided into four steps. Firstly, in order to make the center of distortion (COD) to the original, the COD is estimated accurately and the calibration image is corrected. Then the corresponding points between the model plane and its image are used to compute the homography and distortion coefficients simultaneously based on the radial distortion division model (DM) [21]. In the third step, all the camera intrinsic parameters are analytically calculated from the homographies of calibration images. Finally, the intrinsic parameters of the camera are optimized via the Levenberg-Marquardt algorithm.


### 3.1. Find the COD

In the literatures [9, 14, 23], the researchers assume that the COD is the center of the image or the principal point of the camera. However, we all know that this is not a good assumption. The COD can be displaced by a number of factors from the center of the image, such as an offset from the center of the CCD lens, a tilt with respect to the plane of the lens sensor, or the installation of a combination of lenses and image cropping. In general consumer level cameras (such as several hundred dollar costs), it should not be assumed that the photocenters of the camera are accurate, because these effects are not much different from the subjective image quality [24].Experiments have been done by Hartley and Zisserman [25] to show that the general assumption that the COD is at the center of the image is not exactly right. Besides, they have also proved that the point may be remarkably deviated from the center of the image, or the principal point of the camera. So the accurate estimation of COD is important to obtain optimal results in camera calibration. Here, we adopt the novel method proposed in [24] to accurately estimate the position of COD in this paper. This method is simple, but it produces perfect results. We briefly introduce this approach here. For more details, the reader can refer to the literature [24]. The known point Pw on the calibration model plane and its corresponding image point Pd are related by a radial fundamental relationship, and it can be represented as follows:(7)PdTe×HPw=0,where H is the same as in (3) and [e]× (the skew symmetric matrix) represents the cross-product. Writing Fr=[e]×H (in [24], the matrix Fr is called as radial fundamental matrix), then (7) can be rewritten as(8)PdTFrPw=0.Note thatFr can be calculated in a general way [25] (such as the eight-point algorithm) from some corresponding points. Then the COD can be estimated from its left epipole:(9)eTFr=eTe×H=0.Of course, in the absence of radial deformation, the basic matrix of the above calculation is not stable, and the estimated value ofe is essentially arbitrary and meaningless. Therefore, if there is no radial distortion, then it is not much of a sense to talk about the distorted center [24].


### 3.2. Compute the Distortion Coefficients and Homography

We will compute the distortion coefficients and homography in this section. Based on the reports in the literatures [6, 22, 26], it is likely that the distortion function is completely dominated by the radial components, in particular the dominant first term. In addition, any more detailed modeling not only does not help (ignore when compared with sensor quantization), but also leads to numerical instability [6, 26].After getting the CODe, we can compute the H by factoring the matrix Fr (Fr=[e]×H). Although [e]× is singular, this decomposition is not unique. Here, a novel method is presented to solve this problem. In doing this, we change the coordinates of the image, so that the COD is the original. As is shown in Figure 1, the image plane coordinate system Oxyxy is transformed to OdX⌢dY⌢d. Where Od is the real position of the COD, we represent it as e⌢=001T in homogeneous coordinates. Then the transformed distorted image point P⌢d can be represented as (10)P⌢d=X⌢dY⌢d1=Xd-du0Yd-dv01=10-du001-dv0001XdYd1=10-du001-dv0001Pd.Similarly, the transformed undistorted image pointP⌢u and the original undistorted image point Pu are described as(11)P⌢u=10-du001-dv0001Pu.When the coordinate origin is converted to the center of distortion, the homography matrix is redefined asH⌢ and the corresponding radial fundamental matrix is defined as F⌢r=[e⌢]×H⌢. Then, from (8), we can get(12)P⌢dTF⌢rPw=0.Using (10) into (12), we can get(13)PdT10-du001-dv0001TF⌢rPw=0.Besides, from (8) and (13), we can get(14)Fr=10-du001-dv0001TF⌢r.F ⌢ r can be calculated as(15)F⌢r=100010du0dv01Fr.Note thatF⌢r=e⌢×H⌢ and e⌢=001T; we have(16)F⌢r=e⌢×H⌢=0-10100000H⌢.LetF⌢r=[f⌢1T;f⌢2T;f⌢3T] and H⌢=[h⌢1T;h⌢2T;h⌢3T]. Here, f⌢iT=f⌢i1f⌢i2f⌢i3 is the ith row of the radial fundamental matrix F⌢r, and h⌢iT=h⌢i1h⌢i2h⌢i3 is the ith row of matrix H⌢. Due to the fact that final row of F⌢r is zero, so we only need to solve the first two rows of H⌢ from (16), that is,(17)h⌢1T=f⌢2T,h⌢2T=-f⌢1T.Then by calculating the homographyH⌢, present up to only three unknown parameters h⌢3T=h⌢31h⌢32h⌢33. As the homography H⌢ relates the model point Pw and its undistorted image point P⌢u,(18)sP⌢u=H⌢Pw.Note that P⌢u=X⌢uY⌢u1 is the undistorted image point after the origin of the image coordinates is transferred to the COD. Considering the division model (DM), we have(19)X⌢u=X⌢dLr⌢d,k,Y⌢u=Y⌢dLr⌢d,k,where L(r⌢d,k) is the same as in (4), that is,(20)Lr⌢d,k=1+k1r⌢d2+k2r⌢d4+⋯.We substitute (19) into (18); then we can get(21)
s
X
⌢
d
L
r
⌢
d
,
k
Y
⌢
d
L
r
⌢
d
,
k
1
=
h
⌢
1
T
h
⌢
2
T
h
⌢
3
T
P
w
.So each point can be given the two equations,(22)X⌢dh⌢3TPw-h⌢1TPw·k1r⌢d2+k2r⌢d4+⋯Y⌢dh⌢3TPw-h⌢2TPw·k1r⌢d2+k2r⌢d4+⋯=h⌢1TPwh⌢2TPw.Using (17) in (22), we can get(23)X⌢dPwT-f⌢2TPwr⌢d2r⌢d4⋯Y⌢dPwTf⌢1TPwr⌢d2r⌢d4⋯h⌢3k1k2⋮=f⌢2TPw-f⌢1TPw.Given N corresponding point and that it satisfies that 2N≥n+3, where n is the number of distortion parameters, we can solve the system of equations in the least squares sense. Then the third row of H⌢ (i.e., h⌢3T) is solved. Moreover, the distortion coefficients k1,k2,… are computed at the same time.


### 3.3. Linear Solution of Intrinsic Parameters of the Camera

AfterH⌢ is obtained, the homography matrix H for the original image coordinates can be computed very easily. From (3), (11), and (18), we can get(24)H=10-du001-dv0001-1H⌢=10du001dv0001H⌢.Let us denote homography matrix byH=hc1hc2hc3 (where hci is the ith column of H). According to the orthogonality of the rotation matrix R(r1Tr2=0,r1Tr1=r2Tr2), we get(25)hc1TK-TK-1hc2=0,hc1TK-TK-1hc1=hc2TK-TK-1hc2.Note thatK is the camera intrinsic matrix. 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. Therefore, we need at least three images to calculate all the parameters. After the homography H is obtained, the analytical solution of the intrinsic parameters of the camera can be solved in the literature [11].


### 3.4. Maximum Likelihood Estimation

The intrinsic parameters obtained in the previous section are not accurate enough due to image noise. We can refine it through maximum likelihood estimation. Ahmed and Farag [13] have proved that including the distortion center and the decentering coefficients in the nonlinear optimization step may lead to instability of the estimation algorithm. Therefore, we do not need to optimize the distortion coefficient and distortion center. Meanwhile, it can reduce the search space of the calibration problem without sacrificing the accuracy and produce more stable and noise-robust results. Here, we give 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. Then the maximum likelihood estimate can be obtained by minimizing the following function:(26)∑i=1n∑j=1mPij-P⌢K,Ri,ti,Pj2,where P⌢(K,Ri,ti,Pj) is the projection of point Pj in image i according to (3). 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. Minimizing (26) is a nonlinear minimization problem, which is solved by the Levenberg-Marquardt algorithm [27]. It requires an initial value of K and {Ri,ti∣i=1,…,n}, which can be obtained using the technique described in the previous subsection.


## 3.1. Find the COD

In the literatures [9, 14, 23], the researchers assume that the COD is the center of the image or the principal point of the camera. However, we all know that this is not a good assumption. The COD can be displaced by a number of factors from the center of the image, such as an offset from the center of the CCD lens, a tilt with respect to the plane of the lens sensor, or the installation of a combination of lenses and image cropping. In general consumer level cameras (such as several hundred dollar costs), it should not be assumed that the photocenters of the camera are accurate, because these effects are not much different from the subjective image quality [24].Experiments have been done by Hartley and Zisserman [25] to show that the general assumption that the COD is at the center of the image is not exactly right. Besides, they have also proved that the point may be remarkably deviated from the center of the image, or the principal point of the camera. So the accurate estimation of COD is important to obtain optimal results in camera calibration. Here, we adopt the novel method proposed in [24] to accurately estimate the position of COD in this paper. This method is simple, but it produces perfect results. We briefly introduce this approach here. For more details, the reader can refer to the literature [24]. The known point Pw on the calibration model plane and its corresponding image point Pd are related by a radial fundamental relationship, and it can be represented as follows:(7)PdTe×HPw=0,where H is the same as in (3) and [e]× (the skew symmetric matrix) represents the cross-product. Writing Fr=[e]×H (in [24], the matrix Fr is called as radial fundamental matrix), then (7) can be rewritten as(8)PdTFrPw=0.Note thatFr can be calculated in a general way [25] (such as the eight-point algorithm) from some corresponding points. Then the COD can be estimated from its left epipole:(9)eTFr=eTe×H=0.Of course, in the absence of radial deformation, the basic matrix of the above calculation is not stable, and the estimated value ofe is essentially arbitrary and meaningless. Therefore, if there is no radial distortion, then it is not much of a sense to talk about the distorted center [24].


## 3.2. Compute the Distortion Coefficients and Homography

We will compute the distortion coefficients and homography in this section. Based on the reports in the literatures [6, 22, 26], it is likely that the distortion function is completely dominated by the radial components, in particular the dominant first term. In addition, any more detailed modeling not only does not help (ignore when compared with sensor quantization), but also leads to numerical instability [6, 26].After getting the CODe, we can compute the H by factoring the matrix Fr (Fr=[e]×H). Although [e]× is singular, this decomposition is not unique. Here, a novel method is presented to solve this problem. In doing this, we change the coordinates of the image, so that the COD is the original. As is shown in Figure 1, the image plane coordinate system Oxyxy is transformed to OdX⌢dY⌢d. Where Od is the real position of the COD, we represent it as e⌢=001T in homogeneous coordinates. Then the transformed distorted image point P⌢d can be represented as (10)P⌢d=X⌢dY⌢d1=Xd-du0Yd-dv01=10-du001-dv0001XdYd1=10-du001-dv0001Pd.Similarly, the transformed undistorted image pointP⌢u and the original undistorted image point Pu are described as(11)P⌢u=10-du001-dv0001Pu.When the coordinate origin is converted to the center of distortion, the homography matrix is redefined asH⌢ and the corresponding radial fundamental matrix is defined as F⌢r=[e⌢]×H⌢. Then, from (8), we can get(12)P⌢dTF⌢rPw=0.Using (10) into (12), we can get(13)PdT10-du001-dv0001TF⌢rPw=0.Besides, from (8) and (13), we can get(14)Fr=10-du001-dv0001TF⌢r.F ⌢ r can be calculated as(15)F⌢r=100010du0dv01Fr.Note thatF⌢r=e⌢×H⌢ and e⌢=001T; we have(16)F⌢r=e⌢×H⌢=0-10100000H⌢.LetF⌢r=[f⌢1T;f⌢2T;f⌢3T] and H⌢=[h⌢1T;h⌢2T;h⌢3T]. Here, f⌢iT=f⌢i1f⌢i2f⌢i3 is the ith row of the radial fundamental matrix F⌢r, and h⌢iT=h⌢i1h⌢i2h⌢i3 is the ith row of matrix H⌢. Due to the fact that final row of F⌢r is zero, so we only need to solve the first two rows of H⌢ from (16), that is,(17)h⌢1T=f⌢2T,h⌢2T=-f⌢1T.Then by calculating the homographyH⌢, present up to only three unknown parameters h⌢3T=h⌢31h⌢32h⌢33. As the homography H⌢ relates the model point Pw and its undistorted image point P⌢u,(18)sP⌢u=H⌢Pw.Note that P⌢u=X⌢uY⌢u1 is the undistorted image point after the origin of the image coordinates is transferred to the COD. Considering the division model (DM), we have(19)X⌢u=X⌢dLr⌢d,k,Y⌢u=Y⌢dLr⌢d,k,where L(r⌢d,k) is the same as in (4), that is,(20)Lr⌢d,k=1+k1r⌢d2+k2r⌢d4+⋯.We substitute (19) into (18); then we can get(21)
s
X
⌢
d
L
r
⌢
d
,
k
Y
⌢
d
L
r
⌢
d
,
k
1
=
h
⌢
1
T
h
⌢
2
T
h
⌢
3
T
P
w
.So each point can be given the two equations,(22)X⌢dh⌢3TPw-h⌢1TPw·k1r⌢d2+k2r⌢d4+⋯Y⌢dh⌢3TPw-h⌢2TPw·k1r⌢d2+k2r⌢d4+⋯=h⌢1TPwh⌢2TPw.Using (17) in (22), we can get(23)X⌢dPwT-f⌢2TPwr⌢d2r⌢d4⋯Y⌢dPwTf⌢1TPwr⌢d2r⌢d4⋯h⌢3k1k2⋮=f⌢2TPw-f⌢1TPw.Given N corresponding point and that it satisfies that 2N≥n+3, where n is the number of distortion parameters, we can solve the system of equations in the least squares sense. Then the third row of H⌢ (i.e., h⌢3T) is solved. Moreover, the distortion coefficients k1,k2,… are computed at the same time.


## 3.3. Linear Solution of Intrinsic Parameters of the Camera

AfterH⌢ is obtained, the homography matrix H for the original image coordinates can be computed very easily. From (3), (11), and (18), we can get(24)H=10-du001-dv0001-1H⌢=10du001dv0001H⌢.Let us denote homography matrix byH=hc1hc2hc3 (where hci is the ith column of H). According to the orthogonality of the rotation matrix R(r1Tr2=0,r1Tr1=r2Tr2), we get(25)hc1TK-TK-1hc2=0,hc1TK-TK-1hc1=hc2TK-TK-1hc2.Note thatK is the camera intrinsic matrix. 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. Therefore, we need at least three images to calculate all the parameters. After the homography H is obtained, the analytical solution of the intrinsic parameters of the camera can be solved in the literature [11].


## 3.4. Maximum Likelihood Estimation

The intrinsic parameters obtained in the previous section are not accurate enough due to image noise. We can refine it through maximum likelihood estimation. Ahmed and Farag [13] have proved that including the distortion center and the decentering coefficients in the nonlinear optimization step may lead to instability of the estimation algorithm. Therefore, we do not need to optimize the distortion coefficient and distortion center. Meanwhile, it can reduce the search space of the calibration problem without sacrificing the accuracy and produce more stable and noise-robust results. Here, we give 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. Then the maximum likelihood estimate can be obtained by minimizing the following function:(26)∑i=1n∑j=1mPij-P⌢K,Ri,ti,Pj2,where P⌢(K,Ri,ti,Pj) is the projection of point Pj in image i according to (3). 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. Minimizing (26) is a nonlinear minimization problem, which is solved by the Levenberg-Marquardt algorithm [27]. It requires an initial value of K and {Ri,ti∣i=1,…,n}, which can be obtained using the technique described in the previous subsection.


## 4. Experimental Results

In this section, the proposed calibration method has been tested on both computer simulated data and real data. Experiments mainly consider two calibration methods, Zhang’s calibration method based on planar target [28] and our method (before-optimization and after-optimization).


### 4.1. Simulated Experiments

The simulated camera has the following property:fu=850, fv=850, u0=512, v0=384. The image size is 1024 × 768 pixels. The skew factor is set to zero. Besides, a second-order radial distortion is simulated with the coefficients k1=-6.09×10-7 pixel-2, k2=-1.97×10-13 pixel-4 and the COD is set to (500,366). The plane is a checkerboard image with 70 corners (7 × 10) evenly distributed and the minimum point interval is set to 23 mm, which is the same as the real data experiments. The direction of the plane is denoted by a 3-dimensional vector r, which is parallel to the rotation axis and whose magnitude is the same as the angle of rotation. In addition, the position of the plane is denoted by a 3-dimensional vector t. We employ four planes with r1=[20∘,0∘,0∘]T,  t1=[-80,-60,200]T,  r2=[0∘,0∘,20∘]T,  t2=[-110,-80,250]T, r3=[-40∘,0∘,20∘]T, t3=[-100,-40,330]T, r4=[-10∘,0∘,20∘]T, t4=[-100,-60,280]T in the experiment. Figure 2 presents the simulation of the images, which show the true projection of the image point and the distorted point. The Gauss noise with a mean of 0 and a standard deviation of σ is added to the distorted image.Figure 2
The simulated images with the true projection image points and the distorted points.Then the estimated camera parameters and the actual value of the real were compared. We calculated the relative error forfu and fv and the absolute errors for u0 and v0. Besides, the mean value of the reprojection error Erms is also used as an evaluation index. It is calculated by the discrepancy between the ground true point Pt,i and the reprojection image point Pr,i, (27)Erms=∑i=1nPr,i-Pt,i2n.Noise is added to the projected image points withσ ranging from 0.1 pixels to 2.0 pixels. For each noise level, 50 independent repeated trials were performed, and the results were shown to be average. The relative error and calibration precision of the camera are shown in Figures 3 and 4, respectively. From the two figures, we can see that the errors increase linearly with the noise level. It is worth noting that the error of the two methods is relatively low when the noise is low. For example, for σ=0.5 (it is larger than the normal noise in practical calibration), the absolute errors in u0 and v0 are around 1 pixel, the relative errors in fu, fv are less than 0.3%, and the mean calibration errors are around 0.2 pixels. Besides, the error in v0 is larger than that in u0. The main reason is that there are less data in the v0 direction than in the u0 direction.Figure 3
Effects of pixel coordinates noise on intrinsic parameters using the proposed method (after-optimization: (a) the relative error forfu and fv and (b) the absolute error for u0 and v0; before-optimization: (c) the relative error for fu and fv and (d) the absolute error for u0 and v0).
(a) (b) (c) (d)Figure 4
Effects of pixel coordinates noise on calibration accuracy: (a) after-optimization; (b) before-optimization.
(a) (b)


### 4.2. Real Data Experiments

For the real data experiments, the calibration images provided by The Robotics Institute of Carnegie Mellon University [28] are used to test our approach. The calibration template is a planar checkerboard pattern with 70 corners (7 × 10) evenly distributed and the minimum point interval is 23 mm in both the vertical and the horizontal directions. Ten images of the plane under different orientations were taken, as shown in Figure 5. (It is worth noting that the algorithm proposed in this paper is suitable for the calibration of all kinds of image of the board, including all the pictures provided by Carnegie Mellon University. Due to the limitation of space, this paper only randomly selected ten pictures as our experimental object.) We can observe an obvious lens distortion in the images, particularly in Image 1, Image 3, Image 5, and Image 6. In addition, the resolution of the image is 1024 × 768 pixels.Figure 5
Images used for calibration.First, we use the corner detection method to get subpixel precision of checkerboard angular point position. Next, we compare our proposed calibration method with Zhang’s method. In the experiment, we utilize the first three images to obtain the camera parameters, which are displayed in Table1. In Zhang’s method, he does not consider the distortion center (du0,dv0), while our method does. It is necessary to explain that (du0,dv0) are the mean of the first three images. In fact, the COD of the three images are 507.368, 508.087, and 506.139, respectively. In addition, our approach adopts the DM to represent lens distortion, and Zhang’s method depends on the PM, so the generating distortion coefficients are different. At last, it is worth noting that the before-optimization method also computes the COD. However, the resulting fu,fv,u0, and v0 are almost the same as Zhang’s. We will see later that the accuracy of the before-optimization is not high.Table 1
Comparative result of intrinsic parameters and distortion coefficients.
Method f u f v u 0 v 0 Zhang 839.50 838.56 507.89 367.02 Before-optimization 839.79 838.87 506.17 367.14 After-optimization 838.23 837.27 505.99 366.21  Method k 1 k 2 d u 0 d v 0  Zhang −0.42 mm-2 0.20 mm-4 — — Before-optimization - 6.09 e - 07 pixel-2 - 1.97 e - 13 pixel-4 507.20 367.59 After-optimization - 6.09 e - 7 pixel-2 - 1.97 e - 13 pixel-4 507.20 367.59To investigate the effectiveness of the proposed method, the other seven test images are used to evaluate the calibration accuracy. Here, we once again put the calculated reprojection errorEr-rms as the evaluation index. Unlike the simulation test, the ground true image points are unknown here, so we use the undistorted point Pud,i instead. The reprojection error Er-rms is defined as(28)Er-rms=∑i=1nPrp,i-Pud,i2n.Table2 shows the comparison of the reprojection error of the test data, and the distribution is shown in Figure 6. As shown in Table 2, in both the average accuracy level and standard error, the after-optimization method is better than Zhang’s method and the before-optimization method. From Figure 6, we can see that the distribution of the after-optimization method is more concentrated in the near zero.Table 2
Comparative result of calibration errorEr-rms evaluated by testing data.
Method Image 4 Image 5 Image 6 Image 7 Image 8 Zhang 0.1683 0.1763 0.1319 0.1380 0.1658 Before-optimization 0.1613 0.1716 0.1389 0.1421 0.1746 After-optimization 0.1534 0.1670 0.1374 0.1310 0.1609  Method Image 9 Image 10 Average Standard error —  Zhang 0.1321 0.1580 0.1529 0.0186 — Before-optimization 0.1347 0.1372 0.1515 0.0172 — After-optimization 0.1212 0.1383 0.1442 0.0167 —Figure 6
Reprojection error: (a) Zhang’s method; (b) proposed method (after-optimization). Each point marked as “+” in the figure denotes the reprojection error of corresponding chessboard corner in the testing images.
(a) (b)Besides, in order to further study the stability of the proposed (after-optimization) algorithm, we have also applied it to some combinations of 4 images from the fifth to the ninth images. The results are shown in Table3, where the second column (5678), for example, displays the result with the quadruple of the fifth, sixth, seventh, and eighth image. The last two columns display the mean and deviations. The deviations for all parameters are very small, which means that the proposed algorithm is quite stable.Table 3
Variation of the calibration results among all quadruples of images.
Quadruples (5678) (5689) (5789) (6789) f u 843.39 843.26 843.31 845.70 f v 842.60 842.44 842.65 844.71 u 0 504.50 504.79 504.73 504.92 v 0 365.68 365.48 365.75 366.51 k 1 - 6.0822 e - 07 - 6.1079 e - 07 - 6.0732 e - 07 - 6.0985 e - 07 k 2 - 2.2510 e - 13 - 2.0867 e - 13 - 2.2243 e - 13 - 1.9719 e - 13  Quadruples (5679) Mean Var —  f u 843.07 843.7460 1.2070 — f v 842.23 842.9260 1.0214 — u 0 504.53 504.6940 0.0315 — v 0 365.41 365.7660 0.1924 — k 1 - 6.1882 e - 07 - 6.1100 e - 07 2.0947 e - 17 — k 2 - 1.2948 e - 13 - 1.9657 e - 13 1.5326 e - 27 —Finally, for a few images with a larger distortion, we conducted an experiment separately. We utilize four images (Image 1, Image 3, Image 5, and Image 6) to calibrate the internal parameters. The results are shown in Table4. As can be seen from the results, we obtained similar results compared to Zhang.Table 4
Comparative result of intrinsic parameters and distortion coefficients.
Method f u f v u 0 v 0 Zhang 846.95 845.92 505.66 365.81 Before-optimization 842.57 841.01 504.97 365.33 After-optimization 846.03 844.93 505.81 366.24  Method k 1 k 2 d u 0 d v 0  Zhang −0.43 mm-2 0.21 mm-4 — — Before-optimization - 6.12 e - 07 pixel-2 - 2.01 e - 13 pixel-4 506.40 366.92 After-optimization - 6.12 e - 07 pixel-2 - 2.01 e - 13 pixel-4 506.40 366.92


## 4.1. Simulated Experiments

The simulated camera has the following property:fu=850, fv=850, u0=512, v0=384. The image size is 1024 × 768 pixels. The skew factor is set to zero. Besides, a second-order radial distortion is simulated with the coefficients k1=-6.09×10-7 pixel-2, k2=-1.97×10-13 pixel-4 and the COD is set to (500,366). The plane is a checkerboard image with 70 corners (7 × 10) evenly distributed and the minimum point interval is set to 23 mm, which is the same as the real data experiments. The direction of the plane is denoted by a 3-dimensional vector r, which is parallel to the rotation axis and whose magnitude is the same as the angle of rotation. In addition, the position of the plane is denoted by a 3-dimensional vector t. We employ four planes with r1=[20∘,0∘,0∘]T,  t1=[-80,-60,200]T,  r2=[0∘,0∘,20∘]T,  t2=[-110,-80,250]T, r3=[-40∘,0∘,20∘]T, t3=[-100,-40,330]T, r4=[-10∘,0∘,20∘]T, t4=[-100,-60,280]T in the experiment. Figure 2 presents the simulation of the images, which show the true projection of the image point and the distorted point. The Gauss noise with a mean of 0 and a standard deviation of σ is added to the distorted image.Figure 2
The simulated images with the true projection image points and the distorted points.Then the estimated camera parameters and the actual value of the real were compared. We calculated the relative error forfu and fv and the absolute errors for u0 and v0. Besides, the mean value of the reprojection error Erms is also used as an evaluation index. It is calculated by the discrepancy between the ground true point Pt,i and the reprojection image point Pr,i, (27)Erms=∑i=1nPr,i-Pt,i2n.Noise is added to the projected image points withσ ranging from 0.1 pixels to 2.0 pixels. For each noise level, 50 independent repeated trials were performed, and the results were shown to be average. The relative error and calibration precision of the camera are shown in Figures 3 and 4, respectively. From the two figures, we can see that the errors increase linearly with the noise level. It is worth noting that the error of the two methods is relatively low when the noise is low. For example, for σ=0.5 (it is larger than the normal noise in practical calibration), the absolute errors in u0 and v0 are around 1 pixel, the relative errors in fu, fv are less than 0.3%, and the mean calibration errors are around 0.2 pixels. Besides, the error in v0 is larger than that in u0. The main reason is that there are less data in the v0 direction than in the u0 direction.Figure 3
Effects of pixel coordinates noise on intrinsic parameters using the proposed method (after-optimization: (a) the relative error forfu and fv and (b) the absolute error for u0 and v0; before-optimization: (c) the relative error for fu and fv and (d) the absolute error for u0 and v0).
(a) (b) (c) (d)Figure 4
Effects of pixel coordinates noise on calibration accuracy: (a) after-optimization; (b) before-optimization.
(a) (b)


## 4.2. Real Data Experiments

For the real data experiments, the calibration images provided by The Robotics Institute of Carnegie Mellon University [28] are used to test our approach. The calibration template is a planar checkerboard pattern with 70 corners (7 × 10) evenly distributed and the minimum point interval is 23 mm in both the vertical and the horizontal directions. Ten images of the plane under different orientations were taken, as shown in Figure 5. (It is worth noting that the algorithm proposed in this paper is suitable for the calibration of all kinds of image of the board, including all the pictures provided by Carnegie Mellon University. Due to the limitation of space, this paper only randomly selected ten pictures as our experimental object.) We can observe an obvious lens distortion in the images, particularly in Image 1, Image 3, Image 5, and Image 6. In addition, the resolution of the image is 1024 × 768 pixels.Figure 5
Images used for calibration.First, we use the corner detection method to get subpixel precision of checkerboard angular point position. Next, we compare our proposed calibration method with Zhang’s method. In the experiment, we utilize the first three images to obtain the camera parameters, which are displayed in Table1. In Zhang’s method, he does not consider the distortion center (du0,dv0), while our method does. It is necessary to explain that (du0,dv0) are the mean of the first three images. In fact, the COD of the three images are 507.368, 508.087, and 506.139, respectively. In addition, our approach adopts the DM to represent lens distortion, and Zhang’s method depends on the PM, so the generating distortion coefficients are different. At last, it is worth noting that the before-optimization method also computes the COD. However, the resulting fu,fv,u0, and v0 are almost the same as Zhang’s. We will see later that the accuracy of the before-optimization is not high.Table 1
Comparative result of intrinsic parameters and distortion coefficients.
Method f u f v u 0 v 0 Zhang 839.50 838.56 507.89 367.02 Before-optimization 839.79 838.87 506.17 367.14 After-optimization 838.23 837.27 505.99 366.21  Method k 1 k 2 d u 0 d v 0  Zhang −0.42 mm-2 0.20 mm-4 — — Before-optimization - 6.09 e - 07 pixel-2 - 1.97 e - 13 pixel-4 507.20 367.59 After-optimization - 6.09 e - 7 pixel-2 - 1.97 e - 13 pixel-4 507.20 367.59To investigate the effectiveness of the proposed method, the other seven test images are used to evaluate the calibration accuracy. Here, we once again put the calculated reprojection errorEr-rms as the evaluation index. Unlike the simulation test, the ground true image points are unknown here, so we use the undistorted point Pud,i instead. The reprojection error Er-rms is defined as(28)Er-rms=∑i=1nPrp,i-Pud,i2n.Table2 shows the comparison of the reprojection error of the test data, and the distribution is shown in Figure 6. As shown in Table 2, in both the average accuracy level and standard error, the after-optimization method is better than Zhang’s method and the before-optimization method. From Figure 6, we can see that the distribution of the after-optimization method is more concentrated in the near zero.Table 2
Comparative result of calibration errorEr-rms evaluated by testing data.
Method Image 4 Image 5 Image 6 Image 7 Image 8 Zhang 0.1683 0.1763 0.1319 0.1380 0.1658 Before-optimization 0.1613 0.1716 0.1389 0.1421 0.1746 After-optimization 0.1534 0.1670 0.1374 0.1310 0.1609  Method Image 9 Image 10 Average Standard error —  Zhang 0.1321 0.1580 0.1529 0.0186 — Before-optimization 0.1347 0.1372 0.1515 0.0172 — After-optimization 0.1212 0.1383 0.1442 0.0167 —Figure 6
Reprojection error: (a) Zhang’s method; (b) proposed method (after-optimization). Each point marked as “+” in the figure denotes the reprojection error of corresponding chessboard corner in the testing images.
(a) (b)Besides, in order to further study the stability of the proposed (after-optimization) algorithm, we have also applied it to some combinations of 4 images from the fifth to the ninth images. The results are shown in Table3, where the second column (5678), for example, displays the result with the quadruple of the fifth, sixth, seventh, and eighth image. The last two columns display the mean and deviations. The deviations for all parameters are very small, which means that the proposed algorithm is quite stable.Table 3
Variation of the calibration results among all quadruples of images.
Quadruples (5678) (5689) (5789) (6789) f u 843.39 843.26 843.31 845.70 f v 842.60 842.44 842.65 844.71 u 0 504.50 504.79 504.73 504.92 v 0 365.68 365.48 365.75 366.51 k 1 - 6.0822 e - 07 - 6.1079 e - 07 - 6.0732 e - 07 - 6.0985 e - 07 k 2 - 2.2510 e - 13 - 2.0867 e - 13 - 2.2243 e - 13 - 1.9719 e - 13  Quadruples (5679) Mean Var —  f u 843.07 843.7460 1.2070 — f v 842.23 842.9260 1.0214 — u 0 504.53 504.6940 0.0315 — v 0 365.41 365.7660 0.1924 — k 1 - 6.1882 e - 07 - 6.1100 e - 07 2.0947 e - 17 — k 2 - 1.2948 e - 13 - 1.9657 e - 13 1.5326 e - 27 —Finally, for a few images with a larger distortion, we conducted an experiment separately. We utilize four images (Image 1, Image 3, Image 5, and Image 6) to calibrate the internal parameters. The results are shown in Table4. As can be seen from the results, we obtained similar results compared to Zhang.Table 4
Comparative result of intrinsic parameters and distortion coefficients.
Method f u f v u 0 v 0 Zhang 846.95 845.92 505.66 365.81 Before-optimization 842.57 841.01 504.97 365.33 After-optimization 846.03 844.93 505.81 366.24  Method k 1 k 2 d u 0 d v 0  Zhang −0.43 mm-2 0.21 mm-4 — — Before-optimization - 6.12 e - 07 pixel-2 - 2.01 e - 13 pixel-4 506.40 366.92 After-optimization - 6.12 e - 07 pixel-2 - 2.01 e - 13 pixel-4 506.40 366.92


## 5. Conclusion

In this paper, we proposed a flexible and high-accuracy camera calibration method. Compared with the traditional method, this method overcomes a lot of problems. Firstly, it decouples the estimation of the distortion coefficients and intrinsic parameters, producing more stable and reliable results. Then, the distortion of the center is accurately estimated and it is important to obtain the best results. Moreover, whether it is for a relatively small image distortion or distortion larger image, utilizing our method can get a good result. Finally, the robustness and accuracy of the proposed method are verified by simulation and real data experiments and the experimental results show that this method has the advantages of simple operation, high accuracy, and better flexibility.



---

*Source: 1392832-2016-12-29.xml*