To perform the calibration, the calibration points--internal corners on the checkerboard--must be found and arranged into a grid structure resembling the calibration target. The user should have to merely provide the number of rows and columns of the target. This process proceeds in several steps: possible point identification and refinement, reduction to a probable set, and the arrangement of the points into a grid structure with missing point adaptions.

scorePoint(Point P){
  Let N be a 13x13 window of an image 
  with a center n66  ∈ N 
  such that n66=P

  Let M1,M2 be the matrices
  with entries specified below. 
  Let m1ij ∈ M1, m2ij ∈ M2

  score= MAX(
   (∑i∑j if(m1ij!=0){ | (nij-128)-127*m1ij | }),
   (∑i∑j if(m2ij!=0){ | (nij-128)-127*m2ij | })
  )
  return score
}

The above method produces a set of points that highly correlates to the calibration points. However, the position of the points need to be refined to sub-pixel accuracy and spurious points need to be eliminated. In order to accomplish this, the neighborhood of the point can be thought of as a surface with color denoting height. If the neighborhood contains a calibration point, the surface will be the form of a saddle point. If a quadratic surface is fit to the neighborhood, the critical point of the surface will be the location of the calibration point. In order to achieve better accuracy, the neighborhood is recentered on the updated point, the quadratic surface is refitted, and the critical point is reevaluated. If the position of the point converges, the point is kept as a possible calibration point; if it diverges the point is removed from consideration. Once all the possible points are refined, the set should be checked to remove any double entries.

Let S be a set of Points
Let n be the desired number of points

While( |S|>n ){
  minAvg=(∑(p ∈ S)dist(p,p'))/|S|
     where p' ∈ S | dist(p,p')<=min(t ∈ S|t!=p)dist(p,t)

  for each point p ∈ S define the score as
    s1=∑(t ∈ S)dist(p,t)
    p' ∈ S| dist(p,p')<=min(t ∈ S|t!=p)dist(p,t)
    p'' ∈ S| dist(p,p'')<=min(t ∈ S|t!=p, t!=p')dist(p,t)
    s2= |dist(p,p')-minAvg|/minAvg +1
    s3= |dist(p,p'')-minAvg|/minAvg +1
    score=s1*s2*s3
  eliminate the point with highest score from S
}
    

The set of remaining points needs to be further culled of outliers which may still be present. The set should be reduced to the number of calibration points present on the target. The possibilities in choosing K points from a set of N can be a magnitude of the factorial of N, and therefore an optimal solution is unlikely to be tractable. Instead, as in a method akin to hill climbing, the least fit points will be eliminated one at a time until the desired number of points is reached. To evaluate the points for fitness, the characteristics of the calibration target must be considered. The set of points that corresponds to calibration points will be fairly compact and for each point in the set, the neighboring points of the set will be nearly uniform. For each pair of points the distance is calculated and for each point the distance to the two nearest points is indexed. To measure compactness, the sum of the distance from one point to all others is calculated. To measure uniformity, a more complex measure is used. First the absolute value of the difference between the distance of a neighbor and the average distance to nearest neighbors is computed. The measure is the sum of that value divided by the average distance to nearest neighbors and one. The measure is computed for the two nearest neighbors of a point. The point to be eliminated is the one with the highest product of the three measures. After a point is eliminated, the distances and averages are updated for the set.

stepList=NULL for each pair of points a,b ∈ S a, b are neighbors if !∃ c ∈ S such that dist(a,b) < dist(a,c) && dist(a,b) < dist(b,c) if a,b are neighbors rise=a.y-b.y run=a.x-b.x if rise < 0 rise= -rise run= -run stepList.add( (rise,run) ) stepList.sort_rise_then_run() horiz=stepList[ |stepList|/4 ] vert=stepList[ |stepList|* 3/4 ] find a Point P ∈ S such that for a given error threshold ε |δx| < ε P+horiz+δ1 ∈ S P-horiz+δ2 ∈ S P+vert+δ3 ∈ S P-vert+δ4 ∈ S P+horiz+vert+δ5 ∈ S P+horiz-vert+δ6 ∈ S P-horiz+vert+δ7 ∈ S P-horiz-vert+δ8 ∈ S else fail Start a grid G such that P ∈ G 8 neighboring points of P ∈ G A calibration target has a number of rows and columns Let the lower bound lb=min( #rows, #columns) Let the upper bound ub=max( #rows, #columns) while G.rowSize < lb && G.colSize < lb G.addRow G.addCol if G can add more rows while G.colSize < ub G.addRow else G can add more columns while G.rowSize < ub G.addCol

In order to perform calibration, the set of calibration points must be arranged so that each point is one-to-one and onto with an abstract calibration target model. In order to fit the points to a grid, characteristics of the arrangement must be recovered. The grid will be roughly rectangular (distortion and skew affect the projection) and as such there should be a general distance between rows and columns. For lack of better terms, these distances will be referred to as the vertical and horizontal steps respectively. In order to recover these general directions, the set of points is analyzed for distances between neighboring points. In this context, two points are considered neighboring if the distance between the points is less than a distance between one of the neighboring points and all other points in the set. These distances--recorded as a change in rise and run (y and x) manipulated so that the run is always positive--are sorted in an array so that the run is ascending. The general horizontal and vertical distances will then be at the first and third quartiles of the array (the entries at the indices 1/4 and 3/4 of total magnitude of the array). Once the general directions have been recovered, the set of points is searched to fill a 3x3 block in which entries are separated the horizontal and vertical steps. After the initial 3x3 block is found, rows and columns are added one at a time until the desired dimensions are reached. Proposed points for the side to be added are determined by the nearest point, the difference of the nearest point and its parent, and the rate of change of the difference. Proposed points are then checked for matches within the set of points or are refined via the quadratic surface fitting defined in the refinement section.

Calibration is performed as described by Zhengyou Zhang. Since the methodology is already written in his paper, I will only provide a brief overview as well as note some places of numerical importance.

Before performing calibration, the model and image points need to be normalized for numerical stability. Model points should be offset and scaled such that the centroid of the model is at the origin and that the average magnitude of the coordinates is 1. Assuming a nearly square picture, the image points need to be offset so that the center of the image is at the origin and assuming a nearly square picture the points scaled by the half the smallest side.

It should be noted, while (B^-1/s) should be theoretically be a positive definite matrix and thereby would enable the extraction of A via Cholesky decomposition, in practice, due to finite precision algebra and least squares approximations, this is not necessarily the case. Thus in extracting A, one may have to take liberties with some of the values to avoid complex numbers. Since these values will be optimized in subsequent steps, errors in calculating A can be mitigated.

At this point, the model and image points need to be restored to their unnormalized values. Additionally, the camera parameters and translation vectors need to be updated to reflect the change in coordinates.

   [ alpha/s  skew/s  u0/s-x_offset ]
A`=[ 0        beta/s  v0/s-y_offset ]
   [ 0        0       1             ]

Once the points have been restored from their normalized values and the camera parameters and translation vectors have been adjusted, all the variables need to be optimized via the LM algorithm. If n images are used, there are 5+6n variables (the 5 camera parameters, and for each image three variables representing the translation vector and three variables representing the rotation matrix encoded according to the Rodriguez formula/angle and axis representation).

Once the values are optimized, the variables need to be optimized again with barrel distortion taken into account. If n images are used, there are then 7+6n variables.