From 1182c513e3003b59adae6c12e23293fd2631f55f Mon Sep 17 00:00:00 2001 From: reiter Date: Tue, 17 Feb 2026 17:10:30 +0100 Subject: [PATCH 1/3] refactor: point neighborhood generation uses hash grid --- include/viennaray/rayPointNeighborhood.hpp | 231 +++++++----------- .../pointNeighborhood2D.cpp | 5 + 2 files changed, 93 insertions(+), 143 deletions(-) diff --git a/include/viennaray/rayPointNeighborhood.hpp b/include/viennaray/rayPointNeighborhood.hpp index e60fed1..12e2e1e 100644 --- a/include/viennaray/rayPointNeighborhood.hpp +++ b/include/viennaray/rayPointNeighborhood.hpp @@ -3,6 +3,7 @@ #include #include +#include #include namespace viennaray { @@ -10,9 +11,29 @@ namespace viennaray { using namespace viennacore; template class PointNeighborhood { + static_assert(D == 2 || D == 3, "Only 2D and 3D are supported"); + std::vector> pointNeighborhood_; NumericType distance_ = 0.; + // Cell index type + using CellIndex = std::array; + + struct CellIndexHash { + size_t operator()(const CellIndex &c) const { + // FNV-1a inspired hash combining + size_t h = 2166136261u; + for (int i = 0; i < D; ++i) { + h ^= static_cast(c[i]); + h *= 16777619u; + } + return h; + } + }; + + using GridMap = + std::unordered_map, CellIndexHash>; + public: PointNeighborhood() = default; @@ -20,48 +41,33 @@ template class PointNeighborhood { void init(std::vector> const &points, NumericType distance, Vec3D const &minCoords, Vec3D const &maxCoords) { + static_assert(Dim >= static_cast(D), + "Point dimension must be >= D"); distance_ = distance; const auto numPoints = points.size(); pointNeighborhood_.clear(); - pointNeighborhood_.resize(numPoints, std::vector{}); - if constexpr (D == 3) { - std::vector side1; - std::vector side2; - - // create copy of bounding box - Vec3D min = minCoords; - Vec3D max = maxCoords; - - std::vector dirs; - for (int i = 0; i < 3; ++i) { - if (min[i] != max[i]) { - dirs.push_back(i); - } - } - dirs.shrink_to_fit(); + pointNeighborhood_.resize(numPoints); - int dirIdx = 0; - NumericType pivot = (max[dirs[dirIdx]] + min[dirs[dirIdx]]) / 2; + if (numPoints == 0 || distance_ <= 0) + return; - // divide point data - for (unsigned int idx = 0; idx < numPoints; ++idx) { - if (points[idx][dirs[dirIdx]] <= pivot) { - side1.push_back(idx); - } else { - side2.push_back(idx); - } - } - createNeighborhood(points, side1, side2, min, max, dirIdx, dirs, pivot); - } else { - /// TODO: 2D divide and conquer algorithm - for (unsigned int idx1 = 0; idx1 < numPoints; ++idx1) { - for (unsigned int idx2 = idx1 + 1; idx2 < numPoints; ++idx2) { - if (checkDistance(points[idx1], points[idx2])) { - pointNeighborhood_[idx1].push_back(idx2); - pointNeighborhood_[idx2].push_back(idx1); - } - } - } + const NumericType invCellSize = NumericType(1) / distance_; + + // Build the grid: map cell index -> list of point indices + GridMap grid; + grid.reserve(numPoints); + for (unsigned int idx = 0; idx < numPoints; ++idx) { + CellIndex cell = computeCell(points[idx], minCoords, invCellSize); + grid[cell].push_back(idx); + } + + // For each point, check all neighboring cells + for (unsigned int idx = 0; idx < numPoints; ++idx) { + const auto &point = points[idx]; + CellIndex cell = computeCell(point, minCoords, invCellSize); + + // Iterate over the (2D: 3x3 = 9, 3D: 3x3x3 = 27) neighborhood of cells + iterateNeighborCells(grid, points, idx, point, cell); } } @@ -78,124 +84,63 @@ template class PointNeighborhood { [[nodiscard]] NumericType getDistance() const { return distance_; } private: - void createNeighborhood(const std::vector> &points, - const std::vector &side1, - const std::vector &side2, - const Vec3D &min, - const Vec3D &max, const int &dirIdx, - const std::vector &dirs, - const NumericType &pivot) { - assert(0 <= dirIdx && dirIdx < dirs.size() && "Assumption"); - if (side1.size() + side2.size() <= 1) { - return; + template + CellIndex computeCell(const VectorType &point, + const Vec3D &minCoords, + NumericType invCellSize) const { + CellIndex cell; + for (int i = 0; i < D; ++i) { + cell[i] = + static_cast(std::floor((point[i] - minCoords[i]) * invCellSize)); } + return cell; + } - // Corner case - // The pivot element should actually be between min and max. - if (pivot == min[dirs[dirIdx]] || pivot == max[dirs[dirIdx]]) { - // In this case the points are extremely close to each other (with respect - // to the floating point precision). - assert((min[dirs[dirIdx]] + max[dirs[dirIdx]]) / 2 == pivot && - "Characterization of corner case"); - auto sides = std::vector(side1); - sides.insert(sides.end(), side2.begin(), side2.end()); - // Add each of them to the neighborhoods - for (unsigned int idx1 = 0; idx1 < sides.size() - 1; ++idx1) { - for (unsigned int idx2 = idx1 + 1; idx2 < sides.size(); ++idx2) { - auto const &pi1 = sides[idx1]; - auto const &pi2 = sides[idx2]; - assert(pi1 != pi2 && "Assumption"); - pointNeighborhood_[pi1].push_back(pi2); - pointNeighborhood_[pi2].push_back(pi1); + template + void iterateNeighborCells( + const GridMap &grid, + const std::vector> &points, unsigned int idx, + const VectorType &point, const CellIndex &cell) { + if constexpr (D == 2) { + for (int dx = -1; dx <= 1; ++dx) { + for (int dy = -1; dy <= 1; ++dy) { + CellIndex neighbor = {cell[0] + dx, cell[1] + dy}; + checkCellNeighbors(grid, points, idx, point, neighbor); } } - return; - } - - // sets of candidates - std::vector side1Cand; - std::vector side2Cand; - - int newDirIdx = (dirIdx + 1) % static_cast(dirs.size()); - NumericType newPivot = (max[dirs[newDirIdx]] + min[dirs[newDirIdx]]) / 2; - - // recursion sets - std::vector s1r1set; - std::vector s1r2set; - std::vector s2r1set; - std::vector s2r2set; - - for (unsigned int idx : side1) { - const auto &point = points[idx]; - assert(point[dirs[dirIdx]] <= pivot && "Correctness Assertion"); - if (point[dirs[newDirIdx]] <= newPivot) { - s1r1set.push_back(idx); - } else { - s1r2set.push_back(idx); - } - if (point[dirs[dirIdx]] + distance_ <= pivot) { - continue; - } - side1Cand.push_back(idx); - } - for (unsigned int idx : side2) { - const auto &point = points[idx]; - assert(point[dirs[dirIdx]] > pivot && "Correctness Assertion"); - if (point[dirs[newDirIdx]] <= newPivot) { - s2r1set.push_back(idx); - } else { - s2r2set.push_back(idx); - } - if (point[dirs[dirIdx]] - distance_ >= pivot) { - continue; - } - side2Cand.push_back(idx); - } - - // Iterate over pairs of candidates - if (!side1Cand.empty() && !side2Cand.empty()) { - for (unsigned int &ci1 : side1Cand) { - for (unsigned int &ci2 : side2Cand) { - const auto &point1 = points[ci1]; - const auto &point2 = points[ci2]; - - assert(std::abs(point1[dirs[dirIdx]] - point2[dirs[dirIdx]]) <= - (2 * distance_) && - "Correctness Assertion"); - if (checkDistance(point1, point2)) { - pointNeighborhood_[ci1].push_back(ci2); - pointNeighborhood_[ci2].push_back(ci1); + } else { // D == 3 + for (int dx = -1; dx <= 1; ++dx) { + for (int dy = -1; dy <= 1; ++dy) { + for (int dz = -1; dz <= 1; ++dz) { + CellIndex neighbor = {cell[0] + dx, cell[1] + dy, cell[2] + dz}; + checkCellNeighbors(grid, points, idx, point, neighbor); } } } } - - // Recurse - if (side1.size() > 1) { - auto newS1Max = max; - newS1Max[dirs[dirIdx]] = pivot; // old diridx and old pivot! - createNeighborhood(points, s1r1set, s1r2set, min, newS1Max, newDirIdx, - dirs, newPivot); - } - if (side2.size() > 1) { - auto newS2Min = min; - newS2Min[dirs[dirIdx]] = pivot; // old diridx and old pivot! - createNeighborhood(points, s2r1set, s2r2set, newS2Min, max, newDirIdx, - dirs, newPivot); - } } template - bool checkDistance(const VectorType &p1, - const VectorType &p2) const { - for (int i = 0; i < D; ++i) { - if (std::abs(p1[i] - p2[i]) >= distance_) - return false; - } - if (Distance(p1, p2) < distance_) - return true; + void + checkCellNeighbors(const GridMap &grid, + const std::vector> &points, + unsigned int idx, + const VectorType &point, + const CellIndex &neighborCell) { + auto it = grid.find(neighborCell); + if (it == grid.end()) + return; - return false; + for (unsigned int otherIdx : it->second) { + // Only add each pair once: store neighbor only for otherIdx > idx + if (otherIdx <= idx) + continue; + if (Distance(point, points[otherIdx]) < distance_) { + pointNeighborhood_[idx].push_back(otherIdx); + pointNeighborhood_[otherIdx].push_back(idx); + } + } } }; + } // namespace viennaray diff --git a/tests/pointNeighborhood2D/pointNeighborhood2D.cpp b/tests/pointNeighborhood2D/pointNeighborhood2D.cpp index 7877bb6..009907f 100644 --- a/tests/pointNeighborhood2D/pointNeighborhood2D.cpp +++ b/tests/pointNeighborhood2D/pointNeighborhood2D.cpp @@ -26,6 +26,7 @@ int main() { geometry.initGeometry(device, points, normals, gridDelta); // setup simple 2D plane grid with normal in y-direction with discs only // overlapping at adjacent grid points x - x - x - x - x + // 0 - 1 - 2 - 3 - 4 // assert boundary points have 1 neighbor // assert inner points have 2 neighbors @@ -33,6 +34,10 @@ int main() { for (unsigned int idx = 0; idx < geometry.getNumPrimitives(); ++idx) { auto point = geometry.getPoint(idx); auto neighbors = geometry.getNeighborIndices(idx); + // std::printf("idx=%u point=(%f,%f,%f) neighbors=%zu\n", idx, point[0], + // point[1], point[2], neighbors.size()); + // for (auto n : neighbors) + // std::printf(" -> %u\n", n); if (std::fabs(point[0]) > 1 - eps) { // corner point VC_TEST_ASSERT(neighbors.size() == 1) From 6b11988efa17581ec19bcf7748231af24945a7b6 Mon Sep 17 00:00:00 2001 From: reiter Date: Wed, 18 Feb 2026 10:46:34 +0100 Subject: [PATCH 2/3] refactor: only 2D uses hash map --- include/viennaray/rayGeometryDisk.hpp | 3 +- include/viennaray/rayPointNeighborhood.hpp | 197 ++++++++++++++++-- tests/pointNeighborhood/pointNeighborhood.cpp | 5 + .../pointNeighborhood2D.cpp | 2 +- 4 files changed, 188 insertions(+), 19 deletions(-) diff --git a/include/viennaray/rayGeometryDisk.hpp b/include/viennaray/rayGeometryDisk.hpp index 040a62d..dbb691d 100644 --- a/include/viennaray/rayGeometryDisk.hpp +++ b/include/viennaray/rayGeometryDisk.hpp @@ -198,8 +198,7 @@ class GeometryDisk : public Geometry { return pointNeighborhood_.getNeighborIndices(idx); } - [[nodiscard]] PointNeighborhood const & - getPointNeighborhood() const { + [[nodiscard]] auto const &getPointNeighborhood() const { return pointNeighborhood_; } diff --git a/include/viennaray/rayPointNeighborhood.hpp b/include/viennaray/rayPointNeighborhood.hpp index 12e2e1e..bf49d8f 100644 --- a/include/viennaray/rayPointNeighborhood.hpp +++ b/include/viennaray/rayPointNeighborhood.hpp @@ -4,6 +4,7 @@ #include #include +#include #include namespace viennaray { @@ -15,6 +16,7 @@ template class PointNeighborhood { std::vector> pointNeighborhood_; NumericType distance_ = 0.; + NumericType distance2_ = 0.; // Cell index type using CellIndex = std::array; @@ -44,6 +46,7 @@ template class PointNeighborhood { static_assert(Dim >= static_cast(D), "Point dimension must be >= D"); distance_ = distance; + distance2_ = distance * distance; const auto numPoints = points.size(); pointNeighborhood_.clear(); pointNeighborhood_.resize(numPoints); @@ -51,24 +54,56 @@ template class PointNeighborhood { if (numPoints == 0 || distance_ <= 0) return; - const NumericType invCellSize = NumericType(1) / distance_; + if constexpr (D == 3) { + std::vector side1; + std::vector side2; - // Build the grid: map cell index -> list of point indices - GridMap grid; - grid.reserve(numPoints); - for (unsigned int idx = 0; idx < numPoints; ++idx) { - CellIndex cell = computeCell(points[idx], minCoords, invCellSize); - grid[cell].push_back(idx); - } + // create copy of bounding box + Vec3D min = minCoords; + Vec3D max = maxCoords; - // For each point, check all neighboring cells - for (unsigned int idx = 0; idx < numPoints; ++idx) { - const auto &point = points[idx]; - CellIndex cell = computeCell(point, minCoords, invCellSize); + std::vector dirs; + for (int i = 0; i < 3; ++i) { + if (min[i] != max[i]) { + dirs.push_back(i); + } + } + dirs.shrink_to_fit(); + + int dirIdx = 0; + NumericType pivot = (max[dirs[dirIdx]] + min[dirs[dirIdx]]) / 2; + + // divide point data + for (unsigned int idx = 0; idx < numPoints; ++idx) { + if (points[idx][dirs[dirIdx]] <= pivot) { + side1.push_back(idx); + } else { + side2.push_back(idx); + } + } + createNeighborhood(points, side1, side2, min, max, dirIdx, dirs, pivot); + } else { + const NumericType invCellSize = NumericType(1) / distance_; + + // Build the grid: map cell index -> list of point indices + GridMap grid; + grid.reserve(numPoints); + for (unsigned int idx = 0; idx < numPoints; ++idx) { + CellIndex cell = computeCell(points[idx], minCoords, invCellSize); + grid[cell].push_back(idx); + } + + // For each point, check all neighboring cells + for (unsigned int idx = 0; idx < numPoints; ++idx) { + const auto &point = points[idx]; + CellIndex cell = computeCell(point, minCoords, invCellSize); - // Iterate over the (2D: 3x3 = 9, 3D: 3x3x3 = 27) neighborhood of cells - iterateNeighborCells(grid, points, idx, point, cell); + // Iterate over the (2D: 3x3 = 9, 3D: 3x3x3 = 27) neighborhood of cells + iterateNeighborCells(grid, points, idx, point, cell); + } } + + assert(isUnique() && "Neighborhood contains duplicate entries"); } [[nodiscard]] std::vector const & @@ -135,12 +170,142 @@ template class PointNeighborhood { // Only add each pair once: store neighbor only for otherIdx > idx if (otherIdx <= idx) continue; - if (Distance(point, points[otherIdx]) < distance_) { + if (checkDistance(point, points[otherIdx])) { pointNeighborhood_[idx].push_back(otherIdx); pointNeighborhood_[otherIdx].push_back(idx); } } } -}; + void createNeighborhood(const std::vector> &points, + const std::vector &side1, + const std::vector &side2, + const Vec3D &min, + const Vec3D &max, const int dirIdx, + const std::vector &dirs, + const NumericType pivot) { + assert(0 <= dirIdx && dirIdx < dirs.size() && "Assumption"); + if (side1.size() + side2.size() <= 1) { + return; + } + + // Corner case + // The pivot element should actually be between min and max. + if (pivot == min[dirs[dirIdx]] || pivot == max[dirs[dirIdx]]) { + // In this case the points are extremely close to each other (with respect + // to the floating point precision). + assert((min[dirs[dirIdx]] + max[dirs[dirIdx]]) / 2 == pivot && + "Characterization of corner case"); + auto sides = std::vector(side1); + sides.insert(sides.end(), side2.begin(), side2.end()); + // Add each of them to the neighborhoods + for (unsigned int idx1 = 0; idx1 < sides.size() - 1; ++idx1) { + for (unsigned int idx2 = idx1 + 1; idx2 < sides.size(); ++idx2) { + auto const &pi1 = sides[idx1]; + auto const &pi2 = sides[idx2]; + assert(pi1 != pi2 && "Assumption"); + pointNeighborhood_[pi1].push_back(pi2); + pointNeighborhood_[pi2].push_back(pi1); + } + } + return; + } + + // sets of candidates + std::vector side1Cand; + std::vector side2Cand; + + int newDirIdx = (dirIdx + 1) % static_cast(dirs.size()); + NumericType newPivot = (max[dirs[newDirIdx]] + min[dirs[newDirIdx]]) / 2; + + // recursion sets + std::vector s1r1set; + std::vector s1r2set; + std::vector s2r1set; + std::vector s2r2set; + + for (unsigned int idx : side1) { + const auto &point = points[idx]; + assert(point[dirs[dirIdx]] <= pivot && "Correctness Assertion"); + if (point[dirs[newDirIdx]] <= newPivot) { + s1r1set.push_back(idx); + } else { + s1r2set.push_back(idx); + } + if (point[dirs[dirIdx]] + distance_ <= pivot) { + continue; + } + side1Cand.push_back(idx); + } + for (unsigned int idx : side2) { + const auto &point = points[idx]; + assert(point[dirs[dirIdx]] > pivot && "Correctness Assertion"); + if (point[dirs[newDirIdx]] <= newPivot) { + s2r1set.push_back(idx); + } else { + s2r2set.push_back(idx); + } + if (point[dirs[dirIdx]] - distance_ >= pivot) { + continue; + } + side2Cand.push_back(idx); + } + + // Iterate over pairs of candidates + if (!side1Cand.empty() && !side2Cand.empty()) { + for (unsigned int &ci1 : side1Cand) { + for (unsigned int &ci2 : side2Cand) { + const auto &point1 = points[ci1]; + const auto &point2 = points[ci2]; + + assert(std::abs(point1[dirs[dirIdx]] - point2[dirs[dirIdx]]) <= + (2 * distance_) && + "Correctness Assertion"); + if (checkDistance(point1, point2)) { + pointNeighborhood_[ci1].push_back(ci2); + pointNeighborhood_[ci2].push_back(ci1); + } + } + } + } + + // Recurse + if (side1.size() > 1) { + auto newS1Max = max; + newS1Max[dirs[dirIdx]] = pivot; // old diridx and old pivot! + createNeighborhood(points, s1r1set, s1r2set, min, newS1Max, newDirIdx, + dirs, newPivot); + } + if (side2.size() > 1) { + auto newS2Min = min; + newS2Min[dirs[dirIdx]] = pivot; // old diridx and old pivot! + createNeighborhood(points, s2r1set, s2r2set, newS2Min, max, newDirIdx, + dirs, newPivot); + } + } + + template + bool checkDistance(const VectorType &p1, + const VectorType &p2) const { + for (int i = 0; i < D; ++i) { + if (std::abs(p1[i] - p2[i]) > distance_) + return false; + } + if (Norm2(p1 - p2) <= distance2_) + return true; + + return false; + } + + bool isUnique() const { + for (const auto &neighbors : pointNeighborhood_) { + std::unordered_set uniqueNeighbors(neighbors.begin(), + neighbors.end()); + if (uniqueNeighbors.size() != neighbors.size()) { + return false; + } + } + return true; + } +}; } // namespace viennaray diff --git a/tests/pointNeighborhood/pointNeighborhood.cpp b/tests/pointNeighborhood/pointNeighborhood.cpp index 06182f2..c85e015 100644 --- a/tests/pointNeighborhood/pointNeighborhood.cpp +++ b/tests/pointNeighborhood/pointNeighborhood.cpp @@ -23,7 +23,12 @@ int main() { auto device = rtcNewDevice(""); GeometryDisk geometry; + Timer timer; + timer.start(); geometry.initGeometry(device, points, normals, gridDelta - eps); + timer.finish(); + std::cout << "Point neighborhood initialization time: " + << timer.currentDuration * 1e-6 << " ms\n"; auto bdBox = geometry.getBoundingBox(); for (unsigned int idx = 0; idx < geometry.getNumPrimitives(); ++idx) { diff --git a/tests/pointNeighborhood2D/pointNeighborhood2D.cpp b/tests/pointNeighborhood2D/pointNeighborhood2D.cpp index 009907f..c92663b 100644 --- a/tests/pointNeighborhood2D/pointNeighborhood2D.cpp +++ b/tests/pointNeighborhood2D/pointNeighborhood2D.cpp @@ -23,7 +23,7 @@ int main() { auto device = rtcNewDevice(""); GeometryDisk geometry; - geometry.initGeometry(device, points, normals, gridDelta); + geometry.initGeometry(device, points, normals, gridDelta - eps); // setup simple 2D plane grid with normal in y-direction with discs only // overlapping at adjacent grid points x - x - x - x - x // 0 - 1 - 2 - 3 - 4 From ec443f957fc207b5df5cc685c4a4025794e41310 Mon Sep 17 00:00:00 2001 From: reiter Date: Fri, 20 Feb 2026 09:38:07 +0100 Subject: [PATCH 3/3] chore: bump version --- CMakeLists.txt | 2 +- README.md | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e803622..6121f29 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.20 FATAL_ERROR) project( ViennaRay LANGUAGES CXX - VERSION 3.11.0) + VERSION 3.11.1) # -------------------------------------------------------------------------------------------------------- # Library switches diff --git a/README.md b/README.md index 907815f..19760ce 100644 --- a/README.md +++ b/README.md @@ -63,7 +63,7 @@ We recommend using [CPM.cmake](https://github.com/cpm-cmake/CPM.cmake) to consum * Installation with CPM ```cmake - CPMAddPackage("gh:viennatools/viennaray@3.11.0") # Use the latest release version + CPMAddPackage("gh:viennatools/viennaray@3.11.1") # Use the latest release version ``` * With a local installation