From 0854e3851c677d32b9bcfb55b88e6345d53abbfd Mon Sep 17 00:00:00 2001 From: = <=> Date: Tue, 26 Nov 2024 09:30:03 +0100 Subject: [PATCH] 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); }