Add tile helper functions

CMakeLists.txt

@@ -88,7 +88,7 @@ endif()
 set(LAZPERF_VERSION "3.0.0")
 # First, see if the repo is cloned locally, and if so, compile it
 # this is required for CI build of sdist
-if(EXISTS "${CMAKE_SOURCE_DIR}/libs/laz-perf")
+if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/libs/laz-perf")
     set(WITH_TEST_TEMP ${WITH_TESTS})
     set(WITH_TESTS OFF) # never build lazperf tests
     add_subdirectory(libs/laz-perf)
@@ -106,7 +106,7 @@ else()
         set(EXTRA_EXPORT_TARGETS "")
     else ()
         # assume lazperf is compiled in tandem
-        message("Lazperf package not found, including ${LAZPERF_DIR}/cpp")
+        message("Lazperf package not found, assuming compiled in tandem...")
         set(LAZPERF_LIB_NAME "lazperf")
         set(EXTRA_EXPORT_TARGETS "lazperf")
     endif ()

cpp/CMakeLists.txt

@@ -7,6 +7,7 @@ set(${LIBRARY_TARGET_NAME}_HDR
         include/${LIBRARY_TARGET_NAME}/copc/copc_config.hpp
         include/${LIBRARY_TARGET_NAME}/geometry/box.hpp
         include/${LIBRARY_TARGET_NAME}/geometry/vector3.hpp
+        include/${LIBRARY_TARGET_NAME}/geometry/helpers.hpp
         include/${LIBRARY_TARGET_NAME}/hierarchy/entry.hpp
         include/${LIBRARY_TARGET_NAME}/hierarchy/key.hpp
         include/${LIBRARY_TARGET_NAME}/hierarchy/node.hpp
@@ -37,6 +38,7 @@ set(${LIBRARY_TARGET_NAME}_SRC
         src/copc/extents.cpp
         src/copc/copc_config.cpp
         src/geometry/box.cpp
+        src/geometry/helpers.cpp
         src/hierarchy/key.cpp
         src/hierarchy/page.cpp
         src/io/base_reader.cpp

cpp/include/copc-lib/geometry/helpers.hpp

@@ -0,0 +1,33 @@
+#ifndef CPP_INCLUDE_COPC_LIB_GEOMETRY_HELPERS
+#define CPP_INCLUDE_COPC_LIB_GEOMETRY_HELPERS
+
+#include <cmath>
+#include <cstdint>
+#include <vector>
+
+namespace copc
+{
+
+namespace las
+{
+class LasHeader;
+}
+class Box;
+
+enum RoundStrategy
+{
+    NEAREST,
+    PREFER_LARGER_TILE,
+    PREFER_SMALLER_TILE
+};
+
+// Returns the nearest depth to the requested tile size
+const double GetNearestDepth(double tile_size, const las::LasHeader &header,
+                             RoundStrategy rounding = RoundStrategy::NEAREST);
+
+const std::vector<copc::Box> GetPossibleTilesAtDepth(int32_t target_depth, const las::LasHeader &header);
+const std::vector<copc::Box> GetPossibleTilesWithSize(double target_tile_size, const las::LasHeader &header);
+
+} // namespace copc
+
+#endif /* CPP_INCLUDE_COPC_LIB_GEOMETRY_HELPERS */

cpp/include/copc-lib/hierarchy/key.hpp

@@ -76,6 +76,8 @@ class VoxelKey
 
     double Resolution(const las::LasHeader &header, const CopcInfo &copc_info) const;
     static double GetResolutionAtDepth(int32_t d, const las::LasHeader &header, const CopcInfo &copc_info);
+    double Span(const las::LasHeader &header) const;
+    static double GetSpanAtDepth(int32_t d, const las::LasHeader &header);
 
     int32_t d;
     int32_t x;

cpp/src/geometry/helpers.cpp

@@ -0,0 +1,98 @@
+#include "copc-lib/geometry/helpers.hpp"
+#include "copc-lib/hierarchy/key.hpp"
+#include "copc-lib/las/header.hpp"
+
+#include <limits>
+
+namespace copc
+{
+
+const double GetNearestDepth(double tile_size, const las::LasHeader &header, RoundStrategy rounding)
+{
+    double depth = std::log2(header.Span() / tile_size);
+    int nearest_depth;
+    if (rounding == RoundStrategy::NEAREST)
+        nearest_depth = std::lround(depth);
+    else if (rounding == RoundStrategy::PREFER_LARGER_TILE)
+        nearest_depth = (int)std::floor(depth);
+    else if (rounding == RoundStrategy::PREFER_SMALLER_TILE)
+        nearest_depth = (int)std::ceil(depth);
+    return std::max(0, nearest_depth);
+}
+
+// Generates a vector of Boxes that tile the dataset in its entirety
+// Note that points in the same node may fall within different tiles,
+// due to rounding
+// Additionally, some boxes may contain no points at all, since
+// we cannot know for certain whether a tile has points without actually loading it
+const std::vector<copc::Box> GetPossibleTilesAtDepth(int32_t target_depth, const las::LasHeader &header)
+{
+    std::vector<copc::Box> possible_tiles;
+
+    // We can compute the maximum VoxelKey coordinate based on the depth
+    int max_possible_coord = std::pow(2, target_depth);
+    // But then, we can limit the search space in each dimension,
+    // since not every dimension will use the full span
+    // (since the span is based on the axis with the largest range)
+    double tile_size = header.Span() / max_possible_coord;
+    int max_x_coord = (int)std::ceil((header.max.x - header.min.x) / tile_size);
+    int max_y_coord = (int)std::ceil((header.max.y - header.min.y) / tile_size);
+    int max_z_coord = (int)std::ceil((header.max.z - header.min.z) / tile_size);
+
+    for (int x = 0; x < max_x_coord; x++)
+        for (int y = 0; y < max_y_coord; y++)
+            for (int z = 0; z < max_z_coord; z++)
+            {
+                auto box = copc::Box(copc::VoxelKey(target_depth, x, y, z), header);
+                // To avoid the same point falling into two boxes when the point's coordinates
+                // are the exact same as the box's limit,
+                // we offset the minimum coordinate in each dimension by the smallest possible
+                // floating point value so that two neighboring boxes don't have the exact same max and min
+                // but, don't do that to the first coordinate in each dimension since there's no other box
+                // for that point to fall in
+                if (x != 0)
+                    box.x_min = std::nextafter(box.x_min, std::numeric_limits<double>::max());
+                if (y != 0)
+                    box.y_min = std::nextafter(box.y_min, std::numeric_limits<double>::max());
+                if (z != 0)
+                    box.z_min = std::nextafter(box.z_min, std::numeric_limits<double>::max());
+                possible_tiles.push_back(box);
+            }
+    return possible_tiles;
+}
+
+const std::vector<copc::Box> GetPossibleTilesWithSize(double target_tile_size, const las::LasHeader &header)
+{
+    std::vector<copc::Box> possible_tiles;
+
+    int max_num_tiles = header.Span() / target_tile_size;
+
+    for (int x = 0; x < max_num_tiles; x++)
+        for (int y = 0; y < max_num_tiles; y++)
+            for (int z = 0; z < max_num_tiles; z++)
+            {
+                double x_min = target_tile_size * x + header.min.x;
+                double y_min = target_tile_size * y + header.min.y;
+                double z_min = target_tile_size * z + header.min.z;
+                double x_max = x_min + target_tile_size;
+                double y_max = y_min + target_tile_size;
+                double z_max = z_min + target_tile_size;
+                auto box = copc::Box(x_min, y_min, z_min, x_max, y_max, z_max);
+                // To avoid the same point falling into two boxes when the point's coordinates
+                // are the exact same as the box's limit,
+                // we offset the minimum coordinate in each dimension by the smallest possible
+                // floating point value so that two neighboring boxes don't have the exact same max and min
+                // but, don't do that to the first coordinate in each dimension since there's no other box
+                // for that point to fall in
+                if (x != 0)
+                    box.x_min = std::nextafter(box.x_min, std::numeric_limits<double>::max());
+                if (y != 0)
+                    box.y_min = std::nextafter(box.y_min, std::numeric_limits<double>::max());
+                if (z != 0)
+                    box.z_min = std::nextafter(box.z_min, std::numeric_limits<double>::max());
+                possible_tiles.push_back(box);
+            }
+    return possible_tiles;
+}
+
+} // namespace copc

cpp/src/hierarchy/key.cpp

@@ -98,6 +98,10 @@ double VoxelKey::GetResolutionAtDepth(int32_t d, const las::LasHeader &header, c
     return VoxelKey(d, 0, 0, 0).Resolution(header, copc_info);
 }
 
+double VoxelKey::Span(const las::LasHeader &header) const { return header.Span() * std::pow(2, -d); }
+
+double VoxelKey::GetSpanAtDepth(int32_t d, const las::LasHeader &header) { return VoxelKey(d, 0, 0, 0).Span(header); }
+
 bool VoxelKey::Intersects(const las::LasHeader &header, const Box &box) const
 {
     return Box(*this, header).Intersects(box);