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"
)
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()
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
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"]
uncertainties["spread1"]
uncertainties["num_samples1"]
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.
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.
m3c2.directions_radii()