用 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 代码数百万次，这是有代价的。此外，它不能安全地多线程。该功能适用​​于原型设计，但对于生产，您必须编写自己的代码。例如，您可以使用 point_metrics() 重现函数 segment_shape()，但 segment_shape() 是纯 C++ 和多线程，并且通常快一个数量级。

让我们尝试大约 3 百万点。这两个例子不等价（输出不同），但计算量几乎相同（特征值分解）。

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++ 的包装器。预计没有收益。

作为结论：

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

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