M3C2 Change Detection with py4dgeo

This page runs M3C2 (Multiscale Model to Model Cloud Comparison) between two LAZ epochs with the py4dgeo library — building a regular core-point grid, choosing the normal and cylinder radii from the data rather than by habit, feeding in the measured registration error, and reading the per-point level of detection and sample counts that make the result defensible.

Why you hit this

A stockpile report, a subsidence alert or a claim that a contractor over-excavated all need a distance with an uncertainty attached, and M3C2 is the standard method that provides one. It is also a method with four parameters that interact, and the defaults in most tutorials are tuned for terrestrial scans of cliffs at hundreds of points per square metre. Run with those settings on 20 pts/m² airborne data and nearly every cylinder is under-sampled: distances come back as NaN across half the site, or as confident numbers computed from three points. The broader workflow this sits inside, from registration to change polygons, is in change detection between LiDAR scan epochs.

Prerequisites

  • py4dgeo>=0.7 (pip install py4dgeo), laspy[lazrs]>=2.5, numpy>=1.24.
  • Two epochs in the same CRS, EPSG:32618+5703 in the examples, registered on stable ground, with vegetation and noise classes removed.
  • The stable-ground registration residual from the alignment step, in metres — typically 0.01–0.03 m between two airborne surveys.
  • A rough idea of the point density of each epoch; pdal info --metadata reports it, or see computing point density from LAZ with PDAL.

Step-by-Step

1. Load the epochs and build a regular core-point grid

python
import numpy as np
import py4dgeo

epoch0, epoch1 = py4dgeo.read_from_las("quarry_2025-06_ground.laz", "quarry_2026-06_ground.laz")
print(f"t0 {epoch0.cloud.shape[0]:,} pts, t1 {epoch1.cloud.shape[0]:,} pts")

def grid_corepoints(cloud, spacing=1.0):
    cells = np.floor(cloud[:, :2] / spacing).astype(np.int64)
    _, first = np.unique(cells, axis=0, return_index=True)
    return cloud[np.sort(first)]

corepoints = grid_corepoints(epoch0.cloud, spacing=1.0)
print(f"{len(corepoints):,} core points at 1 m spacing")

A grid of core points, one per square metre of the reference epoch, is better than every n-th point for two reasons. Its results can be rasterised and summed per cell without a density weighting, which matters for volumes. And it spends computation evenly across the site instead of concentrating it where overlapping flight lines tripled the density.

2. Choose the normal radius from surface roughness

The normal radius D sets the scale at which the surface orientation is estimated. Too small and the normal follows every pebble; too large and it smooths over a real edge.

python
from scipy.spatial import cKDTree

def roughness_at(cloud, sample, radius):
    tree = cKDTree(cloud)
    out = []
    for p in sample:
        nbrs = cloud[tree.query_ball_point(p, radius)]
        if len(nbrs) < 10:
            continue
        centred = nbrs - nbrs.mean(axis=0)
        out.append(np.sqrt(np.linalg.eigvalsh(centred.T @ centred / len(nbrs))[0]))
    return np.median(out)

sample = corepoints[np.random.default_rng(7).choice(len(corepoints), 2000, replace=False)]
for r in (0.25, 0.5, 1.0, 2.0):
    print(f"radius {r:4.2f} m → roughness σ ≈ {roughness_at(epoch0.cloud, sample, r) * 1000:.1f} mm")

The square root of the smallest eigenvalue is the standard deviation of the points about their best-fit plane — the local roughness. The guidance from the method’s authors is that the normal scale should be roughly twenty to twenty-five times the roughness, so that surface noise does not tilt the normal. For a quarry floor with 20 mm roughness that points to a normal radius around 0.5 m; for a rubble slope with 80 mm, closer to 2 m. Passing several radii lets py4dgeo pick, per core point, the scale at which the surface is most planar.

The two radii that define an M3C2 measurement At a core point, the normal is estimated from all points within the normal radius. A cylinder of the projection radius is then extended along that normal through both epochs up to the maximum distance. The distance is the separation of the two epochs' mean positions along the normal, and the spread and count of points in each epoch's slice of the cylinder give the level of detection. epoch t₀ epoch t₁ cylinder radius d ≥ 20–30 points per epoch normal radius D ≈ 20–25 × roughness distance max_distance caps the cylinder length
D decides which way the cylinder points; d decides how many points each epoch contributes to the mean. They are tuned to different properties of the data.

3. Choose the cylinder radius from point density

The cylinder has to contain enough points in each epoch for the mean position and its spread to be meaningful.

python
density_t0 = 18.0      # pts/m², from pdal info on the ground-classified tile
density_t1 = 31.0
target = 30
cyl_radius = np.sqrt(target / (np.pi * min(density_t0, density_t1)))
print(f"cylinder radius for ≥{target} pts in the sparser epoch: {cyl_radius:.2f} m")

For 18 pts/m² that gives about 0.73 m. The number is set by the sparser epoch, because the level of detection is dominated by whichever mean is less certain. Round up rather than down: a cylinder with 25 points gives a slightly smoother distance map; one with 8 gives a level of detection that is itself noisy.

4. Run M3C2 with the registration error

python
m3c2 = py4dgeo.M3C2(
    epochs=(epoch0, epoch1),
    corepoints=corepoints,
    normal_radii=(0.5, 1.0, 2.0),
    cyl_radius=0.75,
    max_distance=15.0,
    registration_error=0.018,
)
distances, uncertainties = m3c2.run()

lod = uncertainties["lodetection"]
n0, n1 = uncertainties["num_samples1"], uncertainties["num_samples2"]
valid = np.isfinite(distances)
print(f"valid {valid.mean() * 100:.1f}% | median LoD {np.nanmedian(lod) * 1000:.0f} mm | "
      f"median samples t0 {np.median(n0[valid]):.0f}, t1 {np.median(n1[valid]):.0f}")

