How much does pulse density really matter for Li 2012 segmentation?

More than any parameter. At 4 returns per square metre the mean spacing is about half a metre, so a 1.5 m threshold spans only three points and crowns fragment. At 15 returns per square metre the same threshold spans dozens and the rule behaves as intended. Below roughly 6 returns per square metre a raster method will beat it on the same data.

Point Cloud Tree Segmentation with the Li 2012 Algorithm in Python

Q: Can I get a proper crown polygon instead of a convex hull?

Yes. An alpha shape over each tree's points follows concave crown boundaries far better, at the cost of an extra dependency and one more parameter. For area-based allometry the convex hull is usually adequate because the bias is small and consistent; for crown-shape classification the alpha shape is worth the trouble.

A canopy height model is a projection, and projection discards the understory. When the inventory question is about vertical structure — suppressed stems beneath a closed overstory, multi-layered tropical forest, a two-storey shelterwood — the segmentation has to happen in three dimensions. This guide implements the Li 2012 top-down distance rule against a normalized point cloud, as the 3-D alternative inside Individual Tree Detection and Crown Segmentation, part of Canopy Height Modeling & Terrain Extraction.

When to use this approach

Point-cloud segmentation costs one to two orders of magnitude more compute than a raster watershed and needs far more returns to work. It buys exactly one thing: trees with no visible apex.

Situation	Raster watershed	Li 2012 in the point cloud
Even-aged conifer plantation	Correct choice — fast and accurate	Wasted compute
Multi-layered broadleaf or tropical forest	Misses the whole understory	The reason this algorithm exists
Pulse density below 6–8 returns/m²	Stable	Unreliable — the distance rule needs points
Regional survey, thousands of tiles	Minutes per tile	Hours per tile
Crown polygons required as output	Direct	Needs a convex hull or alpha shape per tree

The honest rule: if the CHM already reaches an F₁ above about 0.85 on your validation plots, stay on the raster. Move to the point cloud when recall is capped by structure rather than by tuning.

How the distance rule works

Li 2012 is a greedy top-down pass with a single geometric decision. Points are sorted by height, descending. The tallest unassigned point becomes the apex of a new tree. Every remaining point is then tested against the trees already open: if its horizontal distance to the current tree’s points is less than a threshold, it joins that tree; otherwise it starts or joins another.

The threshold is not constant, because a tall tree’s crown is wider than a sapling’s. A two-level rule covers most stands:

