Detailed Voxels for Vegetation Modelling

Notebook creator: Hannah Weiser, 2026

This demo demonstrates the different detailedVoxel options for modelling vegetation using a single tree scene. We perform both virtual UAV-borne laser scanning (ULS) and airborne laser scanning (ALS).

[1]:
import helios
import numpy as np
import matplotlib.pyplot as plt

Creating the virtual scene

Our virtual scene consists of a ground plane and five versions of the same tree object. The stem of the tree is always loaded as 3D polygon mesh models (OBJ). The tree crown is loaded using detailedVoxels with different modes and parameters, from left to right:

  • Fixed voxel size (25 cm)

  • Fixed voxel size (5 cm)

  • Transmittive voxels (25 cm)

  • Scaled voxels by plant area density (25 cm)

  • Scaled voxels by plant area density (25 cm) with a random shift of scaled voxels within their voxel space. We use the VOX format by AMAPVox, a software for voxel-based computation of plant area density (PAD).

For reference, see Weiser et al. 2021 and the Scene page in the Wiki.

The scene parts are furthermore geometrically transformed to position them in the scene.

[2]:
# load objs and create transformations
groundplane = helios.ScenePart.from_obj(
    "../data/sceneparts/basic/groundplane/groundplane.obj"
).scale(
    15
)  # scale groundplane by 15
tree_positions = [
    [float(x), 15.0, -241.401] for x in range(-22, 3, 5)
]  # trees positioned 5 m apart along the x-axis

stem1 = helios.ScenePart.from_obj(
    "../data/sceneparts/syssifoss/F_BR08_08_stem_Poisson.obj"
)
crown1 = helios.ScenePart.from_vox(
    "../data/sceneparts/syssifoss/F_BR08_08_crown_250.vox", intersection_mode="fixed"
)

stem2 = helios.ScenePart.from_obj(
    "../data/sceneparts/syssifoss/F_BR08_08_stem_Poisson.obj"
)
crown2 = helios.ScenePart.from_vox(
    "../data/sceneparts/syssifoss/F_BR08_08_crown_50.vox", intersection_mode="fixed"
)

stem3 = helios.ScenePart.from_obj(
    "../data/sceneparts/syssifoss/F_BR08_08_stem_Poisson.obj"
)
crown3 = helios.ScenePart.from_vox(
    "../data/sceneparts/syssifoss/F_BR08_08_merged.vox",
    intersection_mode="transmittive",
)  # default

stem4 = helios.ScenePart.from_obj(
    "../data/sceneparts/syssifoss/F_BR08_08_stem_Poisson.obj"
)
crown4 = helios.ScenePart.from_vox(
    "../data/sceneparts/syssifoss/F_BR08_08_merged.vox",
    intersection_mode="scaled",
    intersection_argument=0.5,
)

stem5 = helios.ScenePart.from_obj(
    "../data/sceneparts/syssifoss/F_BR08_08_stem_Poisson.obj"
)
crown5 = helios.ScenePart.from_vox(
    "../data/sceneparts/syssifoss/F_BR08_08_merged.vox",
    intersection_mode="scaled",
    intersection_argument=0.5,
    random_shift=True,
)

# Position the trees
for i, (stem, crown) in enumerate(
    zip([stem1, stem2, stem3, stem4, stem5], [crown1, crown2, crown3, crown4, crown5])
):
    stem.translate(tree_positions[i])
    crown.translate(tree_positions[i])

# create scene
scene = helios.StaticScene(
    scene_parts=[
        groundplane,
        stem1,
        crown1,
        stem2,
        crown2,
        stem3,
        crown3,
        stem4,
        crown4,
        stem5,
        crown5,
    ]
)

Platform and scanner

[3]:
# scanner and platform for ALS and ULS surveys respectively
scanner_als = helios.scanner_from_name("riegl_lms_q780")
platform_als = helios.platform_from_name("sr22")
scanner_uls = helios.scanner_from_name("riegl_vux_1uav")
platform_uls = helios.platform_from_name("copter_linearpath")

Scanner, platform and full waveform settings

[4]:
scanner_settings_als = helios.ScannerSettings(
    pulse_frequency=100_000 * helios.units.Hz,
    scan_frequency=200 * helios.units.Hz,
    scan_angle=30 * helios.units.deg,
    trajectory_time_interval=0.05 * helios.units.s,
)

