Basic M3C2 algorithm

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 numpy as np
import py4dgeo

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:

epoch1, epoch2 = py4dgeo.read_from_xyz(
    "plane_horizontal_t1.xyz", "plane_horizontal_t2.xyz"
)

[2024-05-14 13:00:45][INFO] Reading point cloud from file '/home/docs/.cache/py4dgeo/./plane_horizontal_t1.xyz'
[2024-05-14 13:00:45][INFO] Reading point cloud from file '/home/docs/.cache/py4dgeo/./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_radii=(2.0,),
    normal_radii=(0.5, 1.0, 2.0),
)

distances, uncertainties = m3c2.run()

[2024-05-14 13:00:45][INFO] Building KDTree structure with leaf parameter 10
[2024-05-14 13:00:45][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.10269298, -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.00565384, 0.00181819, 0.00255157, 0.00332511, 0.00281475,
       0.00299644, 0.00295429, 0.00199006, 0.00360769])

uncertainties["spread1"]

array([0.00417673, 0.00193209, 0.00203436, 0.00368965, 0.00259167,
       0.00286324, 0.0032822 , 0.00247362, 0.00356988])

uncertainties["num_samples1"]

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

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.