Zhang's method from first principles

The previous post ended with a few hundred sub-pixel corner positions and a shrug: they are only coordinates. This one turns them into the camera.

The claim is genuinely strange the first time you meet it. Print a checkerboard. Hold it up to your camera at half a dozen angles you do not measure, at distances you do not measure, in orientations you never write down. Take a photo each time. From nothing but the corner positions in those photos (no ruler, no rig, no knowledge whatsoever of where the board was), you can recover the camera’s focal lengths in pixels, where its optical axis pierces the sensor, and how far its pixel grid is from square.

That is Zhang’s method, and it is why every calibration tutorial on the internet involves a sheet of A4 taped to a cereal box. What follows is the chain that makes it work, written from the ~140 lines of NumPy I wrote for it in 2019, plus the parts I got wrong.

From a lens to three numbers

A real lens obeys the thin-lens equation, and the report opens with it (Eq. 2.1):

$$\frac{1}{f} = \frac{1}{d_0} + \frac{1}{d_1}$$

which is only a reasonable approximation when the lens is much thinner than its radius of curvature. Even so it is too complicated to calibrate, because a single point in the world maps to a region of the sensor whenever it is out of focus: light from one source spreads across several photosites, and once information has overlapped like that there is nothing to undo. So you assume the interesting parts of the scene are in focus and swap the lens for an ideal pinhole: an aperture whose size tends to zero, blocking every ray from a source except the one that goes straight through.

Projection through a pinhole is then two operations. First the perspective divide: push every 3-D point onto the plane one unit in front of the camera by dividing through by its own depth (Eq. 2.2):

$$\mathbf{X}_i = \begin{pmatrix} x_i \\ y_i \\ z_i \end{pmatrix} \;\rightarrow\; \mathbf{X}'_i = \begin{pmatrix} x_i / z_i \\ y_i / z_i \\ 1 \end{pmatrix}$$

Then scale that normalised plane up to pixels. If a lens were perfect and a sensor were perfectly made and perfectly placed, one number would do it (Eq. 2.3):

$$\mathbf{U}_i = \mathbf{K}\cdot\mathbf{X}'_i = \begin{pmatrix} f' & 0 & 0 \\ 0 & f' & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} x'_i \\ y'_i \\ 1 \end{pmatrix}$$

It will not do. The lens does not scale equally in $x$ and $y$, and the sensor applies its own scaling in $u$ and $v$, so $f'$ splits into $f_u$ and $f_v$. The sensor is not perfectly centred on the optical axis either (you would like the shift to be exactly half the resolution each way, and it never is), so two more parameters $c_u, c_v$ move the origin. And if the sensor is not quite square to the imaging plane you get a shear, absorbed into a single parameter $s$. That is the whole camera calibration matrix (Eq. 2.4):

$$\mathbf{K} = \begin{pmatrix} f_u & s & c_u \\ 0 & f_v & c_v \\ 0 & 0 & 1 \end{pmatrix}$$

Five numbers. Every one of them is a property of the camera and none of them is a property of the scene, which is exactly why they are worth recovering once and reusing forever.

Where the camera is takes six more numbers (three of rotation and three of translation) combined into one transformation (Eq. 2.5):

$$\mathbf{T} = \begin{pmatrix} R_{11} & R_{12} & R_{13} & T_1 \\ R_{21} & R_{22} & R_{23} & T_2 \\ R_{31} & R_{32} & R_{33} & T_3 \end{pmatrix}$$

Twelve entries, six degrees of freedom: the nine rotation entries are three angles wearing a disguise, because a rotation matrix is orthonormal. That constraint looks like bookkeeping here and turns out to be the entire trick later.

Why a plane changes everything

The straightforward way to calibrate is the Direct Linear Transform: photograph an object whose 3-D geometry you know exactly, write down the correspondences, and solve directly for the projection. Four correspondences are enough in principle; in practice you use many more, the system is over-determined, and you minimise the residual with an SVD. It works well (the report calls it highly accurate) but it comes with two conditions that are painful in real life. The reference points must not all lie on one plane, and no three of the four may be collinear. Building an accurate non-planar 3-D calibration object is a machining problem, not a printing problem.