scanner_settings_uls = helios.ScannerSettings(
    pulse_frequency=100_000 * helios.units.Hz,
    scan_frequency=50 * helios.units.Hz,
    scan_angle=90 * helios.units.deg,
    trajectory_time_interval=0.01 * helios.units.s,
)

fullwave_settings = helios.FullWaveformSettings(
    beam_sample_quality=3,
    bin_size=0.25 * helios.units.ns,
    win_size=1 * helios.units.ns,
    max_fullwave_range=100 * helios.units.ns,
)

platform_settings_uls = helios.DynamicPlatformSettings(z=50, speed_m_s=5)
platform_settings_als = helios.DynamicPlatformSettings(z=500, speed_m_s=50)

Survey Route

[5]:
survey_uls = helios.Survey(
    scanner=scanner_uls,
    platform=platform_uls,
    scene=scene,
    full_waveform_settings=fullwave_settings,
)
survey_als = helios.Survey(
    scanner=scanner_als,
    platform=platform_als,
    scene=scene,
    full_waveform_settings=fullwave_settings,
)
[6]:
waypoints_uls = [
    [-30, -30, True],
    [-30, 30, False],
    [0, 30, True],
    [0, -30, False],
    [30, -30, True],
    [30, 30, True],
    [-30, 30, False],
    [-30, 0, True],
    [30, 0, False],
    [30, -30, True],
    [-30, -30, False],
]
for x, y, active in waypoints_uls:
    survey_uls.add_leg(
        x=x,
        y=y,
        is_active=active,
        platform_settings=platform_settings_uls,
        scanner_settings=scanner_settings_uls,
    )
[7]:
waypoints_als = [
    [-145, 50, True],
    [-145, -50, False],
    [145, -50, True],
    [145, 50, False],
]
for x, y, active in waypoints_als:
    survey_als.add_leg(
        x=x,
        y=y,
        is_active=active,
        platform_settings=platform_settings_als,
        scanner_settings=scanner_settings_als,
    )

Runninng the survey

[8]:
points_uls, trajectory_uls = survey_uls.run(
    verbosity=helios.LogVerbosity.VERBOSE, format=helios.OutputFormat.NPY
)
CRS bounding box (by vertices): Min: dvec3(-15.000000, -15.000000, 0.000000), Max: dvec3(15.000000, 15.000000, 28.228257)
Shift: dvec3(0.000000, 0.000000, 14.114129)
# vertices to translate: 427052
Actual bounding box (by vertices): Min: dvec3(-15.000000, -15.000000, -14.114129), Max: dvec3(15.000000, 15.000000, 14.114129)
Building KD-Grove...
KDTree (num. primitives 361488) :
        Max. # primitives in leaf: 45
        Min. # primitives in leaf: 1
        Max. depth reached: 37
        KDTree axis-aligned surface area: 5187.31
        Interior nodes: 1405364
        Leaf nodes: 1381307
        Total tree cost: 9.59737
KDGrove stats:
        Number of trees: 1
        Number of static trees: 1
        Number of dynamic trees: 0
        Statistics (min, max, total, mean, stdev):
                Building time: (1.8110, 1.8110, 1.8110, 1.8110, 0.0000)
                Tree primitives: (361488, 361488, 361488, 361488.0000, 0.0000)
                Max primitives in leaf: (45, 45, 45, 45.0000, 0.0000)
                Min primitives in leaf: (1, 1, 1, 1.0000, 0.0000)
                Maximum depth: (37, 37, 37, 37.0000, 0.0000)
                Axis-aligned surface area: (5187.3053, 5187.3053, 5187.3053, 5187.3053, 0.0000)
                Number of interior nodes: (1405364, 1405364, 1405364, 1405364.0000, 0.0000)
                Number of leaf nodes: (1381307, 1381307, 1381307, 1381307.0000, 0.0000)
                Tree cost: (9.5974, 9.5974, 9.5974, 9.5974, 0.0000)

KDG built in 1.814s
Reading Spectral Library...
10 materials found
Warning: material None of primitive 8Triangle (/home/runner/work/helios/helios/example_notebooks/../data/sceneparts/basic/groundplane/groundplane.mtl) has no spectral definition
Number of subsampling rays (riegl_vux-1uav): 19
Simulation: Scanner changed!
Pulse frequency set to 100000
Scan angle set to 90
Applying settings for PolygonMirrorBeamDeflector...
Vertical angle min/max nan/nan degrees
 -- verticalAngleMin not set, using the value of -90 degrees
 -- verticalAngleMax not set, using the value of 90 degrees
