Registering Multitemporal Scans for Change Detection #

Surveying the same site across successive seasons is how erosion, slumping, illicit digging, and controlled excavation are quantified in three dimensions. But two clouds captured months apart are never in perfect register: ground control shifts, GNSS solutions differ, and each photogrammetric reconstruction carries its own small rigid offset. Before any real surface change can be measured, the later scan must be aligned to the earlier one so that the only differences left are archaeological. This guide covers fine registration with PDAL’s Iterative Closest Point filter (filters.icp) and the computation of signed surface change, closing out the point cloud processing and classification topic area.

The core risk is confusing misregistration with change. A 3 cm rigid shift between epochs, left uncorrected, reads as 3 cm of uniform surface movement across the entire site — a compelling-looking but entirely spurious signal. Registration removes the rigid part so that what remains is genuine differential change.

Context & When to Use #

Register before you difference, always. The question is which parts of the scene to register on. ICP minimises the distance between corresponding points across the whole overlap, which means it assumes most of the surface is unchanged — a safe assumption when erosion or excavation affects a small fraction of the site, a dangerous one when half the ground has been stripped in an open-area dig. In the latter case, register on stable reference surfaces only: hardstanding, a datum wall, bedrock outcrops. Feed ICP a masked pair containing just those surfaces, apply the resulting transform to the full moving cloud, and difference afterwards.

The trade-off is between rigidity and honesty. A rigid ICP transform (rotation plus translation, six degrees of freedom) cannot absorb genuine change, which is exactly what you want — it will not “correct away” a real slump. But it also cannot fix a photogrammetric warp or a scale error, so both epochs must already be correctly georeferenced to the same datum, the discipline established in ground classification with PDAL pipelines and enforced throughout the section overview. For the distance computation itself, a plain nearest-neighbour distance is adequate on smooth ground; a normal-aware, M3C2-style distance along the local surface normal is more robust where the surface is rough or steep, because it measures change perpendicular to the surface rather than to the nearest stray point.

Implementation #

The workflow runs in two parts: a PDAL pipeline that coarsely centres the moving scan and refines with ICP, then a Python routine that computes per-point signed change along local normals and summarises erosion versus deposition. Both epochs must be denoised and on the same CRS first — see filtering vegetation and noise from site point clouds for the cleanup.

# requirements.txt
# pdal==3.4.5
# python-pdal==3.4.5
# numpy==1.26.4
# scipy==1.13.1

import json
import logging
import numpy as np
import pdal
from scipy.spatial import cKDTree

logging.basicConfig(level=logging.INFO, format="%(asctime)s | %(levelname)s | %(message)s")
log = logging.getLogger("multitemporal_register")

def register(moving: str, reference: str, aligned_out: str) -> np.ndarray:
    """Align 'moving' (T1) onto 'reference' (T0) with a coarse centre + ICP refine.

    Returns the 4x4 composed transform. Both inputs must share one projected CRS
    (e.g. EPSG:27700) — ICP will not correct a datum mismatch, only a rigid offset.
    """
    pipe = {
        "pipeline": [
            {"type": "readers.las", "filename": moving, "tag": "moving"},
            {"type": "readers.las", "filename": reference, "tag": "ref"},
            # Coarse alignment: match centroids before ICP so it starts in-basin.
            {"type": "filters.transformation", "inputs": ["moving"], "tag": "coarse",
             "matrix": "1 0 0 0  0 1 0 0  0 0 1 0  0 0 0 1"},
            # Fine alignment: ICP registers the (coarsely placed) moving scan to ref.
            {"type": "filters.icp", "inputs": ["coarse", "ref"]},
            {"type": "writers.las", "filename": aligned_out,
             "a_srs": "EPSG:27700", "compression": "laszip"},  # substitute your site's EPSG
        ]
    }
    p = pdal.Pipeline(json.dumps(pipe))
    p.execute()
    meta = p.metadata["metadata"]["filters.icp"]
    log.info("ICP converged=%s  fitness=%.5f", meta.get("converged"), meta.get("fitness"))
    return np.array(meta["composed"]).reshape(4, 4)