Zhang’s method makes the opposite bargain. It insists that all the reference points lie on one plane, and gains a great deal from it. Put the world coordinate frame on the board itself and every board point has $Z = 0$, which kills the third column of the rotation matrix (it is multiplied by zero) and collapses the $3\times 4$ projection into a $3\times 3$ homography (Eq. 3.10):

$$s\begin{pmatrix} u_{m,n} \\ v_{m,n} \\ 1 \end{pmatrix} = \mathbf{K}\begin{pmatrix} R_{11} & R_{12} & T_1 \\ R_{21} & R_{22} & T_2 \\ R_{31} & R_{32} & T_3 \end{pmatrix} \begin{pmatrix} X_n \\ Y_n \\ 1 \end{pmatrix} \qquad\Longrightarrow\qquad s\,\mathbf{x}_{n,m} = \mathbf{H}_m \mathbf{X}_n$$

The scale $s$ is free: a homography is only defined up to scale, which sounds like a loss and is actually a convenience: you solve for $\mathbf{H}$ up to scale and fix the scale afterwards. The price is stated plainly in the report: the extrinsic parameters cannot be fully recovered, and the method is less accurate than a true DLT against a machined target. You need at least four points from at least two views; in practice, as we will see, three views is the real floor.

Step one: one homography per view

Rearranging $s\,\mathbf{x} = \mathbf{H}\mathbf{X}$ to eliminate $s$ gives two linear equations per correspondence in the nine unknowns of $\mathbf{H}$, and $N$ points stack into a $2N \times 9$ system (Eq. 3.11), where $\delta_m$ is whatever residue the noise in the corner detection leaves behind:

$$\begin{pmatrix} -X_n & -Y_n & -1 & 0 & 0 & 0 & x_{n,m}X_n & x_{n,m}Y_n & x_{n,m} \\ 0 & 0 & 0 & -X_n & -Y_n & -1 & y_{n,m}X_n & y_{n,m}Y_n & y_{n,m} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \end{pmatrix} \begin{pmatrix} H_{11} \\ H_{12} \\ H_{13} \\ H_{21} \\ H_{22} \\ H_{23} \\ H_{31} \\ H_{32} \\ H_{33} \end{pmatrix} = \delta_m \qquad\text{i.e.}\qquad \mathbf{A}_m \bar{H}_m = \delta_m$$

Minimising $\|\delta_m\|$ subject to $\|\bar H_m\| = 1$ is the classic homogeneous least-squares problem: decompose $\mathbf{A}_m = \mathbf{U}_m\mathbf{S}_m\mathbf{V}_m^\top$ and take the right singular vector belonging to the smallest singular value. NumPy does that for every view at once, because np.linalg.svd broadcasts over leading axes. From Modules/ZhangAlgV2.py, L56–68:

M = np.empty((Nx.shape[0], 2 * Nx.shape[1], 9), dtype=np.float64)
M[:, 0::2] = np.stack((-NX[:, :, 0], -NX[:, :, 1], nones,
                       zeros, zeros, zeros,
                       NX[:, :, 0] * Nx[:, :, 0], NX[:, :, 1] * Nx[:, :, 0], Nx[:, :, 0]),
                      axis=2)
M[:, 1::2] = np.stack((zeros, zeros, zeros,
                       -NX[:, :, 0], -NX[:, :, 1], nones,
                       NX[:, :, 0] * Nx[:, :, 1], NX[:, :, 1] * Nx[:, :, 1], Nx[:, :, 1]),
                      axis=2)
_, s_1, Vt = np.linalg.svd(M, full_matrices=False)
H = Vt[:, -1].reshape(-1, 3, 3)
Hs = np.dot(np.einsum("ijk,ikm->ijm", iMSx, H), MSX[0])
Hs /= Hs[:, 2, 2][:, None, None]

Thirteen lines for every homography in the capture. The [:, 0::2] / [:, 1::2] slicing is how the two rows per point get interleaved without a loop, and the einsum on the second-to-last line is the un-normalising step, which needs its own section.

Step two: normalisation, and what it is actually for

