Basic M3C2 algorithm

Implemented by

Dominic Kempf (@dokempf, Heidelberg University)

Author(s) of the method Dimitri Lague, Nicolas Brodu, Jérôme Leroux (University of Rennes)

Original publication of the method Lague, D., Brodu, N., & Leroux, J. (2013). Accurate 3D comparison of complex topography with terrestrial laser scanner: Application to the Rangitikei canyon (N-Z). ISPRS Journal of Photogrammetry and Remote Sensing, 82, pp. 10-26. doi: 10.1016/j.isprsjprs.2013.04.009.

This presents how the M3C2 algorithm (Lague et al., 2013) for point cloud distance computation can be run using the py4dgeo package. As a first step, we import the py4dgeo and numpy packages:

import py4dgeo

import numpy as np

from py4dgeo.util import find_file

Next, we need to load two datasets that cover the same scene at two different points in time. Point cloud datasets are represented by numpy arrays of shape n x 3 using a 64 bit floating point type (np.float64). Here, we work with a rather small synthetical data set:

# Fetch original test data
test_filename = "usage_data.zip"
original_file = find_file(test_filename)
print(f"File downloaded to: {original_file}")

File downloaded to: /home/docs/.cache/py4dgeo/usage_data.zip

# Read XYZ files from the extracted directory
epoch1, epoch2 = py4dgeo.read_from_xyz(
    find_file("plane_horizontal_t1.xyz"), find_file("plane_horizontal_t2.xyz")
)

[2026-06-23 06:35:44][INFO] Reading point cloud from file '/home/docs/.cache/py4dgeo/extracted/usage_data/plane_horizontal_t1.xyz'
[2026-06-23 06:35:44][INFO] Reading point cloud from file '/home/docs/.cache/py4dgeo/extracted/usage_data/plane_horizontal_t2.xyz'

The analysis of point cloud distances is executed on so-called core points (cf. Lague et al., 2013). These could be, e.g., one of the input point clouds, a subsampled version thereof, points in an equidistant grid, etc. Here, we choose a subsampling by taking every 50th point of the reference point cloud:

corepoints = epoch1.cloud[::50]

Next, we instantiate the algorithm class and run the distance calculation:

m3c2 = py4dgeo.M3C2(
    epochs=(epoch1, epoch2),
    corepoints=corepoints,
    cyl_radius=2.0,
    normal_radii=[0.5, 1.0, 2.0],
)

distances, uncertainties = m3c2.run()

[2026-06-23 06:35:44][INFO] Building KDTree structure with leaf parameter 10
[2026-06-23 06:35:44][INFO] Building KDTree structure with leaf parameter 10

The calculated result is an array with one distance per core point. The order of distances corresponds exactly to the order of input core points.

distances

array([-0.10269281, -0.09957772, -0.10179986, -0.10063675, -0.10091512,
       -0.10070981, -0.09819643, -0.09926434, -0.09911222])

Corresponding to the derived distances, an uncertainty array is returned which contains several quantities that can be accessed individually: The level of detection lodetection, the spread of the distance across points in either cloud (spread1 and spread2, by default measured as the standard deviation of distances) and the total number of points taken into consideration in either cloud (num_samples1 and num_samples2):

uncertainties["lodetection"]

array([0.0054738 , 0.00181819, 0.00255157, 0.00332511, 0.00281475,
       0.00299644, 0.0029543 , 0.00199006, 0.0036077 ])

uncertainties["spread1"]

array([0.0039509 , 0.00193209, 0.00203436, 0.00368965, 0.00259167,
       0.00286324, 0.00328221, 0.00247362, 0.00356988])

uncertainties["num_samples1"]

array([ 5, 11, 11, 11, 12, 10, 11, 12, 10])

The direction of surface change in the M3C2 algorithm is determined via local normal vectors per core point. The normal vectors used in the calculation can be accessed via the directions() method of the M3C2 algorithm in py4dgeo, which returns an array (Nx3) of length N corresponding to the number of core points with three entries for the normal vector components in x, y, and z direction.

directions = m3c2.directions()

The property directions_radii returns an array (Nx1) of length N corresponding to the number of core points and one entry for the radius used for normal computation at the respective core point. This is relevant for the multi-scale functionality of the M3C2, i.e. the possibility to specify multiple normal radii of which the one with maximized planarity is used for change analysis.

direction_radii = m3c2.directions_radii()

After running the M3C2 algorithm, you can save the results as a LAS file. The write_m3c2_results_to_las function exports the corepoints along with their calculated attributes.

outputfilepath = "m3c2_results.las"
attributes_to_save = {
    "distance": distances,  # The M3C2 distances
    "lodetection": uncertainties["lodetection"],  # Level of detection
    "spread1": uncertainties["spread1"],  # Spread in first cloud
    "spread2": uncertainties["spread2"],  # Spread in second cloud
    "num_points1": uncertainties["num_samples1"],  # Number of points in first cloud
    "num_points2": uncertainties["num_samples2"],  # Number of points in second cloud
    "normal_x": directions[:, 0],  # X component of normal vector
    "normal_y": directions[:, 1],  # Y component of normal vector
    "normal_z": directions[:, 2],  # Z component of normal vector
    "normal_radius": direction_radii,  # Radius used for normal computation
}

py4dgeo.write_m3c2_results_to_las(outputfilepath, m3c2, attributes_to_save)

import os
from pathlib import Path

if os.path.exists(outputfilepath):
    print(f"Results successfully saved to: {os.path.abspath(outputfilepath)}")
else:
    print("Failed to save results.")

Results successfully saved to: /home/docs/checkouts/readthedocs.org/user_builds/py4dgeo/checkouts/latest/doc/m3c2_results.las

References

• Lague, D., Brodu, N., & Leroux, J. (2013). Accurate 3D comparison of complex topography with terrestrial laser scanner: Application to the Rangitikei canyon (N-Z). ISPRS Journal of Photogrammetry and Remote Sensing, 82, pp. 10-26. doi: 10.1016/j.isprsjprs.2013.04.009.