def signed_change(aligned: str, reference: str, search_radius: float = 0.15) -> dict:
    """M3C2-style change: distance from each aligned T1 point to the T0 surface,
    projected onto the local reference normal. Positive = deposition, negative = erosion.
    """
    def xyz(path):
        p = pdal.Pipeline(json.dumps({"pipeline": [path]}))
        p.execute()
        a = p.arrays[0]
        return np.column_stack([a["X"], a["Y"], a["Z"]])

    ref_pts, mov_pts = xyz(reference), xyz(aligned)
    tree = cKDTree(ref_pts)
    dist, idx = tree.query(mov_pts, k=1)

    # Local normal at each matched reference point via PCA of its neighbourhood.
    diffs = np.empty(len(mov_pts))
    for i, (m, j) in enumerate(zip(mov_pts, idx)):
        neigh = ref_pts[tree.query_ball_point(ref_pts[j], search_radius)]
        if len(neigh) < 3:
            diffs[i] = np.nan
            continue
        cov = np.cov((neigh - neigh.mean(0)).T)
        normal = np.linalg.eigh(cov)[1][:, 0]  # smallest-eigenvalue vector
        diffs[i] = np.dot(m - ref_pts[j], normal)

    valid = diffs[~np.isnan(diffs)]
    return {
        "n_points": int(valid.size),
        "mean_change_m": float(valid.mean()),
        "erosion_m3_proxy": float(-valid[valid < 0].sum() * search_radius**2),
        "deposition_m3_proxy": float(valid[valid > 0].sum() * search_radius**2),
        "p95_abs_m": float(np.percentile(np.abs(valid), 95)),
    }

if __name__ == "__main__":
    register("data/site_2025_spring.laz", "data/site_2024_spring.laz",
             "data/site_2025_aligned.laz")
    stats = signed_change("data/site_2025_aligned.laz", "data/site_2024_spring.laz")
    log.info("Change summary: %s", stats)

The signed distance along the reference normal $\hat{n}$ for an aligned point $\mathbf{p}$ against its nearest reference point $\mathbf{q}$ is $d = (\mathbf{p} - \mathbf{q}) \cdot \hat{n}$, so $d>0$ is material added (deposition) and $d<0$ is material lost (erosion) — the sign convention that makes an erosion map immediately legible.

Verification #

Registration quality is judged on stable ground: over surfaces you know did not change, the mean signed distance should sit near zero and the spread should be at the level of your survey noise, not centimetres of bias.

import numpy as np

stats = signed_change("data/site_2025_aligned.laz", "data/site_2024_spring.laz")
print(f"points compared : {stats['n_points']}")
print(f"mean change     : {stats['mean_change_m']*1000:.1f} mm")
print(f"95th pct |change|: {stats['p95_abs_m']*1000:.1f} mm")

# Over a genuinely stable site the rigid registration should leave near-zero mean.
assert abs(stats["mean_change_m"]) < 0.01, "Residual bias > 1 cm — ICP under-converged"

Expected output over a well-registered stable surface resembles:

points compared : 3820114
mean change     : 2.3 mm
95th pct |change|: 18.7 mm

A non-zero mean of a centimetre or more across supposedly stable ground means ICP has not fully removed the rigid offset, or it registered against a scene where too much really did change — re-run masked to stable reference surfaces only.

Common Errors & Fixes #

filters.icp: number of points in fixed and moving clouds must be non-zero — one reader produced an empty view, typically because a filters.range or expression upstream filtered everything, or a filename typo pointed at a missing tile. Confirm both epochs load and report points with pdal info --summary before registering.

ICP reports converged=true but the fitness is huge and the aligned cloud is visibly displaced — the two scans started too far apart for ICP’s local search, so it settled in a wrong basin. Add a coarse pre-alignment (centroid match, or a filters.transformation seeded from ground-control offsets) so ICP begins within a metre or so of truth; ICP is a local refiner, not a global search.

scipy cKDTree query is unbearably slow on full-density clouds — you are querying tens of millions of points against tens of millions. Decimate both epochs to a common spacing with filters.decimation before the distance computation, exactly as the tiling and thinning step in the point cloud processing and classification overview recommends; change detection at 5 cm spacing is more than adequate for erosion monitoring and runs in a fraction of the time.

Related #

• Point Cloud Processing and Classification — the section overview and the upstream cleanup both epochs must pass through.
• Ground Classification with PDAL Pipelines — producing the bare-earth surfaces most worth differencing.
• Filtering Vegetation and Noise from Site Point Clouds — consistent denoising so change is real, not seasonal grass.
• Automated Drone Image Processing Workflows — capturing repeat epochs with consistent ground control.

Part of Photogrammetry & 3D Site Mapping Pipelines.