Fix #732

NEWS.md

@@ -3,6 +3,7 @@ If you are viewing this file on CRAN, please check [the latest news on GitHub](h
 ## lidR v4.0.5 (Release date: )
 
 - Fix: [#726](https://github.com/r-lidar/lidR/issues/726) character palette causes error in plot.
+- Fix: [#732](https://github.com/r-lidar/lidR/issues/732) octree spatial indexes is not working properly. lidR now use voxel partition and no longer support octree until the problem is fixed.
 
 ## lidR v4.0.4 (Release date: 2023-09-07)
 

inst/include/SpatialIndex.h

@@ -54,7 +54,7 @@ inline SpatialIndex::SpatialIndex(const Rcpp::S4 las)
   case GRIDPARTITION:
   case VOXELPARTITION: grid = GridPartition(las); break;
   case QUADTREE: quadtree = QuadTree(las); break;
-  case OCTREE: octree = Octree(las); break;
+  case OCTREE: Rcpp::stop("Error: octree no longer supported."); break; // # nocov
   default: Rcpp::stop("Internal error: spatial index code inccorect."); break; // # nocov
   }
 }
@@ -186,7 +186,7 @@ inline int SpatialIndex::index_selector(const Rcpp::S4 las) const
       if (sensor == UKN || sensor == ALS)
         code = GRIDPARTITION;
       else if (sensor == TLS || sensor == UAV || sensor == DAP)
-        code = OCTREE;
+        code = VOXELPARTITION;
       else
         code = QUADTREE;
     }

inst/include/lidR/GridPartition.h

@@ -184,8 +184,15 @@ template<typename T> void GridPartition::lookup(T& shape, std::vector<PointXYZ>&
   int colmax = std::ceil((xmax - this->xmin) / xres);
   int rowmin = std::floor((this->ymax - ymax) / yres);
   int rowmax = std::ceil((this->ymax - ymin) / yres);
-  int laymin = std::floor((zmin - this->zmin) / zres);
-  int laymax = std::ceil((zmax - this->zmin) / zres);
+
+  int laymin = 0;
+  int laymax = nlayers;
+  if (zmin > this->zmin && zmax < this->zmax)
+  {
+    laymin = std::floor((zmin - this->zmin) / zres);
+    laymax = std::ceil((zmax - this->zmin) / zres);
+  }
+
   int cell;
 
   res.clear();

inst/include/lidR/Octree.h

@@ -390,8 +390,14 @@ template<typename T>  Node::Ocnode* Octree::locate_region(T shape)
   double bbxmax = (shape.xmax - xmin)/(xmax-xmin);
   double bbymin = (shape.ymin - ymin)/(ymax-ymin);
   double bbymax = (shape.ymax - ymin)/(ymax-ymin);
-  double bbzmin = (shape.zmin - zmin)/(zmax-zmin);
-  double bbzmax = (shape.zmax - zmin)/(zmax-zmin);
+
+  double bbzmin = 0;
+  double bbzmax = 1;
+  if (shape.zmin > this->zmin && shape.zmax < this->zmax)
+  {
+    bbzmin = (shape.zmin - zmin)/(zmax-zmin);
+    bbzmax = (shape.zmax - zmin)/(zmax-zmin);
+  }
 
   if (bbxmax < 0 || bbxmin > 1 || bbymax < 0 || bbymin > 1 || bbzmax < 0 || bbzmin > 1)
     return 0;

tests/testthat/test-spatialindex.R

@@ -162,7 +162,7 @@ las = LAS(D)
 
 test_that("circle lookup works for each spatial index", {
 
-  for(index in 0:4)
+  for(index in 0:3)
   {
     index(las) <- index
 
@@ -183,7 +183,7 @@ test_that("circle lookup works for each spatial index", {
 
 test_that("oriented rectangle lookup works for each spatial index", {
 
-  for(index in 0:4)
+  for(index in 0:3)
   {
     index(las) <- index
 
@@ -202,7 +202,7 @@ test_that("oriented rectangle lookup works for each spatial index", {
 
 test_that("knn in 2d works for each spatial indexes", {
 
-  for(index in 0:4)
+  for(index in 0:3)
   {
     index(las) <- index
 
@@ -219,7 +219,7 @@ test_that("knn in 2d works for each spatial indexes", {
 
 test_that("knn in 3d works for each spatial indexes", {
 
-  for(index in 0:4)
+  for(index in 0:3)
   {
     index(las) <- index
 
@@ -248,8 +248,8 @@ test_that("Spatial indexes work with more points (coverage)", {
   y = runif(1, 10, 900)
   z = runif(1, 1, 19)
 
-  u = vector("list", 4)
-  for(index in 0:4)
+  u = vector("list", 3)
+  for(index in 0:3)
   {
     index(las) <- index
     id = lidR:::C_knn3d_lookup(las, x,y,z,10)
@@ -258,8 +258,8 @@ test_that("Spatial indexes work with more points (coverage)", {
 
   expect_true(all(sapply(u, identical, u[[1]])))
 
-  u = vector("list", 4)
-  for(index in 0:4)
+  u = vector("list", 3)
+  for(index in 0:3)
   {
     index(las) <- index
     id = lidR:::C_circle_lookup(las, x, y, 5)
@@ -278,8 +278,8 @@ test_that("Spatial indexes work with 0 point", {
   y = runif(1, 10, 900)
   z = runif(1, 1, 19)
 
-  u = vector("list", 4)
-  for(index in 0:4)
+  u = vector("list", 3)
+  for(index in 0:3)
   {
     index(las) <- index
     id = lidR:::C_knn3d_lookup(las, x, y, z,10)
@@ -288,8 +288,8 @@ test_that("Spatial indexes work with 0 point", {
 
   expect_true(all(sapply(u, identical, integer(0))))
 
-  u = vector("list", 4)
-  for(index in 0:4)
+  u = vector("list", 3)
+  for(index in 0:3)
   {
     index(las) <- index
     id = lidR:::C_circle_lookup(las, x, y, 5)
@@ -308,8 +308,8 @@ test_that("Spatial indexes work with 1 point", {
   y = runif(1, 10, 11)
   z = runif(1, 1, 19)
 
-  u = vector("list", 4)
-  for(index in 0:4)
+  u = vector("list", 3)
+  for(index in 0:3)
   {
     index(las) <- index
     id = lidR:::C_knn3d_lookup(las, x, y, z, 10)
@@ -318,8 +318,8 @@ test_that("Spatial indexes work with 1 point", {
 
   expect_true(all(sapply(u, identical, 1L)))
 
-  u = vector("list", 4)
-  for(index in 0:4)
+  u = vector("list", 3)
+  for(index in 0:3)
   {
     index(las) <- index
     id = lidR:::C_circle_lookup(las, x, y, 5)
@@ -338,8 +338,8 @@ test_that("Spatial indexes work no point to find", {
   y = -100
   z = 0
 
-  u = vector("list", 4)
-  for(index in 0:4)
+  u = vector("list", 3)
+  for(index in 0:3)
   {
     index(las) <- index
     id = lidR:::C_circle_lookup(las, x, y, 5)