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.