eigen vectors (normal vector alias PC3)

[DijoG]

Dear lidR experts, dear Jean Romain,

could there be any way to speed up the computation of the following code that tries to extract eigen vectors (I need the normal vector ~ PC3 ~ on the centroid) and eigen values. I know there is a much faster version in C++ of getting the eigen values, and there is the 4-threaded point_eigenvalues() function, however they do not supply the eigen vectors.
It has been running for more than 2 hours and has just reached 70%:

get_eigen = function(x, y, z) {
  xyz = data.table::data.table(x = x,
                               y = y,
                               z = z)
  pca = prcomp(xyz)
  return(list(eigen_largest        = pca$sdev[1],
              eigen_medium         = pca$sdev[2],
              eigen_smallest       = pca$sdev[3],
              eigen_largest_vecX   = pca$rotation["x", "PC1"],
              eigen_largest_vecY   = pca$rotation["y", "PC1"],
              eigen_largest_vecZ   = pca$rotation["z", "PC1"],
              eigen_medium_vecx    = pca$rotation["x", "PC2"],
              eigen_medium_vecY    = pca$rotation["y", "PC2"],
              eigen_medium_vecZ    = pca$rotation["z", "PC2"],
              eigen_smallest_vecX  = pca$rotation["x", "PC3"],
              eigen_smallest_vecY  = pca$rotation["y", "PC3"],
              eigen_smallest_vecZ  = pca$rotation["z", "PC3"]))
}

start_time = Sys.time()
M = point_metrics(l, ~ get_eigen(X, Y, Z), k = 12, xyz = TRUE)
end_time = Sys.time()
end_time - start_time

Thank you in advance and
best regards,
Gergo Dioszegi

[Jean-Romain]

Your problem is known and documented. In the documentation of point_metrics there is a section performance that explains why it is computationally demanding. This section is not exhaustive and it misses to explain the additional overhead of spatial queries and the overhead of mapping an R function at C++ level. In summary the function allows to do anything because it is not specialized but this has a cost. There are two ways to improve the performances.

In the documetation there is an example with the following comment

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

I don't know if prcomp is similar to cov and eigen but there is maybe a factor 10 to gain with a simple C++ implementation of the eigen decomposition.

There is another factor 10 to gain by coding the function in pure C++. I provided one example on gis.stackexchange here. It does not solve your issue but shows you all the tools you can used from lidR to resolve your problem.

point_metrics() is designed for prototyping. The complexity of the operation implies that, once you have a valid proof of concept, you should convert to pure C++ if you want reasonable performances.

[DijoG]

So it is Rcpp and RcppArmadillo. Thank you, everything is clear!

[Jean-Romain]

45 seconds to 5 seconds (9 folds speed-up) using armadillo

(side node: pca$sdev does not contains eigen values. check the doc)

Rcpp::sourceCpp(code = "
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
Rcpp::List cpp_prcomp(arma::mat A) {
arma::mat coeff;
arma::mat score;
arma::vec latent;
arma::princomp(coeff, score, latent, A);
Rcpp::NumericVector eig = Rcpp::wrap(latent);
Rcpp::NumericMatrix rot = Rcpp::wrap(coeff);
return Rcpp::List::create(eig, rot);
}")

get_eigen = function(x, y, z) {
  xyz = cbind(x, y, z)
  pca = cpp_prcomp(xyz)
  return(list(eigen_largest        = pca[[1]][1,],
              eigen_medium         = pca[[1]][2,],
              eigen_smallest       = pca[[1]][3,],
              eigen_largest_vecX   = pca[[2]][1,1],
              eigen_largest_vecY   = pca[[2]][2,1],
              eigen_largest_vecZ   = pca[[2]][3,1],
              eigen_medium_vecx    = pca[[2]][1,2],
              eigen_medium_vecY    = pca[[2]][2,2],
              eigen_medium_vecZ    = pca[[2]][3,2],
              eigen_smallest_vecX  = pca[[2]][1,3],
              eigen_smallest_vecY  = pca[[2]][2,3],
              eigen_smallest_vecZ  = pca[[2]][3,3]))
}

library(lidR)
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile)
start_time = Sys.time()
M = point_metrics(las, ~ get_eigen(X, Y, Z), k = 12, xyz = TRUE)
end_time = Sys.time()
end_time - start_time
#> Time difference of 5.266334 secs

