Add visibility test and parallelize

ViennaTools · Nov 26, 2024 · cc26206 · cc26206
1 parent 2fc81d7
commit cc26206
Show file tree

Hide file tree

Showing 3 changed files with 122 additions and 60 deletions.
diff --git a/include/viennals/lsCalculateVisibilities.hpp b/include/viennals/lsCalculateVisibilities.hpp
@@ -57,79 +57,98 @@ template <class NumericType, int D> class CalculateVisibilities {
     auto numDefinedPoints = domain.getNumberOfPoints();
     std::vector<NumericType> visibilities(numDefinedPoints);
 
-    // 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;
+#pragma omp parallel num_threads(levelSet->getNumberOfSegments())
+    {
+      int p = 0;
+#ifdef _OPENMP
+      p = omp_get_thread_num();
+#endif
+
+      hrleVectorType<hrleIndexType, D> startVector =
+          (p == 0) ? grid.getMinGridPoint() : domain.getSegmentation()[p - 1];
+
+      hrleVectorType<hrleIndexType, D> endVector =
+          (p != static_cast<int>(domain.getNumberOfSegments() - 1))
+              ? domain.getSegmentation()[p]
+              : grid.incrementIndices(grid.getMaxGridPoint());
+
+      // Calculate the starting index for the visibility vector
+      hrleSizeType id = 0;
+      for (int i = 0; i < p; ++i) {
+        id += domain.getDomainSegment(i).getNumberOfPoints();
       }
 
-      // Starting position of the point
-      Vec3D<NumericType> currentPos;
-      for (int i = 0; i < D; ++i) {
-        currentPos[i] = it.getStartIndices(i);
-      }
+      for (hrleSparseIterator<typename Domain<NumericType, D>::DomainType> it(
+               domain, startVector);
+           it.getStartIndices() < endVector; ++it) {
 
-      // Start tracing the ray
-      NumericType minLevelSetValue = it.getValue(); // Starting level set value
-      Vec3D<NumericType> rayPos = currentPos;
-      bool visibility = true;
+        if (!it.isDefined())
+          continue;
 
-      while (1) {
-        // Update the ray position
+        // Starting position of the point
+        Vec3D<NumericType> currentPos;
         for (int i = 0; i < D; ++i) {
-          rayPos[i] += dir[i];
+          currentPos[i] = it.getStartIndices(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]);
-        }
+        // Start tracing the ray
+        NumericType minLevelSetValue =
+            it.getValue(); // Starting level set value
+        Vec3D<NumericType> rayPos = currentPos;
+        bool visibility = true;
 
-        // // Before adding a cell, check if it's already visited
-        // if (std::find(visitedCells.begin(), visitedCells.end(), nearestCell)
-        // == visitedCells.end()) {
-        //     visitedCells.push_back(nearestCell);
-        // }
+        while (1) {
+          // Update the ray position
+          for (int i = 0; i < D; ++i) {
+            rayPos[i] += dir[i];
+          }
 
-        // 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;
+          // 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]);
           }
-        }
 
-        if (outOfBounds) {
-          break; // Ray is outside the grid
-        }
+          // // 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;
+            }
+          }
 
-        // 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();
-        }
+          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;
+          // Update the minimum value encountered
+          if (neighborValue < minLevelSetValue) {
+            visibility = false;
+            break;
+          }
         }
-      }
 
-      // Update visibility for this point
-      visibilities[id++] = visibility ? 1.0 : 0.0;
-      it.next();
+        // Update visibility for this point
+        visibilities[id++] = visibility ? 1.0 : 0.0;
+      }
     }
 
     auto &pointData = levelSet->getPointData();
@@ -138,8 +157,6 @@ template <class NumericType, int D> class CalculateVisibilities {
       pointData.eraseScalarData(i);
     }
     pointData.insertNextScalarData(visibilities, visibilitiesLabel);
-
-    assert(id == numDefinedPoints);
   }
 };
 

diff --git a/tests/Visibility/CMakeLists.txt b/tests/Visibility/CMakeLists.txt
@@ -0,0 +1,7 @@
+project(Visibility LANGUAGES CXX)
+
+add_executable(${PROJECT_NAME} "${PROJECT_NAME}.cpp")
+target_link_libraries(${PROJECT_NAME} PRIVATE ViennaLS)
+
+add_dependencies(ViennaLS_Tests ${PROJECT_NAME})
+add_test(NAME ${PROJECT_NAME} COMMAND $<TARGET_FILE:${PROJECT_NAME}>)
diff --git a/tests/Visibility/Visibility.cpp b/tests/Visibility/Visibility.cpp
@@ -0,0 +1,38 @@
+#include <chrono>
+#include <iostream>
+
+#include <lsCalculateVisibilities.hpp>
+#include <lsDomain.hpp>
+#include <lsExpand.hpp>
+#include <lsMakeGeometry.hpp>
+#include <lsToMesh.hpp>
+#include <lsVTKWriter.hpp>
+
+namespace ls = viennals;
+
+int main() {
+
+  constexpr int D = 2;
+
+  double gridDelta = 0.4;
+
+  auto sphere1 = ls::SmartPointer<ls::Domain<double, D>>::New(
+      gridDelta); //, boundaryCons);
+
+  double origin[3] = {5., 0., 0.};
+  double radius = 7.3;
+
+  ls::MakeGeometry<double, D>(
+      sphere1, ls::SmartPointer<ls::Sphere<double, D>>::New(origin, radius))
+      .apply();
+
+  ls::Expand<double, D>(sphere1, 5).apply();
+  ls::CalculateVisibilities<double, D>(sphere1, ls::Vec3D<double>{0., -1.0, 0.})
+      .apply();
+
+  auto mesh = ls::SmartPointer<ls::Mesh<>>::New();
+  ls::ToMesh<double, D>(sphere1, mesh).apply();
+  ls::VTKWriter<double>(mesh, "visibility_test.vtp").apply();
+
+  return 0;
+}