FEAT-449: Resolution query functions

.gitignore

@@ -61,8 +61,4 @@ share/python-wheels/
 *.egg
 MANIFEST
 
-/CMakeSettings.json
-/test/data/autzen-classified.copc.laz
-/test/data/writer_test.copc.laz
-/test/data/autzen-trimmed.copc.laz
-/test/data/new-copc.copc.laz
+/CMakeSettings.json

CMakeLists.txt

@@ -122,6 +122,7 @@ include(AddUninstallTarget)
 
 # Build Submodules
 if (WITH_TESTS_AND_EXAMPLES)
+    set(BaseAutzenFilePath "${CMAKE_SOURCE_DIR}/test/data/autzen-classified.copc.laz") # TEMPORARY!!!!!!!!!!!!!! REMOVE ME
     add_subdirectory(test)
     add_subdirectory(example)
 endif()

cmake/DownloadExampleFiles.cmake

@@ -2,6 +2,7 @@ macro(download_test_files OUT_PATH)
 
 	if(NOT EXISTS ${OUT_PATH})
 		if (EXISTS ${BaseAutzenFilePath})
+			file(WRITE ${OUT_PATH}) # this will create any necessary pathes, as CREATE_LINK doens't
 			file(CREATE_LINK ${BaseAutzenFilePath} ${OUT_PATH})
 		else()
 			message(STATUS "Downloading test files to ${OUT_PATH}")

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

@@ -56,6 +56,7 @@ class VoxelKey
     bool Crosses(const las::LasHeader &header, const Box &box) const;
 
     double Resolution(const las::LasHeader &header, const las::CopcVlr &copc_info) const;
+    static double GetDepthResolution(int32_t d, const las::LasHeader &header, const las::CopcVlr &copc_info);
 
     int32_t d;
     int32_t x;

cpp/include/copc-lib/io/reader.hpp

@@ -47,10 +47,10 @@ class Reader : public BaseIO
 
     // Spatial query functions
     // Definitions taken from https://shapely.readthedocs.io/en/stable/manual.html#binary-predicates
-    std::vector<Node> GetNodesWithinBox(const Box &box, double min_resolution = std::numeric_limits<double>::min());
-    std::vector<Node> GetNodesIntersectBox(const Box &box, double min_resolution = std::numeric_limits<double>::min());
-    las::Points GetPointsWithinBox(const Box &box, double min_resolution = std::numeric_limits<double>::min());
-    int32_t GetDepthAtResolution(double resolution) const;
+    std::vector<Node> GetNodesWithinBox(const Box &box, double min_resolution = 0);
+    std::vector<Node> GetNodesIntersectBox(const Box &box, double min_resolution = 0);
+    las::Points GetPointsWithinBox(const Box &box, double min_resolution = 0);
+    int32_t GetDepthAtResolution(double resolution);
     std::vector<Node> GetNodesAtResolution(double resolution);
     std::vector<Node> GetNodesWithinResolution(double resolution);
 

cpp/src/hierarchy/key.cpp

@@ -76,6 +76,11 @@ double VoxelKey::Resolution(const las::LasHeader &header, const las::CopcVlr &co
     return (header.max.x - header.min.x) / copc_info.span / std::pow(2, d);
 }
 
