pointtree.operations

Operations for point cloud processing.

pointtree.operations.cloth_simulation_filtering( coords: ndarray[Any, dtype[float32 | float64]], classification_threshold: float, resolution: float, rigidness: int, correct_steep_slope: bool = False, iterations: int = 100, ) → ndarray[Any, dtype[int64]][source]

Detects ground points using the Cloth Simulation Filtering (CSF) algorithm proposed in Zhang, Wuming, et al. “An Easy-to-Use Airborne LiDAR Data Filtering Method Based on Cloth Simulation.” Remote Sensing 8.6 (2016): 501.

Parameters:

coords (ndarray[Any, dtype[float32 | float64]]) – Point coordinates.
classification_threshold (float) – Maximum height above the cloth a point can have in order to be classified as terrain point. All points whose distance to the cloth is equal or below this threshold are classified as terrain points.
resolution (float) – Resolution of the cloth grid (in meter).
rigidness (int) – Rigidness of the cloth (the three levels 1, 2, and 3 are available, where 1 is the lowest and 3 the highest rigidness).
correct_steep_slope (bool) – Whether the cloth should be corrected for steep slopes in a post-pressing step. (default: False)
iterations (int) – Maximum number of iterations. (default: 100)

Returns:

Class IDs for each point. For terrain points, the class ID is set to 0 and for non-terrain points to 1.

Raises:

ValueError – If rigidness is not 1, 2, or 3.

Shape:

coords: \((N, 3)\)
Output: \((N)\).

where

\(N = \text{ number of points}\)

pointtree.operations.create_digital_terrain_model( terrain_xyz: ndarray[Any, dtype[float32 | float64]], grid_resolution: float, k: int, p: float, voxel_size: float | None = None, num_workers: int = 1, ) → Tuple[ndarray[Any, dtype[float32 | float64]], ndarray[Any, dtype[float32 | float64]]][source]

Constructs a rasterized digital terrain model (DTM) from a set of terrain points. The DTM is constructed by creating a grid of regularly arranged DTM points and interpolating the height of the \(k\) closest terrain points for each DTM point on the grid. In the interpolation, terrain points \(x_t\) are weighted with a factor proportional to a power of \(p\) of their inverse distance to the corresponding DTM point \(x_{dtm}\), i.e., \(\frac{1}{||(x_{dtm} - x_{t})||^p}\). If there are terrain points whose distance to the DTM point is zero, only these points are used to calculate the DTM height and more distant points are ignored. Before constructing the DTM, the terrain points can optionally be downsampled using voxel-based subsampling.

Parameters:

terrain_xyz (ndarray[Any, dtype[float32 | float64]]) – Coordinates of the terrain points from which to construct the DTM.
grid_resolution (float) – Resolution of the DTM grid (in meter).
k (int) – Number of terrain points between which interpolation is performed to obtain the terrain height of a DTM point.
p (float) – Power \(p\) for inverse-distance weighting in the interpolation of terrain points.
voxel_size (float | None) – Voxel size with which the terrain points are downsampled before the DTM is created. If set to None, no downsampling is performed. (default: None)
num_workers (int) – Number of workers to use for parallel processing. If workers is set to -1, all CPU threads are used. (default: 1)

Returns:

Tuple of two arrays:

DTM
X- and y-coordinate of the top left corner of the DTM grid

Shape:

terrain_xyz: \((N, 3)\)
Output: \((H, W)\), \((2)\)

where

\(N = \text{ number of terrain points}\)
\(H = \text{ extent of the DTM grid in y-direction}\)
\(W = \text{ extent of the DTM grid in x-direction}\)

pointtree.operations.distance_to_dtm( xyz: ndarray[Any, dtype[float32 | float64]], dtm: ndarray[Any, dtype[float32 | float64]], dtm_offset: ndarray[Any, dtype[float32 | float64]], dtm_resolution: float, allow_outside_points: bool = True, ) → ndarray[Any, dtype[float32 | float64]][source]

