Camera calibration

1What calibration is

The note on 3D transformations ended with a model of the camera:

It takes a point M′ = [X′ Y′ Z′ 1]^T, given in the coordinates of some object, to its pixel w = [u v 1]^T. It has eleven parameters. Five are intrinsic, α, β, γ, u₀ and v₀, in K and P. Six are extrinsic, the angles ω_X, ω_Y, ω_Z and the three components of t′, in S′. Calibration is estimating them from images.

The idea is to photograph an object whose shape is known exactly, a calibration object, and pick N points on it. Their coordinates M′_i are known from the object, in millimeters, and their pixels w_i are measured in the image. Collect the parameters in a vector θ. For any guess of θ, the model predicts where each point should appear:

where F builds K, P and S′ from θ, multiplies, and divides by λ. A good θ puts every ŵ_i close to the measured w_i, so calibration looks for the θ that makes the average distance between them smallest:

J is measured in pixels, so it says directly how far off the model is, on average.

This is called photogrammetric calibration, since it relies on an object of known size. Self-calibration does without one: it only relates the same points seen in several images. But it cannot tell the true size of the scene, since a small object near the camera gives the same images as a larger copy farther away. When the result must be in millimeters, as for a robot that picks up parts, a calibration object is used.

The table

The example throughout this page is a table, 800 mm deep, 1250 mm wide and 900 mm tall, photographed from the front in an image of 1842 × 1408 pixels. Its eight known points are the top and bottom ends of its four legs, numbered 0 to 7. The table's coordinate system has its origin at the bottom of the back left leg, corner 3, with X′ pointing from the back legs toward the front legs and the camera, Y′ to the right along the table, and Z′ up, a right-handed system. So the bottom of the front left leg, corner 0, is at (800, 0, 0), and the top of the front right leg, corner 5, is at (800, 1250, 900). The pixels were clicked by hand on the photo, so each is off by a few pixels.

2Linear calibration

Multiplied out, KPS′ is a single 3 × 4 matrix. Call it A, the general projection matrix (in code it is usually called P):

Instead of the eleven parameters, the linear method estimates the twelve entries of A. They are not all free: A and kA give the same pixels for any number k ≠ 0, since k only changes λ. So one entry can be fixed, and the usual choice is a₃₄ = 1. That leaves eleven unknowns, as many as the camera has parameters.

Why may a₃₄ be fixed? The third row of A gives λ. For the object's origin, M′ = [0 0 0 1]^T, it gives λ = a₃₄, so a₃₄ is the depth of that origin in front of the camera, times the scale of A. Setting it to 1 just picks the scale. The one exception is a true a₃₄ of 0, when the object's origin is level with the optical center, on the plane through C parallel to the image. No scale turns 0 into 1. An origin somewhere in view, as on the table, avoids this.

Two equations per point

The third row gives λ = a₃₁X′ + a₃₂Y′ + a₃₃Z′ + 1. Putting it into the first two rows removes λ:

Both equations are linear in the unknowns, because u, v, X′, Y′ and Z′ are known numbers. Stack the unknowns in a = [a₁₁ a₁₂ … a₃₃]^T, eleven entries row by row, and move them to the left:

Each point adds these two rows, and N points give a linear system

with 2N equations, where Q has eleven columns and r lists the measured u and v.

How many points

Eleven unknowns need at least eleven equations, so at least six points: five give only ten. And the points must not all lie on one plane. If they did, say all on the floor, where Z′ = 0, then a₁₃, a₂₃ and a₃₃ would multiply 0 in every row, and nothing would fix them. The points of a plane only fix the homography between that plane and the image, as in the note on 2D transformations, and a homography has eight numbers, not eleven. The eight corners of the table lie on two levels, the floor and the tabletop, and they fix A.

Least squares

Clicked pixels are off by a few pixels, so with more than six points no a satisfies every equation. As for homographies, take the a that makes the sum of squares ‖Qa − r‖² smallest:

The columns of Q differ in size by millions: a 1 next to −uY′, which reaches two million here. Forming Q^TQ squares that spread and can lose digits, so numerical libraries, such as np.linalg.lstsq, reach the same â without forming it.

