ViennaLS
Loading...
Searching...
No Matches
lsMarkVoidPoints.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <hrleSparseStarIterator.hpp>
4
5#include <lsDomain.hpp>
6#include <lsGraph.hpp>
7
8namespace viennals {
9
10using namespace viennacore;
11
19enum struct VoidTopSurfaceEnum : unsigned {
24};
25
28template <class T, int D> class MarkVoidPoints {
29 using IndexType = std::size_t;
30
31 SmartPointer<Domain<T, D>> domain = nullptr;
32 bool reverseVoidDetection = false;
33 bool saveComponents = false;
34 bool detectLargestSurface = false;
35
36 // two points are connected if they have the same sign
37 static bool areConnected(const T &value1, const T &value2) {
38 return (value1 >= 0) == (value2 >= 0);
39 }
40
41 std::vector<IndexType> static mergeComponentCounts(
42 const std::vector<IndexType> &components,
43 const std::vector<IndexType> &pointsPerComponent) {
44 // find number of connected components after merge
45 // TODO: the last element in the components vector is very likely
46 // the highest index, so could just set it instead of checking all values
47 unsigned highestIndex = 0;
48 for (auto it : components) {
49 if (highestIndex < it) {
50 highestIndex = it;
51 }
52 }
53
54 // merge pointsPerComponent according to the connected components
55 std::vector<IndexType> pointsPerConnected;
56 pointsPerConnected.resize(highestIndex + 1, 0);
57 for (unsigned i = 0; i < components.size(); ++i) {
58 pointsPerConnected[components[i]] += pointsPerComponent[i];
59 }
60
61 return pointsPerConnected;
62 }
63
64 IndexType
65 calculateTopID(const std::vector<IndexType> &components,
66 const std::vector<IndexType> &pointsPerConnected) const {
67 // check which component has the most points
68 IndexType topId = 0;
69 // use first component, which contains more than 0 points
70 while (topId < pointsPerConnected.size() &&
71 pointsPerConnected[topId] == 0) {
72 ++topId;
73 }
74 for (unsigned i = topId + 1; i < pointsPerConnected.size(); ++i) {
75 if ((pointsPerConnected[topId] < pointsPerConnected[i]) !=
76 reverseVoidDetection) {
77 topId = i;
78 }
79 }
80
81 return topId;
82 }
83
84public:
85 static constexpr char voidPointLabel[] = "VoidPointMarkers";
86
87 MarkVoidPoints() = default;
88
89 MarkVoidPoints(SmartPointer<Domain<T, D>> passedlsDomain,
90 bool passedReverseVoidDetection = false)
91 : domain(passedlsDomain),
92 reverseVoidDetection(passedReverseVoidDetection) {}
93
94 void setLevelSet(SmartPointer<Domain<T, D>> passedlsDomain) {
95 domain = passedlsDomain;
96 }
97
102 void setReverseVoidDetection(bool passedReverseVoidDetection) {
103 reverseVoidDetection = passedReverseVoidDetection;
104 }
105
112 void setDetectLargestSurface(bool passedDetect) {
113 detectLargestSurface = passedDetect;
114 }
115
119 switch (topSurface) {
121 reverseVoidDetection = true;
122 detectLargestSurface = false;
123 break;
125 reverseVoidDetection = false;
126 detectLargestSurface = false;
127 break;
129 reverseVoidDetection = false;
130 detectLargestSurface = true;
131 break;
133 reverseVoidDetection = true;
134 detectLargestSurface = true;
135 break;
136
137 default:
138 Logger::getInstance()
139 .addWarning("MarkVoidPoints: Invalid VoidTopSurfaceEnum set. "
140 "Using default values.")
141 .print();
142 reverseVoidDetection = false;
143 detectLargestSurface = false;
144 break;
145 }
146 }
147
151 void setSaveComponentIds(bool scid) { saveComponents = scid; }
152
153 void apply() {
154 lsInternal::Graph graph;
155
156 std::vector<std::vector<std::vector<IndexType>>> componentList(
157 domain->getNumberOfSegments());
158 for (unsigned segmentId = 0; segmentId < domain->getNumberOfSegments();
159 ++segmentId) {
160 componentList[segmentId].resize(D + 1);
161 for (int dim = -1; dim < D; ++dim) {
162 componentList[segmentId][dim + 1].resize(
163 domain->getDomain().getNumberOfRuns(segmentId, dim), -1);
164 }
165 }
166
167 std::size_t numberOfComponents = 0;
168 std::vector<IndexType> pointsPerComponent;
169
170 // cycle through and set up the graph to get connectivity information
171 for (viennahrle::ConstSparseStarIterator<typename Domain<T, D>::DomainType,
172 1>
173 neighborIt(domain->getDomain());
174 !neighborIt.isFinished(); neighborIt.next()) {
175 auto &center = neighborIt.getCenter();
176
177 // component id of the current run
178 auto &currentComponentId =
179 componentList[center.getSegmentId()][center.getLevel()]
180 [center.getRunTypePosition()];
181
182 // -1 means it is not set yet
183 if (currentComponentId == -1) {
184 for (int k = 0; k < 2 * D; ++k) {
185 auto &neighbor = neighborIt.getNeighbor(k);
186 const auto neighborComponentId =
187 componentList[neighbor.getSegmentId()][neighbor.getLevel()]
188 [neighbor.getRunTypePosition()];
189 // if neighbor is already defined
190 // set current component to neighbor component
191 // if they are connected
192 if (neighborComponentId != -1) {
193 if (areConnected(center.getValue(), neighbor.getValue())) {
194 currentComponentId = neighborComponentId;
195 if (center.getValue() >= 0.)
196 ++pointsPerComponent[currentComponentId];
197 break;
198 }
199 }
200 }
201 }
202
203 // it is still not set, so add new vertex
204 if (currentComponentId == -1) {
205 currentComponentId = numberOfComponents;
206 pointsPerComponent.push_back((center.getValue() > 0.) ? 1 : 0);
207 graph.insertNextVertex();
208 ++numberOfComponents;
209 }
210
211 // check if edge can be set
212 for (int k = 0; k < 2 * D; ++k) {
213 auto &neighbor = neighborIt.getNeighbor(k);
214 auto &neighborComponentId =
215 componentList[neighbor.getSegmentId()][neighbor.getLevel()]
216 [neighbor.getRunTypePosition()];
217 if (areConnected(center.getValue(), neighbor.getValue())) {
218 // if neighbor is already set
219 if (neighborComponentId != -1) {
220 // if neighbor is part of different component
221 if (currentComponentId != neighborComponentId) {
222 graph.insertNextEdge(currentComponentId, neighborComponentId);
223 }
224 } else {
225 neighborComponentId = currentComponentId;
226 if (neighbor.getValue() >= 0.)
227 ++pointsPerComponent[neighborComponentId];
228 }
229 }
230 }
231 }
232
233 auto components = graph.getConnectedComponents();
234
235 int topComponent =
236 (reverseVoidDetection) ? components[0] : components.back();
237
238 // identify which layer to keep
239 {
240 auto pointsPerConnected =
241 mergeComponentCounts(components, pointsPerComponent);
242
243 // find largest connected surface
244 if (detectLargestSurface) {
245 topComponent = calculateTopID(components, pointsPerConnected);
246 } else { // if component does not contain points, take the next one
247 while (pointsPerConnected[topComponent] == 0) {
248 if (reverseVoidDetection)
249 ++topComponent;
250 else
251 --topComponent;
252 }
253 }
254 }
255
256 std::vector<T> voidPointMarkers;
257 voidPointMarkers.resize(domain->getNumberOfPoints());
258
259 std::vector<T> componentMarkers;
260 if (saveComponents)
261 componentMarkers.resize(domain->getNumberOfPoints());
262
263 // cycle through again to set correct voidPointMarkers
264 for (viennahrle::ConstSparseStarIterator<typename Domain<T, D>::DomainType,
265 1>
266 neighborIt(domain->getDomain());
267 !neighborIt.isFinished(); neighborIt.next()) {
268 auto center = neighborIt.getCenter();
269
270 if (!center.isDefined())
271 continue;
272
273 // if it is positive, just check if it is part of the top component
274 if (center.getValue() >= 0) {
275 const int &oldComponentId = componentList[center.getSegmentId()][0]
276 [center.getRunTypePosition()];
277 voidPointMarkers[center.getPointId()] =
278 (components[oldComponentId] != topComponent);
279 } else {
280 // if it is negative, check all neighbours (with different sign),
281 // because they are part of of a positive component, which might
282 // be the top
283 unsigned k;
284 for (k = 0; k < 2 * D; ++k) {
285 auto &neighbor = neighborIt.getNeighbor(k);
286 if (std::signbit(neighbor.getValue()) ==
287 std::signbit(center.getValue()))
288 continue;
289 const int &oldneighborComponentId =
290 componentList[neighbor.getSegmentId()][neighbor.getLevel()]
291 [neighbor.getRunTypePosition()];
292 if (components[oldneighborComponentId] == topComponent) {
293 break;
294 }
295 }
296 voidPointMarkers[center.getPointId()] = (k == 2 * D);
297 }
298
299 if (saveComponents) {
300 const int &oldComponentId = componentList[center.getSegmentId()][0]
301 [center.getRunTypePosition()];
302 componentMarkers[center.getPointId()] = components[oldComponentId];
303 }
304 }
305
306 auto &pointData = domain->getPointData();
307 auto voidMarkersPointer = pointData.getScalarData(voidPointLabel, true);
308 // if vector data does not exist
309 if (voidMarkersPointer == nullptr) {
310 pointData.insertNextScalarData(voidPointMarkers, voidPointLabel);
311 } else {
312 *voidMarkersPointer = std::move(voidPointMarkers);
313 }
314
315 if (saveComponents) {
316 auto componentMarkersPointer =
317 pointData.getScalarData("ConnectedComponentId", true);
318 // if vector data does not exist
319 if (componentMarkersPointer == nullptr) {
320 pointData.insertNextScalarData(componentMarkers,
321 "ConnectedComponentId");
322 } else {
323 *componentMarkersPointer = std::move(componentMarkers);
324 }
325 }
326 }
327};
328
329} // namespace viennals
Definition lsGraph.hpp:13
std::size_t insertNextVertex()
add new vertex
Definition lsGraph.hpp:47
std::vector< IndexType > getConnectedComponents()
Definition lsGraph.hpp:68
void insertNextEdge(std::size_t vertex1, std::size_t vertex2)
add connection of vertex1 to vertex2
Definition lsGraph.hpp:54
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
void setLevelSet(SmartPointer< Domain< T, D > > passedlsDomain)
Definition lsMarkVoidPoints.hpp:94
void setReverseVoidDetection(bool passedReverseVoidDetection)
Set whether the "top" level set should be the most positive(default) connected chain of level set val...
Definition lsMarkVoidPoints.hpp:102
void apply()
Definition lsMarkVoidPoints.hpp:153
void setSaveComponentIds(bool scid)
Set whether the connected component IDs used to generate the void points should be saved....
Definition lsMarkVoidPoints.hpp:151
void setDetectLargestSurface(bool passedDetect)
Set whether the number of points of one connected surface should be used to detect void points....
Definition lsMarkVoidPoints.hpp:112
static constexpr char voidPointLabel[]
Definition lsMarkVoidPoints.hpp:85
void setVoidTopSurface(VoidTopSurfaceEnum topSurface)
Set which connected component to use as the top surface and mark all other components as void points.
Definition lsMarkVoidPoints.hpp:118
MarkVoidPoints(SmartPointer< Domain< T, D > > passedlsDomain, bool passedReverseVoidDetection=false)
Definition lsMarkVoidPoints.hpp:89
Definition lsAdvect.hpp:36
VoidTopSurfaceEnum
Enumeration describing which connected component to use as top surface during void point detection....
Definition lsMarkVoidPoints.hpp:19
@ LEX_HIGHEST
Definition lsMarkVoidPoints.hpp:21
@ SMALLEST
Definition lsMarkVoidPoints.hpp:23
@ LEX_LOWEST
Definition lsMarkVoidPoints.hpp:20
@ LARGEST
Definition lsMarkVoidPoints.hpp:22
constexpr int D
Definition pyWrap.cpp:70
double T
Definition pyWrap.cpp:68