Point Cloud Registration

In order to lower the uncertainty of calculations done with py4dgeo, point cloud registration should be performed to tightly align the point clouds. py4dgeo supports this in two ways:

• It allows you to apply arbitrary affine transformations to epochs and stores the transformation parameters as part of the epoch. You can use the tooling of your choice to calculate these.

• It provides a number of registration algorithms that allow you to calculate the transformation directly in py4dgeo. Currently, this is limited to standard ICP.

This notebook show cases both ways of usage.

Note: Be aware that an initial transformation is required that roughly aligns the two point clouds.

import py4dgeo
import numpy as np

We load two epochs and register them to each other:

epoch1 = py4dgeo.read_from_xyz("plane_horizontal_t1.xyz")  # replace with own data
epoch2 = py4dgeo.read_from_xyz("plane_horizontal_t2.xyz")  # replace with own data

[2024-05-14 13:04:38][INFO] Reading point cloud from file '/home/docs/.cache/py4dgeo/./plane_horizontal_t1.xyz'
[2024-05-14 13:04:38][INFO] Reading point cloud from file '/home/docs/.cache/py4dgeo/./plane_horizontal_t2.xyz'

Registration algorithms use a simple function call interface. It returns an affine transformation as a 3x4 matrix which contains the Rotation matrix \(\mathbf{R}\) and the translation vector \(\mathbf{t}\):

\[\begin{split}\mathbf{T} = \left(\begin{array}{ccc}r_{11}&r_{12}&r_{13}&t_1\\r_{21}&r_{22}&r_{23}&t_2\\r_{31}&r_{32}&r_{33}&t_3\end{array}\right)\end{split}\]

When calculating and applying such transformations in py4dgeo, it is always possible to specify the reduction point \(\mathbf{x_0}\). This allows shifting coordinates towards the origin before applying rotations - leading to a numerically robust operation. The algorithm for a transformation is:

\[\mathbf{Tx} = \left(\mathbf{R}(\mathbf{x}-\mathbf{x_0})\right) + \mathbf{t} + \mathbf{x_0}\]

trafo = py4dgeo.iterative_closest_point(
    epoch1, epoch2, reduction_point=np.array([0, 0, 0])
)

[2024-05-14 13:04:38][INFO] Building KDTree structure with leaf parameter 10

This transformation can then be used to transform epoch2:

epoch2.transform(trafo)

## or use an external transformation matrix
# external_trafo = np.loadtxt("trafo.txt") # replace with own data
# epoch2.transform(external_trafo)

The Epoch class records applied transformations and makes them available through the transformation property:

epoch2.transformation

[Transformation(affine_transformation=array([[ 1.00000000e+00,  6.34089843e-10, -1.42106695e-05,
         -7.59183919e-03],
        [-1.25012920e-10,  9.99999999e-01,  3.58235678e-05,
          7.09922947e-03],
        [ 1.42106695e-05, -3.58235678e-05,  9.99999999e-01,
          2.11030351e+01],
        [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
          1.00000000e+00]]), reduction_point=array([0, 0, 0]))]

Finally, we can export the transformed epoch to inspect the registration quality e.g. in CloudCompare.

# epoch2.save("plane_horizontal_t2_transformed.xyz") # replace with own data

Available registration algorithms

Currently only standard point to point ICP is available, but other algorithms are currently implemented:

?py4dgeo.iterative_closest_point