[Jean-Romain]

5 seconds to 0.5 seconds (10 folds speed-up) using pure C++ copying the example given on gis.stackexchange.

This gives us a 90-100x speed-up compared to your example. You problem is therefore solvable in something like 2 minutes

// [[Rcpp::depends(lidR)]]
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::plugins(openmp)]]

#include <RcppArmadillo.h>
#include <SpatialIndex.h>

using namespace Rcpp;
using namespace lidR;

// [[Rcpp::export]]
NumericMatrix eigen_vector(S4 las, int k, int ncpu = 1)
{
  DataFrame data = as<DataFrame>(las.slot("data"));
  NumericVector X = data["X"];
  NumericVector Y = data["Y"];
  NumericVector Z = data["Z"];
  int npoints = X.size();

  NumericMatrix out(npoints, 12);

  SpatialIndex index(las);

  #pragma omp parallel for num_threads(ncpu)
  for (unsigned int i = 0 ; i < npoints ; i++)
  {
    arma::mat A(k,3);
    arma::mat coeff;  // Principle component matrix
    arma::mat score;
    arma::vec latent; // Eigenvalues in descending order

    PointXYZ p(X[i], Y[i], Z[i]);

    std::vector<PointXYZ> pts;
    index.knn(p, k, pts);

    for (unsigned int j = 0 ; j < pts.size() ; j++)
    {
      A(j,0) = pts[j].x;
      A(j,1) = pts[j].y;
      A(j,2) = pts[j].z;
    }

    arma::princomp(coeff, score, latent, A);

    #pragma omp critical
    {
      out(i, 0) = latent[0];
      out(i, 1) = latent[1];
      out(i, 2) = latent[2];
      out(i, 3) = coeff(0,0);
      out(i, 4) = coeff(1,0);
      out(i, 5) = coeff(2,0);
      out(i, 6) = coeff(1,0);
      out(i, 7) = coeff(1,1);
      out(i, 8) = coeff(1,2);
      out(i, 9) = coeff(2,0);
      out(i, 10) = coeff(2,1);
      out(i, 11) = coeff(2,2);
    }
  }

  return out;
}

library(lidR)
LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile)
start_time = Sys.time()
M = eigen_vector(las, k = 12, ncpu = 4)
end_time = Sys.time()
end_time - start_time
#> Time difference of 0.5432258 secs

[Jean-Romain]

I also want to implement radius based (sphere) selection in order to acquire eigenvalues, eigenvectors and centroids

Change

PointXYZ p(X[i], Y[i], Z[i]);
std::vector<PointXYZ> pts;
index.knn(p, k, pts);

to

Sphere s(X[i], Y[i], Z[i], radius)
std::vector<PointXYZ> pts;
index.lookup(s, pts)

Do not help me, it will a good challenge for me to implement it in C++

The C++ API is not really documented. Unless you read the source code you cannot guess it.

[DijoG]

Well, thank you! I have tried it, left k in C++ code.
If my understanding of the code is correct:
the dimension of the matrix has to be set by k (k needs to be a number of points that has an outer sphere boundary where the radius of the sphere defined by k > radius given by radius), afterwards with Sphere s(X[i], Y[i], Z[i], radius) smaller Sphere is set. The points of the smaller sphere are selected from the A matrix in the loop of for (unsigned int j = 0 ; j < pts.size() ; j++)??
Am I right??