d_{t} (z) = {d_{low} d_{high} z < z_{break} z \geq z_{break}

with $d_{low} \approx 1.5$ m below a break height of about 15 m and $d_{high} \approx 2.0$ m above it. The comparison that matters is between the distance to the current tree and the distance to the nearest other tree — a point closer to a neighbouring crown than to the one being grown is left for that neighbour.

Minimal reproducible example

The implementation below is deliberately self-contained — no LiDAR-specific dependency beyond laspy for reading — so the rule is visible rather than hidden behind a library call. It uses a KD-tree over each tree’s assigned points so the distance test stays sub-linear.

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


def li2012(xyz: np.ndarray, dt_low: float = 1.5, dt_high: float = 2.0,
           z_break: float = 15.0, min_points: int = 12) -> np.ndarray:
    """Li 2012 top-down segmentation. xyz is (n,3) normalized height-above-ground.

    Returns an int array of tree ids (0 = unassigned).
    """
    order = np.argsort(-xyz[:, 2])          # tallest first
    labels = np.zeros(len(xyz), dtype="int32")
    trees: list[list[int]] = []             # index lists, one per tree

    for i in order:
        if labels[i]:
            continue
        p = xyz[i, :2]
        dt = dt_high if xyz[i, 2] >= z_break else dt_low

        best, best_d = -1, np.inf
        for t, members in enumerate(trees):
            pts = xyz[members][:, :2]
            d = float(np.min(np.linalg.norm(pts - p, axis=1)))
            if d < best_d:
                best, best_d = t, d

        if best >= 0 and best_d <= dt:
            trees[best].append(int(i))
            labels[i] = best + 1
        else:
            trees.append([int(i)])
            labels[i] = len(trees)

    # Drop trees too sparse to be real, renumbering the survivors from 1.
    keep = [t for t, m in enumerate(trees) if len(m) >= min_points]
    remap = {old + 1: new + 1 for new, old in enumerate(keep)}
    return np.array([remap.get(v, 0) for v in labels], dtype="int32")

The inner loop over trees is the honest but slow formulation. For production, keep one cKDTree per open tree and query it, rebuilding lazily as members accumulate:

def li2012_fast(xyz, dt_low=1.5, dt_high=2.0, z_break=15.0, min_points=12):
    """Same rule, with a KD-tree per open tree and a spatial pre-filter."""
    order = np.argsort(-xyz[:, 2])
    labels = np.zeros(len(xyz), dtype="int32")
    seeds: list[np.ndarray] = []            # apex xy per tree, for the pre-filter
    members: list[list[int]] = []
    trees_kd: list[cKDTree | None] = []

    for i in order:
        p = xyz[i, :2]
        dt = dt_high if xyz[i, 2] >= z_break else dt_low
        best, best_d = -1, np.inf
        if seeds:
            apex = np.vstack(seeds)
            # Only trees whose apex is within 4*dt can possibly win.
            near = np.where(np.linalg.norm(apex - p, axis=1) <= 4 * dt)[0]
            for t in near:
                if trees_kd[t] is None:
                    trees_kd[t] = cKDTree(xyz[members[t]][:, :2])
                d, _ = trees_kd[t].query(p)
                if d < best_d:
                    best, best_d = int(t), float(d)
        if best >= 0 and best_d <= dt:
            members[best].append(int(i))
            trees_kd[best] = None            # invalidate; rebuilt on next query
            labels[i] = best + 1
        else:
            seeds.append(p)
            members.append([int(i)])
            trees_kd.append(None)
            labels[i] = len(members)

    keep = [t for t, m in enumerate(members) if len(m) >= min_points]
    remap = {old + 1: new + 1 for new, old in enumerate(keep)}
    return np.array([remap.get(v, 0) for v in labels], dtype="int32")

Parameter reference

Parameter	Type	Default	Recommended range	Ecological rationale
`dt_low`	float (m)	1.5	1.0–2.0	Distance threshold below the break height; roughly the crown radius of an understory stem
`dt_high`	float (m)	2.0	1.5–3.0	Threshold in the overstory; scale it to dominant crown radius
`z_break`	float (m)	15.0	10–20	Height at which crowns start behaving like overstory; set near the base of the dominant layer
`min_points`	int	12	8–30	Rejects clusters too sparse to be a stem; raise it as pulse density rises
pulse density	returns/m²	≥ 8	8–40	Below 6 the distance rule fragments crowns; above 40 it is compute-bound, so thin first

Thin very dense clouds before segmenting. Decimating a 60 returns/m² cloud to about 20 changes the result marginally and cuts runtime by roughly an order of magnitude.

Derive per-tree metrics and crowns

Once every point carries a tree id, the metrics that matter come straight out of the groups. A convex hull gives a crown polygon; an alpha shape is closer to reality but rarely worth the dependency.

import geopandas as gpd
from shapely.geometry import MultiPoint


def tree_metrics(xyz: np.ndarray, labels: np.ndarray, crs) -> gpd.GeoDataFrame:
    """Per-tree height, point count, and convex-hull crown from a segmented cloud."""
    rows = []
    for tid in np.unique(labels):
        if tid == 0:
            continue
        pts = xyz[labels == tid]
        hull = MultiPoint([(x, y) for x, y in pts[:, :2]]).convex_hull
        rows.append({
            "tree_id": int(tid),
            "height_m": float(pts[:, 2].max()),
            "n_points": int(len(pts)),
            "crown_area_m2": float(hull.area),
            "base_height_m": float(np.percentile(pts[:, 2], 5)),
            "geometry": hull,
        })
    return gpd.GeoDataFrame(rows, crs=crs)

The base_height_m column is the payoff for working in 3-D: it is the bottom of the live crown, which is exactly the quantity fuel modelling needs and which no height raster can supply. It feeds directly into Deriving Canopy Base Height from LiDAR Vertical Profiles.

Expected output and verification

A correct run assigns the large majority of vegetation points to a tree, produces a height distribution with a visible understory mode, and yields plausible crown areas.

import laspy

las = laspy.read("normalized_tile.laz")
veg = las.classification >= 3            # ASPRS low/med/high vegetation
xyz = np.vstack([las.x[veg], las.y[veg], las.HeightAboveGround[veg]]).T
xyz = xyz[xyz[:, 2] >= 2.0]

labels = li2012_fast(xyz)
assigned = (labels > 0).mean()
n_trees = int(labels.max())

print(f"{n_trees} trees, {assigned:.1%} of vegetation points assigned")
assert assigned > 0.85, "too many orphan points — dt is too small or the cloud is too sparse"
assert n_trees > 0, "no trees segmented — check the height normalization"

An assignment rate below about 85% means the thresholds are too tight for the point spacing. A rate near 100% with implausibly few trees means the opposite: everything merged into a handful of giant clusters, which shows up immediately as crown areas in the hundreds of square metres.

Common pitfalls

Running on absolute elevations. The rule compares heights, so the cloud must already be normalized to height above ground. Feeding raw Z produces one enormous tree per hilltop.
Including ground and noise classes. Ground returns sit below every crown and, processed last, join whichever tree is nearest, inflating crown hulls. Filter to vegetation classes and above the minimum height first.
Leaving the break height at a published default. z_break should sit near the base of your dominant canopy layer. In a 12 m coppice stand a 15 m break means every point uses the understory threshold.
Segmenting a whole tile in one pass. Runtime grows super-linearly with open trees. Chunk into 50–100 m blocks with a 5 m overlap, segment each, and reconcile by apex position.

Frequently Asked Questions

How much does pulse density really matter?

More than any parameter. At 4 returns/m² the mean spacing is about 0.5 m, so a 1.5 m threshold spans only three points and a crown fragments into pieces. At 15 returns/m² the same threshold spans dozens and the rule behaves as intended. Below roughly 6 returns/m² a raster method will beat this one on the same data.

Can I get a proper crown polygon instead of a convex hull?

Yes — an alpha shape over each tree’s points follows concave crown boundaries far better, at the cost of an extra dependency and one more parameter to choose. For area-based allometry the convex hull is usually adequate because the bias is small and consistent; for crown-shape classification it is not, and the alpha shape is worth the trouble.

Why does the algorithm sometimes split one tree into a trunk and a crown?

Because the horizontal distance from a stem return to the nearest crown return can exceed the threshold in a tree with a long clear bole. Raising dt_low fixes it but merges genuinely separate understory stems. The cleaner fix is to run the assignment on crown points only, then attach stem returns to whichever tree owns the column above them.

Individual Tree Detection and Crown Segmentation — the raster-first workflow this is the 3-D alternative to
Normalizing LiDAR Point Clouds with PDAL — produces the height-above-ground field this algorithm requires
Computing Vertical Canopy Complexity Metrics from LiDAR — the plot-level view of the same vertical structure
Deriving Canopy Base Height from LiDAR Vertical Profiles — consumes the per-tree base height this produces

Up: Individual Tree Detection & Segmentation · Canopy Height Modeling & Terrain Extraction