What does least squares minimize here? The first equation of point i leaves the residual (first row of A)M′_i − u_iλ_i. The first term is λ_iû_i, so the residual is λ_i(û_i − u_i): the error in pixels, times λ_i. With a₃₄ = 1, λ_i is the point's depth divided by the depth of the object's origin, so points farther from the camera weigh more. The result is close to the θ with the smallest J, but not the same, and this is one reason for the nonlinear method of section 4.

The choice a₃₄ = 1 can also be avoided, as for homographies. Keep all twelve entries as unknowns, add −u and −v as a twelfth column of Q, and solve Qa = 0 with ‖a‖ = 1. The answer is the right singular vector of Q with the smallest singular value. For the table, the two A agree to three significant digits, and their average errors differ by less than 0.001 pixels.

3Where the camera is, and its eleven parameters

A holds the whole camera, but its entries mix everything together. Two facts from the camera model take it apart. Write A = [B b], with B its left 3 × 3 block and b its last column.

The optical center

The optical center C has no image: every ray starts there, so A sends it to the zero vector. In the object's coordinates, AC = BC + b = 0, and so

This needs no parameters at all: it finds where the camera stood from A alone.

Peeling off K and R′

The model says that A is a multiple of KPS′. Since KP = [K 0], that product is K[R′ t′], and

for some number k. Two properties split B. K is upper triangular, with bottom row [0 0 1]. And the rows of R′, call them r₁, r₂ and r₃, are the camera's axes X, Y and Z written in the object's system: unit vectors, perpendicular to each other. Likewise, call m₁, m₂ and m₃ the rows of B. Row by row from the bottom, B = kKR′ says

and the unknowns come out one at a time.

• r₃ has length 1, so k = ±‖m₃‖ and r₃ = m₃/k.
• A dot product with r₃ picks out the part of m₂/k along r₃, since r₂ is perpendicular to it: v₀ = m₂·r₃/k. What is left, m₂/k − v₀r₃, is βr₂, whose length is β and whose direction is r₂.
• In the same way, u₀ = m₁·r₃/k and γ = m₁·r₂/k, and what is left of m₁/k is αr₁.

This is the Gram–Schmidt process run on the rows of B, from the bottom up. The sign of k is the one that puts the object in front of the camera, with λ positive for its points. Then t′ = K⁻¹b/k, and the angles come from R = R′^T, whose first row is

It gives ω_Y and ω_Z, and ω_X follows from the last column of R.

Lengths are positive, so α and β always come out positive, and a mirror in the data shows up in R′: if det R′ = −1, R′ turns a right hand into a left hand, and no real camera gives those pixels.

Which way the six numbers point

On this page, as in the note on 3D transformations, the camera's coordinate system is the old one and the object's is the new one. The angles ω turn the camera's axes into the object's, and t′ is where the camera sees the object's origin. Many programs describe the same camera the other way round, with the camera's system as the new one: a point M′ of the object is at R′M′ + t′ in the camera's system. Their translation is the same t′, but their rotation matrix is R′, the transpose of R, so their three angles are different numbers for the same camera. Before using angles from elsewhere, check which of the two matrices they build.

For the table, this matters in one more way. Its X′ axis points almost straight back at the camera, so ω_Y comes out near 90°, close to the gimbal lock of the note on 3D transformations: at 90°, ω_X and ω_Z would turn about the same line. The angles are still well defined at 79°, but some small changes are magnified. Moving corner 5 down by 10 pixels turns the camera by half a degree, yet changes ω_X by 2.5° and ω_Z by 2.3°, both in the same direction and about five times as much. This is one reason why calibration programs store a rotation as an axis and an angle instead.

4Nonlinear calibration

The linear method never looks at the parameters. It estimates twelve numbers and minimizes an error in the equations, not in the image. Nonlinear calibration works with the parameters themselves,

and minimizes J(θ), the average distance in pixels. For each θ, F builds K, turns the angles into R = R_XR_YR_Z, forms S′ = [R^T t′; 0^T 1], and projects every point with λ_iŵ_i = KPS′M′_i.

