Brown-Conrady distortion fitting from point pairs alone

You have about 100 pairs of matching points — one set measured through your camera and lens, the other the true positions from a precision stage, a fixture or a CAD drawing — and you need the Brown-Conrady distortion coefficients xc, yc, k1, k2, p1 and p2. A curve-fitting app wants a single input variable and a camera calibrator wants images, so neither accepts raw pairs. This guide shows how to fit the six coefficients with a least-squares solver you already have: the equations, how many pairs you need and where to put them, a worked example with the arithmetic shown, and the point at which better optics beat a better model.

Can you fit Brown-Conrady distortion from point pairs alone?

Yes — the six coefficients are six unknowns and 100 point pairs give 200 equations, so an ordinary nonlinear least-squares solver converges on them in seconds on a laptop.

You do not need a checkerboard or a calibration image. You need three things: a model, a residual (predicted minus measured) and a solver. The lsqcurvefit function and scipy.optimize.least_squares were built for exactly this problem shape. Give each the residual function and a starting guess and it returns the six parameters. What your pairs must supply is coverage: they have to sample the field, not a line through it, and the next two sections put numbers on that.

Which equations connect the two point sets?

The forward Brown-Conrady map: it starts from the undistorted pair (xu, yu), applies radial terms k1 and k2 and tangential terms p1 and p2 about the distortion centre (xc, yc), and predicts the measured pair (xd, yd).

Work in normalized units: subtract the nominal centre and divide by half the image width, so 1.0 is the right edge of a landscape sensor. Scale matters for seeding. The example lens below has k1 = −0.052 normalized, but k1 = 0.052/1224² ≈ 3.5×10⁻⁸ per pixel² in raw pixels, and a solver started at zero crawls toward numbers that small. A k1 of −0.052 also tells you the distortion size at a glance: about 12% radial pull-in at the field corner (0.052×1.70 + 0.011×2.89 = 0.12).

Conventions differ between libraries: the sign of k1, the order of parameters, even which pair counts as input. The interpreting distortion coefficients on an unrectified image walkthrough explains what each term does to an image, and garbage depth from unknown distortion model coefficients is what happens downstream when the assumed model does not match the coefficients at hand. Match conventions before you port anything.

How many point pairs does a stable fit need?

Three pairs is the mathematical floor — six unknowns, two equations per pair — but plan on 100 across the whole field. Coefficient scatter shrinks roughly with √N, so 100 pairs are about 1.8× steadier than 30 (√(100/30) ≈ 1.8).

Placement decides more than count. At least a quarter of the pairs should sit in the outer third of the field, because the corners carry most of the information about k1 and k2. Keep a few near the centre, because they pin xc and yc. Never fit a set that lies along one line: a collinear set cannot separate a centre shift from a k1 change at all. Hold out 10 pairs before fitting so you can test on data the solver never saw. And measure the pairs in the state you will run: the same area scan camera, lens and focus setting as production, because a refocus or a bumped mount changes the geometry the coefficients describe.

How do you set up and run the least-squares fit?

Five steps: normalize both point sets, write the residual from the equations above, seed the solver at zero, bound the parameters, and judge the run by its RMS residual in pixels — the solve itself is one call.

1. Normalize: subtract the nominal centre and divide by half the image width, for both point sets.
2. Write the residual: for each pair, predict (xd, yd) from (xu, yu) with the current parameters and subtract the measured values — 2N residuals for N pairs.
3. Seed: theta0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] for (xc, yc, k1, k2, p1, p2), all float64.
4. Bound: |k1|, |k2| ≤ 0.5; |p1|, |p2| ≤ 0.01; centre inside the image.
5. Run, then report RMS = √(Σresiduals²/N) in normalized units; multiply by the half-width in pixels to state it in pixels.

import numpy as np
from scipy.optimize import least_squares

def residual(theta, xu, yu, xd, yd):
    # theta = [xc, yc, k1, k2, p1, p2], float64, normalized units
    xc, yc, k1, k2, p1, p2 = theta
    xn, yn = xu - xc, yu - yc
    r2 = xn**2 + yn**2
    rad = 1.0 + k1*r2 + k2*r2**2
    dx = p1*(r2 + 2.0*xn**2) + 2.0*p2*xn*yn
    dy = p2*(r2 + 2.0*yn**2) + 2.0*p1*xn*yn
    return np.concatenate([xn*rad + dx - xd, yn*rad + dy - yd])

theta0 = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0])  # start: identity map
sol = least_squares(residual, theta0, args=(xu, yu, xd, yd))
rms_px = np.sqrt(np.mean(sol.fun**2)) * 1224.0      # residual in pixels

