r - 使用 R 尽快计算大型点云的特征值

Question

我有大型点云(每个点云有几百万个点)并且想要计算所有点的几何属性。为此，我需要计算特征值。我怎样才能尽快做到这一点？我的数据存储在 *.las 文件中，我使用 lidR 包读取它们。我还使用这个包来计算点指标。根据this帖子，我实现了这个版本:

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

但是，这对于我的小型测试云来说仍然非常慢(1400 万个点大约需要 15 分钟)。在计算时，它告诉 mit 它只使用一个线程，我不知道如何更改它以及这是否会大大提高性能。无论如何，我不知道如何使这个过程更快。我知道一定有办法，因为当我使用“CloudCompare”软件时，处理速度要快得多，并且在几秒钟内完成。我还尝试使用lidR包中的stdshapemetrics()，该包使用不可访问的函数fast_eigen_values()。这需要同样长的时间，并告诉我它只使用一个线程。

Answer 1

point_metrics() 的问题是它调用用户的 R 代码数百万次，并且这是有成本的。而且它不能安全地进行多线程。该函数非常适合原型(prototype)设计，但对于生产，您必须编写自己的代码。例如，您可以使用 point_metrics() 重现函数 segment_shape()，但 segment_shape() 是纯 C++ 且是多线程的，并且通常是速度快一个数量级。

让我们尝试一下大约 300 万点。这两个示例并不等效(输出不同)，但计算量几乎相同(特征值分解)。

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

此外，建议您使用适当的函数 readLAS() 作为点云的函数。如果使用 readTLSLAS() 读取前面的示例，则分别需要 190 秒和 50 秒，因为我的点云是 ALS。但如果您使用 TLS，则最好使用 readTLSLAS()。

此外，.stdshapemetrics只不过是您已经找到的 C++ 的包装器。预计不会有 yield 。

结论:

• 编写您自己的 C++ 代码。从 segment_shape 的 lidR 源代码中获取灵感
• 使用lidR spatial index C++ API重新生成lidR代码。
• 向lidR询问https://gis.stackexchange.com/questions/tagged/lidr的问题带有标签 lidr

或者继续lidR repo并要求 native 快速特征值分解的新功能。