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
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<NumericType, D> {
     return pointNeighborhood_.getNeighborIndices(idx);
   }
 
-  [[nodiscard]] PointNeighborhood<NumericType, D> const &
-  getPointNeighborhood() const {
+  [[nodiscard]] auto const &getPointNeighborhood() const {
     return pointNeighborhood_;
   }
 
diff --git a/include/viennaray/rayPointNeighborhood.hpp b/include/viennaray/rayPointNeighborhood.hpp
index e60fed1..bf49d8f 100644
--- a/include/viennaray/rayPointNeighborhood.hpp
+++ b/include/viennaray/rayPointNeighborhood.hpp
@@ -3,6 +3,8 @@
 #include <vcVectorType.hpp>
 
 #include <cassert>
+#include <unordered_map>
+#include <unordered_set>
 #include <vector>
 
 namespace viennaray {
@@ -10,8 +12,29 @@ namespace viennaray {
 using namespace viennacore;
 
 template <typename NumericType, int D> class PointNeighborhood {
+  static_assert(D == 2 || D == 3, "Only 2D and 3D are supported");
+
   std::vector<std::vector<unsigned int>> pointNeighborhood_;
   NumericType distance_ = 0.;
+  NumericType distance2_ = 0.;
+
+  // Cell index type
+  using CellIndex = std::array<int, D>;
+
+  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<size_t>(c[i]);
+        h *= 16777619u;
+      }
+      return h;
+    }
+  };
+
+  using GridMap =
+      std::unordered_map<CellIndex, std::vector<unsigned int>, CellIndexHash>;
 
 public:
   PointNeighborhood() = default;
@@ -20,10 +43,17 @@ template <typename NumericType, int D> class PointNeighborhood {
   void init(std::vector<VectorType<NumericType, Dim>> const &points,
             NumericType distance, Vec3D<NumericType> const &minCoords,
             Vec3D<NumericType> const &maxCoords) {
+    static_assert(Dim >= static_cast<size_t>(D),
+                  "Point dimension must be >= D");
     distance_ = distance;
+    distance2_ = distance * distance;
     const auto numPoints = points.size();
     pointNeighborhood_.clear();
-    pointNeighborhood_.resize(numPoints, std::vector<unsigned int>{});
+    pointNeighborhood_.resize(numPoints);
+
+    if (numPoints == 0 || distance_ <= 0)
+      return;
+
     if constexpr (D == 3) {
       std::vector<unsigned int> side1;
       std::vector<unsigned int> side2;
@@ -53,16 +83,27 @@ template <typename NumericType, int D> class PointNeighborhood {
       }
       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);
       }
     }
+
+    assert(isUnique() && "Neighborhood contains duplicate entries");
   }
 
   [[nodiscard]] std::vector<unsigned int> const &
@@ -78,13 +119,71 @@ template <typename NumericType, int D> class PointNeighborhood {
   [[nodiscard]] NumericType getDistance() const { return distance_; }
 
 private:
+  template <size_t Dim>
+  CellIndex computeCell(const VectorType<NumericType, Dim> &point,
+                        const Vec3D<NumericType> &minCoords,
+                        NumericType invCellSize) const {
+    CellIndex cell;
+    for (int i = 0; i < D; ++i) {
+      cell[i] =
+          static_cast<int>(std::floor((point[i] - minCoords[i]) * invCellSize));
+    }
+    return cell;
+  }
+
+  template <size_t Dim>
+  void iterateNeighborCells(
+      const GridMap &grid,
+      const std::vector<VectorType<NumericType, Dim>> &points, unsigned int idx,
+      const VectorType<NumericType, Dim> &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);
+        }
+      }
+    } 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);
+          }
+        }
+      }
+    }
+  }
+
+  template <size_t Dim>
+  void
+  checkCellNeighbors(const GridMap &grid,
+                     const std::vector<VectorType<NumericType, Dim>> &points,
+                     unsigned int idx,
+                     const VectorType<NumericType, Dim> &point,
+                     const CellIndex &neighborCell) {
+    auto it = grid.find(neighborCell);
+    if (it == grid.end())
+      return;
+
+    for (unsigned int otherIdx : it->second) {
+      // Only add each pair once: store neighbor only for otherIdx > idx
+      if (otherIdx <= idx)
+        continue;
+      if (checkDistance(point, points[otherIdx])) {
+        pointNeighborhood_[idx].push_back(otherIdx);
+        pointNeighborhood_[otherIdx].push_back(idx);
+      }
+    }
+  }
+
   void createNeighborhood(const std::vector<Vec3D<NumericType>> &points,
                           const std::vector<unsigned int> &side1,
                           const std::vector<unsigned int> &side2,
                           const Vec3D<NumericType> &min,
-                          const Vec3D<NumericType> &max, const int &dirIdx,
+                          const Vec3D<NumericType> &max, const int dirIdx,
                           const std::vector<int> &dirs,
-                          const NumericType &pivot) {
+                          const NumericType pivot) {
     assert(0 <= dirIdx && dirIdx < dirs.size() && "Assumption");
     if (side1.size() + side2.size() <= 1) {
       return;
@@ -189,13 +288,24 @@ template <typename NumericType, int D> class PointNeighborhood {
   bool checkDistance(const VectorType<NumericType, Dim> &p1,
                      const VectorType<NumericType, Dim> &p2) const {
     for (int i = 0; i < D; ++i) {
-      if (std::abs(p1[i] - p2[i]) >= distance_)
+      if (std::abs(p1[i] - p2[i]) > distance_)
         return false;
     }
-    if (Distance(p1, p2) < distance_)
+    if (Norm2(p1 - p2) <= distance2_)
       return true;
 
     return false;
   }
+
+  bool isUnique() const {
+    for (const auto &neighbors : pointNeighborhood_) {
+      std::unordered_set<unsigned int> 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<NumericType, D> 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 7877bb6..c92663b 100644
--- a/tests/pointNeighborhood2D/pointNeighborhood2D.cpp
+++ b/tests/pointNeighborhood2D/pointNeighborhood2D.cpp
@@ -23,9 +23,10 @@ int main() {
 
   auto device = rtcNewDevice("");
   GeometryDisk<NumericType, D> 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
 
   // 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)