From cc2620659059d2363a70c3ea14394310e20a0477 Mon Sep 17 00:00:00 2001 From: Tobias Reiter Date: Tue, 26 Nov 2024 22:54:58 +0100 Subject: [PATCH] Add visibility test and parallelize --- include/viennals/lsCalculateVisibilities.hpp | 137 +++++++++++-------- tests/Visibility/CMakeLists.txt | 7 + tests/Visibility/Visibility.cpp | 38 +++++ 3 files changed, 122 insertions(+), 60 deletions(-) create mode 100644 tests/Visibility/CMakeLists.txt create mode 100644 tests/Visibility/Visibility.cpp diff --git a/include/viennals/lsCalculateVisibilities.hpp b/include/viennals/lsCalculateVisibilities.hpp index f5d07e76..da46ff65 100644 --- a/include/viennals/lsCalculateVisibilities.hpp +++ b/include/viennals/lsCalculateVisibilities.hpp @@ -57,79 +57,98 @@ template class CalculateVisibilities { auto numDefinedPoints = domain.getNumberOfPoints(); std::vector visibilities(numDefinedPoints); - // std::vector> visitedCells; - - hrleSizeType id = 0; - hrleSparseIterator::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 startVector = + (p == 0) ? grid.getMinGridPoint() : domain.getSegmentation()[p - 1]; + + hrleVectorType endVector = + (p != static_cast(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 currentPos; - for (int i = 0; i < D; ++i) { - currentPos[i] = it.getStartIndices(i); - } + for (hrleSparseIterator::DomainType> it( + domain, startVector); + it.getStartIndices() < endVector; ++it) { - // Start tracing the ray - NumericType minLevelSetValue = it.getValue(); // Starting level set value - Vec3D rayPos = currentPos; - bool visibility = true; + if (!it.isDefined()) + continue; - while (1) { - // Update the ray position + // Starting position of the point + Vec3D 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 nearestCell; - for (int i = 0; i < D; ++i) { - nearestCell[i] = static_cast(rayPos[i]); - } + // Start tracing the ray + NumericType minLevelSetValue = + it.getValue(); // Starting level set value + Vec3D 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 nearestCell; + for (int i = 0; i < D; ++i) { + nearestCell[i] = static_cast(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::max(); - hrleSparseIterator::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::max(); + hrleSparseIterator::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 CalculateVisibilities { pointData.eraseScalarData(i); } pointData.insertNextScalarData(visibilities, visibilitiesLabel); - - assert(id == numDefinedPoints); } }; diff --git a/tests/Visibility/CMakeLists.txt b/tests/Visibility/CMakeLists.txt new file mode 100644 index 00000000..3d61b659 --- /dev/null +++ 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 $) diff --git a/tests/Visibility/Visibility.cpp b/tests/Visibility/Visibility.cpp new file mode 100644 index 00000000..7624aaf5 --- /dev/null +++ b/tests/Visibility/Visibility.cpp @@ -0,0 +1,38 @@ +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace ls = viennals; + +int main() { + + constexpr int D = 2; + + double gridDelta = 0.4; + + auto sphere1 = ls::SmartPointer>::New( + gridDelta); //, boundaryCons); + + double origin[3] = {5., 0., 0.}; + double radius = 7.3; + + ls::MakeGeometry( + sphere1, ls::SmartPointer>::New(origin, radius)) + .apply(); + + ls::Expand(sphere1, 5).apply(); + ls::CalculateVisibilities(sphere1, ls::Vec3D{0., -1.0, 0.}) + .apply(); + + auto mesh = ls::SmartPointer>::New(); + ls::ToMesh(sphere1, mesh).apply(); + ls::VTKWriter(mesh, "visibility_test.vtp").apply(); + + return 0; +}