if (FALSE) {
LASfile <- system.file("extdata", "Topography.laz", package="lidR")
# Read only 0.5 points/m^2 for the purposes of this example
las = readLAS(LASfile, filter = "-thin_with_grid 2")
# Computes the eigenvalues of the covariance matrix of the neighbouring
# points and applies a test on these values. This function simulates the
# 'shp_plane()' algorithm from 'segment_shape()'
plane_metrics1 = function(x,y,z, th1 = 25, th2 = 6) {
xyz <- cbind(x,y,z)
cov_m <- cov(xyz)
eigen_m <- eigen(cov_m)$value
is_planar <- eigen_m[2] > (th1*eigen_m[3]) && (th2*eigen_m[2]) > eigen_m[1]
return(list(planar = is_planar))
}
# Apply a user-defined function
M <- point_metrics(las, ~plane_metrics1(X,Y,Z), k = 25)
#> Computed in 6.3 seconds
# We can verify that it returns the same as 'shp_plane'
las <- segment_shapes(las, shp_plane(k = 25), "planar")
#> Computed in 0.1 seconds
all.equal(M$planar, las$planar)
# At this stage we can be clever and find that the bottleneck is
# the eigenvalue computation. Let's write a C++ version of it with
# Rcpp and RcppArmadillo
Rcpp::sourceCpp(code = "
#include
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
SEXP eigen_values(arma::mat A) {
arma::mat coeff;
arma::mat score;
arma::vec latent;
arma::princomp(coeff, score, latent, A);
return(Rcpp::wrap(latent));
}")
plane_metrics2 = function(x,y,z, th1 = 25, th2 = 6) {
xyz <- cbind(x,y,z)
eigen_m <- eigen_values(xyz)
is_planar <- eigen_m[2] > (th1*eigen_m[3]) && (th2*eigen_m[2]) > eigen_m[1]
return(list(planar = is_planar))
}
M <- point_metrics(las, ~plane_metrics2(X,Y,Z), k = 25)
#> Computed in 0.5 seconds
all.equal(M$planar, las$planar)
# Here we can see that the optimized version is way better but is still 5-times slower
# because of the overhead of calling R functions and switching back and forth from R to C++.
M <- point_eigenvalues(las, k = 25)
is_planar = M$eigen_medium > (25*M$eigen_smallest) & (6*M$eigen_medium) > M$eigen_largest
# Use the filter argument to process only first returns
M1 <- point_metrics(las, ~plane_metrics2(X,Y,Z), k = 25, filter = ~ReturnNumber == 1)
dim(M1) # 13894 instead of 17182 previously.
}
Run the code above in your browser using DataLab