Computing M3C2 Distances on Drone Point Clouds

A quarry face lost several cubic metres of rock between two flights. The difference raster shows almost nothing, because a vertical surface has no vertical difference: the material moved sideways, out of the face, and a top-down comparison is blind to it.

The standard answer is a cloud-to-cloud distance measured along the local surface normal, with an averaging radius that suppresses roughness and a per-point confidence interval that says whether each measurement is significant. This page implements that in Python for photogrammetric drone clouds, where the surfaces are dense, noisy and unevenly sampled — conditions that differ from the terrestrial laser scans the method was designed around.

It is the cloud-based half of change detection between survey epochs.

What the method actually computes

For each point in a reference cloud, the procedure is three steps. Estimate a surface normal from neighbours within a normal scale radius. Project a cylinder of a given projection radius along that normal into both clouds, and take the mean position of the points falling inside each. The signed distance between those two means is the measured change.

The averaging is what makes the method robust. A single nearest-neighbour distance between two noisy clouds measures noise; averaging tens or hundreds of points inside a cylinder reduces that noise by the square root of the count, and — importantly — the scatter inside each cylinder gives a direct estimate of how uncertain the result is at that point.

That per-point uncertainty is the property that distinguishes the method from a plain cloud-to-cloud distance. It means the output is not one map but two: the change, and whether the change is larger than that location’s own noise.

The three scales in an M3C2 measurement on a rough face A cross-section through a rough rock face with two epochs of points. A large circle marks the normal scale, within which neighbours are used to fit a local surface orientation; the resulting normal arrow points out of the face. A narrow cylinder aligned with that normal marks the projection scale, containing about forty points from each epoch. The mean position of each epoch's points inside the cylinder is marked, and the signed distance between them is the measured change of one point one metres. A note states that the scatter inside the cylinder gives the confidence interval directly. rough face, epoch 1 normal scale — orientation only local normal projection cylinder measured change 1.1 m between the two cylinder means scatter inside each cylinder gives the per-point confidence directly

Figure 1 — Three radii, three jobs. Confusing the normal scale with the projection scale is the most common implementation error.

Minimal reproducible solution

import numpy as np
from scipy.spatial import cKDTree


def m3c2(core_pts: np.ndarray, cloud_a: np.ndarray, cloud_b: np.ndarray,
         *, normal_scale: float = 1.0, projection_radius: float = 0.5,
         max_depth: float = 5.0, min_points: int = 8,
         registration_error: float = 0.02) -> dict:
    """Signed change along the local normal, with a per-point confidence.

    core_pts are the locations to measure at — usually a subsample of
    cloud_a, because computing at every point is wasteful and the result is
    spatially smooth anyway.
    """
    tree_a, tree_b = cKDTree(cloud_a), cKDTree(cloud_b)
    tree_n = cKDTree(cloud_a)

    n = len(core_pts)
    change = np.full(n, np.nan)
    uncert = np.full(n, np.nan)
    counts = np.zeros((n, 2), dtype=int)

    for i, p in enumerate(core_pts):
        # 1. local normal from the normal-scale neighbourhood
        idx = tree_n.query_ball_point(p, normal_scale)
        if len(idx) < min_points:
            continue
        nb = cloud_a[idx] - cloud_a[idx].mean(axis=0)
        normal = np.linalg.svd(nb, full_matrices=False)[2][-1]

        # 2. points within the cylinder, in each cloud
        stats = []
        for tree, cloud in ((tree_a, cloud_a), (tree_b, cloud_b)):
            cand = cloud[tree.query_ball_point(p, np.hypot(projection_radius,
                                                           max_depth))]
            if len(cand) == 0:
                stats.append(None)
                continue
            rel = cand - p
            along = rel @ normal
            radial = np.linalg.norm(rel - np.outer(along, normal), axis=1)
            inside = (radial <= projection_radius) & (np.abs(along) <= max_depth)
            if inside.sum() < min_points:
                stats.append(None)
                continue
            stats.append((float(along[inside].mean()),
                          float(along[inside].std(ddof=1)),
                          int(inside.sum())))

        if stats[0] is None or stats[1] is None:
            continue
        (ma, sa, na), (mb, sb, nb_) = stats
        change[i] = mb - ma
        uncert[i] = 1.96 * np.sqrt(sa ** 2 / na + sb ** 2 / nb_) + registration_error
        counts[i] = (na, nb_)

    return {"change": change, "uncertainty": uncert, "counts": counts,
            "significant": np.abs(change) > uncert}

