Skip to content

Leaves on a QSM

A QSM describes the wood. sylva.leaves adds the foliage as flat polygons whose total area, position and orientation follow what was measured, giving a whole-tree model for radiative transfer or visualisation.

from sylva import leaves, qsm

wood = leaves.classify_leaf_wood(tree)                  # bool per point
model = qsm.build_qsm(qsm.wood_points(tree), base_xy=(x, y))
foliage = tree[~wood]

angles = leaves.leaf_angle_distribution(foliage)        # from point normals
print(angles.mean_deg, angles.de_wit, angles.chi, angles.g([0.0, 1.0]))

area = leaves.leaf_area_density(foliage, voxel_size=0.25)
mesh = leaves.add_leaves(model, area, angles, leaf_points=foliage,
                         leaf_length=0.08, leaf_width=0.04)
leaves.write_tree_obj("tree.obj", model, mesh)         # objects "wood" and "leaves"

The four steps

Leaf / wood labels. classify_leaf_wood uses, by default, the graph-based separation of Tian and Li (2022), following their GBSeparation code. Shortest paths from the base give every point a path length and a direction of growth; graph edges are cut where they are long for their neighbourhood or join points growing in different directions, which severs leaves from the branch they hang on before any shape is judged. The cut graph is split into shells of path length at several scales (0.1–1 m for trees under 15 m, 0.5–3 m above), and a connected piece of a shell is wood when it runs the whole shell along its growth direction and is either cylindrical (a circle fits its cross-section) or linear; a piece thicker than the wood below it on its path is rejected. Wood then spreads down every path to the base and to close neighbours. Judging a 10–100 cm segment rather than a 5 cm neighbourhood is what makes it work: over a few centimetres a leaf is as planar as bark.

method="passage" is the QSM's own wood filter (path passage plus local anisotropy, here with a second 10 cm scale). It keeps nearly all the wood and half the leaves with it, which is what a QSM wants (a cylinder fit suffers more from missing wood than from stray leaves) and the wrong trade for leaf work:

30 manually labelled tropical trees (Van den Broeck et al. 2025) accuracy mIoU wood recall / precision leaf recall / precision
anisotropy only 0.75 0.56 0.68 / 0.57 0.77 / 0.84
passage (the QSM filter) 0.61 0.44 0.97 / 0.45 0.44 / 0.97
passage, two scales 0.68 0.51 0.95 / 0.50 0.55 / 0.96
graph-based, 0.1–1 m shells 0.80 0.66 0.93 / 0.64 0.73 / 0.96
graph-based, 0.5–3 m shells (default for tall trees) 0.90 0.79 0.91 / 0.80 0.89 / 0.96

As the QSM input the graph-based labels are worse (destructive-harvest volume rRMSE 26 % against 20 %), so qsm.wood_points keeps the passage filter; wood_points(method="gbs") is there to try. return_scores=True adds a wood confidence per point, a weight for qsm.build_qsm(weights=) (see Wood labels and wood weights).

Leaf angle distribution. Normals from a PCA over each leaf point's neighbours (Vicari et al. 2019) give inclinations (angle between the leaf normal and the vertical), weighted by the area each point stands for so densely scanned leaves do not dominate. The result carries the histogram, mean, a two-parameter beta fit (Goel and Strebel 1984; goel_strebel gives it as the lad_params of lad="twoParamBeta" in sylva.voxels), Campbell's (1990) ellipsoidal χ, the nearest de Wit (1965) type and the projection function g(zenith) (Wilson 1960). LeafAngleDistribution.from_type("planophile") gives the analytic types.

Leaf area density. Without pulses, leaf_area_density counts the surface in the points: thinned to one point per cube of side res, a surface with normal n crosses (|nx| + |ny| + |nz|) / res² cubes per unit area, so each survivor stands for res² / (|nx| + |ny| + |nz|). res defaults to 3.5 × the median point spacing. This is a box count with no plateau (smaller cubes undercount, larger ones overcount at leaf edges), and it only sees foliage the scanner saw. With pulse data prefer a ray-traced grid, LeafAreaGrid.from_voxels(grid, field); with a known total (litterfall, allometry, hemispherical photos) keep the spatial pattern and rescale: area.scaled_to(total_m2). A plain number instead of a grid does that in one step: add_leaves(model, 85.0, angles, leaf_points=foliage).

Insertion. Each voxel receives leaves until its area is met, centred on leaf points of that voxel (uniformly inside it if it has none), with normals drawn from the angle distribution and a uniform azimuth; blades within max_branch_distance of a cylinder point away from it and record that cylinder. Leaves may intersect each other: no collision test is made, unlike the non-intersecting insertion of Åkerblom et al. (2018).

Leaf shape and size

The blade is a LeafShape: a mesh in unit leaf space ((along, across, up) with the base at the origin, the tip at along = 1 and the greatest width 1 across), plus the length and width it is placed at. The default is the six-vertex outline at 8 × 4 cm; shape.area is the area of one leaf, as is single_leaf_area(length, width).

Set the size per call, or once for the session:

mesh = leaves.add_leaves(model, area, angles, leaf_points=foliage,
                         leaf_length=0.15, leaf_width=0.06)

leaves.set_default_leaf(length=0.15, width=0.06)   # every later call
leaves.default_leaf()                              # what is in force

A custom blade takes any triangle mesh of a single leaf, so a scanned or modelled one (lobed, curled, or several leaflets) can be used instead. It is read with the base at the smallest x, the tip along +x, the blade across ±y and any curl in z; along then scales with the length and across and up with the width. A mesh drawn in metres keeps the size it was drawn at:

shape = leaves.LeafShape.from_obj("eucalypt_leaf.obj")   # at its own size
shape = shape.resized(length=0.12)                       # or set one
shape = shape.scaled_to(0.004)                           # or an area (m2)
mesh = leaves.add_leaves(model, area, angles, leaf_points=foliage, shape=shape)
leaves.set_default_leaf(shape)                           # or make it the default

The leaf count follows from the area to be met divided by the area of one leaf, so a bigger blade gives fewer leaves for the same leaf area, and a curled blade counts the area of its triangles, not of its outline.

How well it works

Against synthetic trees with separate wood and leaf meshes, sampled as points without occlusion (so the estimators are isolated from visibility):

result
mean leaf inclination within 0.5° of the mesh (48–59°), histogram overlap 0.95, correct de Wit type, G within 0.01, and the same from the classified leaf points as from the true ones
point-based leaf area, true leaf points 0.9–1.15 × the mesh area from 5 000 to 80 000 points per m²
leaf / wood labels (graph-based) accuracy 0.81–0.96; 96–99 % of leaf points found at 80–98 % precision; wood precision 0.86–0.97, wood recall 0.24–0.94 by point count (most wood surface in these trees is millimetre twigs inside the foliage, which end up as leaf)
leaf area after classification 0.95–1.16 × the mesh area
inserted leaves area met to within one leaf, inclinations as asked, vertical leaf-area profile r = 0.97–1.00 against the mesh

So the labels, the angle distribution and the placement hold up on clouds without occlusion. What a real scan adds is visibility: the point-based area only counts foliage the scanner saw, which is why a ray-traced or independently known total is the better input for dense crowns.