Aligning Two Epochs with ICP Before Differencing

The difference map is one colour. Every cell across two hundred hectares shows between three and five centimetres of apparent subsidence, including the concrete weighbridge, the office roof and the public road along the boundary. Nothing on that site subsided; the two surveys are simply four centimetres apart in the vertical, and the map is showing that offset rather than the site.

Fixing it is a registration problem, and it is one where the obvious approach — run ICP on everything and let it sort itself out — produces a subtly worse result than doing nothing. This page covers doing it correctly, and reading the answer as a diagnostic rather than accepting it silently.

# Why registering on everything is worse than it looks

Iterative closest point finds the rigid transform minimising the distance between two point sets. If a site has genuinely changed — a bench cut back, a stockpile grown, a cell filled — those changed regions are part of the cost function, and the optimiser will happily rotate and translate the whole survey to reduce the disagreement they cause.

The result is a transform that splits the difference: the changed area reports less change than occurred, and the unchanged area acquires a small artificial change in the opposite direction. On a site where five percent of the area changed by two metres, the whole-site fit can absorb several centimetres of that into the transform — which is the same order as the real detection limit.

The fix is simple to state and requires a site decision: register using only ground both surveys agree cannot have moved. A polygon of hardstanding, roads, building roofs and any permanent structures, established once with the site and reused every epoch.

Figure 1 — Why the default is the wrong default. The bias is in the direction that matters most.

# Minimal reproducible solution

import numpy as np
import open3d as o3d


def register_on_stable(moving: np.ndarray, reference: np.ndarray,
                       moving_stable: np.ndarray, reference_stable: np.ndarray,
                       *, max_distance: float = 0.5,
                       voxel: float = 0.25) -> dict:
    """Point-to-plane ICP restricted to stable ground, with diagnostics.

    Point-to-plane rather than point-to-point: photogrammetric surfaces are
    dense and locally planar, and the plane metric converges in far fewer
    iterations without sliding along the surface.
    """
    def to_o3d(arr):
        p = o3d.geometry.PointCloud()
        p.points = o3d.utility.Vector3dVector(arr)
        return p.voxel_down_sample(voxel)

    src = to_o3d(moving[moving_stable])
    tgt = to_o3d(reference[reference_stable])
    tgt.estimate_normals(
        o3d.geometry.KDTreeSearchParamHybrid(radius=voxel * 4, max_nn=30))

    res = o3d.pipelines.registration.registration_icp(
        src, tgt, max_distance, np.eye(4),
        o3d.pipelines.registration.TransformationEstimationPointToPlane(),
        o3d.pipelines.registration.ICPConvergenceCriteria(max_iteration=60))

    t = np.asarray(res.transformation)
    translation = t[:3, 3]
    rot = np.degrees(np.arccos(np.clip((np.trace(t[:3, :3]) - 1) / 2, -1, 1)))
    return {"transform": t,
            "translation_m": translation.tolist(),
            "vertical_shift_m": float(translation[2]),
            "horizontal_shift_m": float(np.linalg.norm(translation[:2])),
            "rotation_deg": float(rot),
            "fitness": float(res.fitness),
            "inlier_rmse": float(res.inlier_rmse)}

Voxel downsampling before the fit is worth doing for a reason beyond speed: it equalises the spatial weighting. A dense patch of hardstanding with ten times the point density of the rest would otherwise dominate the cost function and pull the fit toward matching that one area perfectly.

# Reading the transform as a diagnostic

The transform is not just a correction to apply; it is a measurement of how far apart the two surveys’ georeferencing is, and its magnitude says whether applying it is appropriate at all.

def assess_transform(diag: dict, *, site_extent_m: float,
                     max_vertical: float = 0.15,
                     max_horizontal: float = 0.20,
                     max_rotation_deg: float = 0.05) -> dict:
    """Decide whether this registration should be applied or investigated."""
    verdicts = []
    if abs(diag["vertical_shift_m"]) > max_vertical:
        verdicts.append(
            f"vertical shift {diag['vertical_shift_m']:+.3f} m is larger than two "
            "good surveys should differ — check the vertical datum, not the ICP")
    if diag["horizontal_shift_m"] > max_horizontal:
        verdicts.append(
            f"horizontal shift {diag['horizontal_shift_m']:.3f} m suggests different "
            "ground control between the epochs")
    # A small rotation over a large site is a large displacement at the edges.
    edge_effect = np.radians(diag["rotation_deg"]) * site_extent_m / 2
    if diag["rotation_deg"] > max_rotation_deg:
        verdicts.append(
            f"rotation {diag['rotation_deg']:.4f}° moves the site edge by "
            f"{edge_effect:.3f} m — this is a survey problem, not an alignment one")
    if diag["fitness"] < 0.6:
        verdicts.append(f"only {diag['fitness']:.0%} of stable points found a "
                        "correspondence — the stable polygon may not overlap")
    return {"apply": not verdicts, "findings": verdicts,
            "edge_displacement_m": float(edge_effect)}