The rows of that matrix are a mess of scales. Board coordinates are tens of millimetres; image coordinates are hundreds of pixels; the products $x_{n,m}X_n$ are hundreds of thousands. Columns whose entries differ by five orders of magnitude make for an ill-conditioned matrix, and the standard fix (Hartley’s) is to pre-transform both point sets so they are centred on their own mean with an RMS radius of $\sqrt 2$, solve, and then undo the transform on the answer (Eq. 3.12):

$$\mathbf{H}_m = N_{x,m}^{-1}\,\mathbf{H}'_m\,N_X$$

get_scale_matrix() builds $N$ and $N^{-1}$ in the same call (L6–20):

def get_scale_matrix(pts, inverted_matrix=True):
    avg = np.mean(pts, axis=1)
    std = 1.4142135623730951 / np.std(pts, axis=1)
    mtx = np.zeros((std.shape[0], 3, 3))
    mtx[:, 2, 2] = 1
    if inverted_matrix:
        mtx[:, 0, 0] = 1 / std[:, 0]
        mtx[:, 1, 1] = 1 / std[:, 1]
        mtx[:, 0, 2] = avg[:, 0]
        mtx[:, 1, 2] = avg[:, 1]
    else:
        mtx[:, 0, 0] = std[:, 0]
        mtx[:, 1, 1] = std[:, 1]
        mtx[:, 0, 2] = -std[:, 0] * avg[:, 0]
        mtx[:, 1, 2] = -std[:, 1] * avg[:, 1]

1.4142135623730951 is $\sqrt2$, typed out. One deviation worth flagging: np.std(pts, axis=1) is taken per axis, so the scaling is anisotropic ($\sqrt2/\sigma_x$ across and $\sqrt2/\sigma_y$ down) where Hartley’s prescription uses one isotropic radius. It conditions the matrix just as well; it is simply not quite the textbook transform, and my Rust port copies the Python rather than the textbook so the two agree.

Now the honest part, because I set out to write “watch the recovered K explode when you turn normalisation off” and it does not explode.

Measured, on the real dataset 3 corners with a 25 mm board: the worst per-view condition number is 5.5 normalised and 2.1 × 10⁵ raw, four and a half orders of magnitude, exactly as advertised. And the recovered K moves from $f_u = 4343$ to $f_u = 4340$. Three parts in four thousand. On dataset 2 the condition number goes from 5.3 to 4.4 × 10⁵ and $f_u$ moves by one pixel.

Two things are covering for the bad conditioning, and both are worth knowing:

• Double precision has about sixteen digits and the problem only eats five of them. A condition number of $10^5$ is not large when your arithmetic has $10^{16}$ of headroom.
• One-sided Jacobi never forms $\mathbf{A}^\top\mathbf{A}$. It rotates pairs of columns until they are mutually orthogonal and reads the singular values off the column norms, which buys it high relative accuracy on badly scaled matrices. That is the algorithm in my Rust. NumPy goes a different way: LAPACK’s gesdd bidiagonalises and then divides and conquers, but it is backward stable, so at a condition number of $10^5$ it too loses only about five of its sixteen digits. Either way, the arithmetic absorbs the damage.

Take either away and it falls over immediately. The textbook shortcut for a homogeneous system is to form the $9\times 9$ Gram matrix $\mathbf{A}^\top\mathbf{A}$ and take the eigenvector of its smallest eigenvalue, cheap, and it squares the condition number in the working precision, turning $2\times10^5$ into $4\times10^{10}$. In single precision that is past the end of the number system. Here is dataset 3 solved four ways, with the LM refinement off so you are looking at the algebraic solution alone:

Null space found by	Normalised	$f_u$	$s$	$c_u$	$f_v$	$c_v$
Jacobi SVD, f64	yes	4343	−24.1	406	4303	154
Jacobi SVD, f64	no	4340	−17.8	393	4303	174
$\mathbf{A}^\top\mathbf{A}$, f32	yes	4343	−24.1	406	4303	154
$\mathbf{A}^\top\mathbf{A}$, f32	no	5085	960	4523	2456	18646

A principal point 18 646 pixels down a 576-pixel image. On the synthetic camera in the widget the same combination does not even get that far: the re-projection error climbs to 11 pixels and $\mathbf{B}$ comes back indefinite, so there is no K at all. Switch normalisation back on and that same naive solver lands within three pixels of the truth.

