Segmenting Trees Directly from Points

TL;DR: Normalize heights with filters.hag_nn, keep vegetation above about 2 m, swap HeightAboveGround into Z with filters.ferry, then run filters.litree with min_height: 3, min_points: 50 and a radius of 50–100 m; each point gets a TreeID, and a pandas groupby turns those labels into a tree table.

# Context and Motivation

This guide is part of Individual Tree Segmentation from LiDAR. A canopy height model throws away every return below the top surface, which is most of the returns in a multi-layered forest. Point-based segmentation keeps them. filters.litree implements the region-growing method of Li, Guo, Jakubowski and Kelly (2012): it repeatedly takes the highest unassigned point as a new tree top, then walks down through the remaining points in height order, assigning each to the current tree if it is closer to that tree’s points than to the competing set, using spacing thresholds that relax lower in the crown.

The benefit is twofold. Understorey trees that grow in gaps or beneath taller crowns can be found, because their points exist even when the CHM hides them. And the output is a label on every point, so per-tree structural metrics — height percentiles, crown base height, vertical profile — come straight from the points rather than from a raster.

Growing trees from the top down A side view of two crowns. Points are processed from highest to lowest. The highest point seeds tree 1. The next high point, far from tree 1, seeds tree 2. Lower points are assigned to whichever tree's existing points are nearer, shown by arrows. A horizontal sweep line marks the current height being processed. seed tree 1 seed tree 2 current height grey points below the sweep line are still unassigned

# Prerequisites and Assumptions

  • PDAL 2.3+ with filters.litree (check with pdal --options filters.litree) and the Python bindings.
  • Ground classified to class 2 for height normalization.
  • Density of at least 10 pts/m²; the method relies on vertical structure inside crowns.
  • Tiles of a few hectares to 1 km². The algorithm’s cost grows faster than linearly with points per tree, so very dense tiles benefit from thinning.

# Step-by-Step Implementation

# Step 1 — Normalize heights and keep vegetation

Compute HeightAboveGround, drop ground and noise, and keep points above 2 m so grass and shrubs do not seed trees.

# Step 2 — Put height above ground into Z

The region growing uses Z for ordering and spacing. On sloping terrain, raw elevation would put the uphill side of a crown “above” the downhill tree next to it. Ferrying HeightAboveGround into Z removes the slope; keep a copy of the original if you need it later.

json
{ "type": "filters.ferry", "dimensions": "Z=>Elevation, HeightAboveGround=>Z" }

# Step 3 — Run filters.litree

Set min_height to the lowest tree top you want, min_points to the smallest tree worth keeping, and radius to bound the search for competing trees.

# Step 4 — Summarise per TreeID

Group points by TreeID in pandas: count, maximum height, plan extent and centroid of the highest points.

# Step 5 — Compare with the CHM result

Match tops from both methods within 2 m. Trees found only by filters.litree are candidates for understorey; inspect a few to confirm they are real.

# Complete Working Example

python
"""Point-based tree segmentation with filters.litree and a per-tree summary."""
from __future__ import annotations

import json
from pathlib import Path

import numpy as np
import pandas as pd
import pdal


def litree(src: Path, dst: Path, min_height: float = 3.0, min_points: int = 50,
           radius: float = 60.0) -> np.ndarray:
    stages = [
        {"type": "readers.las", "filename": str(src)},
        {"type": "filters.range", "limits": "Classification![7:7],Classification![18:18]"},
        {"type": "filters.hag_nn", "count": 2},
        {"type": "filters.range",
         "limits": "Classification![2:2],HeightAboveGround[2:70]"},
        {"type": "filters.ferry", "dimensions": "Z=>Elevation, HeightAboveGround=>Z"},
        {"type": "filters.litree", "min_points": min_points, "min_height": min_height,
         "radius": radius},
        {"type": "writers.las", "filename": str(dst), "minor_version": 4,
         "dataformat_id": 6, "extra_dims": "TreeID=uint32,Elevation=double"},
    ]
    p = pdal.Pipeline(json.dumps({"pipeline": stages}))
    n = p.execute()
    print(f"{n} vegetation points processed")
    return p.arrays[0]