Compute the height above the terrain for each point of a point cloud by subtracting the corresponding terrain height from the z-coordinate of the point. The terrain height for a given point is obtained by bilinearly interpolating the terrain heights of the four closest grid points of the digital terrain model.

Parameters:

xyz (ndarray[Any, dtype[float32 | float64]]) – Point coordinates of the point cloud to normalize.
dtm (ndarray[Any, dtype[float32 | float64]]) – Rasterized digital terrain model.
dtm_offset (ndarray[Any, dtype[float32 | float64]]) – X and y-coordinates of the top left corner of the DTM grid.
allow_outside_points (bool) – If this option is set to True and a point in the point cloud to be normalized is not in the area covered by the DTM, the height of the nearest DTM points is still determined and used for normalization. Otherwise, a ValueError is thrown if points are outside the area covered by the DTM. (default: True)

Returns:

Height above the terrain of each point.

Raises:

ValueError – If the point cloud to be normalized covers a larger base area than the DTM.

Shape:

xyz: \((N, 3)\)
dtm: \((H, W)\)
dtm_offset: \((2)\)
Output: \((N)\)

where

\(N = \text{ number of points}\)
\(H = \text{ extent of the DTM in grid in y-direction}\)
\(W = \text{ extent of the DTM in grid in x-direction}\)

pointtree.operations.estimate_with_linear_model( x_train: ndarray[Any, dtype[float32 | float64]], y_train: ndarray[Any, dtype[float32 | float64]], x_predict: ndarray[Any, dtype[float32 | float64]], ) → Tuple[ndarray[Any, dtype[float32 | float64]], LinearRegression][source]

Fits a linear model to the training data x_train and y_train, and computes the predictions of the model for x_predict.

Parameters:

x_train (ndarray[Any, dtype[float32 | float64]]) – Training data.
y_train (ndarray[Any, dtype[float32 | float64]]) – Target values.
x_predict (ndarray[Any, dtype[float32 | float64]]) – Inference data.

Returns:

Tuple of two elements:

An array containing the predictions for x_predict
The fitted linear model

Shape:

x_train: \((N)\) or \((N, D)\)
y_train: \((N)\)
x_predict: \((N')\) or \((N', D)\)
Output: \((N')\)

where

\(N = \text{ number of training samples}\)
\(N' = \text{ number of inference samples}\)
\(D = \text{ number of input features}\)

pointtree.operations.fit_ellipse( xy: ndarray[Any, dtype[float32 | float64]], batch_lengths: ndarray[Any, dtype[int64]], num_workers: int = 1, ) → ndarray[Any, dtype[float32 | float64]][source]

Fits an ellipse to a set of 2D points using the least-squares method described in Halir, Radim, and Jan Flusser. “Numerically Stable Direct Least Squares Fitting of Ellipses.” Proc. 6th International Conference in Central Europe on Computer Graphics and Visualization. WSCG. Vol. 98. Plzen-Bory: Citeseer, 1998. This method supports batch processing, i.e., ellipses can be fitted to separate sets of points (batch items) in parallel. For this purpose, batch_lengths must be set to specify which point belongs to which set.

Parameters:

xy (ndarray[Any, dtype[float32 | float64]]) – X- and y- coordinates of the points to which the ellipses are to be fitted. If the xy array has a row-major storage layout (numpy’s default), a copy of the array is created. To pass xy by reference, xy must be in column-major format.
batch_lengths (ndarray[Any, dtype[int64]]) – Number of points in each item of the input batch. For batch processing, it is expected that all points belonging to the same batch item are stored consecutively in the xy input array. For example, if a batch comprises two batch items with \(N_1\) points and \(N_2\) points, then batch_lengths should be set to [N_1, N_2] and xy[:N_1] should contain the points of the first batch item and xy[N_1:] the points of the second batch item. If batch_lengths is set to None, it is assumed that the input points belong to a single batch item and batch processing is disabled.
num_workers (int) – Number of workers threads to use for parallel processing. If set to -1, all CPU threads are used. (default: 1)

Returns:

Parameters of the fitted ellipses in the following order: X- and y-coordinates of the center, radius along the semi-major and along the semi-minor axis, and the counterclockwise angle of rotation from the x-axis to the semi-major axis of the ellipse. If no ellipse is detected for a batch item, all ellipse parameters for this batch item are set to -1.

Raises:

TypeError – If xy or batch_length have an invalid shape or data type.
ValueError – If the length of xy is not equal to the sum of batch_lengths.

Shape:

xy: \((N, 2)\)
batch_lengths: \((B)\)
Output: \((B, 5)\)

where

\(B = \text{ batch size}\)

\(N = \text{ number of points}\)

pointtree.operations.points_in_ellipse( xy: ndarray[Any, dtype[float32 | float64]], ellipse: ndarray[Any, dtype[float32 | float64]], ) → ndarray[Any, dtype[bool_]][source]

Tests whether 2D points are within the boundaries of an ellipse.

If the input arrays have a row-major storage layout (numpy’s default), a copy of the input arrays is created. To pass them by reference, they must be in column-major format.

Parameters:

xy (ndarray[Any, dtype[float32 | float64]]) – Coordinates of the points to test.
ellipse (ndarray[Any, dtype[float32 | float64]]) – Parameters of the ellipse in the following order: X- and y-coordinates of its center, radius along the semi-major and along the semi-minor axis, and the counterclockwise angle of rotation from the x-axis to the semi-major axis of the ellipse.

Returns:

Boolean array that indicates for each point whether it lies within the ellipse.

Raises:

ValueError – If the shape of xy or ellipse is invalid.

Shape:

xy: \((N, 2)\)
ellipse: \((5)\)
Output: \((N)\)

where

\(N = \text{ number of points}\)

pointtree.operations.polygon_area( x: ndarray[Any, dtype[float32 | float64]], y: ndarray[Any, dtype[float32 | float64]], ) → float[source]

Computes the area of a convex polygon using the Shoelace formula. It is expected that the input vertices are already sorted so that adjacent vertices are neighbors in the input vertex array.

Parameters:

x (ndarray[Any, dtype[float32 | float64]]) – X-coordinates of the polygon vertices.
y (ndarray[Any, dtype[float32 | float64]]) – Y-coordinates of the polygon vertices.

Returns:

Polygon area.

Shape:

x: \((N)\)
y: \((N)\)

where

\(n = \text{ number of points}\)

pointtree.operations.statistical_outlier_removal( xyz: ndarray[Any, dtype[float32 | float64]], k: int, std_multiplier: float, ) → Tuple[ndarray[Any, dtype[float32 | float64]], ndarray[Any, dtype[int64]]][source]

Statistical outlier filter for 3D point clouds: First, it computes the average distance that each point has to its k nearest neighbors. Next, the mean and standard deviation of these average distances are computed to determine a distance threshold. The distance threshold is set to: mean + std_multiplier * stddev. Then, points are classified as outlier and removed if their average neighbor distance is above this threshold.

Parameters:

xyz (ndarray[Any, dtype[float32 | float64]]) – Coordinates of the points to be filtered.
k (int) – Number of neighboring points to consider for the filtering (the point itself is not included).
std_multiplier (float) – Multiplier for the standard deviation in the calculation of the distance threshold.

Returns:

Tuple of two arrays:

Coordinates of the inlier points remaining after filtering.
Indices of the inlier points remaining after filtering with respect to the input array.

Shape:

xyz: \((N, 3)\)
Output: \((N', 3)\), \((N')\)

where

\(N = \text{ number of points before the filtering}\)

\(N' = \text{ number of points after the filtering}\)