The honest summary, then, is that normalisation is insurance rather than a fix. With a careful solver and double precision you will not miss it; it costs four lines, and it makes the problem well conditioned enough that a naive solver, a single-precision pipeline or a future you who reaches for the Gram matrix cannot get hurt. Take it.

The earlier version of the file, Modules/ZhangAlg.py, is the foil. It does not derive the scaling from the data at all: it hard-codes the sensor size (L13–17):

ScaleM = np.array([[2/640.0, 0, -1],
                   [0, 2/480.0, -1],
                   [0,       0,  1]])
iScaleM = np.linalg.inv(ScaleM)
all_detected_points_scaled = 2*all_detected_points / [[480.0], [640.0]] -1

which maps a VGA image to $[-1, 1]^2$ and does nothing at all to the world coordinates. It is better than nothing and it stops being correct the moment you point it at a 768×576 frame. That constant is why V2 exists.

Step three: Levenberg–Marquardt, which Zhang does not do

The SVD minimises $\|\mathbf{A}\bar H\|$, which is an algebraic residual: a quantity with no units and no geometric meaning. What you actually care about is the re-projection error in pixels (Eq. 3.13):

$$E_m = \sum_{n=0}^{N-1}\left\|\mathbf{x}_{n,m} - \mathbf{H}_m\mathbf{X}_n\right\|^2$$

The report is explicit that this step “isn’t done in the original paper by Zhang” and comes from Burger’s tutorial instead. It is non-linear because of the perspective divide, so it wants Levenberg–Marquardt, and LM wants a Jacobian. Writing one out by hand is a small pleasure. With $w = h_{31}X + h_{32}Y + h_{33}$ and $u = (h_{11}X + h_{12}Y + h_{13})/w$:

$$\frac{\partial u}{\partial \bar h} = \left(\tfrac{X}{w},\; \tfrac{Y}{w},\; \tfrac{1}{w},\; 0,\; 0,\; 0,\; -\tfrac{uX}{w},\; -\tfrac{uY}{w},\; -\tfrac{u}{w}\right)$$ $$\frac{\partial v}{\partial \bar h} = \left(0,\; 0,\; 0,\; \tfrac{X}{w},\; \tfrac{Y}{w},\; \tfrac{1}{w},\; -\tfrac{vX}{w},\; -\tfrac{vY}{w},\; -\tfrac{v}{w}\right)$$

The first six entries are the quotient rule’s easy half. The last three are where the perspective divide bites: every parameter in the bottom row of $\mathbf{H}$ moves the denominator, so it drags the projected point by $-1/w^2$ times the numerator, which is $-u/w$, the point’s own position scaled down by its depth. That is the whole content of the $-1/w^2$ terms, and it is why a homography Jacobian looks different from a linear one.

Vectorised over all $N$ points, that is getJ (L33–44):

def getJ(H_vals, X, _):
    H = H_vals.reshape(3, 3)
    proj = np.dot(H, X)
    sw = 1 / proj[2]
    sw2 = -1 / (proj[2] * proj[2])
    swx = proj[0] * sw2
    swy = proj[1] * sw2
    zeros = np.zeros(X.shape[1])
    J = np.empty((9, 2 * X.shape[1]))
    J[:, 0::2] = np.stack((X[0] * sw, X[1] * sw, sw, zeros, zeros, zeros, X[0] * swx, X[1] * swx, swx), axis=0)
    J[:, 1::2] = np.stack((zeros, zeros, zeros, X[0] * sw, X[1] * sw, sw, X[0] * swy, X[1] * swy, swy), axis=0)
    return J.T.reshape(2 * X.shape[1], -1)

swx = proj[0] * sw2 is $-u/w$ and X[0] * swx is $-uX/w$: the derivation above, exactly, in five lines of arithmetic and two np.stacks. It is handed to SciPy with the residual function and MINPACK does the rest (L27–30 and L75–84):

def getError(H_vals, X, x):
    H = H_vals.reshape(3, 3)
    hold = np.dot(H, X)
    return np.square(x - (hold[:2] / hold[2]).T).flatten()

...

