Point Cloud Segmentation with PDAL

Almost every object-level product built from LiDAR — a building footprint, a tree, a pole, a parked car, a power-line span — starts with the same move: take a set of points and decide which of them belong together. That is segmentation, and in PDAL it comes down to two stages. filters.cluster performs Euclidean connected-component labelling: any two points closer than a tolerance end up in the same group. filters.dbscan performs density-based clustering: points join a group only if they sit in a dense enough neighbourhood, and isolated points are labelled noise. Both write a ClusterID dimension. Choosing between them, and choosing their two or three parameters, decides whether the objects that come out are the objects you meant. This topic is the shared foundation for the rest of the classification and feature extraction section.

# Prerequisites

• PDAL 2.3+; filters.dbscan arrived in the 2.x series and filters.cluster has been available much longer. Check both with pdal --options filters.dbscan.
• Python 3.10+ with NumPy, pandas and SciPy for spacing estimates and segment statistics.
• A projected CRS in metres. Tolerances and eps are distances; in a geographic CRS they would be degrees and meaningless.
• Points already restricted to the objects you want. Segmentation is only as good as its input: remove ground with a class filter and cut to a height band with height above ground before grouping.
• A sense of point spacing. Both methods are parameterized by distance; the median nearest-neighbour distance of your input is the number every setting is measured against.

# Core Workflow Architecture

1. Restrict. Filter to the points that could belong to the objects of interest — a class, a height band, a feature threshold.
2. Measure spacing. Compute the median and 90th-percentile nearest-neighbour distance of the restricted points.
3. Choose the method. Euclidean when objects are well separated and you want every point labelled; DBSCAN when objects touch through thin bridges or the input carries isolated clutter.
4. Set distances from spacing. tolerance or eps at two to three times the median spacing; min_points from the smallest object you care about.
5. Segment. Run the stage and read back ClusterID, noting that DBSCAN uses -1 for noise and Euclidean clustering leaves unclustered points at 0 when a where clause excludes them.
6. Summarise and filter. Aggregate each segment’s size, extent and features in pandas, and drop segments that fail size or shape tests.

# Full Implementation

"""Segment points with Euclidean or DBSCAN clustering and summarise each segment."""
from __future__ import annotations

import json
import logging
from pathlib import Path
from typing import Literal

import numpy as np
import pandas as pd
import pdal
from scipy.spatial import cKDTree

log = logging.getLogger("segment")


def restricted(src: Path, where: str, hag_band: tuple[float, float]) -> np.ndarray:
    lo, hi = hag_band
    p = pdal.Pipeline(json.dumps({"pipeline": [
        {"type": "readers.las", "filename": str(src)},
        {"type": "filters.hag_nn", "count": 2},
        {"type": "filters.range", "limits": f"HeightAboveGround[{lo}:{hi}]"},
        {"type": "filters.expression", "expression": where},
    ]}))
    p.execute()
    return p.arrays[0]


def spacing(points: np.ndarray, sample: int = 200_000) -> tuple[float, float]:
    xyz = np.column_stack([points["X"], points["Y"], points["Z"]])
    if len(xyz) > sample:
        xyz = xyz[np.random.default_rng(0).choice(len(xyz), sample, replace=False)]
    d, _ = cKDTree(xyz).query(xyz, k=2)
    return float(np.median(d[:, 1])), float(np.percentile(d[:, 1], 90))


def segment(points: np.ndarray, method: Literal["euclidean", "dbscan"],
            min_points: int, factor: float = 2.5) -> np.ndarray:
    median, p90 = spacing(points)
    dist = round(max(median * factor, p90), 2)
    if method == "euclidean":
        stage = {"type": "filters.cluster", "tolerance": dist,
                 "min_points": min_points, "is3d": True}
    else:
        stage = {"type": "filters.dbscan", "eps": dist,
                 "min_points": min_points, "dimensions": "X,Y,Z"}
    log.info("median spacing %.2f m, p90 %.2f m -> %s distance %.2f m",
             median, p90, method, dist)
    p = pdal.Pipeline(json.dumps({"pipeline": [stage]}), arrays=[points])
    p.execute()
    return p.arrays[0]


