Visibilities calculated in any direction

include/viennals/lsCalculateVisibilities.hpp

@@ -1,99 +1,136 @@
 #pragma once
 
 #include <lsDomain.hpp>
+#include <vcVectorUtil.hpp>
 
 namespace viennals {
 using namespace viennacore;
 
 template <class NumericType, int D> class CalculateVisibilities {
   SmartPointer<Domain<NumericType, D>> levelSet;
-  int openBoundaryDirection;
-  bool isOpenBoundaryNegative;
+  Vec3D<NumericType> direction;
+  const NumericType epsilon = static_cast<NumericType>(1e-6);
 
 public:
   static constexpr char visibilitiesLabel[] = "Visibilities";
 
   CalculateVisibilities(
       const SmartPointer<Domain<NumericType, D>> &passedLevelSet,
-      int passedOpenBoundaryDirection = D - 1,
-      bool passedIsOpenBoundaryNegative = false)
-      : levelSet(passedLevelSet),
-        openBoundaryDirection(passedOpenBoundaryDirection),
-        isOpenBoundaryNegative(passedIsOpenBoundaryNegative) {}
+      const Vec3D<NumericType> passedDirection)
+      : levelSet(passedLevelSet), direction(passedDirection) {}
 
   void apply() {
 
     auto &domain = levelSet->getDomain();
     auto &grid = levelSet->getGrid();
-    const NumericType max = std::numeric_limits<NumericType>::max();
 
-    auto numDefinedPoints = domain.getNumberOfPoints();
-    std::vector<NumericType> visibilities(numDefinedPoints);
-
-    std::vector<hrleIndexType> oldIndices(
-        D - 1 - openBoundaryDirection,
-        std::numeric_limits<hrleIndexType>::max());
-
-    unsigned int size = 1;
-    for (int i = 0; i < openBoundaryDirection; ++i) {
-      assert(!grid.isPosBoundaryInfinite(i));
-      assert(!grid.isNegBoundaryInfinite(i));
-
-      size *= (grid.getMaxIndex(i) - grid.getMinIndex(i) + 1);
+    // *** Determine extents of domain ***
+    Vec3D<NumericType> minDefinedPoint;
+    Vec3D<NumericType> maxDefinedPoint;
+    // Initialize with extreme values
+    for (int i = 0; i < D; ++i) {
+        minDefinedPoint[i] = std::numeric_limits<NumericType>::max();
+        maxDefinedPoint[i] = std::numeric_limits<NumericType>::lowest();
     }
-
-    std::vector<NumericType> minValues(size, max);
-
-    hrleSizeType id = 0;
-
-    hrleSparseIterator<typename Domain<NumericType, D>::DomainType> it(
-        domain, !isOpenBoundaryNegative);
-    while (!it.isFinished()) {
-
-      for (int i = 0; i < D - 1 - openBoundaryDirection; ++i) {
-        bool b = false;
-        if (oldIndices[i] !=
-            it.getStartIndices(i + openBoundaryDirection + 1)) {
-          oldIndices[i] = it.getStartIndices(i + openBoundaryDirection + 1);
-          b = true;
-        }
-        if (b)
-          minValues.assign(size, max);
+    // Iterate through all defined points in the domain
+    for (hrleSparseIterator<typename Domain<NumericType, D>::DomainType> it(domain); !it.isFinished(); it.next()) {
+      if (!it.isDefined()) continue; // Skip undefined points
+
+      // Get the coordinate of the current point
+      auto point = it.getStartIndices();
+      for (int i = 0; i < D; ++i) {
+        // Compare to update min and max defined points
+        NumericType coord = point[i]; // * grid.getGridDelta();
+        minDefinedPoint[i] = std::min(minDefinedPoint[i], coord);
+        maxDefinedPoint[i] = std::max(maxDefinedPoint[i], coord);
       }
+    }
+    //****************************
 
-      unsigned int pos_begin = 0;
-      unsigned int pos_end = 0;
+    // Invert the vector
+    auto dir = Normalize(Inv(direction));
 
-      for (int i = openBoundaryDirection - 1; i >= 0; --i) {
-        pos_begin *= (grid.getMaxIndex(i) - grid.getMinIndex(i) + 1);
-        pos_end *= (grid.getMaxIndex(i) - grid.getMinIndex(i) + 1);
-        pos_begin += (it.getStartIndices(i) - grid.getMinIndex(i));
-        pos_end += (it.getEndIndices(i) - grid.getMinIndex(i));
-      }
-
-      if (it.isDefined()) {
-        visibilities[isOpenBoundaryNegative ? id
-                                            : (numDefinedPoints - 1 - id)] =
-            (it.getValue() < minValues.at(pos_begin)) ? 1. : 0.;
-        ++id;
-      }
+    auto numDefinedPoints = domain.getNumberOfPoints();
+    std::vector<NumericType> visibilities(numDefinedPoints);
 
-      for (unsigned int i = pos_begin; i <= pos_end; ++i)
-        minValues.at(i) = std::min(minValues.at(i), it.getValue());
+      // std::vector<Vec3D<hrleIndexType>> visitedCells;
+
+      hrleSizeType id = 0;
+      hrleSparseIterator<typename Domain<NumericType, D>::DomainType> it(domain);
+
+      while (!it.isFinished()) {
+          if (!it.isDefined()) {
+              it.next();
+              continue;
+          }
+
+          // Starting position of the point
+          Vec3D<NumericType> currentPos;
+          for (int i = 0; i < D; ++i) {
+              currentPos[i] = it.getStartIndices(i);
+          }
+
+          // Start tracing the ray
+          NumericType minLevelSetValue = it.getValue(); // Starting level set value
+          Vec3D<NumericType> rayPos = currentPos;
+          bool visibility = true;
+
+          while(1) {
+            // Update the ray position
+            for (int i = 0; i < D; ++i) {
+              rayPos[i] += dir[i];
+            }
+
+            // Determine the nearest grid cell (round to nearest index)
+            Vec3D<hrleIndexType> nearestCell;
+            for (int i = 0; i < D; ++i) {
+              nearestCell[i] = static_cast<hrleIndexType>(rayPos[i]);
+            }
+
+            // // Before adding a cell, check if it's already visited
+            // if (std::find(visitedCells.begin(), visitedCells.end(), nearestCell) == visitedCells.end()) {
+            //     visitedCells.push_back(nearestCell);
+            // }
+
+            // Check if the nearest cell is within bounds
+            bool outOfBounds = false;
+            for (int i = 0; i < D; ++i) {
+              if (nearestCell[i] < minDefinedPoint[i] || nearestCell[i] > maxDefinedPoint[i]) {
+                outOfBounds = true;
+                break;
+              }
+            }
+
+            if (outOfBounds) {
+              break; // Ray is outside the grid
+            }
+
+            // Access the level set value at the nearest cell
+            NumericType neighborValue = std::numeric_limits<NumericType>::max();
+            hrleSparseIterator<typename Domain<NumericType, D>::DomainType> neighborIt(domain);
+            neighborIt.goToIndices(nearestCell);
+            if (neighborIt.isDefined()) {
+              neighborValue = neighborIt.getValue();
+            }
+
+            // Update the minimum value encountered
+            if (neighborValue < minLevelSetValue) {
+              visibility = false;
+              break;
+            }
+        }
 
-      if (isOpenBoundaryNegative) {
+        // Update visibility for this point
+        visibilities[id++] = visibility ? 1.0 : 0.0;
         it.next();
-      } else {
-        it.previous();
       }
-    }
 
     auto &pointData = levelSet->getPointData();
     // delete if already exists
-    if (int i = pointData.getScalarDataIndex(visibilitiesLabel); i != -1) {
+    if (int i = pointData.getScalarDataIndex("Visibilities"); i != -1) {
       pointData.eraseScalarData(i);
     }
-    pointData.insertNextScalarData(visibilities, visibilitiesLabel);
+    pointData.insertNextScalarData(visibilities, "Visibilities");
 
     assert(id == numDefinedPoints);
   }