Benchmarking SMRF Against Reference Ground Points

TL;DR: Interpolate the classified ground surface at each control point, difference it against the surveyed elevation, and report RMSE, mean bias and the 95th percentile separately for open, sloped and vegetated terrain — a single aggregate number hides exactly the failure you are looking for.

# Context and Motivation

This guide is part of SMRF Ground Classification in PDAL. Tuning without measurement is guesswork with confidence, and the measurement is not difficult — it just has to be set up so it can fail.

The subtlety is that ground classification has two independent error modes and one aggregate statistic cannot express both. Omission is real ground rejected, which shows up as voids and, where interpolation fills them, as a surface pulled toward whatever survived. Commission is objects accepted, which shows up as a surface pulled up. A parameter change that halves one usually increases the other, so an RMSE that improved tells you nothing about which trade you made. Reporting them separately, and by terrain class, is what makes the number actionable.

# Prerequisites and Assumptions

The datum row is where most benchmarks go wrong before they start. Control in orthometric heights against a cloud in ellipsoidal heights produces a uniform thirty-metre bias that looks like a catastrophic classification failure.

# Step-by-Step Implementation

# Step 1 — Classify and keep ground only

{"type": "filters.range", "limits": "Classification[2:2]"}

# Step 2 — Interpolate the surface at each control location

Take the mean elevation of ground returns within a small radius — one metre is typical — of each control point. A control point with no ground returns nearby is an omission failure and must be counted as such rather than dropped.

# Step 3 — Difference and split by class

Signed differences, so the mean is a bias and not another magnitude.

# Step 4 — Report bias, RMSE and the tail, per class

The 95th percentile of absolute error is what a client notices; the mean bias is what a systematic problem looks like.

# Complete Working Example

"""Benchmark a classified ground surface against surveyed control points."""
from __future__ import annotations

import json
from pathlib import Path

import numpy as np
import pdal


def ground_points(src: Path) -> np.ndarray:
    p = pdal.Pipeline(json.dumps({"pipeline": [
        {"type": "readers.las", "filename": str(src)},
        {"type": "filters.range", "limits": "Classification[2:2]"},
    ]}))
    p.execute()
    return p.arrays[0]


def surface_at(ground: np.ndarray, x: float, y: float, radius: float = 1.0) -> float | None:
    dx = ground["X"] - x
    dy = ground["Y"] - y
    near = (dx * dx + dy * dy) <= radius * radius
    if not near.any():
        return None
    return float(ground["Z"][near].mean())


def benchmark(src: Path, control: list[dict], radius: float = 1.0) -> dict:
    ground = ground_points(src)
    per_class: dict[str, list[float]] = {}
    missing: dict[str, int] = {}

    for pt in control:
        cls = pt["terrain"]
        z = surface_at(ground, pt["x"], pt["y"], radius)
        if z is None:
            missing[cls] = missing.get(cls, 0) + 1
            continue
        per_class.setdefault(cls, []).append(z - pt["z"])

    out = {}
    for cls, diffs in sorted(per_class.items()):
        d = np.asarray(diffs)
        out[cls] = {
            "n": int(d.size),
            "no_ground_nearby": missing.get(cls, 0),
            "bias_m": round(float(d.mean()), 3),
            "rmse_m": round(float(np.sqrt((d ** 2).mean())), 3),
            "p95_abs_m": round(float(np.percentile(np.abs(d), 95)), 3),
        }
    all_diffs = np.concatenate([np.asarray(v) for v in per_class.values()])
    out["aggregate"] = {
        "n": int(all_diffs.size),
        "bias_m": round(float(all_diffs.mean()), 3),
        "rmse_m": round(float(np.sqrt((all_diffs ** 2).mean())), 3),
    }
    return out


if __name__ == "__main__":
    control = json.loads(Path("control.json").read_text())
    print(json.dumps(benchmark(Path("classified.laz"), control), indent=2))

# Key Parameter Table

# Verification

No uniform bias. A mean of −30 m across every class is a datum problem, not a classifier problem.

Missing counts are reported. Silently dropping control points where classification found no ground turns the worst failure into a better-looking number.

The comparison is fair. Both parameter sets benchmarked against the same control on the same tile, with only the parameters changed.

# Gotchas and Edge Cases

Control points on hard surfaces are easy. A benchmark made only of road-centreline points reports a flattering RMSE that says nothing about vegetated terrain.

Interpolation radius interacts with slope. On steep ground a one-metre radius spans real elevation change, which appears as error. Shrink it or model the local slope.

A better RMSE can be a worse surface. Halving commission by rejecting more ground improves nothing if it strips ridges — which is precisely what the per-class breakdown reveals.

# Frequently Asked Questions

Why report errors by terrain class?

Because ground classification has two opposing error modes and one aggregate number cannot express both. Two parameter sets can produce almost identical overall RMSE while one is twice as accurate on open ground and far worse on slopes. The breakdown is what makes the comparison actionable.

What should I do with control points where no ground was found?

Count them. A control point with no ground return nearby is the clearest possible omission signal, and dropping it from the statistics converts the worst failure mode into a better-looking RMSE.

My benchmark shows a uniform 30 metre bias — what is wrong?

The datums, almost certainly. Control in orthometric heights compared against a cloud in ellipsoidal heights produces exactly this, and it looks like a catastrophic classification failure. Check the vertical CRS on both sides before touching any parameter.

How many control points do I need?

At least thirty per terrain class for the RMSE to mean anything. A benchmark with ten points in a class reports noise, and a benchmark made entirely of road-centreline points reports a flattering number that says nothing about vegetated ground.

# Related

• SMRF Ground Classification — the classifier being measured
• Tuning SMRF for Forested Terrain — the parameter changes this benchmark evaluates
• SMRF vs PMF for Dense Urban LiDAR — the same comparison between two algorithms
• Setting a Vertical CRS on a Point Cloud — the datum agreement every benchmark depends on
• Ground Filtering and DTM/DSM Generation with PDAL — the section overview