SOURCE Leg with serial ID:0 waypoints:
        Origin: (-30, -30, 35.8859)
        Target: (-30, 30, 35.8859)
        Next: (0, 30, 35.8859)

Starting simulation loop 1 ...
Waypoint reached!
Pulse frequency set to 100000
Scan angle set to 90
Applying settings for PolygonMirrorBeamDeflector...
Vertical angle min/max nan/nan degrees
 -- verticalAngleMin not set, using the value of -90 degrees
 -- verticalAngleMax not set, using the value of 90 degrees
It was not possible to determine attitude with a single computation at MovingPlatform::initLegManual
        angle = 3.14159 but it should be below 0.025
        Using iterative computation instead
Iterative mode was used for manual leg initialization because default one failed for MovingPlatform
SOURCE Leg with serial ID:1 waypoints:
        Origin: (-30, 30, 35.8859)
        Target: (0, 30, 35.8859)
        Next: (0, -30, 35.8859)

Waypoint reached!
Pulse frequency set to 100000
Scan angle set to 90
Applying settings for PolygonMirrorBeamDeflector...
Vertical angle min/max nan/nan degrees
 -- verticalAngleMin not set, using the value of -90 degrees
 -- verticalAngleMax not set, using the value of 90 degrees
It was not possible to determine attitude with a single computation at MovingPlatform::initLegManual
        angle = 3.14159 but it should be below 0.025
        Using iterative computation instead
Iterative mode was used for manual leg initialization because default one failed for MovingPlatform
SOURCE Leg with serial ID:2 waypoints:
        Origin: (0, 30, 35.8859)
        Target: (0, -30, 35.8859)
        Next: (30, -30, 35.8859)

Waypoint reached!
Pulse frequency set to 100000
Scan angle set to 90
Applying settings for PolygonMirrorBeamDeflector...
Vertical angle min/max nan/nan degrees
 -- verticalAngleMin not set, using the value of -90 degrees
 -- verticalAngleMax not set, using the value of 90 degrees
SOURCE Leg with serial ID:3 waypoints:
        Origin: (0, -30, 35.8859)
        Target: (30, -30, 35.8859)
        Next: (30, 30, 35.8859)

Waypoint reached!
Pulse frequency set to 100000
Scan angle set to 90
Applying settings for PolygonMirrorBeamDeflector...
Vertical angle min/max nan/nan degrees
 -- verticalAngleMin not set, using the value of -90 degrees
 -- verticalAngleMax not set, using the value of 90 degrees
SOURCE Leg with serial ID:4 waypoints:
        Origin: (30, -30, 35.8859)
        Target: (30, 30, 35.8859)
        Next: (-30, 30, 35.8859)

Waypoint reached!
Pulse frequency set to 100000
Scan angle set to 90
Applying settings for PolygonMirrorBeamDeflector...
Vertical angle min/max nan/nan degrees
 -- verticalAngleMin not set, using the value of -90 degrees
 -- verticalAngleMax not set, using the value of 90 degrees
SOURCE Leg with serial ID:5 waypoints:
        Origin: (30, 30, 35.8859)
        Target: (-30, 30, 35.8859)
        Next: (-30, 0, 35.8859)

Waypoint reached!
Pulse frequency set to 100000
Scan angle set to 90
Applying settings for PolygonMirrorBeamDeflector...
Vertical angle min/max nan/nan degrees
 -- verticalAngleMin not set, using the value of -90 degrees
 -- verticalAngleMax not set, using the value of 90 degrees
SOURCE Leg with serial ID:6 waypoints:
        Origin: (-30, 30, 35.8859)
        Target: (-30, 0, 35.8859)
        Next: (30, 0, 35.8859)

Waypoint reached!
Pulse frequency set to 100000
Scan angle set to 90
Applying settings for PolygonMirrorBeamDeflector...
Vertical angle min/max nan/nan degrees
 -- verticalAngleMin not set, using the value of -90 degrees
 -- verticalAngleMax not set, using the value of 90 degrees
SOURCE Leg with serial ID:7 waypoints:
        Origin: (-30, 0, 35.8859)
        Target: (30, 0, 35.8859)
        Next: (30, -30, 35.8859)