def summarise(points: np.ndarray) -> pd.DataFrame:
    df = pd.DataFrame({k: points[k] for k in ("X", "Y", "Z", "HeightAboveGround", "ClusterID")})
    noise = int((df.ClusterID < 0).sum())
    unassigned = int((df.ClusterID == 0).sum())
    df = df[df.ClusterID > 0]
    seg = df.groupby("ClusterID").agg(
        n=("X", "size"),
        x=("X", "mean"), y=("Y", "mean"),
        dx=("X", lambda s: s.max() - s.min()),
        dy=("Y", lambda s: s.max() - s.min()),
        height=("HeightAboveGround", "max"),
    )
    log.info("%d segments, %d noise points, %d unassigned", len(seg), noise, unassigned)
    return seg


if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s")
    pts = restricted(Path("block_07.laz"), "Classification == 1", (1.0, 60.0))
    labelled = segment(pts, "dbscan", min_points=10)
    table = summarise(labelled)
    print(table.sort_values("n", ascending=False).head(20))

# Code Breakdown

Restriction happens in PDAL, segmentation on an array. The restricted points are passed back into a second pipeline with arrays=[points]. That split lets you measure spacing on exactly the points that will be segmented, then choose distances from the measurement, without reading the file twice. The mechanics are covered in passing NumPy arrays into a PDAL pipeline.

Spacing on a sample. A k-d tree query over two hundred thousand points is fast and gives a stable median; querying every point of a dense tile would take longer than the segmentation itself.

max(median * factor, p90). The median alone underestimates gaps in sparse parts of an object — the far side of a roof from the scanner, the underside of a crown. Taking at least the 90th percentile of spacing keeps those parts connected.

is3d: true. Euclidean clustering in 3D keeps a bridge deck separate from the road beneath it and a crown separate from a car under it. Set it to false only when you deliberately want plan-view grouping, for example to merge all returns of a building regardless of storey.

dimensions: "X,Y,Z" for DBSCAN. The stage can cluster on any dimensions, which is powerful — clustering on X,Y,HeightAboveGround removes the effect of terrain slope — but it means distances are only meaningful if every listed dimension is in the same units.

Three kinds of label. DBSCAN writes -1 for noise. Euclidean clustering drops groups smaller than min_points back to label 0, and points excluded by a where clause also stay 0. The summary counts them separately because a sudden change in either is the first sign of badly chosen parameters.

# Choosing the Space You Segment In

Distance is only meaningful in the space where you measure it, and the default — raw X,Y,Z — is not always the right one. filters.dbscan makes the choice explicit through its dimensions option, and thinking about it saves a lot of threshold tuning.

Raw coordinates are right when objects are separated in three dimensions: a bridge deck above a road, a crown above a car, stacked pipes in a plant. Here height matters, and removing it would merge things that are genuinely apart.

Plan position plus height above ground — X,Y,HeightAboveGround — is right on slopes. A hedge running down a hillside spans ten metres of elevation but only one metre of height above ground; in raw coordinates its points are far apart vertically and the hedge fragments, while in normalized space it is a compact band. The same applies to terraced houses stepping down a street.

Plan position only — X,Y, or is3d: false on filters.cluster — is right when you want a footprint, not an object: all returns of a building including balconies and eaves, or all returns under one crown including the understorey.

Scaled feature spaces are possible but need care. Clustering on X,Y,Z,Intensity mixes metres with intensity counts, so a difference of 1 in intensity weighs the same as a metre of distance. If you want to include a non-spatial dimension, rescale it first with a Python filter so that one unit means something comparable to the spatial tolerance.

# Parameter Reference Table

# Validation and Integrity Checks

