From 3b7d705c188ac47cbdb6c6c69eafca294cfeb308 Mon Sep 17 00:00:00 2001 From: Jean-Romain Date: Wed, 18 Oct 2023 21:27:07 -0400 Subject: [PATCH] principal component coefficient # --- NEWS.md | 1 + R/RcppExports.R | 4 ++-- R/metrics_point.R | 10 ++++---- man/point_metrics.Rd | 17 +++++++++++--- src/LAS.cpp | 54 +++++++++++++++++++++++++++++++++++++++++++- src/LAS.h | 2 +- src/RcppExports.cpp | 9 ++++---- src/RcppFunction.cpp | 4 ++-- 8 files changed, 84 insertions(+), 17 deletions(-) diff --git a/NEWS.md b/NEWS.md index cd057979..1c9644ba 100644 --- a/NEWS.md +++ b/NEWS.md @@ -3,6 +3,7 @@ If you are viewing this file on CRAN, please check [the latest news on GitHub](h ## lidR v4.1.0 (Release date: ) 1. New functions `head()` and `tail()` for `LAS` objects +2. New: `point_eigenvalues` gained an argument `coeff` to return the principal component coefficients ## lidR v4.0.4 (Release date: 2023-09-07) diff --git a/R/RcppExports.R b/R/RcppExports.R index 7920e481..0a625707 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -101,8 +101,8 @@ C_check_gpstime <- function(t, rn) { .Call(`_lidR_C_check_gpstime`, t, rn) } -C_eigen_metrics <- function(las, k, r, filter, ncpu) { - .Call(`_lidR_C_eigen_metrics`, las, k, r, filter, ncpu) +C_eigen_metrics <- function(las, k, r, coeffs, filter, ncpu) { + .Call(`_lidR_C_eigen_metrics`, las, k, r, coeffs, filter, ncpu) } fast_table <- function(x, size = 5L) { diff --git a/R/metrics_point.R b/R/metrics_point.R index 660cfe0f..80291cf3 100644 --- a/R/metrics_point.R +++ b/R/metrics_point.R @@ -18,7 +18,7 @@ #' is still computationally demanding.\cr\cr #' To help users to get an idea of how computationally demanding this function is, let's compare it to #' \link{pixel_metrics}. Assuming we want to apply `mean(Z)` on a 1 km² tile with 1 point/m² -#' with a resolution of 20 m (400 m² cells), then the function `mean is called roughly 2500 +#' with a resolution of 20 m (400 m² cells), then the function `mean` is called roughly 2500 #' times (once per cell). On the contrary, with `point_metrics`, `mean` is called 1000000 #' times (once per point). So the function is expected to be more than 400 times slower in this specific #' case (but it does not provide the same feature).\cr\cr @@ -33,7 +33,8 @@ #' computes with a sphere neighborhood, if k and r are given computes with the knn and a limit on the #' search distance. #' @param xyz logical. Coordinates of each point are returned in addition to each metric. Otherwise an -#' ID refering to each point. +#' ID referring to each point. +#' @param coeffs logical. Principal component coefficients are returned #' @param filter formula of logical predicates. Enables the function to run only on points of interest #' in an optimized way. See examples. #' @param ... unused. @@ -211,7 +212,7 @@ point_metrics.LAS <- function(las, func, k, r, xyz = FALSE, filter = NULL, ...) #' @export #' @rdname point_metrics -point_eigenvalues = function(las, k, r, xyz = FALSE, metrics = FALSE, filter = NULL) +point_eigenvalues = function(las, k, r, xyz = FALSE, metrics = FALSE, coeffs = FALSE, filter = NULL) { pointID <- NULL @@ -248,7 +249,7 @@ point_eigenvalues = function(las, k, r, xyz = FALSE, metrics = FALSE, filter = N assert_is_a_bool(xyz) filter <- parse_filter(las, filter) - M <- C_eigen_metrics(las, k, r, filter, getThreads()) + M <- C_eigen_metrics(las, k, r, coeffs, filter, getThreads()) data.table::setDT(M) data.table::setorder(M, pointID) @@ -277,3 +278,4 @@ point_eigenvalues = function(las, k, r, xyz = FALSE, metrics = FALSE, filter = N return(M) } + diff --git a/man/point_metrics.Rd b/man/point_metrics.Rd index c0783bfe..4c328cee 100644 --- a/man/point_metrics.Rd +++ b/man/point_metrics.Rd @@ -7,7 +7,15 @@ \usage{ point_metrics(las, func, k, r, xyz = FALSE, filter = NULL, ...) -point_eigenvalues(las, k, r, xyz = FALSE, metrics = FALSE, filter = NULL) +point_eigenvalues( + las, + k, + r, + xyz = FALSE, + metrics = FALSE, + coeffs = FALSE, + filter = NULL +) } \arguments{ \item{las}{An object of class LAS} @@ -21,7 +29,7 @@ computes with a sphere neighborhood, if k and r are given computes with the knn search distance.} \item{xyz}{logical. Coordinates of each point are returned in addition to each metric. Otherwise an -ID refering to each point.} +ID referring to each point.} \item{filter}{formula of logical predicates. Enables the function to run only on points of interest in an optimized way. See examples.} @@ -29,6 +37,8 @@ in an optimized way. See examples.} \item{...}{unused.} \item{metrics}{logical. Compute metrics or not} + +\item{coeffs}{logical. Coefficients of the covariance matrix are returned.} } \description{ Computes a series of user-defined descriptive statistics for a LiDAR dataset for each point based @@ -51,7 +61,8 @@ mapping a user-defined function at the point level using optimized memory manage is still computationally demanding.\cr\cr To help users to get an idea of how computationally demanding this function is, let's compare it to \link{pixel_metrics}. Assuming we want to apply \code{mean(Z)} on a 1 km² tile with 1 point/m² -with a resolution of 20 m (400 m² cells), then the function \verb{mean is called roughly 2500 times (once per cell). On the contrary, with }point_metrics\verb{, }mean` is called 1000000 +with a resolution of 20 m (400 m² cells), then the function \code{mean} is called roughly 2500 +times (once per cell). On the contrary, with \code{point_metrics}, \code{mean} is called 1000000 times (once per point). So the function is expected to be more than 400 times slower in this specific case (but it does not provide the same feature).\cr\cr This is why the user-defined function is expected to be well-optimized, otherwise it might drastically diff --git a/src/LAS.cpp b/src/LAS.cpp index 13716cfe..276c9dc7 100644 --- a/src/LAS.cpp +++ b/src/LAS.cpp @@ -1498,13 +1498,36 @@ List LAS::point_metrics(unsigned int k, double r, DataFrame data, int nalloc, SE return output; } -DataFrame LAS::eigen_decomposition(int k, double r) +DataFrame LAS::eigen_decomposition(int k, double r, bool get_coef) { int n = std::count(skip.begin(), skip.end(), true); IntegerVector pointID(n); + NumericVector eigen_largest(n); NumericVector eigen_medium(n); NumericVector eigen_smallest(n); + NumericVector coeff0; + NumericVector coeff1; + NumericVector coeff2; + NumericVector coeff3; + NumericVector coeff4; + NumericVector coeff5; + NumericVector coeff6; + NumericVector coeff7; + NumericVector coeff8; + + if (get_coef) + { + coeff0 = NumericVector(n); + coeff1 = NumericVector(n); + coeff2 = NumericVector(n); + coeff3 = NumericVector(n); + coeff4 = NumericVector(n); + coeff5 = NumericVector(n); + coeff6 = NumericVector(n); + coeff7 = NumericVector(n); + coeff8 = NumericVector(n); + } bool abort = false; unsigned int j = 0; @@ -1579,9 +1602,24 @@ DataFrame LAS::eigen_decomposition(int k, double r) #pragma omp critical { pointID[j] = i; + eigen_largest[j] = latent[0]; eigen_medium[j] = latent[1]; eigen_smallest[j] = latent[2]; + + if (get_coef) + { + coeff0[j] = coeff(0,0); + coeff1[j] = coeff(0,1); + coeff2[j] = coeff(0,2); + coeff3[j] = coeff(1,0); + coeff4[j] = coeff(1,1); + coeff5[j] = coeff(1,2); + coeff6[j] = coeff(2,0); + coeff7[j] = coeff(2,1); + coeff8[j] = coeff(2,2); + } + j++; } } @@ -1591,6 +1629,20 @@ DataFrame LAS::eigen_decomposition(int k, double r) out.push_back(eigen_largest, "eigen_largest"); out.push_back(eigen_medium, "eigen_medium"); out.push_back(eigen_smallest, "eigen_smallest"); + + if (get_coef) + { + out.push_back(coeff0, "coeff00"); + out.push_back(coeff1, "coeff01"); + out.push_back(coeff2, "coeff02"); + out.push_back(coeff3, "coeff10"); + out.push_back(coeff4, "coeff11"); + out.push_back(coeff5, "coeff12"); + out.push_back(coeff6, "coeff20"); + out.push_back(coeff7, "coeff21"); + out.push_back(coeff8, "coeff22"); + } + return out; } diff --git a/src/LAS.h b/src/LAS.h index 1aa1866f..e3d8916f 100644 --- a/src/LAS.h +++ b/src/LAS.h @@ -47,7 +47,7 @@ class LAS List point_metrics(unsigned int k, double r, DataFrame data, int nalloc, SEXP call, SEXP env); NumericVector fast_knn_metrics(unsigned int k, IntegerVector metrics); NumericVector interpolate_knnidw(NumericVector x, NumericVector y, int k, double p, double rmax); - DataFrame eigen_decomposition(int k, double r); + DataFrame eigen_decomposition(int k, double r, bool get_coef); private: static bool coplanar (arma::vec& latent, arma::mat& coeff, NumericVector& th) { return latent[1] > th[0]*latent[2] && th[1]*latent[1] > latent[0]; } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 459e7de4..878863fb 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -345,16 +345,17 @@ BEGIN_RCPP END_RCPP } // C_eigen_metrics -DataFrame C_eigen_metrics(S4 las, int k, double r, LogicalVector filter, int ncpu); -RcppExport SEXP _lidR_C_eigen_metrics(SEXP lasSEXP, SEXP kSEXP, SEXP rSEXP, SEXP filterSEXP, SEXP ncpuSEXP) { +DataFrame C_eigen_metrics(S4 las, int k, double r, bool coeffs, LogicalVector filter, int ncpu); +RcppExport SEXP _lidR_C_eigen_metrics(SEXP lasSEXP, SEXP kSEXP, SEXP rSEXP, SEXP coeffsSEXP, SEXP filterSEXP, SEXP ncpuSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::traits::input_parameter< S4 >::type las(lasSEXP); Rcpp::traits::input_parameter< int >::type k(kSEXP); Rcpp::traits::input_parameter< double >::type r(rSEXP); + Rcpp::traits::input_parameter< bool >::type coeffs(coeffsSEXP); Rcpp::traits::input_parameter< LogicalVector >::type filter(filterSEXP); Rcpp::traits::input_parameter< int >::type ncpu(ncpuSEXP); - rcpp_result_gen = Rcpp::wrap(C_eigen_metrics(las, k, r, filter, ncpu)); + rcpp_result_gen = Rcpp::wrap(C_eigen_metrics(las, k, r, coeffs, filter, ncpu)); return rcpp_result_gen; END_RCPP } @@ -611,7 +612,7 @@ static const R_CallMethodDef CallEntries[] = { {"_lidR_C_local_maximum", (DL_FUNC) &_lidR_C_local_maximum, 4}, {"_lidR_C_isolated_voxel", (DL_FUNC) &_lidR_C_isolated_voxel, 3}, {"_lidR_C_check_gpstime", (DL_FUNC) &_lidR_C_check_gpstime, 2}, - {"_lidR_C_eigen_metrics", (DL_FUNC) &_lidR_C_eigen_metrics, 5}, + {"_lidR_C_eigen_metrics", (DL_FUNC) &_lidR_C_eigen_metrics, 6}, {"_lidR_fast_table", (DL_FUNC) &_lidR_fast_table, 2}, {"_lidR_fast_countequal", (DL_FUNC) &_lidR_fast_countequal, 2}, {"_lidR_fast_countbelow", (DL_FUNC) &_lidR_fast_countbelow, 2}, diff --git a/src/RcppFunction.cpp b/src/RcppFunction.cpp index 8f52caa1..53b7507b 100644 --- a/src/RcppFunction.cpp +++ b/src/RcppFunction.cpp @@ -211,11 +211,11 @@ int C_check_gpstime(NumericVector t, IntegerVector rn) } //[[Rcpp::export(rng = false)]] -DataFrame C_eigen_metrics(S4 las, int k, double r, LogicalVector filter, int ncpu) +DataFrame C_eigen_metrics(S4 las, int k, double r, bool coeffs, LogicalVector filter, int ncpu) { LAS pt(las, ncpu); pt.new_filter(filter); - return pt.eigen_decomposition(k, r); + return pt.eigen_decomposition(k, r, coeffs); }