Compute eigenvalues for large point clouds as quickly as possible with R

Question

I have large point clouds (each with several million points) and want to compute geometric attributes for all points. In order to do so, I need to calculate eigenvalues. How can I do this as quickly as possible? My data is stored in *.las files and I read them in with the package lidR. I also use this package to calculate point metrics. According to this post, I implemented this version:

# load data
cloud_raw <- readLAS(path_points)

# because eigen() is really slow, use C++
Rcpp::sourceCpp(code = "
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
arma::vec eigen_values(arma::mat A) {
arma::mat coeff, score;
arma::vec latent;
arma::princomp(coeff, score, latent, A);
return(latent);
}")

metric_geometry_features <- function(x, y, z) {
  xyz <- cbind(x, y, z)
  cov_matrix <- cov(xyz)
  eigen_matrix <- eigen_values(cov_matrix)
  geometries = list(
    planarity = (eigen_matrix[2] - eigen_matrix[3]) / eigen_matrix[1],
    linearity = (eigen_matrix[1] - eigen_matrix[2]) / eigen_matrix[1]
  )
  return(geometries)
}

metrics <- point_metrics(cloud_raw, ~metric_geometry_features(X,Y,Z), k = 20)

However, this is still very slow for my small test cloud (14 million points take about 15 minutes). While calculating, it tells mit it only uses one thread and I don't know how to change this and if that would improve the performance a lot. Anyway, I don't know how to make this process faster. I know there has to be a way, since when I use the software "CloudCompare", the processing is much faster and is done within seconds. I also tried to use stdshapemetrics() from the lidR package which uses an inaccessible function fast_eigen_values(). This needs equally long and tells me it uses only one thread.

Answer 1

The problem of point_metrics() is that it calls user's R code millions of times and this have a cost. Moreover it cannot be safely multithreaded. The function is good for prototyping but for production you must write your own code. For example you can reproduce the function segment_shape() with point_metrics() but segment_shape() is pure C++ and multi-threaded and is often an order of magnitude faster.

Let try with ~3 millions points. The two examples are not equivalent (different output) but the computation load is almost the same (eigen value decomposition).

system.time({point_metrics(cloud_raw, .stdshapemetrics, k = 20)})
#> 110 seconds
set_lidr_threads(4L)
system.time({segment_shapes(cloud_raw, shp_plane(k = 20))})
#> 17 seconds

Also you may be advised to use the adequate function readLAS() as a function of your point-cloud. The previous examples, if read with readTLSLAS(), take 190 and 50 seconds respectively because my point-cloud is ALS. But if you are working with TLS you are better to use readTLSLAS().

Also .stdshapemetrics is nothing else than a wrapper around the C++ you already found. No gain expected.

As a conclusion:

• Write your own C++ code. Take inspiration from lidR source code of segment_shape
• Use lidR spatial index C++ API to reprocuce lidR code.
• Ask lidR questions on https://gis.stackexchange.com/questions/tagged/lidr with the tag lidr

Or go on the lidR repo and ask for a new feature of a native fast eigen value decomposition.