FEAT-449: Resolution query functions

example/example-writer.cpp

@@ -3,6 +3,7 @@
 #include <random>
 
 #include <copc-lib/geometry/vector3.hpp>
+#include <copc-lib/hierarchy/key.hpp>
 #include <copc-lib/io/reader.hpp>
 #include <copc-lib/io/writer.hpp>
 #include <copc-lib/las/header.hpp>
@@ -78,6 +79,102 @@ void TrimFileExample(bool compressor_example_flag)
     }
 }
 
+// In this example, we'll filter the points in the autzen dataset based on bounds.
+void BoundsTrimFileExample()
+{
+    // We'll get our point data from this file
+    FileReader reader("test/data/autzen-classified.copc.laz");
+    auto old_header = reader.GetLasHeader();
+
+    // Take horizontal 2D box of [200,200] roughly in the middle of the point cloud.
+    Box box(637190, 851109, 637390, 851309);
+    {
+        // Copy the header to the new file
+        Writer::LasConfig cfg(old_header, reader.GetExtraByteVlr());
+
+        // Now, we can create our actual writer, with an optional `span` and `wkt`:
+        FileWriter writer("test/data/autzen-bounds-trimmed.copc.laz", cfg, reader.GetCopcHeader().span,
+                          reader.GetWkt());
+
+        // The root page is automatically created and added for us
+        Page root_page = writer.GetRootPage();
+
+        for (const auto &node : reader.GetAllChildren())
+        {
+
+            if (node.key.Within(old_header, box))
+            {
+                // If node is within the box then add all points (without decompressing)
+                writer.AddNodeCompressed(root_page, node.key, reader.GetPointDataCompressed(node), node.point_count);
+            }
+            else if (node.key.Intersects(old_header, box))
+            {
+                // If node only crosses the box then decompress points data and get subset of points that are within the
+                // box
+                auto points = reader.GetPoints(node).GetWithin(box);
+                writer.AddNode(root_page, node.key, las::Points(points).Pack());
+            }
+        }
+
+        // Make sure we call close to finish writing the file!
+        writer.Close();
+    }
+
+    // Now, let's test our new file
+    FileReader new_reader("test/data/autzen-bounds-trimmed.copc.laz");
+
+    // Let's go through each point and make sure they fit within the Box
+    for (const auto &node : new_reader.GetAllChildren())
+    {
+        auto points = new_reader.GetPoints(node);
+        assert(points.Within(box));
+    }
+}
+
+// In this example, we'll filter the points in the autzen dataset based on resolution.
+void ResolutionTrimFileExample()
+{
+    // We'll get our point data from this file
+    FileReader reader("test/data/autzen-classified-new.copc.laz");
+    auto old_header = reader.GetLasHeader();
+
+    double resolution = 10;
+    {
+        // Copy the header to the new file
+        Writer::LasConfig cfg(old_header, reader.GetExtraByteVlr());
+
+        // Now, we can create our actual writer, with an optional `span` and `wkt`:
+        FileWriter writer("test/data/autzen-resolution-trimmed.copc.laz", cfg, reader.GetCopcHeader().span,
+                          reader.GetWkt());
+
+        // The root page is automatically created and added for us
+        Page root_page = writer.GetRootPage();
+
+        for (const auto &node : reader.GetAllChildren())
+        {
+            if (node.key.Resolution(writer.GetLasHeader(), writer.GetCopcHeader()) >= resolution)
+            {
+                writer.AddNodeCompressed(root_page, node.key, reader.GetPointDataCompressed(node), node.point_count);
+            }
+        }
+
+        // Make sure we call close to finish writing the file!
+        writer.Close();
+    }
+
+    // Now, let's test our new file
+    FileReader new_reader("test/data/autzen-resolution-trimmed.copc.laz");
+
+    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);
+    }
+}
+
 // constants
 const Vector3 MIN_BOUNDS = {-2000, -5000, 20};
 const Vector3 MAX_BOUNDS = {5000, 1034, 125};
@@ -174,5 +271,7 @@ int main()
 {
     TrimFileExample(false);
     TrimFileExample(true);
+    BoundsTrimFileExample();
+    ResolutionTrimFileExample();
     NewFileExample();
 }

example/example-writer.py

@@ -75,6 +75,102 @@ def TrimFileExample(compressor_example_flag):
     new_reader.Close()
 
 