Adding registration_error to the confidence rather than folding it into the scatter is deliberate. The cylinder statistics capture how well each cloud’s local surface is determined; they say nothing about whether the two clouds are correctly aligned with each other, which is a separate and usually larger term derived from the co-registration in aligning two epochs with ICP before differencing.

Choosing the two scales

The normal scale should span several times the surface roughness and be small compared with the features being measured. On a blocky rock face with 20 cm relief, a normal scale of about 1 m gives stable orientations; at 20 cm the normals follow individual blocks and rotate wildly between neighbouring core points, which produces a noisy and asymmetric change map.

The projection radius trades resolution against precision. Larger radii average more points, so the uncertainty falls as the square root of the count, but a feature smaller than the radius is smoothed away. A useful starting point is a radius equal to the smallest feature that matters, with the normal scale two to three times larger.

import numpy as np
from scipy.spatial import cKDTree


def suggest_scales(cloud: np.ndarray, sample: int = 20_000) -> dict:
    """Data-driven starting values for the two scales.

    Roughness is estimated as the median residual of a local plane fit at a
    range of radii; the normal scale is taken where roughness stops growing,
    which is the scale at which the surface stops looking like a plane.
    """
    rng = np.random.default_rng(0)
    pts = cloud[rng.choice(len(cloud), size=min(sample, len(cloud)), replace=False)]
    tree = cKDTree(cloud)

    out = {}
    for radius in (0.2, 0.5, 1.0, 2.0, 4.0):
        resid = []
        for p in pts[:2000]:
            idx = tree.query_ball_point(p, radius)
            if len(idx) < 10:
                continue
            nb = cloud[idx] - cloud[idx].mean(axis=0)
            normal = np.linalg.svd(nb, full_matrices=False)[2][-1]
            resid.append(np.median(np.abs(nb @ normal)))
        out[radius] = float(np.median(resid)) if resid else float("nan")

    roughness = out
    knee = max((r for r in roughness if np.isfinite(roughness[r])), default=1.0)
    return {"roughness_by_radius": roughness,
            "suggested_normal_scale": knee,
            "suggested_projection_radius": knee / 2}

When the extra cost is justified

M3C2 is slower than differencing two rasters by a wide margin, and on a flat site it will give you the same answer.

Where M3C2 earns its cost and where a raster difference is enough Two columns. The M3C2 column covers steep faces and overhangs where a vertical difference is meaningless, surfaces whose normal direction changes across the scene, clouds of differing density between epochs, and any case needing a per-point confidence interval. The raster difference column covers near-horizontal ground, stockpiles on a flat pad, uniform density across both epochs, and reporting that needs a volume rather than per-point distances. use M3C2 steep faces and overhangs normals varying across the scene epochs of differing density a per-point confidence interval is needed a raster difference suffices near-horizontal ground stockpiles on a flat pad uniform density in both epochs the deliverable is a volume M3C2 measures along the local normal, which is what makes it right where vertical is the wrong direction.

Figure 3 — The choice follows the geometry, not the budget.

Edge-case matrix

Situation Symptom Handling
Normal scale too small Noisy, sign-flipping change Increase until normals stabilise
Normal scale too large Real curvature smoothed, change biased Reduce below the feature size
Projection radius too small Few points, huge uncertainty Increase, or subsample core points less
Projection radius too large Small features vanish Match to the smallest feature of interest
Thin vegetation on the face Change dominated by leaves Filter to ground and rock classes first
Overhanging geometry Cylinder catches two surfaces Reduce max_depth below the overhang separation
Very different point densities Uncertainty asymmetric between epochs Expected; the formula already handles it
Unregistered clouds Uniform apparent change Co-register first; M3C2 cannot fix alignment