NumericMatrix eigen_vector(**S4 las, float radius, int k, int ncpu = 1)
{...
    arma::mat A(k,3);
    arma::mat coeff;  // Principle component matrix
    arma::mat score;
    arma::vec latent; // Eigenvalues in descending order

Sphere s(X[i], Y[i], Z[i], radius);

    std::vector<PointXYZ> pts;
    index.lookup(s, pts);

for (unsigned int j = 0 ; j < pts.size() ; j++)
    {
      A(j,0) = pts[j].x;
      A(j,1) = pts[j].y;
      A(j,2) = pts[j].z;
    }
...

What if I want to kick out int k and select the A matrix beforehand?

NumericMatrix eigen_vector(S4 las, float radius, int ncpu = 1)

So, if my logic is correct:

arma::mat A(k,3);   // k should be the number of points selected by the the radius ????

?

[DijoG]

The solution must be embedded in pts.size() in the A matrix dimension setting:

  #pragma omp parallel for num_threads(ncpu)
  for (unsigned int i = 0 ; i < npoints ; i++)
  {
    Sphere s(X[i], Y[i], Z[i], radius);
    
    std::vector<PointXYZ> pts;
    index.lookup(s, pts);
    
    arma::mat A(pts.size(), 3);  // **_pts.size()_**
    arma::mat coeff;  // Principle component matrix
    arma::mat score;
    arma::vec latent; // Eigenvalues in descending order

    for (unsigned int j = 0 ; j < pts.size() ; j++)
    {
      A(j,0) = pts[j].x;
      A(j,1) = pts[j].y;
      A(j,2) = pts[j].z;
    }

[Jean-Romain]

You are correct.

[DijoG]

Dear Jean-Romain,

your guidance is highly appreciated and I could not have been able to solve the issue without it!

If someone needs the pure C++ code for the eigenvalues, -vectors and centroid extraction (the radius ~ sphere point selection version), here it is:

// [[Rcpp::depends(lidR)]]
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::plugins(openmp)]]

#include <RcppArmadillo.h>
#include <SpatialIndex.h>

using namespace Rcpp;
using namespace lidR;

// [[Rcpp::export]]
NumericMatrix get_eigen(S4 las, float radius, int ncpu = 1)
{
  DataFrame data = as<DataFrame>(las.slot("data"));
  NumericVector X = data["X"];
  NumericVector Y = data["Y"];
  NumericVector Z = data["Z"];
  int npoints = X.size();

  NumericMatrix out(npoints, 15);

  SpatialIndex index(las);

  #pragma omp parallel for num_threads(ncpu)
  for (unsigned int i = 0 ; i < npoints ; i++)
  {
    
    Sphere s(X[i], Y[i], Z[i], radius);
    
    std::vector<PointXYZ> pts;
    index.lookup(s, pts);
    
    arma::mat A(pts.size(), 3);
    arma::mat coeff;      // Principle component matrix
    arma::mat score;
    arma::vec latent;     // Eigenvalues in descending order
       

    for (unsigned int j = 0 ; j < pts.size() ; j++)
    {
      A(j,0) = pts[j].x;
      A(j,1) = pts[j].y;
      A(j,2) = pts[j].z;
    }

    arma::princomp(coeff, score, latent, A);
    arma::mat centro   = arma::mean(A);                                           // Centroid 
    arma::vec centroid = arma::conv_to<arma::vec>::from(centro);  // Centroid vector

    #pragma omp critical
    {
      out(i, 0) = latent[0];
      out(i, 1) = latent[1];
      out(i, 2) = latent[2];
      out(i, 3) = coeff(0,0);
      out(i, 4) = coeff(1,0);
      out(i, 5) = coeff(2,0);
      out(i, 6) = coeff(1,0);
      out(i, 7) = coeff(1,1);
      out(i, 8) = coeff(1,2);
      out(i, 9) = coeff(2,0);
      out(i, 10) = coeff(2,1);
      out(i, 11) = coeff(2,2);
      out(i, 12) = centroid[0];
      out(i, 13) = centroid[1];
      out(i, 14) = centroid[2];
    }
  }

  return out;
}

Thank you very much and
best regards,
Gergo