Computing RMSE for Control Point Alignment in Python

To measure an alignment honestly, split the surveyed points into a control subset used for the fit and a check subset the fit never sees, transform the check points, and report the maximum residual and its index alongside the root-mean-square error. An RMSE computed on the points that were fitted is a measure of the optimiser, not of the transform. This page is part of Scale and Rotation Synchronization.

Why the Split Matters More Than the Statistic

A similarity transform in the plane has four unknowns and in three-space seven. Given exactly the minimum number of points, the solve reproduces them exactly and the residual is zero — whether or not the correspondences were right, whether or not one point was mislabelled, whether or not the two datasets are on the same datum at all.

Adding points reduces that effect without removing it. The solver is still minimising the residual on those points, so their residual remains the quantity that was optimised rather than an independent test. Only a point excluded from the fit measures prediction.

The second half of the argument is about which statistic to act on. The root-mean-square error describes the bulk and is the number usually quoted. The maximum describes the worst case and is the number that identifies problems, because the characteristic failure in this domain is not general imprecision but a single mismatched point — the same monument identified as two different points in the two datasets. A set with an RMSE of 0.03 m and a maximum of 0.42 m is not a slightly imprecise transform; it is a good transform with one bad correspondence, and the average conceals exactly the observation that finds it.

Production-Ready Script

# numpy>=1.24, Python 3.9+
from __future__ import annotations

from dataclasses import dataclass
import numpy as np


@dataclass(frozen=True)
class ResidualReport:
    n_control: int
    n_check: int
    rmse: float
    mean: float
    p95: float
    maximum: float
    worst_index: int
    passed: bool
    tolerance: float

    def summary(self) -> str:
        verdict = "PASS" if self.passed else "FAIL"
        return (f"{verdict}  rmse={self.rmse:.4f} m  max={self.maximum:.4f} m "
                f"(check point {self.worst_index})  tol={self.tolerance:.3f} m")


def split_control_check(points: np.ndarray, *, holdout: float = 0.34, seed: int = 0):
    """Deterministic split — the same points every run, so results compare."""
    n = len(points)
    rng = np.random.default_rng(seed)
    idx = rng.permutation(n)
    n_check = max(2, int(round(n * holdout)))
    if n - n_check < 3:
        raise ValueError(f"{n} points is too few to both fit and check")
    return idx[n_check:], idx[:n_check]        # control, check


def residuals(transformed: np.ndarray, surveyed: np.ndarray) -> np.ndarray:
    a = np.asarray(transformed, dtype=float)
    b = np.asarray(surveyed, dtype=float)
    if a.shape != b.shape:
        raise ValueError(f"shape mismatch {a.shape} vs {b.shape}")
    return np.linalg.norm(a - b, axis=1)


def report(
    transformed_check: np.ndarray,
    surveyed_check: np.ndarray,
    *,
    n_control: int,
    tolerance_m: float,
) -> ResidualReport:
    d = residuals(transformed_check, surveyed_check)
    worst = int(d.argmax())
    return ResidualReport(
        n_control=n_control,
        n_check=int(d.size),
        rmse=float(np.sqrt((d ** 2).mean())),
        mean=float(d.mean()),
        p95=float(np.percentile(d, 95)),
        maximum=float(d[worst]),
        worst_index=worst,
        # Both statistics have to clear: the bulk AND the worst case.
        passed=bool(np.sqrt((d ** 2).mean()) <= tolerance_m and d[worst] <= tolerance_m * 2),
        tolerance=tolerance_m,
    )


if __name__ == "__main__":
    surveyed = np.array([[...]])          # (n, 2) or (n, 3) surveyed positions
    local = np.array([[...]])             # the same points in the source frame
    control_idx, check_idx = split_control_check(surveyed)
    # fit(...) on local[control_idx] -> surveyed[control_idx], then:
    # moved = apply(fit, local[check_idx])
    # print(report(moved, surveyed[check_idx],
    #              n_control=len(control_idx), tolerance_m=0.05).summary())

Key implementation notes:

• The split is seeded, so two runs on the same data compare. An unseeded split makes every re-run a different measurement.
• passed requires both the RMSE and the maximum to clear, with a wider allowance on the maximum. A single criterion on either alone lets through the failure the other catches.
• worst_index is returned rather than just the value, because the index is what an investigation needs.
• The split refuses to leave fewer than three points for the fit. A degenerate fit reports a beautiful residual.

Compatibility Matrix

Fallback Strategies

1. Too few points to split. Fit on all of them and state plainly that the residual is not independent evidence. Do not report it as a validation.

2. Maximum far above the RMSE. One mismatched correspondence. Inspect the point named by worst_index, refit without it, and see whether the maximum drops to the bulk level — if it does, the point was the problem.

3. Every residual is large and similar. Not a correspondence problem but a systematic one: a datum difference, a unit error or an unapplied vertical offset. The uniformity is the diagnosis.

4. Residuals grow with distance from the centroid. A scale error. Check whether the solve was allowed to fit scale, and whether it should have been.

5. The RMSE passes and the model is visibly wrong. Check for a mirrored transform — a reflection fits control points as well as a rotation does. That is the subject of the sibling guide on detecting mirrored transforms.

FAQ

Related Pages

• Scale and Rotation Synchronization — parent reference on the similarity solve these residuals measure
• Aligning BIM Models with GIS Survey Data — the SVD solve this measurement is applied to
• Detecting Mirrored Transforms in Python — a fault that produces a good RMSE and a wrong model