Verification snippet

import numpy as np


def verify_m3c2(result: dict, stable_core_mask: np.ndarray) -> dict:
    """Stable core points must come out insignificant and unbiased."""
    ch = result["change"][stable_core_mask]
    sig = result["significant"][stable_core_mask]
    ch = ch[np.isfinite(ch)]

    bias = float(np.median(ch)) if ch.size else float("nan")
    false_rate = float(np.count_nonzero(sig) / max(sig.size, 1))
    problems = []
    if abs(bias) > 0.03:
        problems.append(f"stable faces show {bias:+.3f} m — registration residual")
    if false_rate > 0.10:
        problems.append(f"{false_rate:.1%} of stable core points are significant — "
                        "the registration error term is understated")
    return {"stable_bias_m": bias, "stable_significant_rate": false_rate,
            "measured_points": int(ch.size), "problems": problems}
Measurement uncertainty and feature resolution against projection radius Two curves against projection radius from ten centimetres to two metres. Measurement uncertainty falls steeply from twenty centimetres at a ten centimetre radius to about three centimetres at one metre, following an inverse square root. The smallest resolvable feature rises linearly with the radius, from twenty centimetres to four metres. The curves cross at about forty centimetres, which is marked as the usable compromise for a rock face where thirty centimetre blocks matter. 0.1 m 0.25 0.5 1.0 1.5 2.0 projection radius uncertainty ∝ 1/√n smallest resolvable feature usable compromise ≈ 0.4 m The right radius is set by the smallest feature that matters, not by the noise you would like.

Figure 2 — The trade the projection radius makes, and where the answer usually sits for a rock face.

Making it fast enough to run on a whole face

The implementation above is a Python loop over core points, which is fine for a few thousand and unusable for a million. Three changes make it practical without leaving Python.

Subsample the core points aggressively. The change field is spatially smooth at the projection radius, so measuring at a spacing of roughly half that radius loses nothing. On a face sampled at 2 cm with a 40 cm projection radius, that is a twenty-fold reduction in core points before any other optimisation.

Batch the neighbour queries. cKDTree.query_ball_point accepts an array of points and a workers argument, so the tree lookups — which dominate the runtime — parallelise across cores with one call instead of one per point.

Keep the cylinder test vectorised. The radial and along-normal decomposition is pure array arithmetic; doing it once per core point over a candidate array is fast, and doing it per candidate point in Python is not.

import numpy as np
from scipy.spatial import cKDTree


def core_points(cloud: np.ndarray, spacing: float) -> np.ndarray:
    """Voxel-downsample a cloud to one representative point per cell."""
    keys = np.floor(cloud / spacing).astype(np.int64)
    _, idx = np.unique(keys, axis=0, return_index=True)
    return cloud[np.sort(idx)]


def batched_neighbours(tree: cKDTree, pts: np.ndarray, radius: float,
                       workers: int = -1) -> list:
    """All neighbour lists in one parallel call rather than one call per point."""
    return tree.query_ball_point(pts, radius, workers=workers)

With those three in place, a face of forty million points measured at a 20 cm core spacing runs in a few minutes on an ordinary workstation, which is comfortably inside the budget of a monthly monitoring job. The memory profile is modest too, because the trees hold indices rather than copies and the candidate arrays are transient.

One caution on subsampling: take the core points from the reference epoch only, and keep the same core set across every survey in a monitoring series. Re-deriving core points each month changes where the measurements are taken, which introduces a small apparent change of its own and makes the series harder to trend. Storing the core set with the site, alongside the stable polygon, removes the problem entirely.

When to escalate

  • The face is too steep for the drone to have seen twice. No distance metric recovers geometry that only one flight captured. Add oblique passes to the flight plan.
  • Significant change appears everywhere at similar magnitude. That is a registration residual, not movement, and the registration error term is understating it. Fix the alignment first.
  • Results must be reconciled with a volume. M3C2 measures distances, not volumes. Converting requires an assumption about the affected area, and it is usually better to compute the volume from surfaces and quote the distances separately.

← Change Detection Between Survey Epochs