for i in range(len(Hs)):
    ret = optAlg(
        getError,
        Hs[i].flatten(),
        jac=getJ,
        args=(h_X, x[i]),
        method='lm'
    )
    if ret.success:
        Hs[i] = ret.x.reshape(3, 3) / ret.x[8]

Step four: the absolute conic, or how a rotation matrix pays you back

Here is the step that makes Zhang’s method feel like a card trick.

Write $\mathbf{H} = \mathbf{K}[\,\mathbf{r}_1\ \mathbf{r}_2\ \mathbf{t}\,]$ with $\mathbf{h}_i$ the $i$-th column of $\mathbf{H}$, so $\mathbf{r}_i \propto \mathbf{K}^{-1}\mathbf{h}_i$. Now use the constraint from the very first section: $\mathbf{R}$ is a rotation, so its columns are orthonormal. Two facts follow, for free, for every single view, whatever the pose was:

$$\mathbf{r}_1^\top\mathbf{r}_2 = 0 \qquad\text{and}\qquad \|\mathbf{r}_1\| = \|\mathbf{r}_2\|$$

Substituting $\mathbf{r}_i \propto \mathbf{K}^{-1}\mathbf{h}_i$ turns both into statements about $\mathbf{h}_i^\top \mathbf{K}^{-\top}\mathbf{K}^{-1}\mathbf{h}_j$, so define (Eq. 3.15):

$$\mathbf{B} = \mathbf{K}^{-\top}\mathbf{K}^{-1}$$

$\mathbf{B}$ is symmetric, so it has six unknowns, and the two constraints are linear in them. Collecting the coefficients gives the vector $v_{p,q,m}$ (Eq. 3.14):

$$v_{p,q,m} = \begin{pmatrix} H_{0,p,m}H_{0,q,m} \\ H_{0,p,m}H_{1,q,m} + H_{1,p,m}H_{0,q,m} \\ H_{1,p,m}H_{1,q,m} \\ H_{2,p,m}H_{0,q,m} + H_{0,p,m}H_{2,q,m} \\ H_{2,p,m}H_{1,q,m} + H_{1,p,m}H_{2,q,m} \\ H_{2,p,m}H_{2,q,m} \end{pmatrix}^{\!\top}$$

and stacking two rows per view gives a $2M\times 6$ homogeneous system for $\bar B = (B_{11}, B_{12}, B_{22}, B_{13}, B_{23}, B_{33})$ (Eq. 3.16):

$$\begin{pmatrix} v_{0,1,0} \\ v_{0,0,0} - v_{1,1,0} \\ v_{0,1,1} \\ v_{0,0,1} - v_{1,1,1} \\ \vdots \end{pmatrix}\cdot \bar B = \bar 0$$

Six unknowns, defined up to scale, so five degrees of freedom; two equations per view. Three views. That is the whole answer to “how many photos do I need”, and it falls straight out of the shapes. More views are noise averaging, not new geometry.

$\mathbf{B}$ is over-determined once you have four or more views, so it goes through the same SVD null-vector machinery as the homographies did. The code builds the three $v$ vectors it needs by hand (L86–92, and two near-identical blocks after it):

v_m_1_2 = np.array((
    Hs[:, 0, 0] * Hs[:, 0, 1],
    Hs[:, 0, 0] * Hs[:, 1, 1] + Hs[:, 1, 0] * Hs[:, 0, 1],
    Hs[:, 0, 0] * Hs[:, 2, 1] + Hs[:, 2, 0] * Hs[:, 0, 1],
    Hs[:, 1, 0] * Hs[:, 1, 1],
    Hs[:, 2, 0] * Hs[:, 1, 1] + Hs[:, 1, 0] * Hs[:, 2, 1],
    Hs[:, 2, 0] * Hs[:, 2, 1]))

Compare that to Eq. 3.14 and the third and fourth entries have swapped places. It is not a bug: the code has simply chosen $\bar B = (B_{11}, B_{12}, B_{13}, B_{22}, B_{23}, B_{33})$ where the report and Zhang’s paper use $(B_{11}, B_{12}, B_{22}, B_{13}, B_{23}, B_{33})$, and the matrix it assembles from $b$ at the end of the function uses the same ordering. Both are self-consistent. Mixing them silently gives you a beautiful, plausible, wrong K, so it is worth checking before you assume you have found a mistake in someone’s code, as I did, in my own, twice.

