Calibration
Intrinsic · Extrinsic · Temporal

Chapter 1: Intrinsics — What the Checkerboard Is Actually Solving For

Forget OpenCV for a chapter. You have a camera and a flat board with a printed pattern, and you can wave the board around in front of the camera as much as you like. What can you solve for, and what is the minimum number of images?

Try it cold, on paper, before reading on. Four minutes.

Count first, derive second

The unknowns. The intrinsic matrix is upper triangular with a 1 in the corner:

Five numbers: two focal lengths in pixels (f_x, f_y — they differ if the sensor's pixels are not square), the principal point (c_x, c_y) where the optical axis pierces the sensor, and the skew γ, which is non-zero only if the sensor rows and columns are not perpendicular. On any sensor manufactured this century, γ is zero to well under a pixel, and most pipelines fix it there. So: 5 unknowns, or 4 with skew fixed.

What one image gives you. This is the whole derivation, and it fits on three lines.

Put the board's own coordinate frame on the board, with Z pointing out of it. Then every board point has Z = 0, so the projection collapses:

The third column of the rotation never gets used, because it is multiplied by zero. What survives is a 3×3 matrix mapping board plane to image plane — a homography, which you fit from the detected corners with a DLT and a least-squares refinement. Call its columns h₁, h₂, h₃. Because homographies are only defined up to scale, there is an unknown λ:

Now the physics. r₁ and r₂ are two columns of a rotation matrix. That is not a decoration; it is two hard facts:

1. They are orthogonal: r₁^T r₂ = 0.
2. They have equal length: r₁^T r₁ = r₂^T r₂ = 1.

Substitute. The unknown λ cancels out of the first because the right-hand side is zero, and cancels out of the second because both sides carry it:

Two equations per image. That is the answer to the counting question — but the next sentence matters just as much, because it is the one that separates counting from reciting.

Count in the right frame. The unknown you actually solve for is not K. It is the 6-vector b holding the distinct entries of B = K^-TK^-1, and b is homogeneous — every constraint is of the form (row)·b = 0, so if b is a solution then 4b and −17b are solutions too. Six entries, one scale you cannot see: 5 degrees of freedom. To pin b down to a single ray you need the null space of V to be exactly one-dimensional, which means

With skew free, the images are your only source of rows: 2n ≥ 5 gives n ≥ 3. That is the answer to "minimum number of images."

How do you actually fix the skew? Everyone says "assume γ = 0" and almost nobody says how the solver is told. It is one line: γ = 0 forces B₁₂ = 0 (you will see why in the expansion two sections down), and B₁₂ is the second entry of b. So you append a literal extra row to V:

That row costs no images. So the count becomes 2n + 1 ≥ 5, and n ≥ 2.

Notice the third row buys you nothing in the count. It buys you a great deal in the conditioning, which is a different quantity and is what the debugging section is about. There is a gap here worth dwelling on: two images is enough, but would you ship two? The answer is no, and the reason is that rank 5 is a statement about the noiseless V, while what you ship depends on how far from rank-deficient the noisy V is. In practice you shoot fifteen to twenty-five.

The change of variable that makes it linear

Both constraints contain the same object, K^-TK^-1. Call it B. It is called the image of the absolute conic, a name from projective geometry that you should recognise but do not need for the derivation. What matters mechanically is: B is symmetric, so it has 6 distinct entries, and the constraints are linear in those 6 entries even though they are horribly nonlinear in K.

Stack them into a vector:

Now expand h_i^T B h_j and collect the coefficient of each entry of b. The only subtlety — and it is the one everyone gets wrong the first time — is that B is symmetric, so each off-diagonal entry appears twice:

with h_ki meaning the k-th element of column i. Then the two constraints are simply v₁₂^Tb = 0 and (v₁₁ − v₂₂)^Tb = 0. Stack 2n such rows into a matrix V, and b is the null vector of V — the last row of V^T from an SVD. One SVD, and you are done.

Deriving the closed form — where those six formulas come from

Every textbook drops six extraction formulas on you at this point, and every reader memorises them and forgets them by Thursday. You do not have to. Write out B in terms of the intrinsics and the formulas invert themselves. Do this once and you will never need the reference card again.

Step 1 — invert K, with γ = 0 for now. An upper-triangular matrix with a 1 in the corner has an upper-triangular inverse with reciprocal diagonal:

Step 2 — read B off as inner products of columns. B = K^-TK^-1 is a Gram matrix: entry (i, j) is just the dot product of column i with column j of K^-1. Name those columns

and six dot products later you have the whole matrix, with no matrix multiplication at all:

Step 3 — now just read the table backwards. Every extraction formula is one line, and every one of them is obvious once B is written this way:

Because B₁₃ = −c_x·B₁₁, dividing kills the focal length and leaves the principal point. Because B₁₁ is literally 1/f_x², one square root gives the focal length. That is the entire "closed form."

Step 4 — and the sixth entry is a redundancy check. Substitute c_x and c_y back into B₃₃:

Six entries, five parameters, one identity left over. So B has exactly the 5 degrees of freedom we counted, and that leftover identity is not waste — it is the ruler that measures the missing homogeneous scale. That is the whole reason Zhang's published formulas have a λ in them.

The λ, finally explained

The SVD does not hand you b. It hands you some multiple of b, because a null vector is only defined up to scale — numpy normalises it to unit length, which is an arbitrary choice with no physical meaning. Write what you actually receive as

Now run the identity from Step 4 on b̃ instead of b. Every term scales by μ (the ratios B₁₃²/B₁₁ carry one factor of μ, not two, because the numerator is quadratic and the denominator is linear):

There it is. λ is not a parameter and not a fudge factor — it is the measured value of the unknown SVD scale, recovered from the one algebraic identity the six entries had left over. Once you have it, the formulas that need an absolute magnitude get corrected and the ones that are ratios do not:

Watch it work on numbers you can check in your head. The true b for f = 800, c = (320, 240) is [1.5625×10^-6, 0, 1.5625×10^-6, −5.0×10^-4, −3.75×10^-4, 1.25] (we derive those six numbers in Worked Example 1 below). Suppose the SVD hands you four times that:

• c_x = −(−2.0×10^-3) / 6.25×10^-6 = 320.0  ✓  — the 4 cancelled, no λ needed
• c_y = 1.5×10^-3 / 6.25×10^-6 = 240.0  ✓
• B̃₁₃²/B̃₁₁ = 4.0×10^-6 / 6.25×10^-6 = 0.640;   B̃₂₃²/B̃₂₂ = 2.25×10^-6 / 6.25×10^-6 = 0.360
• λ = 5.0 − 0.640 − 0.360 = 4.000  ✓  — it recovered the 4
• f_x = √(4.000 / 6.25×10^-6) = √640000 = 800.0  ✓

Had you skipped λ you would have got √(1/6.25×10^-6) = 400 px — exactly half, because μ = 4 and √4 = 2. A calibration that reports half the focal length it should is almost always a missing homogeneous-scale correction, not a bad board. That is a real bug that has shipped, and it is worth recognising by its signature: the principal point is perfect and the focal length is off by a clean multiplicative constant.

And now the published version, with skew

Everything above assumed γ = 0. Zhang's paper does not, which is the only reason his formulas look forbidding. Redo Step 2 with the full upper-triangular inverse

and the same six dot products give

Two things fall out immediately, and both are worth saying at the whiteboard:

1. B₁₂ = 0 if and only if γ = 0. That is the justification for the [0, 1, 0, 0, 0, 0] row we appended to V above — it is not a heuristic, it is exactly the skew-free condition.
2. B₁₁B₂₂ − B₁₂² = 1/(f_x²f_y²) — the γ² terms cancel exactly. That combination is the 2×2 determinant of the upper-left block of B, and it is the skew-immune stand-in for B₂₂. Every place the γ = 0 formulas used B₂₂, the general ones use this determinant divided by B₁₁. That single substitution converts one set of formulas into the other.

Apply that substitution to the four lines from Step 3 and you have literally reproduced Zhang's Appendix B:

Check the collapse yourself in ten seconds: set B₁₂ = 0 in the right column. v₀ becomes −B₁₁B₂₃/(B₁₁B₂₂) = −B₂₃/B₂₂. β becomes √(λB₁₁/(B₁₁B₂₂)) = √(λ/B₂₂). γ becomes 0, so u₀ becomes −B₁₃α²/λ = −B₁₃(λ/B₁₁)/λ = −B₁₃/B₁₁. And λ collapses too, because v₀(−B₁₁B₂₃) / B₁₁ = −v₀B₂₃ = +B₂₃²/B₂₂. Every line lands back on the left column.

Worked example 1 — unpack K from B, every step

The derivation above tells you why the formulas have that shape. Now run them on real digits, forwards then backwards, which is how you sanity-check an implementation at 2 a.m. when you have a K you trust and a solver you do not.

Given f_x = f_y = 800, c = (320, 240), γ = 0. Invert K by inspection — for an upper-triangular K the inverse is upper triangular with reciprocal diagonal:

Form B = K^-TK^-1. Multiply out entry by entry:

• B₁₁ = 0.00125 × 0.00125 = 1.5625 × 10^-6
• B₁₂ = 0.00125 × 0 + 0 × 0.00125 = 0
• B₂₂ = 0.00125 × 0.00125 = 1.5625 × 10^-6
• B₁₃ = 0.00125 × (−0.4) = −5.0 × 10^-4
• B₂₃ = 0.00125 × (−0.3) = −3.75 × 10^-4
• B₃₃ = (−0.4)² + (−0.3)² + 1 = 0.16 + 0.09 + 1 = 1.25

Now run the closed form and watch 800 fall out. First the vertical principal point:

Numerator: 0 × (−5.0×10^-4) − 1.5625×10^-6 × (−3.75×10^-4) = 0 + 5.859375×10^-10.
Denominator: 1.5625×10^-6 × 1.5625×10^-6 − 0 = 2.44140625×10^-12.
Divide: 5.859375×10^-10 / 2.44140625×10^-12 = 240.000.  ✓

Then the overall scale λ (a different λ from the homography one — Zhang unfortunately reuses the letter):

B₁₃² = (5.0×10^-4)² = 2.5×10^-7.
v₀ × numerator = 240 × 5.859375×10^-10 = 1.40625×10^-7.
Sum = 3.90625×10^-7. Divide by B₁₁: 3.90625×10^-7 / 1.5625×10^-6 = 0.250.
So λ = 1.25 − 0.25 = 1.000.

Then the focal lengths:

And the horizontal principal point, with γ = −B₁₂f_x²f_y/λ = 0 because B₁₂ = 0:

Worked example 2 — when three images are not three images

The counting argument says n ≥ 3. But the counting argument is only a necessary condition, and the gap between necessary and sufficient is where real calibrations fail. So do not assert the degeneracy — compute it. Two views, nine numbers each, all of it doable on the back of an envelope.

Set up two fronto-parallel views. Same camera as always: f_x = f_y = 800, c = (320, 240). Board held square-on, so R = I in both, which means r₁ = (1, 0, 0) and r₂ = (0, 1, 0). Slide it sideways between shots: t_A = (0, 0, 0.60) and t_B = (0.05, 0.02, 0.60). Build H = K[r₁ r₂ t] column by column.

The first two columns are trivial — K times a standard basis vector just picks off a column of K:

The third differs between the views, because it is the only place the translation lives:

So the two homographies are genuinely different matrices:

Now form the constraint rows and watch what happens. Recall v_ij = [h_1ih_1j,  h_1ih_2j+h_2ih_1j,  h_2ih_2j,  h_3ih_1j+h_1ih_3j,  h_3ih_2j+h_2ih_3j,  h_3ih_3j]. For view A, column 1 is (800, 0, 0) and column 2 is (0, 800, 0), so plug in slot by slot:

• v₁₂ = [800×0,  800×800 + 0×0,  0×800,  0×0 + 800×0,  0×800 + 0×0,  0×0] = [0,  640000,  0,  0,  0,  0]
• v₁₁ = [800²,  2×800×0,  0²,  2×0×800,  0,  0] = [640000,  0,  0,  0,  0,  0]
• v₂₂ = [0²,  2×0×800,  800²,  0,  0,  0] = [0,  0,  640000,  0,  0,  0]
• v₁₁ − v₂₂ = [640000,  0,  −640000,  0,  0,  0]

And here is the whole point: every one of those numbers came from h₁ and h₂ alone. h₃ appears only in slots 4, 5 and 6, always multiplied by h₃₁ or h₃₂ — and both of those are zero for a fronto-parallel board, because r₁ and r₂ have no z-component. The translation cannot reach the constraint rows. So view B, whose only difference from A is its translation, produces the identical two rows:

Four rows, two of them literal duplicates. You can even get the singular values by hand, because the two distinct rows are orthogonal: stacking a row twice multiplies the squared singular value by 2, so σ = √2 × (row norm).

— from √2×640000√2 = 1280000 for the equal-norm row (norm 640000√2) and √2×640000 = 905097 for the orthogonality row. Rank 2. A four-dimensional null space. Not "poorly conditioned" — genuinely, exactly undetermined, and adding the other eighteen images changes nothing, because they all produce these same two rows. Twenty images, rank 2, confirmed numerically in the lab below.

The second degeneracy, proved rather than asserted

Now the subtle one, and it has a genuinely pretty proof. Tilt the board once, then spin it about its own normal. The board face still points the same way; you have only rotated the printed pattern within its own plane. Intuition says you have added information, because the images look completely different. Intuition is wrong, and here is why in three lines.

Spinning about the board normal by φ post-multiplies R by a z-rotation, which replaces the first two columns with rotated combinations of themselves:

K is linear, so h₁′ and h₂′ are the same combinations of h₁ and h₂. Write h_ij as shorthand for h_i^TBh_j and expand the new orthogonality constraint:

with c = cosφ, s = sinφ. Do the same for the equal-norm constraint:

Using cos2φ = c² − s² and sin2φ = 2cs, the new pair of rows is an exact linear map of the old pair:

And now compute the determinant of that 2×2:

Determinant exactly 1, for every φ. M is invertible always, so the two new rows span precisely the same 2-plane as the two old ones. Spin the board a hundred times at a hundred angles and the row space of V is identical to what one view gave you. This is not near-degeneracy that better numerics could rescue; it is an algebraic identity. Verified numerically to 5×10^-10 at φ = 17°, 40° and 90°.

Near-degenerate is the one that ships

Nobody hands you exactly-parallel boards. They hand you boards tilted by three degrees, because tilting is awkward and the corner detector likes fronto-parallel views. Then V is technically full rank and the solve technically succeeds — with a covariance the size of a house.

The structural quantity first. Take the noiseless V from eight views whose tilt direction is spread evenly around a circle at magnitude T, and look at σ₅/σ₁ — the smallest informative singular value, normalised so it is dimensionless. (σ₆ is the true null direction and is machine-zero; σ₅ is the one that decides whether the null space is really one-dimensional.)

That last column is flat to 7% over a 25× range of tilt. σ₅/σ₁ grows as T². The reason is visible in the hand computation above: tilt puts a z-component of size sin T ≈ T into r₁ and r₂, the constraints are quadratic forms in those columns, and so the information about the three entries that were previously unreachable enters at second order. Halving the tilt quarters the observability. That is the single most useful sentence in this chapter and it takes ten seconds to say.

Now the consequence, measured. Same eight views, 54 corners each, Gaussian corner noise with σ = 0.3 px, 120 independent Monte Carlo trials per row, board at 0.60 m, seed 20250728. You reproduce every cell of this table in Code Lab 2 below — it is not a quoted result, it is the printed output of code that ships with the chapter:

From 20° of tilt to 1° of tilt, the standard deviation of the recovered focal length degrades by a factor of 3172.8 / 5.4 = 590. Same number of images, same corner detector, same noise. All that changed is the geometry of the excitation.

Read the last column too, because it carries a second lesson. At 1.1° the estimator is not merely noisy — it is biased by a factor of 4.6, reporting a mean focal length of 3701 px for a camera whose true focal length is 800 px. Near-degenerate least squares does not scatter symmetrically about the truth; it runs away along the weak direction. Anyone who tells you "averaging more calibrations will fix it" has not looked at that column.

Does shooting more images rescue a badly tilted set? This is the exact question every team asks, so here is the measurement. Fixing the tilt and sweeping the image count instead:

Above about 6° the 1/√N law holds to within 3%, because there the error really is corner noise being averaged down. At 1.1° it fails outright — sixteen images came out worse than eight — because the error is not noise, it is rank deficiency, and near-duplicate rows do not average anything down. You cannot buy observability with sample size. The tilt-bench widget at the end of this chapter enforces exactly this: push tilt below 6° and the images slider stops helping.

Reconciling the table with the 2 px bar

Sharp readers will already have noticed something and it is worth saying before they ask: the best cell in that table is 3.9 px, and the debugging section below calls anything over 2 px unshippable. Is the bar invented? No — the table is deliberately a hard configuration, and the gap between the two is exactly the list of levers you have.

The table is 8 images, 0.3 px corner noise, closed form only. Production is 22 images, sub-pixel-refined corners at 0.2 px, and LM on top. Each of those three is a separate multiplier, and all three are measured:

1.66 px, under the bar, with a mean recovered f_x of 800.05 and a fitted per-coordinate RMS of 0.194 px against an input noise of 0.200 px. So the bar is reachable, and now you know exactly what it costs to reach it.

Two of those four multipliers deserve a second look, because they are the ones people get backwards:

• LM only buys 1.32×. Everyone assumes the nonlinear refinement is where the accuracy comes from. It is not — it is where the distortion comes from, and where the covariance comes from. The accuracy was decided when the operator chose how to wave the board.
• The corner detector is worth as much as the refinement. Going from an unrefined detector (σ ≈ 0.8 px) to an 11×11 sub-pixel window (σ ≈ 0.2 px) is a factor of 4 — strictly better than doubling the image count twice, and it is free. This is why the frontier section's learned-corner-detector work matters more than it first appears: it moves the multiplier that dominates.

Where distortion enters

Everything above assumed a perfect pinhole. Real lenses bend rays. Before writing the model down, derive its shape — because "why are there only even powers of r?" is the single most common follow-up in this whole chapter, and it has a clean two-sentence answer.

Why the radial series has only even powers. A lens ground on a lathe is rotationally symmetric about its optical axis. A symmetric optic cannot push a ray sideways — there is no preferred direction for it to push toward — so all it can do is move the image point along its own radius. That means the distorted point is the undistorted one times a scalar:

Now ask what g is allowed to be. It can only depend on the point through the rotation-invariant combination r² = x² + y², because that is the only function of (x, y) unchanged by rotating the image. And a smooth lens gives a smooth g, so expand it as a Taylor series in that combination:

That is the whole answer. The even powers are not a modelling convention, they are forced: a term in r¹ or r³ would be a power series in √(x²+y²), which has a corner at the origin and is not differentiable there — it would predict a kink in the lens at the exact centre of the image, which no ground glass has. The constant term is 1 because zero distortion must be the identity. Three coefficients is a cubic in r², which is where everyone stops for a normal lens.

And why the tangential terms look the way they do. If the lens element is not exactly parallel to the sensor — a fraction of a degree of decentering from the assembly line — the symmetry argument breaks, because now there is a preferred direction: the axis of the tilt. To first order that adds a displacement field which is not radial. The remarkable thing, and the fact that makes the p-terms memorable, is that the whole field is the gradient of one scalar:

Differentiate it and check: ∂φ/∂x = 2x(p₁y + p₂x) + r²p₂ = 2p₁xy + p₂(r² + 2x²). Exactly the published term, with no fitting and no hand-waving. φ is the lowest-order scalar that is radially weighted (the r²) but picks out a direction (the p₁y + p₂x), so (p₁, p₂) is literally a vector pointing along the decentering axis and √(p₁²+p₂²) is its magnitude. Verified numerically to 1×10^-9.

Put the two together and you have Brown–Conrady, applied in normalised coordinates, before K:

where x = (u − c_x)/f_x, y = (v − c_y)/f_y and r² = x² + y². Note where it is applied: in normalised coordinates, which is why r is dimensionless and why the same k₁ describes the lens no matter what resolution you read the sensor out at. Get that wrong — apply the model in pixels — and your k₁ is off by f².

Distortion breaks the linear derivation, because a distorted board is no longer related to the image by a homography. The standard resolution is bootstrap then refine: solve the closed form assuming zero distortion to get a starting K, then hand everything — K, the distortion coefficients, and all n board poses — to Levenberg–Marquardt minimising true reprojection error. The closed form is never the answer; it is the initialisation that stops LM landing in a bad basin.

Pin the sensor down before you quote a single distortion number

Here is where distortion quietly catches people, and where the YAML from Chapter 0 has a trap in it. Every distortion figure is meaningless until you know the image size, because r is normalised by the focal length and the corner's r depends on how many pixels away the corner is.

The rig used in this chapter, and in both of its Code Labs, is 640 × 480 with f_x = f_y = 800 and c = (320, 240) — which is exactly the centre of a 640 × 480 image, as it should be. Everything follows from those four numbers:

• Field of view. HFOV = 2·atan(320/800) = 2·atan(0.400) = 2 × 21.80° = 43.6°. Vertically, 2·atan(240/800) = 33.4°. Diagonally, 2·atan(400/800) = 53.1°. This is a moderately wide lens — a 28 mm-equivalent, not a fisheye. A genuine 90° horizontal field on a 640-wide sensor would need f_x = 320/tan(45°) = 320 px, not 800.
• The corner's normalised radius. The far corner is √(320² + 240²) = 400 px from the principal point — the 3-4-5 triangle, so this one is exact — and r = 400/800 = 0.500 exactly. That is why every distortion number below is a clean fraction.

Now the magnitudes, for judgment. The Chapter 0 coefficients are k₁ = −0.2814, k₂ = 0.0977, k₃ = −0.0163. Evaluate g(r) = 1 + k₁r² + k₂r⁴ + k₃r⁶ and multiply the pull by the radius to get the displacement in pixels:

Two things to carry out of that table. First, displacement grows as r³, not r — the pull is quadratic and you multiply by the radius. Quadruple the radius from 100 px to 400 px and the displacement goes 0.437 → 25.8 px, a factor of 59, close to 4³ = 64 (the k₂ and k₃ terms bend it back slightly). Second, the k₁-only estimate is 28.1 px at the corner and the full three-term answer is 25.8 px, so k₂ and k₃ together are worth 2.3 px — real, but a 9% correction on a term that is itself only visible in the outer ring. Remember that ratio; it is the whole argument of failure mode B below.

Ignore distortion entirely and your reprojection error at the corners is 129 times the noise floor. Ignore it at the centre and you would never notice. That asymmetry is why a single aggregate RMS is a bad calibration monitor — a board set that never reached the corners will report a beautiful RMS from a model that is badly wrong where it matters.

And the tangential terms, sized honestly. With p₁ = 0.0001 and p₂ = −0.0002 at the corner (x = 0.400, y = 0.300, r² = 0.250):

Magnitude √((9.00×10^-5)² + (5.00×10^-6)²) = 9.01×10^-5 in normalised units, which is 800 × 9.01×10^-5 = 0.072 px. Seventy-two thousandths of a pixel, at the worst point in the image, against a corner detector whose own noise is 0.2–0.3 px. The tangential terms on this lens are a third of the noise floor. That number — not a preference, not a rule of thumb — is the justification for CALIB_ZERO_TANGENT_DIST in the production call below. Compute it and you get to say "I measured it," which is a different sentence from "I usually turn it off."

The calibration service, with numbers

"Where does this live?" is the next question, and it wants rates, sizes and budgets, not boxes. Every number below is for the 640 × 480 mono8 rig pinned down in the previous section — fix the sensor size before you quote a single byte, because every figure in the flow is proportional to it.

The parameter ledger, because "how big is the optimisation?" is a follow-up and the honest answer is "depends what you fixed":

The whole reason a 2376×137 problem solves in 180 ms on one core is that Schur complement: 132 of those 137 parameters live in 22 independent 6×6 blocks that get eliminated analytically, leaving a 5×5 dense system to actually factorise. Marginalising structure out of a bundle is the same trick you will use in Chapter 7's SLAM backend on 10⁵ landmarks; here it is small enough to see whole.

The latency budget that actually matters is not the calibration — it runs once, offline, on a workstation. It is the undistortion that runs forever, per frame, on the robot. You never evaluate the Brown–Conrady polynomial per pixel at runtime; you precompute a remap look-up table once:

• LUT size: 640 × 480 × 2 maps × 4 B (float32) = 2457600 B = 2.46 MB, or 1.23 MB with the fixed-point 16-bit variant (OpenCV's CV_16SC2). Note this is 8× the size of the image it corrects — the LUT, not the frame, is the memory event.
• Per frame, bilinear remap touches 2.46 MB of LUT + 0.307 MB in + 0.307 MB out = 3.07 MB of traffic. On an embedded SoC with ~8 GB/s of usable DRAM bandwidth that is 3.07×10⁶ / 8×10⁹ = 0.38 ms, and it is memory-bound, not compute-bound — the arithmetic is four multiply-adds per pixel, roughly 1.2 Mflop for the whole frame, which any core does in microseconds.
• At 30 Hz and two cameras that is 3.07 × 30 × 2 = 184 MB/s of the memory budget spent on undistortion alone — 2.3% of the bus. Comfortable here; note that the same rig at 1280×960 would be 4× everything (9.83 MB LUT, 12.3 MB/frame, 1.5 ms, 738 MB/s, 9% of the bus), which is the point at which teams stop undistorting images. Then you undistort points instead (a few hundred features, a Newton iteration each, microseconds) and only rectify full frames when something downstream truly needs a rectified image.

Build the solve

The whole of Zhang that is worth being able to write is the constraint row. Everything else is an SVD call and a page of algebra you can look up. So that is what the lab makes you write — and then it runs the experiment from Chapter 0 on your own solver.

And the production form, which is what you would actually ship — say this sentence while you type it:

python
# everything above, plus distortion, plus LM refinement, in one call
rms, K, D, rvecs, tvecs, stdev_int, stdev_ext, per_view = cv2.calibrateCameraExtended(
    objpoints,          # list of (54,3) float32 - board model, IN METRES
    imgpoints,          # list of (54,2) float32 - detected corners
    (640, 480),   # the size these corners were detected at - NOT a header you copied
    None, None,
    flags=cv2.CALIB_FIX_K3 | cv2.CALIB_ZERO_TANGENT_DIST | cv2.CALIB_FIX_ASPECT_RATIO)

# stdev_int is the thing nobody reads and everybody needs:
#   [fx, fy, cx, cy, k1, k2, p1, p2, k3, ...] one-sigma, in pixels
print("fx = %.2f +/- %.2f px" % (K[0,0], stdev_int[0]))
assert stdev_int[0] < 2.0, "focal length is not observable from this image set"

What each argument is deciding. Five things, and they are the difference between "used OpenCV" and "owns calibration". Every one of them is a decision with a cost — anyone can name a flag; the cost is the part worth knowing.

1. "I pass objpoints in metres, and that number is the only metric input in the whole pipeline." Ties back to Chapter 0 unprompted.
2. "The image size argument is 640 by 480 because that is the resolution the corners were detected at. It is not a header I copied — and if it disagrees with where the principal point lands, one of the two is a lie." That sentence is the fix for the trap in the Chapter 0 YAML, and checking it at the moment you type the argument is much cheaper than discovering it in the field.
3. "I use the Extended variant because I want stdev_int. An RMS of 0.28 px with σ(f_x) of 40 px means the fit is tight and the parameter is not observable — those are different failures and the RMS cannot tell them apart."
4. "CALIB_FIX_ASPECT_RATIO forces f_x = f_y. I set it because on a modern sensor the pixels are square to better than 0.1%, so a free f_y is a parameter that buys me under a pixel of real signal and correlates strongly with the board's tilt about the horizontal axis — I would be spending an unknown to absorb a pose error. It costs me: free intrinsics drop from 4 to 3, the total from 141 to 140, and the Schur-reduced block from 9×9 to 8×8. If I ever calibrate an anamorphic lens or a sensor with binning applied on one axis only, this is the first flag I remove."
5. "I fix k₃ and the tangential terms unless I have evidence I need them. On this lens the tangential displacement at the worst point in the image is 0.072 px against a 0.2 px detector noise floor — I computed it — so those two parameters can only fit noise." Which is exactly the next section.

Note what those last two flags do together to the ledger: skew fixed by default, aspect ratio fixed by flag 4, and k₃ plus both tangential terms fixed by flag 5, leaves 3 intrinsics + 2 distortion = 5 free non-pose parameters, a total of 137, and a 5×5 Schur block. If you quote "141" while typing these flags, someone will notice. Quote the ledger row you are actually on.

Regenerate the tilt table yourself

The table in Worked Example 2 is the empirical spine of this chapter, and a number you cannot regenerate is a number you should not quote. So regenerate it. This lab is the exact script that produced every cell — same seed, same poses, same noise — and it prints the five rows to the digits shown above.

The two TODOs are the two ideas the table is about. Everything else is provided.

Three failure modes and the metric that reveals each

Failure mode A: insufficient tilt (the observability failure).

Failure mode B: over-parameterised distortion (the overfit).

Failure mode C: the intrinsics are right and the image size is wrong. This is the one hiding in the Chapter 0 YAML, and it is the most common calibration bug in a codebase that has ever changed resolution.

What is changing, briefly

Zhang's method is from Zhang, "A Flexible New Technique for Camera Calibration," TPAMI 2000, and it is still the default in every toolbox twenty-five years later, which tells you something. Three directions are actually moving:

• Learned corner detection. Hu, DeTone & Malisiewicz, "Deep ChArUco," CVPR 2019 replaces the classical corner refinement with a network, buying robustness under motion blur and extreme lighting where sub-pixel refinement simply fails. It does not change the estimator — it changes the noise σ that feeds it, which by the table above is worth as much as more images.
• Calibration-free geometry. Wang et al., "DUSt3R," CVPR 2024 and Leroy, Cabon & Revaud, "MASt3R," ECCV 2024 regress pointmaps directly from uncalibrated image pairs, and Wang et al., "VGGT: Visual Geometry Grounded Transformer," CVPR 2025 extends it to whole sets while also predicting camera intrinsics. The honest read for a robot: excellent for getting a plausible reconstruction from a phone video, not yet a replacement for a metric calibration with a covariance, because it produces no uncertainty you can gate on and no metres.
• Multi-camera batch toolboxes. Rameau et al., "MC-Calib," CVIU 2022 calibrates a whole rig of non-overlapping cameras against a set of ChArUco boards in one optimisation, which is what you actually need for a six-camera robot and what hand-rolling pairwise calibrations gets wrong.

Practice

Chapter 2: Extrinsics — The Rigid Transform Nobody Measures Directly

A scenario every perception team eventually lives through. A LiDAR is bolted under a camera. The intern drove the robot 60 metres down a straight aisle, logged both sensors, and ran the extrinsic calibration script. It converged. Is that good data?

The answer is no, the rotation is fine and the lever arm is completely unobservable, and you should be able to prove it in ninety seconds. Here is how.

The same motion, seen twice

Let X be the unknown rigid transform from sensor B's frame to sensor A's frame. The robot moves. Sensor A measures its own motion and calls it A. Sensor B measures the same physical motion and calls it B. Both are elements of SE(3).

Follow one physical point through the two routes. Route one: express it in B's old frame, move it with B, then convert to A. Route two: convert to A first, then move it with A. Because the sensors are rigidly attached, the two routes are the same thing:

That is the hand-eye equation, so named because it first appeared for a camera ("eye") bolted to a robot arm ("hand"). It is 4×4, but only 6 of its equations are independent, and the clean way to attack it is to split it.

Split it. Write X = (R_X, t_X), A = (R_A, t_A), B = (R_B, t_B). Multiply out both sides:

Matching the two parts gives two separate problems:

The rotation equation involves only R_X. So solve rotation first, then substitute it into a linear system for the translation. That decoupling is the single most useful structural fact about hand-eye calibration and it should be the first thing out of your mouth.

The rotation half is a Procrustes problem

R_AR_X = R_XR_B rearranges to R_A = R_X R_B R_X^T — a similarity transform. Take the matrix logarithm of both sides. The log of a rotation is its axis-angle vector: the direction is the axis, the length is the angle in radians. And conjugation by R_X just rotates that vector:

Write a_i = log(R_A,i) and b_i = log(R_B,i). The problem has collapsed into something you can see: you have two sets of arrows, and one set is the other set rigidly rotated. Find the rotation. That is the classic orthogonal Procrustes problem, and its solution is three lines:

The diag(1, 1, det) guard exists because U V^T can come back with determinant −1, which is a reflection. A reflection minimises the same least-squares cost but is not a rotation, and shipping one produces a left-handed coordinate frame and a very confusing week.

Worked example 1 — a Procrustes solve entirely by hand

People treat the SVD as a black box. It is not, at this size. Do this one on paper and the method stops feeling like magic.

The setup. The true answer — which the solver does not know — is a rotation of +90° about z:

Motion 1. The LiDAR reports a rotation of 0.40 rad about its own z axis, so b₁ = 0.40 e_z. The camera sees the same motion as a₁ = R_Xb₁. Since R_Xe_z = e_z (the third column), a₁ = 0.40 e_z.

Motion 2. The LiDAR reports 0.55 rad about its own y axis, so b₂ = 0.55 e_y. R_Xe_y is the second column of R_X, which is (−1, 0, 0) = −e_x. So a₂ = −0.55 e_x.

Build M. Two outer products:

e_ze_z^T puts a 1 in the bottom-right; e_xe_y^T puts a 1 in row 1, column 2. So:

The SVD, read off by inspection. Feed M each basis vector and see what comes out — each column of M is M applied to a basis vector:

• M e_y = (−0.3025, 0, 0) = 0.3025 · (−e_x). So σ₁ = 0.3025, v₁ = e_y, u₁ = −e_x.
• M e_z = (0, 0, 0.16) = 0.16 · e_z. So σ₂ = 0.16, v₂ = e_z, u₂ = e_z.
• M e_x = (0, 0, 0) = 0. So σ₃ = 0, v₃ = e_x, and u₃ is only pinned down up to sign: u₃ = ±e_y.

Assemble. U = [−e_x | e_z | s·e_y] and V = [e_y | e_z | e_x], with s = ±1 still undetermined. Then

Fix the sign with the determinant. Expanding along the first row: det = 0 − (−1)·(s·1 − 0·0) + 0 = s. A rotation needs det = +1, so s = +1, and

Exactly recovered, from two motions, by hand. And notice what the determinant guard just did: σ₃ was zero, so the data said nothing about the third singular direction — the guard supplied the missing bit of information from the requirement that the answer be a rotation rather than a reflection. That is not a numerical trick; it is the algorithm using a piece of physics the measurements do not carry.

Worked example 2 — the translation half, and why the corridor drive fails

With R_X in hand, the translation equation is linear:

Let us do a planar case where every number is visible. Robot on a floor, so all rotations are about z and we can work in 2D.

Truth (hidden from the solver): R_X = 90°, lever arm t_X = (0.30, 0.10) m — the camera sits 30 cm forward and 10 cm left of the LiDAR.

The motion. The LiDAR measures its own motion as a 90° turn with translation t_B = (0.50, 0.00) m. Generate what the camera must have measured:

R_Xt_B = [0 −1; 1 0](0.50, 0.00) = (0.00, 0.50).
R_A = R_XR_BR_X^T = 90° as well, since rotations about the same axis commute.
R_At_X = [0 −1; 1 0](0.30, 0.10) = (−0.10, 0.30).
From R_At_X + t_A = R_Xt_B + t_X:
t_A = (0.00, 0.50) + (0.30, 0.10) − (−0.10, 0.30) = (0.40, 0.30) m.

Now solve it back, pretending we only know R_A, t_A, R_X, t_B:

Right-hand side: R_Xt_B − t_A = (0.00, 0.50) − (0.40, 0.30) = (−0.40, 0.20).
Left matrix: R_A − I = [0 −1; 1 0] − [1 0; 0 1] = [−1 −1; 1 −1].
Its determinant: (−1)(−1) − (−1)(1) = 1 + 1 = 2.
Its inverse: (1/2)·[−1  1; −1 −1].
Multiply: (1/2)·[(−1)(−0.40) + (1)(0.20),   (−1)(−0.40) + (−1)(0.20)] = (1/2)·(0.60, 0.20) = (0.30, 0.10) m  ✓

And now the corridor. The intern drove straight. No turning means R_A = I, so

The left side is the zero matrix. Every value of t_X satisfies the equation equally well. Sixty metres of beautiful data and the lever arm is exactly as unknown as it was before you started. The solver "converged" because the least-squares system was rank-deficient and whatever damping or pseudo-inverse it used quietly returned the minimum-norm answer — which for a lever arm is zero, i.e. "both sensors are at the same place."

How much rotation? A conditioning number to quote

"Rotate" is advice. "Rotate 30 degrees" is engineering. Here is where the number comes from.

In 2D, the singular values of R(θ) − I are both equal to 2|sin(θ/2)|. Proof in one line: (R−I)^T(R−I) = 2I − (R + R^T) = 2(1 − cosθ) I, and the square root of 2(1−cosθ) is 2|sin(θ/2)|.

So an error δ in the measured motion translation becomes an error δ / (2 sin(θ/2)) in the recovered lever arm. With a realistic 1 cm error in the per-motion translation:

That table is the whole design of a calibration manoeuvre. It is why the standard IMU-camera dance is a figure-of-eight with aggressive wrist rotation rather than a smooth glide, and why "we collected two hours of driving data" is worth less than ninety seconds of deliberate excitation.

Where the extrinsic lives, and what an error costs in pixels

The extrinsic is not a file that sits on disk. It is a live edge in a transform tree that every consumer queries.

The pixel arithmetic that makes you sound like you have done this. Project a point at depth Z through a camera with focal length f. Then:

A rotation error does not have a Z in it. It shifts everything by the same number of pixels, near and far. A translation error is divided by depth, so it is enormous up close and vanishes at range. With f = 800 px:

They cross at Z = Δt / Δθ = 0.02 / 0.01745 = 1.15 m. Below about a metre you are looking at a translation error; above it, at a rotation error.

Build the rotation solve

The lab gives you the SO(3) log and exponential and the synthetic motions. You write the two lines that are the method — and then the check makes you prove the observability claims from this chapter on your own solver.

The production form, and the two sentences that go with it:

python
# R_gripper2base / t_gripper2base : the "A" motions (n of them)
# R_target2cam  / t_target2cam    : the "B" motions
R_cam2gripper, t_cam2gripper = cv2.calibrateHandEye(
    R_gripper2base, t_gripper2base,
    R_target2cam,   t_target2cam,
    method=cv2.CALIB_HAND_EYE_PARK)     # TSAI, PARK, HORAUD, ANDREFF, DANIILIDIS

# the check nobody writes and everybody needs: is the translation observable?
S = np.vstack([R - np.eye(3) for R in R_gripper2base])
smin = np.linalg.svd(S, compute_uv=False)[-1]
assert smin > 0.5, f"only {smin:.3f} of rotational excitation - lever arm is unobservable"

Say while typing: "PARK is the closed-form rotation-then-translation solve; DANIILIDIS solves both jointly with dual quaternions and is better when the rotation is poorly excited, because it lets translation information inform the rotation. And the assert is the part I care about — every hand-eye library will happily return a number from degenerate data."

Two failure modes and the metric that reveals each

Failure mode A: rotation-starved extrinsic (the corridor drive).

Failure mode B: a temporal error masquerading as an extrinsic error.

Targetless and observability-aware, briefly

• Targetless motion-based extrinsics. Taylor & Nieto, "Motion-based calibration of multimodal sensor extrinsics and timing offset estimation," ICRA 2015 / T-RO 2016 is the canonical statement that you can recover extrinsics and the time offset from ordinary driving, with an observability analysis attached. This is what a modern fleet does instead of a board dance.
• Appearance-based LiDAR-camera. Yuan et al., "Pixel-level Extrinsic Self Calibration of High Resolution LiDAR and Camera in Targetless Environments," RA-L 2021 aligns LiDAR depth discontinuities against image edges. It needs no motion at all, which sidesteps the observability problem entirely — and buys a different one, because edges are scene-dependent.
• Observability-aware data selection. Lv et al., "Observability-Aware Intrinsic and Extrinsic Calibration of LiDAR-IMU Systems," T-RO 2022 scores each segment of a trajectory by how much information it carries about each parameter and keeps only the informative parts. This is the right modern answer to "how much data do I need" — you do not need more data, you need data with a well-conditioned Fisher information matrix, which is exactly the argument in Fisher Information & the CRLB.
• Estimate it online instead. Geneva et al., "OpenVINS," ICRA 2020 ships camera-IMU extrinsics (and t_d) as optional filter states, so they are refined continuously in the field. The tradeoff is honest: online estimation tracks thermal and mechanical drift, but it can also absorb a genuine fault into the calibration and hide it, so you still want an offline reference to compare against.

Practice

Chapter 3: Temporal — Estimating the Time Offset as a State

Your camera and IMU are on different clocks. The hardware team says the offset is "about 15 milliseconds, maybe." What do you do?

The wrong answer is "ask them to measure it." The right answer is "I do not want a constant from the hardware team. I want t_d as a state in my estimator, with a Jacobian, a covariance, and an observability gate — because it is not constant, and because I can estimate it better than they can measure it."

The derivative that makes t_d estimable

Start from what an offset does. Your camera reports that a frame was captured at time t. It was really captured at t + t_d. Everything else in the estimator — the IMU integration, the pose prior, the extrinsic — is evaluated at t. So you predict where a feature should be at t, and you measure where it actually was at t + t_d.

Let u(t) be the image position of a tracked feature as a function of time. Then to first order:

where u̇ is the feature's image velocity in pixels per second — a quantity your tracker already computes for free, because it is the frame-to-frame displacement divided by the frame interval.

So write the residual with t_d in it explicitly:

and differentiate:

That is the entire trick, and it is one line. The Jacobian of the reprojection residual with respect to the time offset is minus the feature's image velocity. Add one column to your Jacobian, one entry to your state vector, and t_d is estimated jointly with everything else by the same Gauss–Newton you already run — see RIP 03 for why every estimator in this field is that same solve.

The degeneracy that is not zero motion

u̇ = 0 is the degeneracy everybody quotes, and it is the easy one — a stationary robot is visibly stationary. The degeneracy that actually costs a team a week is the opposite: plenty of motion, all of it the same.

Rotate at a perfectly constant ω. Then u̇ = fω is a large non-zero constant, the information ∑u̇²/σ² is enormous, every excitation gate you just wrote is green — and the residual an unmodelled offset creates is

Now ask what else shifts every feature by the same constant amount. A camera–IMU extrinsic whose yaw is wrong by Δθ rotates every predicted bearing by Δθ, which moves every predicted pixel by f·Δθ. Same shape, same sign structure, same everything. Put the two Jacobian columns side by side:

Both are constants, and their ratio is ω — the same ratio for every feature and every frame in the window. Two proportional columns are one column. Write the 2×2 information submatrix over (t_d, Δθ_yaw) for N features and read its determinant:

Rank 1. Determinant exactly zero, for every N and every ω. Stacking more features or more frames of the same rate multiplies a singular matrix by a bigger number; it does not make it invertible. This is a strictly harder failure than u̇ = 0, because all of your health metrics look excellent while it is happening.

The number, so it is not an abstraction. At ω = 0.30 rad/s a time offset of t_d = 18.4 ms produces fωt_d = 800 × 0.30 × 0.0184 = 4.416 px — Worked example 1's residual. A fixed yaw extrinsic error of

produces fΔθ = 800 × 0.00552 = 4.416 px. Identical, feature by feature, to every digit you can print. Ask the data which hypothesis it prefers and the data has no opinion. What comes back is whatever mixture the priors and the numerics happen to land on — usually some of each, which is the worst available outcome, because now both parameters are wrong and both report converged.

What breaks the tie is change of rate, not amount of rate. Take two frames, at rates ω₁ and ω₂. The 2×2 Jacobian over (t_d, Δθ_yaw) is

The determinant contains the difference of the rates and nothing else at all. A careful smooth sweep from 0.29 to 0.30 rad/s gives det = 640,000 × 0.01 = 6,400. A reversal between +0.30 and −0.30 rad/s gives 640,000 × 0.60 = 384,000 — sixty times better conditioned from the same amount of shaking. Rate magnitude buys information about t_d in isolation (that was Worked example 2). Rate variation is the only thing that buys the ability to tell t_d apart from everything else that also scales with rate.

Worked example 1 — how loud is an 18 ms offset?

Numbers, so the abstraction becomes a plot you can picture.

Step 1 — convert body rotation into image velocity. A camera rotating at ω rad/s sweeps the image at f·ω pixels per second, because a small rotation Δθ moves a feature by f·Δθ pixels. Take a gentle yaw of ω = 0.30 rad/s (about 17°/s, a slow turn) and f = 800 px:

Step 2 — the residual an unmodelled t_d creates. With t_d = 18.4 ms:

Step 3 — compare it to the noise floor. A good tracker has σ = 0.3 px. So the residual is

Fourteen sigma. This is not a subtle effect — during a gentle turn an unmodelled 18 ms offset is screaming. And yet it hides, because during the 80% of a log when the robot is going straight, u̇ is small and the residual is small, so the average residual looks acceptable. That is the trap: t_d faults are invisible in aggregate statistics and obvious in conditional ones.

Worked example 2 — how well can t_d possibly be known?

The natural follow-up: how precisely can you estimate it? That is a Cramér–Rao question, and it has a closed form here that is worth memorising.

The Fisher information a set of measurements carries about a scalar parameter is the sum over measurements of the squared derivative divided by the noise variance (see Fisher Information & the CRLB for where that comes from):

With N = 200 tracked features all moving at u̇ = 240 px/s and σ = 0.3 px:

u̇² = 240² = 57,600.
σ² = 0.3² = 0.09.
Per feature: 57,600 / 0.09 = 640,000 s^-2.
Times 200 features: I = 1.28 × 10⁸ s^-2.
The bound is σ_td ≥ 1/√I = 1/11,313.7 = 8.84 × 10^-5 s = 88 µs.

Cancel the algebra and the result is memorable:

One frame of a gentle turn pins the time offset to under a hundred microseconds. That is roughly two hundred times better than the hardware team's "about 15 milliseconds, maybe," from data you were collecting anyway. It is also why online estimation wins: the information is free and abundant.

And read the formula backwards for the failure. As u̇ → 0, σ_td → ∞. There is no amount of N that fixes a stationary robot, because N multiplies under a square root while u̇ multiplies directly. Ten times more features buys you a factor of 3.2; ten times more motion buys you a factor of 10.

Worked example 3 — one update, by hand

The Cramér–Rao bound says what is achievable. The next thing to demand of yourself is one actual update, worked in full. So here is a single scalar Kalman update on t_d with every intermediate written down. No matrix in it is bigger than 1×1, which is the whole reason to do it this way first.

The setup. The hardware team said "about 15 milliseconds, maybe." You take them at their word — as a prior, not as a fact: t̂_d⁻ = 15.0 ms with σ⁻ = 5 ms. The truth, which you do not know, is 18.4 ms. One feature is tracked at u̇ = 240 px/s (the gentle yaw of Worked example 1) and the tracker's noise is σ_px = 0.3 px.

Step 1 — the prior variance, in SI. Always SI. Milliseconds inside a covariance are how you earn a factor-of-10⁶ bug that survives three code reviews.

Step 2 — the Jacobian and the measurement noise. Straight from the derivative at the top of the chapter:

Keep the units in view, because they are the sanity check: H is pixels per second-of-offset, so H²P⁻ comes out in px² and is directly comparable with R. That comparison is the update; everything below is bookkeeping around it.

Step 3 — the innovation covariance. How large a residual should we expect, given what we do not yet know?

Read those two terms rather than summing them mechanically. The timing uncertainty is worth 1.44 px² of expected residual — 1.2 px of scatter. The tracker contributes 0.09 px² — 0.3 px. The thing we are trying to learn is sixteen times louder in variance than the noise we have to see through, which is why the update about to happen is decisive rather than incremental.

Step 4 — the gain.

Units check: seconds of clock offset per pixel of unexplained image motion. One pixel of residual is worth 3.92 ms of clock. That number is worth carrying in your head — it is the exchange rate between the two quantities in this problem, and it makes the "4.4 px = 18 ms" claim from Worked example 1 fall out again from the other direction.

Step 5 — the posterior variance.

KH = 0.94118 is not a new quantity: it is exactly 1.44/1.53, the share of the expected residual that the prior uncertainty owns. The measurement gets 94.1% of the vote, the prior keeps 5.9% — from one feature, in one frame, off a 5 ms prior.

Step 6 — the state. Take the residual noise-free so the arithmetic stays visible. With the truth at 18.4 ms and the prediction built at 15.0 ms:

We are carrying H as the derivative of the residual — which is what the chapter derived and what a Gauss–Newton code path actually holds in memory — so the step is t̂⁺ = t̂⁻ − K r. That is identical to the textbook x + Kν with ν = −r; only the sign convention differs, and mixing the two is the single most common way this update ships backwards:

The prior was 3.40 ms wrong. After one feature it is 0.20 ms wrong — and 0.20/3.40 = 0.0588 = 1 − KH, exactly. The gain is not a tuning knob; it is the fraction of your error that this measurement deletes. If you can say that sentence and then show the two numbers that prove it, you are done with this question.

Step 7 — all 200 features at once, in information form. Repeating step 4 two hundred times is arithmetic, not insight. The information form does the whole frame in one line, because independent measurements add information:

That is Worked example 2's Cramér–Rao bound, to the digit. It is not a coincidence, and reproducing it is the reason to do both calculations: for a linear-Gaussian problem the Kalman posterior attains the CRLB, so the bound you quote at the whiteboard is the standard deviation you will actually measure in the log. If your logged spread is much wider than the bound, you have a bug or an unmodelled correlation — not bad luck.

Notice the other number hiding in that line. The prior contributed 4 × 10⁴ out of 1.2804 × 10⁸ — three parts in ten thousand. After a single frame of gentle yaw, the hardware team's estimate is arithmetically irrelevant. That is the quantitative version of "I can estimate it better than they can measure it," and it is a much better answer than the assertion.

Step 8 — the same update below the observability gate. Set u̇ = 5 px/s, under the 20 px/s gate, and run the identical arithmetic with nothing else changed:

One feature moved σ from 5.000 ms to 4.983 ms. Seventeen microseconds of learning. All two hundred together give 1/P⁺ = 4 × 10⁴ + 55,556 = 95,556 s^-2, i.e. σ⁺ = 3.23 ms — after a full frame the answer is still dominated by the hardware team's guess, and 42% of the information in it came from a prior somebody typed into a YAML file. The gate is not conservatism. It is arithmetic: below it, the update is theatre.

The other route: batch correlation

Adding a state to a live estimator is the right long-term answer. But on day one, before the estimator exists, you want a number from a bag file. That is the correlation route, and it is the one you will most likely be asked to implement.

The idea: find two scalar signals that measure the same physical quantity through the two clocks. For camera and IMU the natural pair is angular rate — the gyroscope measures it directly, and the camera gives it as the frame-to-frame rotation divided by the frame interval. Then slide one against the other and find the alignment that maximises the normalised cross-correlation:

Why normalised and not a raw dot product. The two signals are in different units with different gains — one is rad/s from a gyro with a scale-factor error, the other is rad/s derived from pixels through a focal length. A raw dot product's peak location is affected by both signals' amplitude envelopes; the normalised version is invariant to any affine change (gain and bias) of either signal, so a gyro scale error cannot move your answer. That invariance is exactly what the lab's first assert tests.

Then refine below the sample grid. The correlation is only evaluated at integer sample lags, so the raw peak is quantised to one sample — 5 ms at 200 Hz, which is worse than the offset you are chasing. Fit a parabola through the peak and its two neighbours and take its vertex. For three points one unit apart with values y₀, y₁, y₂:

Derive it, do not memorise it. Put the middle sample at d = 0 and fit y = Ad² + Bd + C. Then y₀ = A − B + C, y₁ = C, y₂ = A + B + C. Subtracting: y₀ − y₂ = −2B, so B = (y₂ − y₀)/2. Adding: y₀ + y₂ − 2y₁ = 2A, so A = (y₀ − 2y₁ + y₂)/2. The vertex of a parabola is at d = −B/(2A), which substitutes to the formula above. Ninety seconds at a whiteboard.

Rolling shutter: a time offset that varies down the image

Almost every cheap CMOS sensor exposes rows sequentially. Row k is captured at

where t_line is the line delay, typically 10–30 µs. For a 960-row sensor at 20 µs per line, the bottom row is captured 960 × 20 µs = 19.2 ms after the top row. That is the same order as the whole camera–IMU offset, which is exactly why the two get conflated.

The distinction is sharp once you look at the right axis:

The number for the thirds test. With H = 960 and t_line = 20 µs, the mean row of the top third is 160 and of the bottom third is 800, a difference of 640 rows = 12.8 ms. So the fitted slope should differ by 12.8 ms between top and bottom thirds — larger than most t_d values, and unmistakable if you look.

Where the offset lives, with numbers

The 2×2 that decides whether the extra state helps you or hurts you

Byte counts are the easy half of "what does one more state cost." The half that decides whether the system works is the new row and column of cross-covariance, because a t_d that cannot be observed does not sit politely still — it parks its error in whatever state it is correlated with, and that state then lies to everyone downstream. So take the 2×2 corner of P over (t_d, b_gz) — the time offset and the z gyro bias — and do the update by hand, the way you did the scalar one.

Why those two share a Jacobian row. Over a sliding window of length Δt, an error δb in the z gyro bias integrates into an attitude error δθ = δb·Δt, which moves every predicted pixel by f·δb·Δt. So with Δt = 1.0 s (a typical keyframe window — for a 0.1 s window divide every coupling number below by ten, and the correlation falls with it):

The priors, and the trick of quoting them in pixels. Start from an uncorrelated P⁻ = diag(2.5 × 10^-5 s², (0.01 rad/s)²) — 5 ms of timing doubt and 10 mrad/s of bias doubt. Push each through its own column to see what it is worth as image motion, which is the only currency the measurement understands:

The bias is a 6.7× cheaper place to put a residual than the offset is. Hold on to that ratio; it decides everything that follows.

One measurement, both states free.

Feed it the same 0.816 px residual as Worked example 3 — a residual caused entirely and only by the 3.4 ms timing error:

Account for the 0.816 px. The bias absorbed 800 × 9.96 × 10^-4 = 0.797 px, 97.8% of it; t_d took 240 × 7.47 × 10^-5 = 0.018 px, 2.2%. The identical measurement that moved t_d by 3.20 ms when the bias was known moves it by 0.075 ms when the bias is free — a factor of 43 — and in exchange it has invented a full milliradian per second of gyro bias that the gyro does not have. Your timing fault has become a bias fault, and the bias state has no way of knowing it is lying.

Two hundred features do not fix it. Changing the rate does. Stack a whole frame, add the information, invert the 2×2. The only thing that varies between the rows below is whether ω keeps one sign across the window or reverses inside it:

(ρ is the correlation coefficient of the (t_d, b_gz) block. The bias standard deviations, in the same order, are 10.00, —, 1.484, 0.027, 0.041 and 0.027 mrad/s: note that the bias is estimated beautifully in every single row, including the ones where t_d learns nothing at all. A dashboard that watches bias convergence as a proxy for "calibration is healthy" will be green throughout this entire table.)

Row three is the one to stare at. Two hundred features, a healthy 240 px/s, every excitation gate green — and t_d improved from 5.000 ms to 4.945 ms. About one percent. The information matrix is rank 1: this is the degeneracy from the top of this chapter again, with the gyro bias now playing the part the extrinsic played there, and rank is not something N can buy. Row four is the same 200 features at the same speed, differing only in that ω changed sign halfway through the window: σ_td = 88 µs, a 56× improvement, and ρ collapses to zero.

Now read the correlation column carefully, because the obvious reading of it is wrong. |ρ| is largest in the well-excited-but-unvarying row (−0.9998) and smaller in the starved below-gate row (−0.762). Correlation on its own is therefore not the alarm. The alarm is σ_td refusing to shrink while |ρ| is large — that pair says "these measurements are informative about some combination of the two states, and it is not the one you asked for." Publish both, not either.

And here is the bill for that correlation. Take the posterior from row three — P = [[2.445 × 10^-5, −7.335 × 10^-6], [−7.335 × 10^-6, 2.201 × 10^-6]], ρ = −0.9998 — and now the robot stops. u̇ = 0 exactly, so the t_d column of the Jacobian is exactly zero: H = [0, −800]. A completely ordinary 0.5 px residual (1.7σ of tracker noise, nothing to see) arrives:

A state whose Jacobian column is exactly zero just moved by two milliseconds, on noise, because the cross-covariance dragged it along behind the bias. Equivalently δt_d = ρ(σ_td/σ_b)δb = −0.9998 × 3.333 × 0.5875 mrad/s. At 30 Hz, with the noise sign flipping frame to frame, that is precisely the random walk of tens of milliseconds that failure mode A describes below — now derived, with a number, instead of asserted.

So: freeze, or inflate R? Both are defensible in a design review, so run both through the same arithmetic instead of arguing about them.

Offset versus drift, because they are different faults. An offset is a constant; drift is an offset that grows because two crystals run at different rates. A cheap oscillator is specified to ±40 parts per million, which is 40 µs per second, 2.4 ms per minute, 144 ms per hour. So a system that estimates a single constant t_d at boot and never revisits it is fine for a two-minute mission and catastrophically wrong for an eight-hour shift.

Build the offline estimator

And the online form, which is the answer you give when they ask what you would actually ship. It is four lines added to a Gauss–Newton you already have:

python
# state x = [ ..., td ]  with td the LAST element
for i, f in enumerate(tracked_features):
    u_pred = project(T_wc, f.point_w, K)          # (2,)
    u_dot  = (f.uv - f.uv_prev) / dt_frame          # (2,) px/s - already computed

    r[2*i:2*i+2] = f.uv - (u_pred + td * u_dot)   # residual WITH td in it
    J[2*i:2*i+2, :15] = jac_pose_and_bias(...)  # unchanged
    J[2*i:2*i+2, -1]  = -u_dot                    # THE new column: dr/dtd = -u_dot

# observability gate: do not let an unexcited window move td
info_td = np.sum(u_dots**2) / sigma_px**2
if 1.0 / np.sqrt(info_td) > 2e-3:                 # worse than 2 ms - refuse
    J[:, -1] = 0.0

This is the mechanism behind Qin & Shen's online temporal calibration in VINS-Mono (IROS 2018) and, in filter form, Li & Mourikis (IJRR 2014). Being able to say "it is one Jacobian column, minus the image velocity" is the difference between having read about it and understanding it.

Two failure modes and the metric that reveals each

Failure mode A: t_d estimated with no excitation.

Failure mode B: the offset is really a drift.

The frontier, briefly

• Continuous-time calibration. Furgale, Rehder & Siegwart, "Unified Temporal and Spatial Calibration for Multi-Sensor Systems," IROS 2013 — the Kalibr paper — represents the trajectory as a B-spline so that every sensor can be evaluated at its own true timestamp without interpolating poses. This is the right formulation when your sensors have genuinely different rates, and it is the tool most teams actually run. Follow-up: Rehder & Siegwart, "Camera/IMU calibration revisited," IEEE Sensors Journal 2017, which is honest about what the original got wrong.
• Online in the filter. Li & Mourikis, "Online temporal calibration for camera-IMU systems," IJRR 2014 (MSCKF-style) and Qin & Shen, "Online Temporal Calibration for Monocular Visual-Inertial Systems," IROS 2018 (optimisation-based) are the two references to name. Both are the one-Jacobian-column trick above.
• Jointly with everything else. Taylor & Nieto, ICRA 2015 recovers extrinsics and timing offsets from ordinary motion, and Lv et al., T-RO 2022 adds an observability-aware selection of which trajectory segments to use. The direction of travel is clear: calibration is becoming a continuously-running estimation problem rather than a factory procedure.

Practice

Chapter 4: Reading the Residuals — Symptoms of Bad Calibration in a Running System

Picture a scatter plot with no axis labels and about two thousand points: the reprojection residual from a robot in the field. One calibration parameter is wrong. You get to pick what goes on the x axis. What do you pick, and in what order?

This is the best puzzle in the whole topic, because it separates people who have read about calibration from people who have debugged it. There is a right answer and it is short.

Every fault has a covariate

A residual is a number. A residual plotted against the right thing is a diagnosis. The reason this works is that each calibration parameter enters the projection multiplied by a different physical quantity, so an error in it produces a residual proportional to that quantity — and to nothing else.

Here is the whole table. It is how the knowledge gets used — but do not commit it to memory yet. Every row is one derivative of the projection map, and the section straight after it takes all five in front of you. A slope you have differentiated is a slope you can defend when someone asks "why r³ and not r²?"; a slope you have memorised is one you will misquote under pressure.

Where the table comes from — five derivatives

Nobody should carry that table on trust, least of all in a room where the next question is "why?". Every row is one partial derivative of the projection map, each takes about three lines, and the payoff is direct: the derivative you are about to write is the slope that the regression later in this chapter reports. Differentiate once and the triage procedure stops being a ritual and becomes a consequence.

Everything starts from the two lines Chapter 1 built. A point sitting at (X, Y, Z) in the camera frame lands at

Each fault perturbs exactly one symbol in that map. Perturb it, keep first order, and read off what multiplies the error. That multiplier is the covariate. There is no more to the theory than this sentence; the rest is five substitutions.

(a) Extrinsic translation δt → residual ∝ 1/Z. A lever-arm error means the camera physically sits δt away from where the extrinsic file says it does, so every point's camera-frame X is wrong by δt_x. Nothing else in the map moves — not f, not c_x, not the row, not the clock. So differentiate u with respect to X alone:

Read the right-hand grouping literally: the residual is a constant times 1/Z. Plot Δu against 1/Z and you get a straight line through the origin whose slope is f·δt_x — one number carrying both the focal length and the millimetres you are hunting. With f = 800 px and δt_x = 2 cm the slope is 800 × 0.02 = 16.0, in units of pixels per inverse metre, which is where the strange "px·m" in the threshold table comes from. Two sections down we cash that out arithmetically.

Notice what is absent: no dependence on where in the image the feature sits, and none on how fast it is moving. A lever-arm error is nearly invisible on a far-field-only log (0.8 px at 20 m) and screams on a near-field one (16 px at 1 m), which is why this test is only honest when the log carries range diversity. Chapter 2 tabulated exactly this fall-off; here it is the derivative that produced the table.

(b) Extrinsic rotation δθ → residual ∝ 1, a pure constant. Now rotate the camera by δθ about its Y axis instead of translating it. A rotation acts on the ray direction, and the normalised coordinate x = X/Z is precisely the tangent of the bearing angle θ. So x = tanθ becomes tan(θ + δθ), and

At the principal point x = 0 and the whole expression collapses to Δu = f·δθ. There is no Z anywhere in it: a rotation moves the point at 1 m and the point at 50 m by the same number of pixels. With f = 800 px and δθ = 1.0° = 0.017453 rad, Δu = 800 × 0.017453 = 13.96 px — the constant the practice widget at the bottom of this chapter injects for its "extrinsic Δθ" fault, and the same 13.96 px that runs down every row of the Chapter 2 comparison table.

The (1 + x²) factor deserves one more sentence, because it is a trap waiting in real data. At the image corner of a 1280×960 sensor with f = 800, x ≈ 0.5, so sec²θ = 1.25 and the same 1° error produces 800 × 0.017453 × 1.25 = 17.45 px rather than 13.96 px. So a large rotation error leaves a weak, genuinely radial residual — about 3.5 px of spread across the field here — sitting on top of its constant, and it can be misread as a mild distortion error. The discriminator is the value at r = 0: distortion is exactly zero at the principal point, a rotation is 13.96 px there.

(c) Radial distortion δk₁ → residual ∝ r³ (and δk₂ ∝ r⁵, δk₃ ∝ r⁷). This is the row you cannot reach from Chapter 1's formula without a change of units, and it is where most people fumble, because Brown–Conrady is written in normalised coordinates while the residual you plot is in pixels. Write r_n = r_px/f. The radial model pushes a point along its own radius from r_n to r_n(1 + k₁r_n² + k₂r_n⁴ + k₃r_n⁶), so the displacement is

Multiply by f to convert to pixels, then substitute r_n = r_px/f and let the powers of f collect:

Look hard at the denominators, because they are the whole reason the pixel-space powers are usable at all. Each extra order in r costs two more powers of f, and f is 800. Without the f² you cannot get from the normalised formula in Chapter 1 to the r³ row of the table, and the numbers come out absurd — k₁r³ alone at r = 400 would be eighteen million pixels.

Check it against a number you already have. Chapter 1 quoted k₁ = −0.281 producing a 28 px inward pull at r = 400 px on an f = 800 lens, derived there in normalised coordinates. The pixel form must agree:

It does, to the digit. That 400³/800² = 100 is worth remembering as a sanity anchor: on an 800 px lens, a distortion coefficient converts to pixels at r = 400 by multiplying by 100.

Now the diagnostic version, which is different in an important way. You are never plotting k₁; you are plotting the residual left by the error δk₁ that survived calibration. Take a 5.0% error on that lens, δk₁ = 0.0140, and evaluate at r = 460 px (roughly the corner radius of a 1280×960 frame measured from centre, and the largest radius the practice widget samples):

Two pixels at the corner from a 5% coefficient error — seven times a 0.30 px noise floor, and zero at the centre. That contrast is the entire diagnostic. The practice widget injects a deliberately gross version: 12 px at r = 400, which by the same formula is δk₁ = 12 × 640,000/64,000,000 = 0.120, a 43% error, chosen so the cubic shape is legible on a phone screen rather than realistic.

And the tangential row falls out of the same substitution. The p-terms in Chapter 1's model are 2p₁xy and p₂(r_n² + 2x²): quadratic in normalised coordinates, so in pixels they scale as p·r_px²/f — one power of f in the denominator, not two. That is the r² row. But they carry x and y separately rather than only through r, so the offset is not radially symmetric: it is largest along one diagonal and near zero along the other. That is the "coloured by angle" instruction in the table — a plain residual-vs-radius plot smears a tangential error into a fat band instead of a curve, and only colouring by azimuth separates it.

(d) Time offset t_d → residual ∝ image velocity u̇. The camera reports an image taken at true time t but stamped t − t_d, so the estimator predicts the feature using the pose from the wrong instant. Over that interval the feature has swept across the sensor at its image velocity u̇, so it lands t_d seconds' worth of motion away from the prediction:

The unit check is not decoration — it is the fastest way to remember which way the fraction goes, and it is why the slope of a residual-vs-velocity regression comes out directly in seconds with no focal length in sight. On the Chapter 5 rig (f = 800 px, yaw 0.30 rad/s, so u̇ = f·ω = 240 px/s) an 18.4 ms offset produces 240 × 0.0184 = 4.42 px. Stand still and it produces exactly nothing, which is the observability collapse Chapter 3 spent its length on.

(e) Line delay t_line → residual ∝ u̇ × row. A rolling shutter exposes image row k a delay of k·t_line after row 0. That is the same derivation as (d) with t_d replaced by k·t_line:

Two covariates multiplied together, which is why this row is the only one in the table whose "plot to draw" column has a division in it. And here is the number that makes the confusion concrete. Take the widget's t_line = 20 µs/row. At the bottom of a 960-row sensor, k = 900:

The widget's separate time-offset fault is 18.4 ms. So at the bottom of the image the two faults are numerically indistinguishable — 18.0 ms of effective lag versus 18.4 ms — while at the top row rolling shutter contributes zero and the time offset still contributes 18.4 ms. Nothing about them differs except how the lag varies with row. That is the derivation-level reason the fourth regression divides the residual by u̇ before looking at row index, and it is a far better answer than "because rolling shutter also depends on velocity."

Five derivatives, six rows of the table. The remaining two rows are already derived elsewhere in this lesson: the principal-point row is Worked Example 1 immediately below, and the gauge row is Chapter 0's argument that a uniform scaling of the world model leaves every projection exactly where it was. Here is the whole thing compressed to what each derivative multiplies:

Worked example 1 — a 10-pixel principal-point error, priced out

Why is c_x nearly invisible? Because it changes where you think the optical axis is, and that is almost the same as saying the camera is pointed slightly differently.

Step 1 — the bearing bias. A pixel at u maps to the ray direction (u − c_x)/f. If c_x is 10 px too large, every ray direction is shifted by

and crucially every ray shifts by the same amount, in the same direction. That is exactly what a small rotation of the camera does. Since a single camera cannot observe its own orientation, nothing in the reprojection residual complains — the estimator simply reports a slightly rotated camera pose, and all the residuals stay small.

Step 2 — where it goes. Visual odometry estimates the direction of translation from the image. A 0.7162° boresight bias rotates that direction by 0.7162°. Over 100 m of straight driving:

Step 3 — what it looks like in the field. The robot believes it drove straight down the aisle; it actually curved 1.25 m to one side over 100 m. Loop closure will fix the map, so the map looks fine. Live localisation between loop closures will be off by up to a metre, and the failure will present as "the robot clips the shelf on the long straight, but only on the eastbound run." The direction-dependence is the tell: a boresight bias is fixed in the camera frame, so it produces an error that flips sign when you drive the other way.

Worked example 2 — fitting the covariate, by hand, with a t-statistic

"The residual grows with image velocity" is an impression. "The slope is 18.5 ms with a t-statistic of 339" is a diagnosis. Here is the whole computation on five points, which is all you need to do at a whiteboard.

You have binned your residuals by image velocity and computed the mean residual in each bin:

Step 1 — fit through the origin. A time-offset error has no intercept (zero motion, zero residual), so fit r = b·u̇ with no constant term. The least-squares slope is

Numerator: 0×0.05 + 200×3.75 + 400×7.30 + 600×11.10 + 800×14.80
= 0 + 750 + 2920 + 6660 + 11840 = 22,170.
Denominator: 0 + 40,000 + 160,000 + 360,000 + 640,000 = 1,200,000.
b = 22,170 / 1,200,000 = 0.0184750 s = 18.475 ms.

Step 2 — the fit residuals. Predicted values are 0.018475 × u̇:
0,  3.695,  7.390,  11.085,  14.780.
Errors: 0.050,  0.055,  −0.090,  0.015,  0.020.
Sum of squares: 0.0025 + 0.003025 + 0.0081 + 0.000225 + 0.0004 = 0.014250.

Step 3 — the residual standard deviation. One fitted parameter, so n − 1 = 4 degrees of freedom:

Step 4 — the standard error of the slope. This is the one formula worth pushing on, because it is what turns a slope into a diagnosis, so derive it rather than quote it. It takes three lines.

4a. The estimator is linear in the data. Rewrite the slope with the denominator pulled inside the sum:

The weights w_i are built only from the covariates, which we treat as known. So b is a fixed weighted sum of the observations — and the variance of a weighted sum of independent numbers is the sum of the squared weights times their variance.

4b. Propagate the variance. If the fit residuals are independent with common variance σ²:

One factor of the sum cancels against the square, and what is left is the whole formula: the more your covariate varied, the better you know the slope. That is the same statement as Chapter 3's observability argument, arrived at from the algebra instead of from the information matrix.

4c. Substitute s for σ and put the numbers in. We already have s² = 0.0035625 from Step 3 and ∑u̇² = 1,200,000 from Step 1:

Step 5 — the t-statistic.

A t of 339 on 4 degrees of freedom is not a hint; it is a certainty. And the estimate itself, 18.475 ms, is your first correction — you have diagnosed and measured the fault from binned residuals you were already logging.

The intercept is not a nuisance — it is a fifth channel

I fitted that example through the origin, which is legitimate here because the intercept in this data is 0.05 px, i.e. nothing. In production, always fit the intercept and always read it. The reason is structural.

Look back at the table: every fault except one produces a residual proportional to some covariate. The exception is an extrinsic rotation error, which produces a constant. If you fit through the origin, that constant has nowhere to go, so it leaks into the slope of every regression you run — and you get four confident false positives at once. Fit an intercept and the constant lands where it belongs.

So the reading is two-dimensional, and the two dimensions are independent:

The four regressions, as a triage procedure

Run all four every time. Each is one ordinary least-squares fit over the same residual buffer, so the whole thing is a handful of dot products.

Where "px·m" comes from, in three lines of arithmetic. Regression 3 is the one whose units people cannot reproduce under pressure, and the threshold table below quotes alarms in px·m, so do it once by hand on the chapter's own numbers — f = 800 px, a Δt = 2 cm lever-arm error, which are exactly the constants the practice widget injects. From derivative (a), the residual is f·Δt/Z:

Regress the residual (px) on the covariate (m⁻¹) through those two points:

A pixel divided by an inverse metre is a pixel-metre. That is the whole mystery. And now divide by the focal length, watching the units cancel:

which is the fault we injected, recovered from two binned residuals. So the thresholds in the health-monitor table are not conventions — they are lengths in disguise. The alarm at 4 px·m is 4/800 = 0.005 m = 5 mm of lever-arm error; the warn at 1 px·m is 1/800 = 1.25 mm. An equivalent phrasing that sticks better: 1 px·m means one pixel of residual on a feature one metre away, because at Z = 1 the covariate is 1 and the slope is the residual.

Why regression 4 divides by u̇. Rolling shutter's residual is u̇ × k × t_line — it carries both covariates. Left alone it shows a strong slope against image velocity and looks exactly like a time offset. Divide the residual by u̇ first and the velocity dependence cancels, leaving t_d + k·t_line: a genuine time offset becomes a flat line at height t_d, and rolling shutter becomes a ramp with slope t_line. Same data, one division, two faults separated.

Separating r³ from r⁵ — the only row of the triage table that needs a second column

Regression 1's payoff line says "the fitted power tells you which order," and that sentence hides a problem: regress() above takes one covariate. It cannot fit a power. Neither should it try — fitting an exponent means a nonlinear solve, a starting guess, and a number like "2.87" that you then have to argue is really 3. There is a better move, and it is the one to build into the monitor.

Concretely: stack X = [r³, r⁵], solve the 2×2 normal equations X^TXβ = X^Ty by Cramer's rule, and get a t for each. Two practical points before the arithmetic. First, no intercept in this one: radial distortion is exactly zero at r = 0, so an intercept here would only soak up an extrinsic-rotation constant that belongs to the intercept channel — run the rotation check first, subtract its constant, then fit this. Second, normalise the covariates. Raw r⁵ at r = 460 is 2×10¹³; in float32 the normal equations would be numerically dead. Divide by r_max so the columns are ρ³ and ρ⁵ with ρ = r/r_max ∈ (0, 1], and the coefficients come out in pixels of offset at the corner, which is the unit you want to report anyway.

The five-point example. Synthetic k₁-only data on the lesson's lens: f = 800 px, a 5.0% coefficient error δk₁ = 0.0140, sampled at r = 100, 200, 300, 400, 460 px. From derivative (c), the true offsets are δk₁·r³/f² = r³/45,714,286, so:

The observed row is the true row plus about ±0.0013 px of noise — realistic for a bin mean, because averaging ~3,000 features at a 0.17 px per-feature σ gives 0.17/√3000 = 0.003 px of standard error on the mean.

Step 1 — the two normalised columns, ρ = r/460, c₃ = ρ³, c₅ = ρ⁵:

Step 2 — five dot products. That is the entire fit; there is nothing else to compute.

Step 3 — Cramer's rule on the 2×2. The determinant first, because it is also the collinearity warning:

Note how nearly the two products cancel — that is 3.2% of S₃₃S₅₅ surviving, and it is the first sign of what is coming in the trap below. Now the coefficients:

Step 4 — residuals and the noise estimate. Fitted values b₃c₃ + b₅c₅ are 0.021872, 0.174981, 0.590599, 1.400065, 2.129468, so the errors are +0.001028, −0.001181, +0.001401, −0.000965, +0.000332. Five points, two fitted parameters, three degrees of freedom:

Step 5 — the two standard errors. For a two-column fit the covariance is s²(X^TX)^-1, and inverting a 2×2 swaps the diagonal:

Step 6 — the verdict, which is the whole point:

Two numbers, one sentence: "the cubic term is 2.13 px at the corner with t = 345, the quintic is zero to within ±0.02 px, so this is k₁ and not k₂." And invert the cubic coefficient to get the parameter itself, using derivative (c) backwards:

which recovers the 0.0140 we injected to four figures. That is the difference between "distortion looks wrong" and "k₁ is high by 0.0140, which is 5.0% — re-run the calibration, do not touch k₂."

So what do you do when it does not? Not fit the power anyway, and not report the larger coefficient as the answer. You do three things, in this order. (1) Report the pair jointly: "the radial residual is 2.1 px at the corner; over this log's radius coverage r³ and r⁵ are 99.95% collinear, so I can state the magnitude but not attribute the order." That sentence is a stronger report than a confident wrong exponent. (2) Publish the condition number or the column correlation next to the coefficients, so the dashboard cannot silently degrade into over-claiming when a robot spends a week working close to a wall. (3) Fix it with data, not with statistics — the radius coverage is a collection problem, exactly like the tilt-diversity problem in Chapter 1 and the motion-diversity problem in Chapter 3. Gate the regression on observed radius span the same way you gate the time-offset regression on median image speed, and refuse to publish an order when the span is short.

The calibration health monitor

All of this should run on the robot, continuously, and cost nothing. Here is the shape of it with real sizes.

The thresholds, and where each comes from:

The regressor

Four lines of numpy that will save you more field time than any other four lines in this lesson:

python
import numpy as np

def regress(cov, res):
    """Fit res = a + b*cov. Returns (a, b, t_of_b, n)."""
    x = np.asarray(cov, float); y = np.asarray(res, float)
    xc = x - x.mean(); yc = y - y.mean()
    Sxx = float(xc @ xc)
    if Sxx < 1e-18:                          # covariate never varied -> no information
        return float(y.mean()), 0.0, 0.0, y.size
    b = float(xc @ yc) / Sxx
    a = float(y.mean() - b * x.mean())
    e = y - a - b * x
    s = np.sqrt(float(e @ e) / max(y.size - 2, 1))
    return a, b, b / (s / np.sqrt(Sxx) + 1e-30), y.size

# the four channels. note regression 4 divides the residual by image speed,
# which is what separates rolling shutter from a plain time offset.
CHANNELS = [("radius",    RAD,      RES),
            ("img_speed", UDOT,     RES),
            ("inv_range", 1.0/Z,   RES),
            ("row",       ROW,      RES / np.maximum(UDOT, 20.0))]
for name, cov, y in CHANNELS:
    a, b, t, n = regress(cov, y)
    print("%-10s intercept %+8.3f  slope %+.6g  t %+8.1f  n %d" % (name, a, b, t, n))

Say while typing: "Two things I would defend here. First, the intercept — every one of these faults except an extrinsic rotation is proportional to a covariate, so if I force the fit through the origin the rotation's constant leaks into all four slopes and I get four false positives. The intercept is its own diagnostic channel. Second, the Sxx guard: if the robot never moved, the covariate is constant, and I want the function to report 'no information' rather than divide by nearly zero and publish a confident nonsense slope to a dashboard someone will act on."

The library form is scipy.stats.linregress or numpy.polyfit, and I would use them for exploration — but the shipped monitor is these four lines, because it has to run on the robot with no scipy and a fixed memory budget.

And the one extension it needs, because regression 1 promised to name the distortion order and a single-covariate fit cannot do that. This is the Cramer's-rule arithmetic from the r³-versus-r⁵ section, transcribed, plus the collinearity guard that section argued for:

python
def regress2(c1, c2, y):
    """Through-origin fit y = b1*c1 + b2*c2. Returns (b1, b2, t1, t2, corr)."""
    S11 = c1 @ c1; S12 = c1 @ c2; S22 = c2 @ c2
    S1y = c1 @ y;  S2y = c2 @ y
    det = S11 * S22 - S12 * S12
    corr = S12 / np.sqrt(S11 * S22)          # publish this NEXT TO the coefficients
    if det <= 1e-12 * S11 * S22:               # columns collinear -> refuse to attribute
        return None, None, 0.0, 0.0, corr
    b1 = (S22 * S1y - S12 * S2y) / det
    b2 = (S11 * S2y - S12 * S1y) / det
    e  = y - b1 * c1 - b2 * c2
    s2 = (e @ e) / max(y.size - 2, 1)
    return b1, b2, b1 / np.sqrt(s2 * S22 / det), b2 / np.sqrt(s2 * S11 / det), corr

# normalise so the columns are O(1) and the coefficients are px AT THE CORNER.
rho = RAD / RAD.max()
b3, b5, t3, t5, corr = regress2(rho**3, rho**5, RES)
# -> 2.1289  0.0006  345.1  0.087  0.984   (k1 is wrong; k2 is not)
# dk1 = b3 * f**2 / RAD.max()**3   -> 0.0140

The reasoning behind those two lines: "I am deliberately not fitting the exponent. The candidate orders are 3, 5 and 7 because that is what the lens model offers, so I put them in as separate columns and let each one carry its own t — that stays linear and closed-form, and it hands me a magnitude per order instead of an exponent I would then have to argue about. The corr return is not diagnostics decoration: over a narrow radius band these two columns run to 0.999 correlated, and I want the monitor to say 'magnitude yes, order no' rather than pick whichever coefficient the noise favoured."

Three named modes

Mode A: healthy residuals, wrong world (the gauge).

Mode B: the slow thermal walk.

Mode C: one bad sensor hidden by the aggregate.

The frontier, briefly

• Self-calibration inside the estimator. Geneva et al., "OpenVINS," ICRA 2020 makes the camera intrinsics, the camera–IMU extrinsic and t_d optional filter states with published covariances — so the monitor in this chapter becomes "watch the filter's own uncertainty" instead of a separate regression. The catch is the one from Chapter 3: a state that absorbs faults hides them, so you still want an independent reference.
• Long-lived maps as the reference. Schneider et al., "maplab," RA-L 2018 treats a persistent map as the calibration ground truth: localise against yesterday's map, and any systematic residual is today's calibration drift. This is the practically strongest answer for a fleet, because it needs no extra hardware and gets better as the map matures.
• Observability-driven monitoring. The idea from Lv et al., T-RO 2022 generalises: rather than thresholding residuals, track the Fisher information your data carries about each parameter and alarm when a parameter has been unobservable for too long. It reframes calibration monitoring as "am I still learning about this?" instead of "is this still small?", which is the right question.

Practice

Chapter 5: The Fault Injector — Build the Table by Breaking Things

Everything so far has been derivations and tables. Tables are how the knowledge is stored; they are a terrible way to learn it, because a table read is forgotten by Thursday and a table you built from experiments is not.

So here is the rig, and here are seven ways to break it. Turn one knob at a time and watch which of the six downstream signals moves. By the time you have been through all seven you will have the symptom-to-cause table in your head, derived rather than memorised — and then the blind drill at the bottom will tell you whether that is true.

This is a running monitor, not a geometry exercise

Before any of the physics, be able to say what this thing is in a system. "We project the LiDAR into the image and look at the residual" immediately raises four questions — at what rate, on what buffer, in what array, and inside what budget — and those four answers are the difference between a monitor and a Jupyter notebook.

The monitor is an online gate: it runs at 1 Hz on a 3-second ring buffer and has 40 ms of one core to produce six numbers. Here is every boundary it crosses, with the shape and the byte count on each one.

The 40 ms, spent. Transform and projection over 432,000 points at ~20 flops each is 8.6 Mflop, about 4 ms. The distance-transform lookup is one cache-hostile bilinear fetch per point at ~50 ns, which is 22 ms and is the hot spot. The six reductions over 48,000 survivors are ~1 ms. The two OLS fits are over six aggregated numbers and are free. Typical tick 27 ms, p99 40 ms — 4% of one core at 1 Hz, which is what makes it affordable to leave running on every robot forever.

The scene, and why each structure has its own image velocity

Chapter 3 established that a time offset is only observable where u̇ ≠ 0, and that you estimate it by regressing residual on image velocity. A regression needs spread in its regressor. If every structure in the scene sweeps at the same 240 px/s, the velocity column has zero variance, S_xx = 0, and the slope is undefined — you can no longer measure t_d, only assume it.

So the monitor's 3-second window is taken during a slalom, where the yaw rate is deliberately not constant. Each structure's residual sample is tagged with the image velocity it actually had:

The lateral offset of each structure in metres follows from X = rZ/f, so the six are at 0.10, 0.40, 1.50, 4.00, 10.00 and 28.75 m off the optical axis — worth having, because the depth-squared term in the derivation above is fX/Z², and fX/Z² = r/Z.

Structures 2 and 5 — the two the colouring gauges read — were both sampled at 460 px/s. That is deliberate too: it means a pure time offset produces equal colouring at 2 m and 20 m, so it perfectly imitates a boresight error on the ratio test and can only be unmasked by the velocity column. The bench would be easier and less honest without it.

Those six velocities are not arbitrary. They were chosen so that the velocity column is nearly orthogonal to the other three regressors the bench cares about — r, r³ and 1/Z. Orthogonal means the mean-centred dot product is zero, so a fault that lives in one regressor leaks nothing into the others. Check the worst of the three by hand, u̇ against r:

Compare that to the scale of the two columns: S_uu = 60²+160²+80²+120²+160²+60² = 79,200 and S_rr = 104,333. The correlation is −400 / √(79,200 × 104,333) = −400 / 90,893 = −0.0044. Four parts in a thousand. That is why, in the experiments below, a pure time offset leaves the radius slope at −0.007 px/100 px — visually zero — and a pure focal error leaves the velocity slope at −0.10 ms.

The six signals on the dashboard

Derive the model the bench is built from

The bench is not evidence. It is a simulator, and a simulator cannot prove the equations it was coded from — if you slide the extrinsic-rotation knob and see a depth-independent offset, all you have proved is that somebody typed a depth-independent expression. So derive the expressions first, and then use the bench for what a bench is actually good at: turning three symbolic terms into an experience you will still have in six months.

Take one LiDAR point p_L. The correct extrinsic puts it in the camera frame as p = (X, Y, Z) = R p_L + t. Now suppose the calibration you shipped is wrong by a small rotation δθ and a small translation δt. A small rotation is I + [δθ]_× to first order, so the point you actually compute is

So the 3-D perturbation is δp = δθ × p + δt. Write out the two components that matter for the horizontal pixel coordinate, using the cross-product definition:

Now push it through the projection u = f·X/Z + c_x. This is where the depth structure appears, and it appears because of the quotient rule and nothing else:

Substitute the two lines above and expand every product:

Six terms. Group them by what they do to the depth axis, keeping only the paraxial ones (X, Y small next to Z, which is the regime the bench models):

Read the three terms and the whole chapter is already in your hands:

• Term 1, f·δθ_y. No Z anywhere. A boresight error produces the same pixel offset on a wall at 2 m and a wall at 50 m. Experiment 2 isolates this term.
• Term 2, f·δt_x/Z. Exactly 1/Z. A lever-arm error is enormous up close and vanishes far away. Experiment 3 isolates this term. Terms 1 and 2 are what makes the near/far ratio the extrinsic discriminator: ratio 1 means term 1, ratio equal to the range ratio means term 2.
• Term 3, −(fX/Z²)·δt_z. A range-direction lever-arm error. Note fX/Z² = (fX/Z)/Z = r/Z, so this term is r·δt_z/Z — simultaneously linear in radius and falling as 1/Z. It is the one fault that lights the radius gauge and the near/far ratio at the same time, and the bench does not model it (see the omissions below).

The intrinsic faults are not in that expansion because they act after projection, on the pixel itself. Chapter 1 derived them: a focal-scale error ε = Δf/f moves a pixel by ε·r (linear in radius), a principal-point error moves it by Δc_x (constant), and a radial distortion coefficient moves it by k₁r³/f² (cubic in radius, because the normalised radius is r/f and the Brown term is k₁·x·(r/f)²). Chapter 3 gave the temporal one: u̇·t_d. And Chapter 0 gave the gauge: zero residual, scale s/25.

That is the complete model behind every gauge below. Nothing in the widget is unexplained.

The seven experiments, and what each one proves

Do these in order. Reset between each. Each one takes fifteen seconds and settles a question that costs teams weeks.

Experiment 1 — board square size, 25.0 → 23.0 mm. Watch every gauge except map scale ratio stay exactly where it was. RMS: 0.28. Colouring near: 0.0. Colouring far: 0.0. Both slopes: 0.0. Map scale: 0.920.

Proves: a metric gauge has no residual signature. The five internal signals are blind by construction, because the reconstruction is perfectly self-consistent in the wrong units. Only the sixth signal, which requires an external reference, sees it. This is the Chapter 0 story, reproduced as an experiment.

Experiment 2 — extrinsic rotation, 0 → 1.00°. Colouring at 2 m goes to 13.96 px. Colouring at 20 m goes to 13.96 px. Identical.

Proves: f·Δθ has no depth in it. The two readings are identical because Z does not appear anywhere in the expression — not because the two structures happen to be similar. Two structures a decade apart in range, the same offset. Both slope gauges stay at exactly zero, because a constant is not a slope — it is an intercept, and OLS puts constants in the intercept. Note also that tan(0.0174533) = 0.0174551 differs from the angle itself only in the sixth decimal, which is why everyone writes f·Δθ and nobody notices; at 10° the small-angle version would be off by 1%.

Experiment 3 — extrinsic translation, 0 → 2.0 cm. Colouring at 2 m goes to 8.0 px. Colouring at 20 m goes to 0.8 px. A factor of exactly ten, for a factor of exactly ten in range.

Proves: f·Δt/Z. The numerator f·Δt = 16.0 px·m is the same for both structures; only the division by Z differs, so the ratio of the two readings is the inverse ratio of the two ranges, always, with no dependence on f or on Δt. Put experiments 2 and 3 side by side and you have the entire LiDAR-camera extrinsic diagnostic: the ratio of near offset to far offset is the answer. Ratio 1 means rotation; ratio 10 means translation. And because the ratio is independent of the fault magnitude, it works on a fault too small to alarm the RMS gauge. Now look at the radius slope, which also moved, to −3.81 px per 100 px. Nothing radial is broken. That reading is a confound — in this scene the far structures are also the peripheral ones, so a 1/Z residual is automatically a decreasing-with-r residual too. It is the single most common way a competent engineer misreads this dashboard, and the worked example at the end of the chapter takes it apart properly.

Experiment 4 — time offset, 0 → 18.4 ms. RMS jumps to 5.92 px. The velocity slope reads 18.4 ms. Both colouring gauges move by the same amount — not because velocity is uniform (it is not), but because structures 2 and 5 were both sampled at 460 px/s.

Now do the regression by hand, because this is the gauge whose whole claim is that it measures rather than restates. Deviations of u̇ from ū = 300 are (−60, +160, −80, −120, +160, −60); deviations of the offsets from their mean 5.520 are (−1.104, +2.944, −1.472, −2.208, +2.944, −1.104):

And the RMS, so you can check the top gauge too: mean of the six squared offsets is (4.416² + 8.464² + 4.048² + 3.312² + 8.464² + 4.416²)/6 = 210.03/6 = 35.005, and RMS = √(0.28² + 35.005) = √35.083 = 5.92 px.

Proves: the velocity slope reads out the fault in its own units. You do not infer a time offset from this plot; you measure it, and the number that comes out of Sxy/Sxx is the number you put in, to five digits. Note that colouring near and far both moved equally (8.464 and 8.464), so it masquerades as a rotation error on the ratio test alone — and the thing that unmasks it is the velocity column, which is why the scene was built to have one. Check the cross-talk while you are here: the radius slope reads −0.007 px/100 px, dead flat, because u̇ was designed orthogonal to r.

Experiment 5 — focal length, 0 → +2.0%. RMS climbs to 6.14 px. The radius slope goes to 2.0 px per 100 px of radius. Colouring at 2 m and 20 m both move, but by amounts set by where in the image the structure is, not by its range.

That middle line deserves a beat: because the offset is exactly ε times the regressor, the OLS slope is ε identically — S_re = ∑(r_i−r̄)(εr_i−εr̄) = ε·S_rr, and the S_rr cancels. The gauge is not correlated with the fault; it is the fault, in units of fractional focal error.

Proves: an intrinsic scale error is radial and linear: offset = ε·r. This is why the LiDAR matters — without a metric reference the camera would happily absorb a wrong focal length into a rescaled reconstruction and show you nothing. Cross-talk check: the velocity slope reads −0.10 ms, inside the 1 ms warn band, so the temporal gauge stays quiet.

Experiment 6 — distortion k₁, 0 → +0.02. The radius slope moves to 0.78 px per 100 px, but look at the shape in the dashboard's scene strip: the offset is nearly nothing at 160 px of radius and large at 460 px.

Proves: distortion is cubic in radius where focal length is linear. Both light up the same gauge; the power separates them, and the separation is enormous — a factor of 23.8 against a factor of 2.9 between the same two structures. That is the whole discriminator, and it is why you fit the residual against a design matrix with both an r column and an r³ column and compare the two t-statistics, rather than fitting "radius" and squinting. The code block below does exactly that, because a rule you cannot implement is a rule you will misapply.

Experiment 7 — principal point, 0 → +10 px. Colouring moves by 10 px everywhere. Both slopes stay at zero. RMS climbs to 10.00 px.

Proves: to first order a principal-point error is a constant image shift — indistinguishable, from these six signals alone, from an extrinsic rotation of 0.7162°. Set the c_x slider to +10 and note the six readings; reset, set the Δθ slider to 0.72°, and the six readings are the same to two decimals. Not similar. The same. That degeneracy is real, it is why the two are strongly correlated in any joint calibration, and it has a literature (see Frontier below). To separate them you need a signal these gauges do not have: a second camera (the principal-point error is per-camera, the extrinsic is per-pair), the (1 + (r/f)²) off-axis term from the derivation above if your lens is wide enough to make it measurable, or the trajectory heading against a survey, as in Chapter 4.

Experiment 8 — turn the noise on, and watch the diagnosis rot

Experiments 1–7 were run at σ = 0, which is a lie of omission every bench tells. In the field each gauge is an estimate, and an estimate without an interval is a rumour. The seventh slider adds Gaussian edge-localisation noise of standard deviation σ to each structure's residual sample, and the widget starts printing ±1.96 SE next to every reading and a whisker on every bar.

The first thing to be honest about is the sample size, because this is where most monitoring dashboards lie to themselves. Roughly 48,000 residuals reach the gauges every tick, so the naive standard error of a per-structure mean is σ/√8000 = σ/89.4 — at σ = 1 px that is 0.011 px, and you would conclude that every fault above a hundredth of a pixel is detectable. It is not, because those 48,000 numbers are nowhere near independent. Edge localisation on one wall segment shares an extraction bias across every point that lands on it; the correlated part does not average away. The honest effective sample size is a handful of independent error draws per structure per window:

Now re-run Experiment 3 (extrinsic translation, 2.0 cm) with σ = 1.0 and read the far gauge:

And that is the whole failure, in one line. The near/far ratio test needs a far reading, and there is no longer one. Propagate the two intervals through the delta method and watch the ratio come apart:

An interval that contains 1 means the data cannot tell an extrinsic translation from an extrinsic rotation. The bench refuses to print a ratio at all in this state and says why — "near/far ratio undefined: the far gauge is consistent with zero" — because a dashboard that prints "10.2" here is lying by omission. An engineer reading only the point estimates from such a dashboard, on the next realisation of the same noise, diagnoses an extrinsic rotation on a robot whose lever arm is 2 cm out. Wrong fault, wrong fix, and a confident-sounding write-up.

Try it: set extrinsic Δt to 2.0 cm, then walk σ up from 0 and watch the verdict line. The far reading holds at 0.80 px the whole way; what changes is the whisker, 1.96·σ/2 = 0.98σ, and the moment that exceeds 0.80 — at σ = 0.85, one slider notch past 0.80 — the diagnosis stops being available. Nothing about the robot changed. Nothing about the fault changed. The instrument stopped being able to see it, and the only thing on the dashboard that told you so was the width of a whisker.

Now dial σ back to 0.35 px and the same bench recovers:

What actually buys precision here. Not more points — we just saw that 48,000 of them bought 4 independent draws. Three things move the interval: a longer buffer (n_eff grows with independent windows, so 30 s instead of 3 s is a factor of √10 = 3.2), a nearer far-field structure is useless but a farther near-field one helps because it widens the 1/Z lever, and better edge localisation lowers σ directly. "I would average more points" is the instinct that fails here; "I would extend the window, because the noise is spatially correlated and only independent windows count" is the one that works.

The slope gauges degrade too, and their standard errors are the classic OLS form SE(slope) = SE(point)/√S_xx. Both regressors were designed with real spread, so both survive:

Now imagine the scene had been sampled at a nearly constant yaw rate — velocities (240, 260, 240, 260, 240, 260) instead of the designed spread. Then S_xx = 6×10² = 600, √S_xx = 24.5, and SE(slope) = 0.50/24.5 = 20.4 ms. The reading would be 18.4 ± 40 ms — a t-statistic of 0.9, indistinguishable from no time offset at all. Same noise, same fault, same estimator; the scene design is what decided whether the fault was visible. That is Chapter 3's observability gate reappearing as a property of where you point the robot.

The diagnostic, from scratch

Everything above is a procedure a human runs by eye. The version that ships is forty lines, it runs once a second on every robot, and writing it is the fastest way to find out whether you actually understood the ratio test — because code will not let you wave at "and then you look at the slope."

Three pieces. An OLS slope by hand, as S_xy/S_xx, so nothing is hidden inside polyfit. A two-column radial fit against r and r³ simultaneously, which is the instruction Experiment 6 gave and which nobody ever implements. And a return type that can say "I do not know" — the degeneracy flag, without which the function will confidently label a principal-point error as a boresight error forever.

python
import numpy as np

F, N_EFF = 800.0, 4          # focal px; independent error draws per structure

def ols_slope(x, y):
    """Slope of y on x, by hand: Sxy / Sxx. No polyfit, no lstsq."""
    n  = len(x)
    mx = sum(x) / n                                   # x-bar
    my = sum(y) / n                                   # y-bar
    Sxy = sum((x[i] - mx) * (y[i] - my) for i in range(n))
    Sxx = sum((x[i] - mx) ** 2          for i in range(n))
    return Sxy / Sxx, Sxx                          # slope, and the lever we paid for it

def radial_powers(r, e, se):
    """Fit e = c + a*r + b*r**3 in ONE regression and return both t-stats.
       Linear-in-r is a focal-length error; cubic-in-r is distortion k1.
       Fitting them separately is the classic mistake: r and r**3 are
       correlated 0.94 over this radius range, so a univariate r-fit on a
       pure distortion fault comes back significant and you call it focal."""
    A    = np.column_stack([np.ones_like(r), r, r ** 3])   # (6,3) design matrix
    AtA  = A.T @ A                                        # (3,3) normal matrix
    beta = np.linalg.solve(AtA, A.T @ e)                 # (3,) = [c, a, b]
    cov  = np.linalg.inv(AtA) * se ** 2                  # Cov(e) = se**2 * I
    t    = beta[1:] / np.sqrt(np.diag(cov)[1:])          # (t_r, t_r3)
    return beta, t

def diagnose(offsets, ranges, radii, velocities, map_scale, sigma=0.35):
    """Six gauges in -> (fault_label, magnitude, degeneracy_flag) out.

    offsets    (6,) mean residual per structure          [px]
    ranges     (6,) Z per structure                      [m]
    radii      (6,) image radius per structure           [px]
    velocities (6,) image velocity per structure         [px/s]
    map_scale  scalar: surveyed distance / mapped distance
    """
    e = np.asarray(offsets,    float)
    r = np.asarray(radii,      float)
    Z = np.asarray(ranges,     float)
    u = np.asarray(velocities, float)
    se = sigma / np.sqrt(N_EFF)     # NOT sigma/sqrt(48000): the points are correlated

    # -- 0. the gauge channel: metric error, and NO residual signature -------
    if abs(map_scale - 1.0) > 0.002 and np.max(np.abs(e)) < 1.96 * se:
        return ('metric gauge (square size / wheel radius / baseline)',
                map_scale, False)

    # -- 1. temporal FIRST: it fakes a ratio of 1 and would be read as -------
    #       a boresight error if the ratio test ran before it.
    td, Suu = ols_slope(u, e)                  # units px / (px/s) = SECONDS
    t_td    = td / (se / np.sqrt(Suu))          # SE(slope) = SE(point)/sqrt(Sxx)
    if abs(t_td) > 3.0:
        return ('time offset', td, False)

    # -- 2. extrinsic: the near/far ratio, guarded by the far reading's CI ---
    near, far = e[1], e[4]                  # the 2 m and 20 m structures
    if abs(far) < 1.96 * se:                 # far gauge consistent with zero:
        return ('extrinsic, undetermined - far gauge is not measuring',
                None, True)                     # the ratio is unbounded. say so.
    ratio = near / far
    if abs(ratio - 1.0) < 0.12:               # constant in range
        dtheta = np.arctan(near / F)          # rad, IF it is boresight at all
        return ('extrinsic rotation OR principal point',
                dtheta, True)                     # <-- DEGENERATE. never collapse this.
    if ratio > 4.0:                            # falls as 1/Z
        return ('extrinsic translation', near * Z[1] / F, False)   # metres

    # -- 3. radial: which POWER of r carries the residual? -------------------
    beta, (t_lin, t_cub) = radial_powers(r, e, se)
    if max(abs(t_lin), abs(t_cub)) > 3.0:
        if abs(t_cub) > abs(t_lin):
            return ('distortion k1', beta[2] * F ** 2, False)   # k1 = b * f**2
        return ('focal length', beta[1], False)              # eps = a = df/f

    return ('no fault above noise', 0.0, False)

Run it on the seven experiments and this is what comes back, at σ = 0.35 px so SE = 0.175 px:

Every magnitude comes back to five digits, and the two degenerate cases come back flagged rather than guessed. But look at the translation row, because it is the one that teaches something: t_r = −52.2 and t_r³ = +30.6. A pure lever-arm error produces a screaming radial signature, far larger than the real focal-length fault's 12.3, purely because in this scene the far structures are also the peripheral ones and so r is correlated with 1/Z. If the radial branch ran before the ratio branch, this function would confidently report a focal-length error on every loose bolt in the fleet.

The library one-liner, which is what you would actually ship once you have earned the right to use it:

python
# One joint fit instead of four ordered tests. Columns, left to right:
#   1     -> constant shift   : dcx  OR  f*dtheta   (structurally degenerate)
#   r     -> focal scale      : eps = df/f
#   r**3  -> distortion       : k1 = coef * f**2
#   1/Z   -> lever arm        : dt = coef / f     [px*m -> m]
#   udot  -> time offset      : td = coef         [seconds, directly]
A = np.column_stack([np.ones(6), r, r**3, 1.0/Z, u])          # (6,5)
beta, *_ = np.linalg.lstsq(A, e, rcond=None)                   # (5,)
t = beta[1:] / np.sqrt(np.diag(np.linalg.inv(A.T @ A))[1:]) / se

On the same seven faults this recovers beta = [0, 0.02, 0, 0, 0] for the focal error, [0, 0, 3.125e-8, 0, 0] for distortion (times f² = 0.0200 = k₁), [0, 0, 0, 16.0, 0] for the lever arm (16.0/800 = 0.0200 m), and [0, 0, 0, 0, 0.0184] for the time offset — each fault in exactly one column, no ordering required. And the boresight and principal-point faults both land entirely in the intercept, 13.964 and 10.000, which is the cleanest possible statement of the degeneracy: they are not two hypotheses that are hard to tell apart; they are the same column of the design matrix. No estimator, however clever, separates them from this data.

The blind drill

Press Mystery fault. One knob moves by a random amount; the readouts turn to question marks. You get the six gauges and nothing else. Name the cause.

Two rules that make the drill worth doing:

1. Read the ratio before you read the magnitude. Beginners look at which bar is biggest. The near/far ratio and the slope signs carry the diagnosis; magnitudes only tell you how bad it is.
2. Say your reasoning out loud before you click. "Colouring near equals colouring far, both slopes zero, map scale one — so it is constant in range, constant in velocity, constant in radius. That is either an extrinsic rotation or a principal-point error, and I cannot separate them here, so I would say rotation and name the ambiguity." That sentence is worth more than a lucky correct guess.

What this bench does not model, and why that matters

A showcase that pretends to be complete teaches overconfidence. Four things this bench leaves out, each of which you should be able to name unprompted:

• Cross-terms. Real faults arrive in pairs and partially cancel. A focal-length error plus a compensating distortion error can hold RMS almost flat while both parameters are far from truth — which is exactly the over-parameterised-distortion failure from Chapter 1. Here the knobs are independent; in a joint calibration they are correlated, and the correlation matrix is the thing to look at. The joint design matrix above is the first step toward modelling this: it at least lets the columns compete for the same residual. What it still cannot do is represent a fault that lives outside its five columns.
• Noise and finite samples — now partly modelled. The σ slider and the ±1.96 SE whiskers close most of this gap, and Experiment 8 is where you feel it. What remains unmodelled is the shape of real residual noise: it is not Gaussian and not homoscedastic. Edge localisation is worse on low-contrast structures, worse at grazing incidence, and heavy-tailed wherever the association picked the wrong edge. A single σ is a convenient fiction; the honest version needs a robust cost (Huber, from Chapter 3) and a per-structure σ_i, which turns the OLS into a weighted least squares and the t-statistics into something you have to be more careful about.
• Time-varying faults. These knobs are constants. Thermal drift, clock drift and a loosening bolt all move slowly, and the diagnostic for those is the trend, not the value. A bench that only shows you steady state will not teach you to plot things against time, and plotting things against time is how the slow faults get caught.
• The depth-squared lever-arm term. The derivation at the top of this chapter produced three terms; the bench implements two. Its offset model is eps*r + dcx + k1*r**3/f**2 + f*tan(dtheta) + f*dt/Z + udot*td — the f*dt/Z is term 2, the lateral lever arm. Term 3, −(fX/Z²)·Δt_z = −(r/Z)·Δt_z, is simply absent. That matters because it is the only fault in the whole taxonomy that is simultaneously linear in radius and falling as 1/Z: it lights the radius slope like a focal error and drags the near/far ratio like a translation error, at the same time, and the six gauges have no way to name it. Concretely, a 2 cm error along the camera's optical axis gives 0.02 × 400/20 = 0.40 px at the 20 m structure and 0.02 × 160/2 = 1.60 px at the 2 m structure — ratio 4.0, right on the boundary of the ratio > 4 branch, which the code would report as a lateral translation of the wrong magnitude in the wrong axis. If someone hands you a dashboard where the radius slope and the near/far ratio both moved and neither the focal nor the translation magnitude cross-checks, this is the term to reach for. It is also the reason range-direction lever arms are the hardest extrinsic component to calibrate: they hide inside two other signatures.

Fault injection and calibration monitoring as a research area

Everything on this bench is a threshold on a residual, which is the 2015 answer. Three papers to name if the conversation goes past the bench, all specific to monitoring rather than to calibration itself:

• Rehder & Siegwart, "Camera/IMU Calibration Revisited," IEEE Sensors Journal, 2017. The paper that goes back over the Kalibr formulation and asks what the reported uncertainties actually mean. The lesson for this chapter is sharp: an injected fault very often shows up in the estimator's own reported covariance long before it shows up in RMS, because the optimiser absorbs the fault into a nearby parameter and pays for it with a widened marginal. If your monitor watches only the cost function you are watching the one number the optimiser is actively minimising. Watch the covariance diagonal and the parameter trend instead. Sit with the question "what does Kalibr's reported sigma on t_d mean" and you learn a great deal.
• Lv, Zuo, Hu et al., "Observability-Aware Intrinsic and Extrinsic Calibration of LiDAR-IMU Systems," IEEE Transactions on Robotics, 2022. Replaces "is the residual above 0.5 px" with "is this parameter observable in the data I just collected", by decomposing the Fisher information and selecting which trajectory segments to use. This is the principled version of the warn/alarm bands on the bench: rather than a fixed threshold on a residual, you gate on information, so a gauge that is uninformative in this window is silent rather than confidently wrong. It is the same idea as Experiment 8's confidence intervals, done properly and per-parameter.
• Yang, Geneva, Zuo, Eckenhoff, Liu & Huang, "Degenerate Motion Analysis for Aided INS with Online Spatial and Temporal Sensor Calibration," IEEE Robotics and Automation Letters, 2019. The citation for the degeneracy this chapter keeps running into. It enumerates the motion profiles under which the online calibration states become unobservable — constant velocity, pure translation, one-axis rotation — and shows exactly which parameter directions collapse in each. The c_x-versus-Δθ pairing from Experiment 7 is one instance of a general phenomenon, and being able to say "that is a known degenerate direction, it is characterised in Yang et al. 2019, and the fix is motion excitation not a better solver" is a much better answer than "they are correlated."

Turning the bench into a diagnosis procedure

When someone hands you a symptom, the bench in your head should run backwards. Here is the reverse lookup, which is the one you actually use:

Notice that every row's confirming question is cheap — it costs one plot or one sentence, not a recalibration campaign. Producing the cheap discriminating question, rather than the expensive experiment, is the habit that separates a diagnostician from a guesser.

One worked reverse lookup, out loud

Here is what the answer sounds like when you do it properly. Set the noise slider to σ = 0.35 px so the intervals are real, and suppose the bench hands you this. (The point estimates below are written noise-free so the arithmetic stays checkable; the intervals are exactly the ones σ = 0.35 produces, and a live realisation will jitter the last digit of each reading.)

Step 1 — ratio before magnitude. Near over far is 8.00 / 0.80 = 10.0. The two structures are at 2 m and 20 m, a factor of ten. So the offset is falling as exactly 1/Z. Nothing else in the fault set does that. And say the interval out loud, because it is what makes this a measurement: SE(ratio) = 10.0 × √((0.175/8)² + (0.175/0.8)²) = 2.20, so the ratio is 10.0 with a 95% interval of 5.7 to 14.3 — inconsistent with 1. Had the interval straddled 1, the correct next sentence would have been "this window cannot separate rotation from translation," not a guess.

Step 2 — rule out the other channels, and do not flinch at the radius slope. The velocity slope is −0.24 ± 1.22 ms, a t of 0.4, so it is not temporal. The map scale is 1.000, so it is not a gauge. The offset is not constant with range, so it is not an extrinsic rotation or a principal-point error. That leaves the radius slope, which is −3.81 px per 100 px at t = 70 — five times more significant than the real focal-length fault in Experiment 5 — and which anyone who has only memorised the table will now use to announce a focal-length error.

Step 2b — why that slope is a ghost, said in one sentence. In this scene the far structures are also the peripheral ones: r and 1/Z have a correlation of −0.91. So a residual that falls as 1/Z is automatically also a residual that falls with r, and the marginal regression on radius happily reports it. The tells are all there once you look: the sign is negative (a focal error grows outward, this shrinks outward), and the cubic coefficient came back significant too (t = +30.6), which no pure focal error ever does. The cure is the joint fit from the code block — put r, r³, 1/Z and u̇ in one design matrix and the entire residual lands in the 1/Z column at 16.0 px·m with the other three collapsing to zero.

Step 3 — invert for the magnitude. The model is Δu = f·Δt/Z. Take the 2 m reading:

Check it against the other reading: 0.80 × 20 / 800 = 0.0200 m. Both agree, which is itself confirmation that the 1/Z model is the right one — if the two disagreed, the fault would be a mixture. Attach the interval here too: ±0.34 px at 2 m propagates to ±0.34 × 2 / 800 = ±0.00085 m, so the honest statement is Δt = 20.0 ± 0.9 mm. A number without a millimetre is a number nobody can act on.

Step 4 — name the likely upstream cause and the confirming question. A 2 cm lever-arm error, almost always from a rotation-starved calibration run. "How much did the robot rotate during that calibration? I want the smallest singular value of the stacked (R_A − I) matrix — under 0.5 and I do not trust the translation at all, and at 2° of excitation a 1 cm motion error becomes 29 cm of lever arm."

Step 5 — state what would change your mind. The one thing this dashboard cannot rule out is the depth-direction lever arm from the omissions list, −(r/Z)·Δt_z, which also falls as 1/Z and grows with r. Its near/far ratio would be 4.0, not 10.0, so the measured 10.0 ± [5.7, 14.3] does not exclude it. "I would confirm by driving the same aisle at two different standoff distances: a lateral lever arm keeps the near/far ratio at the range ratio, a range-direction one does not." Naming the alternative you did not choose, and the cheap experiment that would separate them, is the difference between a diagnosis and an opinion.

Chapter 6: The Field Guide

One last dashboard, and then the reference material. A fleet dashboard from a warehouse-robot company: the maps come back 8 % short against the surveyed floor plan, the reprojection RMS is flat at 0.28 px and has not moved in three months, and every loop closure in the last week accepted on the first try. What do you check first?

On a live fleet, you have about ninety seconds of credibility before "let me get back to you" becomes the answer. This chapter is the set of instruments for exactly those ninety seconds — the one derivation you can rebuild from a blank page, the numbers you can do without a calculator, the discriminator that separates two hypotheses that fit the same plot, and the sentence that closes the investigation. It is a reference; do not read it instead of the chapters.

0 · The one derivation to keep in muscle memory

There is one derivation at the centre of camera calibration, and it is always the same one: why three views of a plane are enough to recover K. It rewards practice because it is short enough to finish inside five minutes and structured enough that you cannot bluff it — one substitution, two facts about rotation matrices, and a counting argument. The tables further down this chapter are only useful to somebody who can also do this, so do it here first, out loud, until it is muscle.

Step 1 — put the frame on the board, and the third rotation column dies. Attach the world frame to the target so every corner has Z_board = 0. The full pinhole projection of a corner is

and that zero deletes r₃ before it ever multiplies anything, leaving a 3×3 map from board coordinates to pixels:

That 3×3 map is a homography, and four or more corner correspondences fit it directly from one image with no knowledge of K at all. But the fit only ever determines it up to an overall scale, because projection divides by the third coordinate and a common factor cancels. So the honest statement of what you have is

Step 2 — move K to the other side. Reading off the columns, and writing h₁, h₂, h₃ for the columns of H:

Nothing has been assumed yet. Now use the only thing you know for free: r₁ and r₂ are two columns of a rotation matrix, so they are orthogonal and both have unit length. That is two scalar facts, and they are the entire method.

Step 3 — impose orthogonality. r₁^Tr₂ = 0 becomes

The unknown λ multiplies both sides of an equation whose right-hand side is zero, so it cancels. This is the load-bearing move of the whole derivation and it is worth saying out loud: the constraint is homogeneous, which is exactly why you never needed to know the scale of H, and also — the same fact wearing a different hat — why the absolute size of the board never enters the intrinsic solve. Zhang gives you f_x in pixels whether the squares are 25 mm or 25 km. That is the gauge problem of Chapter 0, visible right here in the algebra.

Step 4 — impose equal length, not unit length. The naive move is to write r₁^Tr₁ = 1. You cannot: that equation contains λ, which you do not know. What survives is the difference of the two norm conditions, where λ cancels again:

So one view of the plane yields exactly two linear equations in the entries of B. Not three, not six. If a derivation claims a view gives more, ask which of the extra equations still contains λ.

Step 5 — count. B = K^-TK^-1 is symmetric 3×3, so it has 6 distinct entries: b₁₁, b₁₂, b₂₂, b₁₃, b₂₃, b₃₃. Both constraints are homogeneous in B (one is "= 0", the other is a difference), so scaling B changes nothing and one degree of freedom is unphysical: 5 free parameters. Each view contributes 2 rows to the stacked system Vb = 0. Therefore

and if you additionally fix the skew at zero — true for every digital sensor made this century — then b₁₂ = 0, B has 4 free parameters, and 2n ≥ 4 gives n = 2. Keep both numbers. The three-view answer with the skew caveat is the complete answer; "three" alone is the memorised one.

Step 6 — get K back. B is positive definite, so B^-1 = K K^T and a Cholesky factorisation of B^-1 hands you K directly (up to a sign convention on the diagonal). No iteration, no initial guess. The LM refinement that follows exists only to add distortion and to reweight by the real corner noise.

Step 7 — and then say why the count is necessary but not sufficient. This is the sentence that separates a good understanding from an excellent one. Three views satisfy 2n ≥ 5, but if all three board planes are parallel, every view produces the same two rows of V up to scale: V has rank 2 no matter how many images you take, and twenty views solve nothing. The count is about the number of rows; observability is about their independence. That is why Chapter 1's tilt bench exists, and why σ(f_x) — not the view count and not the RMS — is the thing to gate on.

The second derivation, and it takes ten seconds. σ_td appears twice in the tables below, so do not let it arrive uninvited. A time offset t_d displaces a tracked feature by exactly the distance it travels in that time, so the residual it induces on one feature is r = u̇·t_d, where u̇ is that feature's image velocity in px/s. Invert it: from one feature, t̂_d = r/u̇, so a corner-localisation noise of σ_px becomes a timing noise of

Averaging N features whose corner noise is independent shrinks the standard deviation by √N, which gives the form in the cheat sheet:

and note what the u̇ in the denominator is telling you before you ever plug numbers in: as the image slows down, the uncertainty diverges. A parked robot cannot calibrate its own time offset at any N. Same fact as ∂r/∂t_d = −u̇ being a column of zeros — unobservability and infinite variance are two descriptions of one thing. The numbers drill in section 7 works this to a figure.

1 · The cheat sheet

2 · System-design patterns

Prompt A — "Design the calibration pipeline for a fleet of 500 robots, each with 2 cameras, a LiDAR and an IMU."

Framework, in this order:

1. Split factory from field. Intrinsics and rigid extrinsics are factory jobs with a controlled target and a fixture that guarantees tilt coverage. Temporal offsets and slow extrinsic drift are field jobs, estimated online.
2. Every artifact is versioned and immutable, keyed by sensor serial number, carrying the raw images, the covariance, the operator, and the target's measured square size with its uncertainty. Never a bare YAML.
3. Gate on covariance, not RMS. The factory line rejects a calibration whose σ(f_x) exceeds 2 px, regardless of how pretty the RMS is. State the number.
4. Budgets: 22 images × 1.23 MB = 27 MB captured, 180 ms to solve 141 parameters against 2376 residuals, ~350 B published per camera_info, and a 9.8 MB undistortion LUT costing 1.5 ms per frame at runtime.
5. Monitoring is the deliverable. Four residual regressions at 0.1 Hz (1.6 Mflop) plus one odometry-ratio watchdog at 1/60 Hz. Without the watchdog the fleet is blind to gauge faults.
6. Rollout: a calibration change is a deploy. Canary it, keep the previous artifact, and be able to roll back in one command.

Prompt B — "We are adding a fourth camera. What changes?"

The trap is answering "run the calibration again." The real answer is about topology: with four cameras you have six possible pairwise extrinsics but only three independent ones, so pairwise calibration will produce an inconsistent loop — go around the ring and you will not get back to identity. Calibrate the rig jointly against a shared board set (MC-Calib does exactly this), report the loop-closure residual as a health metric, and store the extrinsics as a tree with one designated root rather than as a set of pairs. Then: does the new camera overlap any existing one? If not, you need a board large enough to be seen by two cameras at once, or a mirror, or a motion-based method.

Prompt C — "The hardware team wants to save $40 per robot by dropping the hardware trigger and timestamping in software."

Do not say no; price it. Software timestamping on a USB3 driver adds 3–12 ms of receive-side latency with several milliseconds of jitter, and jitter is the part that hurts — a constant offset is estimable, a random one is noise you cannot remove. At 1 m/s, 4 ms of jitter is 4 mm of position noise per frame, and at 60°/s it is 0.24° of pointing noise. Then state the condition: "If our localisation budget is 5 cm and our pointing budget is 1°, I can absorb this by estimating t_d online and modelling the jitter as extra measurement noise — at the cost of roughly a factor of two in effective feature noise. If the budget is 1 cm, the trigger is cheaper than the engineering." Then ask what the 500-unit volume makes the real number.

Prompt D — "How do you know your calibration is still good, six months in?"

Three independent instruments, and say why three. (1) The four residual regressions catch every geometric fault and localise it to a parameter. (2) The odometry-versus-VIO distance ratio catches gauge faults, which regression 1 structurally cannot. (3) Localisation against a long-lived map catches slow drift that both of the others average away, because the map is a reference that does not drift with the robot. Then the trends: RMS against temperature (thermal, reversible), t_d against time (clock drift, in ppm), and per-camera RMS spread (one bad sensor hiding in an aggregate).

3 · Coding drills

Drill 1 — "Recover the board pose from a homography." The most common calibration coding question, because it is short and it has a trap.

python
def pose_from_homography(K, H):
    """K : (3,3) float64, PIXEL units   -- fx, fy, cx, cy, skew
       H : (3,3) float64, board -> image, DEFINED ONLY UP TO SCALE
       returns R : (3,3) float64, det(R) = +1
               t : (3,)  float64, METRES, board origin in the camera frame
       t inherits its unit from the board's square size. Feed findHomography
       a board grid in mm and t comes back in mm. This is the gauge."""
    Ki = np.linalg.inv(K)
    h1, h2, h3 = H[:, 0], H[:, 1], H[:, 2]
    lam = 1.0 / np.linalg.norm(Ki @ h1)   # r1 is UNIT length - that pins lambda
    r1, r2 = lam * (Ki @ h1), lam * (Ki @ h2)
    t      = lam * (Ki @ h3)
    r3     = np.cross(r1, r2)                # not free: fixed by r1 and r2
    U, _, Vt = np.linalg.svd(np.column_stack([r1, r2, r3]))
    R = U @ np.diag([1, 1, np.linalg.det(U @ Vt)]) @ Vt   # back onto SO(3)
    return R, t

Say while writing: "H is only defined up to scale, and the thing that fixes the scale is physics, not algebra — r₁ is a column of a rotation matrix, so it has unit length. If I skip λ, R is still roughly right and t is wrong by a pure scale factor, which is the same signature as a wrong square size and just as hard to notice. And the final SVD is not decoration: with real corner noise, r₁ and r₂ are not exactly orthogonal, and every downstream consumer assumes R^TR = I exactly."

Now run one concrete round trip through it, because "t is wrong by a scale factor" is a claim, and a claim you can produce numbers for is an answer. Take the lesson's camera, f_x = f_y = 800, c_x = 640, c_y = 480, and a board rotated 20° about the camera's y axis at t = [0.12, −0.05, 0.83] m. The homography your corner fit actually returns — already normalised by the library to some arbitrary scale — is

So dropping one line costs you ‖t‖ = 0.83702 m instead of 0.84012 m — 3.1 mm short on an 0.84 m baseline, every component low by the same 0.371 %. Note what did not change: R comes out of an SVD, and svd(cM) = U·(cS)·V^T shares its U and V^T with svd(M), so the rotation is bit-identical either way. You get a perfect rotation, a uniformly shrunken translation, and a reprojection error that never moves — because the board corners still land where they landed. That is the exact signature of a wrong square size from Chapter 0, arriving this time from a missing line of code. That is the sentence that shows you know what the trap costs, not merely that there is one.

Drill 2 — "Solve AX = XB for the rotation."

python
def handeye_rotation(A_list, B_list):
    """A_list, B_list : lists of (4,4) float64 SE(3) RELATIVE motions, same
       length, index-aligned. A[i] is sensor A's motion from pose i to i+1
       in A's own frame; B[i] is the SAME physical motion seen by B. Do NOT
       pass absolute poses -- AX = XB is a statement about motions.
       so3_log(T) : (4,4) -> (3,) axis-angle vector, RADIANS, norm = angle.
       returns R_AB : (3,3) float64, det = +1, rotates B-frame into A-frame.
       Two motions minimum, and their axes must NOT be parallel."""
    M = sum(np.outer(so3_log(A), so3_log(B)) for A, B in zip(A_list, B_list))
    U, S, Vt = np.linalg.svd(M)
    if S[1] < 1e-6 * S[0]:
        raise ValueError("rank 1: all rotation axes are parallel")
    return U @ np.diag([1, 1, np.linalg.det(U @ Vt)]) @ Vt

Say while writing: "The determinant guard is the difference between a rotation and a reflection — both minimise the same cost. And I raise on rank deficiency rather than returning a number, because every library I have used will happily hand you an answer from one motion. The requirement is two motions about non-parallel axes; counting motions is the wrong check."

Drill 3 — "Estimate a time offset between two rate signals."

python
def time_offset(a, b, dt, max_lag=20):
    """a, b : 1-D float64, EQUAL length, uniformly sampled at dt seconds.
       Rate signals, not poses -- e.g. gyro |omega| and image-derived
       rotation rate. Different units and gains are fine (ncc is scale-free).
       dt : float, SECONDS per sample. max_lag : int, SAMPLES either way.
       returns float, SECONDS. POSITIVE means a LAGS b: a's event at
       t + offset corresponds to b's event at t, so ADD it to a's stamps."""
    def ncc(u, v):
        u = u - u.mean(); v = v - v.mean()
        return (u @ v) / np.sqrt((u @ u) * (v @ v))   # gain-invariant
    lags = np.arange(-max_lag, max_lag + 1)
    c = np.array([ncc(a[L:], b[:len(b)-L]) if L > 0 else
                  ncc(a[:len(a)+L], b[-L:]) if L < 0 else ncc(a, b) for L in lags])
    k = int(np.argmax(c)); k = min(max(k, 1), len(lags) - 2)
    y0, y1, y2 = c[k-1], c[k], c[k+1]
    d = 0.5 * (y0 - y2) / (y0 - 2*y1 + y2)          # parabola vertex
    return -(lags[k] + d) * dt

Say while writing: "Normalised, because the two sensors have different gains and I do not want a gyro scale error to move my peak. And the parabola, because the integer grid is 5 ms and the offset I am chasing is 18 — without sub-sample refinement my answer is quantised to worse than the thing I am measuring. On clean data the refinement is about a hundred times better than the grid."

Drill 4 — "Given a residual buffer, tell me which calibration is wrong." This is the drill that separates reading about calibration from owning it, so it does not get to be a pointer at another chapter. Write the whole thing. It is one helper, four calls and a ladder, and the ladder is the part that carries the diagnosis.

python
def triage(res, rad, vel, invZ, row):
    """All five are 1-D float64 arrays of the SAME length N (one entry per
       tracked feature in the buffer; N ~ 200 is plenty).
         res  px      signed reprojection residual along the flow direction
         rad  px      radius from the principal point       (0 .. ~500)
         vel  px/s    image-plane feature velocity          (0 .. ~900)
         invZ 1/m     inverse range from the depth sensor   (0.04 .. 0.8)
         row  index   sensor row the feature was read out on (0 .. H-1)
       returns (name: str, value: float, t: float). The VALUE's unit differs
       per channel and is printed below -- that is the point of the ladder."""
    def ols(x, y):                       # slope + its t-statistic
        x = x - x.mean(); y = y - y.mean()
        Sxx = x @ x
        b   = (x @ y) / Sxx
        e   = y - b * x
        se  = np.sqrt((e @ e) / (len(x) - 2) / Sxx)   # SE(b)
        return b, b / se

    b_r, t_r = ols(rad,  res)        # px per px    -> df/f  = b_r
    b_v, t_v = ols(vel,  res)        # seconds      -> t_d   = b_v
    b_z, t_z = ols(invZ, res)        # px*m         -> dt    = b_z / f
    b_w, t_w = ols(row,  res / vel)  # s per row    -> t_line = b_w
    a  = res.mean()                                     # px  the constant channel
    ta = a / (res.std(ddof=1) / np.sqrt(len(res)))   # its own t

    # ORDER MATTERS. Rolling shutter also lights the velocity channel,
    # so it must be tested BEFORE t_d or it will be misread as one.
    if   abs(t_r) > 4: return "intrinsics: focal or k1", b_r, t_r
    elif abs(t_w) > 4: return "rolling shutter",          b_w, t_w
    elif abs(t_v) > 4: return "time offset t_d",         b_v, t_v
    elif abs(t_z) > 4: return "extrinsic translation",   b_z, t_z
    elif abs(ta)  > 4: return "extrinsic rotation (or c_x)", a, ta
    else:               return "no residual signature: check the gauge", 0.0, 0.0

Say while writing: "Five channels, five faults, and the ladder is ordered by confusability, not by how common the fault is. Radius first because an intrinsic error contaminates every other channel downstream and has to be cleared before anything else means what it says. Rolling shutter before t_d, because rolling shutter is a per-row t_d and therefore lights the velocity channel too — dividing the residual by u̇ before regressing against row is exactly what separates them. And I gate on the t-statistic, not the slope, because a slope of 0.4 ms is decisive with 200 features and meaningless with six. The last branch is the one that matters: it does not say 'clean', it says 'stop looking at residuals'. Five flat channels is a positive result — it has ruled out every fault that has a residual signature, which leaves the class that does not."

Two follow-ups they will ask, so have both. "Why 4 and not 2?" Because you run this on five channels at 0.1 Hz on a fleet: at |t| > 2 you are firing on roughly one channel in twenty every ten seconds and the on-call engineer stops reading the alerts within a day. |t| > 4 is about 6×10^-5 per channel per test, which on 500 robots × 5 channels × 8640 tests a day is still a handful of false pages — so in production the alert also requires the channel to stay hot for three consecutive windows. "What if two faults are present?" The ladder returns the first one, which is correct behaviour: fix it, re-run, and the second one is now the top of the ladder. Reporting both at once from one pass is how you get a confident wrong answer, because a large radius slope biases every other channel's estimate.

4 · Debugging scenarios

5 · Classical versus modern, and when to use which

6 · Recommended reading

The one book. Hartley & Zisserman, Multiple View Geometry in Computer Vision (2nd edition, 2004). Chapter 8 on the absolute conic is where Zhang's B matrix actually comes from, and chapter 10's treatment of projective, affine and metric reconstruction is the rigorous version of this lesson's "the gauge" argument. If you read one chapter, read chapter 8. The runner-up, for the estimator underneath everything here, is Timothy Barfoot, State Estimation for Robotics (2nd edition, 2024).

Five papers, and why each.

1. Zhang, "A Flexible New Technique for Camera Calibration," TPAMI 2000. Because it is still the default in every toolbox and you will be asked to derive it. Read it for the two constraints, not the closed form.
2. Park & Martin, "Robot sensor calibration: solving AX = XB on the Euclidean group," IEEE T-RA 1994. Because it makes the rotation/translation decoupling explicit and treats the observability properly. Tsai & Lenz (1989) is the more-cited ancestor; Daniilidis (IJRR 1999) is the dual-quaternion joint solve worth naming as an alternative.
3. Furgale, Rehder & Siegwart, "Unified Temporal and Spatial Calibration for Multi-Sensor Systems," IROS 2013. The Kalibr paper. Read it for the continuous-time B-spline trajectory, which is the right way to handle sensors at different rates, and for its honesty about what is observable.
4. Qin & Shen, "Online Temporal Calibration for Monocular Visual-Inertial Systems," IROS 2018. Six pages, one idea, and it is the idea from Chapter 3: shift the feature by t_d times its image velocity and let the optimiser have the derivative. Pair with Li & Mourikis (IJRR 2014) for the filter version.
5. Lv et al., "Observability-Aware Intrinsic and Extrinsic Calibration of LiDAR-IMU Systems," T-RO 2022. The modern answer to "how much data do I need": not more, but better conditioned. It formalises the rank arguments this lesson makes by hand.

Five repositories, and what to look at in each.

1. ethz-asl/kalibr — the camera-IMU spatial and temporal calibrator. Look at the report it prints: it gives you parameter covariances and a warning when a direction is poorly excited. Reading that report is a skill in itself.
2. opencv/opencv, module calib3d — read calibrateCameraExtended for how the LM parameterisation is set up and which flags fix which parameters, and calibrateHandEye for all five hand-eye methods side by side in one file. It is the fastest way to see how they differ.
3. rpng/open_vins — look at how the state class carries optional camera intrinsics, camera-IMU extrinsics and calib_camimu_dt, and how each gets its Jacobian block. This is the reference implementation of Chapter 3's one extra column.
4. hku-mars/livox_camera_calib — targetless LiDAR-camera calibration by aligning depth discontinuities to image edges. Look at the edge extraction and the cost function; it is the clearest small example of an appearance-based extrinsic solve.
5. rameau-fr/MC-Calib — multi-camera rig calibration against multiple ChArUco boards, including non-overlapping cameras. Look at how it builds and closes the camera graph; that is the answer to Prompt B above.

7 · The numbers drill

One significant figure, out loud, no calculator. If you cannot do these in about three seconds each, you cannot hold a design conversation at pace.

A one-line "how" column is a reminder for someone who has already done the arithmetic once. If you have not, the reminder is useless, so here are three of those rows worked with every intermediate on the page. Do them by hand before you trust the one-liners.

Worked A — σ on t_d from 200 features. Inputs: N = 200 tracked features, image velocity u̇ = 240 px/s, corner-localisation noise σ_px = 0.30 px. Formula from section 0: σ_td = σ_px / (√N · u̇).

Read step 1 against step 5 before you move on: 200 features bought you a factor of 14, from 1.25 ms to 88 µs. That is the whole economics of the estimate. It also tells you the shape of every follow-up — wanting 10 µs from the same motion needs N ≈ 15,000, which you will not have, so you get it by increasing u̇ instead. Doubling the yaw rate is worth as much as quadrupling the feature count, and it is free.

Worked B — the lever arm you cannot see. Inputs: the hand-eye translation solve has 1 cm of motion error, δ = 0.010 m, and the calibration dance only rotated the rig by θ = 2°. The conditioning factor is ‖(R − I)t‖ = 2 sin(θ/2)·‖t_⊥‖, so the error in the recovered lever arm is δt_X = δ / (2 sin(θ/2)).

The ratio is the punchline: 0.28650 / 0.0070711 = 40.5× — call it a factor of forty. The same 1 cm of motion error produces 7 mm of lever-arm uncertainty after a proper 90° excitation and 29 cm after a timid 2° wrist-flick. Nothing about the data volume changed; only the conditioning did. And notice that at small θ the factor is 2 sin(θ/2) ≈ θ, so the error grows as 1/θ and diverges at θ = 0 — which is the singularity of (R_A − I) from Chapter 2, arriving as a number instead of a rank argument. This is why "we collected 400 images" is not an answer to "is your extrinsic observable".

Worked C — the undistortion LUT, and why it costs 1.5 ms. Inputs: a 1280×960 sensor, undistortion precomputed by initUndistortRectifyMap into two float32 maps (one for the source x, one for the source y), consumed by remap at 30 Hz.

Step 5 is the claim in the system-design section, and now it has a derivation instead of a vibe: 1.5 ms per frame is spent before a single source pixel has been touched, because the map is bigger than any last-level cache and has to come from DRAM every frame. The source-image gather that follows is worse per byte, because it is random access rather than streaming. Two consequences worth having ready: the fixed-point form (CV_16SC2 coordinates plus a CV_16UC1 interpolation-weight table) is 4 + 2 = 6 B/px = 7.37 MB, about 25 % less traffic for a sub-pixel accuracy loss you will not measure; and on four cameras this is 1.18 GB/s of pure map traffic, at which point undistorting on the GPU — or not undistorting at all and pushing the distortion model into the feature tracker — stops being a micro-optimisation and becomes an architecture decision.

8 · The five sentences to have ready

1. The gauge sentence. "The board square size is the only quantity in the dataset carrying units of metres, so it is the gauge — it scales every extrinsic translation, leaves the intrinsics and the residual bit-identical, and can only be caught by an external metric reference."
2. The tilt sentence. "Tilt is not folklore. It is what makes the two homography columns explore different directions, which is what makes the constraint rows independent. Under about five degrees, σ(f_x) goes from a few pixels to a few thousand while the RMS never moves."
3. The lever-arm sentence. "You cannot see a lever arm without rotating about it. Two points on a rigid body translating together move identically, so (R_A − I) is singular and the translation is in its null space."
4. The time-offset sentence. "It is one extra state and one extra Jacobian column, and the column is minus the feature's image velocity. Which also tells me exactly when it is unobservable: whenever nothing is moving in the image."
5. The triage sentence. "Regress the residual against radius, image velocity, inverse range and row index, and read the intercept too. Five channels, five faults. If all five are clean and the system is still wrong, it is a gauge, and I go find a second measurement of a distance."

9 · The questions that X-ray a stack

Calibration questions are unusually revealing to ask of any robotics organisation — a team you are joining, a vendor you are evaluating, or your own — because the answers tell you a great deal about how it actually operates. Four worth having ready:

Asking these also does something for you that no dashboard can: each one points, in the form of a question, at a place where this subject's failures live. When you can predict what a good answer sounds like before you hear it, the lesson has done its job.

And when you want to pressure-test yourself under a clock, the Studio button on this page runs a timed practice session on exactly this material.