The angles enter through sines and cosines, and each ŵ_i is divided by λ_i, which depends on θ too. So, unlike least squares, the minimum of J has no formula. An iterative search finds it instead: start from a guess θ₀, and keep changing θ while J goes down. Some searches, such as the Nelder–Mead simplex, only compare values of J; others follow its slope. Either way, a search finds the lowest point near where it starts, which may not be the lowest overall, so the guess matters. Section 3 gives a good one: the parameters of the linear A.

Why search, when the linear A already fits? For three reasons.

• J measures what matters, distances in the image, with every point counting the same.
• The model can grow. Lens distortion, in the next section, or a manipulator that turns the object between images, adds parameters that no matrix A can hold.
• Knowledge about the camera can be built in. For a camera known to have no skew, γ can be fixed at 0, which the twelve free entries of A cannot express.

An average of distances, or of squares

J averages distances, not squared distances, and the two lead to different cameras. For the table, the smallest J is 2.69 pixels, reached with five of the eight corners fitted exactly and the other three up to 10.6 pixels off. Minimizing the squared distances instead spreads the error over all eight corners, none more than 5 pixels off, with an average distance of 3.40 pixels. An average of distances behaves like a median: it puts up with a few large errors when the rest fit exactly, so one bad click pulls it less. A square makes one large error cost more than several small ones, and most calibration software minimizes squares and reports their root mean square.

The average of distances also has sharp creases, where one distance reaches 0 and its slope jumps. They can stall a search: for the table, the Nelder–Mead simplex stopped anywhere between 2.9 and 3.3 pixels, depending on where it started. Solving a sequence of weighted least-squares problems, each one weighting a point's squared error by one over its current distance, reached 2.69 from each of several starts.

A mirror image that fits

One wrong camera fits the pixels exactly as well as the right one. Turn the camera half a turn about its own Y axis, so that it faces away from the object, and change the signs of β and γ. In the camera's coordinates, the turn changes the signs of X and Z for every point. Dividing −X by −Z changes nothing, but Y is now divided by −Z, and the new signs of β and γ, the two parameters that multiply Y, make up for it. The new camera's matrix is exactly −A, and −A gives the same image as A. But every λ_i is now negative: all the points are behind the camera. A search that starts from such a guess, with a negative β, stays with it. So after calibrating, check that α and β are positive and that every λ_i is positive.

5Lens distortion

The model so far assumes that rays go straight through the optical center. A real lens bends them a little, most near the edges of the image, where its surface is most tilted from the optical axis. Straight lines in the scene then look curved in the image.

Distortion is modeled as one more step between the projection and the conversion to pixels. P gives the ideal image point (x, y), a distortion function moves it to the real image point (x′, y′), the one the sensor records, and K converts that to pixels:

The full model has four steps: S′ takes the point to the camera's coordinates, P projects it onto the ideal image, the distortion moves it to the real image, and K gives its pixel. Only the real image is ever seen.

Radial and tangential

Distortion is usually split into two parts, around a center where it vanishes, most often the principal point. In polar coordinates (r, θ) around that center, a radial part δ_r(r) moves a point along its ray from the center and depends only on its distance r. A tangential part δ_θ(θ) turns it about the center and depends only on its angle θ. Back in x and y, with the center at the origin:

The first terms move the point by δ_r along the direction (cos θ, sin θ); the second terms move it by an arc r δ_θ along the perpendicular direction (−sin θ, cos θ). Both parts are usually polynomials of order higher than one. A δ_r proportional to r would only scale the whole image, which α and β already do, and a constant δ_θ would only turn it, which the angles already do.

A radial distortion that pushes points outward more and more toward the edges bends the straight lines of a square grid inward, like a pincushion. One that pulls them inward bows the lines outward, like a barrel, which is common in wide-angle lenses.

An example

A lens images nine points on a grid, at x and y in {−1, 0, 1}. The center and the four middles of the sides stay where they are, but the four corners move from (±1, ±1) to (±1.2, ±1.2). Every point stays on its own ray from the center, so there is no tangential distortion: δ_θ = 0. The distance from the center changes only at the corners, from √2 to 1.2√2, so δ_r(1) = 0 and δ_r(√2) = 0.2√2. A polynomial that vanishes at r = 0 and r = 1 is δ_r(r) = k(r³ − r), and at r = √2 it gives k√2, so k = 0.2. Since cos θ = x/r and sin θ = y/r,

Check a corner: x′ = 1 + 0.2 · 1 · (1 + 1 − 1) = 1.2. The middle of a side, (1, 0), has x² + y² − 1 = 0 and stays put. Other polynomials fit the nine points too. This one has only odd powers of r, which keeps δ_x and δ_y polynomials in x and y.

Why only the nonlinear method can include it

A distortion bends straight lines, and no matrix can do that: A sends every line to a line. So a linear calibration cannot model distortion. A nonlinear one adds its coefficients to θ and estimates them with the rest.

6Outliers and RANSAC

Least squares trusts every point. If one corner is clicked in the wrong place, say corner 6, 150 pixels to the left, the fit bends toward it and every corner suffers. The errors grow from 2–5 pixels to 9–67 pixels, J grows from 3.69 to 35.7 pixels, and the camera found moves by 22 cm. The wrong corner does not even stand out: its error, 56 pixels, is smaller than that of corner 5, which was clicked correctly. A wrong point in the data is called an outlier.

RANSAC, for random sample consensus, turns the question around. Instead of fitting all the points at once, it repeats four steps:

1. Pick a random sample with as few points as the model needs: six corners for a camera.
2. Fit A to the sample.
3. Count how many of all the points this A reprojects within a tolerance of their pixels. These are its inliers.
4. Keep the sample with the most inliers seen so far.

After enough tries, some sample is likely to hold no outlier, and its A agrees with every good point, which gives it the most inliers. A final least-squares fit on those inliers gives the answer. The sample is kept as small as possible because a small sample is more likely to be free of outliers, so fewer tries are needed.

The same method estimated homographies in the note on 2D transformations, with samples of four point pairs.

RANSAC with only eight corners

Eight corners are few for a model with eleven parameters. Six corners give twelve equations for eleven unknowns, so A fits almost any six corners closely, and a random sample says little about the rest. With so few corners, though, there is no need to sample: the figure in section 3 fits A to every set of eight, seven or six corners, 37 sets in all, and keeps the largest set whose corners all land within the tolerance of the fit to them. That set is the consensus RANSAC looks for with random samples. Switch the figure to RANSAC and drag a corner far away. For corner 6 clicked 150 pixels off, the largest such set is the other seven corners, their fit reprojects each within 3.4 pixels, and corner 6 is left out as an outlier, 162 pixels from where the fit puts it.

The tolerance is 20 pixels, because a camera fitted to seven of the corners predicts the eighth about 7 to 24 pixels off even when every click is right, and a tighter tolerance would call good corners outliers. And a corner dragged only a little can sometimes be explained as well by leaving out a neighboring corner and shifting the principal point. So when two sets of the same size qualify, the figure prefers the camera with its principal point nearest the center of the image and the least skew. Drag the corner farther, and it is the one left out.

Clicking points by hand is slow. Programs find corresponding points on their own, by detecting the corners of a chessboard or by matching keypoints between images, and some of those matches are wrong. RANSAC is what makes that practical. Between two views of one scene, the epipolar constraint gives a quick test for a proposed match.

7Calibration in practice

Calibration objects are usually flat: a printed chessboard, whose inner corners a program finds by itself to a fraction of a pixel. By section 2, one view of a plane is not enough, since it only fixes a homography. But every view of the board, held in a different pose, gives another homography, and all of them share the same K. A handful of views in different poses fix the five intrinsic parameters, and then each view's R′ and t′.

OpenCV does all of this in one call, cv.calibrateCamera. It takes the board's corners in the board's own coordinates, with Z′ = 0, and their pixels in every view. It returns K, the coefficients of its lens distortion model, and, for each view, R′ and t′. OpenCV takes the camera's system as the new one, the other way round described in section 3, and stores R′ as a rotation vector: an axis whose length is the angle turned about it. The matrix A of each view is then [K 0][R′ t′; 0^T 1].

X-ray imaging systems are calibrated the same way, from images of a calibration board. The X-ray source plays the part of the optical center, as the note on 3D transformations showed.

8Summary

The key ideas next to their NumPy equivalents. In code, the 3D points are the columns of M, their pixels the columns of w, and A is P.