Merge pull request #822 from r-lidar/devel

Devel

Author: Jean-Romain

DESCRIPTION

@@ -2,7 +2,7 @@ Package: lidR
 Type: Package
 Title: Airborne LiDAR Data Manipulation and Visualization for Forestry
     Applications
-Version: 4.2.0
+Version: 4.2.2
 Authors@R: c(
     person("Jean-Romain", "Roussel", email = "[email protected]", role = c("aut", "cre", "cph")),
     person("David", "Auty", email = "", role = c("aut", "ctb"), comment = "Reviews the documentation"),
@@ -56,7 +56,7 @@ Suggests:
     knitr,
     rmarkdown
 RoxygenNote: 7.3.2
-LinkingTo: BH (>= 1.72.0),Rcpp,RcppArmadillo
+LinkingTo: BH (>= 1.72.0),Rcpp,RcppArmadillo,Rnanoflann
 Encoding: UTF-8
 ByteCompile: true
 VignetteBuilder: knitr

NEWS.md

@@ -1,6 +1,14 @@
 If you are viewing this file on CRAN, please check [the latest news on GitHub](https://github.com/r-lidar/lidR/blob/master/NEWS.md) where the formatting is also better
 
-## lidR v4.2.0 (Release date: )
+## lidR v4.2.2 (Release date: 2025-06-02)
+
+- Use `nanoflann` for knn computation
+
+## lidR v4.2.1 (Release date: 2025-06-02)
+
+- Fix compilation on MacOS
+
+## lidR v4.2.0 (Release date: 2025-05-25)
 
 **v4.2.0 brings new tools for terrestrial data (TLS, MLS):**
 

R/RcppExports.R

@@ -113,6 +113,10 @@ C_eigen_metrics <- function(las, k, r, coeffs, filter, ncpu) {
     .Call(`_lidR_C_eigen_metrics`, las, k, r, coeffs, filter, ncpu)
 }
 
+C_fast_knn_eigen_decomposition <- function(las, k, coeffs, ncpu) {
+    .Call(`_lidR_C_fast_knn_eigen_decomposition`, las, k, coeffs, ncpu)
+}
+
 C_connected_component <- function(las, res) {
     .Call(`_lidR_C_connected_component`, las, res)
 }

R/st_crs.R

@@ -269,7 +269,7 @@ setGeneric("epsg<-", function(object, value) standardGeneric("epsg<-"))
 #' @rdname st_crs
 setMethod("epsg", "LASheader", function(object, ...)
 {
-  if (use_wktcs(object))  warning("This LAS object stores the CRS as WKT. 0 returned; use 'wkt()' instead.", call. = FALSE)
+  if (use_wktcs(object))  warning("This LAS object stores the CRS as WKT. CRS field might not be correctly populated, yielding uncertain results; use 'wkt()' instead.", call. = FALSE)
   return(rlas::header_get_epsg(as.list(object)))
 })
 

src/LAS.cpp

@@ -2,10 +2,35 @@
 #include "Progress.h"
 #include "myomp.h"
 #include "SpatialIndex.h"
+#include "nanoflann.hpp"
 #include <limits>
 
 using namespace lidR;
 
+class DataFrameAdaptor
+{
+public:
+  std::vector<Rcpp::NumericVector> coords;
+  size_t dim;
+
+  DataFrameAdaptor(const Rcpp::DataFrame& df, std::vector<std::string> col_names)
+  {
+    dim = col_names.size();
+    coords.reserve(dim);
+    for (const auto& name : col_names)
+      coords.push_back(df[name]);
+  }
+
+  inline size_t kdtree_get_point_count() const { return coords[0].size(); }
+  inline double kdtree_get_pt(const size_t idx, const size_t d) const {
+    return coords[d][idx];
+  }
+  template <class BBOX> bool kdtree_get_bbox(BBOX&) const { return false; }
+};
+
+using KDTree = nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Simple_Adaptor<double, DataFrameAdaptor>, DataFrameAdaptor, 3>;
+
+
 bool pnpoly(NumericMatrix polygon, double x, double y)
 {
   int nvert = polygon.nrow();
@@ -1547,6 +1572,112 @@ List LAS::point_metrics(unsigned int k, double r, DataFrame data, int nalloc, SE
 }
 #endif
 
+DataFrame LAS::fast_knn_eigen_decomposition(int k, bool get_coef)
+{
+  int n = npoints;
+
+  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;
+
+  DataFrame df = as<DataFrame>(las.slot("data"));
+  DataFrameAdaptor adaptor(df, {"X", "Y", "Z"});
+  KDTree tree = KDTree(3, adaptor, nanoflann::KDTreeSingleIndexAdaptorParams(10));
+  tree.buildIndex();
+
+  Progress pb(npoints, "Eigen decomposition: ");
+
+  #pragma omp parallel for num_threads(ncpu)
+  for (unsigned int i = 0 ; i < npoints ; i++)
+  {
+    if (abort) continue;
+    if (pb.check_interrupt()) abort = true;
+    pb.increment();
+
+    std::vector<uint32_t> indices(k);
+    std::vector<KDTree::DistanceType> dists(k);
+    double p[3] = { X[i], Y[i], Z[i] };
+    tree.knnSearch(p, k, indices.data(), dists.data());
+
+
+    arma::mat A(k,3);
+    arma::mat coeff;  // Principle component matrix
+    arma::mat score;
+    arma::vec latent; // Eigenvalues in descending order
+
+    for (unsigned int j = 0 ; j < indices.size() ; j++)
+    {
+      A(j,0) = X[indices[j]];
+      A(j,1) = Y[indices[j]];
+      A(j,2) = Z[indices[j]];
+    }
+
+    arma::princomp(coeff, score, latent, A);
+
+    eigen_largest[i] = latent[0];
+    eigen_medium[i] = latent[1];
+    eigen_smallest[i] = latent[2];
+
+    if (get_coef)
+    {
+      coeff0[i] = coeff(0,0);
+      coeff1[i] = coeff(0,1);
+      coeff2[i] = coeff(0,2);
+      coeff3[i] = coeff(1,0);
+      coeff4[i] = coeff(1,1);
+      coeff5[i] = coeff(1,2);
+      coeff6[i] = coeff(2,0);
+      coeff7[i] = coeff(2,1);
+      coeff8[i] = coeff(2,2);
+    }
+  }
+
+  DataFrame out;
+  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;
+}
+
+
 DataFrame LAS::eigen_decomposition(int k, double r, bool get_coef)
 {
   int n = std::count(skip.begin(), skip.end(), true);

src/LAS.h

@@ -50,6 +50,7 @@ class LAS
     NumericVector knn_distance(unsigned int k);
     NumericVector interpolate_knnidw(NumericVector x, NumericVector y, int k, double p, double rmax);
     DataFrame eigen_decomposition(int k, double r, bool get_coef);
+    DataFrame fast_knn_eigen_decomposition(int k, 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

@@ -385,6 +385,19 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// C_fast_knn_eigen_decomposition
+DataFrame C_fast_knn_eigen_decomposition(S4 las, int k, bool coeffs, int ncpu);
+RcppExport SEXP _lidR_C_fast_knn_eigen_decomposition(SEXP lasSEXP, SEXP kSEXP, SEXP coeffsSEXP, 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< bool >::type coeffs(coeffsSEXP);
+    Rcpp::traits::input_parameter< int >::type ncpu(ncpuSEXP);
+    rcpp_result_gen = Rcpp::wrap(C_fast_knn_eigen_decomposition(las, k, coeffs, ncpu));
+    return rcpp_result_gen;
+END_RCPP
+}
 // C_connected_component
 IntegerVector C_connected_component(S4 las, double res);
 RcppExport SEXP _lidR_C_connected_component(SEXP lasSEXP, SEXP resSEXP) {
@@ -690,6 +703,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_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, 6},
+    {"_lidR_C_fast_knn_eigen_decomposition", (DL_FUNC) &_lidR_C_fast_knn_eigen_decomposition, 4},
     {"_lidR_C_connected_component", (DL_FUNC) &_lidR_C_connected_component, 2},
     {"_lidR_C_voxel_id", (DL_FUNC) &_lidR_C_voxel_id, 2},
     {"_lidR_fast_table", (DL_FUNC) &_lidR_fast_table, 2},

src/RcppFunction.cpp

@@ -236,6 +236,13 @@ DataFrame C_eigen_metrics(S4 las, int k, double r, bool coeffs, LogicalVector fi
   return pt.eigen_decomposition(k, r, coeffs);
 }
 
+//[[Rcpp::export(rng = false)]]
+DataFrame C_fast_knn_eigen_decomposition(S4 las, int k, bool coeffs, int ncpu)
+{
+  LAS pt(las, ncpu);
+  return pt.fast_knn_eigen_decomposition(k, coeffs);
+}
+
 
 #include "lidR/Grid3D.h"
 

src/knn.cpp

@@ -1,21 +1,46 @@
-#include "SpatialIndex.h"
+#include "nanoflann.hpp"
 #include "Progress.h"
 
 #include <Rcpp.h>
 #include "myomp.h"
 
 using namespace Rcpp;
-using namespace lidR;
+
+class DataFrameAdaptor
+{
+public:
+  std::vector<Rcpp::NumericVector> coords;
+  size_t dim;
+
+  DataFrameAdaptor(const Rcpp::DataFrame& df, std::vector<std::string> col_names)
+  {
+    dim = col_names.size();
+    coords.reserve(dim);
+    for (const auto& name : col_names)
+      coords.push_back(df[name]);
+  }
+
+  inline size_t kdtree_get_point_count() const { return coords[0].size(); }
+  inline double kdtree_get_pt(const size_t idx, const size_t d) const {
+    return coords[d][idx];
+  }
+  template <class BBOX> bool kdtree_get_bbox(BBOX&) const { return false; }
+};
+
+using KDTree = nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Simple_Adaptor<double, DataFrameAdaptor>, DataFrameAdaptor, 3>;
 
 // [[Rcpp::export]]
 List cpp_knn(S4 data, int k, int ncpu)
 {
-  SpatialIndex tree(data);
 
-  DataFrame tmp = as<DataFrame>(data.slot("data"));
-  NumericVector X = tmp["X"];
-  NumericVector Y = tmp["Y"];
-  NumericVector Z = tmp["Z"];
+  DataFrame df = as<DataFrame>(data.slot("data"));
+  NumericVector X = df["X"];
+  NumericVector Y = df["Y"];
+  NumericVector Z = df["Z"];
+
+  DataFrameAdaptor adaptor(df, {"X", "Y", "Z"});
+  KDTree tree = KDTree(3, adaptor, nanoflann::KDTreeSingleIndexAdaptorParams(10));
+  tree.buildIndex();
 
   int npoints = X.size();
 
@@ -36,19 +61,15 @@ List cpp_knn(S4 data, int k, int ncpu)
     if (pb.check_interrupt()) abort = true;
     pb.increment();
 
-    std::vector<PointXYZ> pts;
-    PointXYZ p(X[i], Y[i], Z[i]);
-    tree.knn(p, k, pts);
+    std::vector<uint32_t> indices(k);
+    std::vector<KDTree::DistanceType> dists(k);
+    double p[3] = { X[i], Y[i], Z[i] };
+    tree.knnSearch(p, k, indices.data(), dists.data());
 
-    for (unsigned int j = 1; j < pts.size(); j++)
+    for (auto j = 1; j < indices.size(); j++)
     {
-      double dx = X[i] - pts[j].x;
-      double dy = Y[i] - pts[j].y;
-      double dz = Z[i] - pts[j].z;
-      double d = std::sqrt(dx*dx+dy*dy+dz*dz);
-
-      knn_idx(i,j-1) = pts[j].id+1;
-      knn_dist(i,j-1) = d;
+      knn_idx(i,j-1) = indices[j]+1;
+      knn_dist(i,j-1) = std::sqrt(dists[j]);
     }
   }
 
@@ -60,7 +81,10 @@ List cpp_knn(S4 data, int k, int ncpu)
 // [[Rcpp::export]]
 List cpp_knnx(S4 data, S4 query, int k, int ncpu)
 {
-  SpatialIndex tree(data);
+  DataFrame temp = as<DataFrame>(data.slot("data"));
+  DataFrameAdaptor adaptor(temp, {"X", "Y", "Z"});
+  KDTree tree = KDTree(3, adaptor, nanoflann::KDTreeSingleIndexAdaptorParams(10));
+  tree.buildIndex();
 
   DataFrame tmp = as<DataFrame>(query.slot("data"));
   NumericVector X = tmp["X"];
@@ -84,19 +108,15 @@ List cpp_knnx(S4 data, S4 query, int k, int ncpu)
     if (pb.check_interrupt()) abort = true;
     pb.increment();
 
-    std::vector<PointXYZ> pts;
-    PointXYZ p(X[i], Y[i], Z[i]);
-    tree.knn(p, k, pts);
+    std::vector<uint32_t> indices(k);
+    std::vector<KDTree::DistanceType> dists(k);
+    double p[3] = { X[i], Y[i], Z[i] };
+    tree.knnSearch(p, k, indices.data(), dists.data());
 
-    for (unsigned int j = 0 ; j < pts.size(); j++)
+    for (auto j = 0; j < indices.size(); j++)
     {
-      double dx = X[i] - pts[j].x;
-      double dy = Y[i] - pts[j].y;
-      double dz = Z[i] - pts[j].z;
-      double d = std::sqrt(dx*dx+dy*dy+dz*dz);
-
-      knn_idx(i,j) = pts[j].id+1;
-      knn_dist(i,j) = d;
+      knn_idx(i,j) = indices[j]+1;
+      knn_dist(i,j) = std::sqrt(dists[j]);
     }
   }