Converting the rotation into an edge displacement is the step that makes it interpretable. A rotation of 0.03° sounds negligible and, across a site 1.5 km wide, moves the far corner by 40 cm. Quoting the angle alone has let more than one survey ship with a real georeferencing fault presented as a successful alignment.

# The sequence in one view

The order below is what keeps an alignment honest, and every stage exists because skipping it produces a specific, recognisable failure.

Figure 3 — The stable subset drives the fit; the whole cloud receives the transform.

# Edge-case matrix

The flat-site row is a genuine limitation. Point-to-plane ICP on a perfectly flat surface has no information about horizontal position — every horizontal translation fits equally well — so the solver will return an arbitrary one. Where that is the case, solve for a vertical shift only:

import numpy as np


def vertical_only_shift(moving_z: np.ndarray, reference_z: np.ndarray) -> float:
    """Median vertical offset over co-located stable cells, no rotation."""
    d = reference_z - moving_z
    return float(np.median(d[np.isfinite(d)]))

# Verification snippet

import numpy as np


def verify_registration(diff_after: np.ndarray, stable_mask: np.ndarray) -> dict:
    """After registration, stable ground must be centred on zero."""
    s = diff_after[stable_mask & np.isfinite(diff_after)]
    bias = float(np.median(s))
    sigma = float(np.median(np.abs(s - bias)) * 1.4826)
    # Residual tilt: fit a plane to the stable residual and check its gradient.
    return {"stable_bias_m": bias, "stable_sigma_m": sigma,
            "ok": abs(bias) < 0.01,
            "note": ("registration successful" if abs(bias) < 0.01
                     else "a residual offset remains — check the vertical datum")}

Figure 2 — Always convert the rotation into a distance before deciding whether the alignment is routine.

# Choosing and maintaining the stable set

Everything above depends on a polygon somebody has to draw, and drawing it well is a site question rather than a processing one. Four properties make a stable set work.

It must be genuinely stable. Concrete hardstanding, adopted roads, building roofs and permanent plant foundations qualify. Compacted haul roads do not — they are regraded constantly. Car parks are borderline: the surface is stable, the vehicles on it are not, so a car park only works if the classification removed the vehicles first.

It must be distributed. A stable set concentrated in one corner constrains the transform there and lets it rotate freely about that point, which produces exactly the edge displacement described above. Three or four separated patches around the site perimeter constrain rotation far better than one large central area of the same total size.

It must be large enough. Below a few thousand points after downsampling, the fit is noisy and the fitness statistic becomes unreliable. A useful floor is roughly one percent of the site area, split across the patches.

It must be reviewed occasionally. Sites change. A building is demolished, a road is resurfaced, a yard is extended. A stable set defined at the start of a three-year monitoring contract and never revisited will eventually include something that moved, and the symptom — a slowly growing registration residual with no obvious cause — is hard to diagnose after the fact.

import numpy as np


def audit_stable_set(diff_after: np.ndarray, stable_labels: np.ndarray,
                     limit: float) -> list[dict]:
    """Per-patch residual, so a patch that stopped being stable is visible."""
    out = []
    for label in np.unique(stable_labels[stable_labels > 0]):
        m = (stable_labels == label) & np.isfinite(diff_after)
        if m.sum() < 100:
            continue
        r = diff_after[m]
        bias = float(np.median(r))
        out.append({"patch": int(label), "cells": int(m.sum()),
                    "bias_m": bias,
                    "suspect": bool(abs(bias) > limit)})
    return sorted(out, key=lambda d: -abs(d["bias_m"]))

Running that audit every epoch and flagging any patch whose bias exceeds the detection limit turns stable-set maintenance from an annual chore into an automatic one. When a patch starts drifting, the pipeline says so in the month it happens rather than after the trend has quietly corrupted a year of differences.

# When to escalate

• The vertical shift exceeds about 15 cm. Two surveys with competent ground control do not differ by that much. The cause is almost always a vertical datum or geoid model difference, covered in geoid models and vertical datum automation, and hiding it inside a transform makes the series wrong from then on.
• Fitness is low and the stable polygon is correct. The two surveys may not overlap as much as assumed, or one may have failed to reconstruct the stable areas. Check coverage before adjusting parameters.
• The residual after registration has a spatial pattern. A rigid transform cannot remove a dome or a bowl, which is a bundle-adjustment fault — see diagnosing doming and bowl effect in reconstructions.

# Related

• Change Detection Between Survey Epochs
• Computing M3C2 Distances on Drone Point Clouds
• Diagnosing Systematic Tilt in GCP Residuals

← Change Detection Between Survey Epochs