TTK
Loading...
Searching...
No Matches
LegacyTopologicalSimplification.h
Go to the documentation of this file.
1
91
92#pragma once
93
94// base code includes
95#include <Debug.h>
96#include <Triangulation.h>
97
98#include <cmath>
99#include <set>
100#include <tuple>
101#include <type_traits>
102
103namespace ttk {
104
105 class SweepCmp {
106 bool isIncreasingOrder_{};
107
108 public:
109 SweepCmp() = default;
110
111 SweepCmp(bool isIncreasingOrder) : isIncreasingOrder_{isIncreasingOrder} {
112 }
113
114 inline void setIsIncreasingOrder(bool isIncreasingOrder) {
115 isIncreasingOrder_ = isIncreasingOrder;
116 }
117
118 template <typename dataType>
119 bool
120 operator()(const std::tuple<dataType, SimplexId, SimplexId> &v0,
121 const std::tuple<dataType, SimplexId, SimplexId> &v1) const {
122 if(isIncreasingOrder_) {
123 return std::get<1>(v0) < std::get<1>(v1);
124 } else {
125 return std::get<1>(v0) > std::get<1>(v1);
126 }
127 }
128 };
129
130 class LegacyTopologicalSimplification : virtual public Debug {
131 public:
133
134 template <typename triangulationType>
135 int getCriticalType(SimplexId vertexId,
136 const SimplexId *const offsets,
137 const triangulationType &triangulation) const;
138
139 template <typename triangulationType>
140 int getCriticalPoints(const SimplexId *const offsets,
141 std::vector<SimplexId> &minList,
142 std::vector<SimplexId> &maxList,
143 const triangulationType &triangulation) const;
144
145 template <typename triangulationType>
146 int getCriticalPoints(const SimplexId *const offsets,
147 std::vector<SimplexId> &minList,
148 std::vector<SimplexId> &maxList,
149 std::vector<bool> &blackList,
150 const triangulationType &triangulation) const;
151
152 template <typename dataType>
153 int addPerturbation(dataType *const scalars,
154 SimplexId *const offsets) const;
155
164 template <typename dataType, typename triangulationType>
165 int execute(const dataType *const inputScalars,
166 dataType *const outputScalars,
167 const SimplexId *const identifiers,
168 const SimplexId *const inputOffsets,
169 SimplexId *const offsets,
170 const SimplexId constraintNumber,
171 const triangulationType &triangulation) const;
172
174 if(triangulation) {
175 vertexNumber_ = triangulation->getNumberOfVertices();
176 triangulation->preconditionVertexNeighbors();
177 }
178 return 0;
179 }
180
181 inline void setConsiderIdentifierAsBlackList(bool onOff) {
183 }
184
185 inline void setAddPerturbation(bool onOff) {
186 addPerturbation_ = onOff;
187 }
188
189 protected:
192 bool addPerturbation_{false};
193 };
194} // namespace ttk
195
196// if the package is a pure template typename, uncomment the following line
197// #include <LegacyTopologicalSimplification.cpp>
198
199template <typename triangulationType>
201 SimplexId vertex,
202 const SimplexId *const offsets,
203 const triangulationType &triangulation) const {
204
205 bool isMinima{true};
206 bool isMaxima{true};
207 SimplexId neighborNumber = triangulation.getVertexNeighborNumber(vertex);
208 for(SimplexId i = 0; i < neighborNumber; ++i) {
209 SimplexId neighbor{-1};
210 triangulation.getVertexNeighbor(vertex, i, neighbor);
211
212 if(offsets[neighbor] < offsets[vertex])
213 isMinima = false;
214 if(offsets[neighbor] > offsets[vertex])
215 isMaxima = false;
216 if(!isMinima and !isMaxima) {
217 return 0;
218 }
219 }
220
221 if(isMinima)
222 return -1;
223 if(isMaxima)
224 return 1;
225
226 return 0;
227}
228
229template <typename triangulationType>
231 const SimplexId *const offsets,
232 std::vector<SimplexId> &minima,
233 std::vector<SimplexId> &maxima,
234 const triangulationType &triangulation) const {
235
236 std::vector<int> type(vertexNumber_, 0);
237
238#ifdef TTK_ENABLE_OPENMP
239#pragma omp parallel for
240#endif
241 for(SimplexId k = 0; k < vertexNumber_; ++k)
242 type[k] = getCriticalType(k, offsets, triangulation);
243
244 for(SimplexId k = 0; k < vertexNumber_; ++k) {
245 if(type[k] < 0)
246 minima.push_back(k);
247 else if(type[k] > 0)
248 maxima.push_back(k);
249 }
250 return 0;
251}
252
253template <typename triangulationType>
255 const SimplexId *const offsets,
256 std::vector<SimplexId> &minima,
257 std::vector<SimplexId> &maxima,
258 std::vector<bool> &extrema,
259 const triangulationType &triangulation) const {
260 std::vector<int> type(vertexNumber_);
261#ifdef TTK_ENABLE_OPENMP
262#pragma omp parallel for num_threads(threadNumber_)
263#endif
264 for(SimplexId k = 0; k < vertexNumber_; ++k) {
265 if(considerIdentifierAsBlackList_ xor extrema[k]) {
266 type[k] = getCriticalType(k, offsets, triangulation);
267 }
268 }
269
270 for(SimplexId k = 0; k < vertexNumber_; ++k) {
271 if(type[k] < 0)
272 minima.push_back(k);
273 else if(type[k] > 0)
274 maxima.push_back(k);
275 }
276 return 0;
277}
278
279template <typename dataType>
281 dataType *const scalars, SimplexId *const offsets) const {
282 dataType epsilon{};
283
284 if(std::is_same<dataType, double>::value)
285 epsilon = Geometry::powIntTen<dataType>(1 - DBL_DIG);
286 else if(std::is_same<dataType, float>::value)
287 epsilon = Geometry::powIntTen<dataType>(1 - FLT_DIG);
288 else
289 return -1;
290
291 std::vector<std::tuple<dataType, SimplexId, SimplexId>> perturbation(
292 vertexNumber_);
293 for(SimplexId i = 0; i < vertexNumber_; ++i) {
294 std::get<0>(perturbation[i]) = scalars[i];
295 std::get<1>(perturbation[i]) = offsets[i];
296 std::get<2>(perturbation[i]) = i;
297 }
298
299 SweepCmp cmp(true);
300 sort(perturbation.begin(), perturbation.end(), cmp);
301
302 for(SimplexId i = 0; i < vertexNumber_; ++i) {
303 if(i) {
304 if(std::get<0>(perturbation[i]) <= std::get<0>(perturbation[i - 1]))
305 std::get<0>(perturbation[i])
306 = std::get<0>(perturbation[i - 1]) + epsilon;
307 }
308 scalars[std::get<2>(perturbation[i])] = std::get<0>(perturbation[i]);
309 }
310
311 return 0;
312}
313
314template <typename dataType, typename triangulationType>
316 const dataType *const inputScalars,
317 dataType *const outputScalars,
318 const SimplexId *const identifiers,
319 const SimplexId *const inputOffsets,
320 SimplexId *const offsets,
321 const SimplexId constraintNumber,
322 const triangulationType &triangulation) const {
323
324 Timer t;
325
326 // pre-processing
327#ifdef TTK_ENABLE_OPENMP
328#pragma omp parallel for num_threads(threadNumber_)
329#endif
330 for(SimplexId k = 0; k < vertexNumber_; ++k) {
331 outputScalars[k] = inputScalars[k];
332 if(std::isnan((double)outputScalars[k]))
333 outputScalars[k] = 0;
334
335 offsets[k] = inputOffsets[k];
336 }
337
338 // get the user extremum list
339 std::vector<bool> extrema(vertexNumber_, false);
340 for(SimplexId k = 0; k < constraintNumber; ++k) {
341 const SimplexId identifierId = identifiers[k];
342
343#ifndef TTK_ENABLE_KAMIKAZE
344 if(identifierId >= 0 and identifierId < vertexNumber_)
345#endif
346 extrema[identifierId] = true;
347 }
348
349 std::vector<SimplexId> authorizedMinima;
350 std::vector<SimplexId> authorizedMaxima;
351 std::vector<bool> authorizedExtrema(vertexNumber_, false);
352
353 getCriticalPoints(
354 offsets, authorizedMinima, authorizedMaxima, extrema, triangulation);
355
356 this->printMsg("Maintaining " + std::to_string(constraintNumber)
357 + " constraints (" + std::to_string(authorizedMinima.size())
358 + " minima and " + std::to_string(authorizedMaxima.size())
359 + " maxima)",
361
362 // declare the tuple-comparison functor
363 SweepCmp cmp;
364
365 // processing
366 for(SimplexId i = 0; i < vertexNumber_; ++i) {
367
368 this->printMsg("Starting simplifying iteration #" + std::to_string(i),
370
371 for(int j = 0; j < 2; ++j) {
372
373 bool isIncreasingOrder = !j;
374
375 if(isIncreasingOrder && authorizedMinima.empty()) {
376 continue;
377 }
378 if(!isIncreasingOrder && authorizedMaxima.empty()) {
379 continue;
380 }
381
382 cmp.setIsIncreasingOrder(isIncreasingOrder);
383 std::set<std::tuple<dataType, SimplexId, SimplexId>, decltype(cmp)>
384 sweepFront(cmp);
385 std::vector<bool> visitedVertices(vertexNumber_, false);
386 std::vector<SimplexId> adjustmentSequence(vertexNumber_);
387
388 // add the seeds
389 if(isIncreasingOrder) {
390 for(SimplexId k : authorizedMinima) {
391 authorizedExtrema[k] = true;
392 sweepFront.emplace(outputScalars[k], offsets[k], k);
393 visitedVertices[k] = true;
394 }
395 } else {
396 for(SimplexId k : authorizedMaxima) {
397 authorizedExtrema[k] = true;
398 sweepFront.emplace(outputScalars[k], offsets[k], k);
399 visitedVertices[k] = true;
400 }
401 }
402
403 // growth by neighborhood of the seeds
404 SimplexId adjustmentPos = 0;
405 do {
406 auto front = sweepFront.begin();
407 if(front == sweepFront.end())
408 return -1;
409
410 SimplexId vertexId = std::get<2>(*front);
411 sweepFront.erase(front);
412
413 SimplexId neighborNumber
414 = triangulation.getVertexNeighborNumber(vertexId);
415 for(SimplexId k = 0; k < neighborNumber; ++k) {
416 SimplexId neighbor{-1};
417 triangulation.getVertexNeighbor(vertexId, k, neighbor);
418 if(!visitedVertices[neighbor]) {
419 sweepFront.emplace(
420 outputScalars[neighbor], offsets[neighbor], neighbor);
421 visitedVertices[neighbor] = true;
422 }
423 }
424 adjustmentSequence[adjustmentPos] = vertexId;
425 ++adjustmentPos;
426 } while(!sweepFront.empty());
427
428 // save offsets and rearrange outputScalars
429 SimplexId offset = (isIncreasingOrder ? 0 : vertexNumber_);
430
431 for(SimplexId k = 0; k < vertexNumber_; ++k) {
432
433 if(isIncreasingOrder) {
434 if(k
435 && outputScalars[adjustmentSequence[k]]
436 <= outputScalars[adjustmentSequence[k - 1]])
437 outputScalars[adjustmentSequence[k]]
438 = outputScalars[adjustmentSequence[k - 1]];
439 ++offset;
440 } else {
441 if(k
442 && outputScalars[adjustmentSequence[k]]
443 >= outputScalars[adjustmentSequence[k - 1]])
444 outputScalars[adjustmentSequence[k]]
445 = outputScalars[adjustmentSequence[k - 1]];
446 --offset;
447 }
448 offsets[adjustmentSequence[k]] = offset;
449 }
450 }
451
452 // test convergence
453 bool needForMoreIterations{false};
454 std::vector<SimplexId> minima;
455 std::vector<SimplexId> maxima;
456 getCriticalPoints(offsets, minima, maxima, triangulation);
457
458 if(maxima.size() > authorizedMaxima.size())
459 needForMoreIterations = true;
460 if(minima.size() > authorizedMinima.size())
461 needForMoreIterations = true;
462
463 this->printMsg(
464 std::vector<std::vector<std::string>>{
465 {"#Minima", std::to_string(minima.size())},
466 {"#Maxima", std::to_string(maxima.size())},
467 },
469
470 if(!needForMoreIterations) {
471 for(SimplexId k : minima) {
472 if(!authorizedExtrema[k]) {
473 needForMoreIterations = true;
474 break;
475 }
476 }
477 }
478 if(!needForMoreIterations) {
479 for(SimplexId k : maxima) {
480 if(!authorizedExtrema[k]) {
481 needForMoreIterations = true;
482 break;
483 }
484 }
485 }
486
487 // optional adding of perturbation
488 if(addPerturbation_)
489 addPerturbation<dataType>(outputScalars, offsets);
490
491 if(!needForMoreIterations)
492 break;
493 }
494
495 this->printMsg(
496 "Simplified scalar field", 1.0, t.getElapsedTime(), this->threadNumber_);
497
498 return 0;
499}
AbstractTriangulation is an interface class that defines an interface for efficient traversal methods...
virtual SimplexId getNumberOfVertices() const
Minimalist debugging class.
Definition: Debug.h:88
TTK processing package for the topological simplification of scalar data.
int preconditionTriangulation(AbstractTriangulation *triangulation)
int execute(const dataType *const inputScalars, dataType *const outputScalars, const SimplexId *const identifiers, const SimplexId *const inputOffsets, SimplexId *const offsets, const SimplexId constraintNumber, const triangulationType &triangulation) const
int addPerturbation(dataType *const scalars, SimplexId *const offsets) const
int getCriticalType(SimplexId vertexId, const SimplexId *const offsets, const triangulationType &triangulation) const
int getCriticalPoints(const SimplexId *const offsets, std::vector< SimplexId > &minList, std::vector< SimplexId > &maxList, const triangulationType &triangulation) const
SweepCmp(bool isIncreasingOrder)
void setIsIncreasingOrder(bool isIncreasingOrder)
SweepCmp()=default
bool operator()(const std::tuple< dataType, SimplexId, SimplexId > &v0, const std::tuple< dataType, SimplexId, SimplexId > &v1) const
double getElapsedTime()
Definition: Timer.h:15
The Topology ToolKit.
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22
printMsg(debug::output::GREEN+"                           "+debug::output::ENDCOLOR+debug::output::GREEN+"▒"+debug::output::ENDCOLOR+debug::output::GREEN+"▒▒▒▒▒▒▒▒▒▒▒▒▒░"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)