Python LiDAR

DBSCAN Segmentation with filters.dbscan

TL;DR: filters.dbscan grows groups only through “core” points that have at least min_points neighbours within eps, and labels everything unreachable as -1. Pick min_points around 2 × the number of dimensions (6–12 for 3D), read eps off the elbow of a k-distance curve, and cluster on X,Y,HeightAboveGround when terrain slope would otherwise stretch objects apart.

# Context and Motivation

This guide is part of Point Cloud Segmentation with PDAL. DBSCAN — density-based spatial clustering of applications with noise — solves the two problems Euclidean clustering cannot. First, a thin bridge of a few points no longer merges two objects, because those bridge points do not have enough neighbours to be core points and cannot pass the group on. Second, isolated clutter is labelled noise rather than becoming thousands of tiny groups.

The cost is a second parameter and a little more thought about density. eps is still a distance; min_points turns it into a density threshold. Together they say: a region is part of an object if it contains at least min_points points within any eps-radius ball.

Core, border, noise A small point set with eps circles drawn around three example points. A core point has five neighbours inside its circle, meeting min_points. A border point has only two neighbours but lies inside a core point's circle, so it joins that group without extending it. A noise point has one neighbour and lies in no core point's circle, so it is labelled -1. core: 5 neighbours ≥ min_points border: joins, cannot extend noise: labelled −1 circles show the eps radius; min_points = 5 here

# Prerequisites and Assumptions

  • PDAL 2.3+ with filters.dbscan and the Python bindings.
  • SciPy for the k-distance curve.
  • Points restricted to candidates (ground and noise removed), in a projected CRS with metres in every clustered dimension.

# Step-by-Step Implementation

# Step 1 — Choose min_points

A common rule is twice the number of dimensions: 6 for X,Y,Z. Raise it to 10–20 on dense data or when you want thin structures treated as noise; lower it to 4 for sparse wires.

# Step 2 — Compute the k-distance curve

For each point, find the distance to its k-th nearest neighbour with k = min_points. Sort those distances. The curve is flat for points inside objects and rises sharply for points on the fringes and in clutter.

# Step 3 — Read eps at the elbow

The elbow — where the sorted curve bends upward — is the natural eps: small enough that clutter fails the density test, large enough that object interiors pass it.

# Step 4 — Pick the clustering space

Use dimensions: "X,Y,Z" for objects separated vertically. Use "X,Y,HeightAboveGround" on slopes, where terrain would otherwise stretch objects.

# Step 5 — Run and account for labels

Count -1 points separately. A noise share of a few percent is normal; above 20 percent usually means eps is too small or min_points too high.

# Complete Working Example

python
"""Choose eps from a k-distance curve, then run filters.dbscan."""
from __future__ import annotations

import json
from pathlib import Path

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


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


def elbow_eps(points: np.ndarray, k: int, dims: tuple[str, ...], sample: int = 150_000) -> float:
    xyz = np.column_stack([points[d] for d in dims])
    if len(xyz) > sample:
        xyz = xyz[np.random.default_rng(1).choice(len(xyz), sample, replace=False)]
    d, _ = cKDTree(xyz).query(xyz, k=k + 1)
    kd = np.sort(d[:, -1])
    # Elbow: point of maximum distance from the chord joining the curve's ends.
    x = np.linspace(0.0, 1.0, len(kd))
    y = (kd - kd[0]) / (kd[-1] - kd[0])
    elbow = int(np.argmax(x - y))
    return float(kd[elbow])


def dbscan(points: np.ndarray, eps: float, min_points: int, dims: str) -> np.ndarray:
    p = pdal.Pipeline(json.dumps({"pipeline": [
        {"type": "filters.dbscan", "eps": round(eps, 2), "min_points": min_points,
         "dimensions": dims}]}), arrays=[points])
    p.execute()
    return p.arrays[0]


if __name__ == "__main__":
    pts = candidates(Path("hillside_22.laz"))
    k = 8
    dims = ("X", "Y", "HeightAboveGround")
    eps = elbow_eps(pts, k, dims)
    out = dbscan(pts, eps, k, ",".join(dims))
    ids = out["ClusterID"]
    print(f"eps={eps:.2f} m  groups={len(np.unique(ids[ids >= 0]))}  "
          f"noise={np.mean(ids < 0):.1%}")
