Known-intrinsics camera calibration: solving a 9-DOF pinhole model

You have one camera, one view, measured 3D world points and the image points they project to, and the geometry is known to be clean: square pixels, zero skew, zero distortion. That is nine unknowns — a shared focal length in pixels, a principal point, three rotations and a translation of three. The planar estimator in a vision toolbox will not fit that model: it wants many views of one flat board and it estimates distortion by default. This guide solves the problem directly, with a direct linear transform for the projection matrix, an RQ decomposition into intrinsics and pose, and a refinement over exactly nine parameters — code, a worked example and the pass thresholds that decide the result.

What you are fitting: nine parameters, not eleven

A 3×4 projection matrix has twelve entries and is fixed only up to scale, so an unconstrained fit carries eleven degrees of freedom. The general calibration matrix K accounts for five of them — fx, fy, cx, cy and skew — and the camera pose adds six more. Your two assumptions delete exactly two: square pixels force fx = fy, and zero skew removes the off-diagonal term. Three intrinsics plus six extrinsics leaves nine numbers, each physically meaningful, and none free to absorb error from the others.

Here f is the focal length in pixels, (cx, cy) the principal point where the optical axis meets the sensor, R a 3×3 rotation and t the world-to-camera translation in millimetres. The scalar s is the depth of the point along the viewing axis; it divides out in the projection. For when the pinhole relation itself is exact enough to use, read undistorted images and the pinhole model: when x = fX/Z holds.

Why the planar toolbox estimator will not do it

The estimator most toolboxes put in front of you — estimateCameraParameters among them — implements Zhang's method: several views of one planar board, a closed-form initialisation, then a non-linear refinement that always includes a distortion model. It exposes no option to drop radial distortion or to fit fewer than four intrinsics, and your data is a single set of non-coplanar points, which a planar-method estimator is never asked to digest.

You have two ways forward. Patch the estimator source to skip its refinement call and accept an eleven-parameter closed form — then you alone maintain a modified library file, and the model still is not the one you specified. Or run the direct solve below: about thirty lines you own, fitting exactly the nine parameters you believe in. For collecting the correspondences themselves, measured scene points rather than a printed board, camera calibration with known points is the companion guide.

The solve: DLT, then a constrained refinement

Stage one is the direct linear transform. Write each correspondence as two linear equations in the twelve entries of P, stack every pair, and take the one-dimensional null space of the stack by singular value decomposition.

Six points give twelve equations for eleven unknowns: that is the algebraic minimum. Stage two imposes your model and polishes the answer. Split P into its left 3×3 block M and last column p4, run an RQ decomposition on M — scipy.linalg.rq returns an upper-triangular K and an orthogonal R — and recover t = K⁻¹ p4. Average the two decomposed focal lengths into one f0, then run a Levenberg–Marquardt refinement over just nine parameters, (f, cx, cy, rvec, tvec), minimising total squared reprojection error. The decomposition algebra is worked through in decomposing a projection matrix into intrinsics.

Fitting nine instead of eleven is not cosmetic. With twenty to thirty points the constrained model has slack where it should — averaging out centroid noise — and none where it should not, because there is no skew coefficient to soak up modelling error your sensor cannot produce.

Step by step

1. Collect at least 20 matched points. Measure or machine world points you can centroid to a tenth of a pixel, and match them by hand or by order. Keep them off one plane: two board positions with a Z separation of at least 5% of the working distance, points pushed into the frame corners.
2. Normalise both sets. Subtract the centroid, then scale: mean distance √2 from the origin for image points, √3 for world points. The solve is badly conditioned without it.
3. Build and solve. Stack the 2n × 12 system and take the right singular vector of the smallest singular value; the code below is the whole of it.
4. Denormalise. With T the image normalisation and U the world normalisation, P = T⁻¹ · P_n · U, where P_n is the matrix you just solved for.
5. Split and sign-fix. RQ-decompose, then flip signs until the diagonal of K is positive and det(R) = +1. Every point must then project with positive depth; a negative one means a sign error, not a broken camera.
6. Constrain and refine. Set f0 to the mean of the decomposed fx and fy, drop the skew, and refine (f, cx, cy, rvec, tvec) by Levenberg–Marquardt on reprojection error.
7. Gate the result against the thresholds two sections down, then stop. Re-measuring beats re-weighting.

import numpy as np

# X: (n, 3) float array of world points in mm, Hartley-normalised
# U: (n, 2) float array of image points in px, Hartley-normalised
A = np.zeros((2 * len(X), 12))
for i, ((x, y, z), (u, v)) in enumerate(zip(X, U)):
    A[2 * i, 0:4] = x, y, z, 1.0
    A[2 * i, 8:12] = -u * x, -u * y, -u * z, -u
    A[2 * i + 1, 4:8] = x, y, z, 1.0
    A[2 * i + 1, 8:12] = -v * x, -v * y, -v * z, -v