max_distance bounds how far along the normal the cylinder looks for the second epoch. Set it just above the largest change you expect — 15 m covers a year of quarrying — because a longer cylinder on a steep face can pass through an unrelated surface and return a plausible, wrong distance. registration_error is added to the level of detection for every core point, so a well-registered pair gains sensitivity everywhere.

5. Classify each core point

python
MIN_SAMPLES = 10
reliable = valid & (n0 >= MIN_SAMPLES) & (n1 >= MIN_SAMPLES)
significant = reliable & (np.abs(distances) > lod)

status = np.full(len(corepoints), "unknown", dtype=object)
status[reliable & ~significant] = "no detectable change"
status[significant & (distances > 0)] = "gain"
status[significant & (distances < 0)] = "loss"

for s in ("gain", "loss", "no detectable change", "unknown"):
    print(f"{s:>22}: {(status == s).mean() * 100:5.1f}%")

Three outcomes are not enough; there have to be four. A core point with no second-epoch points in its cylinder, or too few, is not “no change” — it is unmeasured, and reporting it as stable is how an occluded area behind a new building ends up certified as unchanged.

Four outcomes for every core point A decision tree. If either epoch has too few samples in the cylinder, or the distance is not finite, the core point is unknown. Otherwise, if the absolute distance is at most the level of detection, there is no detectable change. If it exceeds the level of detection, the sign of the distance decides between gain and loss. samples ≥ 10 in both epochs? unknown |distance| > LoD? no change change gain loss no yes
"Unknown" is a result. Folding it into "no change" turns every occlusion shadow into a false certificate of stability.

6. Compute volume change per cell

Because core points sit on a 1 m grid, each significant distance on near-horizontal ground represents one square metre of change.

python
horizontal = np.abs(m3c2.directions()[:, 2]) > 0.9            # normal within ~25° of vertical
cell_area = 1.0
gain_m3 = (distances[significant & horizontal & (distances > 0)]).sum() * cell_area
loss_m3 = -(distances[significant & horizontal & (distances < 0)]).sum() * cell_area
lod_vol = np.sqrt((lod[significant & horizontal] ** 2).sum()) * cell_area
print(f"gain {gain_m3:,.0f} m³, loss {loss_m3:,.0f} m³, ± {lod_vol:,.0f} m³ (quadrature LoD)")

Restricting volume to near-horizontal normals avoids counting a quarry face twice — once as horizontal retreat along its normal and again in the floor below it. Summing levels of detection in quadrature assumes independent cells, which is optimistic for registration error shared by the whole site; for a conservative bound, add registration_error × area linearly.

Expected Output & Verification

text
t0 7,204,611 pts, t1 12,388,092 pts
181,442 core points at 1 m spacing
cylinder radius for ≥30 pts in the sparser epoch: 0.73 m
valid 96.4% | median LoD 61 mm | median samples t0 29, t1 52
                  gain:   3.8%
                  loss:  21.6%
  no detectable change:  70.9%
               unknown:   3.7%
gain 1,420 m³, loss 38,915 m³, ± 312 m³ (quadrature LoD)

Verify against something the pipeline did not produce. Swap the epochs and rerun: gain and loss should exchange, and each core point’s distance should change sign within its LoD. Then compare the loss against the site’s extraction records — weighbridge tonnage divided by the material’s bulk density — and expect agreement within a few percent; a larger gap usually traces back to a stockpile outside the surveyed extent.

Effect of cylinder radius on level of detection and resolution As the cylinder radius grows from 0.25 to 2 metres on an 18 points per square metre survey, the median number of samples per cylinder rises and the level of detection falls quickly at first and then flattens, while the smallest resolvable feature grows linearly. The useful range is where the level of detection has flattened but features are still small, around 0.6 to 1 metre. 0.250.61.0 1.52.0 m level of detection smallest resolvable feature useful range cylinder radius, for an 18 pts/m² airborne survey
Beyond about 30 samples per cylinder, a larger radius buys almost no sensitivity and keeps costing spatial resolution.

Common Errors

Most distances are NaN. The cylinder radius is too small for the density, or max_distance is shorter than the change. Check num_samples2 on a known-changed area: zero samples with a large expected change means max_distance; a handful of samples everywhere means the radius.

The level of detection is implausibly small, a few millimetres. registration_error was left at its default of zero, so only the roughness term is counted. Two independent surveys never agree to a few millimetres; pass the measured stable-ground residual.

Distances have the wrong sign on some slopes. Normals are oriented by default towards positive z, and on overhanging or near-vertical faces that flips between neighbouring core points. For cliffs and walls, orient normals towards the scanner position or a known viewpoint using py4dgeo’s orientation options for your version before comparing signs.

Frequently Asked Questions

How is this different from CloudCompare’s M3C2 plugin?

Same algorithm and the same parameters, run from Python, which is what a pipeline needs: reproducible settings in version control, tiles processed in a batch job and results written straight into the next step. Use CloudCompare to explore a pair of epochs interactively and to sanity-check parameters.

Should core points come from the earlier or later epoch?

Conventionally the earlier one, so distances describe what happened to the reference surface. For a newly built area with no earlier points, use a regular XY grid lifted to the later epoch’s surface; otherwise the new structure has no core points at all.

Can M3C2 run on a whole city at once?

Not in one call on one machine. Tile both epochs on the same grid with an overlap at least the size of the largest radius, run tiles in parallel and keep only each tile’s interior core points.

Back to Change Detection Between LiDAR Scan Epochs.