def tree_table(points: np.ndarray) -> pd.DataFrame:
    df = pd.DataFrame({"x": points["X"], "y": points["Y"], "h": points["Z"],
                       "tree": points["TreeID"]})
    df = df[df.tree > 0]
    top = df.loc[df.groupby("tree").h.idxmax(), ["tree", "x", "y", "h"]].set_index("tree")
    stats = df.groupby("tree").agg(points=("h", "size"),
                                   h95=("h", lambda s: float(np.percentile(s, 95))),
                                   x_span=("x", np.ptp), y_span=("y", np.ptp))
    out = top.join(stats).rename(columns={"h": "height_m"})
    out["crown_diam_m"] = (out.x_span + out.y_span) / 2
    return out.drop(columns=["x_span", "y_span"]).round(2)


if __name__ == "__main__":
    pts = litree(Path("stand_12.laz"), Path("out/stand_12/litree.laz"))
    trees = tree_table(pts)
    print(f"{len(trees)} trees; median height {trees.height_m.median():.1f} m")
    print(trees.sort_values("height_m", ascending=False).head(10))

The label value 0 marks points not assigned to any tree — usually low vegetation or isolated returns below min_height that never reached a seeded tree.

# Key Parameter Table

Option Type Default Guidance
min_points int 10 Minimum points for a tree; 30–80 at 20 pts/m² removes branch fragments
min_height float, m 3.0 Lowest tree top considered; matches the inventory’s height threshold
radius float, m 100.0 Search radius for competing trees; 40–60 m is ample and faster
lower height cut float, m 2.0 Points below never join trees; prevents shrubs seeding
input thinning voxel, m none 0.2–0.3 m voxels speed dense tiles with little effect on results
Where each method finds its trees Grouped bars comparing trees detected by a CHM watershed and by filters.litree on a plot with 140 field-measured trees. For dominant trees both find about 70. For intermediate trees the watershed finds 28 and litree 34. For suppressed understorey trees the watershed finds 4 and litree 17. A note says litree's advantage is almost entirely below the main canopy. dominant intermediate understorey CHM watershed filters.litree illustrative plot of 140 trees — the difference sits almost entirely below the main canopy

# Verification

  • Seed count sanity. The number of trees should be in the range of stand records. Many more usually means min_points is too low and branch clusters are becoming trees.
  • Height agreement. For trees found by both methods, heights should agree within about half a metre; the point method is often slightly higher because it reads the actual top return.
  • Visual check by TreeID. Colour points by TreeID in a viewer with a random palette. Crowns should be contiguous blobs; a tree whose points appear in two separate places indicates a spacing problem.
python
labelled = int((pts["TreeID"] > 0).sum())
share = labelled / len(pts)
assert 0.5 < share <= 1.0, f"only {share:.0%} of vegetation points assigned to trees"

# Gotchas and Edge Cases

Slope without normalization. Skipping the ferry step makes uphill crowns dominate their downhill neighbours on any real slope. It is the most common cause of implausibly large trees on hillsides.

Slope turns equal trees into giants and dwarfs Left: two 18 metre trees on a 20 degree slope drawn in raw elevation; the uphill tree's top is 14 metres higher than the downhill tree's, so region growing treats the downhill crown as lower branches of the uphill tree. Right: the same trees after height normalization stand on a flat baseline with equal tops, and are segmented as two trees. raw Z: merged into one height above ground: two

Very dense data is slow. Run time grows quickly with points per crown. Thinning with filters.voxelcenternearestneighbor at 0.25 m before segmentation typically cuts time several-fold; then transfer labels back to full-density points with a nearest-neighbour join if you need every point labelled.

Leaning and multi-stemmed trees. Top-down growth assumes each tree has one highest point above its crown. Leaning trees whose top overhangs a neighbour, and coppiced trees with several leaders, get split or merged. No setting fixes that; report it as a known limitation for those stand types.

Tile edges. Like every tree method, crowns cut by the tile edge are wrong. Use a buffer and keep trees whose top lies inside the nominal tile.

# Frequently Asked Questions

Is filters.litree better than a CHM watershed?

It finds more trees below the main canopy and gives per-point labels, at the cost of much longer run times and fewer tuning options. For dominant trees in open stands both methods perform similarly; for multi-layered forests the point method has a real advantage.

Why must height above ground replace Z before filters.litree?

The algorithm orders and compares points by Z. On a slope, raw elevation mixes terrain height into that comparison, so trees uphill appear taller than their neighbours downhill. Normalized heights compare trees as they actually stand.

How do I speed up filters.litree on dense data?

Restrict input to vegetation above 2 m, reduce radius to 40 to 60 metres, and voxel-thin to about 0.25 metres. Labels can be transferred back to the full-density cloud afterwards.

What does TreeID 0 mean?

The point was not assigned to any tree, usually because it lies below min_height and never connected to a seeded tree, or belongs to a group smaller than min_points.