Extracting Objects from Segment Labels

TL;DR: Load the labelled points into pandas, groupby("ClusterID") to build one row per object with count, extent, height statistics and median features, filter objects with plain boolean rules, assign a class per object by rule or majority vote and map it back to points with Series.map, then export object points with a PDAL filters.groupby pipeline and object footprints with GeoPandas.

# Context and Motivation

This guide is part of Point Cloud Segmentation with PDAL. Segmentation stages stop at a label per point. Every workflow in this section then needs the same second half: describe each labelled group, decide which groups are real objects, give each object a class, and deliver it in a form someone can use. Doing that consistently — one table, one set of rules, one mapping back to points — is what separates a reproducible extraction from a notebook full of one-off masks.

The pattern also fixes the most visible flaw of per-point classification: speckle. When the class is decided per object and written to every point of that object, a roof cannot contain a handful of “vegetation” points, because the decision was never made per point.

# Prerequisites and Assumptions

• Points with a segment label — ClusterID from filters.cluster or filters.dbscan, or TreeID from filters.litree.
• Any per-point features you want to aggregate, such as HeightAboveGround, Planarity, Scattering.
• Python with pandas, NumPy, Shapely 2.x and GeoPandas; PDAL bindings.

# Step-by-Step Implementation

# Step 1 — Build the object table

Put the dimensions you need into a DataFrame and aggregate by label with named aggregations, so every column has a clear meaning.

# Step 2 — Filter objects

Apply size, extent and height rules as boolean expressions on the table. Keep the rules in one place and log how many objects each removes.

# Step 3 — Decide a class per object

Rules on medians are the simplest; a majority vote of per-point model predictions is the alternative when a classifier exists.

# Step 4 — Map the decision back to points

Series.map from label to class gives a per-point class column in one vectorized step; points of dropped objects keep their original class.

# Step 5 — Export

Write classified points to LAS, each object to its own LAZ with filters.groupby if the client wants object files, and footprints plus attributes to a GeoPackage.

# Complete Working Example

"""Object table, per-object classes, and exports from segment labels."""
from __future__ import annotations

import json
from pathlib import Path

import geopandas as gpd
import numpy as np
import pandas as pd
import pdal
import shapely
from shapely.geometry import MultiPoint


def object_table(a: np.ndarray, label: str = "ClusterID") -> pd.DataFrame:
    df = pd.DataFrame({"x": a["X"], "y": a["Y"], "hag": a["HeightAboveGround"],
                       "planarity": a["Planarity"], "scattering": a["Scattering"],
                       "label": a[label]})
    df = df[df.label > 0]
    return df.groupby("label").agg(
        n=("x", "size"),
        cx=("x", "mean"), cy=("y", "mean"),
        dx=("x", np.ptp), dy=("y", np.ptp),
        h_max=("hag", "max"), h_p50=("hag", "median"),
        planarity=("planarity", "median"), scattering=("scattering", "median"),
    )


def decide(obj: pd.DataFrame) -> pd.Series:
    cls = pd.Series(0, index=obj.index, dtype="uint8")        # 0 = leave unchanged
    big = obj.n >= 50
    cls[big & (obj.planarity > 0.7) & (obj.scattering < 0.1) & (obj.h_max > 2.5)] = 6
    cls[big & (obj.scattering >= 0.2) & (obj.h_max > 5.0)] = 5
    cls[big & (obj.scattering >= 0.2) & obj.h_max.between(2.0, 5.0)] = 4
    print(cls.value_counts().rename({0: "unchanged", 4: "medium veg", 5: "high veg", 6: "building"}))
    return cls


def apply(a: np.ndarray, cls: pd.Series, label: str = "ClusterID") -> np.ndarray:
    out = a.copy()
    per_point = pd.Series(out[label]).map(cls).fillna(0).to_numpy().astype("uint8")
    change = per_point > 0
    out["Classification"][change] = per_point[change]
    return out