_, _, Vt = np.linalg.svd(A)
P = Vt[-1].reshape(3, 4)              # 11-DOF projection matrix, up to scale

# test that proves the solve: reproject point 0 and compare
h = P @ np.append(X[0], 1.0)
print(h[:2] / h[2], U[0])             # must agree to well under one pixel

Worked example: 12 MP camera, 16 mm lens, 250 mm working distance

The inputs, each with its source:

• Camera: 4000 × 3000 px on a 1/1.7-inch sensor, about 7.4 mm wide, so the pixel size is 7.4 mm / 4000 px = 1.85 µm; confirm the exact pitch on the sensor data sheet.
• Lens: 16 mm focal length, set to a middle f-number.
• Working distance: 250 mm from lens to the first target plane.
• Target: a 76 × 76 mm dot grid with 0.5 mm dots, at two Z positions 12.5 mm apart — 5% of the working distance.
• Points: 24 centroided dots, 12 per plane.

The rig, end to end:

The same 12 MP IMX226 sensor comes on GigE as the LS-CMR-0053 when USB3 cable reach is the limit.

First, the focal length you should get. In pixels, f = 16 mm / 0.00185 mm = 8,649 px. Write the acceptance window down before you solve: ±5% is 8,216 to 9,081 px. An answer outside that window is a scale error in the world points or a wrong pixel size far more often than a failed solver.

Second, the field the rig covers, so the board fills it. Ground sample size at 250 mm is s·Z/f = 0.00185 × 250 / 16 = 0.0289 mm per pixel, so 34.6 px per mm. The frame then covers 4000 × 0.0289 = 116 mm across and 3000 × 0.0289 = 87 mm tall — the same relation the working distance and field of view formula guide works from the other end.

The 76 × 76 mm board fills 66% of the width and 88% of the height, and every 0.5 mm dot spans 0.5 × 34.6 = 17 px, comfortable for sub-pixel centroiding. Starting from a target instead of a lens, the calculator asks the question backwards: a 116 mm field at 250 mm on a 7.4 mm sensor needs 250 × 7.4 / 116 ≈ 16 mm.

The solve on this data set returns f = 8,520 px, cx = 2014 px, cy = 1508 px, and an RMS reprojection error of 0.31 px. Reading the gates: f sits 1.5% below the 8,649 px expectation, inside 8,216 to 9,081; cx and cy sit 14 px and 8 px from the 2000, 1500 centre, inside the 200 px and 150 px that 5% of each dimension allows; 0.31 px passes the 0.5 px RMS gate and equals 0.31 × 0.0289 = 9 µm on the object. Decision: accept K and put it to work. If any gate fails, re-measure the board positions before touching the solver.

Pass thresholds before you accept the result

• Reprojection RMS under 0.5 px, and under 0.3 px when the calibration feeds measurement.
• f within ±5% of focal length over pixel size. Outside that, check the world-point scale and the pixel size first.
• Principal point within about 5% of the image centre. It is the noisiest of the nine; treat larger offsets as a symptom, not a finding.
• Positive depth on every point. The third row of P dotted with each homogeneous world point must exceed zero after sign fixing.
• Flat, directionless residuals. Error growing smoothly with radius from (cx, cy) is distortion asking for the parameter you removed.

When nine parameters is the wrong model

Three situations break this model, and each has a tell.

Radial residual pattern. A lens with real distortion bends the residual field smoothly outward from the principal point, and wide-angle glass — fields of view past about 45° — almost always does. Fit a distortion model instead, or characterise the lens first: building a lens distortion profile without an LCP file covers that.

Effectively coplanar points. One plane leaves the depth scale unobservable; the solve turns ill-conditioned or quietly wrong. Add the second plane, or change tools: with intrinsics already trusted, ground-plane pose from four known points recovers the pose from a homography.

A second camera. Do not solve nine parameters twice for a stereo pair. Freeze this calibration and solve the relative pose only — stereo calibration with fixed intrinsics, extrinsics only is the shorter, stiffer problem.

Rolling-shutter cameras are fine here while everything is static; a target that moves during capture smears the geometry in a way no parameter count fixes.

Ordering checks and the next step

Three checks before you order parts for a rig like this. Match the image circle to the sensor diagonal first and focal length second — a 16-22 mm circle leaves a 9.5 mm sensor diagonal in the quiet centre of the field, where distortion runs well below the printed corner figure. Match lens resolution to sensor resolution with margin: a 20 MP lens rating over a 12 MP sensor leaves room to spare. And choose a lens whose data sheet prints a distortion figure in the tenths of a percent, because that number is your zero-distortion argument.

The parts live in machine vision lenses, area scan cameras and calibration targets and lens accessories; the finders at the lens finder and the camera finder filter by sensor format, so the circle-to-diagonal check is two clicks. Size the whole chain — camera, lens, working distance — in the system configurator, and if your point set has a quirk worth a second opinion, contact the applications team.