Hierarchical change analysis

This notebook demonstrates how the hierarchical change analysis algorithm (Tabernig et al., 2025) can be run with py4dgeo using the vapc (Voxel Analysis for Point Clouds) submodule.

Implemented by

Ronald Tabernig (@tabernig)

Author(s) of the method Ronald Tabernig, Bernhard Höfle (Heidelberg University)

Original publication of the method Tabernig, R., Albert, W., Weiser, H., & Höfle, B. (2025). A hierarchical approach for near real-time 3D surface change analysis of permanent laser scanning point clouds. In: 6th Joint International Symposium on Deformation Monitoring (JISDM). doi: 10.5445/IR/1000180377

As a first step, we import the py4dgeo and numpy packages:

import py4dgeo
import numpy as np
import pooch

Numba is not installed.
Use `pip install numba` to install it to enable faster computations in Vapc.
You can run py4dgeo.Vapc without Numba, but it may be slower.

Next, we need to load two point clouds of the same scene taken at different times and specify where to store the output.

# Set up pooch to download data from Zenodo
p = pooch.Pooch(base_url="doi:10.5281/zenodo.18432391/", path=pooch.os_cache("py4dgeo"))
p.load_registry_from_doi()

try:
    # Download and extract the dataset
    p.fetch("trier_sim.zip", processor=pooch.Unzip(members=["trier_sim"]))

    # Define path to the extracted data
    data_path = p.path / "trier_sim.zip.unzip"
    print(f"Data path: {data_path}")

    before_rockfall_file = (
        data_path / "trier_sim_epoch_0.laz"
    )  # Synthetic data of terrain before a simulated rockfall at Trier study site
    after_rockfall_file = (
        data_path / "trier_sim_epoch_1.laz"
    )  # Synthetic data of terrain after a simulated rockfall at Trier study site

    epoch1, epoch2 = py4dgeo.read_from_las(before_rockfall_file, after_rockfall_file)


except Exception as e:
    print(f"Failed to download or extract data: {e}")

Downloading file 'trier_sim.zip' from 'doi:10.5281/zenodo.18432391/trier_sim.zip' to '/home/docs/.cache/py4dgeo'.
Extracting 'trier_sim' from '/home/docs/.cache/py4dgeo/trier_sim.zip' to '/home/docs/.cache/py4dgeo/trier_sim.zip.unzip'

Data path: /home/docs/.cache/py4dgeo/trier_sim.zip.unzip
[2026-06-12 12:50:32][INFO] Reading point cloud from file '/home/docs/.cache/py4dgeo/trier_sim.zip.unzip/trier_sim_epoch_0.laz'
[2026-06-12 12:50:32][INFO] Reading point cloud from file '/home/docs/.cache/py4dgeo/trier_sim.zip.unzip/trier_sim_epoch_1.laz'

outfile_las = str(before_rockfall_file).replace(
    ".laz", "_hierarchical_change_analysis_result.las"
)
outfile_ply = str(before_rockfall_file).replace(
    ".laz", "_hierarchical_change_analysis_result.ply"
)

The hierarchical change analysis is based on the rapid detection of changes in voxelised point clouds, followed by a detailed 3D surface analysis of the points, which is applied only to areas where changes were detected in the first step. Accordingly, we need to define the voxel size (voxel_size), the significance threshold (alpha) and the minimum number of points per voxel required to check for statistically significant changes (min_points).

voxel_size = 6
alpha = 0.999
min_points = 30

We can now check for voxels that have changed significantly. First, we convert our Epoch objects into Vapc objects. We can then compute the bitemporal Mahalanobis distance using one of these Vapc objects.

# Mute vapc function trace and timeit for cleaner output
py4dgeo.enable_trace(False)
py4dgeo.enable_timeit(False)

voxel_epoch1 = py4dgeo.Vapc(epoch1, voxel_size=voxel_size)
voxel_epoch2 = py4dgeo.Vapc(epoch2, voxel_size=voxel_size)

# Compute delta Vapc
mahalanobis_result = voxel_epoch1.compute_bitemporal_mahalanobis(
    voxel_epoch2, alpha=alpha, min_points=min_points
)

This intermediate result indicates whether significant change was detected or not. For a detailed analysis of 3D surface changes, we only need to compute changes in areas where significant changes have been detected. Accordingly, we extract points from voxels with significant changes, reducing our Vapc object to these points.

