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 const 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 const delta = epsilon_ * delta_;
404
405 while(decimationLevel_ > stoppingDecimationLevel_) {
406 decimationLevel_--;
407 multiresTriangulation_.setDecimationLevel(decimationLevel_);
408 // double progress = (double)(startingDecimationLevel_ - decimationLevel_)
409 // / (startingDecimationLevel_);
410 // this->printMsg("decimation level: " + std::to_string(decimationLevel_),
411 // progress, timer.getElapsedTime() - tm_allocation,
412 // threadNumber_);
413
414 int const ret = updateGlobalPolarity(delta, isNew, vertexLinkPolarity,
415 toProcess, toReprocess, fakeScalars,
416 offsets, monotonyOffsets);
417 if(ret == -1) {
418 std::cout << "Found ERROR - aborting" << std::endl;
419 return -1;
420 }
421 } // end while
422
423 computeCriticalPoints(vertexLinkPolarity, toPropagateMin, toPropagateMax,
424 toProcess, link, vertexLink, vertexLinkByBoundaryType,
425 saddleCCMin, saddleCCMax, fakeScalars, offsets,
426 monotonyOffsets);
427
428 updatePropagation(toPropagateMin, toPropagateMax, vertexRepresentativesMin,
429 vertexRepresentativesMax, saddleCCMin, saddleCCMax,
430 vertLockMin, vertLockMax, isUpToDateMin, isUpToDateMax,
431 fakeScalars, offsets, monotonyOffsets);
432
433 computePersistencePairsFromSaddles(
434 CTDiagram_, fakeScalars, offsets, monotonyOffsets, vertexRepresentativesMin,
435 vertexRepresentativesMax, toPropagateMin, toPropagateMax);
436 // ADD GLOBAL MIN-MAX PAIR
437 CTDiagram_.emplace_back(this->globalMin_, this->globalMax_, -1);
438 // fakeScalars[this->globalMax_] - fakeScalars[this->globalMin_], -1);
439
440 // std::cout << "# epsilon " << epsilon_ << " in "
441 // << timer.getElapsedTime() - tm_allocation << " s." << std::endl;
442
443 this->printMsg("Complete", 1.0, timer.getElapsedTime() - tm_allocation,
444 this->threadNumber_);
445
446 // finally sort the diagram
447 // std::cout << "SORTING" << std::endl;
448 sortPersistenceDiagramApproximate(
449 CTDiagram_, fakeScalars, offsets, monotonyOffsets);
450 // std::cout << "Final epsilon " << epsilon_ << std::endl;
451
452 // sort vertices to generate correct output offset order
453 std::vector<SimplexId> sortedVertices{};
454 sortVertices(vertexNumber, sortedVertices, vertsOrder, fakeScalars, offsets,
455 monotonyOffsets);
456 return 0;
457}
458
459template <typename ScalarType, typename OffsetType>
460void ttk::ApproximateTopology::computePersistencePairsFromSaddles(
461 std::vector<PersistencePair> &CTDiagram,
462 const ScalarType *const fakeScalars,
463 const OffsetType *const offsets,
464 const int *const monotonyOffsets,
465 std::vector<std::vector<SimplexId>> &vertexRepresentativesMin,
466 std::vector<std::vector<SimplexId>> &vertexRepresentativesMax,
467 const std::vector<polarity> &toPropagateMin,
468 const std::vector<polarity> &toPropagateMax) const {
469
470 Timer timer{};
471 // CTDiagram.clear();
472 std::vector<triplet> tripletsMax{}, tripletsMin{};
473 const SimplexId nbDecVert = multiresTriangulation_.getDecimatedVertexNumber();
474
475 for(SimplexId localId = 0; localId < nbDecVert; localId++) {
476 const SimplexId globalId
477 = 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 const 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 const 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 sortedVertices.resize(vertexNumber);
644
645 // fill with numbers from 0 to vertexNumber - 1
646 std::iota(sortedVertices.begin(), sortedVertices.end(), 0);
647
648 // sort vertices in ascending order following scalarfield / offsets
649 TTK_PSORT(this->threadNumber_, sortedVertices.begin(), sortedVertices.end(),
650 [&](const SimplexId a, const SimplexId b) {
651 return ((fakeScalars[a] < fakeScalars[b])
652 || (fakeScalars[a] == fakeScalars[b]
653 && ((monotonyOffsets[a] < monotonyOffsets[b])
654 || (monotonyOffsets[a] == monotonyOffsets[b]
655 && offsets[a] < offsets[b]))));
656 });
657
658#ifdef TTK_ENABLE_OPENMP
659#pragma omp parallel for num_threads(threadNumber_)
660#endif // TTK_ENABLE_OPENMP
661 for(size_t i = 0; i < sortedVertices.size(); ++i) {
662 vertsOrder[sortedVertices[i]] = i;
663 }
664}
665
666template <typename scalarType, typename offsetType>
667int ttk::ApproximateTopology::getMonotonyChangeByOldPointCPApproximate(
668 const SimplexId vertexId,
669 double eps,
670 const std::vector<polarity> &isNew,
671 std::vector<polarity> &toProcess,
672 std::vector<polarity> &toReprocess,
673 std::vector<std::pair<polarity, polarity>> &vlp,
674 scalarType *fakeScalars,
675 const offsetType *const offsets,
676 int *monotonyOffsets) const {
677
678 int hasMonotonyChanged = 0;
679 const SimplexId neighborNumber
680 = multiresTriangulation_.getVertexNeighborNumber(vertexId);
681 for(SimplexId i = 0; i < neighborNumber; i++) {
682 SimplexId neighborId = -1;
683 multiresTriangulation_.getVertexNeighbor(vertexId, i, neighborId);
684
685 // check for monotony changes
686 // const bool lowerStatic
687 // = (scalarField[neighborId] < fakeScalars[vertexId])
688 // || ((scalarField[neighborId] == fakeScalars[vertexId])
689 // && (offsets[neighborId] < offsets[vertexId]));
690 const bool lowerDynamic
691 = ((fakeScalars[neighborId] < fakeScalars[vertexId])
692 || (fakeScalars[neighborId] == fakeScalars[vertexId]
693 && ((monotonyOffsets[neighborId] < monotonyOffsets[vertexId])
694 || (monotonyOffsets[neighborId] == monotonyOffsets[vertexId]
695 && offsets[neighborId] < offsets[vertexId]))));
696
697 const polarity isUpperDynamic = lowerDynamic ? 0 : 255;
698 // const polarity isUpperStatic = lowerStatic ? 0 : 255;
699 const polarity isUpperOld = vlp[i].first;
700
701 if(isUpperDynamic != isUpperOld) { // change of monotony
702 SimplexId oldNeighbor = -1;
703 const int oldDecimation = pow(2, decimationLevel_ + 1);
704 multiresTriangulation_.getVertexNeighborAtDecimation(
705 vertexId, i, oldNeighbor, oldDecimation);
706
707 double const replacementValueDynamic
708 = (0.5 * (double)fakeScalars[oldNeighbor]
709 + .5 * (double)fakeScalars[vertexId]); // depends on epsilon
710 double const deltaDynamic
711 = fabs((double)fakeScalars[neighborId] - replacementValueDynamic);
712
713 //=====================
714 SimplexId oldNeighNumber = 0;
715 const SimplexId nnumber
716 = multiresTriangulation_.getVertexNeighborNumber(neighborId);
717 for(SimplexId iii = 0; iii < nnumber; iii++) {
718 SimplexId neighborId2 = -1;
719 multiresTriangulation_.getVertexNeighbor(neighborId, iii, neighborId2);
720 if(!isNew[neighborId2]) {
721 oldNeighNumber++;
722 }
723 }
724
725 if(deltaDynamic > eps or !isNew[neighborId] or oldNeighNumber > 2) {
726 hasMonotonyChanged = 1;
727
728 toReprocess[vertexId] = 255;
729 if(isNew[neighborId]) {
730 toProcess[neighborId] = 255;
731 } else {
732 toReprocess[neighborId] = 255;
733 }
734 const SimplexId neighborNumberNew
735 = multiresTriangulation_.getVertexNeighborNumber(neighborId);
736 for(SimplexId j = 0; j < neighborNumberNew; j++) {
737 SimplexId neighborIdNew = -1;
738 multiresTriangulation_.getVertexNeighbor(
739 neighborId, j, neighborIdNew);
740 if(isNew[neighborIdNew])
741 toProcess[neighborIdNew] = 255;
742 }
743 vlp[i].second = 255;
744 } else {
745 fakeScalars[neighborId] = replacementValueDynamic;
746
747 // corrects rounding error when we process an integer scalar field
748 if(fakeScalars[neighborId] == fakeScalars[oldNeighbor]) {
749 fakeScalars[neighborId] = fakeScalars[vertexId];
750 }
751
752 // change monotony offset
753 // The monotony must be preserved, meaning that
754 // if we had f(vertexId)>f(oldNeighbor)
755 // we should have f(vertexId)>f(neighborId)
756 // which is enforced when
757 // monotonyOffsets[vertexId]>monotonyOffsets[neighborId]
758 if(isUpperOld) { // we should enforce f(vertexId)<f(neighborId)
759 if(offsets[vertexId] > offsets[neighborId]) {
760 monotonyOffsets[neighborId]
761 = monotonyOffsets[vertexId] + pow(2, decimationLevel_);
762 if(monotonyOffsets[vertexId] == monotonyOffsets[oldNeighbor]
763 and fakeScalars[vertexId] == fakeScalars[oldNeighbor]) {
764 std::cout << "THIS IS AN ISSUE" << std::endl;
765 }
766 } else {
767 monotonyOffsets[neighborId] = monotonyOffsets[vertexId];
768 }
769 } else { // we should enforce f(vertexId)>f(neighborId)
770 if(offsets[vertexId] < offsets[neighborId]) {
771 monotonyOffsets[neighborId]
772 = monotonyOffsets[vertexId] - pow(2, decimationLevel_);
773 if(monotonyOffsets[vertexId] == monotonyOffsets[oldNeighbor]
774 and fakeScalars[vertexId] == fakeScalars[oldNeighbor]) {
775 std::cout << "THIS IS AN ISSUE" << std::endl;
776 }
777 } else {
778 monotonyOffsets[neighborId] = monotonyOffsets[vertexId];
779 }
780 }
781 }
782 } // end if change of monotony
783 } // end for neighbors
784 return hasMonotonyChanged;
785}
786
787template <typename scalarType, typename offsetType>
788ttk::SimplexId ttk::ApproximateTopology::propagateFromSaddles(
789 const SimplexId vertexId,
790 std::vector<Lock> &vertLock,
791 std::vector<polarity> &toPropagate,
792 std::vector<std::vector<SimplexId>> &vertexRepresentatives,
793 std::vector<std::vector<SimplexId>> &saddleCC,
794 std::vector<polarity> &isUpdated,
795 std::vector<SimplexId> &globalExtremum,
796 const bool splitTree,
797 const scalarType *fakeScalars,
798 const offsetType *const offsets,
799 const int *const monotonyOffsets) const {
800
801 auto &toProp = toPropagate[vertexId];
802 auto &reps = vertexRepresentatives[vertexId];
803 auto &updated = isUpdated[vertexId];
804
805 if(updated) {
806 return reps[0];
807 }
808
809 // const auto gt = [=](const SimplexId v1, const SimplexId v2) {
810 // return ((fakeScalars[v1] > fakeScalars[v2])
811 // || (fakeScalars[v1] == fakeScalars[v2]
812 // && offsets[v1] > offsets[v2]))
813 // == splitTree;
814 // };
815 const auto gt = [=](const SimplexId v1, const SimplexId v2) {
816 return ((fakeScalars[v1] > fakeScalars[v2])
817 || (fakeScalars[v1] == fakeScalars[v2]
818 && ((monotonyOffsets[v1] > monotonyOffsets[v2])
819 || (monotonyOffsets[v1] == monotonyOffsets[v2]
820 && offsets[v1] > offsets[v2]))))
821 == splitTree;
822 };
823
824 if(this->threadNumber_ > 1) {
825 vertLock[vertexId].lock();
826 }
827 if(saddleCC[vertexId].size()
828 and !toProp) { // tis a saddle point, should have to propagate on it
829 printErr("ERRRROR");
830 }
831 if(toProp) { // SADDLE POINT
832 if(debugLevel_ > 5)
833 printMsg("to saddle " + std::to_string(vertexId) + " "
834 + std::to_string(saddleCC[vertexId].size()));
835 // for(auto ccid : saddleCC[vertexId]) {
836 // SimplexId ccV = -1;
837 // multiresTriangulation_.getVertexNeighbor(vertexId, ccid, ccV);
838 // std::cout << " " << ccid << " (" << ccV << ") , ";
839 // }
840 // std::cout << std::endl;
841 const auto &CC = saddleCC[vertexId];
842 reps.clear();
843 reps.reserve(CC.size());
844 for(size_t r = 0; r < CC.size(); r++) {
845 SimplexId neighborId = -1;
846 const SimplexId localId = CC[r];
847 multiresTriangulation_.getVertexNeighbor(vertexId, localId, neighborId);
848 // printMsg(" CC " + std::to_string(CC[r]) + " "
849 // + std::to_string(neighborId) + " ("
850 // + std::to_string(fakeScalars[neighborId]) + " , "
851 // + std::to_string(offsets[neighborId]) + ")");
852 SimplexId const ret = propagateFromSaddles(
853 neighborId, vertLock, toPropagate, vertexRepresentatives, saddleCC,
854 isUpdated, globalExtremum, splitTree, fakeScalars, offsets,
855 monotonyOffsets);
856 reps.emplace_back(ret);
857 }
858
859 if(reps.size() > 1) {
860 // sort & remove duplicate elements
861 std::sort(reps.begin(), reps.end(), gt);
862 const auto last = std::unique(reps.begin(), reps.end());
863 reps.erase(last, reps.end());
864 }
865
866 updated = 255;
867 if(this->threadNumber_ > 1) {
868 vertLock[vertexId].unlock();
869 }
870
871 return reps[0];
872
873 } else {
874 if(debugLevel_ > 5)
875 printMsg("to non saddle " + std::to_string(vertexId) + " "
876 + std::to_string(saddleCC[vertexId].size()));
877
878 SimplexId ret = vertexId;
879 const SimplexId neighborNumber
880 = multiresTriangulation_.getVertexNeighborNumber(vertexId);
881 SimplexId maxNeighbor = vertexId;
882 // std::cout << "neigh number" << neighborNumber << std::endl;
883 for(SimplexId i = 0; i < neighborNumber; i++) {
884 SimplexId neighborId = -1;
885 multiresTriangulation_.getVertexNeighbor(vertexId, i, neighborId);
886 if(gt(neighborId, maxNeighbor)) {
887 maxNeighbor = neighborId;
888 }
889 }
890 if(maxNeighbor != vertexId) { // not an extremum
891 ret = propagateFromSaddles(maxNeighbor, vertLock, toPropagate,
892 vertexRepresentatives, saddleCC, isUpdated,
893 globalExtremum, splitTree, fakeScalars,
894 offsets, monotonyOffsets);
895
896 } else { // needed to find the globalExtremum per thread
897#ifdef TTK_ENABLE_OPENMP
898 const auto tid = omp_get_thread_num();
899#else
900 const auto tid = 0;
901#endif // TTK_ENABLE_OPENMP
902 if(gt(vertexId, globalExtremum[tid])) {
903 // if(splitTree)
904 // std::cout << "new global max " << std::endl;
905 // else
906 // std::cout << "new global min " << std::endl;
907 globalExtremum[tid] = vertexId;
908 }
909 }
910 reps.resize(1);
911 reps[0] = ret;
912 updated = 255;
913 if(this->threadNumber_ > 1) {
914 vertLock[vertexId].unlock();
915 }
916 return ret;
917 }
918}
919
920template <typename scalarType, typename offsetType>
921void ttk::ApproximateTopology::updatePropagation(
922 std::vector<polarity> &toPropagateMin,
923 std::vector<polarity> &toPropagateMax,
924 std::vector<std::vector<SimplexId>> &vertexRepresentativesMin,
925 std::vector<std::vector<SimplexId>> &vertexRepresentativesMax,
926 std::vector<std::vector<SimplexId>> &saddleCCMin,
927 std::vector<std::vector<SimplexId>> &saddleCCMax,
928 std::vector<Lock> &vertLockMin,
929 std::vector<Lock> &vertLockMax,
930 std::vector<polarity> &isUpdatedMin,
931 std::vector<polarity> &isUpdatedMax,
932 const scalarType *fakeScalars,
933 const offsetType *const offsets,
934 const int *const monotonyOffsets) {
935
936 Timer tm{};
937 const size_t nDecVerts = multiresTriangulation_.getDecimatedVertexNumber();
938
939 if(debugLevel_ > 5) {
940 const auto pred = [](const polarity a) { return a > 0; };
941 const auto numberOfCandidatesToPropagateMax
942 = std::count_if(toPropagateMax.begin(), toPropagateMax.end(), pred);
943 std::cout << " sad-max we have " << numberOfCandidatesToPropagateMax
944 << " vertices to propagate from outta " << nDecVerts << std::endl;
945 const auto numberOfCandidatesToPropagateMin
946 = std::count_if(toPropagateMin.begin(), toPropagateMin.end(), pred);
947 std::cout << " min-sad we have " << numberOfCandidatesToPropagateMin
948 << " vertices to propagate from outta " << nDecVerts << std::endl;
949 }
950
951 std::vector<SimplexId> globalMaxThr(threadNumber_, 0);
952 std::vector<SimplexId> globalMinThr(threadNumber_, 0);
953
954 // reset updated flag
955#ifdef TTK_ENABLE_OPENMP
956#pragma omp parallel for num_threads(threadNumber_)
957#endif // TTK_ENABLE_OPENMP
958 for(size_t i = 0; i < nDecVerts; i++) {
959 const SimplexId v = multiresTriangulation_.localToGlobalVertexId(i);
960 isUpdatedMin[v] = 0;
961 isUpdatedMax[v] = 0;
962 }
963
964 // propagate along split tree
965#ifdef TTK_ENABLE_OPENMP
966#pragma omp parallel for num_threads(threadNumber_)
967#endif // TTK_ENABLE_OPENMP
968 for(size_t i = 0; i < nDecVerts; i++) {
969 SimplexId const v = multiresTriangulation_.localToGlobalVertexId(i);
970 if(toPropagateMin[v]) {
971 propagateFromSaddles(v, vertLockMin, toPropagateMin,
972 vertexRepresentativesMin, saddleCCMin, isUpdatedMin,
973 globalMinThr, false, fakeScalars, offsets,
974 monotonyOffsets);
975 }
976 if(toPropagateMax[v]) {
977 propagateFromSaddles(v, vertLockMax, toPropagateMax,
978 vertexRepresentativesMax, saddleCCMax, isUpdatedMax,
979 globalMaxThr, true, fakeScalars, offsets,
980 monotonyOffsets);
981 }
982 }
983
984 const auto lt = [=](const SimplexId a, const SimplexId b) -> bool {
985 return ((fakeScalars[a] < fakeScalars[b])
986 || (fakeScalars[a] == fakeScalars[b]
987 && ((monotonyOffsets[a] < monotonyOffsets[b])
988 || (monotonyOffsets[a] == monotonyOffsets[b]
989 && offsets[a] < offsets[b]))));
990 };
991
992 globalMin_ = *std::min_element(globalMinThr.begin(), globalMinThr.end(), lt);
993 globalMax_ = *std::max_element(globalMaxThr.begin(), globalMaxThr.end(), lt);
994
995 if(globalMin_ == 0 or globalMax_ == 0) {
996#ifdef TTK_ENABLE_OPENMP
997#pragma omp parallel for num_threads(threadNumber_)
998#endif // TTK_ENABLE_OPENMP
999 for(size_t i = 0; i < nDecVerts; i++) {
1000 SimplexId const v = multiresTriangulation_.localToGlobalVertexId(i);
1001
1002#ifdef TTK_ENABLE_OPENMP
1003 const auto tid = omp_get_thread_num();
1004#else
1005 const auto tid = 0;
1006#endif // TTK_ENABLE_OPENMP
1007
1008 if(lt(globalMaxThr[tid], v)) {
1009 globalMaxThr[tid] = v;
1010 }
1011 if(lt(v, globalMinThr[tid])) {
1012 globalMinThr[tid] = v;
1013 }
1014 }
1015 globalMin_
1016 = *std::min_element(globalMinThr.begin(), globalMinThr.end(), lt);
1017 globalMax_
1018 = *std::max_element(globalMaxThr.begin(), globalMaxThr.end(), lt);
1019 // printMsg("Explicitly found global extremas");
1020 }
1021 if(debugLevel_ > 3) {
1022 printMsg("Propagation Update", 1, tm.getElapsedTime(), threadNumber_);
1023 }
1024}
1025
1026template <typename scalarType, typename offsetType>
1027void ttk::ApproximateTopology::buildVertexLinkPolarityApproximate(
1028 const SimplexId vertexId,
1029 std::vector<std::pair<polarity, polarity>> &vlp,
1030 const scalarType *fakeScalars,
1031 const offsetType *const offsets,
1032 const int *const monotonyOffsets) const {
1033
1034 const SimplexId neighborNumber
1035 = multiresTriangulation_.getVertexNeighborNumber(vertexId);
1036 vlp.resize(neighborNumber);
1037
1038 for(SimplexId i = 0; i < neighborNumber; i++) {
1039 SimplexId neighborId0 = 0;
1040 multiresTriangulation_.getVertexNeighbor(vertexId, i, neighborId0);
1041
1042 // const bool lower0 = vertsOrder_[neighborId0] < vertsOrder_[vertexId];
1043 const bool lower0
1044 = ((fakeScalars[neighborId0] < fakeScalars[vertexId])
1045 || (fakeScalars[neighborId0] == fakeScalars[vertexId]
1046 && ((monotonyOffsets[neighborId0] < monotonyOffsets[vertexId])
1047 || (monotonyOffsets[neighborId0] == monotonyOffsets[vertexId]
1048 && offsets[neighborId0] < offsets[vertexId]))));
1049
1050 const polarity isUpper0 = static_cast<polarity>(!lower0) * 255;
1051 vlp[i] = std::make_pair(isUpper0, 0);
1052 }
1053}
1054
1055template <typename scalarType, typename offsetType>
1056bool ttk::ApproximateTopology::printPolarity(
1057 std::vector<polarity> &isNew,
1058 const SimplexId vertexId,
1059 std::vector<std::vector<std::pair<polarity, polarity>>> &vertexLinkPolarity,
1060 const scalarType *scalars,
1061 const scalarType *fakeScalars,
1062 const offsetType *const offsets,
1063 const int *const monotonyOffsets,
1064 bool verbose) {
1065
1066 bool error = false;
1067 std::stringstream mymsg;
1068 std::vector<std::pair<polarity, polarity>> &vlp
1069 = vertexLinkPolarity[vertexId];
1070 mymsg << "POLARITY PRINT"
1071 << "\n";
1072 mymsg << "vertex " << vertexId << " has "
1073 << multiresTriangulation_.getVertexNeighborNumber(vertexId)
1074 << " neighbors"
1075 << "\n";
1076 mymsg << "\tself f:" << fakeScalars[vertexId] << " s:" << scalars[vertexId]
1077 << " o:" << offsets[vertexId] << " m:" << monotonyOffsets[vertexId]
1078 << " isnew: " << (int)isNew[vertexId] << "\n";
1079 if(vlp.empty()) {
1080 if(verbose) {
1081 std::cout << "\tpolarity not initialized for " << vertexId << std::endl;
1082 std::cout << mymsg.str() << std::endl;
1083 }
1084 return false;
1085 }
1086 for(SimplexId i = 0;
1087 i < multiresTriangulation_.getVertexNeighborNumber(vertexId); i++) {
1088 SimplexId nId = -1;
1089 multiresTriangulation_.getVertexNeighbor(vertexId, i, nId);
1090 // find reverse pol
1091 std::vector<std::pair<polarity, polarity>> &vlp2 = vertexLinkPolarity[nId];
1092 bool init = false;
1093 polarity rpol;
1094 if(!vlp2.empty()) {
1095 for(SimplexId j = 0;
1096 j < multiresTriangulation_.getVertexNeighborNumber(nId); j++) {
1097 SimplexId nId2 = -1;
1098 multiresTriangulation_.getVertexNeighbor(nId, j, nId2);
1099 if(nId2 == vertexId) {
1100 rpol = vlp2[j].first;
1101 init = true;
1102 }
1103 }
1104 }
1105
1106 const bool lower
1107 = ((fakeScalars[nId] < fakeScalars[vertexId])
1108 || (fakeScalars[nId] == fakeScalars[vertexId]
1109 && ((monotonyOffsets[nId] < monotonyOffsets[vertexId])
1110 || (monotonyOffsets[nId] == monotonyOffsets[vertexId]
1111 && offsets[nId] < offsets[vertexId]))));
1112
1113 const polarity isUpper = lower ? 0 : 255;
1114
1115 mymsg << " " << i << "th: " << nId << " f:" << fakeScalars[nId]
1116 << " s:" << scalars[nId] << " o:" << offsets[nId]
1117 << " m:" << monotonyOffsets[nId] << " , pol:" << (bool)vlp[i].first
1118 << "(" << (bool)vlp[i].second << ")"
1119 << " rpol:" << (bool)rpol << " true pol:" << (bool)isUpper
1120 << " init " << init << " isnew: " << (int)isNew[nId] << "\n";
1121 if((rpol == isUpper and !vlp2.empty())
1122 or (isUpper != vlp[i].first and !vlp[i].second)) {
1123 mymsg << "POLARITY ERROR "
1124 << "\n";
1125 error = true;
1126 }
1127 }
1128 if(error or verbose) {
1129 std::cout << mymsg.str() << std::endl;
1130 }
1131 return error;
1132}
1133
1134template <typename scalarType, typename offsetType>
1135void ttk::ApproximateTopology::getCriticalTypeApproximate(
1136 const SimplexId &vertexId,
1137 std::vector<std::pair<polarity, polarity>> &vlp,
1138 uint8_t &vertexLink,
1139 DynamicTree &link,
1140 VLBoundaryType &vlbt,
1141 const scalarType *fakeScalars,
1142 const offsetType *const offsets,
1143 const int *const monotonyOffsets) const {
1144
1145 if(vlp.empty()) {
1146 buildVertexLinkPolarityApproximate(
1147 vertexId, vlp, fakeScalars, offsets, monotonyOffsets);
1148 }
1149 const SimplexId neighborNumber
1150 = multiresTriangulation_.getVertexNeighborNumber(vertexId);
1151 link.alloc(neighborNumber);
1152
1153 // associate vertex link boundary
1154 vertexLink = multiresTriangulation_.getVertexBoundaryIndex(vertexId);
1155
1156 // update the link polarity for old points that are processed for
1157 // the first time
1158 const auto &vl = vlbt[vertexLink];
1159
1160 for(size_t edgeId = 0; edgeId < vl.size(); edgeId++) {
1161 const SimplexId n0 = vl[edgeId].first;
1162 const SimplexId n1 = vl[edgeId].second;
1163 if(vlp[n0].first == vlp[n1].first) {
1164 // the smallest id (n0) becomes the parent of n1
1165 link.insertEdge(n1, n0);
1166 }
1167 }
1168}
1169
1170template <typename scalarType, typename offsetType>
1171void ttk::ApproximateTopology::initGlobalPolarity(
1172 std::vector<polarity> &isNew,
1173 std::vector<std::vector<std::pair<polarity, polarity>>> &vertexLinkPolarity,
1174 std::vector<polarity> &toProcess,
1175 const scalarType *fakeScalars,
1176 const offsetType *const offsets,
1177 const int *const monotonyOffsets) const {
1178
1179 Timer timer{};
1180 const size_t nDecVerts = multiresTriangulation_.getDecimatedVertexNumber();
1181
1182 // computes the critical types of all points
1183#ifdef TTK_ENABLE_OPENMP
1184#pragma omp parallel for num_threads(threadNumber_)
1185#endif // TTK_ENABLE_OPENMP
1186 for(size_t i = 0; i < nDecVerts; i++) {
1187 SimplexId const globalId = multiresTriangulation_.localToGlobalVertexId(i);
1188 buildVertexLinkPolarityApproximate(globalId, vertexLinkPolarity[globalId],
1189 fakeScalars, offsets, monotonyOffsets);
1190 toProcess[globalId] = 255;
1191 isNew[globalId] = 0;
1192 }
1193 printMsg("Polarity Init", 1, timer.getElapsedTime(), threadNumber_,
1194 debug::LineMode::NEW, debug::Priority::DETAIL);
1195}
1196
1197template <typename ScalarType, typename offsetType>
1198int ttk::ApproximateTopology::updateGlobalPolarity(
1199 double eps,
1200 std::vector<polarity> &isNew,
1201 std::vector<std::vector<std::pair<polarity, polarity>>> &vertexLinkPolarity,
1202 std::vector<polarity> &toProcess,
1203 std::vector<polarity> &toReprocess,
1204 ScalarType *fakeScalars,
1205 const offsetType *const offsets,
1206 int *monotonyOffsets) const {
1207
1208 const auto nDecVerts = multiresTriangulation_.getDecimatedVertexNumber();
1209
1210#ifdef TTK_ENABLE_OPENMP
1211#pragma omp parallel for num_threads(threadNumber_)
1212#endif // TTK_ENABLE_OPENMP
1213 for(SimplexId localId = 0; localId < nDecVerts; localId++) {
1214 SimplexId const globalId
1215 = multiresTriangulation_.localToGlobalVertexId(localId);
1216 if(!isNew[globalId]) {
1217 getMonotonyChangeByOldPointCPApproximate(
1218 globalId, eps, isNew, toProcess, toReprocess,
1219 vertexLinkPolarity[globalId], fakeScalars, offsets, monotonyOffsets);
1220 }
1221 }
1222
1223 // second Loop process or reprocess
1224#ifdef TTK_ENABLE_OPENMP
1225#pragma omp parallel for num_threads(threadNumber_)
1226#endif
1227 for(int i = 0; i < multiresTriangulation_.getDecimatedVertexNumber(); i++) {
1228 SimplexId const globalId = multiresTriangulation_.localToGlobalVertexId(i);
1229 if(isNew[globalId]) { // new point
1230 if(decimationLevel_ > stoppingDecimationLevel_) {
1231 buildVertexLinkPolarityApproximate(
1232 globalId, vertexLinkPolarity[globalId], fakeScalars, offsets,
1233 monotonyOffsets);
1234 }
1235 isNew[globalId] = 0;
1236
1237 } else { // old point
1238 if(toReprocess[globalId]) {
1239 updateLinkPolarityPonctual(vertexLinkPolarity[globalId]);
1240
1241 toProcess[globalId] = 255; // mark for processing
1242 toReprocess[globalId] = 0;
1243 }
1244 }
1245 } // end for openmp
1246 return 0;
1247}
1248
1249template <typename ScalarType, typename offsetType>
1250void ttk::ApproximateTopology::computeCriticalPoints(
1251 std::vector<std::vector<std::pair<polarity, polarity>>> &vertexLinkPolarity,
1252 std::vector<polarity> &toPropagateMin,
1253 std::vector<polarity> &toPropagateMax,
1254 std::vector<polarity> &toProcess,
1255 std::vector<DynamicTree> &link,
1256 std::vector<uint8_t> &vertexLink,
1257 VLBoundaryType &vertexLinkByBoundaryType,
1258 std::vector<std::vector<SimplexId>> &saddleCCMin,
1259 std::vector<std::vector<SimplexId>> &saddleCCMax,
1260 ScalarType *fakeScalars,
1261 const offsetType *const offsets,
1262 int *monotonyOffsets) {
1263
1264 Timer tm;
1265
1266#ifdef TTK_ENABLE_OPENMP
1267#pragma omp parallel for num_threads(threadNumber_)
1268#endif
1269 for(int i = 0; i < multiresTriangulation_.getDecimatedVertexNumber(); i++) {
1270 SimplexId const globalId = multiresTriangulation_.localToGlobalVertexId(i);
1271
1272 if(toProcess[globalId]) {
1273 // nbComputations++;
1274 getCriticalTypeApproximate(globalId, vertexLinkPolarity[globalId],
1275 vertexLink[globalId], link[globalId],
1276 vertexLinkByBoundaryType, fakeScalars, offsets,
1277 monotonyOffsets);
1278 getValencesFromLink(globalId, vertexLinkPolarity[globalId],
1279 link[globalId], toPropagateMin, toPropagateMax,
1280 saddleCCMin, saddleCCMax);
1281 }
1282 } // end for openmp
1283
1284 if(debugLevel_ > 3) {
1285 printMsg(
1286 "Critical Points Computation", 1, tm.getElapsedTime(), threadNumber_);
1287 }
1288}
1289
1290template <typename scalarType>
1291int ttk::ApproximateTopology::sortPersistenceDiagramApproximate(
1292 std::vector<PersistencePair> &diagram,
1293 scalarType *fakeScalars,
1294 const SimplexId *const offsets,
1295 int *monotonyOffsets) const {
1296 auto cmp = [fakeScalars, offsets, monotonyOffsets](
1297 const PersistencePair &pA, const PersistencePair &pB) {
1298 const ttk::SimplexId a = pA.birth;
1299 const ttk::SimplexId b = pB.birth;
1300
1301 return ((fakeScalars[a] < fakeScalars[b])
1302 || (fakeScalars[a] == fakeScalars[b]
1303 && ((monotonyOffsets[a] < monotonyOffsets[b])
1304 || (monotonyOffsets[a] == monotonyOffsets[b]
1305 && offsets[a] < offsets[b]))));
1306 };
1307
1308 std::sort(diagram.begin(), diagram.end(), cmp);
1309
1310 return 0;
1311}
1312
1313template <typename scalarType>
1314int ttk::ApproximateTopology::computeApproximatePD(
1315 std::vector<PersistencePair> &CTDiagram,
1316 const scalarType *scalars,
1317 scalarType *const fakeScalars,
1318 SimplexId *const outputOffsets,
1319 int *const outputMonotonyOffsets) {
1320
1321 int ret = -1;
1322 std::stringstream ss;
1323 ss << "Approximate Persistence Diagram computation with "
1324 << debug::output::UNDERLINED << debug::output::YELLOW << epsilon_ * 100
1325 << "%" << debug::output::ENDCOLOR << debug::output::ENDCOLOR << " error";
1326 printMsg(ss.str());
1327
1328 ret = executeApproximateTopology(
1329 scalars, fakeScalars, outputOffsets, outputMonotonyOffsets);
1330 CTDiagram = std::move(CTDiagram_);
1331 CTDiagram_.clear();
1332
1333 return ret;
1334}
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:1027
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:667
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:1171
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:1291
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:1314
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:921
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:1056
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:1135
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:788
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::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:1250
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:1198
ttk::Debug::setDebugMsgPrefix
void setDebugMsgPrefix(const std::string &prefix)
Definition Debug.h:364
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::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)
Generated on Thu Sep 26 2024 06:50:17 for TTK by doxygen 1.9.8