def footprints(a: np.ndarray, obj: pd.DataFrame, cls: pd.Series, crs: str) -> gpd.GeoDataFrame:
    rows = []
    for lab in cls.index[cls == 6]:
        sel = a["ClusterID"] == lab
        hull = shapely.concave_hull(MultiPoint(np.column_stack([a["X"][sel], a["Y"][sel]])),
                                    ratio=0.3)
        rows.append({"object": int(lab), "height_m": round(float(obj.at[lab, "h_max"]), 2),
                     "points": int(obj.at[lab, "n"]), "geometry": hull})
    return gpd.GeoDataFrame(rows, crs=crs)


def per_object_files(src_with_labels: Path, out_dir: Path) -> None:
    out_dir.mkdir(parents=True, exist_ok=True)
    pdal.Pipeline(json.dumps({"pipeline": [
        str(src_with_labels),
        {"type": "filters.range", "limits": "Classification[6:6]"},
        {"type": "filters.groupby", "dimension": "ClusterID"},
        {"type": "writers.las", "filename": str(out_dir / "building_#.laz")},
    ]})).execute()


if __name__ == "__main__":
    p = pdal.Pipeline(json.dumps({"pipeline": ["block_07_features.laz"]}))
    p.execute()
    arr = p.arrays[0]
    objects = object_table(arr)
    classes = decide(objects)
    classified = apply(arr, classes)
    pdal.Writer.las(filename="block_07_classified.laz", minor_version=4, dataformat_id=6,
                    forward="all", extra_dims="ClusterID=int64").pipeline(classified).execute()
    footprints(classified, objects, classes, "EPSG:6347").to_file(
        "block_07_objects.gpkg", layer="buildings", driver="GPKG")
    per_object_files(Path("block_07_classified.laz"), Path("out/buildings"))

In writers.las, the # in the filename is replaced by a running number for each group that filters.groupby produces, which is how one pipeline writes one file per object.

# Key Parameter Table

# Verification

• Point conservation. The classified array has the same length as the input, and only points in decided objects changed class.
• Table and points agree. For a sample of objects, the class in the table equals the class of every point with that label.
• Rules removed what you expect. Log how many objects each rule accepted, and look at the rejected objects closest to each threshold.

before, after = arr["Classification"], classified["Classification"]
changed = before != after
labels_changed = np.unique(arr["ClusterID"][changed])
assert set(labels_changed) <= set(classes.index[classes > 0]), "points outside decided objects changed"

# Gotchas and Edge Cases

Label 0 and -1. Euclidean clustering uses 0 for unassigned points and DBSCAN uses -1 for noise. Exclude both before aggregating, or they become one giant “object” that dominates every statistic.

Labels are per tile. The same number means different objects in different tiles. Build a global key such as f"{tile}_{label}" before merging tables across tiles.

Mapping with missing keys. Series.map returns NaN for labels not in the decision table; the fillna(0) is what keeps those points unchanged. Forgetting it turns them into class 0 after the cast.

Objects on tile edges. A building cut by the tile edge becomes two partial objects in two tiles. Process buffered tiles and keep objects whose centroid lies inside the nominal tile.

# Frequently Asked Questions

Why decide classes per object instead of per point?

Because many properties that identify an object, such as its footprint area, height range or median planarity, only exist at object level, and because one decision per object removes the scattered misclassified points that per-point decisions leave behind.

How do I write one LAS file per object?

Use filters.groupby on the label dimension followed by writers.las with a # in the filename. Each group is written to its own file with the placeholder replaced by a running number.

What should happen to objects that no rule accepts?

Leave their points’ classification unchanged, usually class 1. Record them in the object table with a reason so you can review the rejected objects nearest the thresholds.

How do I carry object IDs into the delivered LAS?

Write the label as an extra dimension with writers.las extra_dims, for example ClusterID as int64. Readers such as PDAL and laspy return it by name.

# Related

• Point Cloud Segmentation — producing the labels
• Euclidean Segmentation with filters.cluster — connected-component labels
• DBSCAN Segmentation with filters.dbscan — density-based labels and noise
• Extracting Building Footprints from LiDAR — the dedicated footprint route
• Training a Random Forest Point Classifier — per-point predictions to vote over