Reading eps off the k-distance curve A curve of the distance to the eighth nearest neighbour, sorted from smallest to largest over all points. It stays nearly flat around 0.4 to 0.6 metres for the first ninety percent of points, then bends sharply upward to over 3 metres for the last few percent. The elbow at about 0.8 metres is marked as the chosen eps, with a dashed chord drawn from the curve's start to its end. elbow: eps ≈ 0.8 m object interiors: flat fringes and clutter points sorted by k-distance 8th-NN distance

# Interpreting the Output

The run prints three numbers worth reading together. The chosen eps should sit within a small multiple of the point spacing you expect from the flight — if it comes out at 3 m on 20 pts/m² data, the elbow detection latched onto something odd, usually a sparse region such as water or a swath edge dominating the sample. The group count should be in the range of objects you expect in the tile. And the noise share tells you how much the density test threw away.

It helps to look at noise in space, not only as a percentage. Write the result with ClusterID as an extra dimension and colour noise points distinctly. Noise concentrated on crown fringes and roof edges is healthy — those are exactly the sparse regions DBSCAN is meant to reject. Noise covering whole objects means those objects were sampled more sparsely than the rest of the tile, and a single eps cannot serve both; that is the cue to thin the dense areas or segment them separately.

Finally, remember that noise is not deleted. The points are still in the array with label -1, and a later step can reassign them — for example, giving each noise point the label of its nearest grouped neighbour within a short distance, which recovers crown fringes after the groups themselves have been found cleanly.

# Key Parameter Table

Option Type Default Guidance
eps float 1.0 Elbow of the k-distance curve; roughly 2–4 × median spacing
min_points int 6 2 × dimensions as a floor; higher treats thin structures as noise
dimensions string X,Y,Z Any dimensions in consistent units; X,Y,HeightAboveGround on slopes
where expression none Run on a subset while keeping all points (unprocessed get the default label)

# Verification

  • Noise share. Report the fraction of -1 points; compare across tiles for consistency.
  • Thin-bridge test. Find a known case — a hedge touching a house, a branch over a roof — and confirm the two objects receive different IDs.
  • Group count stability. As with Euclidean clustering, rerun at ±20 percent eps. DBSCAN is usually more stable than Euclidean clustering around its chosen setting.
python
ids = out["ClusterID"]
noise = float(np.mean(ids < 0))
assert noise < 0.2, f"{noise:.0%} noise: eps too small or min_points too high"

# Gotchas and Edge Cases

Varying density defeats one eps. DBSCAN assumes one density defines objects. Near nadir and in overlap zones density doubles; at swath edges it halves. A single eps then over-segments sparse areas or merges dense ones. Normalizing density with voxel thinning first, or running per flightline, helps.

Border points are order-dependent. A border point within eps of core points in two groups joins whichever group reaches it first. The groups themselves are stable; a few border assignments can change between runs or versions.

Mixed units in dimensions. Adding Intensity or ReturnNumber to dimensions mixes metres with counts. Rescale such dimensions before clustering, or leave them out.

One eps, two densities Two identical shrubs. The left one, near nadir in an overlap zone, is densely sampled and forms one group. The right one, at the swath edge, is sampled at half the density; with the same eps, many of its points fail the core test and are labelled noise, splitting the shrub into a small group and scattered noise. overlap zone: one group swath edge: fragments every point has ≥ min_points within eps grey: labelled −1 at the same eps

# Frequently Asked Questions

How do I choose eps for filters.dbscan?

Compute the distance from each point to its k-th nearest neighbour, where k equals min_points, sort those distances and pick the value at the elbow where the curve bends upward. That separates object interiors from fringes and clutter.

What does ClusterID -1 mean?

Noise: the point is not a core point and is not within eps of any core point. It belongs to no group. Treat it as unassigned, and track its share as a diagnostic.

When should I use DBSCAN instead of filters.cluster?

When objects touch through thin connections, or when the input contains scattered clutter you want discarded. On clean, well-separated objects Euclidean clustering is simpler and faster.

Can DBSCAN cluster on attributes as well as coordinates?

Yes, through the dimensions option, but every dimension must be on a comparable scale because eps is a single distance. Rescale non-spatial dimensions first.

Related in Point Cloud Segmentation with PDAL

↑ Point Cloud Segmentation with PDAL overview --- title: "Point Cloud Segmentation with PDAL" description: "Group points into objects with filters.cluster and filters.dbscan: how... Extracting Objects from Segment Labels --- title: "Extracting Objects from Segment Labels" description: "From ClusterID to usable objects: build a per-object table with pandas,... Euclidean Segmentation with filters.cluster --- title: "Euclidean Segmentation with filters.cluster" description: "Group points into connected objects with PDAL filters.cluster:...