+# In this example, we'll filter the points in the autzen dataset based on bounds.
+def BoundsTrimFileExample():
+    # We'll get our point data from this file
+    reader = copc.FileReader("../test/data/autzen-classified.copc.laz")
+    old_header = reader.GetLasHeader()
+
+    # Take horizontal 2D box of [200,200] roughly in the middle of the point cloud.
+    box = copc.Box(637190, 851109, 637390, 851309)
+
+    # Copy the header to the new file
+    cfg = copc.LasConfig(old_header, reader.GetExtraByteVlr())
+
+    # Now, we can create our actual writer, with an optional `span` and `wkt`:
+    writer = copc.FileWriter(
+        "../test/data/autzen-bounds-trimmed.copc.laz",
+        cfg,
+        reader.GetCopcHeader().span,
+        reader.GetWkt(),
+    )
+
+    # The root page is automatically created and added for us
+    root_page = writer.GetRootPage()
+
+    for node in reader.GetAllChildren():
+        if node.key.Within(old_header, box):
+            # If node is within the box then add all points (without decompressing)
+            writer.AddNodeCompressed(
+                root_page,
+                node.key,
+                reader.GetPointDataCompressed(node),
+                node.point_count,
+            )
+        elif node.key.Intersects(old_header, box):
+            # If node only crosses the box then decompress points data and get subset of points that are within the box
+            points = reader.GetPoints(node).GetWithin(box)
+            writer.AddNode(root_page, node.key, copc.Points(points).Pack())
+
+    # Make sure we call close to finish writing the file!
+    writer.Close()
+
+    # Now, let's test our new file
+    new_reader = copc.FileReader("../test/data/autzen-bounds-trimmed.copc.laz")
+
+    # Let's go through each point and make sure they fit in the within the Box
+    for node in new_reader.GetAllChildren():
+        points = new_reader.GetPoints(node)
+        assert points.Within(box)
+
+
+# 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("../test/data/autzen-classified-new.copc.laz")
+    old_header = reader.GetLasHeader()
+
+    resolution = 10
+
+    # Copy the header to the new file
+    cfg = copc.LasConfig(old_header, reader.GetExtraByteVlr())
+
+    # Now, we can create our actual writer, with an optional `span` and `wkt`:
+    writer = copc.FileWriter(
+        "../test/data/autzen-resolution-trimmed.copc.laz",
+        cfg,
+        reader.GetCopcHeader().span,
+        reader.GetWkt(),
+    )
+
+    # The root page is automatically created and added for us
+    root_page = writer.GetRootPage()
+
+    for node in reader.GetAllChildren():
+        if (
+            node.key.Resolution(writer.GetLasHeader(), writer.GetCopcHeader())
+            >= resolution
+        ):
+            writer.AddNodeCompressed(
+                root_page,
+                node.key,
+                reader.GetPointDataCompressed(node),
+                node.point_count,
+            )
+
+    # Make sure we call close to finish writing the file!
+    writer.Close()
+
+    # Now, let's test our new file
+    new_reader = copc.FileReader("../test/data/autzen-resolution-trimmed.copc.laz")
+
+    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
+
+
 # constants
 MIN_BOUNDS = copc.Vector3(-2000, -5000, 20)  # Argument Constructor
 MAX_BOUNDS = copc.Vector3([5000, 1034, 125])  # List Constructor
@@ -185,4 +281,6 @@ def NewFileExample():
 if __name__ == "__main__":
     TrimFileExample(False)
     TrimFileExample(True)
+    BoundsTrimFileExample()
+    ResolutionTrimFileExample()
     NewFileExample()

include/copc-lib/hierarchy/key.hpp

@@ -10,6 +10,7 @@
 
 #include <copc-lib/geometry/box.hpp>
 #include <copc-lib/geometry/vector3.hpp>