Step five: Cholesky, and one line that hides a failure mode

$\mathbf{B} = \mathbf{K}^{-\top}\mathbf{K}^{-1}$ has the shape of a Cholesky factorisation already. Cholesky gives $\mathbf{B} = \mathbf{L}\mathbf{L}^\top$ with $\mathbf{L}$ lower triangular; $\mathbf{K}^{-1}$ is upper triangular because $\mathbf{K}$ is; and the factorisation is unique. So $\mathbf{L} = \mathbf{K}^{-\top}$, and

$$\mathbf{K} = \left(\mathbf{L}^\top\right)^{-1}$$

which is where the function ends (L124–132 and L136, the three commented-out lines in between hold the closed-form $\alpha, \beta, \gamma, u_c, v_c$ expressions, the alternative to a Cholesky, still sitting in the source where I abandoned them):

B = np.array([[b[0], b[1], b[2]],
              [b[1], b[3], b[4]],
              [b[2], b[4], b[5]]])
B = B/B[2,2]
try:
    KnT = np.linalg.cholesky(B)
except LinAlgError:
    return "Did not Converge"
K = np.linalg.inv(KnT.T)
return K/K[2,2]

That return "Did not Converge" is my favourite line in the repository, and not because it is good. A function whose contract is “returns a 3×3 array” returning a string on failure is exactly the kind of thing that gets TypeErrord in a batch job at 2 a.m. But it is pointing at something real, and it is worth understanding what.

A null vector from an SVD has an arbitrary sign. $\mathbf{V}$’s columns are only determined up to $\pm 1$, so $b$ and $-b$ are equally valid answers and which one you get depends on the arithmetic. If you get $-b$, then $\mathbf{B}$ is negative definite, Cholesky finds a negative pivot on the very first row, and NumPy raises LinAlgError.

Except that B = B/B[2,2] on the line above quietly fixes it. Multiply $\mathbf{K}^{-\top}\mathbf{K}^{-1}$ out and $B_{33}$ comes to $\left(\tfrac{s c_v - c_u f_v}{f_u f_v}\right)^2 + \left(\tfrac{c_v}{f_v}\right)^2 + 1$, a sum of squares plus one, so strictly positive for any camera whatsoever. Dividing by it restores the sign whichever way the SVD went. V2 is immune. V1 is not: ZhangAlg.py L68–72 calls np.linalg.cholesky(B) on the raw B with no scaling and no try, so a coin-flip in LAPACK’s sign convention crashes it outright. That is the sign ambiguity in the wild, and my Rust exposes both paths so the tests can pin the behaviour down:

// V1: straight to Cholesky. Fails.
assert_eq!(k_from_b(&bmat, false).unwrap_err(), Status::NotPositiveDefinite);
// V2: divide by B[2][2] first. Same K as the unflipped run.
let (k, _, _) = k_from_b(&bmat, true).expect("sign fixed");

What is left of "Did not Converge" in V2, then, is the genuinely interesting failure: $\mathbf{B}$ that is not positive definite on its own merits, because the views did not constrain it. Hold the board parallel to the sensor in every shot and each view repeats the constraints of the last one; the $2M\times6$ system has an effectively two-dimensional null space; the vector you pull out of it is an arbitrary mixture and corresponds to no real camera at all. You can do that on purpose in the widget below: pick “all fronto-parallel”, set the noise slider to zero, and there it is.

Calibrate a camera

Two modes. Synthetic builds a camera you control, projects a board through it at poses you choose, adds pixel noise, and hands the corners (and nothing else) to the same zhang.rs the tests exercise. The K on the right is what the algorithm recovers from those dots alone; the error column is against a truth it never saw. Real corners loads the corner sets my detector actually extracted from the calibration datasets, so you can reproduce a row of the results table below in your browser.

With JavaScript on, this becomes a live calibrator: sliders for a ground-truth camera’s fu, fv, cu, cv and skew, a pixel-noise control, a choice of board poses, and the recovered K beside them with the per-parameter error. Switches turn Hartley normalisation and the Levenberg–Marquardt refinement on and off and change how the null space is computed. It also reads the ordered-corner JSON exported by the detector post.

