A Decoupled Calibration Method for Camera Intrinsic Parameters and Distortion Coefficients

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 divisionmodel, 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.


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][5][6][7][8][9][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 2 Mathematical Problems in Engineering to real sensor array plane.Ahmed and Farag [13] proposed a robust approach to distortion calibration by using a leastmedian-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][16][17][18][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.Section 2 gives a brief description of the camera model and the lens distortion model.Section 3 describes the detail procedure of the Figure 1: The model of camera.
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.

Camera Model
The camera model is shown in Figure 1.Given one homogeneous coordinates point P  = (  ,   ,   , 1)  in 3D space.P  = (  ,   , 1)  is the perspective projection of P  in the image plane based on pinhole model without lens distortion.P  = (  ,   , 1)  is the real projection of P  in the image plane of camera considering the lens distortion.Without consideration of lens distortion, the mapping between the 3D point P  and 2D image point P  is given by where  represents the nonzero scale factors; K is the camera intrinsic matrix, with ( 0 , V 0 ) being the coordinates of the principal point,   and  V being the effective focal length in pixels, and   being the parameter describing the skew of the two image axes; [R t] 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 plane   = 0 with no loss of generality.From (1), we have where r  is the th column of the rotation matrix R. By abuse of notation, we still denote the point on the model plane by P  = (  ,   , 1)  .Therefore, 3D point P  and its image point P  are related by a homography H, With the influence of lens distortion, as shown in Figure 1, the actual image point P  is not the point P  , which is the image plane intersection of the connection of the 3D point P  and the optical center O  , 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: where  1 ,  2 , . . .are the distortion coefficients, e = ( 0 , V 0 , 1)  is the homogeneous coordinate of the center of distortion (COD), and   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, where (  , ) is the same as in (4).Combining ( 4) and ( 5) we can get a more generic Rational Model (RM): 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.

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.

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 P  on the calibration model plane and its corresponding image point P  are related by a radial fundamental relationship, and it can be represented as follows: where H is the same as in (3) and [e] × (the skew symmetric matrix) represents the cross-product.Writing F  = [e] × H (in [24], the matrix F  is called as radial fundamental matrix), then (7) can be rewritten as Note that F  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: Of course, in the absence of radial deformation, the basic matrix of the above calculation is not stable, and the estimated value of e 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].

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 COD e, we can compute the H by factoring the matrix F  (F  = [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 Similarly, the transformed undistorted image point  P  and the original undistorted image point P  are described as When the coordinate origin is converted to the center of distortion, the homography matrix is redefined as  H and the corresponding radial fundamental matrix is defined as H.Then, from (8), we can get Using ( 10) into (12), we can get Besides, from ( 8) and ( 13), we can get F  can be calculated as Note that  F  = [ e] ×  H and  e = [0 0 1]  ; we have Let is the th row of the radial fundamental matrix  F  , and is the th row of matrix  H. Due to the fact that final row of  F  is zero, so we only need to solve the first two rows of  H from (16), that is, Then where (   , ) is the same as in (4), that is, We substitute ( 19) into (18); then we can get So each point can be given the two equations, Using ( 17) in (22), we can get Given  corresponding point and that it satisfies that 2 ≥  + 3, where  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  3 ) is solved.Moreover, the distortion coefficients { 1 ,  2 , . ..} are computed at the same time.

Linear Solution of Intrinsic Parameters of the Camera.
After  H is obtained, the homography matrix H for the original image coordinates can be computed very easily.From (3), (11), and (18), we can get Let us denote homography matrix by H = [h 1 h 2 h 3 ] (where h  is the th column of H).According to the orthogonality of the rotation matrix R (r  1 r 2 = 0, r  1 r 1 = r  2 r 2 ), we get Note that K 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].

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  images of a model plane and there are  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: where  P(K, R  , t  , P  ) is the projection of point P  in image  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 {R  ,   |  = 1, . . ., }, which can be obtained using the technique described in the previous subsection.

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).

Simulated Experiments.
The simulated camera has the following property:   = 850,  V = 850,  0 = 512, V 0 = 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  1 = −6.09× 10 −7 pixel-2,  2 = −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 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.
Then the estimated camera parameters and the actual value of the real were compared.We calculated the relative error for   and  V and the absolute errors for  0 and V 0 .Besides, the mean value of the reprojection error  rms is also used as an evaluation index.It is calculated by the discrepancy between the ground true point P , and the reprojection image point P , , 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  0 and V 0 are around 1 pixel, the relative errors in   ,  V are less than 0.3%, and the mean calibration errors are around 0.2 pixels.Besides, the error in V 0 is larger than that in  0 .The main reason is that there are less data in the V 0 direction than in the  0 direction.

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.
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 Table 1.In Zhang's method, he does not consider the distortion center ( 0 , V 0 ), while our method does.It is necessary to explain that ( 0 , V 0 ) 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   ,  V ,  0 , and V 0 are almost the same as Zhang's.We will see later that the accuracy of the before-optimization is not high.
To investigate the effectiveness of the proposed method, the other seven test images are used to evaluate the calibration Table 2 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 afteroptimization method is more concentrated in the near zero.
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 Table 3, 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.
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 Table 4.As can be seen from the results, we obtained similar results compared to Zhang.

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  Relative error (%) Noise level (pixels) Noise level (pixels)      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.
, the image plane coordinate system O   is transformed to O        .Where O  is the real position of the COD, we represent it as  e = [0 0 1]  in homogeneous coordinates.Then the transformed distorted image point  P  can be represented as  P

Figure 2 :
Figure 2: The simulated images with the true projection image points and the distorted points.

Figure 3 :
Figure3: Effects of pixel coordinates noise on intrinsic parameters using the proposed method (after-optimization: (a) the relative error for   and  V and (b) the absolute error for  0 and V 0 ; before-optimization: (c) the relative error for   and  V and (d) the absolute error for  0 and V 0 ).

Table 1 :
Comparative result of intrinsic parameters and distortion coefficients.Here, we once again put the calculated reprojection error  -rms as the evaluation index.Unlike the simulation test, the ground true image points are unknown here, so we use the undistorted point P , instead.The reprojection error  -rms is defined as =1 √      P , − P ,      2 ) .

Table 2 :
Comparative result of calibration error  -rms evaluated by testing data.