Test the plumbing before you trust it. 100 synthetic pairs built from the worked-example coefficients with 0.05 px of random noise round-trip to k1 = −0.052 ± 0.001 and an RMS residual of 0.05 px — the fit recovers what went in. The data, variable names and noise figure here are illustrative; only the solver call is a library interface.

What starting values and bounds make the solver converge?

Seed the centre at the nominal image centre with every coefficient at zero, and bound |k1| and |k2| to 0.5 and |p1| and |p2| to 0.01; from that start a field-wide pair set converges on the first call.

At zero the model is the identity map, so the first residual is just your data and every gradient is well defined. The bounds encode physics: tangential terms from element decentering are typically 100× smaller than the radial ones, so a p that runs to its 0.01 bound means mislabeled pairs, not a strange lens. If a run stalls, suspect the points before the seed. A line of pairs, a swapped (x, y) column or duplicated targets cause far more failures than starting values do.

How do you check a fit before you trust it?

Judge it on pairs the solver never saw. In the worked example below the held-out RMS is 0.68 px, which is 0.033 mm at a pixel size of 0.049 mm — comfortable for a ±0.1 mm gauge and far too coarse for ±0.02 mm.

The setup: a 2448 × 2048 image over a 120 mm field of view, so one pixel is 120/2448 = 0.049 mm and the normalization half-width is 1224 px. (The object-space pixel size calculation generalizes this step.) Say the fit over 90 pairs returns k1 = −0.052, k2 = 0.011, p1 = −0.0002, p2 = 0.0001 and a centre at (0.000, 0.000) — example values chosen to make the arithmetic easy to follow. For one held-out pair with true position (0.600, 0.350):

• r² = 0.600² + 0.350² = 0.4825
• radial factor = 1 + (−0.052)(0.4825) + (0.011)(0.4825²) = 1 − 0.0251 + 0.0026 = 0.9775
• predicted xd = 0.9775 × 0.600 − 0.0002 = 0.5863 and yd = 0.9775 × 0.350 − 0.00001 = 0.3421
• the bench measured (0.5858, 0.3424), so the residuals are −0.0005 and +0.0003 normalized
• in pixels: −0.0005 × 1224 = −0.61 px and +0.0003 × 1024 = +0.31 px
• point RMS = √(0.61² + 0.31²) = √(0.372 + 0.096) = 0.68 px, and 0.68 px × 0.049 mm/px = 0.033 mm

The decision rule: keep the held-out RMS below about a fifth of the tolerance band. A ±0.1 mm gauge has a 0.2 mm band, so the limit is 0.04 mm and 0.033 mm passes. A ±0.02 mm gauge limits you to 0.008 mm and the same fit fails — for that part, software is the wrong tool, as the last section covers.

Two more checks take a minute. Convert k1 and k2 to corner distortion — 0.052 × 1.70 + 0.011 × 2.89 = 0.12, about 12% — and ask whether your lens family produces that. A corrected FA lens prints under 0.3%, so a fit that returns −0.052 on such a lens means bad pairs, not barrel distortion. Then plot residuals against radius: structure that grows toward one side means the centre is absorbing a decentered element. For where your RMS should land next to calibration targets, the calibration residual average and maximum targets article gives the values to aim for.

How do you undistort new points with the fitted coefficients?

Run the model backwards: the equations map ideal to measured, so correcting a measurement means inverting them — a fixed-point update that typically settles in three to five passes per point.

Start each point from xu = xd, yu = yd and iterate until the update moves it by less than 0.001 px. For |k1| below about 0.2 that takes three to five passes and a handful of multiplies each, so even a million points undistort in well under a second. Then feed the corrected coordinates downstream exactly as before.

The fisheye undistortImage model-matching failure covers the same trap from the function side: an undistort call only does what the model it was built for allows.

When is fitting from point pairs the wrong fix?

Two cases: you can take a calibration image after all, in which case a standard calibrator solves intrinsics and distortion together; or the fitted error exceeds a fifth of your tolerance band, in which case you buy optics that do not produce the error.

A checkerboard beats a point-pair fit whenever it is physically possible: it also returns focal length, principal point and pixel scale, which the pair fit cannot see. But when the tolerance rules software out, change the optics. The LS-FFL-0567 prints distortion under 0.3%. On the worked example's 120 mm field that bounds corner error near 0.36 mm before correction, and the fit then cleans up whatever remains. For dimensional gauging, telecentric lenses remove the whole problem class: constant magnification means radial distortion is designed out, not fitted around, and the LS-TCL-0574 — 0.4x at a 227 mm working distance — is that kind of part.

If the numbers say software is enough, you are done: archive the pairs together with the coefficients, and re-validate the fit whenever the lens is refocused, swapped or knocked. If they say optics, size the field of view and tolerance in the lens finder, compare machine vision lenses by their printed distortion, and work the full camera-to-light chain in the lens configurator.