Things worth doing to it:

• Select “AᵀA in f32: the naive way” and untick both normalisation and LM. The re-projection error jumps from 0.40 px to 10.9 px and the solver gives up: B is not positive definite. Tick normalisation back on, leaving the naive solver selected, and it recovers K to within three pixels. That is the entire argument for the four lines of get_scale_matrix, and it is why the LM box has to come off too. With the refinement running, LM starts from the wrecked homography and walks back to the right answer anyway, which is a nice illustration of what a re-projection-error minimiser is worth.
• Choose “all fronto-parallel” with the noise slider at zero. Every board parallel to the sensor, so every view repeats the constraints of the last one, the conic system is genuinely rank-deficient, and you get the real thing: B is not positive definite, the literal "Did not Converge" of the Python. Now push the noise up to 0.3 px and it is worse, not better: the noise breaks the tie, a K appears, and it claims $f_u = 44670$ for a 640-pixel sensor. The tell is the conic system’s condition number, which goes from about 3 × 10³ on a good capture to 8 × 10⁵ here; the widget says so above the table. Then try “barely tilted”, which is not degenerate and still moves $c_u$ by five pixels. Tilting the board is not a nicety. It is the measurement.
• Drop to two views. The solver refuses rather than returning a number, which is the behaviour I wish more numerical code had.
• Turn LM off with noise at 1 px, everything else sane. K moves by a pixel or two and the RMS barely twitches. The refinement is worth having and it is not the difference between working and not: Zhang shipped without it.
• Load dataset 3 and pull “views used” down from 15. Watch $c_v$ (always the least determined parameter) swing around until there are enough views to pin it.
• Load your own corners. The detector post’s Download corners (JSON) button produces a checkerboard-corners/v1 file, and this widget reads it. One file is one view, so you will want several; the shipped bundles show the {"board": …, "views": […]} shape that holds a whole capture.

Results

The benchmark harness runs four combinations (OpenCV’s corner detector or mine, OpenCV’s calibrateCamera or mine) over each dataset, with every OpenCV distortion coefficient fixed at zero so the two are solving the same problem. This is Table 2 of the report, in full. The bottom row of every matrix is $(0, 0, 1)$ and is omitted:

#	OpenCV corners, OpenCV calibration	OpenCV corners, this calibration	Detected corners, OpenCV calibration	Detected corners, this calibration
1	2740, 0, 567 / 0, 2734, 371	2746, 4, 565 / 0, 2738, 370	2743, 0, 567 / 0, 2737, 370	2746, 4, 565 / 0, 2739, 370
2	692, 0, 318 / 0, 688, 268	684, −2.8, 321 / 0, 678, 279	671, 0, 319 / 0, 681, 246	685, −2.8, 321 / 0, 678, 279
3	4343, 0, 394 / 0, 4311, 157	4361, −10, 407 / 0, 4328, 117	4389, 0, 387 / 0, 4353, 164	4361, −10, 407 / 0, 4328, 117
5	738, 0, 325 / 0, 741, 225	741, 0, 323 / 0, 744, 229	736, 0, 324 / 0, 740, 225	741, 0, 323 / 0, 744, 229
6	755, 0, 392 / 0, 762, 235	753, 3, 391 / 0, 759, 229	NA	NA
7	1312, 0, 791 / 0, 1230, 514	1336, −5, 792 / 0, 1249, 524	1319, 0, 793 / 0, 1235, 517	1336, −5, 792 / 0, 1249, 524

Dataset 4 is absent, and dataset 6 has no “detected corners” columns because the detector scored 0 out of 10 on it: the previous post has that story.

Focal lengths agree to within about half a percent across the board: 2740 against 2746, 4343 against 4361, 1312 against 1336. The report’s own summary of the discrepancy is the one I would still give:

That is the right explanation and it is worth spelling out. OpenCV’s model has no $s$; it pins the skew to zero and lets its tangential distortion terms absorb the same physical effect. This implementation has an $s$ and no distortion at all. So when this code reports $s = -10$ for dataset 3, that number has to come out of somewhere, and it comes out of the other four, which is exactly where you see dataset 3’s 13-pixel disagreement in $c_u$ and its 40-pixel one in $c_v$. Two models, both under-specified in different directions, meeting a few pixels apart. Neither column is “the truth”.

