ViennaLS
Loading...
Searching...
No Matches
lsCalculateVisibilities.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <hrleSparseIterator.hpp>
4#include <lsDomain.hpp>
5#include <vcVectorType.hpp>
6
7namespace viennals {
8using namespace viennacore;
9
10template <class T, int D> class CalculateVisibilities {
11 using hrleDomainType = typename Domain<T, D>::DomainType;
12
13 SmartPointer<Domain<T, D>> levelSet;
14 Vec3D<T> direction;
15 const T epsilon = static_cast<T>(1e-6);
16 const std::string visibilitiesLabel;
17
18public:
19 CalculateVisibilities(const SmartPointer<Domain<T, D>> &passedLevelSet,
20 const Vec3D<T> &passedDirection,
21 std::string label = "Visibilities")
22 : levelSet(passedLevelSet), direction(passedDirection),
23 visibilitiesLabel(std::move(label)) {}
24
25 void apply() {
26
27 auto &domain = levelSet->getDomain();
28 auto &grid = levelSet->getGrid();
29
30 // *** Determine extents of domain ***
31 Vec3D<T> minDefinedPoint;
32 Vec3D<T> maxDefinedPoint;
33 // Initialize with extreme values
34 for (int i = 0; i < D; ++i) {
35 minDefinedPoint[i] = std::numeric_limits<T>::max();
36 maxDefinedPoint[i] = std::numeric_limits<T>::lowest();
37 }
38 // Iterate through all defined points in the domain
39 for (viennahrle::SparseIterator<hrleDomainType> it(domain);
40 !it.isFinished(); it.next()) {
41 if (!it.isDefined())
42 continue; // Skip undefined points
43
44 // Get the coordinate of the current point
45 auto point = it.getStartIndices();
46 for (int i = 0; i < D; ++i) {
47 // Compare to update min and max defined points
48 T coord = point[i]; // * grid.getGridDelta();
49 minDefinedPoint[i] = std::min(minDefinedPoint[i], coord);
50 maxDefinedPoint[i] = std::max(maxDefinedPoint[i], coord);
51 }
52 }
53 //****************************
54
55 // Invert the vector
56 auto dir = Normalize(Inv(direction)) *
57 static_cast<T>(domain.getGrid().getGridDelta());
58
59 auto numDefinedPoints = domain.getNumberOfPoints();
60 std::vector<T> visibilities(numDefinedPoints);
61
62#pragma omp parallel num_threads(levelSet->getNumberOfSegments())
63 {
64 int p = 0;
65#ifdef _OPENMP
66 p = omp_get_thread_num();
67#endif
68
69 const viennahrle::Index<D> startVector =
70 (p == 0) ? grid.getMinGridPoint() : domain.getSegmentation()[p - 1];
71
72 const viennahrle::Index<D> endVector =
73 (p != static_cast<int>(domain.getNumberOfSegments() - 1))
74 ? domain.getSegmentation()[p]
75 : grid.incrementIndices(grid.getMaxGridPoint());
76
77 // Calculate the starting index for the visibility vector
78 viennahrle::SizeType id = 0;
79 for (int i = 0; i < p; ++i) {
80 id += domain.getDomainSegment(i).getNumberOfPoints();
81 }
82
83 for (viennahrle::SparseIterator<hrleDomainType> it(domain, startVector);
84 it.getStartIndices() < endVector; ++it) {
85
86 if (!it.isDefined())
87 continue;
88
89 // Starting position of the point
90 Vec3D<T> currentPos;
91 for (int i = 0; i < D; ++i) {
92 currentPos[i] = it.getStartIndices(i);
93 }
94
95 // Start tracing the ray
96 T minLevelSetValue = it.getValue(); // Starting level set value
97 Vec3D<T> rayPos = currentPos;
98 bool visibility = true;
99
100 while (true) {
101 // Update the ray position
102 for (int i = 0; i < D; ++i) {
103 rayPos[i] += dir[i];
104 }
105
106 // Determine the nearest grid cell (round to nearest index)
107 viennahrle::Index<D> nearestCell;
108 for (int i = 0; i < D; ++i) {
109 nearestCell[i] = static_cast<viennahrle::IndexType>(rayPos[i]);
110 }
111
112 // // Before adding a cell, check if it's already visited
113 // if (std::find(visitedCells.begin(), visitedCells.end(),
114 // nearestCell)
115 // == visitedCells.end()) {
116 // visitedCells.push_back(nearestCell);
117 // }
118
119 // Check if the nearest cell is within bounds
120 bool outOfBounds = false;
121 for (int i = 0; i < D; ++i) {
122 if (nearestCell[i] < minDefinedPoint[i] ||
123 nearestCell[i] > maxDefinedPoint[i]) {
124 outOfBounds = true;
125 break;
126 }
127 }
128
129 if (outOfBounds) {
130 break; // Ray is outside the grid
131 }
132
133 // Access the level set value at the nearest cell
134 T value =
135 viennahrle::SparseIterator<hrleDomainType>(domain, nearestCell)
136 .getValue();
137
138 // Update the minimum value encountered
139 if (value < minLevelSetValue) {
140 visibility = false;
141 break;
142 }
143 }
144
145 // Update visibility for this point
146 visibilities[id++] = visibility ? 1.0 : 0.0;
147 }
148 }
149
150 auto &pointData = levelSet->getPointData();
151 // delete if already exists
152 if (int i = pointData.getScalarDataIndex(visibilitiesLabel); i != -1) {
153 pointData.eraseScalarData(i);
154 }
155 pointData.insertNextScalarData(visibilities, visibilitiesLabel);
156 }
157};
158
159} // namespace viennals
CalculateVisibilities(const SmartPointer< Domain< T, D > > &passedLevelSet, const Vec3D< T > &passedDirection, std::string label="Visibilities")
Definition lsCalculateVisibilities.hpp:19
void apply()
Definition lsCalculateVisibilities.hpp:25
Class containing all information about the level set, including the dimensions of the domain,...
Definition lsDomain.hpp:27
viennahrle::Domain< T, D > DomainType
Definition lsDomain.hpp:32
Definition lsAdvect.hpp:36
constexpr int D
Definition pyWrap.cpp:70
double T
Definition pyWrap.cpp:68