• Label accounting. Segmented plus noise plus unassigned must equal the input count. Segmentation stages do not drop points, so any mismatch means a stage between reading and summarising did.
• Size distribution. Plot a histogram of segment sizes on a log scale. A healthy result has a long tail of small segments and a clear population of object-sized ones. A single enormous segment holding most of the points means the distance is far too large.
• Spot checks by colour. Write the labelled points with ClusterID as an extra dimension and colour by it in a viewer. Two minutes looking at a handful of blocks catches merges that no statistic will.
• Stability across settings. Rerun at 0.8 and 1.2 times the chosen distance. If the number of object-sized segments changes by more than about ten percent, you are on the slope of the curve above, not on the plateau.

def label_accounting(before: int, labelled: np.ndarray) -> None:
    ids = labelled["ClusterID"]
    parts = int((ids > 0).sum()) + int((ids < 0).sum()) + int((ids == 0).sum())
    assert parts == before == len(labelled), "segmentation changed the point count"

# Performance Tuning

Both stages build a k-d tree and then do one radius search per point, so run time grows with point count and with the number of neighbours inside the distance. Two consequences follow. First, restricting the input is by far the biggest saving — segmenting only elevated, non-ground points typically handles a third of the tile. Second, an overly large eps or tolerance is slow as well as wrong, because every search returns hundreds of neighbours.

Neither stage is multithreaded, so parallelism comes from running tiles in separate processes; see parallel tile processing. Use buffered tiles and keep only segments whose centroid lies in the unbuffered area, so objects on edges are counted once.

# Common Errors and Troubleshooting

One segment contains most of the tile. Ground or low vegetation was not removed, and it connects everything. Check the restriction step; a single filters.range on HeightAboveGround usually fixes it.

Thousands of one-point segments. min_points is left at its default of 1 for filters.cluster. Set it to the smallest meaningful object size.

ClusterID is -1 everywhere after DBSCAN. min_points is higher than the number of neighbours within eps anywhere in the data. Either the data is sparser than you think or eps is in the wrong units; measure spacing first.

Segments differ from run to run. DBSCAN assigns border points to whichever core point reaches them first, which can depend on order. The objects are the same, but a few border points move. If exact reproducibility matters, sort the input with filters.sort on GpsTime before segmenting.

Buildings on a slope segment in pieces. Terrain slope stretches objects vertically. Cluster on X,Y,HeightAboveGround with DBSCAN instead of raw Z, which flattens the terrain out of the distance calculation.

# Frequently Asked Questions

What is the difference between filters.cluster and filters.dbscan?

filters.cluster links any two points within the tolerance, so every point ends up in some group and thin bridges merge objects. filters.dbscan only grows groups through dense neighbourhoods and labels sparse points as noise, which separates objects joined by thin connections and discards clutter.

How do I choose eps or tolerance?

Measure the median nearest-neighbour distance of the points you will segment and start at two to three times that value. Then check that the number of object-sized segments is stable when you change the distance by twenty percent in either direction.

Why are some points labelled 0 and others -1?

DBSCAN labels noise as -1. filters.cluster uses 0 for points in groups smaller than min_points or larger than max_points, and any stage leaves points excluded by a where clause at 0. Treat both as unassigned, but count them separately when diagnosing.

When should I set max_points on filters.cluster?

Set it whenever there is a known upper bound on object size, such as the largest building in the project at your point density. Groups above the limit are relabelled 0 instead of silently becoming one merged super-object, which turns a hidden failure into a visible count you can alert on in a batch run.

Can I segment in 2D and 3D in the same pipeline?

Yes. Run filters.cluster once with is3d set to false to get footprint groups, copy ClusterID to another dimension with filters.ferry, then run it again in 3D. Each point then carries both labels, which is useful when a building’s storeys must be grouped but its separate towers must not.

Can segmentation run in streaming mode?

No. Both stages need a spatial index over all input points, so they force standard execution. Restrict the input first and use tiles small enough to fit in memory.

# Related

• Euclidean Segmentation with filters.cluster — tolerance, size limits and 2D versus 3D
• DBSCAN Segmentation with filters.dbscan — eps, min_points and the noise label
• Extracting Objects from Segment Labels — from ClusterID to per-object tables and files
• Building Extraction from LiDAR — segmentation applied to roofs
• Power Line Detection — DBSCAN applied to gappy wires