使用 R 尽快计算大点云的特征值
Compute eigenvalues for large point clouds as quickly as possible with R
我有大量的点云(每个点云有数百万个点)并且想要计算所有点的几何属性。为此,我需要计算特征值。我怎样才能尽快做到这一点?我的数据存储在 *.las 文件中,我用包 lidR 读取它们。我还使用这个包来计算点指标。根据thispost,我实现了这个版本:
# 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()
。这需要同样长的时间并告诉我它只使用一个线程。
point_metrics()
的问题是它调用用户的 R 代码数百万次,这是有代价的。此外,它不能安全地多线程。该功能适用于原型制作,但对于生产,您必须编写自己的代码。例如,您可以使用 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++ 的包装器。预期没有收益。
作为结论:
- 编写您自己的 C++ 代码。从 segment_shape
的 lidR
源代码中获取灵感
- 使用 lidR spatial index C++ API 重新生成 lidR 代码。
- 在 https://gis.stackexchange.com/questions/tagged/lidr 上使用标签
lidr
问 lidR 问题
或者继续 lidR repo 并要求新特性的原生快速特征值分解。
我有大量的点云(每个点云有数百万个点)并且想要计算所有点的几何属性。为此,我需要计算特征值。我怎样才能尽快做到这一点?我的数据存储在 *.las 文件中,我用包 lidR 读取它们。我还使用这个包来计算点指标。根据thispost,我实现了这个版本:
# 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()
。这需要同样长的时间并告诉我它只使用一个线程。
point_metrics()
的问题是它调用用户的 R 代码数百万次,这是有代价的。此外,它不能安全地多线程。该功能适用于原型制作,但对于生产,您必须编写自己的代码。例如,您可以使用 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++ 的包装器。预期没有收益。
作为结论:
- 编写您自己的 C++ 代码。从 segment_shape 的
- 使用 lidR spatial index C++ API 重新生成 lidR 代码。
- 在 https://gis.stackexchange.com/questions/tagged/lidr 上使用标签
lidr
问 lidR 问题
lidR
源代码中获取灵感
或者继续 lidR repo 并要求新特性的原生快速特征值分解。