# Filter significant changes
sig_filter = mahalanobis_result.out["significance"] == 1
# Apply the filter to the delta Vapc object - this will keep only the significant changes in the Vapc object
mahalanobis_result.filter(sig_filter, overwrite=True)
# Select points with significant changes in epoch 1
# This will return a new Vapc object with only the points that have significant changes
voxel_epoch_1_with_significant_change, mask = voxel_epoch1.select_by_mask(
    mahalanobis_result
)

Now we use only these points with significant changes as (core-)points for subsequent change analysis, e.g., with the well-known M3C2 algorithm.

m3c2 = py4dgeo.M3C2(
    epochs=(epoch1, epoch2),
    corepoints=voxel_epoch_1_with_significant_change.epoch.cloud,
    cyl_radius=1.0,
    normal_radii=[1.0],
    max_distance=10.0,
)

distances, uncertainties = m3c2.run()

[2026-06-12 12:50:33][INFO] Building KDTree structure with leaf parameter 10
[2026-06-12 12:50:33][INFO] Building KDTree structure with leaf parameter 10

In the final step, we update the Vapc object with the computed M3C2 results and add a field to indicate whether the detected M3C2 change is significant. Then, we save the file to the specified output path.

voxel_epoch_1_with_significant_change.out["M3C2_distance"] = distances
d = {name: uncertainties[name] for name in uncertainties.dtype.names}
d_filtered = {k: v for k, v in d.items() if k != "dtype"}
voxel_epoch_1_with_significant_change.out.update(d)
voxel_epoch_1_with_significant_change.out["significant_change"] = (
    np.abs(distances) > uncertainties["lodetection"]
)

# Save the Vapc object with M3C2 results
voxel_epoch_1_with_significant_change.save_as_las(outfile_las)
print(f"Results saved to folder: {data_path}")

Adding extra dimension 'M3C2_distance' of type float64
Adding extra dimension 'lodetection' of type float64
Adding extra dimension 'spread1' of type float64
Adding extra dimension 'num_samples1' of type int64
Adding extra dimension 'spread2' of type float64
Adding extra dimension 'num_samples2' of type int64
Adding extra dimension 'significant_change' of type <class 'numpy.uint8'>
Results saved to folder: /home/docs/.cache/py4dgeo/trier_sim.zip.unzip

We may also wish to save the voxels as 3D cubes. The save_as_ply function accomplishes this by saving occupied voxels as cubes in a triangle mesh. The edge length of these cubes is defined by the voxel size set before. The features to be stored with each voxel must be listed. In this example, we select all available features. The mode option allows us to define the center of each sube. It has the following options:

• “closest_to_centroids”

• “closest_to_voxel_centers”

• “centroid” and

• “voxel_center”

We use the “voxel_center” which results in cubes centered on the voxel grid. Before saving, we reduce the Vapc object to one point per voxel using the reduce_to_feature method (see the notebook on spatial subsampling), which takes the same mode options as described above. This ensures that we du not write duplicate voxels.

# Let's reduce our point cloud to one point per voxel to ensure that we don't write duplicate voxels.
try:
    reduce_to_mode = "closest_to_voxel_centers"  # other options are "closest_to_centroids", "closest_to_voxel_centers", "centroid", "voxel_center"
    reduced_vapc = voxel_epoch_1_with_significant_change.reduce_to_feature(
        reduce_to_mode
    )

    reduced_vapc.save_as_ply(
        outfile=outfile_ply, features=reduced_vapc.out.keys(), mode=reduce_to_mode
    )
    print(f"Results saved to folder: {data_path}")
except:
    print("Failed to save PLY file. Check if 'plyfile' is installed.")
    print(
        "You can try installing it by uncommenting and running the following lines of code in a new cell:"
    )
    print("import sys")
    print("!conda install --yes --prefix {sys.prefix} conda-forge::plyfile")

plyfile is not installed.
Use `pip install plyfile` to install it to enable the save_as_ply functionality in Vapc.
Failed to save PLY file. Check if 'plyfile' is installed.
You can try installing it by uncommenting and running the following lines of code in a new cell:
import sys
!conda install --yes --prefix {sys.prefix} conda-forge::plyfile

References

Tabernig, R., Albert, W., Weiser, H., & Höfle, B. (2025). A hierarchical approach for near real-time 3D surface change analysis of permanent laser scanning point clouds. In: 6th Joint International Symposium on Deformation Monitoring (JISDM). doi: 10.5445/IR/1000180377