diff --git a/CMakeLists.txt b/CMakeLists.txt
index cbd9d30..27c682e 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.1.0)
+  VERSION 3.1.3)
 
 # --------------------------------------------------------------------------------------------------------
 # Library switches
diff --git a/README.md b/README.md
index 0e645db..3d1303b 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.1.1")
+  CPMAddPackage("gh:viennatools/viennaray@3.1.3")
   ```
 
 * With a local installation
diff --git a/include/viennaray/rayGeometry.hpp b/include/viennaray/rayGeometry.hpp
index fb6f4a6..26b4202 100644
--- a/include/viennaray/rayGeometry.hpp
+++ b/include/viennaray/rayGeometry.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <rayPointNeighborhood.hpp>
 #include <rayUtil.hpp>
 
 namespace viennaray {
@@ -95,10 +96,13 @@ template <typename NumericType, int D> class Geometry {
     assert(rtcGetDeviceError(device) == RTC_ERROR_NONE &&
            "RTC Error: rtcCommitGeometry");
 
-    initPointNeighborhood(points);
     if (materialIds_.empty()) {
       materialIds_.resize(numPoints_, 0);
     }
+
+    // Initialize point neighborhood
+    pointNeighborhood_ = PointNeighborhood<NumericType, D>(
+        points, 2 * discRadii_, minCoords_, maxCoords_);
   }
 
   template <typename MatIdType>
@@ -124,8 +128,13 @@ template <typename NumericType, int D> class Geometry {
 
   [[nodiscard]] std::vector<unsigned int> const &
   getNeighborIndicies(const unsigned int idx) const {
-    assert(idx < numPoints_ && "Geometry: Index out of bounds");
-    return pointNeighborhood_[idx];
+    assert(pointNeighborhood_.getDistance() > 0.); // check if initialized
+    return pointNeighborhood_.getNeighborIndicies(idx);
+  }
+
+  [[nodiscard]] PointNeighborhood<NumericType, D> const &
+  getPointNeighborhood() const {
+    return pointNeighborhood_;
   }
 
   [[nodiscard]] size_t getNumPoints() const { return numPoints_; }
@@ -186,175 +195,6 @@ template <typename NumericType, int D> class Geometry {
     }
   }
 
-private:
-  template <size_t Dim>
-  void initPointNeighborhood(
-      std::vector<std::array<NumericType, Dim>> const &points) {
-    pointNeighborhood_.clear();
-    pointNeighborhood_.resize(numPoints_, std::vector<unsigned int>{});
-
-    if constexpr (D == 3) {
-      std::vector<unsigned int> side1;
-      std::vector<unsigned int> side2;
-
-      // create copy of bounding box
-      Vec3D<NumericType> min = minCoords_;
-      Vec3D<NumericType> max = maxCoords_;
-
-      std::vector<int> 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 {
-      // 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], 2 * discRadii_)) {
-            pointNeighborhood_[idx1].push_back(idx2);
-            pointNeighborhood_[idx2].push_back(idx1);
-          }
-        }
-      }
-    }
-  }
-
-  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 std::vector<int> &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<unsigned int>(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<unsigned int> side1Cand;
-    std::vector<unsigned int> side2Cand;
-
-    int newDirIdx = (dirIdx + 1) % dirs.size();
-    NumericType newPivot = (max[dirs[newDirIdx]] + min[dirs[newDirIdx]]) / 2;
-
-    // recursion sets
-    std::vector<unsigned int> s1r1set;
-    std::vector<unsigned int> s1r2set;
-    std::vector<unsigned int> s2r1set;
-    std::vector<unsigned int> s2r2set;
-
-    for (unsigned int idx = 0; idx < side1.size(); ++idx) {
-      const auto &point = points[side1[idx]];
-      assert(point[dirs[dirIdx]] <= pivot && "Correctness Assertion");
-      if (point[dirs[newDirIdx]] <= newPivot) {
-        s1r1set.push_back(side1[idx]);
-      } else {
-        s1r2set.push_back(side1[idx]);
-      }
-      if (point[dirs[dirIdx]] + 2 * discRadii_ <= pivot) {
-        continue;
-      }
-      side1Cand.push_back(side1[idx]);
-    }
-    for (unsigned int idx = 0; idx < side2.size(); ++idx) {
-      const auto &point = points[side2[idx]];
-      assert(point[dirs[dirIdx]] > pivot && "Correctness Assertion");
-      if (point[dirs[newDirIdx]] <= newPivot) {
-        s2r1set.push_back(side2[idx]);
-      } else {
-        s2r2set.push_back(side2[idx]);
-      }
-      if (point[dirs[dirIdx]] - 2 * discRadii_ >= pivot) {
-        continue;
-      }
-      side2Cand.push_back(side2[idx]);
-    }
-
-    // Iterate over pairs of candidates
-    if (side1Cand.size() > 0 && side2Cand.size() > 0) {
-      for (unsigned int ci1 = 0; ci1 < side1Cand.size(); ++ci1) {
-        for (unsigned int ci2 = 0; ci2 < side2Cand.size(); ++ci2) {
-          const auto &point1 = points[side1Cand[ci1]];
-          const auto &point2 = points[side2Cand[ci2]];
-
-          assert(std::abs(point1[dirs[dirIdx]] - point2[dirs[dirIdx]]) <=
-                     (4 * discRadii_) &&
-                 "Correctness Assertion");
-          if (checkDistance(point1, point2, 2 * discRadii_)) {
-            pointNeighborhood_[side1Cand[ci1]].push_back(side2Cand[ci2]);
-            pointNeighborhood_[side2Cand[ci2]].push_back(side1Cand[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 <size_t Dim>
-  static bool checkDistance(const std::array<NumericType, Dim> &p1,
-                            const std::array<NumericType, Dim> &p2,
-                            const NumericType &dist) {
-    for (int i = 0; i < D; ++i) {
-      if (std::abs(p1[i] - p2[i]) >= dist)
-        return false;
-    }
-    if (Distance(p1, p2) < dist)
-      return true;
-
-    return false;
-  }
-
 private:
   // "RTC_GEOMETRY_TYPE_POINT:
   // The vertex buffer stores each control vertex in the form of a single
@@ -382,8 +222,9 @@ template <typename NumericType, int D> class Geometry {
   NumericType discRadii_;
   Vec3D<NumericType> minCoords_;
   Vec3D<NumericType> maxCoords_;
-  pointNeighborhoodType pointNeighborhood_;
   std::vector<int> materialIds_;
+
+  PointNeighborhood<NumericType, D> pointNeighborhood_;
 };
 
 } // namespace viennaray
diff --git a/include/viennaray/rayParticle.hpp b/include/viennaray/rayParticle.hpp
index c38f7d9..1e17375 100644
--- a/include/viennaray/rayParticle.hpp
+++ b/include/viennaray/rayParticle.hpp
@@ -134,12 +134,14 @@ class TestParticle : public Particle<TestParticle<NumericType>, NumericType> {
                         const unsigned int primID, const int materialId,
                         TracingData<NumericType> &localData,
                         const TracingData<NumericType> *globalData,
-                        RNG &rngState) override final {}
+                        RNG &rngState) override final {
+    localData.getVectorData(0)[primID] += rayWeight;
+  }
 
   NumericType getSourceDistributionPower() const override final { return 1.; }
 
   std::vector<std::string> getLocalDataLabels() const override final {
-    return {};
+    return {"testFlux"};
   }
 
   void logData(DataLog<NumericType> &log) override final {}
diff --git a/include/viennaray/rayPointNeighborhood.hpp b/include/viennaray/rayPointNeighborhood.hpp
new file mode 100644
index 0000000..741db32
--- /dev/null
+++ b/include/viennaray/rayPointNeighborhood.hpp
@@ -0,0 +1,201 @@
+#pragma once
+
+#include <vcVectorUtil.hpp>
+
+#include <cassert>
+#include <vector>
+
+namespace viennaray {
+
+using namespace viennacore;
+
+template <typename NumericType, int D> class PointNeighborhood {
+  std::vector<std::vector<unsigned int>> pointNeighborhood_;
+  NumericType distance_ = 0.;
+
+public:
+  PointNeighborhood() = default;
+
+  template <size_t Dim>
+  PointNeighborhood(std::vector<std::array<NumericType, Dim>> const &points,
+                    NumericType distance, Vec3D<NumericType> const &minCoords,
+                    Vec3D<NumericType> const &maxCoords)
+      : distance_(distance) {
+
+    const auto numPoints = points.size();
+    pointNeighborhood_.resize(numPoints, std::vector<unsigned int>{});
+    if constexpr (D == 3) {
+      std::vector<unsigned int> side1;
+      std::vector<unsigned int> side2;
+
+      // create copy of bounding box
+      Vec3D<NumericType> min = minCoords;
+      Vec3D<NumericType> max = maxCoords;
+
+      std::vector<int> 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 {
+      /// 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);
+          }
+        }
+      }
+    }
+  }
+
+  [[nodiscard]] std::vector<unsigned int> const &
+  getNeighborIndicies(const unsigned int idx) const {
+    assert(idx < pointNeighborhood_.size() && "Index out of bounds");
+    return pointNeighborhood_[idx];
+  }
+
+  [[nodiscard]] size_t getNumPoints() const {
+    return pointNeighborhood_.size();
+  }
+
+  [[nodiscard]] NumericType getDistance() const { return distance_; }
+
+private:
+  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 std::vector<int> &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<unsigned int>(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<unsigned int> side1Cand;
+    std::vector<unsigned int> side2Cand;
+
+    int newDirIdx = (dirIdx + 1) % dirs.size();
+    NumericType newPivot = (max[dirs[newDirIdx]] + min[dirs[newDirIdx]]) / 2;
+
+    // recursion sets
+    std::vector<unsigned int> s1r1set;
+    std::vector<unsigned int> s1r2set;
+    std::vector<unsigned int> s2r1set;
+    std::vector<unsigned int> s2r2set;
+
+    for (unsigned int idx = 0; idx < side1.size(); ++idx) {
+      const auto &point = points[side1[idx]];
+      assert(point[dirs[dirIdx]] <= pivot && "Correctness Assertion");
+      if (point[dirs[newDirIdx]] <= newPivot) {
+        s1r1set.push_back(side1[idx]);
+      } else {
+        s1r2set.push_back(side1[idx]);
+      }
+      if (point[dirs[dirIdx]] + distance_ <= pivot) {
+        continue;
+      }
+      side1Cand.push_back(side1[idx]);
+    }
+    for (unsigned int idx = 0; idx < side2.size(); ++idx) {
+      const auto &point = points[side2[idx]];
+      assert(point[dirs[dirIdx]] > pivot && "Correctness Assertion");
+      if (point[dirs[newDirIdx]] <= newPivot) {
+        s2r1set.push_back(side2[idx]);
+      } else {
+        s2r2set.push_back(side2[idx]);
+      }
+      if (point[dirs[dirIdx]] - distance_ >= pivot) {
+        continue;
+      }
+      side2Cand.push_back(side2[idx]);
+    }
+
+    // Iterate over pairs of candidates
+    if (side1Cand.size() > 0 && side2Cand.size() > 0) {
+      for (unsigned int ci1 = 0; ci1 < side1Cand.size(); ++ci1) {
+        for (unsigned int ci2 = 0; ci2 < side2Cand.size(); ++ci2) {
+          const auto &point1 = points[side1Cand[ci1]];
+          const auto &point2 = points[side2Cand[ci2]];
+
+          assert(std::abs(point1[dirs[dirIdx]] - point2[dirs[dirIdx]]) <=
+                     (2 * distance_) &&
+                 "Correctness Assertion");
+          if (checkDistance(point1, point2)) {
+            pointNeighborhood_[side1Cand[ci1]].push_back(side2Cand[ci2]);
+            pointNeighborhood_[side2Cand[ci2]].push_back(side1Cand[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 <size_t Dim>
+  bool checkDistance(const std::array<NumericType, Dim> &p1,
+                     const std::array<NumericType, Dim> &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;
+
+    return false;
+  }
+};
+} // namespace viennaray
\ No newline at end of file
diff --git a/include/viennaray/rayTrace.hpp b/include/viennaray/rayTrace.hpp
index 9b2614d..6ee0f9a 100644
--- a/include/viennaray/rayTrace.hpp
+++ b/include/viennaray/rayTrace.hpp
@@ -259,17 +259,46 @@ template <class NumericType, int D> class Trace {
 
   /// Helper function to smooth the recorded flux by averaging over the
   /// neighborhood in a post-processing step.
-  void smoothFlux(std::vector<NumericType> &flux) {
+  void smoothFlux(std::vector<NumericType> &flux, int numNeighbors = 1) {
     assert(flux.size() == geometry_.getNumPoints() &&
            "Unequal number of points in smoothFlux");
     auto oldFlux = flux;
+    PointNeighborhood<NumericType, D> pointNeighborhood;
+    if (numNeighbors == 1) {
+      // re-use the neighborhood from the geometry
+      pointNeighborhood = geometry_.getPointNeighborhood();
+    } else {
+      // create a new neighborhood with a larger radius
+      auto boundingBox = geometry_.getBoundingBox();
+      std::vector<Vec3D<NumericType>> points(geometry_.getNumPoints());
+#pragma omp parallel for
+      for (int idx = 0; idx < geometry_.getNumPoints(); idx++) {
+        points[idx] = geometry_.getPoint(idx);
+      }
+      pointNeighborhood = PointNeighborhood<NumericType, D>(
+          points, numNeighbors * 2 * diskRadius_, boundingBox[0],
+          boundingBox[1]);
+    }
+
 #pragma omp parallel for
     for (int idx = 0; idx < geometry_.getNumPoints(); idx++) {
-      auto neighborhood = geometry_.getNeighborIndicies(idx);
+
+      NumericType vv = oldFlux[idx];
+
+      auto const &neighborhood = pointNeighborhood.getNeighborIndicies(idx);
+      NumericType sum = 1.;
+      auto const normal = geometry_.getPrimNormal(idx);
+
       for (auto const &nbi : neighborhood) {
-        flux[idx] += oldFlux[nbi];
+        auto nnormal = geometry_.getPrimNormal(nbi);
+        auto weight = DotProduct(normal, nnormal);
+        if (weight > 0.) {
+          vv += oldFlux[nbi] * weight;
+          sum += weight;
+        }
       }
-      flux[idx] /= (neighborhood.size() + 1);
+
+      flux[idx] = vv / sum;
     }
   }
 
diff --git a/tests/diskAreas/diskAreas.cpp b/tests/diskAreas/diskAreas.cpp
index cca4128..c99be55 100644
--- a/tests/diskAreas/diskAreas.cpp
+++ b/tests/diskAreas/diskAreas.cpp
@@ -47,6 +47,9 @@ int main() {
 
   TestParticle<NumericType> particle;
   auto cp = particle.clone();
+  localData.setNumberOfVectorData(cp->getLocalDataLabels().size());
+  auto numPoints = geometry.getNumPoints();
+  localData.resizeAllVectorData(numPoints, 0.);
 
   DataLog<NumericType> log;
   TraceInfo info;
diff --git a/tests/traceInterface/traceInterface.cpp b/tests/traceInterface/traceInterface.cpp
index 6060d08..5dcb5c4 100644
--- a/tests/traceInterface/traceInterface.cpp
+++ b/tests/traceInterface/traceInterface.cpp
@@ -55,6 +55,13 @@ int main() {
 
   rayTracer.apply();
 
+  auto flux = rayTracer.getLocalData().getVectorData(0);
+  VC_TEST_ASSERT(flux.size() == points.size());
+
+  rayTracer.normalizeFlux(flux);
+  rayTracer.smoothFlux(flux, 2);
+  VC_TEST_ASSERT(flux.size() == points.size());
+
   auto info = rayTracer.getRayTraceInfo();
   VC_TEST_ASSERT(info.numRays == 4410);