FEAT-449: Resolution query functions

.gitignore

@@ -62,7 +62,3 @@ share/python-wheels/
 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

CHANGELOG.md

@@ -7,6 +7,19 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Added
+
+- **\[Python/C++\]** Add spatial function `Crosses` to `VoxelKey`
+- **\[Python/C++\]** Add resolution function `Resolution` and `GetResolutionAtDepth` to `VoxelKey`
+- **\[Python/C++\]** Add `Vector3` operators to `Vector3`
+- **\[Python/C++\]** Add resolution query functions `GetDepthAtResolution`, `GetNodesAtResolution` and `GetNodesWithinResolution` to `Reader`
+- **\[Python/C++\]** Add `BoundsTrimFileExample` and `ResolutionTrimFileExample` to example-writer files
+### Changed
+
+- **\[Python/C++\]** Change order of arguments in `VoxelKey` spatial functions `Intersects`, `Contains`, and `Within`
+- **\[Python/C++\]** Add optional `resolution` argument to `Reader` spatial query functions `GetNodesWithinBox`, `GetNodesIntersectBox`, and `GetPointsWithinBox`. `resolution` can be used to limit the resolution during spatial queries
+- **\[Python/C++\]** Update `span` of `autzen-classified.copc.laz` test file from 0 to 128
+
 ## [1.3.0] - 2021-10-04
 
 ### Added

cmake/DownloadExampleFiles.cmake

@@ -2,11 +2,12 @@ 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}")
 			file (DOWNLOAD
-				"https://github.com/PDAL/data/raw/62e514b6484ec59cd48bb48d5c6fe8a00216a6ac/autzen/autzen-classified.copc.laz"
+				"https://cloudfront.rockrobotic.com/public/laz-examples/copc/drafts/span/autzen-classified.copc.laz"
 				${OUT_PATH}
 			)
 

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

@@ -29,6 +29,11 @@ struct Vector3
     Vector3 operator+(const double &d) const { return Vector3(x + d, y + d, z + d); }
     Vector3 operator-(const double &d) const { return Vector3(x - d, y - d, z - d); }
 
+    Vector3 operator*(const Vector3 &other) const { return Vector3(x * other.x, y * other.y, z * other.z); }
+    Vector3 operator/(const Vector3 &other) const { return Vector3(x / other.x, y / other.y, z / other.z); }
+    Vector3 operator+(const Vector3 &other) const { return Vector3(x + other.x, y + other.y, z + other.z); }
+    Vector3 operator-(const Vector3 &other) const { return Vector3(x - other.x, y - other.y, z - other.z); }
+
     std::string ToString()
     {
         std::stringstream ss;

cpp/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,16 @@ 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;
+    static double GetResolutionAtDepth(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

@@ -2,6 +2,7 @@
 #define COPCLIB_IO_READER_H_
 
 #include <istream>
+#include <limits>
 #include <string>
 
 #include "copc-lib/copc/file.hpp"
@@ -44,10 +45,19 @@ class Reader : public BaseIO
     // Helper function to get all points from the root
     las::Points GetAllPoints();
 
+    // Resolution query functions
+    // The resulting resolution may not be exactly this value: the minimum possible resolution that is at least as
+    // precise as the requested resolution will be selected. Therefore the result may be a bit more precise than
+    // requested.
+    int32_t GetDepthAtResolution(double resolution);
+    std::vector<Node> GetNodesAtResolution(double resolution);
+    std::vector<Node> GetNodesWithinResolution(double resolution);
+
     // 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 resolution = 0);
+    std::vector<Node> GetNodesIntersectBox(const Box &box, double resolution = 0);
+    las::Points GetPointsWithinBox(const Box &box, double resolution = 0);
 
   protected:
     Reader() = default;

cpp/src/hierarchy/key.cpp

@@ -69,15 +69,31 @@ bool VoxelKey::ChildOf(VoxelKey parent_key) const
     return false;
 }
 
-bool VoxelKey::Intersects(const Box &box, const las::LasHeader &header) const
+double VoxelKey::Resolution(const las::LasHeader &header, const las::CopcVlr &copc_info) const
+{
+    if (copc_info.span <= 0)
+        throw std::runtime_error("VoxelKey::Resolution: Octree span must be greater than 0.");
+    return (header.max.x - header.min.x) / copc_info.span / std::pow(2, d);
+}
+
+double VoxelKey::GetResolutionAtDepth(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);
 }
-bool VoxelKey::Contains(const Box &box, const las::LasHeader &header) const { return Box(*this, header).Contains(box); }
-bool VoxelKey::Contains(const Vector3 &point, const las::LasHeader &header) const
+bool VoxelKey::Contains(const las::LasHeader &header, const Box &box) const { return Box(*this, header).Contains(box); }
+bool VoxelKey::Contains(const las::LasHeader &header, const Vector3 &point) const
 {
     return Box(*this, header).Contains(point);
 }
-bool VoxelKey::Within(const Box &box, const las::LasHeader &header) const { return Box(*this, header).Within(box); }
+bool VoxelKey::Within(const las::LasHeader &header, const Box &box) const { return Box(*this, header).Within(box); }
+bool VoxelKey::Crosses(const las::LasHeader &header, const Box &box) const
+{
+    return Box(*this, header).Intersects(box) && !Box(*this, header).Within(box);
+}
 
 } // namespace copc

cpp/src/io/reader.cpp

@@ -1,5 +1,8 @@
-#include "copc-lib/io/reader.hpp"
+#include <cmath>
+#include <stdexcept>
+
 #include "copc-lib/hierarchy/internal/hierarchy.hpp"
+#include "copc-lib/io/reader.hpp"
 #include "copc-lib/las/header.hpp"
 #include "copc-lib/laz/decompressor.hpp"
 
@@ -186,53 +189,128 @@ las::Points Reader::GetAllPoints()
     return out;
 }
 
-std::vector<Node> Reader::GetNodesWithinBox(const Box &box)
+std::vector<Node> Reader::GetNodesWithinBox(const Box &box, double resolution)
 {
     std::vector<Node> out;
 
     auto header = GetLasHeader();
+    auto copc_info = GetCopcHeader();
+    auto max_depth = GetDepthAtResolution(resolution);
+
     for (const auto &node : GetAllChildren())
     {
-        if (node.key.Within(box, header))
+        if (node.key.Within(header, box) && node.key.d <= max_depth)
             out.push_back(node);
     }
 
     return out;
 }
 
-std::vector<Node> Reader::GetNodesIntersectBox(const Box &box)
+std::vector<Node> Reader::GetNodesIntersectBox(const Box &box, double resolution)
 {
     std::vector<Node> out;
 
     auto header = GetLasHeader();
+    auto copc_info = GetCopcHeader();
+    auto max_depth = GetDepthAtResolution(resolution);
+
+    // Get all nodes in octree
     for (const auto &node : GetAllChildren())
     {
-        if (node.key.Intersects(box, header))
+        if (node.key.Intersects(header, box) && node.key.d <= max_depth)
             out.push_back(node);
     }
 
     return out;
 }
 
-las::Points Reader::GetPointsWithinBox(const Box &box)
+las::Points Reader::GetPointsWithinBox(const Box &box, double resolution)
 {
     auto header = GetLasHeader();
+    auto copc_info = GetCopcHeader();
+    auto max_depth = GetDepthAtResolution(resolution);
     auto out = las::Points(header);
 
     // Get all nodes in octree
     for (const auto &node : GetAllChildren())
     {
-        // If node fits in Box
-        if (node.key.Intersects(box, header))
+        if (node.key.d <= max_depth)
         {
-            // Add points that fit in the box
-            auto points = GetPoints(node);
-            out.AddPoints(points.GetWithin(box));
+            // If node fits in Box
+            if (node.key.Within(header, box))
+            {
+                // If the node is within the box add all points
+                out.AddPoints(GetPoints(node));
+            }
+            else if (node.key.Intersects(header, box))
+            {
+                // If the node only crosses the box then get subset of points within box
+                auto points = GetPoints(node);
+                out.AddPoints(points.GetWithin(box));
+            }
         }
     }
     return out;
 }
 
+int32_t Reader::GetDepthAtResolution(double resolution)
+{
+    // Compute max depth
+    int32_t max_depth = -1;
+    // Get all nodes in octree
+    for (const auto &node : GetAllChildren())
+    {
+        if (node.key.d > max_depth)
+            max_depth = node.key.d;
+    }
+
+    // If query resolution is <=0 return the octree's max depth
+    if (resolution <= 0.0)
+        return max_depth;
+    if (GetCopcHeader().span <= 0)
+        throw std::runtime_error("Reader::GetDepthAtResolution: Octree span must be greater than 0.");
+
+    auto current_resolution = (GetLasHeader().max.x - GetLasHeader().min.x) / GetCopcHeader().span;
+
+    for (int32_t i = 0; i <= max_depth; i++)
+    {
+        if (current_resolution <= resolution)
+            return i;
+        current_resolution /= 2;
+    }
+    return max_depth;
+}
+
+std::vector<Node> Reader::GetNodesAtResolution(double resolution)
+{
+    auto target_depth = GetDepthAtResolution(resolution);
+
+    std::vector<Node> out;
+
+    for (const auto &node : GetAllChildren())
+    {
+        if (node.key.d == target_depth)
+            out.push_back(node);
+    }
+
+    return out;
+}
+
+std::vector<Node> Reader::GetNodesWithinResolution(double resolution)
+{
+    auto target_depth = GetDepthAtResolution(resolution);
+
+    std::vector<Node> out;
+
+    for (const auto &node : GetAllChildren())
+    {
+        if (node.key.d <= target_depth)
+            out.push_back(node);
+    }
+
+    return out;
+}
+
 las::EbVlr Reader::ReadExtraByteVlr(std::map<uint64_t, las::VlrHeader> &vlrs)
 {
     for (auto &[offset, vlr_header] : vlrs)