From 301ba079d1184eea453d2e10eb47e5570f5adc95 Mon Sep 17 00:00:00 2001 From: Tobias Reiter Date: Thu, 21 Nov 2024 15:43:28 +0100 Subject: [PATCH 01/12] Add CalculateVisibilities function --- include/viennals/lsCalculateVisibilities.hpp | 97 ++++++++++++++++++++ 1 file changed, 97 insertions(+) create mode 100644 include/viennals/lsCalculateVisibilities.hpp diff --git a/include/viennals/lsCalculateVisibilities.hpp b/include/viennals/lsCalculateVisibilities.hpp new file mode 100644 index 00000000..d5b5d60a --- /dev/null +++ b/include/viennals/lsCalculateVisibilities.hpp @@ -0,0 +1,97 @@ +#pragma once + +#include + +namespace viennals { +using namespace viennacore; + +template class CalculateVisibilities { + SmartPointer> levelSet; + int openBoundaryDirection; + bool isOpenBoundaryNegative; + +public: + static constexpr char visibilitiesLabel[] = "Visibilities"; + + CalculateVisibilities( + const SmartPointer> &passedLevelSet, + int passedOpenBoundaryDirection, bool passedIsOpenBoundaryNegative) + : levelSet(passedLevelSet), + openBoundaryDirection(passedOpenBoundaryDirection), + isOpenBoundaryNegative(passedIsOpenBoundaryDirection) {} + + void apply() { + + auto &domain = levelSet->getDomain(); + auto &grid = levelSet->getGrid(); + const NumericType max = std::numeric_limits::max(); + + auto numDefinedPoints = domain.getNumberOfPoints(); + std::vector visibilities(numDefinedPoints); + + std::vector oldIndices( + D - 1 - openBoundaryDirection, + std::numeric_limits::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); + } + + std::vector minValues(size, max); + + hrleSizeType id = 0; + + hrleSparseIterator::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); + } + + unsigned int pos_begin = 0; + unsigned int pos_end = 0; + + 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; + } + + for (unsigned int i = pos_begin; i <= pos_end; ++i) + minValues.at(i) = std::min(minValues.at(i), it.getValue()); + + if (isOpenBoundaryNegative) { + it.next(); + } else { + it.previous(); + } + } + + levelSet->getPointData().insertNextScalarData(visibilities, + visibilitiesLabel); + + assert(id == numDefinedPoints); + } +}; + +} // namespace viennals \ No newline at end of file From fab9760c20751a514f05d0693979ec6ee02f674c Mon Sep 17 00:00:00 2001 From: Tobias Reiter Date: Thu, 21 Nov 2024 15:45:49 +0100 Subject: [PATCH 02/12] Fix typo --- include/viennals/lsCalculateVisibilities.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/viennals/lsCalculateVisibilities.hpp b/include/viennals/lsCalculateVisibilities.hpp index d5b5d60a..81bb0a5c 100644 --- a/include/viennals/lsCalculateVisibilities.hpp +++ b/include/viennals/lsCalculateVisibilities.hpp @@ -18,7 +18,7 @@ template class CalculateVisibilities { int passedOpenBoundaryDirection, bool passedIsOpenBoundaryNegative) : levelSet(passedLevelSet), openBoundaryDirection(passedOpenBoundaryDirection), - isOpenBoundaryNegative(passedIsOpenBoundaryDirection) {} + isOpenBoundaryNegative(passedIsOpenBoundaryNegative) {} void apply() { From 76b3dcfbb6b9341c7990448d2065301161c75686 Mon Sep 17 00:00:00 2001 From: Tobias Reiter Date: Thu, 21 Nov 2024 15:46:47 +0100 Subject: [PATCH 03/12] Fix iterator template parameter --- include/viennals/lsCalculateVisibilities.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/viennals/lsCalculateVisibilities.hpp b/include/viennals/lsCalculateVisibilities.hpp index 81bb0a5c..a05d4301 100644 --- a/include/viennals/lsCalculateVisibilities.hpp +++ b/include/viennals/lsCalculateVisibilities.hpp @@ -45,7 +45,7 @@ template class CalculateVisibilities { hrleSizeType id = 0; - hrleSparseIterator::DomainType> it( + hrleSparseIterator::DomainType> it( domain, !isOpenBoundaryNegative); while (!it.isFinished()) { From 476b2743ed6bfc6247487a6b4ff64d1d48da5cd2 Mon Sep 17 00:00:00 2001 From: Tobias Reiter Date: Thu, 21 Nov 2024 16:14:29 +0100 Subject: [PATCH 04/12] Default params --- include/viennals/lsCalculateVisibilities.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/viennals/lsCalculateVisibilities.hpp b/include/viennals/lsCalculateVisibilities.hpp index a05d4301..8d6bfe09 100644 --- a/include/viennals/lsCalculateVisibilities.hpp +++ b/include/viennals/lsCalculateVisibilities.hpp @@ -15,7 +15,8 @@ template class CalculateVisibilities { CalculateVisibilities( const SmartPointer> &passedLevelSet, - int passedOpenBoundaryDirection, bool passedIsOpenBoundaryNegative) + int passedOpenBoundaryDirection = D - 1, + bool passedIsOpenBoundaryNegative = false) : levelSet(passedLevelSet), openBoundaryDirection(passedOpenBoundaryDirection), isOpenBoundaryNegative(passedIsOpenBoundaryNegative) {} From 8200717d453c7a8f558ab4d922aa70ed0e344587 Mon Sep 17 00:00:00 2001 From: Tobias Reiter Date: Fri, 22 Nov 2024 09:44:55 +0100 Subject: [PATCH 05/12] Delete visibility scalar data if it already exists --- include/viennals/lsCalculateVisibilities.hpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/include/viennals/lsCalculateVisibilities.hpp b/include/viennals/lsCalculateVisibilities.hpp index 8d6bfe09..a1d685b3 100644 --- a/include/viennals/lsCalculateVisibilities.hpp +++ b/include/viennals/lsCalculateVisibilities.hpp @@ -88,8 +88,12 @@ template class CalculateVisibilities { } } - levelSet->getPointData().insertNextScalarData(visibilities, - visibilitiesLabel); + auto &pointData = levelSet->getPointData(); + // delete if already exists + if (int i = pointData.getScalarDataIndex(visibilitiesLabel); i != -1) { + pointData.eraseScalarData(i); + } + pointData.insertNextScalarData(visibilities, visibilitiesLabel); assert(id == numDefinedPoints); } From 0854e3851c677d32b9bcfb55b88e6345d53abbfd Mon Sep 17 00:00:00 2001 From: = <=> Date: Tue, 26 Nov 2024 09:30:03 +0100 Subject: [PATCH 06/12] Visibilities calculated in any direction --- include/viennals/lsCalculateVisibilities.hpp | 161 ++++++++++++------- 1 file changed, 99 insertions(+), 62 deletions(-) diff --git a/include/viennals/lsCalculateVisibilities.hpp b/include/viennals/lsCalculateVisibilities.hpp index a1d685b3..dc007542 100644 --- a/include/viennals/lsCalculateVisibilities.hpp +++ b/include/viennals/lsCalculateVisibilities.hpp @@ -1,99 +1,136 @@ #pragma once #include +#include namespace viennals { using namespace viennacore; template class CalculateVisibilities { SmartPointer> levelSet; - int openBoundaryDirection; - bool isOpenBoundaryNegative; + Vec3D direction; + const NumericType epsilon = static_cast(1e-6); public: static constexpr char visibilitiesLabel[] = "Visibilities"; CalculateVisibilities( const SmartPointer> &passedLevelSet, - int passedOpenBoundaryDirection = D - 1, - bool passedIsOpenBoundaryNegative = false) - : levelSet(passedLevelSet), - openBoundaryDirection(passedOpenBoundaryDirection), - isOpenBoundaryNegative(passedIsOpenBoundaryNegative) {} + const Vec3D passedDirection) + : levelSet(passedLevelSet), direction(passedDirection) {} void apply() { auto &domain = levelSet->getDomain(); auto &grid = levelSet->getGrid(); - const NumericType max = std::numeric_limits::max(); - auto numDefinedPoints = domain.getNumberOfPoints(); - std::vector visibilities(numDefinedPoints); - - std::vector oldIndices( - D - 1 - openBoundaryDirection, - std::numeric_limits::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 minDefinedPoint; + Vec3D maxDefinedPoint; + // Initialize with extreme values + for (int i = 0; i < D; ++i) { + minDefinedPoint[i] = std::numeric_limits::max(); + maxDefinedPoint[i] = std::numeric_limits::lowest(); } - - std::vector minValues(size, max); - - hrleSizeType id = 0; - - hrleSparseIterator::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::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 visibilities(numDefinedPoints); - for (unsigned int i = pos_begin; i <= pos_end; ++i) - minValues.at(i) = std::min(minValues.at(i), it.getValue()); + // std::vector> visitedCells; + + hrleSizeType id = 0; + hrleSparseIterator::DomainType> it(domain); + + while (!it.isFinished()) { + if (!it.isDefined()) { + it.next(); + continue; + } + + // Starting position of the point + Vec3D 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 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 nearestCell; + for (int i = 0; i < D; ++i) { + nearestCell[i] = static_cast(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::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; + } + } - 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); } From 2fc81d7a4ec44c78e05215438e30dd7452dbeffe Mon Sep 17 00:00:00 2001 From: Tobias Reiter Date: Tue, 26 Nov 2024 21:56:19 +0100 Subject: [PATCH 07/12] Change visibility label --- include/viennals/lsCalculateVisibilities.hpp | 155 ++++++++++--------- 1 file changed, 81 insertions(+), 74 deletions(-) diff --git a/include/viennals/lsCalculateVisibilities.hpp b/include/viennals/lsCalculateVisibilities.hpp index dc007542..f5d07e76 100644 --- a/include/viennals/lsCalculateVisibilities.hpp +++ b/include/viennals/lsCalculateVisibilities.hpp @@ -10,14 +10,15 @@ template class CalculateVisibilities { SmartPointer> levelSet; Vec3D direction; const NumericType epsilon = static_cast(1e-6); + const std::string visibilitiesLabel; public: - static constexpr char visibilitiesLabel[] = "Visibilities"; - CalculateVisibilities( const SmartPointer> &passedLevelSet, - const Vec3D passedDirection) - : levelSet(passedLevelSet), direction(passedDirection) {} + const Vec3D passedDirection, + const std::string label = "Visibilities") + : levelSet(passedLevelSet), direction(passedDirection), + visibilitiesLabel(std::move(label)) {} void apply() { @@ -29,12 +30,15 @@ template class CalculateVisibilities { Vec3D maxDefinedPoint; // Initialize with extreme values for (int i = 0; i < D; ++i) { - minDefinedPoint[i] = std::numeric_limits::max(); - maxDefinedPoint[i] = std::numeric_limits::lowest(); + minDefinedPoint[i] = std::numeric_limits::max(); + maxDefinedPoint[i] = std::numeric_limits::lowest(); } // Iterate through all defined points in the domain - for (hrleSparseIterator::DomainType> it(domain); !it.isFinished(); it.next()) { - if (!it.isDefined()) continue; // Skip undefined points + for (hrleSparseIterator::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(); @@ -53,84 +57,87 @@ template class CalculateVisibilities { auto numDefinedPoints = domain.getNumberOfPoints(); std::vector visibilities(numDefinedPoints); - // std::vector> visitedCells; - - hrleSizeType id = 0; - hrleSparseIterator::DomainType> it(domain); + // std::vector> visitedCells; - while (!it.isFinished()) { - if (!it.isDefined()) { - it.next(); - continue; - } + hrleSizeType id = 0; + hrleSparseIterator::DomainType> it(domain); + + while (!it.isFinished()) { + if (!it.isDefined()) { + it.next(); + continue; + } + + // Starting position of the point + Vec3D 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 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 nearestCell; + for (int i = 0; i < D; ++i) { + nearestCell[i] = static_cast(rayPos[i]); + } - // Starting position of the point - Vec3D currentPos; - for (int i = 0; i < D; ++i) { - currentPos[i] = it.getStartIndices(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; } + } - // Start tracing the ray - NumericType minLevelSetValue = it.getValue(); // Starting level set value - Vec3D 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 nearestCell; - for (int i = 0; i < D; ++i) { - nearestCell[i] = static_cast(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::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; - } + if (outOfBounds) { + break; // Ray is outside the grid } - // Update visibility for this point - visibilities[id++] = visibility ? 1.0 : 0.0; - it.next(); + // 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 visibility for this point + visibilities[id++] = visibility ? 1.0 : 0.0; + it.next(); + } + auto &pointData = levelSet->getPointData(); // delete if already exists - if (int i = pointData.getScalarDataIndex("Visibilities"); i != -1) { + if (int i = pointData.getScalarDataIndex(visibilitiesLabel); i != -1) { pointData.eraseScalarData(i); } - pointData.insertNextScalarData(visibilities, "Visibilities"); + pointData.insertNextScalarData(visibilities, visibilitiesLabel); assert(id == numDefinedPoints); } From cc2620659059d2363a70c3ea14394310e20a0477 Mon Sep 17 00:00:00 2001 From: Tobias Reiter Date: Tue, 26 Nov 2024 22:54:58 +0100 Subject: [PATCH 08/12] 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; +} From c447ffeeea6419201cba3184fce6257587f1bf30 Mon Sep 17 00:00:00 2001 From: = <=> Date: Thu, 28 Nov 2024 09:33:17 +0100 Subject: [PATCH 09/12] Fixed visibility calculation --- include/viennals/lsCalculateVisibilities.hpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/include/viennals/lsCalculateVisibilities.hpp b/include/viennals/lsCalculateVisibilities.hpp index da46ff65..a1872f12 100644 --- a/include/viennals/lsCalculateVisibilities.hpp +++ b/include/viennals/lsCalculateVisibilities.hpp @@ -135,9 +135,7 @@ template class CalculateVisibilities { hrleSparseIterator::DomainType> neighborIt(domain); neighborIt.goToIndices(nearestCell); - if (neighborIt.isDefined()) { - neighborValue = neighborIt.getValue(); - } + neighborValue = neighborIt.getValue(); // Update the minimum value encountered if (neighborValue < minLevelSetValue) { From 49384f57d0fb83b1780ddeb976b9a2c575da1f35 Mon Sep 17 00:00:00 2001 From: Tobias Reiter Date: Thu, 28 Nov 2024 15:03:06 +0100 Subject: [PATCH 10/12] Scale direction vector by grid delta --- include/viennals/lsCalculateVisibilities.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/viennals/lsCalculateVisibilities.hpp b/include/viennals/lsCalculateVisibilities.hpp index a1872f12..c1a7ee1e 100644 --- a/include/viennals/lsCalculateVisibilities.hpp +++ b/include/viennals/lsCalculateVisibilities.hpp @@ -52,7 +52,7 @@ template class CalculateVisibilities { //**************************** // Invert the vector - auto dir = Normalize(Inv(direction)); + auto dir = Normalize(Inv(direction)) * domiain.getGrid().getGridDelta(); auto numDefinedPoints = domain.getNumberOfPoints(); std::vector visibilities(numDefinedPoints); From 5a94f29861be1f1857d44e8602c014a350ba934f Mon Sep 17 00:00:00 2001 From: Tobias Reiter Date: Thu, 28 Nov 2024 15:32:31 +0100 Subject: [PATCH 11/12] Fix typo --- include/viennals/lsCalculateVisibilities.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/viennals/lsCalculateVisibilities.hpp b/include/viennals/lsCalculateVisibilities.hpp index c1a7ee1e..105688b1 100644 --- a/include/viennals/lsCalculateVisibilities.hpp +++ b/include/viennals/lsCalculateVisibilities.hpp @@ -52,7 +52,7 @@ template class CalculateVisibilities { //**************************** // Invert the vector - auto dir = Normalize(Inv(direction)) * domiain.getGrid().getGridDelta(); + auto dir = Normalize(Inv(direction)) * domain.getGrid().getGridDelta(); auto numDefinedPoints = domain.getNumberOfPoints(); std::vector visibilities(numDefinedPoints); From 916579b35d23c2bc0b3da70441c3937b6e2b3d31 Mon Sep 17 00:00:00 2001 From: Tobias Reiter Date: Thu, 28 Nov 2024 16:08:35 +0100 Subject: [PATCH 12/12] Cast grid delta to correct type --- include/viennals/lsCalculateVisibilities.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/viennals/lsCalculateVisibilities.hpp b/include/viennals/lsCalculateVisibilities.hpp index 105688b1..a1fea5c5 100644 --- a/include/viennals/lsCalculateVisibilities.hpp +++ b/include/viennals/lsCalculateVisibilities.hpp @@ -52,7 +52,8 @@ template class CalculateVisibilities { //**************************** // Invert the vector - auto dir = Normalize(Inv(direction)) * domain.getGrid().getGridDelta(); + auto dir = Normalize(Inv(direction)) * + static_cast(domain.getGrid().getGridDelta()); auto numDefinedPoints = domain.getNumberOfPoints(); std::vector visibilities(numDefinedPoints);