Waypoint reached!
Pulse frequency set to 100000
Scan angle set to 90
Applying settings for PolygonMirrorBeamDeflector...
Vertical angle min/max nan/nan degrees
 -- verticalAngleMin not set, using the value of -90 degrees
 -- verticalAngleMax not set, using the value of 90 degrees
It was not possible to determine attitude with a single computation at MovingPlatform::initLegManual
        angle = 3.14159 but it should be below 0.025
        Using iterative computation instead
Iterative mode was used for manual leg initialization because default one failed for MovingPlatform
SOURCE Leg with serial ID:8 waypoints:
        Origin: (30, 0, 35.8859)
        Target: (30, -30, 35.8859)
        Next: (-30, -30, 35.8859)

Waypoint reached!
Pulse frequency set to 100000
Scan angle set to 90
Applying settings for PolygonMirrorBeamDeflector...
Vertical angle min/max nan/nan degrees
 -- verticalAngleMin not set, using the value of -90 degrees
 -- verticalAngleMax not set, using the value of 90 degrees
It was not possible to determine attitude with a single computation at MovingPlatform::initLegManual
        angle = 3.14158 but it should be below 0.025
        Using iterative computation instead
Iterative mode was used for manual leg initialization because default one failed for MovingPlatform
SOURCE Leg with serial ID:9 waypoints:
        Origin: (30, -30, 35.8859)
        Target: (-30, -30, 35.8859)
        Next: (-30, -30, 35.8859)

Waypoint reached!
Pulse frequency set to 100000
Scan angle set to 90
Applying settings for PolygonMirrorBeamDeflector...
Vertical angle min/max nan/nan degrees
 -- verticalAngleMin not set, using the value of -90 degrees
 -- verticalAngleMax not set, using the value of 90 degrees
Waypoint reached!
Finishing simulation loop 1 ...
Finished simulation loop 1.
Elapsed simulation steps = 9600011
Elapsed virtual time = 96.0001 sec.
Main thread simulation loop finished in 16.9173 sec.
Waiting for completion of pulse computation tasks...
Pulse computation tasks finished in 16.9173 sec.
[9]:
points_als, trajectory_als = survey_als.run(
    verbosity=helios.LogVerbosity.VERBOSE, format=helios.OutputFormat.NPY
)
Reading Spectral Library...
10 materials found
Number of subsampling rays (riegl_lms-q780): 19
Simulation: Scanner changed!
Pulse frequency set to 100000
Scan angle set to 30
Applying settings for PolygonMirrorBeamDeflector...
Vertical angle min/max nan/nan degrees
 -- verticalAngleMin not set, using the value of -30 degrees
 -- verticalAngleMax not set, using the value of 30 degrees
SOURCE Leg with serial ID:0 waypoints:
        Origin: (-145, 50, 485.886)
        Target: (-145, -50, 485.886)
        Next: (145, -50, 485.886)

Starting simulation loop 1 ...
Waypoint reached!
Pulse frequency set to 100000
Scan angle set to 30
Applying settings for PolygonMirrorBeamDeflector...
Vertical angle min/max nan/nan degrees
 -- verticalAngleMin not set, using the value of -30 degrees
 -- verticalAngleMax not set, using the value of 30 degrees
SOURCE Leg with serial ID:1 waypoints:
        Origin: (-145, -50, 485.886)
        Target: (145, -50, 485.886)
        Next: (145, 50, 485.886)

Waypoint reached!
Pulse frequency set to 100000
Scan angle set to 30
Applying settings for PolygonMirrorBeamDeflector...
Vertical angle min/max nan/nan degrees
 -- verticalAngleMin not set, using the value of -30 degrees
 -- verticalAngleMax not set, using the value of 30 degrees
SOURCE Leg with serial ID:2 waypoints:
        Origin: (145, -50, 485.886)
        Target: (145, 50, 485.886)
        Next: (145, 50, 485.886)

Waypoint reached!
Pulse frequency set to 100000
Scan angle set to 30
Applying settings for PolygonMirrorBeamDeflector...
Vertical angle min/max nan/nan degrees
 -- verticalAngleMin not set, using the value of -30 degrees
 -- verticalAngleMax not set, using the value of 30 degrees