What the Rust gives

The widget is not running the NumPy. It is running my port, on corners produced by my port of the detector, so it is a second implementation of both halves, and that makes it a check. Detection first: 15/15 frames on dataset 3 and 10/25 on dataset 2, which is Table 1 exactly. Then, feeding those corners through zhang.rs with normalisation and LM on:

Dataset	Rust port, its own corners	Report, detected corners + OpenCV calibration	Report, this calibration	RMS
3	4337, −17.4, 400 / 4300, 178	4389, 0, 387 / 4353, 164	4361, −10, 407 / 4328, 117	0.569 px
2	672, −1.2, 321 / 681, 260	671, 0, 319 / 681, 246	685, −2.8, 321 / 678, 279	1.339 px

Dataset 2 is the striking one: $f_u = 672$ against the report’s 671, $f_v = 681$ against 681, $c_u = 321$ against 319. That column is my 2019 corners fed to OpenCV’s calibration, and the Rust here reproduces it to within a pixel on three of five parameters: two independent implementations of a detector and two independent implementations of a calibrator, agreeing. Dataset 3 is looser, about 1% on the focal lengths, and $c_v$ is the usual offender.

What this does not do

Worth saying plainly, because a calibration tool that overstates itself is worse than none:

• No lens distortion. Not radial, not tangential, no coefficients at all. A wide-angle or cheap lens will not be well described by this K at the edges of the frame.
• No extrinsics. get_camera_calib_m returns K and stops. Zhang’s method can recover $\mathbf{R}$ and $\mathbf{t}$ per view (they are one $\mathbf{K}^{-1}\mathbf{H}$ away) but this implementation never does, and the report lists that as a known drawback of the planar approach as implemented here.
• No uncertainty. There is no covariance on the recovered parameters, so "$c_v = 117$" and "$c_v = 178$" look equally confident and are not.

The Rust

The maths lives in two new modules of calib-wasm, the crate the detector post established, and it stayed dependency-free: the shapes here are so small that pulling in a linear-algebra library would have cost more bytes than writing the two decompositions.

• linalg.rs: a one-sided Jacobi SVD (rotate pairs of columns until they are mutually orthogonal; the column norms are then the singular values and the accumulated rotation is $\mathbf{V}$), the same loop again in f32, the $\mathbf{A}^\top\mathbf{A}$ shortcut, a textbook Cholesky that returns None on the first non-positive pivot, and a 9×9 Gaussian solve for the LM normal equations.
• zhang.rs: normalisation, the DLT, LM, the conic solve, K, and the C ABI.

Where it differs from the NumPy, and why:

	NumPy	Rust
SVD	np.linalg.svd, LAPACK gesdd	hand-rolled one-sided Jacobi; a f32 variant and an $\mathbf{A}^\top\mathbf{A}$ variant for the widget’s solver switch
LM	scipy.optimize.root(method='lm'), MINPACK lmder, on the squared residual with a mismatched Jacobian	damped normal equations on the unsquared residual with the matching Jacobian
Normalisation	always on	switchable
Failure	the string "Did not Converge"	a Status enum the JS side reads as an integer
Diagnostics	discarded	condition numbers, the conic system’s singular spectrum, B, L and per-view RMS all reported

cargo test -p calib-wasm covers the SVD against a closed-form decomposition built from two explicit rotations, Cholesky against the standard integer hand case, a synthetic camera recovered from noiseless projections to within $10^{-4}$ of a K with non-zero skew, the condition-number gap with normalisation on and off, that the naive solver only survives when normalised, that LM never makes the RMS worse, and that a sign-flipped null vector breaks Cholesky unless the B[2,2] division is applied. The crate grew from 40.7 kB to 60.0 kB.

The next thing to build on this is obvious and I did not build it: the detector and the calibrator on one page, pointed at a webcam, so you print a board and walk away with your own camera’s numbers. That needs a capture protocol, a pose-diversity indicator so you are not quietly collecting eight fronto-parallel views, and the honesty section above stapled to the bottom of it.