+double VoxelKey::GetDepthResolution(int32_t d, const las::LasHeader &header, const las::CopcVlr &copc_info)
+{
+    return VoxelKey(d, 0, 0, 0).Resolution(header, copc_info);
+}
+
 bool VoxelKey::Intersects(const las::LasHeader &header, const Box &box) const
 {
     return Box(*this, header).Intersects(box);

cpp/src/io/reader.cpp

@@ -195,9 +195,11 @@ std::vector<Node> Reader::GetNodesWithinBox(const Box &box, double min_resolutio
 
     auto header = GetLasHeader();
     auto copc_info = GetCopcHeader();
+    auto max_depth = GetDepthAtResolution(min_resolution);
+
     for (const auto &node : GetAllChildren())
     {
-        if (node.key.Within(header, box) && node.key.Resolution(header, copc_info) >= min_resolution)
+        if (node.key.Within(header, box) && node.key.d <= max_depth)
             out.push_back(node);
     }
 
@@ -210,9 +212,12 @@ std::vector<Node> Reader::GetNodesIntersectBox(const Box &box, double min_resolu
 
     auto header = GetLasHeader();
     auto copc_info = GetCopcHeader();
+    auto max_depth = GetDepthAtResolution(min_resolution);
+
+    // Get all nodes in octree
     for (const auto &node : GetAllChildren())
     {
-        if (node.key.Intersects(header, box) && node.key.Resolution(header, copc_info) >= min_resolution)
+        if (node.key.Intersects(header, box) && node.key.d <= max_depth)
             out.push_back(node);
     }
 
@@ -223,12 +228,13 @@ las::Points Reader::GetPointsWithinBox(const Box &box, double min_resolution)
 {
     auto header = GetLasHeader();
     auto copc_info = GetCopcHeader();
+    auto max_depth = GetDepthAtResolution(min_resolution);
     auto out = las::Points(header);
 
     // Get all nodes in octree
     for (const auto &node : GetAllChildren())
     {
-        if (node.key.Resolution(header, copc_info) >= min_resolution)
+        if (node.key.d <= max_depth)
         {
             // If node fits in Box
             if (node.key.Within(header, box))
@@ -247,21 +253,29 @@ las::Points Reader::GetPointsWithinBox(const Box &box, double min_resolution)
     return out;
 }
 
-int32_t Reader::GetDepthAtResolution(double resolution) const
+int32_t Reader::GetDepthAtResolution(double resolution)
 {
+    // If query resolution is <=0 return the octree's max depth
     if (resolution <= 0)
-        throw std::runtime_error("Reader::GetDepthAtResolution: Query resolution must be greater than 0.");
+    {
+        int32_t max_depth = 0;
+        // Get all nodes in octree
+        for (const auto &node : GetAllChildren())
+        {
+            if (node.key.d > max_depth)
+                max_depth++;
+        }
+        return max_depth;
+    }
     if (GetCopcHeader().span <= 0)
         throw std::runtime_error("Reader::GetDepthAtResolution: Octree span must be greater than 0.");
 
     auto root_resolution = (GetLasHeader().max.x - GetLasHeader().min.x) / GetCopcHeader().span;
-    return static_cast<int32_t>(std::floor(std::max(0.0, std::log(root_resolution / resolution) / std::log(2))));
+    return static_cast<int32_t>(std::ceil(std::max(0.0, std::log(root_resolution / resolution) / std::log(2))));
 }
 
 std::vector<Node> Reader::GetNodesAtResolution(double resolution)
 {
-    if (resolution <= 0)
-        throw std::runtime_error("Query resolution must be greater than 0.");
     auto target_depth = GetDepthAtResolution(resolution);
 
     std::vector<Node> out;
@@ -277,8 +291,6 @@ std::vector<Node> Reader::GetNodesAtResolution(double resolution)
 
 std::vector<Node> Reader::GetNodesWithinResolution(double resolution)
 {
-    if (resolution <= 0)
-        throw std::runtime_error("Query resolution must be greater than 0.");
     auto target_depth = GetDepthAtResolution(resolution);
 
     std::vector<Node> out;

example/example-writer.cpp

@@ -137,7 +137,7 @@ void BoundsTrimFileExample()
 void ResolutionTrimFileExample()
 {
     // We'll get our point data from this file
-    FileReader reader("autzen-classified-new.copc.laz");
+    FileReader reader("autzen-classified.copc.laz");
     auto old_header = reader.GetLasHeader();
 
     double resolution = 10;
@@ -151,9 +151,11 @@ void ResolutionTrimFileExample()
         // The root page is automatically created and added for us
         Page root_page = writer.GetRootPage();
 
+        auto max_depth = reader.GetDepthAtResolution(resolution);
+
         for (const auto &node : reader.GetAllChildren())
         {
-            if (node.key.Resolution(writer.GetLasHeader(), writer.GetCopcHeader()) >= resolution)
+            if (node.key.d <= max_depth)
             {
                 writer.AddNodeCompressed(root_page, node.key, reader.GetPointDataCompressed(node), node.point_count);
             }
@@ -169,11 +171,9 @@ void ResolutionTrimFileExample()
     auto new_header = new_reader.GetLasHeader();
     auto new_copc_info = new_reader.GetCopcHeader();
 
-    // Let's go through each node we've written and make sure the resolution is correct
-    for (const auto &node : new_reader.GetAllChildren())
-    {
-        assert(node.key.Resolution(new_header, new_copc_info) >= resolution);
-    }
+    // Let's make sure the max resolution is at least as much as we requested
+    auto max_depth = new_reader.GetDepthAtResolution(0);
+    assert(VoxelKey::GetDepthResolution(max_depth, new_header, new_copc_info) <= resolution);
 }
 
 // constants

example/example-writer.py

@@ -127,7 +127,7 @@ def BoundsTrimFileExample():
 # In this example, we'll filter the points in the autzen dataset based on resolution.
 def ResolutionTrimFileExample():
     # We'll get our point data from this file
-    reader = copc.FileReader("autzen-classified-new.copc.laz")
+    reader = copc.FileReader("autzen-classified.copc.laz")
     old_header = reader.GetLasHeader()
 
     resolution = 10
@@ -145,12 +145,10 @@ def ResolutionTrimFileExample():
 
     # The root page is automatically created and added for us
     root_page = writer.GetRootPage()
+    max_depth = reader.GetDepthAtResolution(resolution)
 
     for node in reader.GetAllChildren():
-        if (
-            node.key.Resolution(writer.GetLasHeader(), writer.GetCopcHeader())
-            >= resolution
-        ):
+        if node.key.d <= max_depth:
             writer.AddNodeCompressed(
                 root_page,
                 node.key,
@@ -166,9 +164,12 @@ def ResolutionTrimFileExample():
 
     new_las_header = new_reader.GetLasHeader()
     new_copc_header = new_reader.GetCopcHeader()
-    # Let's go through each node we've written and make sure the resolution is correct
-    for node in new_reader.GetAllChildren():
-        assert node.key.Resolution(new_las_header, new_copc_header) >= resolution
+    # Let's make sure the max resolution is at least as much as we requested
+    max_depth = new_reader.GetDepthAtResolution(0)
+    assert (
+        copc.VoxelKey.GetDepthResolution(max_depth, new_las_header, new_copc_header)
+        <= resolution
+    )
 
 
 # constants

python/bindings.cpp

@@ -58,6 +58,8 @@ PYBIND11_MODULE(copclib, m)
         .def("GetParents", &VoxelKey::GetParents, py::arg("include_current"))
         .def("ChildOf", &VoxelKey::ChildOf, py::arg("parent_key"))
         .def("Resolution", &VoxelKey::Resolution, py::arg("las_header"), py::arg("copc_header"))
+        .def_static("GetDepthResolution", &VoxelKey::GetDepthResolution, py::arg("depth"), py::arg("las_header"),
+                    py::arg("copc_header"))
         .def("Intersects", &VoxelKey::Intersects)
         .def("Contains", py::overload_cast<const las::LasHeader &, const Box &>(&VoxelKey::Contains, py::const_))
         .def("Contains", py::overload_cast<const las::LasHeader &, const Vector3 &>(&VoxelKey::Contains, py::const_))
@@ -372,12 +374,9 @@ PYBIND11_MODULE(copclib, m)
         .def("GetAllChildren", py::overload_cast<const VoxelKey &>(&Reader::GetAllChildren), py::arg("key"))
         .def("GetAllChildren", py::overload_cast<>(&Reader::GetAllChildren))
         .def("GetAllPoints", &Reader::GetAllPoints)
-        .def("GetNodesWithinBox", &Reader::GetNodesWithinBox, py::arg("box"),
-             py::arg("min_resolution") = std::numeric_limits<double>::min())
-        .def("GetNodesIntersectBox", &Reader::GetNodesIntersectBox, py::arg("box"),
-             py::arg("min_resolution") = std::numeric_limits<double>::min())
-        .def("GetPointsWithinBox", &Reader::GetPointsWithinBox, py::arg("box"),
-             py::arg("min_resolution") = std::numeric_limits<double>::min())
+        .def("GetNodesWithinBox", &Reader::GetNodesWithinBox, py::arg("box"), py::arg("min_resolution") = 0)
+        .def("GetNodesIntersectBox", &Reader::GetNodesIntersectBox, py::arg("box"), py::arg("min_resolution") = 0)
+        .def("GetPointsWithinBox", &Reader::GetPointsWithinBox, py::arg("box"), py::arg("min_resolution") = 0)
         .def("GetDepthAtResolution", &Reader::GetDepthAtResolution, py::arg("resolution"))
         .def("GetNodesAtResolution", &Reader::GetNodesAtResolution, py::arg("resolution"))
         .def("GetNodesWithinResolution", &Reader::GetNodesWithinResolution, py::arg("resolution"));

test/reader_test.cpp

@@ -17,7 +17,7 @@ TEST_CASE("Reader tests", "[Reader]")
         SECTION("GetCopc Test")
         {
             auto copc = reader.GetCopcHeader();
-            REQUIRE(copc.span == 0);
+            REQUIRE(copc.span == 128);
             REQUIRE(copc.root_hier_offset == 93169718);
             REQUIRE(copc.root_hier_size == 8896);
             REQUIRE(copc.laz_vlr_offset == 643);
@@ -56,7 +56,7 @@ TEST_CASE("Reader tests", "[Reader]")
         SECTION("GetCopc Test")
         {
             auto copc = reader.GetCopcHeader();
-            REQUIRE(copc.span == 0);
+            REQUIRE(copc.span == 128);
             REQUIRE(copc.root_hier_offset == 93169718);
             REQUIRE(copc.root_hier_size == 8896);
             REQUIRE(copc.laz_vlr_offset == 643);
@@ -191,7 +191,11 @@ TEST_CASE("Point Error Handling Test", "[Reader]")
 
 TEST_CASE("Spatial Query Functions", "[Reader]")
 {
-    FileReader reader("autzen-classified-new.copc.laz");
+    FileReader reader("autzen-classified.copc.laz");
+
+    // Make horizontal 2D box of [200,200] roughly in the middle of the point cloud.
+    auto middle = (reader.GetLasHeader().max + reader.GetLasHeader().min) / 2;
+    Box middle_box(middle.x - 200, middle.y - 200, middle.x + 200, middle.y + 200);
 
     SECTION("GetNodesWithinBox")
     {
@@ -207,8 +211,7 @@ TEST_CASE("Spatial Query Functions", "[Reader]")
 
     SECTION("GetNodesIntersectBox")
     {
-        // Take horizontal 2D box of [200,200] roughly in the middle of the point cloud.
-        auto subset_nodes = reader.GetNodesIntersectBox(Box(637190, 851109, 637390, 851309));
+        auto subset_nodes = reader.GetNodesIntersectBox(middle_box);
         REQUIRE(subset_nodes.size() == 13);
     }
 
@@ -233,28 +236,29 @@ TEST_CASE("Spatial Query Functions", "[Reader]")
         }
         {
             // Take horizontal 2D box of [200,200] roughly in the middle of the point cloud.
-            auto subset_points = reader.GetPointsWithinBox(Box(637190, 851109, 637390, 851309));
-            REQUIRE(subset_points.Get().size() == 22902);
+            auto subset_points = reader.GetPointsWithinBox(middle_box);
+            REQUIRE(subset_points.Get().size() == 91178);
         }
     }
 
     SECTION("GetDepthAtResolution")
     {
-        REQUIRE(reader.GetDepthAtResolution(1) == 3);
-        REQUIRE_THROWS(reader.GetDepthAtResolution(0));
+        REQUIRE(reader.GetDepthAtResolution(1) == 4);
+        REQUIRE(reader.GetDepthAtResolution(0) == 5);
     }
 
     SECTION("GetNodesAtResolution")
     {
         auto subset_nodes = reader.GetNodesAtResolution(1);
-        REQUIRE(subset_nodes.size() == 48);
-        REQUIRE_THROWS(reader.GetDepthAtResolution(0));
+        REQUIRE(subset_nodes.size() == 192);
+        for (const auto &node : reader.GetNodesAtResolution(0))
+            REQUIRE(node.key.d == 5);
     }
 
     SECTION("GetNodesWithinResolution")
     {
         auto subset_nodes = reader.GetNodesWithinResolution(1);
-        REQUIRE(subset_nodes.size() == 65);
-        REQUIRE_THROWS(reader.GetDepthAtResolution(0));
+        REQUIRE(subset_nodes.size() == 257);
+        REQUIRE(reader.GetNodesWithinResolution(0).size() == reader.GetAllChildren().size());
     }
 }