Basic M3C2-EP algorithm

This presents how the M3C2-EP algorithm (Winiwarter et al., 2021) 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).

Please ensure two datasets include scan positions which are specified by attribute name sp_name and scan positions configuration information in sp_file.

Here, we work with a rather small data set:

epoch1, epoch2 = py4dgeo.read_from_las(
    "ahk_2017_652900_5189100_gnd_subarea.laz",
    "ahk_2018A_652900_5189100_gnd_subarea.laz",
    additional_dimensions={"point_source_id": "scanpos_id"},
)

[2024-05-14 13:00:50][INFO] Reading point cloud from file '/home/docs/.cache/py4dgeo/./ahk_2017_652900_5189100_gnd_subarea.laz'
[2024-05-14 13:00:50][INFO] Reading point cloud from file '/home/docs/.cache/py4dgeo/./ahk_2018A_652900_5189100_gnd_subarea.laz'

M3C2-EP leverages knowledge regarding the data acquisition sensor. We extract the sensor orientation details from a JSON configuration file and assign them to a dictionary. These settings are then applied to each epoch of the data since both epochs share the same sensor configuration.

with open(py4dgeo.find_file("sps.json"), "r") as load_f:
    scanpos_info_dict = eval(load_f.read())
epoch1.scanpos_info = scanpos_info_dict
epoch2.scanpos_info = scanpos_info_dict

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, or something else. Here, we choose a subsampling by taking every 50th point of the reference point cloud:

corepoints = py4dgeo.read_from_las("ahk_cp_652900_5189100_subarea.laz").cloud

[2024-05-14 13:00:51][INFO] Reading point cloud from file '/home/docs/.cache/py4dgeo/./ahk_cp_652900_5189100_subarea.laz'

The algorithm needs an alignment covariance matrix of shape 12 x 12, an affine transformation matrix \(T\) of shape 3 x 4 and a reduction point \((x_0, y_0, z_0)^T\) (rotation origin, 3 parameters) obtained from aligning the two point clouds. The transformation follows this notation:

\[\begin{split}T = \begin{pmatrix} t_1 & t_2 & t_3 & t_4 \\ t_5 & t_6 & t_7 & t_8 \\ t_9 & t_{10} & t_{11} & t_{12} \\ \end{pmatrix}\end{split}\]

Where the transformation is applied as follows:

\[\begin{split}y = \begin{pmatrix} t_1 & t_2 & t_3 \\ t_5 & t_6 & t_7 \\ t_9 & t_{10} & t_{11} \\ \end{pmatrix} \left( x-\begin{pmatrix} x_0 \\ y_0 \\ z_0 \\ \end{pmatrix} \right) + \begin{pmatrix} t_4 \\ t_8 \\ t_{12} \\ \end{pmatrix} + \begin{pmatrix} x_0 \\ y_0 \\ z_0 \\ \end{pmatrix}\end{split}\]

The order of the elements in the covariance matrix is:

\[t_1, t_2, t_3, t_4, t_5, t_6, t_7, t_8, t_9, t_{10}, t_{11}, t_{12}\]

, meaning that transformation and rotation/scaling parameters are interleaved.

We can decide whether to perform the transformation by a boolean flag ‘perform_trans’ and it is performed by default.

Cxx = np.loadtxt(py4dgeo.find_file("Cxx.csv"), dtype=np.float64, delimiter=",")
tfM = np.loadtxt(py4dgeo.find_file("tfM.csv"), dtype=np.float64, delimiter=",")
refPointMov = np.loadtxt(
    py4dgeo.find_file("redPoint.csv"), dtype=np.float64, delimiter=","
)

Next, we instantiate the algorithm class and run the distance calculation. The parameters are very similar to the base M3C2 implementation, but extended to work for M3C2-EP.

m3c2_ep = py4dgeo.M3C2EP(
    epochs=(epoch1, epoch2),
    corepoints=corepoints,
    normal_radii=(0.5, 1.0, 2.0),
    cyl_radii=(0.5,),
    max_distance=3.0,
    Cxx=Cxx,
    tfM=tfM,
    refPointMov=refPointMov,
)

distances, uncertainties, covariance = m3c2_ep.run()

M3C2EP running
[2024-05-14 13:00:51][INFO] Building KDTree structure with leaf parameter 10
[2024-05-14 13:00:52][INFO] Building KDTree structure with leaf parameter 10
M3C2EP end

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. In contrast to the base M3C2, additional covariance information is returned as a third element in the tuple.

distances[0:8]

array([-0.17164396, -0.62896535, -0.5648335 , -0.15689298, -0.12769975,
       -0.46240109, -0.18725634, -1.08448984])

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"][0:8]

array([0.02308354, 0.02319384, 0.02437987, 0.02349072, 0.02414121,
       0.02387164, 0.02578069, 0.03018877])

uncertainties["spread1"][0:8]

array([0.00689291, 0.00695656, 0.00744346, 0.00719902, 0.0074957 ,
       0.00757502, 0.00818117, 0.0084822 ])

uncertainties["num_samples1"][0:8]

array([119, 122, 163,  11, 281,  41, 199, 172])

Corresponding to the derived distances, a 3D covariance information for the point cloud is returned.

covariance["cov1"][0, :, :]

array([[ 5.99898020e-05,  1.54621578e-05, -5.85061284e-06],
       [ 1.54621578e-05,  5.78886408e-05, -5.46655748e-06],
       [-5.85061284e-06, -5.46655748e-06,  4.95475375e-05]])

Finally we could visualize our distances results.

import matplotlib.cm as cm
import matplotlib.pyplot as plt


def plt_3d(corepoints, distances):
    fig, ax = plt.subplots(figsize=(10, 10), subplot_kw={"projection": "3d"})

    # add axis labels
    ax.set_xlabel("X [m]")
    ax.set_ylabel("Y [m]")
    ax.set_zlabel("Z [m]")

    # plot the corepoints colored by their distance
    x, y, z = np.transpose(corepoints)
    vmin = np.min(distances)
    vmax = np.max(distances)
    pts = ax.scatter(
        x, y, z, s=10, c=distances, vmin=vmin, vmax=vmax, cmap=cm.seismic_r
    )

    # add colorbar
    cmap = plt.colorbar(pts, shrink=0.5, label="Distance [m]", ax=ax)

    # add title
    ax.set_title("Visualize Changes")

    ax.set_aspect("equal")
    ax.view_init(22, 113)
    plt.show()

plt_3d(corepoints, distances)

References

• Winiwarter, L., Anders, K., & Höfle, B. (2021). M3C2-EP: Pushing the limits of 3D topographic point cloud change detection by error propagation. ISPRS Journal of Photogrammetry and Remote Sensing, 178, 240-258. doi: 10.1016/j.isprsjprs.2021.06.011.

• 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.