A Flexible New Technique for Camera Calibration

Zhengyou Zhang of Microsoft Research showed that a camera can be calibrated from a few photographs of a flat printed pattern, with neither the camera motion nor the plane motion known. He presented it as “Flexible Camera Calibration by Viewing a Plane from Unknown Orientations” at the Seventh International Conference on Computer Vision (ICCV) in Corfu in September 1999, and as “A Flexible New Technique for Camera Calibration” in IEEE Transactions on Pattern Analysis and Machine Intelligence (TPAMI) in 2000. The technical report MSR-TR-98-71, first issued in December 1998, is the fullest account; the five-page journal paper keeps the derivation but refers to the report for distortion, degeneracies, and most experiments, and reports a different real experiment. This article treats the three as one publication and notes where they differ. The ICCV paper received the 2013 Helmholtz Prize. The procedure as used today is described in camera calibration.

Problem

Zhang divided existing techniques into two groups. Photogrammetric calibration observed an object of precisely known 3D geometry, usually two or three orthogonal planes or a plane under a precisely known translation as in Tsai’s 1987 method; it was efficient but needed expensive apparatus and an elaborate setup. Self-calibration used no object, relying on the rigidity of a static scene as Maybank and Faugeras had shown; Zhang judged it flexible but not yet mature, with too many parameters to estimate reliably. Aiming at a desktop vision system for non-experts with cheap cameras, he placed his method between the two: it uses 2D metric information, the known layout of points on a plane.

Triggs (ECCV 1998) had autocalibrated from five or more views of an unknown planar scene, a method Zhang called more flexible than his but hard to initialize. Liebowitz and Zisserman (CVPR 1998) rectified planes from metric information such as known angles and, as Zhang summarized, noted that three such planes would calibrate the intrinsics, without showing experiments. Zhang credits Zisserman with pointing out that this work uses the same constraint in a different form.

Contribution

• A known plane is enough. A pattern printed on a laser printer and attached to a flat surface, such as a hard book cover, determines the intrinsics from at least two orientations, with the camera or the pattern moved by hand.
• Two linear constraints per view on the image of the absolute conic, from each view’s homography: three views fix the five intrinsics linearly, or two if the skew is zero.
• A complete procedure: closed-form initialization, radial distortion, maximum-likelihood refinement by nonlinear least squares, degenerate configurations, and tests on simulated and real data.

Method

Constraints from one plane. Zhang writes the intrinsic matrix as $A$; this article uses $K$, as the neighboring articles do. With the pattern on the plane $Z = 0$, model points map to the image by $H = [\mathbf{h}_1\;\mathbf{h}_2\;\mathbf{h}_3] = \lambda K[\mathbf{r}_1\;\mathbf{r}_2\;\mathbf{t}]$. A homography has eight degrees of freedom and the pose six, so orthonormality of $\mathbf{r}_1$ and $\mathbf{r}_2$ leaves two constraints per view:

$$\mathbf{h}_1^\top B\,\mathbf{h}_2 = 0, \qquad \mathbf{h}_1^\top B\,\mathbf{h}_1 = \mathbf{h}_2^\top B\,\mathbf{h}_2, \qquad B = K^{-\top}K^{-1}.$$

The report and the journal paper derive them geometrically, from the condition that the images $\mathbf{h}_1 \pm i\,\mathbf{h}_2$ of the plane’s two circular points lie on the image of the absolute conic $B$; the ICCV paper only notes that $B$ is that image.

Closed-form solution. Writing the symmetric $B$ as the 6-vector $\mathbf{b} = [B_{11}, B_{12}, B_{22}, B_{13}, B_{23}, B_{33}]^\top$ makes $\mathbf{h}_i^\top B\,\mathbf{h}_j = \mathbf{v}_{ij}^\top\mathbf{b}$, so each view contributes the rows $\mathbf{v}_{12}^\top$ and $(\mathbf{v}_{11} - \mathbf{v}_{22})^\top$, and $n$ views give

$$V\mathbf{b} = \mathbf{0}, \qquad V \in \mathbb{R}^{2n \times 6},$$

solved as the right singular vector of $V$ for the smallest singular value. Three views determine $\mathbf{b}$ up to scale in general; with two, the zero-skew row $[0, 1, 0, 0, 0, 0]\,\mathbf{b} = 0$ is added. The report and the journal paper note that one view gives two parameters, such as the focal lengths when the principal point is known and the skew zero. Closed-form formulas then give the intrinsics from $B$; in the one for $u_0$, the ICCV and TPAMI texts print $\gamma v_0/\alpha$ where the report, revised in 2002 for a typo in its Appendix B, has the correct $\gamma v_0/\beta$. Each rotation computed from $K^{-1}H$ is replaced by the nearest rotation in the Frobenius norm, which the report derives as $UW^\top$ from the singular value decomposition $USW^\top$.

Refinement and distortion. Because the closed form minimizes an algebraic distance with no physical meaning, Zhang refines all parameters by minimizing the reprojection error, the maximum-likelihood estimate under independent, identically distributed image noise, using Rodrigues vectors for rotations and Minpack’s Levenberg–Marquardt. The ICCV paper and the report model two radial lens distortion terms, citing reports that richer models gain little and can be numerically unstable; the journal paper refers to the report:

$$\breve{x} = x + x\left[k_1 r^2 + k_2 r^4\right], \qquad \breve{y} = y + y\left[k_1 r^2 + k_2 r^4\right], \qquad r^2 = x^2 + y^2,$$

where $(x, y)$ are ideal and $(\breve{x}, \breve{y})$ distorted normalized coordinates, with the distortion centered at the principal point. Given the other parameters, $k_1$ and $k_2$ are linear, but alternating that step with the refinement converged slowly, so the recommended procedure estimates them linearly once and then refines all parameters jointly.

Degenerate configurations. A pure translation of the plane adds no constraint, and neither does a plane parallel to an earlier position, which the ICCV paper states and the report proves: parallel planes share the same circular points. The report adds that a known pure translation still allows calibration; the journal paper only mentions the parallel case.

Example

The code below uses the camera, plane orientations, and lateral offsets of the paper’s simulations, but not their depths: the listed depths of 500 to 525, stated in centimeters, would shrink the $18 \times 25$ cm pattern to about 45 pixels, so the code uses one tenth of them. It estimates each homography by the direct linear transform, solves $V\mathbf{b} = \mathbf{0}$, and applies the paper’s formulas for the intrinsics, then tries two views with the zero-skew row and a pair of parallel planes.

import numpy as np
from scipy.spatial.transform import Rotation

# Camera and plane orientations from Zhang's simulations; depths 1/10 of the listed ones (cm).
A = np.array([[1250.0, 1.09083, 255], [0, 900, 255], [0, 0, 1]])
X, Y = np.meshgrid(np.linspace(0, 18, 10), np.linspace(0, 25, 14))
M = np.c_[X.ravel(), Y.ravel(), np.ones(X.size)]            # model points [X, Y, 1]
deg = np.pi / 180
poses = [([20 * deg, 0, 0], [-9, -12.5, 50]),
         ([0, 20 * deg, 0], [-9, -12.5, 51]),
         (np.array([-30, -30, -15]) * deg / np.sqrt(5), [-10.5, -12.5, 52.5])]


def image_of_plane(rvec, t, sigma, rng):
    R = Rotation.from_rotvec(rvec).as_matrix()
    H = A @ np.c_[R[:, 0], R[:, 1], t]                         # H = A [r1 r2 t]
    m = M @ H.T
    return m[:, :2] / m[:, 2:] + sigma * rng.standard_normal((len(M), 2))


def homography(m):
    """DLT with Hartley normalization of the image points (L x = 0)."""
    c, s = m.mean(0), np.sqrt(2) / np.linalg.norm(m - m.mean(0), axis=1).mean()
    T = np.array([[s, 0, -s * c[0]], [0, s, -s * c[1]], [0, 0, 1]])
    u, v = ((m - c) * s).T
    Z = np.zeros_like(M)
    L = np.vstack([np.c_[M, Z, -u[:, None] * M], np.c_[Z, M, -v[:, None] * M]])
    return np.linalg.inv(T) @ np.linalg.svd(L)[2][-1].reshape(3, 3)


def v_ij(H, i, j):
    hi, hj = H[:, i], H[:, j]
    return np.array([hi[0] * hj[0], hi[0] * hj[1] + hi[1] * hj[0], hi[1] * hj[1],
                     hi[2] * hj[0] + hi[0] * hj[2], hi[2] * hj[1] + hi[1] * hj[2],
                     hi[2] * hj[2]])


def V_matrix(Hs):
    return np.vstack([[v_ij(H, 0, 1), v_ij(H, 0, 0) - v_ij(H, 1, 1)] for H in Hs])


def intrinsics(V):
    """b = smallest right singular vector of V; then Zhang's formulas (u0 as corrected)."""
    B11, B12, B22, B13, B23, B33 = np.linalg.svd(V)[2][-1]
    v0 = (B12 * B13 - B11 * B23) / (B11 * B22 - B12**2)
    lam = B33 - (B13**2 + v0 * (B12 * B13 - B11 * B23)) / B11
    alpha = np.sqrt(lam / B11)
    beta = np.sqrt(lam * B11 / (B11 * B22 - B12**2))
    gamma = -B12 * alpha**2 * beta / lam
    u0 = gamma * v0 / beta - B13 * alpha**2 / lam
    return np.array([alpha, beta, gamma, u0, v0])


def show(label, x):
    print(f"{label:22s}" + "".join(f"{a:10.2f}" for a in x))


rng = np.random.default_rng(1)
print(" " * 22 + "     alpha      beta     gamma        u0        v0")
show("truth", [1250, 900, 1.09083, 255, 255])
for sigma in [0.0, 0.5]:
    est = [intrinsics(V_matrix([homography(image_of_plane(r, t, sigma, rng))
                                for r, t in poses])) for _ in range(100)]
    show(f"3 views, noise {sigma}", np.mean(est, 0))
    show("  std (100 runs)", np.std(est, 0))

# Two views: add the zero-skew row [0, 1, 0, 0, 0, 0] b = 0 (the camera above has gamma ~ 1).
V2 = np.vstack([V_matrix([homography(image_of_plane(r, t, 0.0, rng)) for r, t in poses[:2]]),
                [0, 1, 0, 0, 0, 0]])
show("2 views + zero skew", intrinsics(V2))

# Degenerate: the second plane is parallel to the first (rotated about its own normal only).
r1 = Rotation.from_rotvec(poses[0][0])
r2 = (r1 * Rotation.from_euler("z", 40, degrees=True)).as_rotvec()
Hp = [homography(image_of_plane(r, t, 0.0, rng)) for r, t in [poses[0], (r2, [-9, -12.5, 56])]]
s = np.linalg.svd(V_matrix(Hp) / np.linalg.norm(V_matrix(Hp), axis=1, keepdims=True))[1]
print("parallel planes: singular values of V", np.array2string(s, precision=2))

Output:

                           alpha      beta     gamma        u0        v0
truth                    1250.00    900.00      1.09    255.00    255.00
3 views, noise 0.0       1250.00    900.00      1.09    255.00    255.00
  std (100 runs)            0.00      0.00      0.00      0.00      0.00
3 views, noise 0.5       1250.09    900.12      1.20    254.96    254.91
  std (100 runs)            5.36      3.95      0.61      1.81      1.13
2 views + zero skew      1241.73    894.05     -0.00    258.00    252.84
parallel planes: singular values of V [1.45e+00 1.37e+00 1.57e-15 1.99e-17]

Three exact views recover all five parameters, including the small skew; with 0.5-pixel noise the focal lengths scatter by about 0.4% from run to run. The zero-skew row is exact only if the skew is zero: forcing the true $\gamma \approx 1.09$ to zero leaves $\alpha$ about 8 pixels off even with exact data. A parallel second plane leaves $V$ with rank 2, adding nothing to the first view.

Results

All three versions report that the refinement converged in three to five Levenberg–Marquardt iterations.

Simulations (ICCV paper and report). The simulated camera is the one in the example, its skew equivalent to an 89.95° angle between the image axes. With three planes and Gaussian noise from 0.1 to 1.5 pixels, averaged over 100 trials, the errors grew linearly with the noise. At 0.5 pixels, which Zhang considered larger than typical, the focal-length errors were below 0.3% and the principal-point errors around 1 pixel. Adding views from 2 to 16 reduced the errors, most sharply from two to three. With three views tilted by an angle from 5° to 75°, accuracy was best near 45°; at 5°, 40% of trials failed because the planes were nearly parallel.

Real data (ICCV paper and report). An off-the-shelf PULNiX CCD camera with a 6 mm lens and $640 \times 480$ images viewed a 17 cm pattern of $8 \times 8$ squares with 256 corners. Calibrating from the first 2, 3, 4, and all 5 of five images gave final focal lengths between about 830 and 833 pixels. With all five images, refinement reduced the RMS reprojection error from 0.881 to 0.335 pixels. The closed-form distortion estimate had the wrong sign, $k_1 = 0.136$ (pincushion) against a final $-0.228$ (barrel), which Zhang attributed to fitting the intrinsics without distortion; the joint refinement recovered the correct shape. Recalibrating from each of the five subsets of four images gave $\alpha$ with a sample deviation of 2.90 pixels, an aspect ratio of 0.99995 ± 0.00012, and a skew not significantly different from zero. Two images of a tea tin with 16 hand-picked matches gave a reconstruction whose two visible, nearly perpendicular sides met at 94.7°.

Real data (journal paper). Citing space limits, the journal paper reports a different experiment: a PULNiX camera with a 6 mm lens as before, a pattern of $9 \times 9$ squares laser-printed on A4 paper and attached to cardboard, and ten images, five close (set A) and five farther away (set B). Sets A, B, and A+B gave $\alpha$ between 834.0 and 836.2 pixels and an angle between the image axes within 0.1° of 90°. Only the distortion disagreed: from set B, where the pattern covered just the image center and distortion is weak, $k_2$ was 1.31 against about 0.37 from the other sets.

Imprecise patterns (report only). Random noise in the model points up to 15% of the square size (about 2 mm) barely affected the focal lengths, below 0.02%, but moved the principal point by about 20 pixels and made $k_1$ unreliable. Systematic bending mattered more, cylindrical more than spherical, and Zhang judged results still usable for non-planarity below about 3% of the pattern size.

Impact

The method let non-specialists calibrate with a printer and a few photographs; the ICCV paper and the report note it was already in routine use in Microsoft Research’s vision and graphics groups. The journal paper is widely cited. The PAMI Technical Committee lists the ICCV 1999 paper among the 2013 winners of the Helmholtz Prize, which recognizes ICCV papers with significant impact on computer vision research.

Limitations

• Algebraic initialization. The closed form ignores distortion and minimizes an algebraic residual, so it can be biased, as the initial $k_1$ on real data showed.
• A simple distortion model. Two radial terms centered at the principal point leave out tangential and higher-order distortion.
• An exact plane assumed. The model points are taken as exactly known and coplanar. The report’s own tests show that the principal point and the distortion coefficients are the parameters most sensitive to errors in the pattern.
• One camera. The paper calibrates a single camera with fixed intrinsics; multi-camera rigs, stereo pairs, and changing zoom are outside its scope.

What Came After

The journal paper notes that Sturm and Maybank (CVPR 1999) independently published similar plane-based calibration, with a simplified camera model whose image axes are orthogonal and an exhaustive study of the degenerate configurations for one and two planes.

According to its documentation, Jean-Yves Bouguet’s Camera Calibration Toolbox for Matlab took its initialization partly from the ICCV paper, with identical homography and final maximum-likelihood estimation but a different closed form for the intrinsics, and adopted Heikkilä and Silvén’s intrinsic model with two tangential distortion terms; one of its examples recalibrates Zhang’s five published images. OpenCV’s source describes its older calibration engine as a port of Bouguet’s toolbox based largely on the journal paper. Since OpenCV 4.14 that engine is an option (CALIB_DISABLE_SCHUR_COMPLEMENT), and the default refinement is a Levenberg–Marquardt solver that eliminates the per-view poses with a Schur complement, as in sparse bundle adjustment. Both engines start from a reduced form of Zhang’s initialization for planar targets: the principal point is fixed at the image center and the skew at zero, and the two focal lengths are solved by least squares from each view’s homography.