+#include <copc-lib/las/vlr.hpp>
 
 namespace copc
 {
@@ -46,10 +47,15 @@ class VoxelKey
     // Tests whether the current key is a child of a given key
     bool ChildOf(VoxelKey parent_key) const;
 
-    bool Intersects(const Box &box, const las::LasHeader &header) const;
-    bool Contains(const Box &vec, const las::LasHeader &header) const;
-    bool Contains(const Vector3 &point, const las::LasHeader &header) const;
-    bool Within(const Box &box, const las::LasHeader &header) const;
+    // Spatial query functions
+    // Definitions taken from https://shapely.readthedocs.io/en/stable/manual.html#binary-predicates
+    bool Intersects(const las::LasHeader &header, const Box &box) const;
+    bool Contains(const las::LasHeader &header, const Box &vec) const;
+    bool Contains(const las::LasHeader &header, const Vector3 &point) const;
+    bool Within(const las::LasHeader &header, const Box &box) const;
+    bool Crosses(const las::LasHeader &header, const Box &box) const;
+
+    double Resolution(const las::LasHeader &header, const las::CopcVlr &copc_info) const;
 
     int32_t d;
     int32_t x;

include/copc-lib/io/reader.hpp

@@ -2,6 +2,7 @@
 #define COPCLIB_IO_READER_H_
 
 #include <istream>
+#include <limits>
 #include <string>
 
 #include <copc-lib/copc/file.hpp>
@@ -41,9 +42,13 @@ class Reader : public BaseIO
     las::Points GetAllPoints();
 
     // Spatial query functions
-    std::vector<Node> GetNodesWithinBox(const Box &box);
-    std::vector<Node> GetNodesIntersectBox(const Box &box);
-    las::Points GetPointsWithinBox(const Box &box);
+    // 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 GetDepthWithResolution(double resolution) const;
+    std::vector<Node> GetNodesWithResolution(double resolution);
+    std::vector<Node> GetNodesDownToResolution(double resolution);
 
   protected:
     Reader() = default;

python/bindings.cpp

@@ -1,5 +1,6 @@
 #include <cmath>
 #include <iostream>
+#include <limits>
 #include <string>
 #include <vector>
 
@@ -56,10 +57,12 @@ PYBIND11_MODULE(copclib, m)
         .def("GetParent", &VoxelKey::GetParent)
         .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("Intersects", &VoxelKey::Intersects)
-        .def("Contains", py::overload_cast<const Box &, const las::LasHeader &>(&VoxelKey::Contains, py::const_))
-        .def("Contains", py::overload_cast<const Vector3 &, const las::LasHeader &>(&VoxelKey::Contains, py::const_))
+        .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_))
         .def("Within", &VoxelKey::Within)
+        .def("Crosses", &VoxelKey::Crosses)
         .def("__str__", &VoxelKey::ToString)
         .def("__repr__", &VoxelKey::ToString);
     py::implicitly_convertible<py::tuple, VoxelKey>();
@@ -129,9 +132,8 @@ PYBIND11_MODULE(copclib, m)
             "__truediv__", [](const Vector3 &vec, const double &d) { return vec / d; }, py::is_operator())
         .def(
             "__floordiv__",
-            [](const Vector3 &vec, const double &d) {
-                return Vector3(std::floor(vec.x / d), std::floor(vec.y / d), std::floor(vec.z / d));
-            },
+            [](const Vector3 &vec, const double &d)
+            { return Vector3(std::floor(vec.x / d), std::floor(vec.y / d), std::floor(vec.z / d)); },
             py::is_operator())
         .def(
             "__add__", [](const Vector3 &vec, const double &d) { return vec + d; }, py::is_operator())
@@ -240,7 +242,8 @@ PYBIND11_MODULE(copclib, m)
     using DiffType = ssize_t;
     using SizeType = size_t;
 
-    auto wrap_i = [](DiffType i, SizeType n) {
+    auto wrap_i = [](DiffType i, SizeType n)
+    {
         if (i < 0)
             i += n;
         if (i < 0 || (SizeType)i >= n)
@@ -268,6 +271,7 @@ PYBIND11_MODULE(copclib, m)
         .def("CreatePoint", &las::Points::CreatePoint)
         .def("ToPointFormat", &las::Points::ToPointFormat, py::arg("point_format_id"))
         .def("Within", &las::Points::Within, py::arg("box"))
+        .def("GetWithin", &las::Points::GetWithin, py::arg("box"))
         .def("Pack", py::overload_cast<>(&las::Points::Pack))
         .def("Unpack", py::overload_cast<const std::vector<char> &, const las::LasHeader &>(&las::Points::Unpack))
         .def("Unpack", py::overload_cast<const std::vector<char> &, const int8_t &, const uint16_t &, const Vector3 &,
@@ -353,9 +357,15 @@ 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)
-        .def("GetNodesIntersectBox", &Reader::GetNodesIntersectBox)
-        .def("GetPointsWithinBox", &Reader::GetPointsWithinBox);
+        .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("GetDepthWithResolution", &Reader::GetDepthWithResolution, py::arg("resolution"))
+        .def("GetNodesWithResolution", &Reader::GetNodesWithResolution, py::arg("resolution"))
+        .def("GetNodesDownToResolution", &Reader::GetNodesDownToResolution, py::arg("resolution"));
 
     py::class_<FileWriter>(m, "FileWriter")
         .def(py::init<const std::string &, Writer::LasConfig const &, const int &, const std::string &>(),