Get eigenvalues of large point cloud using lidR

Question

I have some very big TLS point clouds (up to 400 million points) and want to assign the eigenvalues to each point in order to later calculate geometric features like planarity or sphericity. I use R and the package lidR following lidR book guidances. I realized that the function point_metrics() in combination with fast_eigen_values() is too slow. It takes ~15 minutes to compute

# 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)

The function segment_shape() is way faster and also computes eigenvalues in the background. However, looking at the C++ Code behind lidR - although I can't write C++ - it seems to be that it is written in a way that it only able to return booleans, not doubles.

Does someone know how to solve this problem? Was there already a version somewhere implemented, which enables me to use the spatial indexing to return eigenvalues without needing to learn C++ just for this? My next resort would be to try using voxel_metrics(). However, I would prefer to have one set of eigenvalues for each point.

Answer 1

Your question is complex but luckily has an answer. For future readers I copy here (parts of) the answer I gave you on SO but that is lost in the sea in absence of appropriated tags on SO.

The problem of point_metrics() is that it calls user's R code millions of times making back and forth between R and C++. This have a cost. Moreover it cannot be safely multi-threaded. 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. In absence of native C++ function in lidR you must write your own C++ code.

The good new is that I'm currently (2021-05-06) working on a native C++ eigen_decomposition() function that will likely be added in lidR 3.2.0. The other good new is that R and lidR provide every tools you need to create your own C++ function exactly like the one we can create in lidR. In the following is a simplified (and standalone) version of the eigen_decomposition() C++ function I'm currently working on. Some features are missing such as progress estimation, user abortion, on-the-fly subset computation, max radius search and so on.

I'm not going through details. Rcpp has its own book, RcppArmadillo is for linear algebra and lidR spatial indexing C++ API is described in the book, opemmp is for multi-threading. You can create a file named eigen_decomposition.cpp:

// [[Rcpp::depends(lidR)]]
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::plugins(openmp)]]

#include <RcppArmadillo.h>
#include <SpatialIndex.h>

using namespace Rcpp;
using namespace lidR;

// [[Rcpp::export]]
NumericMatrix eigen_decomposition(S4 las, int k, int ncpu = 1)
{
  DataFrame data = as<DataFrame>(las.slot("data"));
  NumericVector X = data["X"];
  NumericVector Y = data["Y"];
  NumericVector Z = data["Z"];
  int npoints = X.size();

  NumericMatrix out(npoints, 3);

  SpatialIndex index(las);

  #pragma omp parallel for num_threads(ncpu)
  for (unsigned int i = 0 ; i < npoints ; i++)
  {
    arma::mat A(k,3);
    arma::mat coeff;  // Principle component matrix
    arma::mat score;
    arma::vec latent; // Eigenvalues in descending order

    PointXYZ p(X[i], Y[i], Z[i]);

    std::vector<PointXYZ> pts;
    index.knn(p, k, pts);

    for (unsigned int j = 0 ; j < pts.size() ; j++)
    {
      A(j,0) = pts[j].x;
      A(j,1) = pts[j].y;
      A(j,2) = pts[j].z;
    }

    arma::princomp(coeff, score, latent, A);

    #pragma omp critical
    {
      out(i, 0) = latent[0];
      out(i, 1) = latent[1];
      out(i, 2) = latent[2];
    }
  }

  return out;
}

Then in R

 Rcpp::sourceCpp('eigen_decomposition.cpp')

Lets benchmark it on ALS data. We can see with 4 cores it is 10-fold faster.

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las <- readLAS(LASfile)

microbenchmark::microbenchmark(
u1 = eigen_decomposition(las, 10, 1),
u2 = eigen_decomposition(las, 10, 4),
u3 = point_metrics(las, .stdshapemetrics, k = 10),
times = 3)
#> expr       min        lq      mean    median        uq       max neval
#>   u1  749.5519  755.0064  780.7645  760.4609  796.3708  832.2807     3
#>   u2  339.2228  339.8316  342.9141  340.4403  344.7597  349.0792     3
#>   u3 3173.9732 3188.8633 3370.6251 3203.7533 3468.9511 3734.1488     3

Last but not least, your are working with TLS not ALS. Again it is important to read the documentation(help("lidR-spatial-index"), book chapter). To make it simple you can use readTLSLAS() to say "hey! it is TLS data, please use appropriated spatial index." Lets benchmark it. We found a 7-fold speed-up. Not bad. It could be more, it it could be less, it depends on the point cloud.

file <- system.file("extdata", "pine_plot.laz", package="TreeLS")
als <- readLAS(file, select='xyz')
tls <- readTLSLAS(file, select='xyz')

microbenchmark::microbenchmark(
  u1 = eigen_decomposition(als, 10, 4),
  u2 = eigen_decomposition(tls, 10, 4),
  u3 = point_metrics(als, .stdshapemetrics, k = 10),
  u4 = point_metrics(tls, .stdshapemetrics, k = 10),
  times = 3)
#> expr       min       lq      mean   median       uq      max neval
#>   u1 1373.9824 1394.448 1462.4316 1414.914 1506.656 1598.398     3
#>   u2  693.2717  864.446  929.9681 1035.620 1048.316 1061.012     3
#>   u3 6441.4701 6564.546 6971.8443 6687.621 7237.031 7786.441     3
#>   u4 4792.9041 5118.473 5444.4808 5444.041 5770.269 6096.497     3

To finish lets check that both eigen_decomposition() and .stdshapemetrics return the same (not exactly because .stdshapemetrics returns more metrics).

as.numeric(u2[1,2:4])
#> [1] 8.186695e-04 6.370588e-04 3.090013e-05
u1[1,]
#> [1] 8.186695e-04 6.370588e-04 3.090013e-05