diff --git a/CMakeLists.txt b/CMakeLists.txt index f229b6a4..cc43e0e8 100644 --- a/CMakeLists.txt +++ b/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 () diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index c521cd18..6190d40e 100644 --- a/cpp/CMakeLists.txt +++ b/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 diff --git a/cpp/include/copc-lib/geometry/helpers.hpp b/cpp/include/copc-lib/geometry/helpers.hpp new file mode 100644 index 00000000..77bdefd6 --- /dev/null +++ b/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 +#include +#include + +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 int GetNearestDepth(double tile_size, const las::LasHeader &header, + RoundStrategy rounding = RoundStrategy::NEAREST); + +const std::vector GetPossibleTilesAtDepth(int32_t target_depth, const las::LasHeader &header); +const std::vector GetPossibleTilesWithSize(double target_tile_size, const las::LasHeader &header); + +} // namespace copc + +#endif /* CPP_INCLUDE_COPC_LIB_GEOMETRY_HELPERS */ diff --git a/cpp/include/copc-lib/hierarchy/key.hpp b/cpp/include/copc-lib/hierarchy/key.hpp index 3b6d6ef9..bb29c143 100644 --- a/cpp/include/copc-lib/hierarchy/key.hpp +++ b/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; diff --git a/cpp/include/copc-lib/las/point.hpp b/cpp/include/copc-lib/las/point.hpp index 9e3c53f2..9bc81cfe 100644 --- a/cpp/include/copc-lib/las/point.hpp +++ b/cpp/include/copc-lib/las/point.hpp @@ -22,13 +22,28 @@ class Point #pragma region XYZ // XYZ Scaled double X() const { return x_scaled_; } - void X(const double &x) { x_scaled_ = x; } + void X(const double &x) + { + if (!std::isfinite(x)) + throw std::runtime_error("Tried to set non-finite X value!"); + x_scaled_ = x; + } double Y() const { return y_scaled_; } - void Y(const double &y) { y_scaled_ = y; } + void Y(const double &y) + { + if (!std::isfinite(y)) + throw std::runtime_error("Tried to set non-finite Y value!"); + y_scaled_ = y; + } double Z() const { return z_scaled_; } - void Z(const double &z) { z_scaled_ = z; } + void Z(const double &z) + { + if (!std::isfinite(z)) + throw std::runtime_error("Tried to set non-finite Z value!"); + z_scaled_ = z; + } #pragma endregion XYZ diff --git a/cpp/include/copc-lib/las/utils.hpp b/cpp/include/copc-lib/las/utils.hpp index 75af069c..c6ad9742 100644 --- a/cpp/include/copc-lib/las/utils.hpp +++ b/cpp/include/copc-lib/las/utils.hpp @@ -17,8 +17,14 @@ bool FormatHasNir(const uint8_t &point_format_id); template double ApplyScale(T value, double scale, double offset) { return (value * scale) + offset; } template T RemoveScale(double value, double scale, double offset) { + if (!std::isfinite(value)) + throw std::runtime_error("The input value " + std::to_string(value) + " is not finite."); + double val = std::round((value - offset) / scale); + if (!std::isfinite(val)) + throw std::runtime_error("The output value " + std::to_string(val) + " is not finite."); + if (val < std::numeric_limits::min() || val > std::numeric_limits::max()) throw std::runtime_error("The value " + std::to_string(value) + " is too large to save into the requested format." + diff --git a/cpp/src/geometry/helpers.cpp b/cpp/src/geometry/helpers.cpp new file mode 100644 index 00000000..9ae05f98 --- /dev/null +++ b/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 + +namespace copc +{ + +const int 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 (int)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 GetPossibleTilesAtDepth(int32_t target_depth, const las::LasHeader &header) +{ + std::vector 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::max()); + if (y != 0) + box.y_min = std::nextafter(box.y_min, std::numeric_limits::max()); + if (z != 0) + box.z_min = std::nextafter(box.z_min, std::numeric_limits::max()); + possible_tiles.push_back(box); + } + return possible_tiles; +} + +const std::vector GetPossibleTilesWithSize(double target_tile_size, const las::LasHeader &header) +{ + std::vector 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::max()); + if (y != 0) + box.y_min = std::nextafter(box.y_min, std::numeric_limits::max()); + if (z != 0) + box.z_min = std::nextafter(box.z_min, std::numeric_limits::max()); + possible_tiles.push_back(box); + } + return possible_tiles; +} + +} // namespace copc diff --git a/cpp/src/hierarchy/key.cpp b/cpp/src/hierarchy/key.cpp index 9c34bec7..1a4124fd 100644 --- a/cpp/src/hierarchy/key.cpp +++ b/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);