TTK
Loading...
Searching...
No Matches
ApproximateTopology.h
Go to the documentation of this file.
1
19
20#pragma once
21
22// base code includes
23#include <DynamicTree.h>
24#include <MultiresTopology.h>
25#include <MultiresTriangulation.h>
26
27#include <limits>
28#include <numeric>
29#include <tuple>
30
31namespace ttk {
32
33 using triplet = std::tuple<ttk::SimplexId, ttk::SimplexId, ttk::SimplexId>;
34 using polarity = unsigned char;
35
36 class ApproximateTopology : public MultiresTopology {
37
38 public:
39 ApproximateTopology() {
40 this->setDebugMsgPrefix("ApproximateTopology");
41 }
42 void setEpsilon(double data) {
43 epsilon_ = data;
44 }
45 void setDelta(double data) {
46 delta_ = data;
47 }
48 template <typename scalarType>
49 int computeApproximatePD(std::vector<PersistencePair> &CTDiagram,
50 const scalarType *scalars,
51 scalarType *const fakeScalars,
52 SimplexId *const outputOffsets,
53 int *const outputMonotonyOffsets);
54
55 template <typename scalarType>
56 int executeApproximateTopology(const scalarType *scalars,
57 scalarType *fakeScalars,
58 SimplexId *outputOffsets,
59 int *outputMonotonyOffsets);
60
61 template <typename scalarType, typename offsetType>
62 void
63 initGlobalPolarity(std::vector<polarity> &isNew,
64 std::vector<std::vector<std::pair<polarity, polarity>>>
65 &vertexLinkPolarity,
66 std::vector<polarity> &toProcess,
67 const scalarType *fakeScalars,
68 const offsetType *const offsets,
69 const int *const monotonyOffsets) const;
70
71 template <typename scalarType, typename offsetType>
72 void updatePropagation(
73 std::vector<polarity> &toPropagateMin,
74 std::vector<polarity> &toPropagateMax,
75 std::vector<std::vector<SimplexId>> &vertexRepresentativesMin,
76 std::vector<std::vector<SimplexId>> &vertexRepresentativesMax,
77 std::vector<std::vector<SimplexId>> &saddleCCMin,
78 std::vector<std::vector<SimplexId>> &saddleCCMax,
79 std::vector<Lock> &vertLockMin,
80 std::vector<Lock> &vertLockMax,
81 std::vector<polarity> &isUpdatedMin,
82 std::vector<polarity> &isUpdatedMax,
83 const scalarType *fakeScalars,
84 const offsetType *const offsetField,
85 const int *const monotonyOffsets);
86
87 template <typename scalarType, typename offsetType>
88 void buildVertexLinkPolarityApproximate(
89 const SimplexId vertexId,
90 std::vector<std::pair<polarity, polarity>> &vlp,
91 const scalarType *fakeScalars,
92 const offsetType *const offsetField,
93 const int *const monotonyOffsets) const;
94
95 template <typename ScalarType, typename OffsetType>
96 void sortTripletsApproximate(std::vector<triplet> &triplets,
97 const ScalarType *const scalars,
98 const ScalarType *const fakeScalars,
99 const OffsetType *const offsets,
100 const int *const monotonyOffsets,
101 const bool splitTree) const;
102
103 template <typename scalarType, typename offsetType>
104 void getCriticalTypeApproximate(
105 const SimplexId &vertexId,
106 std::vector<std::pair<polarity, polarity>> &vlp,
107 uint8_t &vertexLink,
108 DynamicTree &link,
109 VLBoundaryType &vlbt,
110 const scalarType *fakeScalars,
111 const offsetType *const offsets,
112 const int *const monotonyOffsets) const;
113
114 template <typename ScalarType, typename offsetType>
115 void computeCriticalPoints(
116 std::vector<std::vector<std::pair<polarity, polarity>>>
117 &vertexLinkPolarity,
118 std::vector<polarity> &toPropagateMin,
119 std::vector<polarity> &toPropagateMax,
120 std::vector<polarity> &toProcess,
121 std::vector<DynamicTree> &link,
122 std::vector<uint8_t> &vertexLink,
123 VLBoundaryType &vertexLinkByBoundaryType,
124 std::vector<std::vector<SimplexId>> &saddleCCMin,
125 std::vector<std::vector<SimplexId>> &saddleCCMax,
126 ScalarType *fakeScalars,
127 const offsetType *const offsets,
128 int *monotonyOffsets);
129
130 template <typename scalarType, typename offsetType>
131 ttk::SimplexId propagateFromSaddles(
132 const SimplexId vertexId,
133 std::vector<Lock> &vertLock,
134 std::vector<polarity> &toPropagate,
135 std::vector<std::vector<SimplexId>> &vertexRepresentatives,
136 std::vector<std::vector<SimplexId>> &saddleCC,
137 std::vector<polarity> &isUpdated,
138 std::vector<SimplexId> &globalExtremum,
139 const bool splitTree,
140 const scalarType *fakeScalars,
141 const offsetType *const offsetField,
142 const int *const monotonyOffsets) const;
143
144 template <typename ScalarType, typename offsetType>
145 void computePersistencePairsFromSaddles(
146 std::vector<PersistencePair> &CTDiagram,
147 const ScalarType *const fakeScalars,
148 const offsetType *const offsets,
149 const int *const monotonyOffsets,
150 std::vector<std::vector<SimplexId>> &vertexRepresentativesMin,
151 std::vector<std::vector<SimplexId>> &vertexRepresentativesMax,
152 const std::vector<polarity> &toPropagateMin,
153 const std::vector<polarity> &toPropagateMax) const;
154
155 template <typename scalarType, typename offsetType>
156 void sortVertices(const SimplexId vertexNumber,
157 std::vector<SimplexId> &sortedVertices,
158 SimplexId *vertsOrder,
159 const scalarType *const fakeScalars,
160 const offsetType *const offsetField,
161 const int *const monotonyOffsets);
162
163 template <typename scalarType>
164 int sortPersistenceDiagramApproximate(std::vector<PersistencePair> &diagram,
165 scalarType *fakeScalars,
166 const SimplexId *const offsets,
167 int *monotonyOffsets) const;
168
169 protected:
170 // maximum link size in 3D
171 static const size_t nLink_ = 27;
172 using VLBoundaryType
173 = std::array<std::vector<std::pair<SimplexId, SimplexId>>, nLink_>;
174
175 void
176 initCriticalPoints(std::vector<polarity> &isNew,
177 std::vector<std::vector<std::pair<polarity, polarity>>>
178 &vertexLinkPolarity,
179 std::vector<polarity> &toProcess,
180 std::vector<DynamicTree> &link,
181 std::vector<uint8_t> &vertexLink,
182 VLBoundaryType &vertexLinkByBoundaryType,
183 std::vector<char> &vertexTypes,
184 const SimplexId *const offsets) const;
185
186 void initSaddleSeeds(std::vector<polarity> &isNew,
187 std::vector<std::vector<std::pair<polarity, polarity>>>
188 &vertexLinkPolarity,
189 std::vector<polarity> &toPropagateMin,
190 std::vector<polarity> &toPropagateMax,
191 std::vector<polarity> &toProcess,
192 std::vector<DynamicTree> &link,
193 std::vector<uint8_t> &vertexLink,
194 VLBoundaryType &vertexLinkByBoundaryType,
195 std::vector<std::vector<SimplexId>> &saddleCCMin,
196 std::vector<std::vector<SimplexId>> &saddleCCMax,
197 const SimplexId *const offsets) const;
198
199 template <typename scalarType, typename offsetType>
200 void updatePropagation(
201 std::vector<polarity> &toPropagateMin,
202 std::vector<polarity> &toPropagateMax,
203 std::vector<std::vector<SimplexId>> &vertexRepresentativesMin,
204 std::vector<std::vector<SimplexId>> &vertexRepresentativesMax,
205 std::vector<std::vector<SimplexId>> &saddleCCMin,
206 std::vector<std::vector<SimplexId>> &saddleCCMax,
207 std::vector<Lock> &vertLockMin,
208 std::vector<Lock> &vertLockMax,
209 std::vector<polarity> &isUpdatedMin,
210 std::vector<polarity> &isUpdatedMax,
211 const scalarType *fakeScalars,
212 const offsetType *const offsets,
213 const int *const monotonyOffsets) const;
214
215 template <typename scalarType, typename offsetType>
216 void
217 buildVertexLinkPolarity(const SimplexId vertexId,
218 std::vector<std::pair<polarity, polarity>> &vlp,
219 const scalarType *fakeScalars,
220 const offsetType *const offsets,
221 const int *const monotonyOffsets) const;
222
223 template <typename ScalarType, typename OffsetType>
224 void sortTriplets(std::vector<triplet> &triplets,
225 const ScalarType *const fakeScalars,
226 const OffsetType *const offsets,
227 const int *const monotonyOffsets,
228 const bool splitTree) const;
229
230 template <typename scalarType, typename offsetType>
231 void tripletsToPersistencePairs(
232 std::vector<PersistencePair> &pairs,
233 std::vector<std::vector<SimplexId>> &vertexRepresentatives,
234 std::vector<triplet> &triplets,
235 const scalarType *const fakeScalars,
236 const offsetType *const offsets,
237 const int *const monotonyOffsets,
238 const bool splitTree) const;
239
240 template <typename scalarType, typename offsetType>
241 int getMonotonyChangeByOldPointCPApproximate(
242 const SimplexId vertexId,
243 double eps,
244 const std::vector<polarity> &isNew,
245 std::vector<polarity> &toProcess,
246 std::vector<polarity> &toReprocess,
247 std::vector<std::pair<polarity, polarity>> &vlp,
248 scalarType *fakeScalars,
249 const offsetType *const offsets,
250 int *monotonyOffsets) const;
251
252 template <typename ScalarType, typename offsetType>
253 int updateGlobalPolarity(
254 double eps,
255 std::vector<polarity> &isNew,
256 std::vector<std::vector<std::pair<polarity, polarity>>>
257 &vertexLinkPolarity,
258 std::vector<polarity> &toProcess,
259 std::vector<polarity> &toReprocess,
260 ScalarType *fakeScalars,
261 const offsetType *const offsets,
262 int *monotonyOffsets) const;
263
264 ttk::SimplexId propagateFromSaddles(
265 const SimplexId vertexId,
266 std::vector<Lock> &vertLock,
267 std::vector<polarity> &toPropagate,
268 std::vector<std::vector<SimplexId>> &vertexRepresentatives,
269 std::vector<std::vector<SimplexId>> &saddleCC,
270 std::vector<polarity> &isUpdated,
271 std::vector<SimplexId> &globalExtremum,
272 const SimplexId *const offsets,
273 const bool splitTree) const;
274
275 void computePersistencePairsFromSaddles(
276 std::vector<PersistencePair> &CTDiagram,
277 const SimplexId *const offsets,
278 std::vector<std::vector<SimplexId>> &vertexRepresentativesMin,
279 std::vector<std::vector<SimplexId>> &vertexRepresentativesMax,
280 const std::vector<polarity> &toPropagateMin,
281 const std::vector<polarity> &toPropagateMax) const;
282
283 template <typename scalarType, typename offsetType>
284 bool printPolarity(std::vector<polarity> &isNew,
285 const SimplexId vertexId,
286 std::vector<std::vector<std::pair<polarity, polarity>>>
287 &vertexLinkPolarity,
288 const scalarType *scalars,
289 const scalarType *fakeScalars,
290 const offsetType *const offsets,
291 const int *const monotonyOffsets,
292 bool verbose = false);
293
294 double epsilon_{};
295 double delta_{};
296 };
297} // namespace ttk
298
299template <typename scalarType>
300int ttk::ApproximateTopology::executeApproximateTopology(
301 const scalarType *ttkNotUsed(scalars),
302 scalarType *fakeScalars,
303 SimplexId *outputOffsets,
304 int *outputMonotonyOffsets) {
305
306 Timer timer;
307
308 SimplexId *const vertsOrder = outputOffsets;
309
310 decimationLevel_ = startingDecimationLevel_;
311 multiresTriangulation_.setTriangulation(triangulation_);
312 const SimplexId vertexNumber = multiresTriangulation_.getVertexNumber();
313
314 int *monotonyOffsets = outputMonotonyOffsets;
315
316#ifdef TTK_ENABLE_KAMIKAZE
317 if(vertexNumber == 0) {
318 this->printErr("No points in triangulation");
319 return 1;
320 }
321#endif // TTK_ENABLE_KAMIKAZE
322
323 double tm_allocation = timer.getElapsedTime();
324
325 const auto dim = multiresTriangulation_.getDimensionality();
326 const size_t maxNeigh = dim == 3 ? 14 : (dim == 2 ? 6 : 0);
327
328 std::vector<std::vector<SimplexId>> saddleCCMin(vertexNumber),
329 saddleCCMax(vertexNumber);
330 std::vector<std::vector<SimplexId>> vertexRepresentativesMin(vertexNumber),
331 vertexRepresentativesMax(vertexNumber);
332
333 std::vector<std::vector<std::pair<polarity, polarity>>> vertexLinkPolarity(
334 vertexNumber);
335
336 std::vector<polarity> isNew(vertexNumber, 255);
337 std::vector<polarity> toPropagateMin(vertexNumber, 0),
338 toPropagateMax(vertexNumber, 0);
339 std::vector<polarity> isUpToDateMin(vertexNumber, 0),
340 isUpToDateMax(vertexNumber, 0);
341
342 // index in vertexLinkByBoundaryType
343 std::vector<uint8_t> vertexLink(vertexNumber);
344 VLBoundaryType vertexLinkByBoundaryType{};
345 std::vector<DynamicTree> link(vertexNumber);
346 std::vector<polarity> toProcess(vertexNumber, 0), toReprocess{};
347
348 std::vector<SimplexId> offsetsVec(vertexNumber);
349 std::iota(offsetsVec.begin(), offsetsVec.end(), 0);
350 const SimplexId *const offsets = offsetsVec.data();
351
352 if(this->startingDecimationLevel_ > this->stoppingDecimationLevel_) {
353 // only needed for progressive computation
354 toReprocess.resize(vertexNumber, 0);
355 }
356
357 std::vector<Lock> vertLockMin(vertexNumber), vertLockMax(vertexNumber);
358
359 if(preallocateMemory_) {
360 double tm_prealloc = timer.getElapsedTime();
361 printMsg("Pre-allocating data structures", 0, 0, threadNumber_,
362 ttk::debug::LineMode::REPLACE);
363 for(SimplexId i = 0; i < vertexNumber; ++i) {
364 // for(int d = 0; d < vertexLinkPolarity.size(); d++) {
365 // vertexLinkPolarity[d][i].reserve(maxNeigh);
366 // }
367 vertexLinkPolarity[i].reserve(maxNeigh);
368 link[i].alloc(maxNeigh);
369 }
370 printMsg("Pre-allocating data structures", 1,
371 timer.getElapsedTime() - tm_prealloc, threadNumber_);
372 }
373
374 tm_allocation = timer.getElapsedTime() - tm_allocation;
375 printMsg("Total memory allocation", 1, tm_allocation, threadNumber_);
376
377 // computation of implicit link
378 std::vector<SimplexId> boundReps{};
379 multiresTriangulation_.findBoundaryRepresentatives(boundReps);
380
381#ifdef TTK_ENABLE_OPENMP
382#pragma omp parallel for num_threads(threadNumber_)
383#endif
384 for(size_t i = 0; i < boundReps.size(); i++) {
385 if(boundReps[i] != -1) {
386 buildVertexLinkByBoundary(boundReps[i], vertexLinkByBoundaryType);
387 }
388 }
389
390 // if(debugLevel_ > 4) {
391 // std::cout << "boundary representatives : ";
392 // for(auto bb : boundReps) {
393 // std::cout << ", " << bb;
394 // }
395 // std::cout << std::endl;
396 // }
397
398 multiresTriangulation_.setDecimationLevel(decimationLevel_);
399
400 initGlobalPolarity(isNew, vertexLinkPolarity, toProcess, fakeScalars, offsets,
401 monotonyOffsets);
402
403 double delta = epsilon_ * delta_;
404
405 while(decimationLevel_ > stoppingDecimationLevel_) {
406 Timer tmIter{};
407 decimationLevel_--;
408 multiresTriangulation_.setDecimationLevel(decimationLevel_);
409 // double progress = (double)(startingDecimationLevel_ - decimationLevel_)
410 // / (startingDecimationLevel_);
411 // this->printMsg("decimation level: " + std::to_string(decimationLevel_),
412 // progress, timer.getElapsedTime() - tm_allocation,
413 // threadNumber_);
414
415 int ret = updateGlobalPolarity(delta, isNew, vertexLinkPolarity, toProcess,
416 toReprocess, fakeScalars, offsets,
417 monotonyOffsets);
418 if(ret == -1) {
419 std::cout << "Found ERROR - aborting" << std::endl;
420 return -1;
421 }
422 } // end while
423
424 computeCriticalPoints(vertexLinkPolarity, toPropagateMin, toPropagateMax,
425 toProcess, link, vertexLink, vertexLinkByBoundaryType,
426 saddleCCMin, saddleCCMax, fakeScalars, offsets,
427 monotonyOffsets);
428
429 updatePropagation(toPropagateMin, toPropagateMax, vertexRepresentativesMin,
430 vertexRepresentativesMax, saddleCCMin, saddleCCMax,
431 vertLockMin, vertLockMax, isUpToDateMin, isUpToDateMax,
432 fakeScalars, offsets, monotonyOffsets);
433
434 computePersistencePairsFromSaddles(
435 CTDiagram_, fakeScalars, offsets, monotonyOffsets, vertexRepresentativesMin,
436 vertexRepresentativesMax, toPropagateMin, toPropagateMax);
437 // ADD GLOBAL MIN-MAX PAIR
438 CTDiagram_.emplace_back(this->globalMin_, this->globalMax_, -1);
439 // fakeScalars[this->globalMax_] - fakeScalars[this->globalMin_], -1);
440
441 // std::cout << "# epsilon " << epsilon_ << " in "
442 // << timer.getElapsedTime() - tm_allocation << " s." << std::endl;
443
444 this->printMsg("Complete", 1.0, timer.getElapsedTime() - tm_allocation,
445 this->threadNumber_);
446
447 // finally sort the diagram
448 // std::cout << "SORTING" << std::endl;
449 sortPersistenceDiagramApproximate(
450 CTDiagram_, fakeScalars, offsets, monotonyOffsets);
451 // std::cout << "Final epsilon " << epsilon_ << std::endl;
452
453 // sort vertices to generate correct output offset order
454 std::vector<SimplexId> sortedVertices{};
455 sortVertices(vertexNumber, sortedVertices, vertsOrder, fakeScalars, offsets,
456 monotonyOffsets);
457 return 0;
458}
459
460template <typename ScalarType, typename OffsetType>
461void ttk::ApproximateTopology::computePersistencePairsFromSaddles(
462 std::vector<PersistencePair> &CTDiagram,
463 const ScalarType *const fakeScalars,
464 const OffsetType *const offsets,
465 const int *const monotonyOffsets,
466 std::vector<std::vector<SimplexId>> &vertexRepresentativesMin,
467 std::vector<std::vector<SimplexId>> &vertexRepresentativesMax,
468 const std::vector<polarity> &toPropagateMin,
469 const std::vector<polarity> &toPropagateMax) const {
470
471 Timer timer{};
472 // CTDiagram.clear();
473 std::vector<triplet> tripletsMax{}, tripletsMin{};
474 const SimplexId nbDecVert = multiresTriangulation_.getDecimatedVertexNumber();
475
476 for(SimplexId localId = 0; localId < nbDecVert; localId++) {
477 SimplexId globalId = multiresTriangulation_.localToGlobalVertexId(localId);
478 if(toPropagateMin[globalId]) {
479 getTripletsFromSaddles(globalId, tripletsMin, vertexRepresentativesMin);
480 }
481 if(toPropagateMax[globalId]) {
482 getTripletsFromSaddles(globalId, tripletsMax, vertexRepresentativesMax);
483 }
484 }
485
486 sortTriplets(tripletsMax, fakeScalars, offsets, monotonyOffsets, true);
487 sortTriplets(tripletsMin, fakeScalars, offsets, monotonyOffsets, false);
488
489 const auto tm_sort = timer.getElapsedTime();
490
491 typename std::remove_reference<decltype(CTDiagram)>::type CTDiagramMin{},
492 CTDiagramMax{};
493
494#ifdef TTK_ENABLE_OPENMP
495#pragma omp parallel sections num_threads(threadNumber_)
496#endif // TTK_ENABLE_OPENMP
497 {
498#ifdef TTK_ENABLE_OPENMP
499#pragma omp section
500#endif // TTK_ENABLE_OPENMP
501 tripletsToPersistencePairs(CTDiagramMin, vertexRepresentativesMax,
502 tripletsMax, fakeScalars, offsets,
503 monotonyOffsets, true);
504#ifdef TTK_ENABLE_OPENMP
505#pragma omp section
506#endif // TTK_ENABLE_OPENMP
507 tripletsToPersistencePairs(CTDiagramMax, vertexRepresentativesMin,
508 tripletsMin, fakeScalars, offsets,
509 monotonyOffsets, false);
510 }
511 CTDiagram = std::move(CTDiagramMin);
512 CTDiagram.insert(CTDiagram.end(), CTDiagramMax.begin(), CTDiagramMax.end());
513
514 if(debugLevel_ > 3) {
515 std::cout << "PAIRS " << timer.getElapsedTime() - tm_sort << std::endl;
516 }
517}
518
519template <typename ScalarType, typename OffsetType>
520void ttk::ApproximateTopology::sortTriplets(std::vector<triplet> &triplets,
521 const ScalarType *const fakeScalars,
522 const OffsetType *const offsets,
523 const int *const monotonyOffsets,
524 const bool splitTree) const {
525 if(triplets.empty())
526 return;
527
528 const auto lt = [=](const SimplexId a, const SimplexId b) -> bool {
529 return ((fakeScalars[a] < fakeScalars[b])
530 || (fakeScalars[a] == fakeScalars[b]
531 && ((monotonyOffsets[a] < monotonyOffsets[b])
532 || (monotonyOffsets[a] == monotonyOffsets[b]
533 && offsets[a] < offsets[b]))));
534 };
535
536 // Sorting step
537 const auto cmp = [=](const triplet &t1, const triplet &t2) {
538 const SimplexId s1 = std::get<0>(t1);
539 const SimplexId s2 = std::get<0>(t2);
540 const SimplexId m1 = std::get<2>(t1);
541 const SimplexId m2 = std::get<2>(t2);
542 if(s1 != s2)
543 return lt(s1, s2) != splitTree;
544 else // s1 == s2
545 return lt(m1, m2) == splitTree;
546 };
547
548 TTK_PSORT(this->threadNumber_, triplets.begin(), triplets.end(), cmp);
549}
550
551template <typename scalarType, typename offsetType>
552void ttk::ApproximateTopology::tripletsToPersistencePairs(
553 std::vector<PersistencePair> &pairs,
554 std::vector<std::vector<SimplexId>> &vertexRepresentatives,
555 std::vector<triplet> &triplets,
556 const scalarType *const fakeScalars,
557 const offsetType *const offsets,
558 const int *const monotonyOffsets,
559
560 const bool splitTree) const {
561 Timer tm;
562 if(triplets.empty())
563 return;
564
565 // // accelerate getRep lookup?
566 // std::vector<SimplexId> firstRep(vertexRepresentatives.size());
567
568 // #ifdef TTK_ENABLE_OPENMP
569 // #pragma omp parallel for num_threads(threadNumber_)
570 // #endif // TTK_ENABLE_OPENMP
571 // for(size_t i = 0; i < firstRep.size(); ++i) {
572 // if(!vertexRepresentatives[i].empty()) {
573 // firstRep[i] = vertexRepresentatives[i][0];
574 // }
575 // }
576
577 const auto lt = [=](const SimplexId a, const SimplexId b) -> bool {
578 return ((fakeScalars[a] < fakeScalars[b])
579 || (fakeScalars[a] == fakeScalars[b]
580 && ((monotonyOffsets[a] < monotonyOffsets[b])
581 || (monotonyOffsets[a] == monotonyOffsets[b]
582 && offsets[a] < offsets[b]))));
583 };
584
585 const auto getRep = [&](SimplexId v) -> SimplexId {
586 auto r = vertexRepresentatives[v][0];
587 while(r != v) {
588 v = r;
589 r = vertexRepresentatives[v][0];
590 }
591 return r;
592 };
593
594 for(const auto &t : triplets) {
595 SimplexId r1 = getRep(std::get<1>(t));
596 SimplexId r2 = getRep(std::get<2>(t));
597 if(r1 != r2) {
598 SimplexId s = std::get<0>(t);
599
600 // Add pair
601 if(splitTree) {
602 // r1 = min(r1, r2), r2 = max(r1, r2)
603 if(lt(r2, r1)) {
604 std::swap(r1, r2);
605 }
606 // pair saddle-max: s -> min(r1, r2);
607 pairs.emplace_back(s, r1, 2);
608 // fakeScalars[r1] - fakeScalars[s], 2);
609
610 } else {
611 // r1 = max(r1, r2), r2 = min(r1, r2)
612 if(lt(r1, r2)) {
613 std::swap(r1, r2);
614 }
615 // pair min-saddle: max(r1, r2) -> s;
616 pairs.emplace_back(r1, s, 0);
617 // fakeScalars[s] - fakeScalars[r1], 0);
618 }
619
620 vertexRepresentatives[std::get<1>(t)][0] = r2;
621 vertexRepresentatives[r1][0] = r2;
622 // vertexRepresentatives[std::get<1>(t)][0] = r2;
623 // vertexRepresentatives[r1][0] = r2;
624 }
625 }
626
627 if(debugLevel_ > 3) {
628 std::string prefix = splitTree ? "[sad-max]" : "[min-sad]";
629 std::cout << prefix << " found all pairs in " << tm.getElapsedTime()
630 << " s." << std::endl;
631 }
632}
633
634template <typename scalarType, typename offsetType>
635void ttk::ApproximateTopology::sortVertices(
636 const SimplexId vertexNumber,
637 std::vector<SimplexId> &sortedVertices,
638 SimplexId *vertsOrder,
639 const scalarType *const fakeScalars,
640 const offsetType *const offsets,
641 const int *const monotonyOffsets) {
642
643 Timer tm;
644
645 sortedVertices.resize(vertexNumber);
646
647 // fill with numbers from 0 to vertexNumber - 1
648 std::iota(sortedVertices.begin(), sortedVertices.end(), 0);
649
650 // sort vertices in ascending order following scalarfield / offsets
651 TTK_PSORT(this->threadNumber_, sortedVertices.begin(), sortedVertices.end(),
652 [&](const SimplexId a, const SimplexId b) {
653 return ((fakeScalars[a] < fakeScalars[b])
654 || (fakeScalars[a] == fakeScalars[b]
655 && ((monotonyOffsets[a] < monotonyOffsets[b])
656 || (monotonyOffsets[a] == monotonyOffsets[b]
657 && offsets[a] < offsets[b]))));
658 });
659
660#ifdef TTK_ENABLE_OPENMP
661#pragma omp parallel for num_threads(threadNumber_)
662#endif // TTK_ENABLE_OPENMP
663 for(size_t i = 0; i < sortedVertices.size(); ++i) {
664 vertsOrder[sortedVertices[i]] = i;
665 }
666
667 // if(debugLevel_ > 2) {
668 // std::cout << "SORT " << tm.getElapsedTime() << std::endl;
669 // }
670}
671
672template <typename scalarType, typename offsetType>
673int ttk::ApproximateTopology::getMonotonyChangeByOldPointCPApproximate(
674 const SimplexId vertexId,
675 double eps,
676 const std::vector<polarity> &isNew,
677 std::vector<polarity> &toProcess,
678 std::vector<polarity> &toReprocess,
679 std::vector<std::pair<polarity, polarity>> &vlp,
680 scalarType *fakeScalars,
681 const offsetType *const offsets,
682 int *monotonyOffsets) const {
683
684 int hasMonotonyChanged = 0;
685 const SimplexId neighborNumber
686 = multiresTriangulation_.getVertexNeighborNumber(vertexId);
687 for(SimplexId i = 0; i < neighborNumber; i++) {
688 SimplexId neighborId = -1;
689 multiresTriangulation_.getVertexNeighbor(vertexId, i, neighborId);
690
691 // check for monotony changes
692 // const bool lowerStatic
693 // = (scalarField[neighborId] < fakeScalars[vertexId])
694 // || ((scalarField[neighborId] == fakeScalars[vertexId])
695 // && (offsets[neighborId] < offsets[vertexId]));
696 const bool lowerDynamic
697 = ((fakeScalars[neighborId] < fakeScalars[vertexId])
698 || (fakeScalars[neighborId] == fakeScalars[vertexId]
699 && ((monotonyOffsets[neighborId] < monotonyOffsets[vertexId])
700 || (monotonyOffsets[neighborId] == monotonyOffsets[vertexId]
701 && offsets[neighborId] < offsets[vertexId]))));
702
703 const polarity isUpperDynamic = lowerDynamic ? 0 : 255;
704 // const polarity isUpperStatic = lowerStatic ? 0 : 255;
705 const polarity isUpperOld = vlp[i].first;
706
707 if(isUpperDynamic != isUpperOld) { // change of monotony
708 SimplexId oldNeighbor = -1;
709 int oldDecimation = pow(2, decimationLevel_ + 1);
710 multiresTriangulation_.getVertexNeighborAtDecimation(
711 vertexId, i, oldNeighbor, oldDecimation);
712
713 double replacementValueDynamic
714 = (0.5 * (double)fakeScalars[oldNeighbor]
715 + .5 * (double)fakeScalars[vertexId]); // depends on epsilon
716 double deltaDynamic
717 = fabs((double)fakeScalars[neighborId] - replacementValueDynamic);
718
719 //=====================
720 SimplexId oldNeighNumber = 0;
721 SimplexId nnumber
722 = multiresTriangulation_.getVertexNeighborNumber(neighborId);
723 for(SimplexId iii = 0; iii < nnumber; iii++) {
724 SimplexId neighborId2 = -1;
725 multiresTriangulation_.getVertexNeighbor(neighborId, iii, neighborId2);
726 if(!isNew[neighborId2]) {
727 oldNeighNumber++;
728 }
729 }
730
731 if(deltaDynamic > eps or !isNew[neighborId] or oldNeighNumber > 2) {
732 hasMonotonyChanged = 1;
733
734 toReprocess[vertexId] = 255;
735 if(isNew[neighborId]) {
736 toProcess[neighborId] = 255;
737 } else {
738 toReprocess[neighborId] = 255;
739 }
740 const SimplexId neighborNumberNew
741 = multiresTriangulation_.getVertexNeighborNumber(neighborId);
742 for(SimplexId j = 0; j < neighborNumberNew; j++) {
743 SimplexId neighborIdNew = -1;
744 multiresTriangulation_.getVertexNeighbor(
745 neighborId, j, neighborIdNew);
746 if(isNew[neighborIdNew])
747 toProcess[neighborIdNew] = 255;
748 }
749 vlp[i].second = 255;
750 } else {
751 fakeScalars[neighborId] = replacementValueDynamic;
752
753 // corrects rounding error when we process an integer scalar field
754 if(fakeScalars[neighborId] == fakeScalars[oldNeighbor]) {
755 fakeScalars[neighborId] = fakeScalars[vertexId];
756 }
757
758 // change monotony offset
759 // The monotony must be preserved, meaning that
760 // if we had f(vertexId)>f(oldNeighbor)
761 // we should have f(vertexId)>f(neighborId)
762 // which is enforced when
763 // monotonyOffsets[vertexId]>monotonyOffsets[neighborId]
764 if(isUpperOld) { // we should enforce f(vertexId)<f(neighborId)
765 if(offsets[vertexId] > offsets[neighborId]) {
766 monotonyOffsets[neighborId]
767 = monotonyOffsets[vertexId] + pow(2, decimationLevel_);
768 if(monotonyOffsets[vertexId] == monotonyOffsets[oldNeighbor]
769 and fakeScalars[vertexId] == fakeScalars[oldNeighbor]) {
770 std::cout << "THIS IS AN ISSUE" << std::endl;
771 }
772 } else {
773 monotonyOffsets[neighborId] = monotonyOffsets[vertexId];
774 }
775 } else { // we should enforce f(vertexId)>f(neighborId)
776 if(offsets[vertexId] < offsets[neighborId]) {
777 monotonyOffsets[neighborId]
778 = monotonyOffsets[vertexId] - pow(2, decimationLevel_);
779 if(monotonyOffsets[vertexId] == monotonyOffsets[oldNeighbor]
780 and fakeScalars[vertexId] == fakeScalars[oldNeighbor]) {
781 std::cout << "THIS IS AN ISSUE" << std::endl;
782 }
783 } else {
784 monotonyOffsets[neighborId] = monotonyOffsets[vertexId];
785 }
786 }
787 }
788 } // end if change of monotony
789 } // end for neighbors
790 return hasMonotonyChanged;
791}
792
793template <typename scalarType, typename offsetType>
794ttk::SimplexId ttk::ApproximateTopology::propagateFromSaddles(
795 const SimplexId vertexId,
796 std::vector<Lock> &vertLock,
797 std::vector<polarity> &toPropagate,
798 std::vector<std::vector<SimplexId>> &vertexRepresentatives,
799 std::vector<std::vector<SimplexId>> &saddleCC,
800 std::vector<polarity> &isUpdated,
801 std::vector<SimplexId> &globalExtremum,
802 const bool splitTree,
803 const scalarType *fakeScalars,
804 const offsetType *const offsets,
805 const int *const monotonyOffsets) const {
806
807 auto &toProp = toPropagate[vertexId];
808 auto &reps = vertexRepresentatives[vertexId];
809 auto &updated = isUpdated[vertexId];
810
811 if(updated) {
812 return reps[0];
813 }
814
815 // const auto gt = [=](const SimplexId v1, const SimplexId v2) {
816 // return ((fakeScalars[v1] > fakeScalars[v2])
817 // || (fakeScalars[v1] == fakeScalars[v2]
818 // && offsets[v1] > offsets[v2]))
819 // == splitTree;
820 // };
821 const auto gt = [=](const SimplexId v1, const SimplexId v2) {
822 return ((fakeScalars[v1] > fakeScalars[v2])
823 || (fakeScalars[v1] == fakeScalars[v2]
824 && ((monotonyOffsets[v1] > monotonyOffsets[v2])
825 || (monotonyOffsets[v1] == monotonyOffsets[v2]
826 && offsets[v1] > offsets[v2]))))
827 == splitTree;
828 };
829
830 if(this->threadNumber_ > 1) {
831 vertLock[vertexId].lock();
832 }
833 if(saddleCC[vertexId].size()
834 and !toProp) { // tis a saddle point, should have to propagate on it
835 printErr("ERRRROR");
836 }
837 if(toProp) { // SADDLE POINT
838 if(debugLevel_ > 5)
839 printMsg("to saddle " + std::to_string(vertexId) + " "
840 + std::to_string(saddleCC[vertexId].size()));
841 // for(auto ccid : saddleCC[vertexId]) {
842 // SimplexId ccV = -1;
843 // multiresTriangulation_.getVertexNeighbor(vertexId, ccid, ccV);
844 // std::cout << " " << ccid << " (" << ccV << ") , ";
845 // }
846 // std::cout << std::endl;
847 const auto &CC = saddleCC[vertexId];
848 reps.clear();
849 reps.reserve(CC.size());
850 for(size_t r = 0; r < CC.size(); r++) {
851 SimplexId neighborId = -1;
852 SimplexId localId = CC[r];
853 multiresTriangulation_.getVertexNeighbor(vertexId, localId, neighborId);
854 // printMsg(" CC " + std::to_string(CC[r]) + " "
855 // + std::to_string(neighborId) + " ("
856 // + std::to_string(fakeScalars[neighborId]) + " , "
857 // + std::to_string(offsets[neighborId]) + ")");
858 SimplexId ret = propagateFromSaddles(
859 neighborId, vertLock, toPropagate, vertexRepresentatives, saddleCC,
860 isUpdated, globalExtremum, splitTree, fakeScalars, offsets,
861 monotonyOffsets);
862 reps.emplace_back(ret);
863 }
864
865 if(reps.size() > 1) {
866 // sort & remove duplicate elements
867 std::sort(reps.begin(), reps.end(), gt);
868 const auto last = std::unique(reps.begin(), reps.end());
869 reps.erase(last, reps.end());
870 }
871
872 updated = 255;
873 if(this->threadNumber_ > 1) {
874 vertLock[vertexId].unlock();
875 }
876
877 return reps[0];
878
879 } else {
880 if(debugLevel_ > 5)
881 printMsg("to non saddle " + std::to_string(vertexId) + " "
882 + std::to_string(saddleCC[vertexId].size()));
883
884 SimplexId ret = vertexId;
885 SimplexId neighborNumber
886 = multiresTriangulation_.getVertexNeighborNumber(vertexId);
887 SimplexId maxNeighbor = vertexId;
888 // std::cout << "neigh number" << neighborNumber << std::endl;
889 for(SimplexId i = 0; i < neighborNumber; i++) {
890 SimplexId neighborId = -1;
891 multiresTriangulation_.getVertexNeighbor(vertexId, i, neighborId);
892 if(gt(neighborId, maxNeighbor)) {
893 maxNeighbor = neighborId;
894 }
895 }
896 if(maxNeighbor != vertexId) { // not an extremum
897 ret = propagateFromSaddles(maxNeighbor, vertLock, toPropagate,
898 vertexRepresentatives, saddleCC, isUpdated,
899 globalExtremum, splitTree, fakeScalars,
900 offsets, monotonyOffsets);
901
902 } else { // needed to find the globalExtremum per thread
903#ifdef TTK_ENABLE_OPENMP
904 const auto tid = omp_get_thread_num();
905#else
906 const auto tid = 0;
907#endif // TTK_ENABLE_OPENMP
908 if(gt(vertexId, globalExtremum[tid])) {
909 // if(splitTree)
910 // std::cout << "new global max " << std::endl;
911 // else
912 // std::cout << "new global min " << std::endl;
913 globalExtremum[tid] = vertexId;
914 }
915 }
916 reps.resize(1);
917 reps[0] = ret;
918 updated = 255;
919 if(this->threadNumber_ > 1) {
920 vertLock[vertexId].unlock();
921 }
922 return ret;
923 }
924}
925
926template <typename scalarType, typename offsetType>
927void ttk::ApproximateTopology::updatePropagation(
928 std::vector<polarity> &toPropagateMin,
929 std::vector<polarity> &toPropagateMax,
930 std::vector<std::vector<SimplexId>> &vertexRepresentativesMin,
931 std::vector<std::vector<SimplexId>> &vertexRepresentativesMax,
932 std::vector<std::vector<SimplexId>> &saddleCCMin,
933 std::vector<std::vector<SimplexId>> &saddleCCMax,
934 std::vector<Lock> &vertLockMin,
935 std::vector<Lock> &vertLockMax,
936 std::vector<polarity> &isUpdatedMin,
937 std::vector<polarity> &isUpdatedMax,
938 const scalarType *fakeScalars,
939 const offsetType *const offsets,
940 const int *const monotonyOffsets) {
941
942 Timer tm{};
943 const size_t nDecVerts = multiresTriangulation_.getDecimatedVertexNumber();
944
945 if(debugLevel_ > 5) {
946 const auto pred = [](const polarity a) { return a > 0; };
947 const auto numberOfCandidatesToPropagateMax
948 = std::count_if(toPropagateMax.begin(), toPropagateMax.end(), pred);
949 std::cout << " sad-max we have " << numberOfCandidatesToPropagateMax
950 << " vertices to propagate from outta " << nDecVerts << std::endl;
951 const auto numberOfCandidatesToPropagateMin
952 = std::count_if(toPropagateMin.begin(), toPropagateMin.end(), pred);
953 std::cout << " min-sad we have " << numberOfCandidatesToPropagateMin
954 << " vertices to propagate from outta " << nDecVerts << std::endl;
955 }
956
957 std::vector<SimplexId> globalMaxThr(threadNumber_, 0);
958 std::vector<SimplexId> globalMinThr(threadNumber_, 0);
959
960 // reset updated flag
961#ifdef TTK_ENABLE_OPENMP
962#pragma omp parallel for num_threads(threadNumber_)
963#endif // TTK_ENABLE_OPENMP
964 for(size_t i = 0; i < nDecVerts; i++) {
965 SimplexId v = multiresTriangulation_.localToGlobalVertexId(i);
966 isUpdatedMin[v] = 0;
967 isUpdatedMax[v] = 0;
968 }
969
970 // propagate along split tree
971#ifdef TTK_ENABLE_OPENMP
972#pragma omp parallel for num_threads(threadNumber_)
973#endif // TTK_ENABLE_OPENMP
974 for(size_t i = 0; i < nDecVerts; i++) {
975 SimplexId v = multiresTriangulation_.localToGlobalVertexId(i);
976 if(toPropagateMin[v]) {
977 propagateFromSaddles(v, vertLockMin, toPropagateMin,
978 vertexRepresentativesMin, saddleCCMin, isUpdatedMin,
979 globalMinThr, false, fakeScalars, offsets,
980 monotonyOffsets);
981 }
982 if(toPropagateMax[v]) {
983 propagateFromSaddles(v, vertLockMax, toPropagateMax,
984 vertexRepresentativesMax, saddleCCMax, isUpdatedMax,
985 globalMaxThr, true, fakeScalars, offsets,
986 monotonyOffsets);
987 }
988 }
989
990 const auto lt = [=](const SimplexId a, const SimplexId b) -> bool {
991 return ((fakeScalars[a] < fakeScalars[b])
992 || (fakeScalars[a] == fakeScalars[b]
993 && ((monotonyOffsets[a] < monotonyOffsets[b])
994 || (monotonyOffsets[a] == monotonyOffsets[b]
995 && offsets[a] < offsets[b]))));
996 };
997
998 globalMin_ = *std::min_element(globalMinThr.begin(), globalMinThr.end(), lt);
999 globalMax_ = *std::max_element(globalMaxThr.begin(), globalMaxThr.end(), lt);
1000
1001 if(globalMin_ == 0 or globalMax_ == 0) {
1002#ifdef TTK_ENABLE_OPENMP
1003#pragma omp parallel for num_threads(threadNumber_)
1004#endif // TTK_ENABLE_OPENMP
1005 for(size_t i = 0; i < nDecVerts; i++) {
1006 SimplexId v = multiresTriangulation_.localToGlobalVertexId(i);
1007
1008#ifdef TTK_ENABLE_OPENMP
1009 const auto tid = omp_get_thread_num();
1010#else
1011 const auto tid = 0;
1012#endif // TTK_ENABLE_OPENMP
1013
1014 if(lt(globalMaxThr[tid], v)) {
1015 globalMaxThr[tid] = v;
1016 }
1017 if(lt(v, globalMinThr[tid])) {
1018 globalMinThr[tid] = v;
1019 }
1020 }
1021 globalMin_
1022 = *std::min_element(globalMinThr.begin(), globalMinThr.end(), lt);
1023 globalMax_
1024 = *std::max_element(globalMaxThr.begin(), globalMaxThr.end(), lt);
1025 // printMsg("Explicitly found global extremas");
1026 }
1027 if(debugLevel_ > 3) {
1028 printMsg("Propagation Update", 1, tm.getElapsedTime(), threadNumber_);
1029 }
1030}
1031
1032template <typename scalarType, typename offsetType>
1033void ttk::ApproximateTopology::buildVertexLinkPolarityApproximate(
1034 const SimplexId vertexId,
1035 std::vector<std::pair<polarity, polarity>> &vlp,
1036 const scalarType *fakeScalars,
1037 const offsetType *const offsets,
1038 const int *const monotonyOffsets) const {
1039
1040 const SimplexId neighborNumber
1041 = multiresTriangulation_.getVertexNeighborNumber(vertexId);
1042 vlp.resize(neighborNumber);
1043
1044 for(SimplexId i = 0; i < neighborNumber; i++) {
1045 SimplexId neighborId0 = 0;
1046 multiresTriangulation_.getVertexNeighbor(vertexId, i, neighborId0);
1047
1048 // const bool lower0 = vertsOrder_[neighborId0] < vertsOrder_[vertexId];
1049 const bool lower0
1050 = ((fakeScalars[neighborId0] < fakeScalars[vertexId])
1051 || (fakeScalars[neighborId0] == fakeScalars[vertexId]
1052 && ((monotonyOffsets[neighborId0] < monotonyOffsets[vertexId])
1053 || (monotonyOffsets[neighborId0] == monotonyOffsets[vertexId]
1054 && offsets[neighborId0] < offsets[vertexId]))));
1055
1056 const polarity isUpper0 = static_cast<polarity>(!lower0) * 255;
1057 vlp[i] = std::make_pair(isUpper0, 0);
1058 }
1059}
1060
1061template <typename scalarType, typename offsetType>
1062bool ttk::ApproximateTopology::printPolarity(
1063 std::vector<polarity> &isNew,
1064 const SimplexId vertexId,
1065 std::vector<std::vector<std::pair<polarity, polarity>>> &vertexLinkPolarity,
1066 const scalarType *scalars,
1067 const scalarType *fakeScalars,
1068 const offsetType *const offsets,
1069 const int *const monotonyOffsets,
1070 bool verbose) {
1071
1072 bool error = false;
1073 std::stringstream mymsg;
1074 std::vector<std::pair<polarity, polarity>> &vlp
1075 = vertexLinkPolarity[vertexId];
1076 mymsg << "POLARITY PRINT"
1077 << "\n";
1078 mymsg << "vertex " << vertexId << " has "
1079 << multiresTriangulation_.getVertexNeighborNumber(vertexId)
1080 << " neighbors"
1081 << "\n";
1082 mymsg << "\tself f:" << fakeScalars[vertexId] << " s:" << scalars[vertexId]
1083 << " o:" << offsets[vertexId] << " m:" << monotonyOffsets[vertexId]
1084 << " isnew: " << (int)isNew[vertexId] << "\n";
1085 if(vlp.empty()) {
1086 if(verbose) {
1087 std::cout << "\tpolarity not initialized for " << vertexId << std::endl;
1088 std::cout << mymsg.str() << std::endl;
1089 }
1090 return false;
1091 }
1092 for(SimplexId i = 0;
1093 i < multiresTriangulation_.getVertexNeighborNumber(vertexId); i++) {
1094 SimplexId nId = -1;
1095 multiresTriangulation_.getVertexNeighbor(vertexId, i, nId);
1096 // find reverse pol
1097 std::vector<std::pair<polarity, polarity>> &vlp2 = vertexLinkPolarity[nId];
1098 bool init = false;
1099 polarity rpol;
1100 if(!vlp2.empty()) {
1101 for(SimplexId j = 0;
1102 j < multiresTriangulation_.getVertexNeighborNumber(nId); j++) {
1103 SimplexId nId2 = -1;
1104 multiresTriangulation_.getVertexNeighbor(nId, j, nId2);
1105 if(nId2 == vertexId) {
1106 rpol = vlp2[j].first;
1107 init = true;
1108 }
1109 }
1110 }
1111
1112 const bool lower
1113 = ((fakeScalars[nId] < fakeScalars[vertexId])
1114 || (fakeScalars[nId] == fakeScalars[vertexId]
1115 && ((monotonyOffsets[nId] < monotonyOffsets[vertexId])
1116 || (monotonyOffsets[nId] == monotonyOffsets[vertexId]
1117 && offsets[nId] < offsets[vertexId]))));
1118
1119 const polarity isUpper = lower ? 0 : 255;
1120
1121 mymsg << " " << i << "th: " << nId << " f:" << fakeScalars[nId]
1122 << " s:" << scalars[nId] << " o:" << offsets[nId]
1123 << " m:" << monotonyOffsets[nId] << " , pol:" << (bool)vlp[i].first
1124 << "(" << (bool)vlp[i].second << ")"
1125 << " rpol:" << (bool)rpol << " true pol:" << (bool)isUpper
1126 << " init " << init << " isnew: " << (int)isNew[nId] << "\n";
1127 if((rpol == isUpper and !vlp2.empty())
1128 or (isUpper != vlp[i].first and !vlp[i].second)) {
1129 mymsg << "POLARITY ERROR "
1130 << "\n";
1131 error = true;
1132 }
1133 }
1134 if(error or verbose) {
1135 std::cout << mymsg.str() << std::endl;
1136 }
1137 return error;
1138}
1139
1140template <typename scalarType, typename offsetType>
1141void ttk::ApproximateTopology::getCriticalTypeApproximate(
1142 const SimplexId &vertexId,
1143 std::vector<std::pair<polarity, polarity>> &vlp,
1144 uint8_t &vertexLink,
1145 DynamicTree &link,
1146 VLBoundaryType &vlbt,
1147 const scalarType *fakeScalars,
1148 const offsetType *const offsets,
1149 const int *const monotonyOffsets) const {
1150
1151 if(vlp.empty()) {
1152 buildVertexLinkPolarityApproximate(
1153 vertexId, vlp, fakeScalars, offsets, monotonyOffsets);
1154 }
1155 const SimplexId neighborNumber
1156 = multiresTriangulation_.getVertexNeighborNumber(vertexId);
1157 link.alloc(neighborNumber);
1158
1159 // associate vertex link boundary
1160 vertexLink = multiresTriangulation_.getVertexBoundaryIndex(vertexId);
1161
1162 // update the link polarity for old points that are processed for
1163 // the first time
1164 const auto &vl = vlbt[vertexLink];
1165
1166 for(size_t edgeId = 0; edgeId < vl.size(); edgeId++) {
1167 const SimplexId n0 = vl[edgeId].first;
1168 const SimplexId n1 = vl[edgeId].second;
1169 if(vlp[n0].first == vlp[n1].first) {
1170 // the smallest id (n0) becomes the parent of n1
1171 link.insertEdge(n1, n0);
1172 }
1173 }
1174}
1175
1176template <typename scalarType, typename offsetType>
1177void ttk::ApproximateTopology::initGlobalPolarity(
1178 std::vector<polarity> &isNew,
1179 std::vector<std::vector<std::pair<polarity, polarity>>> &vertexLinkPolarity,
1180 std::vector<polarity> &toProcess,
1181 const scalarType *fakeScalars,
1182 const offsetType *const offsets,
1183 const int *const monotonyOffsets) const {
1184
1185 Timer timer{};
1186 const size_t nDecVerts = multiresTriangulation_.getDecimatedVertexNumber();
1187
1188 // computes the critical types of all points
1189#ifdef TTK_ENABLE_OPENMP
1190#pragma omp parallel for num_threads(threadNumber_)
1191#endif // TTK_ENABLE_OPENMP
1192 for(size_t i = 0; i < nDecVerts; i++) {
1193 SimplexId globalId = multiresTriangulation_.localToGlobalVertexId(i);
1194 buildVertexLinkPolarityApproximate(globalId, vertexLinkPolarity[globalId],
1195 fakeScalars, offsets, monotonyOffsets);
1196 toProcess[globalId] = 255;
1197 isNew[globalId] = 0;
1198 }
1199 printMsg("Polarity Init", 1, timer.getElapsedTime(), threadNumber_,
1200 debug::LineMode::NEW, debug::Priority::DETAIL);
1201}
1202
1203template <typename ScalarType, typename offsetType>
1204int ttk::ApproximateTopology::updateGlobalPolarity(
1205 double eps,
1206 std::vector<polarity> &isNew,
1207 std::vector<std::vector<std::pair<polarity, polarity>>> &vertexLinkPolarity,
1208 std::vector<polarity> &toProcess,
1209 std::vector<polarity> &toReprocess,
1210 ScalarType *fakeScalars,
1211 const offsetType *const offsets,
1212 int *monotonyOffsets) const {
1213
1214 Timer tm;
1215 const auto nDecVerts = multiresTriangulation_.getDecimatedVertexNumber();
1216
1217#ifdef TTK_ENABLE_OPENMP
1218#pragma omp parallel for num_threads(threadNumber_)
1219#endif // TTK_ENABLE_OPENMP
1220 for(SimplexId localId = 0; localId < nDecVerts; localId++) {
1221 SimplexId globalId = multiresTriangulation_.localToGlobalVertexId(localId);
1222 if(!isNew[globalId]) {
1223 getMonotonyChangeByOldPointCPApproximate(
1224 globalId, eps, isNew, toProcess, toReprocess,
1225 vertexLinkPolarity[globalId], fakeScalars, offsets, monotonyOffsets);
1226 }
1227 }
1228
1229 // second Loop process or reprocess
1230#ifdef TTK_ENABLE_OPENMP
1231#pragma omp parallel for num_threads(threadNumber_)
1232#endif
1233 for(int i = 0; i < multiresTriangulation_.getDecimatedVertexNumber(); i++) {
1234 SimplexId globalId = multiresTriangulation_.localToGlobalVertexId(i);
1235 if(isNew[globalId]) { // new point
1236 if(decimationLevel_ > stoppingDecimationLevel_) {
1237 buildVertexLinkPolarityApproximate(
1238 globalId, vertexLinkPolarity[globalId], fakeScalars, offsets,
1239 monotonyOffsets);
1240 }
1241 isNew[globalId] = 0;
1242
1243 } else { // old point
1244 if(toReprocess[globalId]) {
1245 updateLinkPolarityPonctual(vertexLinkPolarity[globalId]);
1246
1247 toProcess[globalId] = 255; // mark for processing
1248 toReprocess[globalId] = 0;
1249 }
1250 }
1251 } // end for openmp
1252 return 0;
1253}
1254
1255template <typename ScalarType, typename offsetType>
1256void ttk::ApproximateTopology::computeCriticalPoints(
1257 std::vector<std::vector<std::pair<polarity, polarity>>> &vertexLinkPolarity,
1258 std::vector<polarity> &toPropagateMin,
1259 std::vector<polarity> &toPropagateMax,
1260 std::vector<polarity> &toProcess,
1261 std::vector<DynamicTree> &link,
1262 std::vector<uint8_t> &vertexLink,
1263 VLBoundaryType &vertexLinkByBoundaryType,
1264 std::vector<std::vector<SimplexId>> &saddleCCMin,
1265 std::vector<std::vector<SimplexId>> &saddleCCMax,
1266 ScalarType *fakeScalars,
1267 const offsetType *const offsets,
1268 int *monotonyOffsets) {
1269
1270 Timer tm;
1271
1272#ifdef TTK_ENABLE_OPENMP
1273#pragma omp parallel for num_threads(threadNumber_)
1274#endif
1275 for(int i = 0; i < multiresTriangulation_.getDecimatedVertexNumber(); i++) {
1276 SimplexId globalId = multiresTriangulation_.localToGlobalVertexId(i);
1277
1278 if(toProcess[globalId]) {
1279 // nbComputations++;
1280 getCriticalTypeApproximate(globalId, vertexLinkPolarity[globalId],
1281 vertexLink[globalId], link[globalId],
1282 vertexLinkByBoundaryType, fakeScalars, offsets,
1283 monotonyOffsets);
1284 getValencesFromLink(globalId, vertexLinkPolarity[globalId],
1285 link[globalId], toPropagateMin, toPropagateMax,
1286 saddleCCMin, saddleCCMax);
1287 }
1288 } // end for openmp
1289
1290 if(debugLevel_ > 3) {
1291 printMsg(
1292 "Critical Points Computation", 1, tm.getElapsedTime(), threadNumber_);
1293 }
1294}
1295
1296template <typename scalarType>
1297int ttk::ApproximateTopology::sortPersistenceDiagramApproximate(
1298 std::vector<PersistencePair> &diagram,
1299 scalarType *fakeScalars,
1300 const SimplexId *const offsets,
1301 int *monotonyOffsets) const {
1302 auto cmp = [fakeScalars, offsets, monotonyOffsets](
1303 const PersistencePair &pA, const PersistencePair &pB) {
1304 const ttk::SimplexId a = pA.birth;
1305 const ttk::SimplexId b = pB.birth;
1306
1307 return ((fakeScalars[a] < fakeScalars[b])
1308 || (fakeScalars[a] == fakeScalars[b]
1309 && ((monotonyOffsets[a] < monotonyOffsets[b])
1310 || (monotonyOffsets[a] == monotonyOffsets[b]
1311 && offsets[a] < offsets[b]))));
1312 };
1313
1314 std::sort(diagram.begin(), diagram.end(), cmp);
1315
1316 return 0;
1317}
1318
1319template <typename scalarType>
1320int ttk::ApproximateTopology::computeApproximatePD(
1321 std::vector<PersistencePair> &CTDiagram,
1322 const scalarType *scalars,
1323 scalarType *const fakeScalars,
1324 SimplexId *const outputOffsets,
1325 int *const outputMonotonyOffsets) {
1326
1327 int ret = -1;
1328 std::stringstream ss;
1329 ss << "Approximate Persistence Diagram computation with "
1330 << debug::output::UNDERLINED << debug::output::YELLOW << epsilon_ * 100
1331 << "%" << debug::output::ENDCOLOR << debug::output::ENDCOLOR << " error";
1332 printMsg(ss.str());
1333
1334 ret = executeApproximateTopology(
1335 scalars, fakeScalars, outputOffsets, outputMonotonyOffsets);
1336 CTDiagram = std::move(CTDiagram_);
1337 CTDiagram_.clear();
1338
1339 return ret;
1340}
ttkNotUsed
#define ttkNotUsed(x)
Mark function/method parameters that are not used in the function body at all.
Definition: BaseClass.h:47
TTK_PSORT
#define TTK_PSORT(NTHREADS,...)
Parallel sort macro.
Definition: OpenMP.h:46
ttk::ApproximateTopology
TTK processing package for progressive Topological Data Analysis.
Definition: ApproximateTopology.h:36
ttk::ApproximateTopology::initSaddleSeeds
void initSaddleSeeds(std::vector< polarity > &isNew, std::vector< std::vector< std::pair< polarity, polarity > > > &vertexLinkPolarity, std::vector< polarity > &toPropagateMin, std::vector< polarity > &toPropagateMax, std::vector< polarity > &toProcess, std::vector< DynamicTree > &link, std::vector< uint8_t > &vertexLink, VLBoundaryType &vertexLinkByBoundaryType, std::vector< std::vector< SimplexId > > &saddleCCMin, std::vector< std::vector< SimplexId > > &saddleCCMax, const SimplexId *const offsets) const
ttk::ApproximateTopology::updatePropagation
void updatePropagation(std::vector< polarity > &toPropagateMin, std::vector< polarity > &toPropagateMax, std::vector< std::vector< SimplexId > > &vertexRepresentativesMin, std::vector< std::vector< SimplexId > > &vertexRepresentativesMax, std::vector< std::vector< SimplexId > > &saddleCCMin, std::vector< std::vector< SimplexId > > &saddleCCMax, std::vector< Lock > &vertLockMin, std::vector< Lock > &vertLockMax, std::vector< polarity > &isUpdatedMin, std::vector< polarity > &isUpdatedMax, const scalarType *fakeScalars, const offsetType *const offsets, const int *const monotonyOffsets) const
ttk::ApproximateTopology::buildVertexLinkPolarityApproximate
void buildVertexLinkPolarityApproximate(const SimplexId vertexId, std::vector< std::pair< polarity, polarity > > &vlp, const scalarType *fakeScalars, const offsetType *const offsetField, const int *const monotonyOffsets) const
Definition: ApproximateTopology.h:1033
ttk::ApproximateTopology::buildVertexLinkPolarity
void buildVertexLinkPolarity(const SimplexId vertexId, std::vector< std::pair< polarity, polarity > > &vlp, const scalarType *fakeScalars, const offsetType *const offsets, const int *const monotonyOffsets) const
ttk::ApproximateTopology::sortTriplets
void sortTriplets(std::vector< triplet > &triplets, const ScalarType *const fakeScalars, const OffsetType *const offsets, const int *const monotonyOffsets, const bool splitTree) const
Definition: ApproximateTopology.h:520
ttk::ApproximateTopology::getMonotonyChangeByOldPointCPApproximate
int getMonotonyChangeByOldPointCPApproximate(const SimplexId vertexId, double eps, const std::vector< polarity > &isNew, std::vector< polarity > &toProcess, std::vector< polarity > &toReprocess, std::vector< std::pair< polarity, polarity > > &vlp, scalarType *fakeScalars, const offsetType *const offsets, int *monotonyOffsets) const
Definition: ApproximateTopology.h:673
ttk::ApproximateTopology::computePersistencePairsFromSaddles
void computePersistencePairsFromSaddles(std::vector< PersistencePair > &CTDiagram, const ScalarType *const fakeScalars, const offsetType *const offsets, const int *const monotonyOffsets, std::vector< std::vector< SimplexId > > &vertexRepresentativesMin, std::vector< std::vector< SimplexId > > &vertexRepresentativesMax, const std::vector< polarity > &toPropagateMin, const std::vector< polarity > &toPropagateMax) const
ttk::ApproximateTopology::initGlobalPolarity
void initGlobalPolarity(std::vector< polarity > &isNew, std::vector< std::vector< std::pair< polarity, polarity > > > &vertexLinkPolarity, std::vector< polarity > &toProcess, const scalarType *fakeScalars, const offsetType *const offsets, const int *const monotonyOffsets) const
Definition: ApproximateTopology.h:1177
ttk::ApproximateTopology::VLBoundaryType
std::array< std::vector< std::pair< SimplexId, SimplexId > >, nLink_ > VLBoundaryType
Definition: ApproximateTopology.h:173
ttk::ApproximateTopology::initCriticalPoints
void initCriticalPoints(std::vector< polarity > &isNew, std::vector< std::vector< std::pair< polarity, polarity > > > &vertexLinkPolarity, std::vector< polarity > &toProcess, std::vector< DynamicTree > &link, std::vector< uint8_t > &vertexLink, VLBoundaryType &vertexLinkByBoundaryType, std::vector< char > &vertexTypes, const SimplexId *const offsets) const
ttk::ApproximateTopology::tripletsToPersistencePairs
void tripletsToPersistencePairs(std::vector< PersistencePair > &pairs, std::vector< std::vector< SimplexId > > &vertexRepresentatives, std::vector< triplet > &triplets, const scalarType *const fakeScalars, const offsetType *const offsets, const int *const monotonyOffsets, const bool splitTree) const
Definition: ApproximateTopology.h:552
ttk::ApproximateTopology::sortTripletsApproximate
void sortTripletsApproximate(std::vector< triplet > &triplets, const ScalarType *const scalars, const ScalarType *const fakeScalars, const OffsetType *const offsets, const int *const monotonyOffsets, const bool splitTree) const
ttk::ApproximateTopology::executeApproximateTopology
int executeApproximateTopology(const scalarType *scalars, scalarType *fakeScalars, SimplexId *outputOffsets, int *outputMonotonyOffsets)
ttk::ApproximateTopology::sortPersistenceDiagramApproximate
int sortPersistenceDiagramApproximate(std::vector< PersistencePair > &diagram, scalarType *fakeScalars, const SimplexId *const offsets, int *monotonyOffsets) const
Definition: ApproximateTopology.h:1297
ttk::ApproximateTopology::sortVertices
void sortVertices(const SimplexId vertexNumber, std::vector< SimplexId > &sortedVertices, SimplexId *vertsOrder, const scalarType *const fakeScalars, const offsetType *const offsetField, const int *const monotonyOffsets)
Definition: ApproximateTopology.h:635
ttk::ApproximateTopology::computeApproximatePD
int computeApproximatePD(std::vector< PersistencePair > &CTDiagram, const scalarType *scalars, scalarType *const fakeScalars, SimplexId *const outputOffsets, int *const outputMonotonyOffsets)
Definition: ApproximateTopology.h:1320
ttk::ApproximateTopology::updatePropagation
void updatePropagation(std::vector< polarity > &toPropagateMin, std::vector< polarity > &toPropagateMax, std::vector< std::vector< SimplexId > > &vertexRepresentativesMin, std::vector< std::vector< SimplexId > > &vertexRepresentativesMax, std::vector< std::vector< SimplexId > > &saddleCCMin, std::vector< std::vector< SimplexId > > &saddleCCMax, std::vector< Lock > &vertLockMin, std::vector< Lock > &vertLockMax, std::vector< polarity > &isUpdatedMin, std::vector< polarity > &isUpdatedMax, const scalarType *fakeScalars, const offsetType *const offsetField, const int *const monotonyOffsets)
Definition: ApproximateTopology.h:927
ttk::ApproximateTopology::printPolarity
bool printPolarity(std::vector< polarity > &isNew, const SimplexId vertexId, std::vector< std::vector< std::pair< polarity, polarity > > > &vertexLinkPolarity, const scalarType *scalars, const scalarType *fakeScalars, const offsetType *const offsets, const int *const monotonyOffsets, bool verbose=false)
Definition: ApproximateTopology.h:1062
ttk::ApproximateTopology::getCriticalTypeApproximate
void getCriticalTypeApproximate(const SimplexId &vertexId, std::vector< std::pair< polarity, polarity > > &vlp, uint8_t &vertexLink, DynamicTree &link, VLBoundaryType &vlbt, const scalarType *fakeScalars, const offsetType *const offsets, const int *const monotonyOffsets) const
Definition: ApproximateTopology.h:1141
ttk::ApproximateTopology::propagateFromSaddles
ttk::SimplexId propagateFromSaddles(const SimplexId vertexId, std::vector< Lock > &vertLock, std::vector< polarity > &toPropagate, std::vector< std::vector< SimplexId > > &vertexRepresentatives, std::vector< std::vector< SimplexId > > &saddleCC, std::vector< polarity > &isUpdated, std::vector< SimplexId > &globalExtremum, const bool splitTree, const scalarType *fakeScalars, const offsetType *const offsetField, const int *const monotonyOffsets) const
Definition: ApproximateTopology.h:794
ttk::ApproximateTopology::propagateFromSaddles
ttk::SimplexId propagateFromSaddles(const SimplexId vertexId, std::vector< Lock > &vertLock, std::vector< polarity > &toPropagate, std::vector< std::vector< SimplexId > > &vertexRepresentatives, std::vector< std::vector< SimplexId > > &saddleCC, std::vector< polarity > &isUpdated, std::vector< SimplexId > &globalExtremum, const SimplexId *const offsets, const bool splitTree) const
ttk::ApproximateTopology::nLink_
static const size_t nLink_
Definition: ApproximateTopology.h:171
ttk::ApproximateTopology::setEpsilon
void setEpsilon(double data)
Definition: ApproximateTopology.h:42
ttk::ApproximateTopology::computeCriticalPoints
void computeCriticalPoints(std::vector< std::vector< std::pair< polarity, polarity > > > &vertexLinkPolarity, std::vector< polarity > &toPropagateMin, std::vector< polarity > &toPropagateMax, std::vector< polarity > &toProcess, std::vector< DynamicTree > &link, std::vector< uint8_t > &vertexLink, VLBoundaryType &vertexLinkByBoundaryType, std::vector< std::vector< SimplexId > > &saddleCCMin, std::vector< std::vector< SimplexId > > &saddleCCMax, ScalarType *fakeScalars, const offsetType *const offsets, int *monotonyOffsets)
Definition: ApproximateTopology.h:1256
ttk::ApproximateTopology::computePersistencePairsFromSaddles
void computePersistencePairsFromSaddles(std::vector< PersistencePair > &CTDiagram, const SimplexId *const offsets, std::vector< std::vector< SimplexId > > &vertexRepresentativesMin, std::vector< std::vector< SimplexId > > &vertexRepresentativesMax, const std::vector< polarity > &toPropagateMin, const std::vector< polarity > &toPropagateMax) const
ttk::ApproximateTopology::updateGlobalPolarity
int updateGlobalPolarity(double eps, std::vector< polarity > &isNew, std::vector< std::vector< std::pair< polarity, polarity > > > &vertexLinkPolarity, std::vector< polarity > &toProcess, std::vector< polarity > &toReprocess, ScalarType *fakeScalars, const offsetType *const offsets, int *monotonyOffsets) const
Definition: ApproximateTopology.h:1204
ttk::BaseClass::threadNumber_
int threadNumber_
Definition: BaseClass.h:95
ttk::Debug::setDebugMsgPrefix
void setDebugMsgPrefix(const std::string &prefix)
Definition: Debug.h:364
ttk::Debug::printMsg
int printMsg(const std::string &msg, const debug::Priority &priority=debug::Priority::INFO, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cout) const
Definition: Debug.h:118
ttk::Debug::printErr
int printErr(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition: Debug.h:149
ttk::DynamicTree
Implements the Dynamic Tree data-structure (or ST-Tree)
Definition: DynamicTree.h:66
ttk::DynamicTree::alloc
void alloc(const std::size_t nbNodes)
Definition: DynamicTree.h:73
ttk::DynamicTree::insertEdge
bool insertEdge(DynTreeNode *const n1, DynTreeNode *const n2)
Definition: DynamicTree.h:93
ttk::MultiresTopology
TTK processing package for progressive Topological Data Analysis.
Definition: MultiresTopology.h:50
ttk::MultiresTopology::multiresTriangulation_
MultiresTriangulation multiresTriangulation_
Definition: MultiresTopology.h:134
ttk::MultiresTopology::CTDiagram_
std::vector< PersistencePair > CTDiagram_
Definition: MultiresTopology.h:144
ttk::MultiresTopology::buildVertexLinkByBoundary
void buildVertexLinkByBoundary(const SimplexId vertexId, VLBoundaryType &vlbt) const
Definition: MultiresTopology.cpp:53
ttk::MultiresTopology::triangulation_
ImplicitTriangulation * triangulation_
Definition: MultiresTopology.h:132
ttk::MultiresTriangulation::setTriangulation
void setTriangulation(ImplicitTriangulation *triangulation)
Definition: MultiresTriangulation.h:809
ttk::MultiresTriangulation::setDecimationLevel
void setDecimationLevel(int decimationLevel)
Definition: MultiresTriangulation.h:882
ttk::MultiresTriangulation::findBoundaryRepresentatives
void findBoundaryRepresentatives(std::vector< SimplexId > &boundaryRepresentatives)
Definition: MultiresTriangulation.cpp:4314
ttk::Timer::getElapsedTime
double getElapsedTime()
Definition: Timer.h:15
ttk::debug::output::ENDCOLOR
const std::string ENDCOLOR
Definition: Debug.h:82
ttk::debug::output::YELLOW
const std::string YELLOW
Definition: Debug.h:77
ttk::debug::output::UNDERLINED
const std::string UNDERLINED
Definition: Debug.h:70
ttk
The Topology ToolKit.
Definition: AbstractTriangulation.h:51
ttk::triplet
std::tuple< ttk::SimplexId, ttk::SimplexId, ttk::SimplexId > triplet
Persistence pair type (with persistence in double)
Definition: ApproximateTopology.h:33
ttk::polarity
unsigned char polarity
Definition: ApproximateTopology.h:34
ttk::SimplexId
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22
printMsg
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)
Generated on Tue Jul 11 2023 08:20:52 for TTK by doxygen 1.9.6