Waypoint reached!
Finishing simulation loop 1 ...
Finished simulation loop 1.
Elapsed simulation steps = 980004
Elapsed virtual time = 9.80004 sec.
Main thread simulation loop finished in 0.463468 sec.
Waiting for completion of pulse computation tasks...
Pulse computation tasks finished in 0.463476 sec.

Visualizing the results

Now we can display the resulting point clouds.

1) ULS point cloud

[10]:
coords = points_uls["position"]
nor = points_uls["number_of_returns"]

# Matplotlib figure
fig = plt.figure(figsize=(15, 10))


# settings for a discrete colorbar
N = int(np.max(nor))
cmap = plt.get_cmap("RdYlBu_r", N)
# Scatter plot (coloured by number of returns).
ax = fig.add_subplot(projection="3d")
sc = ax.scatter(coords[:, 0], coords[:, 1], coords[:, 2], c=nor, cmap=cmap, s=0.02)

# Add axis labels.
ax.set_xlabel("$X$")
ax.set_ylabel("$Y$")
ax.set_zlabel("$Z$")

# set equal axes
box = (np.ptp(coords[:, 0]), np.ptp(coords[:, 1]), np.ptp(coords[:, 2]))
ax.set_box_aspect(box)

# Set title.
ax.set_title(label="ULS point cloud", fontsize=18)

tick_locs = (np.arange(1, N + 1) + 0.5) * (N) / (N + 1)
cbar = plt.colorbar(sc, ticks=tick_locs)

cbar.set_label("Number of returns", fontsize=15)
tick_labels = [str(i) for i in range(1, N + 1)]
cbar.ax.set_yticklabels(tick_labels)

# Display results
plt.show()
_images/08-als_uls_detailed_voxel_16_0.png

On the tree on the very left, we clearly see the outlines of the fairly large (25 cm) opaque voxels. Due to the large voxel size, the beams can barely penetrate into the canopy and generate mostly single returns at the surface of the outer voxels.

The second tree from the left already looks more realistic due to the finer voxel grid. Still, there are relatively little multiple returns.

The tree in the middle is using the “transmittive” mode, where we assume that each voxel is filled with randomly distributed infinitely small sized leaf scatterers (“turbid medium”). Each voxel has a leaf area density and a leaf angle distribution is defined for the entire crown. From this, the extinction coefficient \(\sigma\) is calculated. Based on \(\sigma\), HELIOS++ determines where and if a return is generated within the voxel using a stochastic model. In this mode, we get the highest number of returns.

The two trees on the right use opaque voxels again, but voxels are scaled, so that the voxel leaf area density is proportional to the scaled voxel cube area. This means that smaller voxel cubes are generated where leaf area density is low. While generated using a voxel grid of 25 cm (like for the first tree), the scaling results in a voxel model which the laser beam can penetrate into and hence more realistic point clouds and number of multiple returns.

For the scaled and the transmittive mode, it is necessary to first compute voxel-based leaf area density. HELIOS++ supports the output format of AMAPVox (Vincent et al. 2017), a LiDAR Voxelisation Software for the computation of transmittance and estimation of leaf area density.

2) ALS point cloud

[11]:
coords = points_als["position"]
nor = points_als["number_of_returns"]

# Matplotlib figure
fig = plt.figure(figsize=(15, 10))

# settings for a discrete colorbar
N = int(np.max(nor))
cmap = plt.get_cmap("jet", N)
# Scatter plot (coloured by number of returns).
ax = fig.add_subplot(projection="3d")
sc = ax.scatter(
    coords[:, 0],
    coords[:, 1],
    coords[:, 2],
    c=nor,
    cmap=cmap,
    s=1.2,
    label="scene",
)

# Add axis labels.
ax.set_xlabel("$X$")
ax.set_ylabel("$Y$")
ax.set_zlabel("$Z$")

# set equal axes
box = (np.ptp(coords[:, 0]), np.ptp(coords[:, 1]), np.ptp(coords[:, 2]))
ax.set_box_aspect(box)

# Set title.
ax.set_title(label="ALS point cloud", fontsize=18)

tick_locs = (np.arange(1, N + 1) + 0.5) * (N) / (N + 1)
cbar = plt.colorbar(sc, ticks=tick_locs)

cbar.set_label("Number of returns", fontsize=15)
tick_labels = [str(i) for i in range(1, N + 1)]
cbar.ax.set_yticklabels(tick_labels)

# Display results
plt.show()
_images/08-als_uls_detailed_voxel_19_0.png
[ ]: