Calibrating a camera from a pile of corners

Two posts ago I derived Zhang’s method end to end: the homography DLT, Hartley normalisation, the absolute-conic solve, a Cholesky that hands back K. That post has no distortion model at all, on purpose: the report it was drawn from says plainly that radial and tangential distortion “can become very complex” and leaves them out. This post is the other half, drawn from a different project’s report and a different repo: a real board, a real Canon 500D, twenty real photographs, and the two questions Zhang’s method alone can’t answer: what does a lens actually do to a ray, and how much of that should you even try to fit.

The board

The project calibrates against a ChArUco board (the same idea post 38 covers), generated by CameraCalibration/BoardInfo.py:

dpmm = 40
A4_shape = 210, 280
desired_block_size_mm = 20
desired_aurco_size_mm = 14

blocksx = A4_shape[0]//desired_block_size_mm   # 10
blocksy = A4_shape[1]//desired_block_size_mm   # 14

aurcoDict = aruco.getPredefinedDictionary(aruco.DICT_4X4_250)
charucoBoard = aruco.CharucoBoard_create(blocksx, blocksy, desired_block_size_mm,
                                          desired_aurco_size_mm, aurcoDict)

Ten blocks by fourteen, 20 mm squares, printed on A4. The interior chessboard corners (the ones actually used for calibration, via charucoBoard.chessboardCorners) are one fewer in each direction: 9 × 13 = 117 corners per fully-visible frame. That is the board the widget below uses too.

What a camera matrix does to a straight line, and why it lies

Post 41 already covers K = [[fu,s,cu],[0,fv,cv],[0,0,1]] and the pinhole model in full, so I won’t re-derive it: go read that post first if you haven’t. The pinhole model it builds assumes every straight line in the world stays straight in the image, and no real lens manages that at the edges of the frame. A real lens bends rays more the further they are from the optical axis, and the sensor’s stack of filters over the photosites shifts the effective centre of projection by a fraction of a pixel in a way that depends on direction, not just distance. Two families of coefficients absorb both effects, applied to the point after the perspective divide and before K, normalised coordinates $(x, y) = (X/Z,\, Y/Z)$ , with $r^2 = x^2+y^2$ :

\begin{aligned} x_d &= x\,(1 + k_1 r^2 + k_2 r^4 + k_3 r^6) + 2 p_1 x y + p_2 (r^2 + 2x^2) \\ y_d &= y\,(1 + k_1 r^2 + k_2 r^4 + k_3 r^6) + p_1 (r^2 + 2y^2) + 2 p_2 x y \end{aligned}

then $(u, v) = K \cdot (x_d,\ y_d,\ 1)$ . This is OpenCV’s model: a polynomial in even powers of the radius for the lens’s radial bending ( $k_1, k_2, k_3$ ), plus a smaller term for the sensor not sitting quite square to the lens ( $p_1, p_2$ , “tangential” or “decentering” distortion). It is not derived from first principles anywhere in the report; it is the standard model every camera calibration tool uses, and the report’s own contribution is deciding how much of it to trust.

Deciding how much of it to trust

Report p.11, discussing the projector-camera calibration (docs/pngs/report-11.png, §4.2.3):

Only around twenty photographs went into this calibration (p.14, §5.3.2: “A set of 20 images are scanned containing various angles of the calibration pattern… It is also ensured that there exist large changes in the angle of the pattern between images”). Five distortion coefficients plus five intrinsics plus six numbers per view is a lot of freedom to hand twenty photographs, and the report’s answer is to hand back some of that freedom: fit the coefficient that matters most ( $k_1$ ) and the two tangential terms, and hold the higher-order radial terms at exactly zero rather than let them chase whatever twenty images’ worth of corner-detection noise happens to look like.

The actual flags, CameraCalibration/main.py lines 72–74:

#CalibrationFlags=cv2.CALIB_ZERO_TANGENT_DIST + cv2.CALIB_FIX_K1 + cv2.CALIB_FIX_K2 + cv2.CALIB_FIX_K3
rep_err_camera, mtx_camera, dist_camera, rvecs_camera, tvecs_camera = cv2.aruco.calibrateCameraCharuco(
    all_charco_corners_camera, all_charco_ids_camera, BoardInfo.charucoBoard, camera_resolution,
    None, None, flags=cv2.CALIB_FIX_K2+cv2.CALIB_FIX_K3+cv2.CALIB_FIX_K4+cv2.CALIB_FIX_K5+cv2.CALIB_FIX_K6)

and, thirty-five lines later with the flags argument dropped entirely, the second calibration (lines 108–109):

rep_err_camera, mtx_camera, dist_camera, rvecs_camera, tvecs_camera = cv2.aruco.calibrateCameraCharuco(
    all_charco_corners_camera, all_charco_ids_camera, BoardInfo.charucoBoard, camera_resolution, None, None)

Both runs happen, back to back, on the same corners, and are saved to two different files: calculated_cams_matrix_less_distortion.npz for the fixed run, calculated_cams_matrix.npz for the unrestricted one. So the report can compare them. It only ever quotes the first.

Caveat

The commented-out line above the fixed call is stale, and the five active flags are not quite what the prose says. The comment describes an earlier intention (CALIB_ZERO_TANGENT_DIST plus fixing $k_1$ – $k_3$ ) that the actual flags= argument does not implement at all: it fixes $K_2$ through $K_6$ , never $K_1$ , and never touches the tangential terms. That matches the prose (“K₂ to K₅… set to 0”) much more closely than the dead comment does, modulo an off-by-one at the top end ( $K_6$ , not $K_5$ ). But $K_4$ , $K_5$ and $K_6$ are OpenCV’s rational model coefficients (a numerator/denominator extension for wide-angle lenses) which only exist in the fit at all if CALIB_RATIONAL_MODEL is also passed. It never is, here. So three of the five flags in that call are inert: fixing a coefficient that was never going to be estimated changes nothing. The two flags that do anything are CALIB_FIX_K2 and CALIB_FIX_K3, and that is exactly the pair the widget below lets you toggle.

The real numbers

Report p.21–22 (docs/pngs/report-21.png, report-22.png), §6.3, Eqs 6.1–6.5: the fixed run, the one the report treats as final:

\mathbf{K}_{camera} = \begin{pmatrix} 6192.0 & 0.0 & 2347.3 \\ 0.0 & 6190.0 & 1730.7 \\ 0.0 & 0.0 & 1.0 \end{pmatrix} \qquad \begin{aligned} k_1 &= 5.520\times10^{-2} \\ p_1 &= 6.015\times10^{-3} \\ p_2 &= -2.028\times10^{-4} \end{aligned}

No $k_2$ or $k_3$ is reported, because both are pinned at zero. A Canon 500D at roughly 4752 × 3168, so $f \approx 6191$ px is a little over 1.3× the sensor width in pixels, a long lens, or a fairly tight crop; the principal point at $(2347, 1731)$ is close to centred horizontally and about 150 px above centre vertically. This is what the widget below uses as its default ground-truth camera, and it is the number the “load corners” bundled sample was generated from.

Whose pose is it anyway: the parameterisation problem

Zhang’s method as post 41 built it stops at K. It never recovers where the camera actually was for each view, even though the information is sitting right there in each homography, post 41 says so explicitly, as a known gap. Recovering it is the first new piece here.

Write $\mathbf{H} = \mathbf{K}[\mathbf{r}_1\ \mathbf{r}_2\ \mathbf{t}]$ , exactly as Zhang’s derivation of the homography does (post 41’s Eq. 3.10). Then $\mathbf{r}_1, \mathbf{r}_2, \mathbf{t}$ are one $\mathbf{K}^{-1}$ away from $\mathbf{H}$ ‘s columns, up to a shared scale fixed by $\|\mathbf{r}_1\| = \|\mathbf{r}_2\| = 1$ , and $\mathbf{r}_3 = \mathbf{r}_1 \times \mathbf{r}_2$ completes the rotation. In noiseless arithmetic that is the whole answer; with real corners $\mathbf{r}_1$ and $\mathbf{r}_2$ come out close to but not exactly orthonormal, so the actual implementation (bundle::nearest_rotation) snaps the approximate matrix to the nearest proper rotation by an SVD/Procrustes step: decompose $M = USV^\top$ and take $R = UV^\top$ , flipping the sign of $U$ ‘s last column if that makes $\det(R) = -1$ .

That gives one pose per view to start from. Refining it is where the real difficulty is. A rotation matrix has nine numbers for three degrees of freedom, and an unconstrained Levenberg–Marquardt step on all nine would immediately walk off the manifold of valid rotations: nothing stops it turning $\mathbf{R}$ into a matrix that isn’t orthogonal at all. The fix is the standard one: parameterise each pose by a three-number axis-angle vector $\mathbf{w}$ (rotate $|\mathbf{w}|$ radians about $\mathbf{w}/|\mathbf{w}|$ ) and rebuild the full rotation from it on every evaluation via Rodrigues’ formula,

R(\mathbf{w}) = \cos\theta\, I + \sin\theta\, [\hat{\mathbf{w}}]_\times + (1-\cos\theta)\, \hat{\mathbf{w}}\hat{\mathbf{w}}^\top, \qquad \theta = \|\mathbf{w}\|,\ \hat{\mathbf{w}} = \mathbf{w}/\theta

three free numbers in, always a valid rotation out, whatever LM does to them. Going the other way, a fitted rotation matrix back to $\mathbf{w}$ (needed once to seed the optimiser from nearest_rotation’s output), goes through a quaternion (Shepperd’s method) rather than the textbook $\theta = \arccos\!\big((\mathrm{tr}\,R - 1)/2\big)$ formula, because that formula loses essentially all its precision for rotations near 180°. None of this is in the report; OpenCV does it inside cv2.aruco.calibrateCameraCharuco and never shows its work.

The bundle: everything, at once

With every view’s rotation reduced to three numbers, the full parameter vector is small: five intrinsics ( $f_u, s, c_u, f_v, c_v$ ), up to five distortion coefficients, and six numbers (axis-angle plus translation) per view. Twenty views is $10 + 120 = 130$ numbers, comfortably inside the “~100 parameters” the plan for this post estimated. Levenberg–Marquardt minimises total reprojection error over all of them jointly, not one homography at a time the way zhang::refine does, but every view’s pose and the shared camera and the shared distortion, simultaneously, because a change to $k_1$ moves every view’s residuals at once and a solver that refined one view in isolation would have no way to feel that.

The damped-normal-equations loop is the same shape as zhang::refine: accept a step if it lowers the cost, otherwise raise $\lambda$ and retry, linalg::solve for the linear system each try, reused deliberately, per the crate’s own README: “post 61 should reuse linalg::solve and the LM structure in zhang::refine… expect a bundle.rs rather than an edit to zhang.rs.” What is genuinely new is the Jacobian, and here this post departs from zhang::refine on purpose.

Caveat

The Jacobian is numerical, not analytic, and that is a deliberate departure from zhang::refine. Post 41’s homography refinement has a Jacobian you can write on one line per entry: a quotient rule and a $-1/w^2$ term. The bundle’s residual is a rotation exponential, composed with a perspective divide, composed with a five-coefficient distortion polynomial, composed with $K$ ; differentiating that by hand through nine parameter directions per view is exactly the kind of derivation zhangs-method’s own report warns about: the getError/getJ mismatch that post found was precisely a hand-differentiation slip nobody caught for years. So bundle::levenberg_marquardt builds each Jacobian column by central differences instead: perturb one free parameter by $10^{-6}$ of its own scale, re-run every residual, twice, subtract, divide. At the sizes this problem runs at (a few hundred points, a few tens of views, well under two hundred free parameters) that costs a few hundred extra reprojections per iteration, which is nothing next to forming and solving the normal equations. The whole calibration in the widget below, sixteen views and up to sixty LM iterations, runs in well under a second. It is measurably slower than a hand-derived Jacobian would be and it cannot be silently wrong the way one can.

Detecting garbage before it becomes a fake K

The plan for this post is explicit that LM “must converge on user garbage and detect degenerate all-coplanar input rather than emitting a nonsense K.” The bundle doesn’t re-invent this check: it calls zhang::calibrate first, purely for its initial K and per-view homographies, and if that fails (fewer than three views, a singular homography, or, the interesting case, every view close enough to fronto-parallel that the absolute-conic system is rank-deficient) the bundle refuses outright rather than seeding LM with nonsense and hoping it recovers. cargo test -p calib-wasm pins this down directly: six views, all pure in-plane rotation with no tilt at all, and bundle::calibrate reports Degenerate with every entry of K still zero, never a fabricated matrix.

Testing the report’s own argument

The most direct way to check “fixing k2/k3 helps with few images” is to reproduce it as a held-out generalisation test rather than trust it by eye. cargo test -p calib-wasm’s with_few_views_fixing_distortion_beats_fitting_it_on_held_out_views fits two models: one with all five distortion coefficients free, one with $k_2$ and $k_3$ fixed at zero, both against the report’s own true camera ( $k_1, p_1, p_2$ non-zero, $k_2 = k_3 = 0$ ), on three training views of a small sixteen-corner board, then evaluates both on views neither fit ever saw, averaged over six independent noise draws:

With sixteen good views instead of three, the gap mostly closes: the bundled sample below (16 synthetic views, the report’s own camera, 0.3 px of pixel noise) recovers $k_2 = -0.0130$ and $k_3 = 0.0602$ when both are left free: non-zero, because LM will always find something nonzero to fit noise with, but the training RMS is 0.294 px either way. The extra freedom doesn’t hurt visibly until you check a view it wasn’t fit on, which is exactly why the report’s own §5.3.2 procedure, “large changes in the angle of the pattern between images”, matters as much as the flag does. More, well-spread views is the other way to buy back the freedom the report chose not to spend on twenty images.

My numbers against the report’s

Feeding the bundled sample (16 synthetic views, projected through the report’s own K and distortion, 0.3 px noise) through bundle_calibrate with $k_2, k_3$ fixed: the same restriction as the report’s own fixed run:

	$f_u$	$f_v$	$c_u$	$c_v$	$k_1$	$p_1$	$p_2$
Report (Eqs 6.1–6.4)	6192.0	6190.0	2347.3	1730.7	0.05520	0.006015	−0.0002028
This implementation	6191.1	6189.4	2344.9	1734.8	0.05475	0.006226	−0.0004383

Every intrinsic agrees to within a few pixels and $k_1$ , $p_1$ to three significant figures, unsurprising, since the sample was generated from the report’s numbers and this is mostly confirming the pipeline round-trips its own ground truth through 0.3 px of injected noise. The one real outlier is $p_2$ : off by a factor of two, though both sides agree it’s small and negative. That’s the least surprising place for a mismatch: $p_2$ is the smallest-magnitude coefficient in the whole model by two orders of magnitude, so a handful of pixels of detector noise pushes it around proportionally far more than it moves $f_u$ or $k_1$ . The report doesn’t state an uncertainty on any of its numbers, so there’s nothing to compare that gap against except the general shape: small coefficients are the ones twenty images buys you the least confidence in, which is the report’s own argument again, one level down.

The report’s real reprojection error is 1.2779 px on real photographs (§6.3, p.22); the synthetic sample above lands at 0.294 px because it only has the 0.3 px of pixel noise I put there and none of a real lens’s departure from a five-coefficient model, none of a real detector’s failure modes, and none of JPEG compression’s effect on sub-pixel corner localisation. 1.2779 px on a 4752-pixel-wide sensor is still good, well under a pixel of error per axis on average, and the report’s own §7.3 makes the same point about the Kinect’s 0.6409 px on a much smaller 640-pixel sensor: a small number in pixels means more when the sensor supplying those pixels is small.

Calibrate a camera

InteractiveCalibrate a camera

With JavaScript on, this becomes a live calibrator. Pick a ground-truth camera and distortion with sliders and click “Add a view” to sample a random tilted pose and project noisy corners through it. A coverage heat-map shows which part of the sensor the captured views have actually touched. Or load the bundled sample (16 synthetic views through the report’s own recovered camera) or your own checkerboard-corners/v1 export. Press Calibrate: a Levenberg–Marquardt bundle refinement runs in WASM and streams its convergence into a small chart. Per-view RMS bars are clickable: pick one to see its detected corners (green) against its reprojected corners (magenta), or drop it and re-solve. Checkboxes fix k2/k3/p1/p2 at zero, so you can reproduce “fewer coefficients, better result with limited data” yourself.

Things worth trying:

Add three or four views, all with similar tilt, then Calibrate. Watch $c_v$ especially (the least determined intrinsic, same as in post 41) swing around between runs. Add ten more views at varied angles and it settles.
Load the bundled sample, calibrate with nothing fixed, then tick k2 and k3 and calibrate again. The recovered $f_u, f_v, c_u, c_v, k_1, p_1, p_2$ barely move and the RMS barely moves either: sixteen well-spread views is enough that the extra freedom mostly goes unused rather than actively hurting. This is the report’s argument in the direction it doesn’t emphasise: fixing coefficients is insurance against too few images, not a free lunch you should always take.
Set the noise slider high (1–2 px) and capture only three or four views, tightly clustered in tilt. Now toggling k2/k3 visibly moves the recovered $k_1$ itself: the restricted model is forced to explain the same residual pattern with fewer knobs, and it picks the more physically real one.
Click through the per-view RMS bars after a calibration with one deliberately extreme view included (tilt near 90°, or very close to the camera). Its magenta reprojected corners visibly diverge from the green detected ones before you even read the number. Drop it and re-solve.
Drag the distortion sliders to zero one at a time and watch the grid-warp preview straighten out: $k_1$ barrels or pincushions the whole frame, $p_1$ / $p_2$ shear it asymmetrically, which is a much smaller effect at these coefficient magnitudes and easy to miss unless you’re looking for it directly.
Try to make it lie. Every view fronto-parallel, tiny noise: the calibration refuses outright rather than returning a number, the same “Did not Converge” behaviour post 41 demonstrated for Zhang’s method alone.

What this does not do

Five coefficients, not eight. This implements OpenCV’s default radial-tangential model, the one the report’s own calibration actually fits (see the flags note above), not the rational ( $k_4$ – $k_6$ ), thin-prism or tilted-sensor extensions OpenCV also offers behind other flags, none of which this project or this post uses.
No covariance. Like post 41’s K, there is no uncertainty estimate on any recovered number. " $p_2 = -0.00044$ " and " $p_2 = -0.0002$ " look equally confident above and are not, see the discussion of that gap.
The Jacobian is numerical. Said above, worth repeating here: this is slower per iteration than a hand-derived analytic Jacobian and, at this problem’s size, the difference is not something you would notice without a profiler.

The Rust

Everything lives in one new module, bundle.rs, added to the calib-wasm crate post 38 and post 41 built. lib.rs gained exactly the one line the crate’s own README called for, pub mod bundle;. zhang.rs, linalg.rs, corners.rs, image.rs, grid.rs and traj.rs are untouched. It reuses zhang::calibrate and zhang::homography directly (for the initial K, per-view homographies and the coplanar-input check) and linalg::{svd, solve, inv3, mul3, transpose3} from post 41’s module. No new linear-algebra primitives were needed, only new code to assemble them into a bundle: the distortion model, Rodrigues forward and inverse, nearest_rotation, pose_from_homography, and the LM loop itself.

cargo test -p calib-wasm: 10 new tests for this module (61 in the crate total): the distortion model round-trips within $10^{-6}$ over a range of normalised coordinates and coefficients; the Rodrigues conversion round-trips including rotations close to 180°, where the naive trace-based formula loses precision; a homography-derived pose is a genuine rotation ( $R^\top R = I$ , $\det = 1$ ); LM’s RMS-per-iteration history never increases; a synthetic camera with non-zero radial and tangential distortion is recovered from noiseless correspondences to within $10^{-2}$ px on the intrinsics and $10^{-3}$ on every distortion coefficient; fixing a coefficient leaves it at exactly zero, not merely small, while the still-free $k_1$ visibly moves away from its own zero start; six fronto-parallel views are reported as Degenerate with a K of all zeros rather than a fabricated answer; too few views is reported the same way; the held-out generalisation test above; and a C-ABI round trip including null pointers and zero-sized inputs.

Size: the crate carried linalg.rs and zhang.rs at 60.0 kB after post 41. bundle.rs alone adds 33.4 kB (measured before post 60’s concurrently-developed board.rs module landed in the same crate). The committed public/blog/wasm/calib-wasm.wasm (which now also carries post 60’s ChArUco board-detection code, built in the same pass) is 145.6 kB, still comfortably under the crate’s ~200 kB budget.