TTK
Loading...
Searching...
No Matches
FTRGraphPrivate_Template.h
Go to the documentation of this file.
1#pragma once
2
3// local includes
4#include "FTRGraph.h"
5#include "FTRGraph_Template.h"
6#include "FTRTasks.h"
7
8// c++ includes
9#include <unordered_map>
10
11// trick to print a full line in an atomic operation, avoid mixed up redulsts in
12// parallel
13#ifndef NDEBUG
14// #define PRINT(msg) {std::stringstream s; s << msg << std::endl; std::cout <<
15// s.str();}
16#define PRINT(msg)
17#else
18#define PRINT(msg)
19#endif
20
21template <typename ScalarType, typename triangulationType>
23 const idVertex seed, Propagation *localProp, idSuperArc currentArc) {
24
25 TTK_FORCE_USE(seed);
26 PRINT(seed << " :: " << localProp->getNbArcs());
27 if(!localProp->getNbArcs()) {
28 PRINT(seed << " Stop direct");
29#ifdef TTK_ENABLE_FTR_TASK_STATS
30 {
31 // This only when tasks grows exclusively from min
32 // This propagation is dying here
33 idVertex curProp;
34 float curTime;
35#pragma omp atomic capture seq_cst
36 {
37 curProp = nbProp_;
38 --nbProp_;
39 }
40#pragma omp critical(stats)
41 { curTime = sweepStart_.getElapsedTime(); }
42 propTimes_[curProp - 1] = curTime;
43 }
44#endif
45 return;
46 }
47
48 // topology
49 bool isJoinLast = false;
50 bool isJoin = false, isSplit = false;
51
52 // containers
53 // vitsit
54 Star star;
55 Comp comp;
56
57 while(!isJoin && !isSplit && !localProp->empty()) {
58 localProp->nextVertex();
59 const idVertex curVert = localProp->getCurVertex();
60 idSuperArc mergeIn = nullSuperArc;
61
62 if(localProp->getNbArcs() == 0) {
63#ifdef TTK_ENABLE_FTR_TASK_STATS
64 {
65 // This only when tasks grows exclusively from min
66 // This propagation is dying here
67 idVertex curProp;
68 float curTime;
69#pragma omp atomic capture seq_cst
70 {
71 curProp = nbProp_;
72 --nbProp_;
73 }
74#pragma omp critical(stats)
75 { curTime = sweepStart_.getElapsedTime(); }
76 propTimes_[curProp - 1] = curTime;
77 }
78#endif
79 return;
80 }
81
82 PRINT("<" << curVert << " " << localProp->goUp() << " a "
83 << localProp->getNbArcs());
84
85#ifdef TTK_ENABLE_FTR_VERT_STATS
86 graph_.incTouch(curVert);
87 graph_.setNbArcActive(curVert, localProp->getNbArcs());
88#endif
89
90 visitStar(localProp, star);
91
92 if(propagations_.hasVisitedOpposite(
93 curVert, localProp) /*&& !graph_.isNode(curVert)*/) {
94 bool ignoreVert = false;
95 for(auto edge : star.lower) {
96 const idSuperArc tmpLowArc = dynGraph(localProp).getSubtreeArc(edge);
97 if(tmpLowArc != nullSuperArc && !graph_.getArc(tmpLowArc).isVisible()
98 && graph_.getArc(tmpLowArc).getPropagation() == localProp) {
99 PRINT("-" << curVert << " " << graph_.printArc(tmpLowArc));
100#ifdef TTK_ENABLE_FTR_VERT_STATS
101 graph_.incAvoid();
102#endif
103 ignoreVert = true;
104 break;
105 } else if(tmpLowArc != nullSuperArc) {
106 PRINT("!" << graph_.printArc(tmpLowArc));
107 }
108 }
109 if(ignoreVert)
110 continue;
111 }
112
113#ifndef TTK_DISABLE_FTR_LAZY
114 if(valences_.lower[curVert] < 2 && valences_.upper[curVert] < 2) {
115
116 // not a local min (for local min, currentArc is already set)
117 if(star.lower.size()) {
118 // not a min nor a saddle: 1 CC below
119 currentArc = dynGraph(localProp).getSubtreeArc(star.lower[0]);
120 if(currentArc == nullSuperArc) {
121 PRINT("n-" << curVert);
122 continue;
123 }
124 if(valences_.upper[curVert] && valences_.lower[curVert]) {
125 // not saddle neither extrema
126 graph_.getArc(currentArc).visit(curVert);
127 }
128 }
129
130 mergeIn = visit(localProp, currentArc);
131
132 lazyUpdatePreimage(localProp, currentArc);
133
134 // ensure we will always recover this arc from the upper neighbors
135 for(const idEdge dgNode : star.upper) {
136 dynGraph(localProp).setCorArc(dgNode, currentArc);
137 }
138 } else {
139 // locally apply the lazy one the current growing arc
140 for(const idEdge e : star.lower) {
141 const auto node{dynGraph(localProp).getNode(e)};
142 const idSuperArc a
143 = node != nullptr ? node->findRootArc() : nullSuperArc;
144 if(a != nullSuperArc && graph_.getArc(a).isVisible()
145 && graph_.getArc(a).getPropagation()->getId()
146 == localProp->getId()) {
147 lazyApply(localProp, a);
148 }
149 }
150#else
151 {
152#endif
153 comp.lower = lowerComps(star.lower, localProp);
154
155 if(comp.lower.size() > 1) {
156 isJoin = true;
157 isJoinLast = checkLast(localProp, star.lower);
158 break;
159 } else {
160 if(comp.lower.size()) {
161 currentArc = (*comp.lower.begin())->getCorArc();
162 if(currentArc == nullSuperArc) {
163 PRINT("n--" << curVert);
164 continue;
165 }
166 }
167 mergeIn = visit(localProp, currentArc);
168 }
169 updatePreimage(localProp, currentArc);
170 comp.upper = upperComps(star.upper, localProp);
171 if(comp.upper.size() > 1) {
172 isSplit = true;
173 }
174
175 if(!isJoin && !isSplit && comp.upper.size()) {
176 // this arc is not empty (not saddle not max)
177 graph_.getArc(currentArc).visit(curVert);
178 }
179 }
180
181 // merge current arc
182 if(mergeIn != nullSuperArc && graph_.getArc(currentArc).isVisible()) {
183 graph_.getArc(currentArc).merge(mergeIn);
184 localProp->lessArc();
185 PRINT(">" << graph_.printArc(currentArc) << " "
186 << graph_.printArc(mergeIn));
187 // Can't stop here
188 // try interleaved launch sequential on dragon.vtu
189 // if (localProp->getNbArcs() == 0) return;
190 }
191
192 // stop on leaves
193 if(!star.upper.size()) {
194 // We have reached a local extrema (max from this propagation)
195 const idNode leafNode = graph_.makeNode(curVert);
196 graph_.closeArc(currentArc, leafNode);
197 PRINT("^" << graph_.printArc(currentArc));
198 if(graph_.getArc(currentArc).isVisible()) {
199 // do not decrease on merged arcs
200 localProp->lessArc();
201 }
202 if(localProp->getNbArcs() == 0) {
203 // no more active arcs
204 PRINT(curVert << " stop");
205#ifdef TTK_ENABLE_FTR_TASK_STATS
206 {
207 // This only when tasks grows exclusively from min
208 // This propagation is dying here
209 idVertex curProp;
210 float curTime;
211#pragma omp atomic capture seq_cst
212 {
213 curProp = nbProp_;
214 --nbProp_;
215 }
216#pragma omp critical(stats)
217 { curTime = sweepStart_.getElapsedTime(); }
218 propTimes_[curProp - 1] = curTime;
219 }
220#endif
221 return;
222 }
223 }
224
225 // add upper star for futur visit
226 localGrowth(localProp, star.upper);
227 } // end propagation while
228
229 // get the corresponding critical point on which
230 // the propagation has stopped (join, split, max)
231 const idVertex upVert = localProp->getCurVertex();
232
233 if(localProp->getNbArcs() == 0 || localProp->empty()) {
234 PRINT(upVert << " stop " << localProp->getNbArcs());
235#ifdef TTK_ENABLE_FTR_TASK_STATS
236 {
237 // This only when tasks grows exclusively from min
238 // This propagation is dying here
239 idVertex curProp;
240 float curTime;
241#pragma omp atomic capture seq_cst
242 {
243 curProp = nbProp_;
244 --nbProp_;
245 }
246#pragma omp critical(stats)
247 { curTime = sweepStart_.getElapsedTime(); }
248 propTimes_[curProp - 1] = curTime;
249 }
250#endif
251 return;
252 }
253
254 PRINT(upVert << " active " << localProp->getNbArcs());
255
256 // reached node id and whether it has been created by this task or already
257 // existed
258 idNode saddleNode;
259 idSuperArc joinParentArc{};
260 bool hideFromHere = false; // if true, new arc are hidden to stop propagation.
261
262#ifdef TTK_ENABLE_OPENMP
263#pragma omp critical
264#endif
265 {
266 bool alreadyNode = graph_.isNode(upVert);
267 hideFromHere = alreadyNode;
268 if(isJoin) {
269 PRINT(upVert << " join " << isJoinLast << " h " << hideFromHere);
270 if(isJoinLast) {
271 graph_.makeNode(upVert);
272 }
273 }
274 if(isSplit) {
275 PRINT(upVert << " split " << hideFromHere);
276 graph_.makeNode(upVert);
277 }
278 }
279
280 // arc management
281 if(isJoinLast) {
282 // ensure we have the good values here, even if other tasks were doing
283 // stuff
284 {
285 // here to solve a 1 over thousands execution bug in parallel
286 // TODO Still required ??
287 visitStar(localProp, star);
288 comp.lower = lowerComps(star.lower, localProp);
289 }
290 saddleNode = graph_.getNodeId(upVert);
291 idSuperArc visibleMerged = mergeAtSaddle(saddleNode, localProp, comp.lower);
292 localProp->lessArc(visibleMerged - 1);
293
294 localGrowth(localProp, star.upper);
295
296 joinParentArc = graph_.openArc(saddleNode, localProp);
297 visit(localProp, joinParentArc);
298 updatePreimage(localProp, joinParentArc);
299 comp.upper = upperComps(star.upper, localProp);
300
301 // do not propagate
302 if(hideFromHere) {
303 if(graph_.getArc(joinParentArc).hide()) {
304 localProp->lessArc();
305 }
306 }
307
308 // split detection required after the merge
309 if(comp.upper.size() > 1) {
310 // this node is both join and split
311 isSplit = true;
312 // will be replaced be new arcs of the split
313 if(graph_.getArc(joinParentArc).hide()) {
314 localProp->lessArc();
315 }
316 }
317
318 PRINT("+" << graph_.printArc(joinParentArc));
319 } else if(!isJoin) {
320 // only one arc coming here
321 saddleNode = graph_.getNodeId(upVert);
322 graph_.closeArc(currentArc, saddleNode);
323 if(localProp->getNbArcs() == 0) {
324 // no more active arcs
325 PRINT(upVert << " stop");
326#ifdef TTK_ENABLE_FTR_TASK_STATS
327 {
328 // This only when tasks grows exclusively from min
329 // This propagation is dying here
330 idVertex curProp;
331 float curTime;
332#pragma omp atomic capture seq_cst
333 {
334 curProp = nbProp_;
335 --nbProp_;
336 }
337#pragma omp critical(stats)
338 { curTime = sweepStart_.getElapsedTime(); }
339 propTimes_[curProp - 1] = curTime;
340 }
341#endif
342 return;
343 }
344 if(graph_.getArc(currentArc).isVisible()) {
345 // if not visible, already decremented the counter
346 localProp->lessArc();
347 }
348 PRINT("|" << graph_.printArc(currentArc));
349 }
350
351 if(isSplit) {
352 splitAtSaddle(localProp, comp.upper, hideFromHere);
353 if(!hideFromHere) {
354 localProp->moreArc(comp.upper.size());
355 }
356 }
357
358 // starting from the saddle
359 if(isSplit && (!isJoin || isJoinLast)) {
360#ifdef TTK_ENABLE_OPENMP
361#pragma omp task OPTIONAL_PRIORITY(PriorityLevel::Low)
362#endif
363 growthFromSeed(upVert, localProp);
364
365 } else if(isJoinLast) {
366
367#ifdef TTK_ENABLE_OPENMP
368#pragma omp task OPTIONAL_PRIORITY(PriorityLevel::Average)
369#endif
370 growthFromSeed(upVert, localProp, joinParentArc);
371 }
372#ifdef TTK_ENABLE_FTR_TASK_STATS
373 else {
374 // This only when tasks grows exclusively from min
375 // This propagation is dying here
376 idVertex curProp;
377 float curTime;
378#pragma omp atomic capture seq_cst
379 {
380 curProp = nbProp_;
381 --nbProp_;
382 }
383#pragma omp critical(stats)
384 { curTime = sweepStart_.getElapsedTime(); }
385 propTimes_[curProp - 1] = curTime;
386 }
387#endif
388}
389
390template <typename ScalarType, typename triangulationType>
392 const idVertex begin, const idVertex stop) {
393 Star star;
394 Comp comp;
395
396 const bool fromMin = begin < stop;
397 Propagation *localProp
398 = newPropagation(scalars_.getSortedVert(begin), fromMin);
399 const idVertex incr = fromMin ? 1 : -1;
400 for(idVertex idv = begin; idv < stop; idv = idv + incr) {
401 const idVertex curVert = scalars_.getSortedVert(idv);
402 localProp->setCurvert(curVert); // never use Fibo Heap here
403
404 idSuperArc currentArc = nullSuperArc;
405 visitStar(localProp, star);
406
407 if(valences_.lower[curVert] == 0) { // min
408 const idNode minNode = graph_.makeNode(curVert);
409 currentArc = graph_.openArc(minNode, localProp);
410 }
411
412#ifndef TTK_DISABLE_FTR_LAZY
413 if(valences_.lower[curVert] < 2 && valences_.upper[curVert] < 2) {
414
415 // simple reeb regular, laziness
416 if(!star.lower.empty()) {
417 // not a min nor a saddle: 1 CC below (need findSubtree)
418 currentArc = dynGraph(localProp).getSubtreeArc(star.lower[0]);
419 if(valences_.upper[curVert] && valences_.lower[curVert]) {
420 // not saddle neither extrema
421 graph_.getArc(currentArc).visit(curVert);
422 }
423 }
424 if(star.upper.size()) {
425 for(const idEdge dgNode : star.upper) {
426 dynGraph(localProp).setCorArc(dgNode, currentArc);
427 }
428 } else { // max
429 const idNode maxNode = graph_.makeNode(curVert);
430 graph_.closeArc(currentArc, maxNode);
431 }
432
433 visit(localProp, currentArc);
434
435 lazyUpdatePreimage(localProp, currentArc);
436 } else {
437
438 // locally aply the lazy one the current growing arc
439 for(const idEdge e : star.lower) {
440 const idSuperArc a = dynGraph(localProp).getNode(e)->findRootArc();
441 if(!lazy_.isEmpty(a)) {
442 // process lazy
443 // sort both list
444 // for each:
445 // if del < add: delete in real RG
446 // if add < del: add in real RG
447 // else: drop both, computation avoided
448 lazyApply(localProp, a);
449 }
450 }
451#else
452 {
453#endif
454 bool isJoin = false;
455 comp.lower = lowerComps(star.lower, localProp);
456 if(comp.lower.size() == 1) { // regular
457 currentArc = (*comp.lower.begin())->getCorArc();
458 } else if(comp.lower.size() > 1) { // join saddle
459 const idNode sadNode = graph_.makeNode(curVert);
460 currentArc = graph_.openArc(sadNode, localProp);
461 mergeAtSaddle(sadNode, comp.lower);
462 isJoin = true;
463 }
464
465 graph_.visit(curVert, currentArc);
466 propagations_.visit(curVert, localProp);
467 updatePreimage(localProp, currentArc);
468
469 comp.upper = upperComps(star.upper, localProp);
470 if(!comp.upper.size()) { // max
471 const idNode maxNode = graph_.makeNode(curVert);
472 graph_.closeArc(currentArc, maxNode);
473 } else if(comp.upper.size() < 2) {
474 if(!isJoin) {
475 // this arc is not empty
476 graph_.getArc(currentArc).visit(curVert);
477 }
478 } else {
479 if(isJoin) {
480 graph_.getArc(currentArc).hide();
481 } else { // split
482 const idNode splitNode = graph_.makeNode(curVert);
483 graph_.closeArc(currentArc, splitNode);
484 }
485 splitAtSaddle(localProp, comp.upper);
486 }
487 }
488 } // end for each vertex
489}
490
491template <typename ScalarType, typename triangulationType>
493 const Propagation *const localProp, Star &star) const {
494
495 star.lower.clear();
496 star.upper.clear();
497
498 const idEdge nbAdjEdges
499 = mesh_.getVertexEdgeNumber(localProp->getCurVertex());
500 star.lower.reserve(nbAdjEdges);
501 star.upper.reserve(nbAdjEdges);
502
503 for(idEdge e = 0; e < nbAdjEdges; ++e) {
504 idEdge edgeId;
505 mesh_.getVertexEdge(localProp->getCurVertex(), e, edgeId);
506 idVertex edgeLowerVert, edgeUpperVert;
507 std::tie(edgeLowerVert, edgeUpperVert)
508 = mesh_.getOrderedEdge(edgeId, localProp->goUp());
509 if(edgeLowerVert == localProp->getCurVertex()) {
510 star.upper.emplace_back(edgeId);
511 } else {
512 star.lower.emplace_back(edgeId);
513 }
514 }
515}
516
517template <typename ScalarType, typename triangulationType>
518std::set<ttk::ftr::DynGraphNode<ttk::ftr::idVertex> *>
520 const std::vector<idEdge> &finishingEdges,
521 const Propagation *const localProp) {
522 return dynGraph(localProp).findRoot(finishingEdges);
523}
524
525template <typename ScalarType, typename triangulationType>
526std::set<ttk::ftr::DynGraphNode<ttk::ftr::idVertex> *>
528 const std::vector<idEdge> &startingEdges,
529 const Propagation *const localProp) {
530 return dynGraph(localProp).findRoot(startingEdges);
531}
532
533template <typename ScalarType, typename triangulationType>
535 const std::vector<DynGraphNode<idVertex> *> &compVect) {
536 for(const auto *dgNode : compVect) {
537 const idSuperArc arc = dgNode->getCorArc();
538 if(arc != nullSuperArc && !graph_.getArc(arc).isVisible()) {
539 return true;
540 }
541 }
542 return false;
543}
544
545template <typename ScalarType, typename triangulationType>
546std::pair<ttk::ftr::valence, ttk::ftr::valence>
548 const idVertex curVert,
549 LocalForests &localForests,
550 const VertCompFN &comp) {
551 // traduce edge id in a local id for the forests
552 std::unordered_map<idEdge, std::size_t> mapNeighDown, mapNeighUp;
553 std::size_t nextId = 0;
554
555 localForests.up.reset();
556 localForests.down.reset();
557 const idVertex oldUpCC = localForests.up.getNbCC();
558 const idVertex oldDownCC = localForests.down.getNbCC();
559 const idCell nbTri = mesh_.getVertexTriangleNumber(curVert);
560 for(idCell t = 0; t < nbTri; ++t) {
561 idCell curTri;
562 mesh_.getVertexTriangle(curVert, t, curTri);
563 // edges found, max 2 per triangles.
564 std::size_t downSide[2], upSide[2];
565 unsigned char curDownSide = 0, curUpSide = 0;
566 for(idVertex v = 0; v < 3; ++v) {
567 idEdge triEdge;
568 mesh_.getTriangleEdge(curTri, v, triEdge);
569 idVertex v0, v1;
570 mesh_.getEdgeVertex(triEdge, 0, v0);
571 mesh_.getEdgeVertex(triEdge, 1, v1);
572
573 // only consider the edges in the star
574 if(v0 == curVert) {
575 if(comp(curVert, v1)) {
576 if(mapNeighUp.count(triEdge)) {
577 upSide[curUpSide++] = mapNeighUp[triEdge];
578 } else {
579 upSide[curUpSide++] = nextId;
580 mapNeighUp[triEdge] = nextId++;
581 }
582 } else {
583 if(mapNeighDown.count(triEdge)) {
584 downSide[curDownSide++] = mapNeighDown[triEdge];
585 } else {
586 downSide[curDownSide++] = nextId;
587 mapNeighDown[triEdge] = nextId++;
588 }
589 }
590 } else if(v1 == curVert) {
591 if(comp(curVert, v0)) {
592 if(mapNeighUp.count(triEdge)) {
593 upSide[curUpSide++] = mapNeighUp[triEdge];
594 } else {
595 upSide[curUpSide++] = nextId;
596 mapNeighUp[triEdge] = nextId++;
597 }
598 } else {
599 if(mapNeighDown.count(triEdge)) {
600 downSide[curDownSide++] = mapNeighDown[triEdge];
601 } else {
602 downSide[curDownSide++] = nextId;
603 mapNeighDown[triEdge] = nextId++;
604 }
605 }
606 }
607 }
608
609 // if both edges of the triangle were up or down, we add a link btwn
610 // them This is how the number of component is reduced
611 if(curDownSide == 2) {
612 localForests.down.insertEdge(downSide[0], downSide[1], 0, nullSuperArc);
613 } else if(curUpSide == 2) {
614 localForests.up.insertEdge(upSide[0], upSide[1], 0, nullSuperArc);
615 }
616 }
617
618 const valence down
619 = mapNeighDown.size() - (oldDownCC - localForests.down.getNbCC());
620 const valence up = mapNeighUp.size() - (oldUpCC - localForests.up.getNbCC());
621 return {down, up};
622}
623
624template <typename ScalarType, typename triangulationType>
626 const Propagation *const localProp, const idSuperArc curArc) {
627 const idCell nbAdjTriangles
628 = mesh_.getVertexTriangleNumber(localProp->getCurVertex());
629
630 orderedTriangle oTriangle;
631
632 // TODO SORT TRANGLES HERE ?
633
634 for(idCell t = 0; t < nbAdjTriangles; ++t) {
635 // Classify current cell
636 idCell curTriangleid;
637 mesh_.getVertexTriangle(localProp->getCurVertex(), t, curTriangleid);
638
639 mesh_.getOrderedTriangle(curTriangleid, localProp->goUp(), oTriangle);
640 vertPosInTriangle curVertPos = getVertPosInTriangle(oTriangle, localProp);
641
642 // Update DynGraph
643 // We can have an end pos on an unvisited triangle
644 // in case of saddle points
645 switch(curVertPos) {
646 case vertPosInTriangle::Start:
647 updatePreimageStartCell(oTriangle, localProp, curArc);
648 break;
649 case vertPosInTriangle::Middle:
650 updatePreimageMiddleCell(oTriangle, localProp, curArc);
651 break;
652 case vertPosInTriangle::End:
653 updatePreimageEndCell(oTriangle, localProp, curArc);
654 break;
655 default:
656 std::cout << "[FTR]: update preimage error, unknown vertPos type"
657 << std::endl;
658 break;
659 }
660 }
661}
662
663template <typename ScalarType, typename triangulationType>
665 const orderedTriangle &oTriangle,
666 const Propagation *const localProp,
667 const idSuperArc curArc) {
668 const orderedEdge e0
669 = mesh_.getOrderedEdge(std::get<0>(oTriangle), localProp->goUp());
670 const orderedEdge e1
671 = mesh_.getOrderedEdge(std::get<1>(oTriangle), localProp->goUp());
672 const idVertex w = getWeight(e0, e1, localProp);
673
674 // this order for history
675 dynGraph(localProp).insertEdge(
676 std::get<1>(oTriangle), std::get<0>(oTriangle), w, curArc);
677}
678
679template <typename ScalarType, typename triangulationType>
681 updatePreimageMiddleCell(const orderedTriangle &oTriangle,
682 const Propagation *const localProp,
683 const idSuperArc curArc) {
684 // Check if exist ?
685 // If not, the triangle will be visited again once a merge have occurred.
686 // So we do not add the edge now
687 dynGraph(localProp).removeEdge(
688 std::get<0>(oTriangle), std::get<1>(oTriangle));
689
690 // keep history inside the dyngraph structure
691 dynGraph(localProp).setCorArc(std::get<0>(oTriangle), curArc);
692 // dynGraph(localProp).setCorArc(std::get<1>(oTriangle), curArc);
693
694 const orderedEdge e1
695 = mesh_.getOrderedEdge(std::get<1>(oTriangle), localProp->goUp());
696 const orderedEdge e2
697 = mesh_.getOrderedEdge(std::get<2>(oTriangle), localProp->goUp());
698 const idVertex w = getWeight(e1, e2, localProp);
699
700 // this order for history
701 dynGraph(localProp).insertEdge(
702 std::get<1>(oTriangle), std::get<2>(oTriangle), w, curArc);
703}
704
705template <typename ScalarType, typename triangulationType>
707 const orderedTriangle &oTriangle,
708 const Propagation *const localProp,
709 const idSuperArc curArc) {
710 dynGraph(localProp).removeEdge(
711 std::get<1>(oTriangle), std::get<2>(oTriangle));
712 dynGraph(localProp).setCorArc(std::get<1>(oTriangle), curArc);
713 // dynGraph(localProp).setCorArc(std::get<2>(oTriangle), curArc);
714}
715
716template <typename ScalarType, typename triangulationType>
718 const idVertex seed,
719 const idEdge neigEdge,
720 const idSuperArc curArc,
721 const Propagation *const localProp) {
722 idVertex v0;
723 idVertex v1;
724 mesh_.getEdgeVertex(neigEdge, 0, v0);
725 mesh_.getEdgeVertex(neigEdge, 1, v1);
726
727 const idVertex other = (v0 == seed) ? v1 : v0;
728
729 if(localProp->compare(seed, other)) {
730 dynGraph(localProp).setSubtreeArc(neigEdge, curArc);
731 }
732}
733
734#ifndef TTK_DISABLE_FTR_LAZY
735template <typename ScalarType, typename triangulationType>
737 Propagation *const localProp, const idSuperArc curArc) {
738 const idCell nbAdjTriangles
739 = mesh_.getVertexTriangleNumber(localProp->getCurVertex());
740
741 orderedTriangle oTriangle;
742
743 for(idCell t = 0; t < nbAdjTriangles; ++t) {
744 // Classify current cell
745 idCell curTriangleid;
746 mesh_.getVertexTriangle(localProp->getCurVertex(), t, curTriangleid);
747
748 mesh_.getOrderedTriangle(curTriangleid, localProp->goUp(), oTriangle);
749 vertPosInTriangle curVertPos = getVertPosInTriangle(oTriangle, localProp);
750
751 // Update DynGraph
752 // We can have an end pos on an unvisited triangle
753 // in case of saddle points
754 switch(curVertPos) {
755 case vertPosInTriangle::Start:
756 updateLazyStart(oTriangle, localProp, curArc);
757 break;
758 case vertPosInTriangle::Middle:
759 updateLazyMiddle(oTriangle, localProp, curArc);
760 break;
761 case vertPosInTriangle::End:
762 updateLazyEnd(oTriangle, localProp, curArc);
763 break;
764 default:
765 std::cout << "[FTR]: lazy update preimage error, unknown vertPos type"
766 << std::endl;
767 break;
768 }
769 }
770}
771
772template <typename ScalarType, typename triangulationType>
774 const orderedTriangle &oTriangle,
775 Propagation *const ttkNotUsed(localProp),
776 const idSuperArc curArc) {
777 lazy_.addEmplace(std::get<0>(oTriangle), std::get<1>(oTriangle), curArc);
778}
779
780template <typename ScalarType, typename triangulationType>
782 const orderedTriangle &oTriangle,
783 Propagation *const localProp,
784 const idSuperArc curArc) {
785 lazy_.delEmplace(std::get<0>(oTriangle), std::get<1>(oTriangle), curArc);
786 dynGraph(localProp).removeEdge(
787 std::get<0>(oTriangle), std::get<1>(oTriangle));
788 dynGraph(localProp).setCorArc(std::get<0>(oTriangle), curArc);
789 dynGraph(localProp).setCorArc(std::get<1>(oTriangle), curArc);
790 lazy_.addEmplace(std::get<1>(oTriangle), std::get<2>(oTriangle), curArc);
791}
792
793template <typename ScalarType, typename triangulationType>
795 const orderedTriangle &oTriangle,
796 Propagation *const localProp,
797 const idSuperArc curArc) {
798 lazy_.delEmplace(std::get<1>(oTriangle), std::get<2>(oTriangle), curArc);
799 dynGraph(localProp).removeEdge(
800 std::get<1>(oTriangle), std::get<2>(oTriangle));
801 dynGraph(localProp).setCorArc(std::get<1>(oTriangle), curArc);
802 dynGraph(localProp).setCorArc(std::get<2>(oTriangle), curArc);
803}
804
805template <typename ScalarType, typename triangulationType>
807 const Propagation *const localProp,
808 const linkEdge edge,
809 const idSuperArc arc) {
810 const orderedEdge e0
811 = mesh_.getOrderedEdge(std::get<0>(edge), localProp->goUp());
812 const orderedEdge e1
813 = mesh_.getOrderedEdge(std::get<1>(edge), localProp->goUp());
814 const idVertex w = getWeight(e0, e1, localProp);
815 dynGraph(localProp).insertEdge(std::get<1>(edge), std::get<0>(edge), w, arc);
816}
817
818template <typename ScalarType, typename triangulationType>
820 const Propagation *const localProp,
821 const linkEdge edge,
822 const idSuperArc ttkNotUsed(arc)) {
823 dynGraph(localProp).removeEdge(std::get<0>(edge), std::get<1>(edge));
824}
825
826template <typename ScalarType, typename triangulationType>
828 Propagation *const localProp, const idSuperArc a) {
829 auto add = lazy_.addGetNext(a);
830 while(add != nullLink) {
831 updateLazyAdd(localProp, add, a);
832 add = lazy_.addGetNext(a);
833 }
834}
835
836#endif
837
838template <typename ScalarType, typename triangulationType>
840 const idVertex seed,
841 const idSuperArc curArc,
842 const Propagation *const localProp) {
843 const idVertex nbEdgesNeigh = mesh_.getVertexEdgeNumber(seed);
844 for(idVertex nid = 0; nid < nbEdgesNeigh; ++nid) {
845 idEdge edgeId;
846 mesh_.getVertexEdge(seed, nid, edgeId);
847
848 updateDynGraphCurArc(seed, edgeId, curArc, localProp);
849 }
850}
851
852template <typename ScalarType, typename triangulationType>
854 Propagation *const localProp, const std::vector<idEdge> &upperEdges) {
855 const idVertex curVert = localProp->getCurVertex();
856 for(const idEdge e : upperEdges) {
857 idVertex v0, v1;
858 mesh_.getEdgeVertex(e, 0, v0);
859 mesh_.getEdgeVertex(e, 1, v1);
860 const idVertex other = (v0 == curVert) ? v1 : v0;
861 if(localProp->compare(curVert, other)) {
862 if(!propagations_.willVisit(other, localProp)) {
863 localProp->addNewVertex(other);
864 propagations_.toVisit(other, localProp);
865 }
866 }
867 }
868}
869
870template <typename ScalarType, typename triangulationType>
872 Propagation *const localProp, const std::vector<idEdge> &starVect) {
873 const idVertex curSaddle = localProp->getCurVertex();
874 AtomicUF *curId = localProp->getId();
875 valence decr = 0;
876
877 // NOTE:
878 // Using propagation id allows to decrement by the number of time this
879 // propagation has reached the saddle, even if the propagation take care
880 // of several of these arcs (after a Hole-split).
881 for(idEdge edgeId : starVect) {
882 const idSuperArc edgeArc = dynGraph(localProp).getSubtreeArc(edgeId);
883 if(edgeArc == nullSuperArc) {
884 continue;
885 }
886 AtomicUF *tmpId = graph_.getArc(edgeArc).getPropagation()->getId();
887 if(tmpId == curId) {
888 graph_.getArc(edgeArc).setEnd(curSaddle);
889 ++decr;
890 }
891 }
892
893#if defined(__GNUC__)
894#pragma GCC diagnostic push
895#pragma GCC diagnostic ignored "-Wunused-value"
896#endif // __GNUC__
897
898 valence oldVal = 0;
899 if(localProp->goUp()) {
900 // for gcc 4.8 and old openMP
901 valence *const vd = &graph_.valDown_[curSaddle];
902#ifdef TTK_ENABLE_OPENMP
903#pragma omp atomic capture
904#endif
905 {
906 oldVal = *vd;
907 *vd -= decr;
908 }
909
910 } else {
911 valence *const vu = &graph_.valUp_[curSaddle];
912#ifdef TTK_ENABLE_OPENMP
913#pragma omp atomic capture
914#endif
915 {
916 oldVal = *vu;
917 *vu -= decr;
918 }
919 }
920
921 if(oldVal == -1) {
922 // First task to touch this saddle, compute the valence
923 idVertex totalVal = starVect.size();
924 valence newVal = 0;
925 if(localProp->goUp()) {
926 valence *const vd = &graph_.valDown_[curSaddle];
927#ifdef TTK_ENABLE_OPENMP
928#pragma omp atomic capture
929#endif
930 {
931 newVal = *vd;
932 *vd += (totalVal + 1);
933 }
934
935 } else {
936 valence *const vu = &graph_.valUp_[curSaddle];
937#ifdef TTK_ENABLE_OPENMP
938#pragma omp atomic capture
939#endif
940 {
941 newVal = *vu;
942 *vu += (totalVal + 1);
943 }
944 }
945 oldVal = decr + newVal + (totalVal + 1);
946 }
947
948#if defined(__GNUC__)
949#pragma GCC diagnostic pop
950#endif // __GNUC__
951
952 return oldVal == decr;
953}
954
955template <typename ScalarType, typename triangulationType>
958 const idNode saddleId,
959 Propagation *localProp,
960 const std::set<DynGraphNode<idVertex> *> &compVect) {
961
962#ifndef TTK_ENABLE_KAMIKAZE
963 if(compVect.size() < 2) {
964 std::cerr << "[FTR]: merge at saddle with only one lower CC" << std::endl;
965 }
966#endif
967
968 idSuperArc visibleClosed = 0;
969 for(auto *dgNode : compVect) {
970 // read in the history (lower comp already contains roots)
971 const idSuperArc endingArc = dgNode->getCorArc();
972 graph_.closeArc(endingArc, saddleId);
973 PRINT("/" << graph_.printArc(endingArc));
974 if(graph_.getArc(endingArc).isVisible()) {
975 ++visibleClosed;
976 }
977 Propagation *arcProp = graph_.getArc(endingArc).getPropagation();
978 localProp->merge(*arcProp);
979 }
980 return visibleClosed;
981}
982
983template <typename ScalarType, typename triangulationType>
986 const idNode saddleId, const std::set<DynGraphNode<idVertex> *> &compVect) {
987 // version for the sequential arc growth, do not merge the propagations
988
989#ifndef TTK_ENABLE_KAMIKAZE
990 if(compVect.size() < 2) {
991 std::cerr << "[FTR]: merge at saddle with only one lower CC" << std::endl;
992 }
993#endif
994
995 idSuperArc visibleClosed = 0;
996 for(auto *dgNode : compVect) {
997 const idSuperArc endingArc = dgNode->getCorArc();
998 graph_.closeArc(endingArc, saddleId);
999 PRINT("/" << graph_.printArc(endingArc));
1000 if(graph_.getArc(endingArc).isVisible()) {
1001 ++visibleClosed;
1002 }
1003 }
1004 return visibleClosed;
1005}
1006
1007template <typename ScalarType, typename triangulationType>
1009 Propagation *const localProp,
1010 const std::set<DynGraphNode<idVertex> *> &compVect,
1011 const bool hidden) {
1012 const idVertex curVert = localProp->getCurVertex();
1013 const idNode curNode = graph_.getNodeId(curVert);
1014
1015 for(auto *dgNode : compVect) {
1016 const idSuperArc newArc = graph_.openArc(curNode, localProp);
1017 dgNode->setRootArc(newArc);
1018 visit(localProp, newArc);
1019
1020 if(hidden)
1021 graph_.getArc(newArc).hide();
1022
1023 PRINT("v" << graph_.printArc(newArc));
1024 }
1025}
1026
1027template <typename ScalarType, typename triangulationType>
1029 Propagation *const localProp, const idSuperArc curArc) {
1030 const idVertex curVert = localProp->getCurVertex();
1031 Visit opposite;
1032 idSuperArc retArc;
1033 // #pragma omp critical
1034 {
1035 propagations_.visit(curVert, localProp);
1036 opposite = propagations_.visitOpposite(curVert, localProp);
1037 bool done;
1038#ifdef TTK_ENABLE_OPENMP
1039#pragma omp atomic read seq_cst
1040#endif
1041 done = opposite.done;
1042 if(!done) {
1043 graph_.visit(curVert, curArc);
1044 retArc = nullSuperArc;
1045 } else {
1046 retArc = graph_.getArcId(curVert);
1047 }
1048 }
1049 return retArc;
1050}
1051
1053
1054template <typename ScalarType, typename triangulationType>
1057 const idVertex leaf, const bool fromMin) {
1058 VertCompFN comp;
1059 if(fromMin)
1060 comp = [&](idVertex a, idVertex b) { return scalars_.isHigher(a, b); };
1061 else
1062 comp = [&](idVertex a, idVertex b) { return scalars_.isLower(a, b); };
1063 return propagations_.newPropagation(leaf, comp, fromMin);
1064}
1065
1066template <typename ScalarType, typename triangulationType>
1068 const orderedEdge &e0,
1069 const orderedEdge &e1,
1070 const Propagation *const localProp) {
1071 const idVertex end0 = std::get<1>(e0);
1072 const idVertex end1 = std::get<1>(e1);
1073
1074 if(localProp->compare(end0, end1)) {
1075 if(localProp->goDown()) {
1076 return -scalars_.getMirror(end0);
1077 }
1078 return scalars_.getMirror(end0);
1079 }
1080
1081 if(localProp->goDown()) {
1082 return -scalars_.getMirror(end1);
1083 }
1084 return scalars_.getMirror(end1);
1085}
1086
1087template <typename ScalarType, typename triangulationType>
1090 const orderedTriangle &oTriangle,
1091 const Propagation *const localProp) const {
1092 orderedEdge firstEdge
1093 = mesh_.getOrderedEdge(std::get<0>(oTriangle), localProp->goUp());
1094 if(std::get<0>(firstEdge) == localProp->getCurVertex()) {
1095 return vertPosInTriangle::Start;
1096 } else if(std::get<1>(firstEdge) == localProp->getCurVertex()) {
1097 return vertPosInTriangle::Middle;
1098 } else {
1099 return vertPosInTriangle::End;
1100 }
1101}
1102
1103template <typename ScalarType, typename triangulationType>
1106 const orderedTriangle &oTriangle,
1107 const Propagation *const localProp) const {
1108 const orderedEdge &higherEdge
1109 = mesh_.getOrderedEdge(std::get<1>(oTriangle), localProp->goUp());
1110 return std::get<1>(higherEdge);
1111}
1112
1113template <typename ScalarType, typename triangulationType>
1116 const orderedTriangle oTri, const idVertex v0, const idVertex v1) {
1117 idVertex edge0Vert, edge1Vert;
1118
1119 mesh_.getEdgeVertex(std::get<0>(oTri), 0, edge0Vert);
1120 mesh_.getEdgeVertex(std::get<0>(oTri), 1, edge1Vert);
1121 if((edge0Vert == v0 && edge1Vert == v1)
1122 || (edge0Vert == v1 && edge1Vert == v0)) {
1123 return std::get<0>(oTri);
1124 }
1125
1126 mesh_.getEdgeVertex(std::get<1>(oTri), 0, edge0Vert);
1127 mesh_.getEdgeVertex(std::get<1>(oTri), 1, edge1Vert);
1128 if((edge0Vert == v0 && edge1Vert == v1)
1129 || (edge0Vert == v1 && edge1Vert == v0)) {
1130 return std::get<1>(oTri);
1131 }
1132
1133#ifndef TTK_ENABLE_KAMIKAZE
1134 mesh_.getEdgeVertex(std::get<2>(oTri), 0, edge0Vert);
1135 mesh_.getEdgeVertex(std::get<2>(oTri), 1, edge1Vert);
1136 if((edge0Vert == v0 && edge1Vert == v1)
1137 || (edge0Vert == v1 && edge1Vert == v0)) {
1138 return std::get<2>(oTri);
1139 }
1140
1141 std::cout << "[FTR]: edge not found in triangle " << v0 << " " << v1
1142 << std::endl;
1143 return nullEdge;
1144#else
1145 return std::get<2>(oTri);
1146#endif
1147}
#define TTK_FORCE_USE(x)
Force the compiler to use the function/method parameter.
Definition: BaseClass.h:57
#define ttkNotUsed(x)
Mark function/method parameters that are not used in the function body at all.
Definition: BaseClass.h:47
#define PRINT(msg)
TTK processing package that efficiently computes the contour tree of scalar data and more (data segme...
TTK FTRGraph processing package.
Definition: FTRGraph.h:77
TTK fTRGraph propagation management with Fibonacci heaps.
vertPosInTriangle
position of a vertex in a triangle
Definition: FTRDataTypes.h:49
long unsigned int idSuperArc
SuperArc index in vect_superArcs_.
Definition: FTRDataTypes.h:25
SimplexId idVertex
Vertex index in scalars_.
Definition: FTRDataTypes.h:29
SimplexId idEdge
Edge index in vect_edgeList_.
Definition: FTRDataTypes.h:31