principal component coefficient #

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)
 

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

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)
 }
+

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;
 }
 

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]; }

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},

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);
 }