TTK
Loading...
Searching...
No Matches
ContourTree.cpp
Go to the documentation of this file.
1#include <ContourTree.h>
2
3#include <cmath>
4#include <iomanip>
5#include <queue>
6#include <set>
7
8using namespace std;
9using namespace ttk;
10
11struct MyCmp {
12 const vector<double> *vertexScalars_;
13 const vector<int> *vertexOffsets_;
14 const vector<Node> *nodeList_;
16
17 MyCmp(const vector<double> *vertexScalars,
18 const vector<int> *vertexOffsets,
19 const vector<Node> *nodeList,
20 bool isAscendingOrder)
21 : vertexScalars_(vertexScalars), vertexOffsets_(vertexOffsets),
22 nodeList_(nodeList), isAscendingOrder_(isAscendingOrder) {
23 }
24
25 bool operator()(int node1, int node2) {
26 int vertex1 = (*nodeList_)[node1].getVertexId();
27 int vertex2 = (*nodeList_)[node2].getVertexId();
28 bool cmp = ((*vertexScalars_)[vertex1] < (*vertexScalars_)[vertex2])
29 || ((*vertexScalars_)[vertex1] == (*vertexScalars_)[vertex2]
30 && (*vertexOffsets_)[vertex1] < (*vertexOffsets_)[vertex2]);
32 return cmp;
33 else
34 return !cmp;
35 }
36};
37
39 bool operator()(const pair<bool, pair<double, pair<int, int>>> &v0,
40 const pair<bool, pair<double, pair<int, int>>> &v1) const {
41
42 if(v0.first) {
43 return ((v0.second.first < v1.second.first)
44 || ((v0.second.first == v1.second.first)
45 && (v0.second.second.first < v1.second.second.first)));
46 } else {
47 return ((v0.second.first > v1.second.first)
48 || ((v0.second.first == v1.second.first)
49 && (v0.second.second.first > v1.second.second.first)));
50 }
51 }
53
55 bool operator()(const pair<pair<int, int>, double> &p0,
56 const pair<pair<int, int>, double> &p1) {
57
58 return p0.second < p1.second;
59 }
60
62
64 _persistenceCmp3(const vector<double> *vertexScalars)
65 : vertexScalars_{vertexScalars} {
66 }
67
68 bool operator()(const pair<pair<int, int>, double> &p0,
69 const pair<pair<int, int>, double> &p1) {
70 return (*vertexScalars_)[p0.first.first]
71 > (*vertexScalars_)[p1.first.first];
72 }
73
74 const vector<double> *vertexScalars_;
75};
76
78 bool operator()(const pair<double, double> &p0,
79 const pair<double, double> &p1) {
80 return p0.first < p1.first;
81 }
83
84void SuperArc::smooth(const vector<Node> &nodeList,
85 const vector<vector<double>> *vertexPositions,
86 bool order) {
87 int N = getNumberOfBarycenters();
88 if(N) {
90 vector<vector<double>> barycenterList(barycenterList_.size());
91 for(unsigned int i = 0; i < barycenterList.size(); ++i)
92 barycenterList[i].resize(3);
93
94 int up_vId;
95 int down_vId;
96 if(order) {
97 up_vId = nodeList[getUpNodeId()].getVertexId();
98 down_vId = nodeList[getDownNodeId()].getVertexId();
99 } else {
100 up_vId = nodeList[getDownNodeId()].getVertexId();
101 down_vId = nodeList[getUpNodeId()].getVertexId();
102 }
103
104 double p0[3];
105 double p1[3];
106 for(unsigned int k = 0; k < 3; ++k) {
107 p0[k] = (*vertexPositions)[down_vId][k];
108 p1[k] = (*vertexPositions)[up_vId][k];
109 }
110
112 if(N > 1) {
113 // first
114 for(unsigned int k = 0; k < 3; ++k)
115 barycenterList[0][k] = (p0[k] + barycenterList_[1][k]) * 0.5;
116
117 // main
118 for(int i = 1; i < N - 1; ++i) {
119 for(unsigned int k = 0; k < 3; ++k)
120 barycenterList[i][k]
121 = (barycenterList_[i - 1][k] + barycenterList_[i + 1][k]) * 0.5;
122 }
123 // last
124 for(unsigned int k = 0; k < 3; ++k)
125 barycenterList[N - 1][k] = (p1[k] + barycenterList_[N - 1][k]) * 0.5;
126 } else {
127 for(unsigned int k = 0; k < 3; ++k)
128 barycenterList[0][k] = (p0[k] + p1[k]) * 0.5;
129 }
130
132 for(int i = 0; i < N; ++i) {
133 for(unsigned int k = 0; k < 3; ++k)
134 barycenterList_[i][k] = barycenterList[i][k];
135 }
136 }
137}
138
139void SuperArc::sortRegularNodes(const vector<double> *vertexScalars,
140 const vector<int> *vertexOffsets,
141 const vector<Node> *nodeList,
142 bool order) {
143 MyCmp cmp(vertexScalars, vertexOffsets, nodeList, order);
144
145 if(order)
146 std::sort(regularNodeList_.begin(), regularNodeList_.end(), cmp);
147 else
148 std::sort(regularNodeList_.rbegin(), regularNodeList_.rend(), cmp);
149}
150
152 const int &downVertexId,
153 const int &upVertexId,
154 const vector<int> &ttkNotUsed(interiorNodeIds)) {
155
156 if(!tree_)
157 return -DBL_MAX;
158
159 double downScalar = 0, upScalar = 0;
160
161 tree_->getVertexScalar(downVertexId, downScalar);
162 tree_->getVertexScalar(upVertexId, upScalar);
163
164 return fabs(upScalar - downScalar);
165}
166
168 this->setDebugMsgPrefix("SubLevelSetTree");
169}
170
171int SubLevelSetTree::appendRegularNode(const int &superArcId,
172 const int &nodeId) {
173
174 superArcList_[superArcId].appendRegularNode(nodeId);
175 vertex2superArc_[nodeList_[nodeId].getVertexId()] = superArcId;
176 vertex2superArcNode_[nodeList_[nodeId].getVertexId()]
177 = superArcList_[superArcId].getNumberOfRegularNodes() - 1;
178
179 // create a low level arc with the previous node
180 int nId = -1;
181 if(superArcList_[superArcId].getNumberOfRegularNodes() == 1) {
182 // the down neighbor is the arc extremity
183 nId = superArcList_[superArcId].getDownNodeId();
184 } else {
185 nId = superArcList_[superArcId].getRegularNodeId(
186 superArcList_[superArcId].getNumberOfRegularNodes() - 2);
187 }
188
189 makeArc(nId, nodeId);
190
191 return 0;
192}
193
195
196 Timer timer;
197
198 if(!vertexNumber_)
199 return -1;
200 if((!vertexScalars_) || ((int)vertexScalars_->size() != vertexNumber_))
201 return -2;
202 if((!vertexSoSoffsets_) || ((int)vertexSoSoffsets_->size() != vertexNumber_))
203 return -3;
204 // TODO uncomment after triangulation templatization
205 // if(triangulation_->isEmpty())
206 // return -4;
207 if((!minimumList_) && (!maximumList_))
208 return -5;
209 if((minimumList_) && (maximumList_))
210 return -6;
211
212 vector<UnionFind> seeds;
213 vector<vector<int>> seedSuperArcs;
214 vector<UnionFind *> vertexSeeds(vertexNumber_, (UnionFind *)nullptr);
215 vector<UnionFind *> starSets;
216 vector<bool> visitedVertices(vertexNumber_, false);
217
218 SimplexId vertexId = -1, nId = -1;
219 UnionFind *seed = nullptr, *firstUf = nullptr;
220
221 const vector<int> *extremumList = nullptr;
222
223 bool isMergeTree = true;
224
225 if(minimumList_) {
226 extremumList = minimumList_;
227 } else {
228 isMergeTree = false;
229 extremumList = maximumList_;
230 }
231
232 set<
233 // is merging
234 pair<bool,
235 // function value
236 pair<double,
237 // vertex offset
238 pair<int,
239 // vertex id
240 int>>>,
242 filtrationFront;
243
244 seeds.resize(extremumList->size());
245 seedSuperArcs.resize(seeds.size());
246
247 for(int i = 0; i < (int)extremumList->size(); i++) {
248 // link each minimum to a union find seed
249 vertexSeeds[(*extremumList)[i]] = &(seeds[i]);
250
251 // open an arc
252 seedSuperArcs[i].push_back(openSuperArc(makeNode((*extremumList)[i])));
253
254 // add each minimum to the filtration front
255 pair<bool, pair<double, pair<int, int>>> v;
256 v.first = isMergeTree;
257 v.second.first = (*vertexScalars_)[(*extremumList)[i]];
258 v.second.second.first = (*vertexSoSoffsets_)[(*extremumList)[i]];
259 v.second.second.second = (*extremumList)[i];
260
261 filtrationFront.insert(v);
262 }
263
264 bool merge = false;
265
266 // filtration loop
267 do {
268
269 vertexId = filtrationFront.begin()->second.second.second;
270 filtrationFront.erase(filtrationFront.begin());
271
272 starSets.clear();
273
274 merge = false;
275 firstUf = nullptr;
276
277 SimplexId neighborNumber
279 for(SimplexId i = 0; i < neighborNumber; i++) {
280 triangulation_->getVertexNeighbor(vertexId, i, nId);
281
282 if(vertexSeeds[nId]) {
283 seed = vertexSeeds[nId]->find();
284 starSets.push_back(seed);
285
286 // is it merging things?
287 if(!firstUf)
288 firstUf = seed;
289 else if(seed != firstUf)
290 merge = true;
291 }
292
293 if(!visitedVertices[nId]) {
294
295 pair<bool, pair<double, pair<int, int>>> v;
296 v.first = isMergeTree;
297 v.second.first = (*vertexScalars_)[nId];
298 v.second.second.first = (*vertexSoSoffsets_)[nId];
299 v.second.second.second = nId;
300
301 filtrationFront.insert(v);
302 visitedVertices[nId] = true;
303 }
304 }
305
306 if(!vertexSeeds[vertexId]) {
307
308 vertexSeeds[vertexId] = UnionFind::makeUnion(starSets);
309
310 int newNodeId = makeNode(vertexId);
311
312 if(merge) {
313
314 vector<int> seedIds;
315 for(int i = 0; i < (int)starSets.size(); i++) {
316 int seedId = starSets[i] - &(seeds[0]);
317 bool found = false;
318 for(int j = 0; j < (int)seedIds.size(); j++) {
319 if(seedIds[j] == seedId) {
320 found = true;
321 break;
322 }
323 }
324 if(!found)
325 seedIds.push_back(seedId);
326 }
327
328 for(int i = 0; i < (int)seedIds.size(); i++) {
329 int superArcId
330 = seedSuperArcs[seedIds[i]][seedSuperArcs[seedIds[i]].size() - 1];
331 closeSuperArc(superArcId, newNodeId);
332 }
333
334 int seedId = vertexSeeds[vertexId] - &(seeds[0]);
335 if(!filtrationFront.empty())
336 seedSuperArcs[seedId].push_back(openSuperArc(newNodeId));
337 } else if(starSets.size()) {
338 // we're dealing with a degree-2 node
339 int seedId = starSets[0] - &(seeds[0]);
340 int superArcId
341 = seedSuperArcs[seedId][seedSuperArcs[seedId].size() - 1];
342
343 if(filtrationFront.empty()) {
344 // last vertex
345 closeSuperArc(superArcId, newNodeId);
346 } else {
348 appendRegularNode(superArcId, newNodeId);
349 }
350 }
351 }
352
353 visitedVertices[vertexId] = true;
354
355 } while(!filtrationFront.empty());
356
357 const std::string tree = isMergeTree ? "JoinTree" : "SplitTree";
358 this->printMsg(
359 tree + " computed", 1.0, timer.getElapsedTime(), this->threadNumber_);
360 this->printMsg(std::vector<std::vector<std::string>>{
361 {"#Nodes", std::to_string(getNumberOfNodes())},
362 {"#Arcs", std::to_string(getNumberOfArcs())}});
363
364 print();
365
366 return 0;
367}
368
369int SubLevelSetTree::clearArc(const int &vertexId0, const int &vertexId1) {
370
371 if((vertexId0 < 0) || (vertexId0 >= vertexNumber_))
372 return -1;
373 if((vertexId1 < 0) || (vertexId1 >= vertexNumber_))
374 return -2;
375
376 int nodeId0 = vertex2node_[vertexId0], nodeId1 = vertex2node_[vertexId1];
377
378 for(int i = 0; i < nodeList_[nodeId0].getNumberOfUpArcs(); i++) {
379 if(nodeList_[arcList_[nodeList_[nodeId0].getUpArcId(i)].getUpNodeId()]
380 .getVertexId()
381 == vertexId1) {
382 nodeList_[nodeId0].removeUpArcId(i);
383 break;
384 }
385 }
386
387 for(int i = 0; i < nodeList_[nodeId1].getNumberOfDownArcs(); i++) {
388 if(nodeList_[arcList_[nodeList_[nodeId1].getDownArcId(i)].getDownNodeId()]
389 .getVertexId()
390 == vertexId0) {
391 nodeList_[nodeId1].removeDownArcId(i);
392 break;
393 }
394 }
395
396 return 0;
397}
398
399int SubLevelSetTree::clearRegularNode(const int &vertexId) {
400
401 if((vertexId < 0) || (vertexId >= vertexNumber_))
402 return -1;
403
404 int nodeId = vertex2node_[vertexId];
405
406 if(!((nodeList_[nodeId].getNumberOfDownArcs() == 1)
407 && (nodeList_[nodeId].getNumberOfUpArcs() == 1)))
408 return -2;
409
410 int downArcId = nodeList_[nodeId].getDownArcId(0),
411 upArcId = nodeList_[nodeId].getUpArcId(0);
412
413 int upNodeId = arcList_[upArcId].getUpNodeId();
414
415 // removes the arcs from the current node
416 nodeList_[nodeId].removeDownArcId(0);
417 nodeList_[nodeId].removeUpArcId(0);
418
419 // update the arc from the down node
420 arcList_[downArcId].setUpNodeId(upNodeId);
421
422 // remove the arc from the up node
423 for(int i = 0; i < nodeList_[upNodeId].getNumberOfDownArcs(); i++) {
424 if(nodeList_[upNodeId].getDownArcId(i) == upArcId) {
425 nodeList_[upNodeId].removeDownArcId(i);
426 break;
427 }
428 }
429 // and add the updated arc
430 nodeList_[upNodeId].addDownArcId(downArcId);
431
432 return 0;
433}
434
435int SubLevelSetTree::clearRoot(const int &vertexId) {
436
437 if((vertexId < 0) || (vertexId >= vertexNumber_))
438 return -1;
439
440 int nodeId = vertex2node_[vertexId];
441
442 if((nodeList_[nodeId].getNumberOfDownArcs() == 1)
443 && (!nodeList_[nodeId].getNumberOfUpArcs())) {
444 clearArc(
445 nodeList_[arcList_[nodeList_[nodeId].getDownArcId(0)].getDownNodeId()]
446 .getVertexId(),
447 vertexId);
448 }
449
450 return 0;
451}
452
454 const int &pointId,
455 vector<int> &vertexIds,
456 const vector<float> *origin,
457 const vector<float> *voxelSize,
458 ofstream &o) {
459
460 vector<float> myOrigin(3), myVoxelSize(3);
461
462 if(origin) {
463 myOrigin = *(origin);
464 } else {
465 myOrigin[0] = myOrigin[1] = myOrigin[2] = 0;
466 }
467
468 if(voxelSize) {
469 myVoxelSize = *(voxelSize);
470 } else {
471 myVoxelSize[0] = myVoxelSize[1] = myVoxelSize[2] = 1;
472 }
473
474 double offset = myVoxelSize[0];
475
476 vector<double> v;
477
478 int downNodeId = superArcList_[arcId].downNodeId_;
479 int upNodeId = superArcList_[arcId].upNodeId_;
480
481 v = (*vertexPositions_)[nodeList_[downNodeId].vertexId_];
482 // fix a bug in paraview
483 for(int i = 0; i < 3; i++) {
484 v[i] -= myOrigin[i];
485 v[i] /= myVoxelSize[i];
486 }
487 v[0] -= offset;
488 v[2] += offset;
489 o << v[0] << " " << v[1] << " " << v[2] << endl;
490 vertexIds.push_back(pointId);
491
492 v = (*vertexPositions_)[nodeList_[downNodeId].vertexId_];
493 // fix a bug in paraview
494 for(int i = 0; i < 3; i++) {
495 v[i] -= myOrigin[i];
496 v[i] /= myVoxelSize[i];
497 }
498 v[0] += offset;
499 v[2] += offset;
500 o << v[0] << " " << v[1] << " " << v[2] << endl;
501 vertexIds.push_back(pointId + 1);
502
503 v = (*vertexPositions_)[nodeList_[downNodeId].vertexId_];
504 // fix a bug in paraview
505 for(int i = 0; i < 3; i++) {
506 v[i] -= myOrigin[i];
507 v[i] /= myVoxelSize[i];
508 }
509 v[0] += offset;
510 v[2] -= offset;
511 o << v[0] << " " << v[1] << " " << v[2] << endl;
512 vertexIds.push_back(pointId + 2);
513
514 v = (*vertexPositions_)[nodeList_[downNodeId].vertexId_];
515 // fix a bug in paraview
516 for(int i = 0; i < 3; i++) {
517 v[i] -= myOrigin[i];
518 v[i] /= myVoxelSize[i];
519 }
520 v[0] -= offset;
521 v[2] -= offset;
522 o << v[0] << " " << v[1] << " " << v[2] << endl;
523 vertexIds.push_back(pointId + 3);
524
525 v = (*vertexPositions_)[nodeList_[upNodeId].vertexId_];
526 // fix a bug in paraview
527 for(int i = 0; i < 3; i++) {
528 v[i] -= myOrigin[i];
529 v[i] /= myVoxelSize[i];
530 }
531 v[0] -= offset;
532 v[2] += offset;
533 o << v[0] << " " << v[1] << " " << v[2] << endl;
534 vertexIds.push_back(pointId + 4);
535
536 v = (*vertexPositions_)[nodeList_[upNodeId].vertexId_];
537 // fix a bug in paraview
538 for(int i = 0; i < 3; i++) {
539 v[i] -= myOrigin[i];
540 v[i] /= myVoxelSize[i];
541 }
542 v[0] += offset;
543 v[2] += offset;
544 o << v[0] << " " << v[1] << " " << v[2] << endl;
545 vertexIds.push_back(pointId + 5);
546
547 v = (*vertexPositions_)[nodeList_[upNodeId].vertexId_];
548 // fix a bug in paraview
549 for(int i = 0; i < 3; i++) {
550 v[i] -= myOrigin[i];
551 v[i] /= myVoxelSize[i];
552 }
553 v[0] += offset;
554 v[2] -= offset;
555 o << v[0] << " " << v[1] << " " << v[2] << endl;
556 vertexIds.push_back(pointId + 6);
557
558 v = (*vertexPositions_)[nodeList_[upNodeId].vertexId_];
559 // fix a bug in paraview
560 for(int i = 0; i < 3; i++) {
561 v[i] -= myOrigin[i];
562 v[i] /= myVoxelSize[i];
563 }
564 v[0] -= offset;
565 v[2] -= offset;
566 o << v[0] << " " << v[1] << " " << v[2] << endl;
567 vertexIds.push_back(pointId + 7);
568
569 return 0;
570}
571
572int SubLevelSetTree::exportNodeColorToVtk(const int &nodeId, ofstream &o) {
573
574 for(int i = 0; i < 8; i++) {
575
576 if(!nodeList_[nodeId].downSuperArcList_.size()) {
577 // extremum
578 if(minimumList_) {
579 // minimum
580 o << "0 ";
581 } else {
582 o << "1 ";
583 }
584 } else if(!nodeList_[nodeId].upSuperArcList_.size()) {
585 if(minimumList_) {
586 // maximum
587 o << "1 ";
588 } else {
589 o << "0 ";
590 }
591 } else {
592 o << "0.5 ";
593 }
594 }
595
596 return 0;
597}
598
600 const int &pointId,
601 vector<int> &vertexIds,
602 const vector<float> *origin,
603 const vector<float> *voxelSize,
604 ofstream &o) {
605
606 vector<float> myOrigin(3), myVoxelSize(3);
607
608 if(origin) {
609 myOrigin = *(origin);
610 } else {
611 myOrigin[0] = myOrigin[1] = myOrigin[2] = 0;
612 }
613
614 if(voxelSize) {
615 myVoxelSize = *(voxelSize);
616 } else {
617 myVoxelSize[0] = myVoxelSize[1] = myVoxelSize[2] = 1;
618 }
619
620 double offset = myVoxelSize[0];
621
622 vector<double> v;
623
624 v = (*vertexPositions_)[nodeList_[nodeId].vertexId_];
625 // fix a bug in paraview
626 for(int i = 0; i < 3; i++) {
627 v[i] -= myOrigin[i];
628 v[i] /= myVoxelSize[i];
629 }
630
631 v[0] -= offset;
632 v[1] -= offset;
633 v[2] -= offset;
634 o << v[0] << " " << v[1] << " " << v[2] << endl;
635 vertexIds.push_back(pointId);
636
637 v = (*vertexPositions_)[nodeList_[nodeId].vertexId_];
638 // fix a bug in paraview
639 for(int i = 0; i < 3; i++) {
640 v[i] -= myOrigin[i];
641 v[i] /= myVoxelSize[i];
642 }
643 v[0] += offset;
644 v[1] -= offset;
645 v[2] -= offset;
646 o << v[0] << " " << v[1] << " " << v[2] << endl;
647 vertexIds.push_back(pointId + 1);
648
649 v = (*vertexPositions_)[nodeList_[nodeId].vertexId_];
650 // fix a bug in paraview
651 for(int i = 0; i < 3; i++) {
652 v[i] -= myOrigin[i];
653 v[i] /= myVoxelSize[i];
654 }
655 v[0] += offset;
656 v[1] += offset;
657 v[2] -= offset;
658 o << v[0] << " " << v[1] << " " << v[2] << endl;
659 vertexIds.push_back(pointId + 2);
660
661 v = (*vertexPositions_)[nodeList_[nodeId].vertexId_];
662 // fix a bug in paraview
663 for(int i = 0; i < 3; i++) {
664 v[i] -= myOrigin[i];
665 v[i] /= myVoxelSize[i];
666 }
667 v[0] -= offset;
668 v[1] += offset;
669 v[2] -= offset;
670 o << v[0] << " " << v[1] << " " << v[2] << endl;
671 vertexIds.push_back(pointId + 3);
672
673 v = (*vertexPositions_)[nodeList_[nodeId].vertexId_];
674 // fix a bug in paraview
675 for(int i = 0; i < 3; i++) {
676 v[i] -= myOrigin[i];
677 v[i] /= myVoxelSize[i];
678 }
679 v[0] -= offset;
680 v[1] -= offset;
681 v[2] += offset;
682 o << v[0] << " " << v[1] << " " << v[2] << endl;
683 vertexIds.push_back(pointId + 4);
684
685 v = (*vertexPositions_)[nodeList_[nodeId].vertexId_];
686 // fix a bug in paraview
687 for(int i = 0; i < 3; i++) {
688 v[i] -= myOrigin[i];
689 v[i] /= myVoxelSize[i];
690 }
691 v[0] += offset;
692 v[1] -= offset;
693 v[2] += offset;
694 o << v[0] << " " << v[1] << " " << v[2] << endl;
695 vertexIds.push_back(pointId + 5);
696
697 v = (*vertexPositions_)[nodeList_[nodeId].vertexId_];
698 // fix a bug in paraview
699 for(int i = 0; i < 3; i++) {
700 v[i] -= myOrigin[i];
701 v[i] /= myVoxelSize[i];
702 }
703 v[0] += offset;
704 v[1] += offset;
705 v[2] += offset;
706 o << v[0] << " " << v[1] << " " << v[2] << endl;
707 vertexIds.push_back(pointId + 6);
708
709 v = (*vertexPositions_)[nodeList_[nodeId].vertexId_];
710 // fix a bug in paraview
711 for(int i = 0; i < 3; i++) {
712 v[i] -= myOrigin[i];
713 v[i] /= myVoxelSize[i];
714 }
715 v[0] -= offset;
716 v[1] += offset;
717 v[2] += offset;
718 o << v[0] << " " << v[1] << " " << v[2] << endl;
719 vertexIds.push_back(pointId + 7);
720
721 return 0;
722}
723
724int SubLevelSetTree::exportPersistenceCurve(const string &fileName) const {
725 vector<pair<double, int>> persistencePlot;
726
727 getPersistencePlot(persistencePlot);
728
729 ofstream file(fileName.data(), ios::out);
730
731 if(!file) {
732 return -1;
733 }
734
735 for(int i = 0; i < (int)persistencePlot.size(); i++) {
736
737 file << setprecision(REAL_SIGNIFICANT_DIGITS) << persistencePlot[i].first
738 << " " << persistencePlot[i].second << endl;
739 }
740
741 file.close();
742
743 return 0;
744}
745
746int SubLevelSetTree::exportPersistenceDiagram(const string &fileName) const {
747
748 vector<pair<double, double>> diagram;
749
750 getPersistenceDiagram(diagram);
751
752 ofstream file(fileName.data(), ios::out);
753
754 if(!file) {
755 return -1;
756 }
757
758 for(int i = 0; i < (int)diagram.size(); i++) {
759 file << setprecision(REAL_SIGNIFICANT_DIGITS) << diagram[i].first << " "
760 << diagram[i].first << endl;
761 file << setprecision(REAL_SIGNIFICANT_DIGITS) << diagram[i].first << " "
762 << diagram[i].second << endl;
763 file << setprecision(REAL_SIGNIFICANT_DIGITS) << diagram[i].first << " "
764 << diagram[i].first << endl;
765 }
766
767 file.close();
768
769 return 0;
770}
771
772int SubLevelSetTree::exportToSvg(const string &fileName,
773 const double &scaleX,
774 const double &scaleY) {
775
776 Timer t;
777
778 if(!vertexScalars_)
779 return -1;
780
781 bool hasLayout = buildPlanarLayout(scaleX, scaleY);
782
783 // put the root right above the highest saddle
784 double minY = -1, maxY = -1;
785 for(int i = 0; i < (int)superArcList_.size(); i++) {
786 const SuperArc *a = &(superArcList_[i]);
787 const Node *n = &(nodeList_[a->getDownNodeId()]);
788 if((!i) || (minY > n->layoutY_)) {
789 minY = n->layoutY_;
790 }
791 if((!i) || (maxY < n->layoutY_)) {
792 maxY = n->layoutY_;
793 }
794
795 n = &(nodeList_[a->getUpNodeId()]);
796 if((!i) || (minY > n->layoutY_)) {
797 minY = n->layoutY_;
798 }
799 if((!i) || (maxY < n->layoutY_)) {
800 maxY = n->layoutY_;
801 }
802 }
803
804 for(int i = 0; i < (int)superArcList_.size(); i++) {
805 const SuperArc *a = &(superArcList_[i]);
806 Node *upNode = &(nodeList_[a->getUpNodeId()]);
807 Node *downNode = &(nodeList_[a->getDownNodeId()]);
808
809 if(!upNode->getNumberOfUpSuperArcs()) {
810 // root
811 if(minimumList_) {
812 // join tree
813 upNode->layoutY_
814 = 0.4 * (downNode->layoutY_ - minY) + downNode->layoutY_;
815 break;
816 }
817 if(maximumList_) {
818 // split tree
819 upNode->layoutY_
820 = downNode->layoutY_ - 0.4 * (maxY - downNode->layoutY_);
821 break;
822 }
823 }
824 }
825
826 string dotFileName = fileName + ".dot";
827
828 ofstream dotFile(dotFileName.data(), ios::out);
829
830 if(!dotFile) {
831 this->printErr("Could not open file `" + dotFileName + "'!");
832 return -2;
833 }
834
835 dotFile << "digraph \"";
836 if(minimumList_)
837 dotFile << "Join";
838 else
839 dotFile << "Split";
840 dotFile << " Tree\"{" << endl;
841
842 double minValue = 0, maxValue = 0;
843 for(int i = 0; i < (int)vertexScalars_->size(); i++) {
844 if((!i) || (minValue > (*vertexScalars_)[i])) {
845 minValue = (*vertexScalars_)[i];
846 }
847 if((!i) || (maxValue < (*vertexScalars_)[i])) {
848 maxValue = (*vertexScalars_)[i];
849 }
850 }
851
852 // super-arc first!
853 for(int i = 0; i < (int)superArcList_.size(); i++) {
854
855 if(!superArcList_[i].pruned_) {
856
857 int downNodeId = superArcList_[i].getDownNodeId();
858 int upNodeId = superArcList_[i].getUpNodeId();
859
860 dotFile << " \"f = ";
861 dotFile << (*vertexScalars_)[nodeList_[downNodeId].vertexId_];
862 if(vertexPositions_) {
863 dotFile << "\\n p = ("
864 << (*vertexPositions_)[nodeList_[downNodeId].vertexId_][0]
865 << " "
866 << (*vertexPositions_)[nodeList_[downNodeId].vertexId_][1]
867 << " "
868 << (*vertexPositions_)[nodeList_[downNodeId].vertexId_][2]
869 << ")";
870 }
871 dotFile << "\" -> \"f = ";
872 dotFile << (*vertexScalars_)[nodeList_[upNodeId].vertexId_];
873 if(vertexPositions_) {
874 dotFile << "\\n p = ("
875 << (*vertexPositions_)[nodeList_[upNodeId].vertexId_][0] << " "
876 << (*vertexPositions_)[nodeList_[upNodeId].vertexId_][1] << " "
877 << (*vertexPositions_)[nodeList_[upNodeId].vertexId_][2] << ")";
878 }
879 dotFile << "\"" << endl;
880 }
881 }
882
883 for(int i = 0; i < (int)nodeList_.size(); i++) {
884
885 if((nodeList_[i].downSuperArcList_.size())
886 || (nodeList_[i].upSuperArcList_.size())) {
887 // not a regular node
888 if(!nodeList_[i].pruned_) {
889 dotFile << " \"f = ";
890 dotFile << (*vertexScalars_)[nodeList_[i].vertexId_];
891 if(vertexPositions_) {
892 dotFile << "\\n p = ("
893 << (*vertexPositions_)[nodeList_[i].vertexId_][0] << " "
894 << (*vertexPositions_)[nodeList_[i].vertexId_][1] << " "
895 << (*vertexPositions_)[nodeList_[i].vertexId_][2] << ")";
896 }
897 dotFile << "\" [";
898
899 // properties
900 if(minimumList_) {
901 // join tree
902
903 // color
904 if(!nodeList_[i].upSuperArcList_.size()) {
905 // global maximum --> green
906 dotFile << "fillcolor=green,";
907 } else if(!nodeList_[i].downSuperArcList_.size()) {
908 // minimum --> blue
909 dotFile << "fillcolor=blue,fontcolor=white,";
910 } else {
911 dotFile << "fillcolor=lightyellow,shape=diamond,";
912 }
913
914 } else {
915 // split tree
916
917 // color
918 if(!nodeList_[i].upSuperArcList_.size()) {
919 // global minimum --> blue
920 dotFile << "fillcolor=blue,fontcolor=white,";
921 } else if(!nodeList_[i].downSuperArcList_.size()) {
922 // maximum --> green
923 dotFile << "fillcolor=green,";
924 } else {
925 dotFile << "fillcolor=lightyellow,shape=diamond,";
926 }
927 }
928
929 dotFile << "style=filled,";
930
931 if(hasLayout) {
932 dotFile << "pos=\"" << nodeList_[i].layoutX_ << ", "
933 << nodeList_[i].layoutY_ << "!\"";
934 }
935
936 dotFile << "]";
937
938 dotFile << endl;
939 }
940 }
941 }
942
943 dotFile << "}" << endl;
944
945 dotFile.close();
946
947 stringstream commandLine;
948 commandLine << "dot -Kneato -Tsvg " << dotFileName << " -o " << fileName
949 << " &> /dev/null";
950 this->printMsg("Calling GraphViz to generate the SVG file.");
951 this->printMsg("This may take a long time...");
952
953 int cmdRet = system(commandLine.str().data());
954
955 if((!cmdRet) || (cmdRet == 256)) {
956 this->printMsg(
957 "Output file `" + fileName + "' generated", 1.0, t.getElapsedTime(), 1);
958 } else {
959 this->printErr("Could not find the `GraphViz' package!");
960 this->printErr("Please install it and re-run this program.");
961 }
962
963 // OsCall::rmFile(dotFileName);
964
965 return 0;
966}
967
968int SubLevelSetTree::exportToVtk(const string &fileName,
969 const vector<float> *origin,
970 const vector<float> *voxelSize) {
971
972 if(!vertexScalars_)
973 return -1;
974
976 return -2;
977
978 ofstream o(fileName.data(), ios::out);
979
980 if(!o) {
981 this->printErr("Could not open file `" + fileName + "'!");
982 return -3;
983 }
984
985 int superArcNumber = 0;
986 int nodeNumber = 0;
987 for(int i = 0; i < (int)superArcList_.size(); i++) {
988 if(!superArcList_[i].pruned_) {
989 superArcNumber++;
990 }
991 }
992 nodeNumber = superArcNumber + 1;
993
994 int pointNumber = 0, cellNumber = 0;
995 pointNumber = 8 * nodeNumber + 8 * superArcNumber;
996 cellNumber = nodeNumber + superArcNumber;
997
998 vector<bool> nodeColorOut(nodeList_.size(), false);
999 vector<bool> nodePosOut(nodeList_.size(), false);
1000 vector<bool> nodeMeshOut(nodeList_.size(), false);
1001 vector<vector<int>> nodeIds(nodeList_.size());
1002 vector<vector<int>> arcIds(superArcList_.size());
1003
1004 o << "<?xml version=\"1.0\"?>" << endl;
1005 o << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\""
1006 << " byte_order=\"LittleEndian\">" << endl;
1007 o << " <UnstructuredGrid>" << endl;
1008 o << " <Piece NumberOfPoints=\"" << pointNumber << "\" NumberOfCells=\""
1009 << cellNumber << "\">";
1010 o << " <PointData Scalars=\"Color Code\">" << endl;
1011 o << " <DataArray type=\"Float32\" Name=\"Critical Color\""
1012 << " format=\"ascii\" NumberOfComponents=\"1\">" << endl;
1013
1014 // node color
1015 for(int i = 0; i < (int)superArcList_.size(); i++) {
1016 if(!superArcList_[i].pruned_) {
1017 int downNodeId = superArcList_[i].downNodeId_;
1018 int upNodeId = superArcList_[i].upNodeId_;
1019
1020 if(!nodeColorOut[downNodeId]) {
1021 o << " ";
1022 exportNodeColorToVtk(downNodeId, o);
1023 o << endl;
1024 nodeColorOut[downNodeId] = true;
1025 }
1026 if(!nodeColorOut[upNodeId]) {
1027 o << " ";
1028 exportNodeColorToVtk(upNodeId, o);
1029 o << endl;
1030 nodeColorOut[upNodeId] = true;
1031 }
1032 }
1033 }
1034
1035 // arc color
1036 for(int i = 0; i < (int)superArcList_.size(); i++) {
1037 if(!superArcList_[i].pruned_) {
1038 o << "\t\t\t\t\t";
1039 for(int j = 0; j < 8; j++) {
1040 o << "2 ";
1041 }
1042 o << endl;
1043 }
1044 }
1045 o << " </DataArray>" << endl;
1046 o << " </PointData>" << endl;
1047
1048 // node position
1049
1050 o << " <Points>" << endl;
1051 o << " <DataArray type=\"Float32\" NumberOfComponents=\"3\""
1052 << " format=\"ascii\">" << endl;
1053
1054 int pointId = 0;
1055 for(int i = 0; i < (int)superArcList_.size(); i++) {
1056 if(!superArcList_[i].pruned_) {
1057 int downNodeId = superArcList_[i].downNodeId_;
1058 int upNodeId = superArcList_[i].upNodeId_;
1059
1060 if(!nodePosOut[downNodeId]) {
1061 o << " ";
1063 downNodeId, pointId, nodeIds[downNodeId], origin, voxelSize, o);
1064 nodePosOut[downNodeId] = true;
1065 o << endl;
1066 pointId += 8;
1067 }
1068 if(!nodePosOut[upNodeId]) {
1069 o << " ";
1071 upNodeId, pointId, nodeIds[upNodeId], origin, voxelSize, o);
1072 nodePosOut[upNodeId] = true;
1073 o << endl;
1074 pointId += 8;
1075 }
1076 }
1077 }
1078 for(int i = 0; i < (int)superArcList_.size(); i++) {
1079 if(!superArcList_[i].pruned_) {
1080 o << " ";
1081 exportArcPosToVtk(i, pointId, arcIds[i], origin, voxelSize, o);
1082 o << endl;
1083 pointId += 8;
1084 }
1085 }
1086 o << " </DataArray>" << endl;
1087 o << " </Points>" << endl;
1088
1089 // cells now
1090 o << " <Cells>" << endl;
1091
1092 o << " <DataArray type=\"Int32\" Name=\"connectivity\""
1093 << " format=\"ascii\">" << endl;
1094
1095 for(int i = 0; i < (int)superArcList_.size(); i++) {
1096 if(!superArcList_[i].pruned_) {
1097 int downNodeId = superArcList_[i].downNodeId_;
1098 int upNodeId = superArcList_[i].upNodeId_;
1099
1100 if(!nodeMeshOut[downNodeId]) {
1101 o << " " << nodeIds[downNodeId][0] << " "
1102 << nodeIds[downNodeId][1] << " " << nodeIds[downNodeId][2] << " "
1103 << nodeIds[downNodeId][3] << " " << nodeIds[downNodeId][4] << " "
1104 << nodeIds[downNodeId][5] << " " << nodeIds[downNodeId][6] << " "
1105 << nodeIds[downNodeId][7] << endl;
1106 nodeMeshOut[downNodeId] = true;
1107 }
1108 if(!nodeMeshOut[upNodeId]) {
1109 o << " " << nodeIds[upNodeId][0] << " " << nodeIds[upNodeId][1]
1110 << " " << nodeIds[upNodeId][2] << " " << nodeIds[upNodeId][3] << " "
1111 << nodeIds[upNodeId][4] << " " << nodeIds[upNodeId][5] << " "
1112 << nodeIds[upNodeId][6] << " " << nodeIds[upNodeId][7] << endl;
1113 nodeMeshOut[upNodeId] = true;
1114 }
1115 }
1116 }
1117 for(int i = 0; i < (int)superArcList_.size(); i++) {
1118 if(!superArcList_[i].pruned_) {
1119 o << " " << arcIds[i][0] << " " << arcIds[i][1] << " "
1120 << arcIds[i][2] << " " << arcIds[i][3] << " " << arcIds[i][4] << " "
1121 << arcIds[i][5] << " " << arcIds[i][6] << " " << arcIds[i][7] << endl;
1122 }
1123 }
1124
1125 o << " </DataArray>" << endl;
1126
1127 o << " <DataArray type=\"Int32\" Name=\"offsets\""
1128 << " format=\"ascii\">" << endl;
1129 pointId = 8;
1130 for(int i = 0; i < nodeNumber; i++) {
1131 o << " " << pointId << endl;
1132 pointId += 8;
1133 }
1134 for(int i = 0; i < superArcNumber; i++) {
1135 o << " " << pointId << endl;
1136 pointId += 8;
1137 }
1138 o << " </DataArray>" << endl;
1139
1140 o << " <DataArray type=\"Int32\" Name=\"types\""
1141 << " format=\"ascii\">" << endl;
1142 for(int i = 0; i < nodeNumber; i++) {
1143 o << " 12" << endl;
1144 }
1145 for(int i = 0; i < superArcNumber; i++) {
1146 o << " 12" << endl;
1147 }
1148 o << " </DataArray>" << endl;
1149
1150 o << " </Cells>" << endl;
1151 o << " </Piece>" << endl;
1152 o << " </UnstructuredGrid>" << endl;
1153 o << "</VTKFile>" << endl;
1154 o.close();
1155
1156 return 0;
1157}
1158
1160
1161 nodeList_.clear();
1162 arcList_.clear();
1163 superArcList_.clear();
1164 vertex2node_.clear();
1165 vertex2superArc_.clear();
1166 vertex2superArcNode_.clear();
1167
1168 vertex2node_.resize(vertexNumber_, -1);
1169 vertex2superArc_.resize(vertexNumber_, -1);
1171
1172 return 0;
1173}
1174
1175int SubLevelSetTree::buildExtremumList(vector<int> &extremumList,
1176 const bool &isSubLevelSet) {
1177
1178 // TODO uncomment after triangulation templatization
1179 if((!triangulation_) /*|| triangulation_->isEmpty()*/)
1180 return -1;
1181
1182 if((!vertexScalars_) || (!vertexScalars_->size()))
1183 return -2;
1184
1185 if((!vertexSoSoffsets_) || (!vertexSoSoffsets_->size()))
1186 return -3;
1187
1188 vector<pair<bool, pair<double, pair<int, int>>>> tmpList;
1189
1190 for(SimplexId i = 0; i < triangulation_->getNumberOfVertices(); i++) {
1191
1192 bool isExtremum = true;
1194 for(SimplexId j = 0; j < neighborNumber; j++) {
1195 SimplexId otherId;
1196 triangulation_->getVertexNeighbor(i, j, otherId);
1197
1198 if(isSubLevelSet) {
1199 // looking for minima
1200 if(isSosHigherThan(i, otherId)) {
1201 // not a minimum then
1202 isExtremum = false;
1203 break;
1204 }
1205 } else {
1206 // looking for maxima
1207 if(!isSosHigherThan(i, otherId)) {
1208 // not a maximum then
1209 isExtremum = false;
1210 break;
1211 }
1212 }
1213 }
1214 if(isExtremum) {
1215 pair<bool, pair<double, pair<int, int>>> entry;
1216 entry.first = isSubLevelSet;
1217 entry.second.first = (*vertexScalars_)[i];
1218 entry.second.second.first = (*vertexSoSoffsets_)[i];
1219 entry.second.second.second = i;
1220 tmpList.push_back(entry);
1221 }
1222 }
1223
1224 sort(tmpList.begin(), tmpList.end(), filtrationCmp);
1225
1226 extremumList.resize(tmpList.size());
1227 for(int i = 0; i < (int)extremumList.size(); i++) {
1228 extremumList[i] = tmpList[i].second.second.second;
1229 }
1230
1231 if(isSubLevelSet) {
1232 minimumList_ = &(extremumList);
1233 } else {
1234 maximumList_ = &(extremumList);
1235 }
1236
1237 return 0;
1238}
1239
1240int SubLevelSetTree::closeSuperArc(const int &superArcId, const int &nodeId) {
1241
1242 if((superArcId < 0) || (superArcId >= (int)superArcList_.size()))
1243 return -1;
1244
1245 if((nodeId < 0) || (nodeId >= (int)nodeList_.size()))
1246 return -2;
1247
1248 superArcList_[superArcId].setUpNodeId(nodeId);
1249
1250 nodeList_[nodeId].addDownSuperArcId(superArcId);
1251
1252 // create a low level arc with the previous node
1253 int nId = -1;
1254 if(!superArcList_[superArcId].getNumberOfRegularNodes()) {
1255 // the down neighbor is the arc extremity
1256 nId = superArcList_[superArcId].getDownNodeId();
1257 } else {
1258 nId = superArcList_[superArcId].getRegularNodeId(
1259 superArcList_[superArcId].getNumberOfRegularNodes() - 1);
1260 }
1261
1262 makeArc(nId, nodeId);
1263
1264 return 0;
1265}
1266
1268 vector<pair<double, double>> &diagram,
1269 vector<pair<pair<int, int>, double>> *pairs) const {
1270
1271 vector<pair<pair<int, int>, double>> defaultPairs{};
1272 if(!pairs) {
1273 pairs = &defaultPairs;
1274 }
1275
1276 if(!pairs->size())
1277 getPersistencePairs(*pairs);
1278
1279 // fast fix :(
1280 diagram.resize(pairs->size());
1281
1282 for(int i = 0; i < (int)pairs->size(); i++) {
1283 if((maximumList_) && (!minimumList_)) {
1284 // split tree
1285 diagram[i].second = (*vertexScalars_)[(*pairs)[i].first.first];
1286 diagram[i].first = (*vertexScalars_)[(*pairs)[i].first.second];
1287 } else {
1288 // join tree or contour tree
1289 diagram[i].first = (*vertexScalars_)[(*pairs)[i].first.first];
1290 diagram[i].second = (*vertexScalars_)[(*pairs)[i].first.second];
1291 }
1292 }
1293
1294 std::sort(diagram.begin(), diagram.end(), _pPairCmp);
1295
1296 return 0;
1297}
1298
1300 vector<pair<pair<int, int>, double>> &pairs,
1301 std::vector<std::pair<std::pair<int, int>, double>> *ttkNotUsed(mergePairs),
1302 std::vector<std::pair<std::pair<int, int>, double>> *ttkNotUsed(
1303 splitPairs)) const {
1304
1305 Timer t;
1306
1307 if(!superArcList_.size())
1308 return -1;
1309
1310 bool isMergeTree = true;
1311
1312 vector<int> *extremumList = minimumList_;
1313 if(!extremumList) {
1314 extremumList = maximumList_;
1315 isMergeTree = false;
1316 }
1317 if(!extremumList)
1318 return -2;
1319
1320 pairs.resize(extremumList->size());
1321
1322 vector<int> min2arc(vertexNumber_, -1);
1323
1324 // compute the initial persistence
1325 int leafId = 0;
1326 for(int i = 0; i < getNumberOfSuperArcs(); ++i) {
1327 const SuperArc *a = getSuperArc(i);
1328 const Node *n = getNode(a->getDownNodeId());
1329 int nodeId = n->getVertexId();
1330
1331 if(!n->getNumberOfDownSuperArcs()) {
1332 // internal, split trees and merges trees have their leaves at the bottom
1333 double extremumScalar = 0;
1334 double saddleScalar = 0;
1335 int saddleId = getNode(a->getUpNodeId())->getVertexId();
1336
1337 extremumScalar = (*vertexScalars_)[nodeId];
1338 saddleScalar = (*vertexScalars_)[saddleId];
1339
1340 pairs[leafId].first.first = nodeId;
1341 pairs[leafId].first.second = saddleId;
1342 pairs[leafId].second = fabs(saddleScalar - extremumScalar);
1343 float persistence = pairs[leafId].second;
1344 if(isnan(persistence))
1345 pairs[leafId].second = 0;
1346
1347 // save that in the map
1348 min2arc[nodeId] = i;
1349 leafId++;
1350 }
1351 }
1352
1353 // sort the pairs and start the processing
1354 sort(pairs.begin(), pairs.end(), _pCmp);
1355
1356 vector<int> pairedSaddles(vertexNumber_, 0);
1357 // now start the pairing...
1358 for(unsigned int i = 0; i < pairs.size(); ++i) {
1359 pairedSaddles[pairs[i].first.second]++;
1360
1361 // now update the persistence of all others
1362 for(unsigned int j = i + 1; j < pairs.size(); ++j) {
1363 if((pairs[j].first.second == pairs[i].first.second)
1364 && (pairedSaddles[pairs[i].first.second] + 1
1365 == (getVertexNode(pairs[i].first.second))
1366 ->getNumberOfDownSuperArcs())) {
1367 // that pair contains a minimum competing with our guy
1368 // and the saddle is fully paired
1369
1370 const SuperArc *a = getSuperArc(min2arc[pairs[j].first.first]);
1371 do {
1373 break;
1374
1375 int nextArcId = getNode(a->getUpNodeId())->getUpSuperArcId(0);
1376 a = getSuperArc(nextArcId);
1377
1378 if(a) {
1379 int saddleId = getNode(a->getUpNodeId())->getVertexId();
1381 break;
1382
1383 if(pairedSaddles[saddleId]
1384 < (getVertexNode(saddleId))->getNumberOfDownSuperArcs() - 1) {
1385
1386 double extremumScalar = 0;
1387 double saddleScalar = 0;
1388
1389 extremumScalar = (*vertexScalars_)[pairs[j].first.first];
1390 saddleScalar = (*vertexScalars_)[saddleId];
1391 pairs[j].first.second = saddleId;
1392 pairs[j].second = fabs(saddleScalar - extremumScalar);
1393 float persistence = pairs[j].second;
1394 if(isnan(persistence))
1395 pairs[j].second = 0;
1396 break;
1397 }
1398 }
1399 } while(a);
1400
1401 if(a) {
1402 // we indeed update the persistence of the competitor
1403 sort(pairs.begin() + j, pairs.end(), _pCmp);
1404 }
1405
1406 break;
1407 }
1408 }
1409 }
1410
1411 // the last entry should have infinite persistence
1412 // pairs.erase(pairs.end() - 1);
1413
1414 // let's pair the global maximum with the global min
1415 for(int i = 0; i < getNumberOfSuperArcs(); i++) {
1416 const SuperArc *a = getSuperArc(i);
1417 const Node *up = getNode(a->getUpNodeId());
1418 if(!up->getNumberOfUpSuperArcs()) {
1419 // global min
1420 pairs.back().first.second = up->getVertexId();
1421 pairs.back().second
1422 = fabs((*vertexScalars_)[pairs.back().first.first]
1423 - (*vertexScalars_)[pairs.back().first.second]);
1424 float persistence = pairs.back().second;
1425 if(isnan(persistence))
1426 pairs.back().second = 0;
1427 break;
1428 }
1429 }
1430
1431 std::string pairtype = isMergeTree ? "(0-1)" : "(1-2)";
1432 std::string pairextr = isMergeTree ? "min" : "max";
1433
1434 this->printMsg(
1435 std::vector<std::vector<std::string>>{
1436 {"#" + pairtype + " pairs", std::to_string(pairs.size())}},
1438 for(const auto &p : pairs) {
1439 this->printMsg(
1440 std::vector<std::vector<std::string>>{
1441 {"#" + pairextr, std::to_string(p.first.first)},
1442 {"#saddles", std::to_string(p.first.second)},
1443 {"Persistence", std::to_string(p.second)}},
1445 }
1446
1447 this->printMsg(std::to_string(pairs.size()) + " persistence pairs computed",
1448 1.0, t.getElapsedTime(), this->threadNumber_);
1449
1450 return 0;
1451}
1452
1454 vector<pair<double, int>> &plot,
1455 vector<pair<pair<int, int>, double>> *persistencePairs) const {
1456
1457 vector<pair<pair<int, int>, double>> defaultPersistencePairs{};
1458 if(!persistencePairs) {
1459 persistencePairs = &defaultPersistencePairs;
1460 }
1461
1462 if(!persistencePairs->size())
1463 getPersistencePairs(*persistencePairs);
1464
1465 plot.resize(persistencePairs->size());
1466
1467 for(int i = 0; i < (int)plot.size(); i++) {
1468 plot[i].first = (*persistencePairs)[i].second;
1469 if(plot[i].first < Geometry::powIntTen(-REAL_SIGNIFICANT_DIGITS)) {
1471 }
1472 plot[i].second = persistencePairs->size() - i;
1473 }
1474
1475 return 0;
1476}
1477
1478bool SubLevelSetTree::buildPlanarLayout(const double &scaleX,
1479 const double &scaleY) {
1480
1481 if((minimumList_) && (maximumList_)) {
1482 this->printErr("Contour tree planar layout not implemented.");
1483 this->printErr("Planar layout is only implemented for merge-trees.");
1484 return false;
1485 }
1486
1487 if((!vertexScalars_) || (!vertexScalars_->size()))
1488 return false;
1489
1490 // global dimansions [0,1]x[0,1]
1491
1492 // 1) set the y coordinate for everyone
1493 for(int i = 0; i < (int)superArcList_.size(); i++) {
1494 int downNodeId = superArcList_[i].downNodeId_;
1495 int upNodeId = superArcList_[i].upNodeId_;
1496
1497 nodeList_[downNodeId].layoutY_
1498 = ((*vertexScalars_)[nodeList_[downNodeId].vertexId_] - minScalar_)
1499 / (maxScalar_ - minScalar_);
1500
1501 nodeList_[upNodeId].layoutY_
1502 = ((*vertexScalars_)[nodeList_[upNodeId].vertexId_] - minScalar_)
1503 / (maxScalar_ - minScalar_);
1504 }
1505
1506 // 2) for each node count the number of leaves under it
1507 // --> start a front from the leaves and increment a table
1508 vector<bool> inFront(nodeList_.size(), false);
1509 vector<int> node2LeafNumber(nodeList_.size(), 0);
1510 // use the same data-structure for front propagation
1511 set<pair<bool, pair<double, pair<int, int>>>, filtrationCtCmp> front;
1512
1513 // 2a) start from the leaves
1514 for(int i = 0; i < (int)superArcList_.size(); i++) {
1515
1516 if(!superArcList_[i].pruned_) {
1517
1518 int downNodeId = superArcList_[i].downNodeId_;
1519 if(!nodeList_[downNodeId].downSuperArcList_.size()) {
1520 pair<bool, pair<double, pair<int, int>>> n;
1521 n.first = (minimumList_ ? true : false);
1522 n.second.first = (*vertexScalars_)[nodeList_[downNodeId].vertexId_];
1523 n.second.second.first
1524 = (*vertexSoSoffsets_)[nodeList_[downNodeId].vertexId_];
1525 n.second.second.second = nodeList_[downNodeId].vertexId_;
1526 front.insert(n);
1527
1528 node2LeafNumber[downNodeId] = 1;
1529 inFront[downNodeId] = true;
1530 }
1531 }
1532 }
1533
1534 // filtration loop
1535 do {
1536
1537 int vertexId = front.begin()->second.second.second;
1538 front.erase(front.begin());
1539
1540 int nodeId = vertex2node_[vertexId];
1541
1542 // retrieve the number of leaves given by its children
1543 for(int i = 0; i < (int)nodeList_[nodeId].downSuperArcList_.size(); i++) {
1544 int arcId = nodeList_[nodeId].downSuperArcList_[i];
1545 int childId = superArcList_[arcId].downNodeId_;
1546
1547 node2LeafNumber[nodeId] += node2LeafNumber[childId];
1548 }
1549
1550 // add its parents to the front
1551 for(int i = 0; i < (int)nodeList_[nodeId].upSuperArcList_.size(); i++) {
1552 int arcId = nodeList_[nodeId].upSuperArcList_[i];
1553 int parentId = superArcList_[arcId].upNodeId_;
1554
1555 if(!inFront[parentId]) {
1556 pair<bool, pair<double, pair<int, int>>> n;
1557 n.first = (minimumList_ ? true : false);
1558 n.second.first = (*vertexScalars_)[nodeList_[parentId].vertexId_];
1559 n.second.second.first
1560 = (*vertexSoSoffsets_)[nodeList_[parentId].vertexId_];
1561 n.second.second.second = nodeList_[parentId].vertexId_;
1562 front.insert(n);
1563
1564 node2LeafNumber[parentId] = 0;
1565 inFront[parentId] = true;
1566 }
1567 }
1568
1569 } while(!front.empty());
1570
1571 // 3) start from the node below the root, these two guys have the same x
1572 // each node has its X coordinate given by a ratio
1573
1574 // let's find the root
1575 int rootId = -1;
1576 int subRootId = -1;
1577
1578 for(int i = 0; i < (int)superArcList_.size(); i++) {
1579 if(!superArcList_[i].pruned_) {
1580 int upNodeId = superArcList_[i].upNodeId_;
1581 if(!nodeList_[upNodeId].upSuperArcList_.size()) {
1582 // that's it
1583 rootId = upNodeId;
1584 break;
1585 }
1586 }
1587 }
1588
1589 // in theory, the root has exactly one arc below it
1590 subRootId = superArcList_[nodeList_[rootId].downSuperArcList_[0]].downNodeId_;
1591
1592 // now start the recursion
1593 vector<pair<double, double>> nodeInterval(
1594 nodeList_.size(), pair<double, double>(0, 1));
1595 queue<int> nodeQueue;
1596
1597 // printf("pushing root: %f\n",
1598 // (*vertexScalars_)[nodeList_[rootId].vertexId_]);
1599 nodeQueue.push(rootId);
1600
1601 do {
1602
1603 int nodeId = nodeQueue.front();
1604 nodeQueue.pop();
1605
1606 // get the number of leaves of the first child
1607 double ratio = 0.5;
1608
1609 if((nodeList_[nodeId].downSuperArcList_.size()) && (nodeId != subRootId)
1610 && (nodeId != rootId)) {
1611
1612 int arcId = nodeList_[nodeId].downSuperArcList_[0];
1613 int childId = superArcList_[arcId].downNodeId_;
1614
1615 ratio = node2LeafNumber[childId] / ((double)node2LeafNumber[nodeId]);
1616 if(ratio > 1)
1617 ratio = 1;
1618 }
1619
1620 // printf("node %f: [%f-%f] r=%f l=%d\n",
1621 // (*vertexScalars_)[nodeList_[nodeId].vertexId_],
1622 // nodeInterval[nodeId].first, nodeInterval[nodeId].second, ratio,
1623 // node2LeafNumber[nodeId]);
1624 nodeList_[nodeId].layoutX_
1625 = nodeInterval[nodeId].first
1626 + ratio * (nodeInterval[nodeId].second - nodeInterval[nodeId].first);
1627 if(nodeList_[nodeId].layoutX_ > 1) {
1628 nodeList_[nodeId].layoutX_ = 1;
1629 }
1630
1631 // pass the right interval to the children and add them to the queue
1632 double xOffset = 0;
1633 for(int i = 0; i < (int)nodeList_[nodeId].downSuperArcList_.size(); i++) {
1634 int arcId = nodeList_[nodeId].downSuperArcList_[i];
1635 int childId = superArcList_[arcId].downNodeId_;
1636
1637 ratio = node2LeafNumber[childId] / ((double)node2LeafNumber[nodeId]);
1638 if(ratio > 1)
1639 ratio = 1;
1640
1641 nodeInterval[childId].first = nodeInterval[nodeId].first + xOffset;
1642 nodeInterval[childId].second
1643 = nodeInterval[childId].first
1644 + ratio * (nodeInterval[nodeId].second - nodeInterval[nodeId].first);
1645 xOffset = nodeInterval[childId].second;
1646
1647 nodeQueue.push(childId);
1648 }
1649
1650 } while(!nodeQueue.empty());
1651
1652 // test
1653 // for(int i = 0; i < (int) superArcList_.size(); i++){
1654 // if(!superArcList_[i].pruned_){
1655 // Node *downNode = &(nodeList_[superArcList_[i].getDownNodeId()]);
1656 // printf("node %f --> %fx%f\n",
1657 // (*vertexScalars_)[downNode->vertexId_],
1658 // downNode->layoutX_, downNode->layoutY_);
1659 // }
1660 // }
1661 //
1662
1663 for(int i = 0; i < (int)nodeList_.size(); i++) {
1664 nodeList_[i].layoutX_ *= 30 * scaleX;
1665 nodeList_[i].layoutY_ *= 1000 * scaleY;
1666 }
1667
1668 return true;
1669}
1670
1671int SubLevelSetTree::buildSaddleList(vector<int> &vertexList) const {
1672
1673 vertexList.clear();
1674
1675 for(int i = 0; i < (int)superArcList_.size(); i++) {
1676 const Node *down = &(nodeList_[superArcList_[i].getDownNodeId()]);
1677
1678 if(down->getNumberOfDownSuperArcs()) {
1679 vertexList.push_back(down->getVertexId());
1680 }
1681 }
1682
1683 return 0;
1684}
1685
1686int SubLevelSetTree::makeArc(const int &nodeId0, const int &nodeId1) {
1687
1688 if((nodeId0 < 0) || (nodeId0 >= (int)nodeList_.size()))
1689 return -1;
1690 if((nodeId1 < 0) || (nodeId1 >= (int)nodeList_.size()))
1691 return -2;
1692
1693 arcList_.resize(arcList_.size() + 1);
1694 arcList_[arcList_.size() - 1].setDownNodeId(nodeId0);
1695 arcList_[arcList_.size() - 1].setUpNodeId(nodeId1);
1696
1697 nodeList_[nodeId0].addUpArcId((int)arcList_.size() - 1);
1698 nodeList_[nodeId1].addDownArcId((int)arcList_.size() - 1);
1699
1700 return (int)arcList_.size() - 1;
1701}
1702
1703int SubLevelSetTree::makeNode(const int &vertexId) {
1704
1705 if((vertexId < 0) || (vertexId >= vertexNumber_))
1706 return -1;
1707
1708 if(vertex2node_[vertexId] == -1) {
1709 nodeList_.resize(nodeList_.size() + 1);
1710 nodeList_[nodeList_.size() - 1].setVertexId(vertexId);
1711
1712 vertex2node_[vertexId] = (int)nodeList_.size() - 1;
1713
1714 return (int)nodeList_.size() - 1;
1715 }
1716
1717 return vertex2node_[vertexId];
1718}
1719
1720int SubLevelSetTree::openSuperArc(const int &nodeId) {
1721
1722 if((nodeId < 0) || (nodeId >= (int)nodeList_.size()))
1723 return -1;
1724
1725 superArcList_.resize(superArcList_.size() + 1);
1726 superArcList_[superArcList_.size() - 1].setDownNodeId(nodeId);
1727
1728 nodeList_[nodeId].addUpSuperArcId((int)superArcList_.size() - 1);
1729
1730 return (int)superArcList_.size() - 1;
1731}
1732
1733bool SubLevelSetTree::isSosHigherThan(const int &vertexId0,
1734 const int &vertexId1) const {
1735
1736 return (((*vertexScalars_)[vertexId0] > (*vertexScalars_)[vertexId1])
1737 || (((*vertexScalars_)[vertexId0] == (*vertexScalars_)[vertexId1])
1738 && ((*vertexSoSoffsets_)[vertexId0]
1739 > (*vertexSoSoffsets_)[vertexId1])));
1740}
1741
1742bool SubLevelSetTree::isSosLowerThan(const int &vertexId0,
1743 const int &vertexId1) const {
1744
1745 return (((*vertexScalars_)[vertexId0] < (*vertexScalars_)[vertexId1])
1746 || (((*vertexScalars_)[vertexId0] == (*vertexScalars_)[vertexId1])
1747 && ((*vertexSoSoffsets_)[vertexId0]
1748 < (*vertexSoSoffsets_)[vertexId1])));
1749}
1750
1752 const Node *oldDown,
1753 const Node *oldUp,
1754 const Node *newDown,
1755 const Node *newUp) {
1756
1757 int arcId = -1;
1758 int nodeId = n - &(nodeList_[0]);
1759
1760 // remove n from oldDown and oldUp and update n
1761 for(int i = 0; i < oldUp->getNumberOfDownArcs(); i++) {
1762 arcId = oldUp->getDownArcId(i);
1763 if(arcList_[arcId].getDownNodeId() == nodeId) {
1764 arcList_[arcId].setDownNodeId(oldDown - &(nodeList_[0]));
1765 break;
1766 }
1767 }
1768 for(int i = 0; i < n->getNumberOfUpArcs(); i++) {
1769 if(arcId == n->getUpArcId(i)) {
1770 nodeList_[nodeId].removeUpArcId(i);
1771 break;
1772 }
1773 }
1774
1775 for(int i = 0; i < oldDown->getNumberOfUpArcs(); i++) {
1776 if(arcList_[oldDown->getUpArcId(i)].getUpNodeId() == nodeId) {
1777 nodeList_[oldDown - &(nodeList_[0])].removeUpArcId(i);
1778 nodeList_[oldDown - &(nodeList_[0])].addUpArcId(arcId);
1779 break;
1780 }
1781 }
1782 for(int i = 0; i < n->getNumberOfDownArcs(); i++) {
1783 if(arcList_[n->getDownArcId(i)].getDownNodeId()
1784 == (oldDown - &(nodeList_[0]))) {
1785 nodeList_[nodeId].removeDownArcId(i);
1786 break;
1787 }
1788 }
1789
1790 // disconnect newDown and newUp
1791 if(newUp) {
1792 for(int i = 0; i < newDown->getNumberOfUpArcs(); i++) {
1793 arcId = newDown->getUpArcId(i);
1794 if(arcList_[arcId].getUpNodeId() == (newUp - &(nodeList_[0]))) {
1795 nodeList_[newDown - &(nodeList_[0])].removeUpArcId(i);
1796 break;
1797 }
1798 }
1799 for(int i = 0; i < newUp->getNumberOfDownArcs(); i++) {
1800 if(newUp->getDownArcId(i) == arcId) {
1801 nodeList_[newUp - &(nodeList_[0])].removeDownArcId(i);
1802 break;
1803 }
1804 }
1805 }
1806
1807 makeArc(newDown - &(nodeList_[0]), n - &(nodeList_[0]));
1808 if(newUp)
1809 makeArc(n - &(nodeList_[0]), newUp - &(nodeList_[0]));
1810
1811 return 0;
1812}
1813
1815
1816 stringstream msg;
1817
1818 this->printMsg(
1819 "Node list (" + std::to_string(getNumberOfNodes()) + " nodes):",
1821
1822 int minCount = 0, saddleCount = 0, maxCount = 0, regularCount = 0;
1823
1824 for(int i = 0; i < getNumberOfNodes(); i++) {
1825 const Node *n = getNode(i);
1826
1827 if(vertex2superArc_[n->getVertexId()] == -1) {
1828 // regular node
1829
1830 this->printMsg(
1831 std::vector<std::vector<std::string>>{
1832 {"Id", std::to_string(i), "VertId", std::to_string(n->getVertexId()),
1833 "D", std::to_string(n->getNumberOfDownSuperArcs()), "U",
1834 std::to_string(n->getNumberOfUpSuperArcs())}},
1836
1837 if(!n->getNumberOfDownSuperArcs())
1838 minCount++;
1839 else if(!n->getNumberOfUpSuperArcs())
1840 maxCount++;
1841 else
1842 saddleCount++;
1843 }
1844 }
1845
1846 this->printMsg(
1847 "Arc list (" + std::to_string(getNumberOfSuperArcs()) + " arcs):",
1849
1850 for(int i = 0; i < getNumberOfSuperArcs(); i++) {
1851 const SuperArc *a = getSuperArc(i);
1852
1853 const Node *down = getNode(a->getDownNodeId()),
1854 *up = getNode(a->getUpNodeId());
1855 if((up) && (down)) {
1856
1857 this->printMsg(
1858 std::vector<std::vector<std::string>>{
1859 {"Id", std::to_string(i), "D",
1860 std::to_string(down->getNumberOfDownSuperArcs()), "U",
1861 std::to_string(up->getNumberOfUpSuperArcs()), "V",
1862 std::to_string(a->getNumberOfRegularNodes())}},
1864 } else {
1865 this->printErr("Arc inconsistency! " + std::to_string(a->getDownNodeId())
1866 + "->" + std::to_string(a->getUpNodeId()));
1867 }
1868
1869 for(int j = 0; j < a->getNumberOfRegularNodes(); j++) {
1870 this->printMsg(
1871 std::vector<std::vector<std::string>>{
1872 {"Regular vertex",
1873 std::to_string(nodeList_[a->getRegularNodeId(j)].getVertexId())}},
1875 }
1876
1877 regularCount += a->getNumberOfRegularNodes();
1878 }
1879
1880 this->printMsg(
1881 std::vector<std::vector<std::string>>{
1882 {"#Minima", std::to_string(minCount)},
1883 {"#Saddles", std::to_string(saddleCount)},
1884 {"#Maxima", std::to_string(maxCount)},
1885 {"#Regular", std::to_string(regularCount)},
1886 {"#Sum",
1887 std::to_string(minCount + saddleCount + maxCount + regularCount)},
1888 {"#Input vertices", std::to_string(vertexNumber_)},
1889 },
1891
1892 return 0;
1893}
1894
1895int SubLevelSetTree::simplify(const double &simplificationThreshold,
1897
1898 Timer t;
1899
1900 // if((simplificationThreshold < 0)||(simplificationThreshold > 1))
1901 // return -1;
1902
1903 PersistenceMetric defaultMetric;
1904
1905 if(metric == nullptr) {
1906 metric = &defaultMetric;
1907 }
1908
1909 metric->tree_ = this;
1910
1911 // double unNormalizedThreshold = simplificationThreshold
1912 // * (maxScalar_ - minScalar_);
1913
1914 if(!nodeList_.size())
1915 return -2;
1916
1917 if(!superArcList_.size())
1918 return -3;
1919
1920 if((minimumList_) && (maximumList_)) {
1921 this->printErr("Contour tree simplification not implemented.");
1922 this->printErr("Simplification is only implemented for merge-trees.");
1923 return -4;
1924 }
1925
1926 int simplifiedArcNumber = 0;
1927 double maximumMetricScore = 0;
1928
1929 if(!originalNodeList_.size()) {
1930 // first time we simplify
1933 } else {
1934 // not the first time
1937 }
1938
1939 if(!simplificationThreshold)
1940 return 0;
1941
1942 // if(!unNormalizedThreshold) return -4;
1943
1944 vector<pair<real, int>> regularNodeList;
1945
1946 int arcToSimplify = -1;
1947 double currentScore = 0, minScore;
1948
1949 do {
1950
1951 arcToSimplify = -1;
1952 minScore = 1.1 * (maxScalar_ - minScalar_);
1953
1954 // search for an elligible arc.
1955 for(int i = 0; i < (int)superArcList_.size(); i++) {
1956
1957 if(!superArcList_[i].pruned_) {
1958
1959 int downNodeId = superArcList_[i].downNodeId_;
1960 int upNodeId = superArcList_[i].upNodeId_;
1961
1962 if(!nodeList_[downNodeId].downSuperArcList_.size()) {
1963
1964 currentScore = metric->computeSuperArcMetric(
1965 nodeList_[downNodeId].vertexId_, nodeList_[upNodeId].vertexId_,
1966 superArcList_[i].regularNodeList_);
1967
1968 if((currentScore < 0) || (currentScore > (maxScalar_ - minScalar_))) {
1969 this->printWrn("Out-of-range score! (user-defined metric)");
1970 } else {
1971 if(currentScore < minScore) {
1972 arcToSimplify = i;
1973 minScore = currentScore;
1974 }
1975 }
1976 }
1977 }
1978 }
1979
1980 if(minScore > simplificationThreshold)
1981 break;
1982
1983 if(minScore > maximumMetricScore)
1984 maximumMetricScore = minScore;
1985
1986 if(arcToSimplify == -1)
1987 break;
1988
1989 // perform the simplification
1990 int downNodeId = superArcList_[arcToSimplify].downNodeId_;
1991 int upNodeId = superArcList_[arcToSimplify].upNodeId_;
1992
1993 int pivotId = upNodeId;
1994 int leafId = downNodeId;
1995
1996 if(!nodeList_[downNodeId].downSuperArcList_.size()) {
1997 // remove an extremum
1998 pivotId = upNodeId;
1999 leafId = downNodeId;
2000 }
2001
2002 bool isSimpleSaddle = (nodeList_[pivotId].downSuperArcList_.size()
2003 + nodeList_[pivotId].upSuperArcList_.size()
2004 == 3);
2005
2006 int brotherId = 0;
2007
2008 for(int i = 0; i < (int)nodeList_[pivotId].downSuperArcList_.size(); i++) {
2009 if(nodeList_[pivotId].downSuperArcList_[i] != arcToSimplify) {
2010 brotherId = nodeList_[pivotId].downSuperArcList_[i];
2011 break;
2012 }
2013 }
2014
2015 int brotherExtremityId = superArcList_[brotherId].downNodeId_;
2016
2017 int parentId = nodeList_[pivotId].upSuperArcList_[0];
2018
2019 // 1) update the parent extremity (only if simple saddles)
2020 if(isSimpleSaddle) {
2021 superArcList_[parentId].downNodeId_ = brotherExtremityId;
2022 }
2023
2024 // 2) tell the brotherExtremityId that it's no longer linked to the
2025 // brother but to the parent
2026 if(isSimpleSaddle) {
2027 for(int i = 0;
2028 i < (int)nodeList_[brotherExtremityId].upSuperArcList_.size(); i++) {
2029 if(nodeList_[brotherExtremityId].upSuperArcList_[i] == brotherId) {
2030 nodeList_[brotherExtremityId].upSuperArcList_[i] = parentId;
2031 break;
2032 }
2033 }
2034 } else {
2035 // multi-saddle, we need to erase arcToSimplify from the pivot
2036 for(int i = 0; i < (int)nodeList_[pivotId].downSuperArcList_.size();
2037 i++) {
2038 if(nodeList_[pivotId].downSuperArcList_[i] == arcToSimplify) {
2039 nodeList_[pivotId].downSuperArcList_.erase(
2040 nodeList_[pivotId].downSuperArcList_.begin() + i);
2041 break;
2042 }
2043 }
2044 }
2045
2046 // 3) do the merging
2047 if(isSimpleSaddle) {
2048 // copy the regular nodes of the brother into the parent
2049 superArcList_[parentId].regularNodeList_.insert(
2050 superArcList_[parentId].regularNodeList_.end(),
2051 superArcList_[brotherId].regularNodeList_.begin(),
2052 superArcList_[brotherId].regularNodeList_.end());
2053
2054 // copy the regular nodes of the arcToSimplify into the parent
2055 superArcList_[parentId].regularNodeList_.insert(
2056 superArcList_[parentId].regularNodeList_.end(),
2057 superArcList_[arcToSimplify].regularNodeList_.begin(),
2058 superArcList_[arcToSimplify].regularNodeList_.end());
2059
2060 // add the pivot
2061 superArcList_[parentId].regularNodeList_.push_back(pivotId);
2062
2063 // add the leaf
2064 superArcList_[parentId].regularNodeList_.push_back(leafId);
2065
2066 // NOTE: optionally, the list of regular nodes could be sorted
2067
2068 // mark the brother as pruned
2069 superArcList_[brotherId].pruned_ = true;
2070 nodeList_[pivotId].pruned_ = true;
2071 simplifiedArcNumber++;
2072 } else {
2073 // just merge with the first brother we found....
2074 superArcList_[brotherId].regularNodeList_.insert(
2075 superArcList_[brotherId].regularNodeList_.end(),
2076 superArcList_[arcToSimplify].regularNodeList_.begin(),
2077 superArcList_[arcToSimplify].regularNodeList_.end());
2078
2079 // add the leaf
2080 superArcList_[brotherId].regularNodeList_.push_back(leafId);
2081 }
2082 superArcList_[arcToSimplify].pruned_ = true;
2083 nodeList_[leafId].pruned_ = true;
2084 simplifiedArcNumber++;
2085
2086 } while(minScore < simplificationThreshold);
2087
2088 // sort the regular nodes
2089 Timer sortTimer;
2090 for(int i = 0; i < getNumberOfSuperArcs(); ++i) {
2091 if(!superArcList_[i].pruned_) {
2092 superArcList_[i].sortRegularNodes(
2094 }
2095 }
2096
2097 this->printMsg("Regular node sorting", 1.0, sortTimer.getElapsedTime(),
2098 this->threadNumber_);
2099
2100 this->printMsg(std::to_string(simplifiedArcNumber)
2101 + " super-arcs pruned (threshold="
2102 + std::to_string(simplificationThreshold) + ")",
2103 1.0, t.getElapsedTime(), this->threadNumber_);
2104
2105 this->printMsg(
2106 "Biggest simplification metric score: " + std::to_string(maximumMetricScore)
2107 + " ("
2108 + std::to_string(maximumMetricScore * 100 / (maxScalar_ - minScalar_))
2109 + "%)");
2110
2111 return 0;
2112}
2113
2115 vector<int> sample;
2116 vector<double> barycenter(3);
2117 int vertexId;
2118
2119 const SuperArc *a;
2120 for(int i = 0; i < getNumberOfSuperArcs(); ++i) {
2121 a = getSuperArc(i);
2122 if(!a->isPruned()) {
2123 for(int j = 0; j < a->getNumberOfSamples(); ++j) {
2124 a->getSample(j, sample);
2125
2126 for(unsigned int k = 0; k < 3; ++k)
2127 barycenter[k] = 0;
2128
2129 for(unsigned int k = 0; k < sample.size(); ++k) {
2130 vertexId = getNode(sample[k])->getVertexId();
2131
2132 for(unsigned int l = 0; l < 3; ++l)
2133 barycenter[l] += (*vertexPositions_)[vertexId][l];
2134 }
2135 if(sample.size()) {
2136 for(unsigned int k = 0; k < 3; ++k)
2137 barycenter[k] /= sample.size();
2138
2139 // update the arc
2140 superArcList_[i].appendBarycenter(barycenter);
2141 }
2142 }
2143 }
2144 }
2145
2146 return 0;
2147}
2148
2150 const vector<double> &scalars,
2151 vector<vector<double>> &skeletonScalars) const {
2152 skeletonScalars.clear();
2153 skeletonScalars.resize(getNumberOfSuperArcs());
2154
2155 vector<int> sample;
2156 int nodeId;
2157 int vertexId;
2158
2159 double f;
2160 double f0;
2161 double f1;
2162 double fmin;
2163 double fmax;
2164 int nodeMinId;
2165 int nodeMaxId;
2166 int nodeMinVId;
2167 int nodeMaxVId;
2168 const SuperArc *a;
2169 for(int i = 0; i < getNumberOfSuperArcs(); ++i) {
2170 a = getSuperArc(i);
2171
2172 if(!a->isPruned()) {
2173 if(maximumList_ && !minimumList_) {
2174 nodeMinId = a->getUpNodeId();
2175 nodeMaxId = a->getDownNodeId();
2176 } else {
2177 nodeMaxId = a->getUpNodeId();
2178 nodeMinId = a->getDownNodeId();
2179 }
2180
2181 nodeMaxVId = getNode(nodeMaxId)->getVertexId();
2182 nodeMinVId = getNode(nodeMinId)->getVertexId();
2183
2184 fmax = scalars[nodeMaxVId];
2185 fmin = scalars[nodeMinVId];
2186
2187 // init: min
2188 f0 = fmin;
2189
2190 // iteration
2191 for(int j = 0; j < a->getNumberOfSamples(); ++j) {
2192 a->getSample(j, sample);
2193
2194 f = 0;
2195 for(unsigned int k = 0; k < sample.size(); ++k) {
2196 nodeId = sample[k];
2197 vertexId = getNode(nodeId)->getVertexId();
2198 f += scalars[vertexId];
2199 }
2200 if(sample.size()) {
2201 f /= sample.size();
2202
2203 f1 = f;
2204 // update the arc
2205 skeletonScalars[i].push_back((f0 + f1) / 2);
2206 f0 = f1;
2207 }
2208 }
2209
2210 // end: max
2211 f1 = fmax;
2212
2213 // update the arc
2214 skeletonScalars[i].push_back((f0 + f1) / 2);
2215 }
2216 }
2217
2218 return 0;
2219}
2220
2221int SubLevelSetTree::computeSkeleton(unsigned int arcResolution) {
2222 sample(arcResolution);
2224
2225 isSkeletonComputed_ = true;
2226 return 0;
2227}
2228
2229int SubLevelSetTree::smoothSkeleton(unsigned int skeletonSmoothing) {
2230 for(unsigned int i = 0; i < skeletonSmoothing; ++i) {
2231 for(int j = 0; j < getNumberOfSuperArcs(); ++j) {
2232 if(!superArcList_[j].isPruned()) {
2233 if(minimumList_)
2234 superArcList_[j].smooth(nodeList_, vertexPositions_, true);
2235 else
2236 superArcList_[j].smooth(nodeList_, vertexPositions_, false);
2237 }
2238 }
2239 }
2240
2241 return 0;
2242}
2243
2244int SubLevelSetTree::sample(unsigned int samplingLevel) {
2245 vector<vector<int>> sampleList(samplingLevel);
2246
2247 const SuperArc *a;
2248 for(int i = 0; i < getNumberOfSuperArcs(); ++i) {
2249 a = getSuperArc(i);
2250
2251 if(!a->isPruned()) {
2252
2253 for(unsigned int j = 0; j < samplingLevel; ++j)
2254 sampleList[j].clear();
2255
2256 double fmax, fmin;
2257 int nodeMaxId, nodeMinId;
2258 int nodeMaxVId, nodeMinVId;
2259 double delta;
2260 if(a->getNumberOfRegularNodes()) {
2261 if(minimumList_) {
2262 nodeMaxId = a->getUpNodeId();
2263 nodeMinId = a->getDownNodeId();
2264 } else {
2265 nodeMaxId = a->getDownNodeId();
2266 nodeMinId = a->getUpNodeId();
2267 }
2268
2269 nodeMaxVId = getNode(nodeMaxId)->getVertexId();
2270 nodeMinVId = getNode(nodeMinId)->getVertexId();
2271
2272 fmax = (*vertexScalars_)[nodeMaxVId];
2273 fmin = (*vertexScalars_)[nodeMinVId];
2274
2275 delta = (fmax - fmin) / samplingLevel;
2276
2277 double f;
2278 int nodeId;
2279 int vertexId;
2280 for(int j = 0; j < a->getNumberOfRegularNodes(); ++j) {
2281 nodeId = a->getRegularNodeId(j);
2282 vertexId = getNode(nodeId)->getVertexId();
2283 f = (*vertexScalars_)[vertexId];
2284
2285 for(unsigned int k = 0; k < samplingLevel; ++k) {
2286 if(f <= (k + 1) * delta + fmin) {
2287 sampleList[k].push_back(nodeId);
2288 break;
2289 }
2290 }
2291 }
2292
2293 // update the arc
2294 for(unsigned int j = 0; j < sampleList.size(); ++j)
2295 superArcList_[i].appendSample(sampleList[j]);
2296 }
2297 }
2298 }
2299
2300 return 0;
2301}
2302
2304 for(int j = 0; j < getNumberOfSuperArcs(); ++j) {
2305 superArcList_[j].clearBarycenters();
2306 superArcList_[j].clearSamples();
2307 }
2308
2309 return 0;
2310}
2311
2313 this->setDebugMsgPrefix("ContourTree");
2314}
2315
2317
2318 Timer timer;
2319
2320 if(!vertexNumber_)
2321 return -1;
2322 if((!vertexScalars_) || ((int)vertexScalars_->size() != vertexNumber_))
2323 return -2;
2325 return -3;
2326
2329
2330 std::vector<int> localVertSoSoffsets{};
2331
2332 // 0) init data structures
2333 if(!externalOffsets_) {
2334 vertexSoSoffsets_ = &localVertSoSoffsets;
2336 for(int i = 0; i < (int)vertexSoSoffsets_->size(); i++)
2337 (*vertexSoSoffsets_)[i] = i;
2338 }
2339
2340 // build the actual extrema list
2341
2342 std::vector<int> minimumVec, maximumVec;
2343 minimumList_ = &minimumVec;
2344 maximumList_ = &maximumVec;
2345
2346 for(int i = 0; i < vertexNumber_; i++) {
2347
2348 bool isMin = true, isMax = true;
2350 for(SimplexId j = 0; j < neighborNumber; j++) {
2351 SimplexId nId;
2353
2354 if(((*vertexScalars_)[nId] > (*vertexScalars_)[i])
2355 || (((*vertexScalars_)[nId] == (*vertexScalars_)[i])
2356 && ((*vertexSoSoffsets_)[nId] > (*vertexSoSoffsets_)[i])))
2357 isMax = false;
2358
2359 if(((*vertexScalars_)[nId] < (*vertexScalars_)[i])
2360 || (((*vertexScalars_)[nId] == (*vertexScalars_)[i])
2361 && ((*vertexSoSoffsets_)[nId] < (*vertexSoSoffsets_)[i])))
2362 isMin = false;
2363
2364 if((!isMin) && (!isMax))
2365 break;
2366 }
2367
2368 if((isMin) && (!isMax))
2369 minimumList_->push_back(i);
2370 if((isMax) && (!isMin))
2371 maximumList_->push_back(i);
2372 }
2373
2374#ifdef TTK_ENABLE_OPENMP
2375#pragma omp parallel sections
2376#endif
2377 {
2378
2379#ifdef TTK_ENABLE_OPENMP
2380#pragma omp section
2381#endif
2382 {
2383 // 1) build the merge tree
2390 mergeTree_.build();
2391 }
2392
2393#ifdef TTK_ENABLE_OPENMP
2394#pragma omp section
2395#endif
2396 {
2397 // 2) build the split tree
2404 splitTree_.build();
2405 }
2406 }
2407
2408 // note: at this point, the split tree is laid out upside down.
2409
2410 // 3) merge the two trees into the contour tree
2411 combineTrees();
2412
2413 // 4) update the high level super arc structure
2414 finalize();
2415
2416 this->printMsg("ContourTree computed", 1.0, timer.getElapsedTime());
2417
2418 this->printMsg(std::vector<std::vector<std::string>>{
2419 {"#Nodes", std::to_string(getNumberOfNodes())},
2420 {"#Arcs", std::to_string(getNumberOfArcs())}});
2421
2422 print();
2423
2424 return 0;
2425}
2426
2428
2429 Timer timer;
2430
2432 return -1;
2433
2434 queue<const Node *> nodeQueue;
2435 const Node *mergeNode = nullptr, *splitNode = nullptr;
2436
2437 do {
2438
2439 int initQueueSize = (int)nodeQueue.size();
2440
2441 for(int i = 0; i < mergeTree_.getNumberOfNodes(); i++) {
2442 mergeNode = mergeTree_.getNode(i);
2443 if(isNodeEligible(mergeNode)) {
2444 nodeQueue.push(mergeNode);
2445 }
2446 }
2447 for(int i = 0; i < splitTree_.getNumberOfNodes(); i++) {
2448 splitNode = splitTree_.getNode(i);
2449 if(isNodeEligible(splitNode)) {
2450 nodeQueue.push(splitNode);
2451 }
2452 }
2453
2454 // no more eligible nodes
2455 if((int)nodeQueue.size() == initQueueSize)
2456 break;
2457
2458 const Node *n0 = nullptr, *n1 = nullptr, *other = nullptr;
2459
2460 do {
2461
2462 n0 = nodeQueue.front();
2463 nodeQueue.pop();
2464
2465#ifndef TTK_ENABLE_KAMIKAZE
2466 if(n0 == nullptr) {
2467 continue;
2468 }
2469#endif // TTK_ENABLE_KAMIKAZE
2470
2471 if((mergeTree_.getNode(n0 - mergeTree_.getNode(0)) == n0)
2472 && (isNodeEligible(n0))) {
2473 // merge node
2474 // get the corresponding split node
2475 other = splitTree_.getVertexNode(n0->getVertexId());
2476
2477 if(!((other->getNumberOfUpArcs())
2478 && (other->getNumberOfDownArcs() > 1))) {
2479
2480 n1 = nullptr;
2481
2482 if(n0->getNumberOfUpArcs()) {
2483
2484 n1 = mergeTree_.getNodeUpNeighbor(n0, 0);
2485
2486 int newNode0 = makeNode(n0->getVertexId()),
2487 newNode1 = makeNode(n1->getVertexId());
2488
2489 // "we move v and its incident edge from JT to the contour tree".
2490 makeArc(newNode0, newNode1);
2491
2492 mergeTree_.clearArc(n0->getVertexId(), n1->getVertexId());
2493 }
2494 if((other->getNumberOfDownArcs() == 1)
2495 && (other->getNumberOfUpArcs() == 1)) {
2496 // "in ST, either v is a degree 2 node or..."
2497 // "in the former case, we suppress v in ST while maintaining the
2498 // connection"
2499 splitTree_.clearRegularNode(other->getVertexId());
2500 } else {
2501 // "or it's the root"
2502 // "in the latter, we delete v from ST"
2503 splitTree_.clearRoot(other->getVertexId());
2504 }
2505
2506 // update n1 if it's eligible
2507 if((n1) && (isNodeEligible(n1))) {
2508
2509 nodeQueue.push(n1);
2510 }
2511 } else {
2512 // unsure
2513 if(nodeQueue.empty())
2514 break;
2515 nodeQueue.push(n0);
2516 }
2517 }
2518
2519 // symmetric case with nodes coming from the split tree
2520 else if((splitTree_.getNode(n0 - splitTree_.getNode(0)) == n0)
2521 && (isNodeEligible(n0))) {
2522
2523 // merge node
2524 // get the corresponding split node
2525 other = mergeTree_.getVertexNode(n0->getVertexId());
2526
2527 if(!((other->getNumberOfUpArcs())
2528 && (other->getNumberOfDownArcs() > 1))) {
2529
2530 n1 = nullptr;
2531
2532 if(n0->getNumberOfUpArcs()) {
2533
2534 n1 = splitTree_.getNodeUpNeighbor(n0, 0);
2535
2536 int newNode0 = makeNode(n0->getVertexId()),
2537 newNode1 = makeNode(n1->getVertexId());
2538
2539 makeArc(newNode1, newNode0);
2540 splitTree_.clearArc(n0->getVertexId(), n1->getVertexId());
2541 }
2542 if((other->getNumberOfDownArcs() == 1)
2543 && (other->getNumberOfUpArcs() == 1)) {
2544 // degree 2 node
2545 mergeTree_.clearRegularNode(other->getVertexId());
2546 } else {
2547 // other has to be the root
2548 mergeTree_.clearRoot(other->getVertexId());
2549 }
2550
2551 // update n1 if it's eligible
2552 if((n1) && (isNodeEligible(n1))) {
2553 nodeQueue.push(n1);
2554 }
2555 } else {
2556 if(nodeQueue.empty())
2557 break;
2558 nodeQueue.push(n0);
2559 }
2560 }
2561 } while(nodeQueue.size());
2562
2563 if((int)nodeList_.size() == vertexNumber_)
2564 break;
2565
2566 // if one of the two trees became a line, we need to re-iterate the process
2567 } while((int)nodeList_.size() < vertexNumber_);
2568
2569 if((int)nodeList_.size() != vertexNumber_) {
2570 this->printErr("Incomplete contour tree! ("
2571 + std::to_string(nodeList_.size()) + " vs. "
2572 + std::to_string(vertexNumber_) + ")");
2573 }
2574
2575 this->printMsg(
2576 "MergeTree and SplitTree combined", 1.0, timer.getElapsedTime(), 1);
2577
2578 return 0;
2579}
2580
2582
2583 vector<bool> inQueue(nodeList_.size(), false);
2584 queue<int> nodeIdQueue;
2585
2586 for(int i = 0; i < (int)nodeList_.size(); i++) {
2587 if(!nodeList_[i].getNumberOfDownArcs()) {
2588 nodeIdQueue.push(i);
2589 inQueue[i] = true;
2590 }
2591 }
2592
2593 while(nodeIdQueue.size()) {
2594
2595 int nodeId = nodeIdQueue.front();
2596 nodeIdQueue.pop();
2597
2598 for(int i = 0; i < nodeList_[nodeId].getNumberOfUpArcs(); i++) {
2599 int nextNodeId = finalizeSuperArc(nodeId, i);
2600
2601 if(!inQueue[nextNodeId]) {
2602 nodeIdQueue.push(nextNodeId);
2603 inQueue[nextNodeId] = true;
2604 }
2605 }
2606 }
2607
2608 return 0;
2609}
2610
2611int ContourTree::finalizeSuperArc(const int &nodeId, const int &arcId) {
2612
2613 if((nodeId < 0) || (nodeId >= (int)nodeList_.size()))
2614 return -1;
2615 if((arcId < 0) || (arcId >= nodeList_[nodeId].getNumberOfUpArcs()))
2616 return -2;
2617
2618 int superArcId = openSuperArc(nodeId);
2619
2620 int currentNodeId = nodeId;
2621
2622 do {
2623
2624 if(nodeList_[currentNodeId].getNumberOfUpArcs()) {
2625 if(currentNodeId != nodeId) {
2626
2627 superArcList_[superArcId].appendRegularNode(currentNodeId);
2628 vertex2superArc_[nodeList_[currentNodeId].getVertexId()] = superArcId;
2629 vertex2superArcNode_[nodeList_[currentNodeId].getVertexId()]
2630 = superArcList_[superArcId].getNumberOfRegularNodes() - 1;
2631
2632 currentNodeId
2633 = arcList_[nodeList_[currentNodeId].getUpArcId(0)].getUpNodeId();
2634 } else {
2635 currentNodeId
2636 = arcList_[nodeList_[currentNodeId].getUpArcId(arcId)].getUpNodeId();
2637 }
2638 }
2639
2640 } while((currentNodeId != nodeId)
2641 && (nodeList_[currentNodeId].getNumberOfUpArcs() == 1)
2642 && (nodeList_[currentNodeId].getNumberOfDownArcs() == 1));
2643
2644 if(currentNodeId != nodeId) {
2645 superArcList_[superArcId].setUpNodeId(currentNodeId);
2646 nodeList_[currentNodeId].addDownSuperArcId(superArcId);
2647 }
2648
2649 return currentNodeId;
2650}
2651
2652bool ContourTree::isNodeEligible(const Node *n) const {
2653
2654#ifndef TTK_ENABLE_KAMIKAZE
2655 if(n == nullptr) {
2656 return false;
2657 }
2658#endif // TTK_ENABLE_KAMIKAZE
2659
2660 const Node *merge = nullptr, *split = nullptr;
2661
2662 if(mergeTree_.getNode(n - mergeTree_.getNode(0)) == n) {
2663 merge = n;
2664 split = splitTree_.getVertexNode(merge->getVertexId());
2665
2666 // "choose a non-root leaf that is not a split in ST"
2667 // non-root
2668 if((!merge->getNumberOfDownArcs())
2669 && (merge->getNumberOfUpArcs())
2670 // it is not a split in ST
2671 && (((split->getNumberOfDownArcs() < 2) && (split->getNumberOfUpArcs()))
2672 // or it's the root
2673 || ((split->getNumberOfDownArcs() == 1)
2674 && (!split->getNumberOfUpArcs())))) {
2675
2676 return true;
2677 }
2678 }
2679
2680 if(splitTree_.getNode(n - splitTree_.getNode(0)) == n) {
2681 split = n;
2682 merge = mergeTree_.getVertexNode(split->getVertexId());
2683
2684 if((!split->getNumberOfDownArcs()) && (split->getNumberOfUpArcs())
2685 && (((merge->getNumberOfDownArcs() < 2) && (merge->getNumberOfUpArcs()))
2686 || ((merge->getNumberOfDownArcs() == 1)
2687 && (!merge->getNumberOfUpArcs())))) {
2688
2689 return true;
2690 }
2691 }
2692
2693 return false;
2694}
2695
2696int ContourTree::computeSkeleton(unsigned int arcResolution) {
2697#ifdef TTK_ENABLE_OPENMP
2698#pragma omp parallel sections
2699#endif
2700 {
2701#ifdef TTK_ENABLE_OPENMP
2702#pragma omp section
2703#endif
2704 { SubLevelSetTree::computeSkeleton(arcResolution); }
2705
2706#ifdef TTK_ENABLE_OPENMP
2707#pragma omp section
2708#endif
2709 { mergeTree_.computeSkeleton(arcResolution); }
2710
2711#ifdef TTK_ENABLE_OPENMP
2712#pragma omp section
2713#endif
2714 { splitTree_.computeSkeleton(arcResolution); }
2715 }
2716
2717 return 0;
2718}
2719
2720int ContourTree::smoothSkeleton(unsigned int skeletonSmoothing) {
2721#ifdef TTK_ENABLE_OPENMP
2722#pragma omp parallel sections
2723#endif
2724 {
2725#ifdef TTK_ENABLE_OPENMP
2726#pragma omp section
2727#endif
2728 { SubLevelSetTree::smoothSkeleton(skeletonSmoothing); }
2729
2730#ifdef TTK_ENABLE_OPENMP
2731#pragma omp section
2732#endif
2733 { mergeTree_.smoothSkeleton(skeletonSmoothing); }
2734
2735#ifdef TTK_ENABLE_OPENMP
2736#pragma omp section
2737#endif
2738 { splitTree_.smoothSkeleton(skeletonSmoothing); }
2739 }
2740
2741 return 0;
2742}
2743
2745#ifdef TTK_ENABLE_OPENMP
2746#pragma omp parallel sections
2747#endif
2748 {
2749#ifdef TTK_ENABLE_OPENMP
2750#pragma omp section
2751#endif
2753
2754#ifdef TTK_ENABLE_OPENMP
2755#pragma omp section
2756#endif
2758
2759#ifdef TTK_ENABLE_OPENMP
2760#pragma omp section
2761#endif
2763 }
2764
2765 return 0;
2766}
2767
2769 vector<pair<pair<int, int>, double>> &pairs,
2770 vector<pair<pair<int, int>, double>> *mergePairs,
2771 vector<pair<pair<int, int>, double>> *splitPairs) const {
2772 if(pairs.size())
2773 return 0;
2774
2775 vector<pair<pair<int, int>, double>> defaultMergePairs{};
2776 if(!mergePairs) {
2777 mergePairs = &defaultMergePairs;
2778 }
2779 vector<pair<pair<int, int>, double>> defaultSplitPairs{};
2780 if(!splitPairs) {
2781 splitPairs = &defaultSplitPairs;
2782 }
2783
2784 if(!mergePairs->size() || !splitPairs->size()) {
2785#ifdef TTK_ENABLE_OPENMP
2786#pragma omp parallel sections
2787#endif
2788 {
2789
2790#ifdef TTK_ENABLE_OPENMP
2791#pragma omp section
2792#endif
2793 {
2794 if(!mergePairs->size())
2795 mergeTree_.getPersistencePairs(*mergePairs);
2796 }
2797
2798#ifdef TTK_ENABLE_OPENMP
2799#pragma omp section
2800#endif
2801 {
2802 if(!splitPairs->size())
2803 splitTree_.getPersistencePairs(*splitPairs);
2804 }
2805 }
2806 }
2807
2808 pairs.resize(mergePairs->size() + splitPairs->size());
2809
2810 for(unsigned int i = 0; i < mergePairs->size(); ++i)
2811 pairs[i] = (*mergePairs)[i];
2812
2813 unsigned int shift = mergePairs->size();
2814 for(unsigned int i = 0; i < splitPairs->size(); ++i)
2815 pairs[shift + i] = (*splitPairs)[i];
2816
2817 std::sort(pairs.begin(), pairs.end(), _pCmp);
2818
2819 return 0;
2820}
2821
2823 vector<pair<double, int>> &plot,
2824 vector<pair<pair<int, int>, double>> *mergePairs,
2825 vector<pair<pair<int, int>, double>> *splitPairs,
2826 vector<pair<pair<int, int>, double>> *pairs) const {
2827
2828 vector<pair<pair<int, int>, double>> defaultPairs{};
2829 if(!pairs) {
2830 pairs = &defaultPairs;
2831 }
2832
2833 if(!pairs->size())
2834 getPersistencePairs(*pairs, mergePairs, splitPairs);
2835
2836 plot.resize(pairs->size());
2837
2838 for(int i = 0; i < (int)plot.size(); i++) {
2839 plot[i].first = (*pairs)[i].second;
2840 if(plot[i].first < Geometry::powIntTen(-REAL_SIGNIFICANT_DIGITS)) {
2842 }
2843 plot[i].second = pairs->size() - i;
2844 }
2845
2846 return 0;
2847}
2848
2850 vector<pair<double, double>> &diagram,
2851 vector<pair<pair<int, int>, double>> *mergePairs,
2852 vector<pair<pair<int, int>, double>> *splitPairs,
2853 vector<pair<pair<int, int>, double>> *pairs) const {
2854
2855 vector<pair<pair<int, int>, double>> defaultPairs{};
2856 if(!pairs) {
2857 pairs = &defaultPairs;
2858 }
2859
2860 if(!pairs->size())
2861 getPersistencePairs(*pairs, mergePairs, splitPairs);
2862
2863 // fast fix :(
2864 diagram.resize(pairs->size());
2865
2866 for(int i = 0; i < (int)pairs->size(); i++) {
2867 if((maximumList_) && (!minimumList_)) {
2868 // split tree
2869 diagram[i].second = (*vertexScalars_)[(*pairs)[i].first.first];
2870 diagram[i].first = (*vertexScalars_)[(*pairs)[i].first.second];
2871 } else {
2872 // join tree or contour tree
2873 diagram[i].first = (*vertexScalars_)[(*pairs)[i].first.first];
2874 diagram[i].second = (*vertexScalars_)[(*pairs)[i].first.second];
2875 }
2876 }
2877
2878 std::sort(diagram.begin(), diagram.end(), _pPairCmp);
2879
2880 return 0;
2881}
2882
2883int ContourTree::simplify(const double &simplificationThreshold,
2885#ifdef TTK_ENABLE_OPENMP
2886#pragma omp parallel sections
2887#endif
2888 {
2889
2890#ifdef TTK_ENABLE_OPENMP
2891#pragma omp section
2892#endif
2893 { mergeTree_.simplify(simplificationThreshold, metric); }
2894
2895#ifdef TTK_ENABLE_OPENMP
2896#pragma omp section
2897#endif
2898 { splitTree_.simplify(simplificationThreshold, metric); }
2899 }
2900
2901 return 0;
2902}
#define ttkNotUsed(x)
Mark function/method parameters that are not used in the function body at all.
Definition: BaseClass.h:47
struct filtrationCtCmp filtrationCmp
struct _persistenceCmp _pCmp
struct _persistencePairCmp _pPairCmp
#define REAL_SIGNIFICANT_DIGITS
Definition: Os.h:46
virtual SimplexId getVertexNeighborNumber(const SimplexId &vertexId) const
virtual SimplexId getNumberOfVertices() const
virtual int getVertexNeighbor(const SimplexId &vertexId, const int &localNeighborId, SimplexId &neighborId) const
int getUpNodeId() const
Definition: ContourTree.h:38
int getDownNodeId() const
Definition: ContourTree.h:34
virtual double computeSuperArcMetric(const int &downVertexId, const int &upVertexId, const std::vector< int > &interiorNodeIds)=0
int getPersistenceDiagram(std::vector< std::pair< double, double > > &diagram, std::vector< std::pair< std::pair< int, int >, double > > *mergePairs=nullptr, std::vector< std::pair< std::pair< int, int >, double > > *splitPairs=nullptr, std::vector< std::pair< std::pair< int, int >, double > > *pairs=nullptr) const
int computeSkeleton(unsigned int arcResolution=3) override
int simplify(const double &simplificationThreshold, ContourTreeSimplificationMetric *metric=nullptr) override
int smoothSkeleton(unsigned int skeletonSmoothing) override
int getPersistencePlot(std::vector< std::pair< double, int > > &plot, std::vector< std::pair< std::pair< int, int >, double > > *mergePairs=nullptr, std::vector< std::pair< std::pair< int, int >, double > > *splitPairs=nullptr, std::vector< std::pair< std::pair< int, int >, double > > *pairs=nullptr) const
bool isNodeEligible(const Node *n) const
SubLevelSetTree mergeTree_
Definition: ContourTree.h:653
int finalizeSuperArc(const int &nodeId, const int &arcId)
SubLevelSetTree splitTree_
Definition: ContourTree.h:653
int clearSkeleton() override
int getPersistencePairs(std::vector< std::pair< std::pair< int, int >, double > > &pairs, std::vector< std::pair< std::pair< int, int >, double > > *mergePairs=nullptr, std::vector< std::pair< std::pair< int, int >, double > > *splitPairs=nullptr) const override
int debugLevel_
Definition: Debug.h:379
int printWrn(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition: Debug.h:159
void setDebugMsgPrefix(const std::string &prefix)
Definition: Debug.h:364
virtual int setDebugLevel(const int &debugLevel)
Definition: Debug.cpp:147
int printMsg(const std::string &msg, const debug::Priority &priority=debug::Priority::INFO, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cout) const
Definition: Debug.h:118
int printErr(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition: Debug.h:149
int getUpArcId(const int &neighborId) const
Definition: ContourTree.h:166
int getNumberOfUpArcs() const
Definition: ContourTree.h:186
int getUpSuperArcId(const int &neighborId) const
Definition: ContourTree.h:172
int getDownArcId(const int &neighborId) const
Definition: ContourTree.h:154
double layoutY_
Definition: ContourTree.h:221
int getNumberOfDownArcs() const
Definition: ContourTree.h:178
int getNumberOfDownSuperArcs() const
Definition: ContourTree.h:182
int getVertexId() const
Definition: ContourTree.h:194
int getNumberOfUpSuperArcs() const
Definition: ContourTree.h:190
double computeSuperArcMetric(const int &downVertexId, const int &upVertexId, const std::vector< int > &interiorNodeIds) override
int clearRoot(const int &vertexId)
int exportPersistenceDiagram(const std::string &fileName="output.plot") const
bool isSosHigherThan(const int &vertexId0, const int &vertexId1) const
int getSkeletonScalars(const std::vector< double > &scalars, std::vector< std::vector< double > > &skeletonScalars) const
void setTriangulation(const AbstractTriangulation *const triangulation)
Definition: ContourTree.h:508
int moveRegularNode(const Node *n, const Node *oldDown, const Node *oldUp, const Node *newDown, const Node *newUp)
void setMaximumList(std::vector< int > &maximumList)
Definition: ContourTree.h:492
void setNumberOfVertices(const int &vertexNumber)
Definition: ContourTree.h:500
std::vector< Node > originalNodeList_
Definition: ContourTree.h:583
int exportToSvg(const std::string &fileName, const double &scaleX=1, const double &scaleY=1)
std::vector< int > * minimumList_
Definition: ContourTree.h:582
int clearRegularNode(const int &vertexId)
int appendRegularNode(const int &superArcId, const int &nodeId)
void setMinimumList(std::vector< int > &minimumList)
Definition: ContourTree.h:496
int getPersistenceDiagram(std::vector< std::pair< double, double > > &diagram, std::vector< std::pair< std::pair< int, int >, double > > *pairs=nullptr) const
const Node * getNode(const int &nodeId) const
Definition: ContourTree.h:305
int exportPersistenceCurve(const std::string &fileName="output.plot") const
virtual int clearSkeleton()
void setVertexScalars(const std::vector< real > *const vertexScalars)
Definition: ContourTree.h:517
const AbstractTriangulation * triangulation_
Definition: ContourTree.h:581
virtual int computeSkeleton(unsigned int arcResolution=3)
int buildSaddleList(std::vector< int > &vertexList) const
std::vector< int > vertex2superArcNode_
Definition: ContourTree.h:586
std::vector< int > * vertexSoSoffsets_
Definition: ContourTree.h:579
void setVertexPositions(std::vector< std::vector< double > > *vertexPositions)
Definition: ContourTree.h:513
std::vector< Arc > arcList_
Definition: ContourTree.h:584
std::vector< SuperArc > superArcList_
Definition: ContourTree.h:585
std::vector< SuperArc > originalSuperArcList_
Definition: ContourTree.h:585
int getNumberOfArcs() const
Definition: ContourTree.h:365
const Node * getVertexNode(const int &vertexId) const
Definition: ContourTree.h:446
int getNumberOfNodes() const
Definition: ContourTree.h:373
std::vector< int > * maximumList_
Definition: ContourTree.h:582
std::vector< Node > nodeList_
Definition: ContourTree.h:583
int getNumberOfSuperArcs() const
Definition: ContourTree.h:369
int exportArcPosToVtk(const int &arcId, const int &pointId, std::vector< int > &vertexIds, const std::vector< float > *origin, const std::vector< float > *voxelSize, std::ofstream &o)
int makeNode(const int &vertexId)
std::vector< int > vertex2node_
Definition: ContourTree.h:586
const SuperArc * getSuperArc(const int &superArcId) const
Definition: ContourTree.h:398
std::vector< std::vector< double > > * vertexPositions_
Definition: ContourTree.h:587
int sample(unsigned int samplingLevel=3)
int exportNodePosToVtk(const int &nodeId, const int &pointId, std::vector< int > &vertexIds, const std::vector< float > *origin, const std::vector< float > *voxelSize, std::ofstream &o)
int makeArc(const int &nodeId0, const int &nodeId1)
virtual int getPersistencePairs(std::vector< std::pair< std::pair< int, int >, double > > &pairs, std::vector< std::pair< std::pair< int, int >, double > > *mergePairs=nullptr, std::vector< std::pair< std::pair< int, int >, double > > *splitPairs=nullptr) const
int buildExtremumList(std::vector< int > &extremumList, const bool &isSubLevelSet=true)
virtual int simplify(const double &simplificationThreshold, ContourTreeSimplificationMetric *metric=nullptr)
int openSuperArc(const int &nodeId)
int exportToVtk(const std::string &fileName, const std::vector< float > *origin=nullptr, const std::vector< float > *voxelSize=nullptr)
int getPersistencePlot(std::vector< std::pair< double, int > > &plot, std::vector< std::pair< std::pair< int, int >, double > > *persistencePairs=nullptr) const
void setVertexSoSoffsets(std::vector< int > *vertexSoSoffsets)
Definition: ContourTree.h:530
const Node * getNodeUpNeighbor(const Node *n, const int &neighborId) const
Definition: ContourTree.h:336
int exportNodeColorToVtk(const int &nodeId, std::ofstream &o)
std::vector< int > vertex2superArc_
Definition: ContourTree.h:586
int getVertexScalar(const int &vertexId, double &scalar)
Definition: ContourTree.h:408
int closeSuperArc(const int &superArcId, const int &nodeId)
bool isSosLowerThan(const int &vertexId0, const int &vertexId1) const
bool buildPlanarLayout(const double &scaleX, const double &scaleY)
const std::vector< real > * vertexScalars_
Definition: ContourTree.h:578
virtual int smoothSkeleton(unsigned int skeletonSmoothing)
int clearArc(const int &vertexId0, const int &vertexId1)
int getNumberOfRegularNodes() const
Definition: ContourTree.h:70
int getNumberOfBarycenters() const
Definition: ContourTree.h:74
void getSample(const int &id, std::vector< int > &sample) const
Definition: ContourTree.h:100
int getRegularNodeId(const int &arcNodeId) const
Definition: ContourTree.h:82
bool isPruned() const
Definition: ContourTree.h:116
std::vector< std::vector< double > > barycenterList_
Definition: ContourTree.h:132
void smooth(const std::vector< Node > &nodeList, const std::vector< std::vector< double > > *vertexPositions, bool order=true)
Definition: ContourTree.cpp:84
std::vector< int > regularNodeList_
Definition: ContourTree.h:131
void sortRegularNodes(const std::vector< double > *vertexScalars, const std::vector< int > *vertexOffsets, const std::vector< Node > *nodeList, bool order=true)
int getNumberOfSamples() const
Definition: ContourTree.h:78
double getElapsedTime()
Definition: Timer.h:15
Union Find implementation for connectivity tracking.
Definition: UnionFind.h:17
static UnionFind * makeUnion(UnionFind *uf0, UnionFind *uf1)
Definition: UnionFind.h:40
UnionFind * find()
Definition: UnionFind.h:99
T powIntTen(const int n)
Compute the nth power of ten.
Definition: Geometry.h:360
The Topology ToolKit.
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22
MyCmp(const vector< double > *vertexScalars, const vector< int > *vertexOffsets, const vector< Node > *nodeList, bool isAscendingOrder)
Definition: ContourTree.cpp:17
bool isAscendingOrder_
Definition: ContourTree.cpp:15
const vector< double > * vertexScalars_
Definition: ContourTree.cpp:12
const vector< Node > * nodeList_
Definition: ContourTree.cpp:14
const vector< int > * vertexOffsets_
Definition: ContourTree.cpp:13
bool operator()(int node1, int node2)
Definition: ContourTree.cpp:25
const vector< double > * vertexScalars_
Definition: ContourTree.cpp:74
_persistenceCmp3(const vector< double > *vertexScalars)
Definition: ContourTree.cpp:64
bool operator()(const pair< pair< int, int >, double > &p0, const pair< pair< int, int >, double > &p1)
Definition: ContourTree.cpp:68
bool operator()(const pair< pair< int, int >, double > &p0, const pair< pair< int, int >, double > &p1)
Definition: ContourTree.cpp:55
bool operator()(const pair< double, double > &p0, const pair< double, double > &p1)
Definition: ContourTree.cpp:78
bool operator()(const pair< bool, pair< double, pair< int, int > > > &v0, const pair< bool, pair< double, pair< int, int > > > &v1) const
Definition: ContourTree.cpp:39