lsCalculateVisibilities.hpp

Go to the documentation of this file.

#pragma once
 
#include <hrleSparseIterator.hpp>
#include <lsDomain.hpp>
#include <vcVectorType.hpp>
 
namespace viennals {
using namespace viennacore;
 
template <class T, int D> class CalculateVisibilities {
  using hrleDomainType = typename Domain<T, D>::DomainType;
 
  SmartPointer<Domain<T, D>> levelSet;
  Vec3D<T> direction;
  const T epsilon = static_cast<T>(1e-6);
  const std::string visibilitiesLabel;
 
public:
  CalculateVisibilities(const SmartPointer<Domain<T, D>> &passedLevelSet,
                        const Vec3D<T> &passedDirection,
                        std::string label = "Visibilities")
      : levelSet(passedLevelSet), direction(passedDirection),
        visibilitiesLabel(std::move(label)) {}
 
  void apply() {
    auto &domain = levelSet->getDomain();
    auto &grid = levelSet->getGrid();
 
    // Invert the vector
    auto dir = Normalize(Inv(direction));
    std::vector<T> visibilities(domain.getNumberOfPoints(), static_cast<T>(-1));
 
    // Check if the direction vector has a component in the simulation
    // dimensions
    T dirNorm = 0;
    for (int i = 0; i < D; ++i) {
      dirNorm += dir[i] * dir[i];
    }
 
    if (dirNorm < epsilon) {
      std::fill(visibilities.begin(), visibilities.end(), static_cast<T>(1.0));
    } else {
      // *** Determine extents of domain ***
      Vec3D<T> minDefinedPoint{};
      Vec3D<T> maxDefinedPoint{};
      // Initialize with extreme values
      for (int i = 0; i < D; ++i) {
        minDefinedPoint[i] = std::numeric_limits<T>::max();
        maxDefinedPoint[i] = std::numeric_limits<T>::lowest();
      }
      // Iterate through all defined points in the domain
      for (viennahrle::SparseIterator<hrleDomainType> 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
          T coord = point[i]; // * grid.getGridDelta();
          minDefinedPoint[i] = std::min(minDefinedPoint[i], coord);
          maxDefinedPoint[i] = std::max(maxDefinedPoint[i], coord);
        }
      }
      //****************************
 
#pragma omp parallel num_threads(levelSet->getNumberOfSegments())
      {
        int p = 0;
#ifdef _OPENMP
        p = omp_get_thread_num();
#endif
 
        const viennahrle::Index<D> startVector =
            (p == 0) ? grid.getMinGridPoint() : domain.getSegmentation()[p - 1];
 
        const viennahrle::Index<D> endVector =
            (p != static_cast<int>(domain.getNumberOfSegments() - 1))
                ? domain.getSegmentation()[p]
                : grid.incrementIndices(grid.getMaxGridPoint());
 
        for (viennahrle::SparseIterator<hrleDomainType> it(domain, startVector);
             it.getStartIndices() < endVector; ++it) {
 
          if (!it.isDefined())
            continue;
 
          // Starting position of the point
          Vec3D<T> currentPos{};
          for (int i = 0; i < D; ++i) {
            currentPos[i] = it.getStartIndices(i);
          }
 
          // Start tracing the ray
          T minLevelSetValue = it.getValue(); // Starting level set value
          Vec3D<T> rayPos = currentPos;
 
          // Step once to skip immediate neighbor
          for (int i = 0; i < D; ++i)
            rayPos[i] += dir[i];
 
          bool visibility = true;
 
          while (true) {
            // Update the ray position
            for (int i = 0; i < D; ++i)
              rayPos[i] += dir[i];
 
            // Determine the nearest grid cell (round to nearest index)
            viennahrle::Index<D> nearestCell;
            for (int i = 0; i < D; ++i)
              nearestCell[i] = static_cast<viennahrle::IndexType>(rayPos[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;
              }
            }
 
            if (outOfBounds)
              break; // Ray is outside the grid
 
            // Access the level set value at the nearest cell
            T value =
                viennahrle::SparseIterator<hrleDomainType>(domain, nearestCell)
                    .getValue();
 
            // Update the minimum value encountered
            if (value < minLevelSetValue) {
              visibility = false;
              break;
            }
          }
 
          // Update visibility for this point
          visibilities[it.getPointId()] = visibility ? 1.0 : 0.0;
        }
      }
    }
 
    int unassignedCount = 0;
    for (size_t i = 0; i < visibilities.size(); ++i) {
      if (visibilities[i] < 0) {
        std::cerr << "Unassigned visibility at point ID: " << i << std::endl;
        ++unassignedCount;
      }
    }
    if (unassignedCount > 0) {
      std::cerr << "[Error] Total unassigned points: " << unassignedCount
                << std::endl;
    }
 
    auto &pointData = levelSet->getPointData();
    // delete if already exists
    if (int i = pointData.getScalarDataIndex(visibilitiesLabel); i != -1) {
      pointData.eraseScalarData(i);
    }
    pointData.insertNextScalarData(visibilities, visibilitiesLabel);
  }
};
 
} // namespace viennals