75 std::vector<int> &mdtTreePointComponentId,
76 std::vector<PointType> &mdtTreePointType,
77 std::vector<double> &mdtTreePointLowInterval,
78 std::vector<double> &mdtTreePointUpInterval,
79 std::vector<int> &mdtTreeEdgeSwitchable,
80 const std::vector<int> &mdtExtremumParentSaddle,
81 const std::vector<int> &mdtSaddleParentSaddle,
82 const std::vector<bool> &isExtremumSimplified,
83 const std::vector<bool> &isSaddleSimplified,
84 const std::vector<std::pair<double, double>> &extremumInterval,
85 const std::vector<std::pair<int, int>> &mandatorySaddleVertices,
86 const int extremaNumber,
87 const int saddleNumber,
91 const double globalOtherExtremumValue)
const {
93 const std::string treeName
98 mdtTreePointComponentId.clear();
99 mdtTreePointComponentId.reserve(extremaNumber + saddleNumber + 1);
100 mdtTreePointType.clear();
101 mdtTreePointType.reserve(extremaNumber + saddleNumber + 1);
102 mdtTreePointLowInterval.clear();
103 mdtTreePointLowInterval.reserve(extremaNumber + saddleNumber + 1);
104 mdtTreePointUpInterval.clear();
105 mdtTreePointUpInterval.reserve(extremaNumber + saddleNumber + 1);
106 mdtTreeEdgeSwitchable.clear();
109 std::vector<int> saddleGraphVertex(saddleNumber, -1);
111 for(
int i = 0; i < extremaNumber; i++) {
113 if(isExtremumSimplified[i])
116 int extremumGraphPoint = mdtTree.
addVertex();
117 mdtTreePointComponentId.push_back(i);
118 mdtTreePointType.push_back(extremumType);
120 mdtTreePointLowInterval.push_back(extremumInterval[i].first);
121 mdtTreePointUpInterval.push_back(extremumInterval[i].second);
123 int parentSaddle = mdtExtremumParentSaddle[i];
125 if(parentSaddle == -1)
128 if(saddleGraphVertex[parentSaddle] == -1) {
129 saddleGraphVertex[parentSaddle] = mdtTree.
addVertex();
130 mdtTreePointComponentId.push_back(parentSaddle);
131 mdtTreePointType.push_back(saddleType);
132 int lowerVertex = mandatorySaddleVertices[parentSaddle].first;
133 int upperVertex = mandatorySaddleVertices[parentSaddle].second;
138 mdtTree.
addEdge(saddleGraphVertex[parentSaddle], extremumGraphPoint);
140 mdtTreeEdgeSwitchable.push_back(0);
145 this->
printWrn(
"Mandatory " + treeName +
" empty.");
153 mdtTreePointComponentId.push_back(-1);
154 mdtTreePointType.push_back(otherExtremumType);
155 mdtTreePointLowInterval.push_back(globalOtherExtremumValue);
156 mdtTreePointUpInterval.push_back(globalOtherExtremumValue);
158 mdtTreeEdgeSwitchable.push_back(0);
159 this->
printWrn(
"Mandatory " + treeName +
" with only one minimum.");
164 for(
int i = 0; i < saddleNumber; i++) {
167 if(isSaddleSimplified[i]) {
171 if(saddleGraphVertex[i] == -1) {
172 saddleGraphVertex[i] = mdtTree.
addVertex();
173 mdtTreePointComponentId.push_back(i);
174 mdtTreePointType.push_back(saddleType);
175 int lowerVertex = mandatorySaddleVertices[i].first;
176 int upperVertex = mandatorySaddleVertices[i].second;
181 int parentSaddle = mdtSaddleParentSaddle[i];
184 if(parentSaddle != i) {
186 if(saddleGraphVertex[parentSaddle] == -1) {
187 saddleGraphVertex[parentSaddle] = mdtTree.
addVertex();
188 mdtTreePointComponentId.push_back(parentSaddle);
189 mdtTreePointType.push_back(saddleType);
190 int lowerVertex = mandatorySaddleVertices[parentSaddle].first;
191 int upperVertex = mandatorySaddleVertices[parentSaddle].second;
196 mdtTree.
addEdge(saddleGraphVertex[parentSaddle], saddleGraphVertex[i]);
198 mdtTreeEdgeSwitchable.push_back(
201 int globalOtherExtremumGraphPoint = mdtTree.
addVertex();
202 mdtTreePointComponentId.push_back(-1);
203 mdtTreePointType.push_back(otherExtremumType);
204 mdtTreePointLowInterval.push_back(globalOtherExtremumValue);
205 mdtTreePointUpInterval.push_back(globalOtherExtremumValue);
206 mdtTree.
addEdge(globalOtherExtremumGraphPoint, saddleGraphVertex[i]);
207 mdtTreeEdgeSwitchable.push_back(0);
213 this->
printMsg(
"Building of the mandatory " + treeName +
"");
217 this->
printMsg(
"List of mandatory " + treeName +
" graph points:",
220 std::stringstream msg{};
221 msg <<
" (" << std::setw(3) << std::right << i <<
") ";
222 msg << std::setw(12) << std::right;
223 switch(mdtTreePointType[i]) {
231 msg <<
"Join Saddle";
234 msg <<
"Split Saddle";
239 msg <<
" id = " << std::setw(3) << mdtTreePointComponentId[i];
240 msg <<
" I = [ " << mdtTreePointLowInterval[i];
241 msg <<
" ; " << mdtTreePointUpInterval[i];
249 std::stringstream msg{};
250 msg <<
" (" << std::setw(3) << std::right << i <<
") ";
252 msg << std::setw(3) << std::right
255 msg << std::setw(3) << std::left
265 const std::vector<std::pair<int, int>> &saddleList,
266 const std::vector<std::vector<int>> &mergedExtrema,
267 const std::vector<std::pair<double, double>> &extremumInterval,
270 std::vector<std::pair<std::pair<int, int>,
double>> &extremaSaddlePair)
273#ifndef TTK_ENABLE_KAMIKAZE
282 extremaSaddlePair.clear();
284 for(
size_t i = 0; i < mergedExtrema.size(); i++) {
285 for(
size_t j = 0; j < mergedExtrema[i].size(); j++) {
287 std::pair<std::pair<int, int>,
double> constructedPair{
288 {mergedExtrema[i][j], i}, 0.0};
290 double lowerValue = 0;
291 double upperValue = 0;
293 lowerValue = extremumInterval[constructedPair.first.first].first;
301 upperValue = extremumInterval[constructedPair.first.first].second;
304 constructedPair.second = fabs(upperValue - lowerValue);
306 extremaSaddlePair.push_back(constructedPair);
311 std::sort(extremaSaddlePair.begin(), extremaSaddlePair.end(),
pairComparison);
314 ?
"minimum - join saddle"
315 :
"maximum - split saddle";
318 for(
size_t i = 0; i < extremaSaddlePair.size(); i++) {
319 std::stringstream msg;
320 msg <<
" (" << std::setw(3) << extremaSaddlePair[i].first.first <<
";";
321 msg << std::setw(3) << extremaSaddlePair[i].first.second <<
")";
322 msg <<
" -> d = " << extremaSaddlePair[i].second;
331 const Graph &mdtTree,
332 const std::vector<PointType> &mdtTreePointType,
333 const std::vector<double> &mdtTreePointLowInterval,
334 const std::vector<double> &mdtTreePointUpInterval,
335 std::vector<double> &xCoord,
336 std::vector<double> &yCoord)
const {
342 const double range = rangeMax - rangeMin;
345 int rootGraphPointId = -1;
347 for(
size_t i = 0; i < mdtTreePointType.size(); i++) {
349 rootGraphPointId = i;
354 for(
size_t i = 0; i < mdtTreePointType.size(); i++) {
356 rootGraphPointId = i;
361 if(rootGraphPointId == -1) {
366 int downRootPointId = -1;
373 if(downRootPointId == -1) {
378 xCoord.resize(numberOfPoints);
379 yCoord.resize(numberOfPoints);
382 std::vector<std::pair<double, double>> xInterval(numberOfPoints);
383 xInterval[rootGraphPointId].first = -0.5 * (double)numberOfPoints;
384 xInterval[rootGraphPointId].second = 0.5 * (double)numberOfPoints;
386 std::vector<std::pair<int, double>> xOrder(numberOfPoints);
388 std::queue<int> graphPointQueue;
389 graphPointQueue.push(rootGraphPointId);
391 while(!graphPointQueue.empty()) {
393 int graphPoint = graphPointQueue.front();
394 graphPointQueue.pop();
398 int numberOfDownEdges = 0;
399 for(
int i = 0; i < numberOfEdges; i++) {
406 xOrder[graphPoint].first = graphPoint;
407 xOrder[graphPoint].second
408 = (xInterval[graphPoint].second + xInterval[graphPoint].first) * 0.5;
410 if(numberOfDownEdges == 0)
413 int downPointCount = 0;
415 for(
int i = 0; i < numberOfEdges; i++) {
417 int downGraphPoint = -1;
422 if(downGraphPoint == -1)
425 xInterval[downGraphPoint].first
426 = xInterval[graphPoint].first
427 + ((xInterval[graphPoint].second - xInterval[graphPoint].first)
428 * (
double)downPointCount / (double)numberOfDownEdges);
429 xInterval[downGraphPoint].second
430 = xInterval[graphPoint].first
431 + ((xInterval[graphPoint].second - xInterval[graphPoint].first)
432 * (
double)(downPointCount + 1) / (
double)numberOfDownEdges);
435 graphPointQueue.push(downGraphPoint);
440 for(
int i = 0; i < numberOfPoints; i++) {
442 = ((0.5 * (mdtTreePointLowInterval[i] + mdtTreePointUpInterval[i]))
450 std::sort(xOrder.begin(), xOrder.end(), xCoordCmp);
453 for(
int i = 0; i < numberOfPoints; i++) {
454 if(xOrder[i].first != rootGraphPointId) {
455 xCoord[xOrder[i].first]
456 = (double)pointCount * (1.0 / ((
int)numberOfPoints - 2.0));
461 xCoord[rootGraphPointId] = xCoord[downRootPointId];
463 const std::string treeName
466 this->
printMsg(
"Planar layout for mandatory " + treeName +
": root point ("
467 + std::to_string(downRootPointId) +
")",
470 for(
int i = 0; i < numberOfPoints; i++) {
471 std::stringstream msg;
472 msg <<
" (" << i <<
") :"
473 <<
" x = " << xCoord[i] <<
" y = " << yCoord[i];
483 const int seedVertexId,
484 const std::vector<double> &vertexScalars,
485 std::vector<int> &componentVertexList)
const {
487 const double value = vertexScalars[seedVertexId];
490 componentVertexList.clear();
496 std::vector<int> subTreeSuperArcId;
501 for(
int i = 0; i < (int)subTreeSuperArcId.size(); i++) {
505 if(subTreeSuperArcId[i] == rootSuperArcId) {
507 for(
int j = 0; j < numberOfRegularNodes; j++) {
510 if(!((sign * nodeScalar) > (sign * value))) {
512 componentVertexList.push_back(vertexId);
518 if(!((sign * nodeScalar) > (sign * value))) {
520 componentVertexList.push_back(vertexId);
524 for(
int j = 0; j < numberOfRegularNodes; j++) {
527 componentVertexList.push_back(vertexId);
532 componentVertexList.push_back(vertexId);
542 std::vector<int> &mandatoryExtremum,
543 std::vector<std::pair<double, double>> &criticalInterval)
const {
547#ifndef TTK_ENABLE_KAMIKAZE
560 size_t extremumNumber = extremumList.size();
561 std::vector<int> vertexId(extremumNumber);
562 std::vector<double> vertexValue(extremumNumber);
563 std::vector<int> superArcId(extremumNumber);
565 std::vector<int> subTreeSuperArcId;
568 mandatoryExtremum.clear();
569 criticalInterval.clear();
572 std::vector<bool> isSuperArcAlreadyVisited(
578 for(
size_t i = 0; i < extremumNumber; i++) {
580 bool isMandatory =
true;
582 vertexId[i] = extremumList[i];
589 int rootSuperArcId = -1;
590 if(!isSuperArcAlreadyVisited[secondTreeSuperArcId]) {
592 &secondTree, secondTreeSuperArcId, vertexValue[i]);
599 subTreeSuperArcId.clear();
600 std::queue<int> superArcQueue;
601 superArcQueue.push(rootSuperArcId);
602 while(isMandatory && (!superArcQueue.empty())) {
603 int spaId = superArcQueue.front();
605 if(isSuperArcAlreadyVisited[spaId]) {
609 int numberOfDownSuperArcs
611 if(numberOfDownSuperArcs > 0) {
612 for(
int j = 0; j < numberOfDownSuperArcs; j++) {
617 subTreeSuperArcId.push_back(spaId);
625 mandatoryExtremum.push_back(vertexId[i]);
628 criticalInterval.emplace_back(vertexValue[i], vertexValue[i]);
629 for(
int j = 0; j < (int)subTreeSuperArcId.size(); j++) {
630 isSuperArcAlreadyVisited[subTreeSuperArcId[j]] =
true;
635 if(downNodeValue < criticalInterval.back().first) {
636 criticalInterval.back().first = downNodeValue;
639 if(downNodeValue > criticalInterval.back().second) {
640 criticalInterval.back().second = downNodeValue;
647 while((spaId != -1) && !(isSuperArcAlreadyVisited[spaId])) {
648 isSuperArcAlreadyVisited[spaId] =
true;
690#ifdef TTK_ENABLE_OPENMP
697 this->
printMsg(
"Computed " + std::to_string(mandatoryExtremum.size())
698 +
" mandatory " + pt,
702 for(
size_t i = 0; i < mandatoryExtremum.size(); i++) {
703 std::stringstream msg;
704 msg <<
" -> " << pt <<
" (" << std::setw(3) << std::right << i <<
") ";
706 <<
"Vertex " << std::setw(9) << std::left << mandatoryExtremum[i];
707 msg <<
" \tInterval [ " << std::setw(12) << std::right
708 << criticalInterval[i].first <<
" ; " << std::setprecision(9)
709 << std::setw(11) << std::left << criticalInterval[i].second <<
" ]";
721 const std::vector<int> &mandatoryExtremumVertex,
722 std::vector<std::pair<int, int>> &mandatorySaddleVertex,
723 std::vector<std::vector<int>> &mandatoryMergedExtrema) {
742 const unsigned int extremaNumber = mandatoryExtremumVertex.size();
744 std::vector<int> upperSaddleList;
745 std::vector<int> lowerSaddleList;
750 std::vector<std::vector<int>> lowerToUpperLinks;
751 std::vector<std::vector<int>> upperToLowerLinks;
753 std::vector<std::vector<int>> mergedExtrema;
761#ifdef TTK_ENABLE_OPENMP
764 for(
int i = 0; i < (int)mandatoryExtremumVertex.size(); i++) {
772 bool multipleSaddleFound;
775 multipleSaddleFound =
false;
778#ifdef TTK_ENABLE_OPENMP
779#pragma omp critical(upperTreeTransverse)
783 if(++upperTransverse[upperSuperArcId]) {
787 if(upperTransverse[upperSuperArcId] > 1) {
788 multipleSaddleFound =
true;
796 if(upperSuperArcId == -1) {
800 if(!multipleSaddleFound) {
803#ifdef TTK_ENABLE_OPENMP
806 upperSaddleList.push_back(saddleId);
809 }
while(!(saddleFound || rootReached));
812 multipleSaddleFound =
false;
815#ifdef TTK_ENABLE_OPENMP
816#pragma omp critical(lowerTreeTransverse)
820 if(++lowerTransverse[lowerSuperArcId]) {
824 if(lowerTransverse[lowerSuperArcId] > 1) {
825 multipleSaddleFound =
true;
833 if(lowerSuperArcId == -1) {
837 if(!multipleSaddleFound) {
840#ifdef TTK_ENABLE_OPENMP
843 lowerSaddleList.push_back(saddleId);
846 }
while(!(saddleFound || rootReached));
853#ifdef TTK_ENABLE_OPENMP
854#pragma omp parallel num_threads(threadNumber_)
859#ifdef TTK_ENABLE_OPENMP
862 for(
int i = 0; i < (int)upperTransverse.size(); i++) {
863 upperTransverse[i] = -1;
866#ifdef TTK_ENABLE_OPENMP
869 for(
int i = 0; i < (int)lowerTransverse.size(); i++) {
870 lowerTransverse[i] = -1;
874#ifdef TTK_ENABLE_OPENMP
877 for(
int i = 0; i < (int)upperSaddleList.size(); i++) {
880 upperTransverse[superArcId] = i;
883#ifdef TTK_ENABLE_OPENMP
886 for(
int i = 0; i < (int)lowerSaddleList.size(); i++) {
889 lowerTransverse[superArcId] = i;
893#ifdef TTK_ENABLE_OPENMP
898#ifdef TTK_ENABLE_OPENMP
902 upperLca.
addNodes(mandatoryExtremumVertex.size()
903 + upperSaddleList.size());
906#ifdef TTK_ENABLE_OPENMP
910 lowerLca.
addNodes(mandatoryExtremumVertex.size()
911 + lowerSaddleList.size());
916#ifdef TTK_ENABLE_OPENMP
919 for(
int i = 0; i < (int)mandatoryExtremumVertex.size(); i++) {
924 while((superArcId != -1) && (upperTransverse[superArcId] == -1)) {
928 if(superArcId != -1) {
929 lcaSaddleId = extremaNumber + upperTransverse[superArcId];
933#ifdef TTK_ENABLE_OPENMP
934#pragma omp critical(upperLca)
941 while((superArcId != -1) && (lowerTransverse[superArcId] == -1)) {
945 if(superArcId != -1) {
946 lcaSaddleId = extremaNumber + lowerTransverse[superArcId];
950#ifdef TTK_ENABLE_OPENMP
951#pragma omp critical(lowerLca)
958#ifdef TTK_ENABLE_OPENMP
961 for(
int i = 0; i < (int)upperSaddleList.size(); i++) {
967 }
while((superArcId != -1) && (upperTransverse[superArcId] == -1));
968 if(superArcId != -1) {
969 int lcaSuccessorSaddleId = extremaNumber + i;
970 int lcaAncestorSaddleId = extremaNumber + upperTransverse[superArcId];
974#ifdef TTK_ENABLE_OPENMP
975#pragma omp critical(upperLca)
977 upperLca.
getNode(lcaAncestorSaddleId)
982#ifdef TTK_ENABLE_OPENMP
985 for(
int i = 0; i < (int)lowerSaddleList.size(); i++) {
991 }
while((superArcId != -1) && (lowerTransverse[superArcId] == -1));
992 if(superArcId != -1) {
993 int lcaSuccessorSaddleId = extremaNumber + i;
994 int lcaAncestorSaddleId = extremaNumber + lowerTransverse[superArcId];
998#ifdef TTK_ENABLE_OPENMP
999#pragma omp critical(lowerLca)
1001 lowerLca.
getNode(lcaAncestorSaddleId)
1007#ifdef TTK_ENABLE_OPENMP
1011#ifdef TTK_ENABLE_OPENMP
1015#ifdef TTK_ENABLE_OPENMP
1022 std::vector<std::vector<int>> localLowerToUpperLinks;
1023 std::vector<std::vector<int>> localUpperToLowerLinks;
1025 std::vector<std::vector<int>> localMergedExtrema;
1027 localLowerToUpperLinks.resize(lowerSaddleList.size());
1028 localUpperToLowerLinks.resize(upperSaddleList.size());
1029 localMergedExtrema.resize(lowerSaddleList.size());
1031 unsigned int kmax = (extremaNumber * (extremaNumber - 1)) / 2;
1033#ifdef TTK_ENABLE_OPENMP
1036 for(
int k = 0; k < (int)kmax; k++) {
1037 unsigned int i = k / extremaNumber;
1038 unsigned int j = k % extremaNumber;
1040 i = extremaNumber - i - 2;
1041 j = extremaNumber - j - 1;
1043 int uppLCA = upperLca.
query(i, j) - extremaNumber;
1044 int lowLCA = lowerLca.
query(i, j) - extremaNumber;
1045 localLowerToUpperLinks[lowLCA].push_back(uppLCA);
1046 localUpperToLowerLinks[uppLCA].push_back(lowLCA);
1047 localMergedExtrema[lowLCA].push_back(i);
1048 localMergedExtrema[lowLCA].push_back(j);
1052#ifdef TTK_ENABLE_OPENMP
1056 lowerToUpperLinks.resize(lowerSaddleList.size());
1057 upperToLowerLinks.resize(upperSaddleList.size());
1058 mergedExtrema.resize(lowerSaddleList.size());
1061 for(
unsigned int i = 0; i < localLowerToUpperLinks.size(); i++) {
1062 std::vector<int>::iterator newEnd;
1064 localLowerToUpperLinks[i].begin(), localLowerToUpperLinks[i].end());
1065#ifdef TTK_ENABLE_OPENMP
1066#pragma omp critical(lowerToUpperFusion)
1069 lowerToUpperLinks[i].insert(
1070 lowerToUpperLinks[i].end(),
1071 make_move_iterator(localLowerToUpperLinks[i].begin()),
1072 make_move_iterator(newEnd));
1074 localLowerToUpperLinks[i].clear();
1076 localLowerToUpperLinks.clear();
1078 for(
unsigned int i = 0; i < localUpperToLowerLinks.size(); i++) {
1079 std::vector<int>::iterator newEnd;
1081 localUpperToLowerLinks[i].begin(), localUpperToLowerLinks[i].end());
1082#ifdef TTK_ENABLE_OPENMP
1083#pragma omp critical(upperToLowerFusion)
1086 upperToLowerLinks[i].insert(
1087 upperToLowerLinks[i].end(),
1088 make_move_iterator(localUpperToLowerLinks[i].begin()),
1089 make_move_iterator(newEnd));
1091 localUpperToLowerLinks[i].clear();
1093 localUpperToLowerLinks.clear();
1095 for(
unsigned int i = 0; i < localMergedExtrema.size(); i++) {
1096 std::vector<int>::iterator newEnd;
1097 std::sort(localMergedExtrema[i].begin(), localMergedExtrema[i].end());
1099 = unique(localMergedExtrema[i].begin(), localMergedExtrema[i].end());
1100#ifdef TTK_ENABLE_OPENMP
1101#pragma omp critical(mergedExtremaFusion)
1104 mergedExtrema[i].insert(
1105 mergedExtrema[i].end(),
1106 make_move_iterator(localMergedExtrema[i].begin()),
1107 make_move_iterator(newEnd));
1109 localMergedExtrema[i].clear();
1111 localMergedExtrema.clear();
1119 std::vector<bool> isLowerVisited(lowerSaddleList.size(),
false);
1120 std::vector<bool> isUpperVisited(upperSaddleList.size(),
false);
1121 std::vector<std::vector<int>> lowComponent;
1122 std::vector<std::vector<int>> uppComponent;
1123 for(
unsigned int i = 0; i < lowerSaddleList.size(); i++) {
1124 if(!isLowerVisited[i]) {
1126 lowComponent.emplace_back();
1127 uppComponent.emplace_back();
1128 int componentId =
static_cast<int>(lowComponent.size()) - 1;
1129 lowComponent[componentId].push_back(i);
1130 isLowerVisited[i] =
true;
1132 std::vector<int> nextList;
1133 std::vector<int> currentList;
1134 currentList.push_back(i);
1136 std::vector<bool> *isVisited = &(isUpperVisited);
1137 std::vector<std::vector<int>> *component = &(uppComponent);
1138 std::vector<std::vector<int>> *linkList = &(lowerToUpperLinks);
1139 while(!currentList.empty()) {
1141 for(
unsigned int j = 0; j < currentList.size(); j++) {
1143 for(
unsigned int k = 0; k < (*linkList)[currentList[j]].size(); k++) {
1145 int neighbor = (*linkList)[currentList[j]][k];
1147 if(!(*isVisited)[neighbor]) {
1149 (*isVisited)[neighbor] =
true;
1150 (*component)[componentId].push_back(neighbor);
1151 nextList.push_back(neighbor);
1156 isVisited = (isVisited == &(isUpperVisited)) ? &(isLowerVisited)
1157 : &(isUpperVisited);
1159 = (component == &(uppComponent)) ? &(lowComponent) : &(uppComponent);
1160 linkList = (linkList == &(lowerToUpperLinks)) ? &(upperToLowerLinks)
1161 : &(lowerToUpperLinks);
1162 currentList.swap(nextList);
1169 int numberOfComponents =
static_cast<int>(lowComponent.size());
1170 mandatorySaddleVertex.resize(numberOfComponents, std::pair<int, int>(-1, -1));
1171 std::vector<std::vector<int>> mandatoryMergedExtrema_tmp;
1172 mandatoryMergedExtrema_tmp.resize(numberOfComponents, std::vector<int>());
1174 for(
int i = 0; i < numberOfComponents; i++) {
1177 nodeId = lowerSaddleList[lowComponent[i][0]];
1179 for(
unsigned int j = 0; j < lowComponent[i].size(); j++) {
1181 int nId = lowerSaddleList[lowComponent[i][j]];
1183 double refScalar = 0;
1185 if(nodeScalar < refScalar) {
1190 for(
unsigned int j = 0; j < lowComponent[i].size(); j++) {
1191 mandatoryMergedExtrema_tmp[i].insert(
1192 mandatoryMergedExtrema_tmp[i].end(),
1193 make_move_iterator(mergedExtrema[lowComponent[i][j]].begin()),
1194 make_move_iterator(mergedExtrema[lowComponent[i][j]].end()));
1195 mergedExtrema[lowComponent[i][j]].clear();
1199 nodeId = upperSaddleList[uppComponent[i][0]];
1201 for(
unsigned int j = 0; j < uppComponent[i].size(); j++) {
1203 int nId = upperSaddleList[uppComponent[i][j]];
1205 double refScalar = 0;
1206 upperTree.
getVertexScalar(mandatorySaddleVertex[i].second, refScalar);
1207 if(nodeScalar > refScalar) {
1214 std::vector<std::pair<int, int>> order;
1215 for(
unsigned int i = 0; i < mandatoryMergedExtrema_tmp.size(); i++) {
1216 std::sort(mandatoryMergedExtrema_tmp[i].begin(),
1217 mandatoryMergedExtrema_tmp[i].end());
1218 std::vector<int>::iterator newEnd
1219 = unique(mandatoryMergedExtrema_tmp[i].begin(),
1220 mandatoryMergedExtrema_tmp[i].end());
1221 mandatoryMergedExtrema_tmp[i].resize(
1222 distance(mandatoryMergedExtrema_tmp[i].begin(), newEnd));
1223 mandatoryMergedExtrema_tmp[i].shrink_to_fit();
1224 order.emplace_back(i, mandatoryMergedExtrema_tmp[i].size());
1228 std::sort(order.begin(), order.end(), cmp);
1230 mandatoryMergedExtrema.clear();
1231 mandatoryMergedExtrema.resize(mandatoryMergedExtrema_tmp.size());
1232 for(
unsigned int i = 0; i < order.size(); i++) {
1233 mandatoryMergedExtrema_tmp[order[i].first].swap(mandatoryMergedExtrema[i]);
1237#ifdef TTK_ENABLE_OPENMP
1238#pragma omp parallel sections
1241#ifdef TTK_ENABLE_OPENMP
1254#ifdef TTK_ENABLE_OPENMP
1269 const std::string st
1272 this->
printMsg(
"Computed " + std::to_string(mandatorySaddleVertex.size())
1273 +
" mandatory " + st,
1277 for(
size_t i = 0; i < mandatorySaddleVertex.size(); i++) {
1278 std::stringstream msg;
1279 msg <<
" -> " << st <<
" (";
1280 msg << std::setw(3) << std::right << i <<
")";
1281 msg <<
" Vertices <";
1282 msg << std::setw(9) << std::right << mandatorySaddleVertex[i].first
1284 msg << std::setw(9) << std::left << mandatorySaddleVertex[i].second <<
">";
1285 msg <<
" Interval [";
1286 double lowerValue = 0;
1287 double upperValue = 0;
1288 lowerTree.
getVertexScalar(mandatorySaddleVertex[i].first, lowerValue);
1289 upperTree.
getVertexScalar(mandatorySaddleVertex[i].second, upperValue);
1290 msg << std::setw(12) << std::right << lowerValue <<
" ; ";
1291 msg << std::setw(12) << std::left << upperValue <<
"]";
1300 const int &startingSuperArcId,
1301 const double &targetValue)
const {
1303 int superArcId = startingSuperArcId;
1305 if(superArcId != -1) {
1313 while(!((sign * targetValue) < (sign * upNodeValue))) {
1314 int numberOfUpSuperArcs
1316 if(numberOfUpSuperArcs > 0) {
1331 const int &vertexId0,
1332 const int &vertexId1)
const {
1334 std::vector<bool> isSuperArcVisited(numberOfSuperArcs,
false);
1337 int superArcId = superArcId0;
1339 isSuperArcVisited[superArcId] =
true;
1342 }
while(superArcId != -1);
1343 superArcId = superArcId1;
1344 while(!isSuperArcVisited[superArcId]) {
1353 const int &rootSuperArcId,
1354 std::vector<int> &subTreeSuperArcId)
const {
1355 std::vector<int> listOfSuperArcId;
1356 std::queue<int> superArcIdsToCompute;
1357 superArcIdsToCompute.push(rootSuperArcId);
1358 while(!superArcIdsToCompute.empty()) {
1359 int superArcId = superArcIdsToCompute.front();
1361 int numberOfUpSuperArcs
1363 if(numberOfUpSuperArcs > 0) {
1364 for(
int i = 0; i < numberOfUpSuperArcs; i++)
1365 superArcIdsToCompute.push(
1368 listOfSuperArcId.push_back(superArcId);
1369 superArcIdsToCompute.pop();
1371 listOfSuperArcId.swap(subTreeSuperArcId);
1375 const double normalizedThreshold,
1377 const std::vector<std::pair<std::pair<int, int>,
double>> &extremaSaddlePair,
1378 const std::vector<std::vector<int>> &mergedExtrema,
1379 const int numberOfExtrema,
1380 std::vector<bool> &extremumSimplified,
1381 std::vector<bool> &saddleSimplified,
1382 std::vector<int> &extremumParentSaddle,
1383 std::vector<int> &saddleParentSaddle)
const {
1386 double simplificationThreshold;
1387 if(normalizedThreshold > 1.0)
1389 else if(normalizedThreshold < 0.0)
1390 simplificationThreshold = 0.0;
1392 simplificationThreshold
1395 int extremaSimplifiedNumber = 0;
1396 int saddlesSimplifiedNumber = 0;
1399 extremumSimplified.resize(numberOfExtrema);
1400 fill(extremumSimplified.begin(), extremumSimplified.end(),
false);
1401 for(
size_t i = 0; i < extremaSaddlePair.size(); i++) {
1402 if(extremaSaddlePair[i].second < simplificationThreshold) {
1403 extremumSimplified[extremaSaddlePair[i].first.first] =
true;
1404 extremaSimplifiedNumber++;
1410 saddleSimplified.resize(mergedExtrema.size());
1411 fill(saddleSimplified.begin(), saddleSimplified.end(),
false);
1412 std::vector<int> nonSimplifiedExtremaNumber(mergedExtrema.size(), 0);
1413 extremumParentSaddle.resize(numberOfExtrema);
1414 fill(extremumParentSaddle.begin(), extremumParentSaddle.end(), -1);
1415 saddleParentSaddle.resize(mergedExtrema.size());
1416 fill(saddleParentSaddle.begin(), saddleParentSaddle.end(), -1);
1417 for(
size_t i = 0; i < mergedExtrema.size(); i++) {
1418 saddleParentSaddle[i] = i;
1419 for(
size_t j = 0; j < mergedExtrema[i].size(); j++) {
1420 if(!extremumSimplified[mergedExtrema[i][j]])
1421 nonSimplifiedExtremaNumber[i]++;
1423 if(nonSimplifiedExtremaNumber[i] > 1) {
1426 while(extremumSimplified[mergedExtrema[i][extremum]]) {
1429 if(!(extremumParentSaddle[mergedExtrema[i][extremum]] == -1)) {
1431 int rootSaddleId = extremumParentSaddle[mergedExtrema[i][extremum]];
1432 int lastSaddleId = -1;
1433 while(rootSaddleId != lastSaddleId) {
1434 lastSaddleId = rootSaddleId;
1435 rootSaddleId = saddleParentSaddle[rootSaddleId];
1438 if(!(nonSimplifiedExtremaNumber[i]
1439 > nonSimplifiedExtremaNumber[rootSaddleId])) {
1440 saddleSimplified[i] =
true;
1441 saddlesSimplifiedNumber++;
1445 saddleSimplified[i] =
true;
1446 saddlesSimplifiedNumber++;
1448 if(!saddleSimplified[i]) {
1449 for(
size_t j = 0; j < mergedExtrema[i].size(); j++) {
1450 int extremaId = mergedExtrema[i][j];
1451 if(!extremumSimplified[extremaId]) {
1453 if(extremumParentSaddle[extremaId] == -1) {
1454 extremumParentSaddle[extremaId] = i;
1457 int rootSaddleId = extremumParentSaddle[extremaId];
1458 int lastSaddleId = -1;
1459 while(rootSaddleId != lastSaddleId) {
1460 lastSaddleId = rootSaddleId;
1461 rootSaddleId = saddleParentSaddle[rootSaddleId];
1464 saddleParentSaddle[rootSaddleId] = i;
1474 "#Simplified " + pt +
": " + std::to_string(extremaSimplifiedNumber)
1475 +
" (threshold value = " + std::to_string(simplificationThreshold) +
")");
1477 if(extremaSimplifiedNumber > 0) {
1479 for(
size_t i = 0; i < extremumSimplified.size(); i++) {
1480 if(extremumSimplified[i]) {
1481 std::stringstream msg;
1482 msg <<
" (" << i <<
")";
1488 const std::string st
1492 "#Simplified " + st +
": " + std::to_string(saddlesSimplifiedNumber)
1493 +
" (threshold value = " + std::to_string(simplificationThreshold) +
")");
1495 if(saddlesSimplifiedNumber > 0) {
1497 for(
size_t i = 0; i < saddleSimplified.size(); i++) {
1498 if(saddleSimplified[i]) {
1499 std::stringstream msg;
1500 msg <<
" (" << i <<
")";
int getDownNodeId() const
int printWrn(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
void setDebugMsgPrefix(const std::string &prefix)
virtual int setDebugLevel(const int &debugLevel)
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
const std::pair< int, int > & getVertexIdx() const
int getNumberOfEdges() const
int getEdgeIdx(const int &connectedEdge) const
const Edge & getEdge(const int idx) const
int addEdge(const int &start, const int &end)
Add an edge between the vertex start and end, returns it's index.
const Vertex & getVertex(const int idx) const
Get a pointer to the vertex idx.
int addVertex()
Add a vertex, returns it's index.
int getNumberOfVertices() const
int getNumberOfEdges() const
void addSuccessor(const int &id)
void setAncestor(const int &id)
Class to answer the lowest common ancestor requests of pairs of nodes in a tree in constant time afte...
int query(int i, int j) const
Node & getNode(const unsigned int &id)
void addNodes(const unsigned int &number)
std::vector< std::pair< int, int > > mandatorySplitSaddleVertex_
Pair of mandatory vertices for each split saddle component.
std::vector< std::pair< std::pair< int, int >, double > > mdtMaxSplitSaddlePair_
Pairs ( (M,S), d(M,S) ) Of maxima and split saddles.
std::vector< double > mdtJoinTreePointXCoord_
int vertexNumber_
Number of vertices.
std::vector< std::vector< int > > mandatoryMinimumComponentVertices_
std::vector< double > mdtSplitTreePointLowInterval_
std::vector< double > mdtJoinTreePointLowInterval_
int computePlanarLayout(const TreeType &treeType, const Graph &mdtTree, const std::vector< PointType > &mdtTreePointType, const std::vector< double > &mdtTreePointLowInterval, const std::vector< double > &mdtTreePointUpInterval, std::vector< double > &xCoord, std::vector< double > &yCoord) const
double getGlobalMaximum() const
int enumerateMandatorySaddles(const PointType pointType, SubLevelSetTree &lowerTree, SubLevelSetTree &upperTree, const std::vector< int > &mandatoryExtremumVertex, std::vector< std::pair< int, int > > &mandatorySaddleVertex, std::vector< std::vector< int > > &mandatoryMergedExtrema)
TODO : Multiplicity.
std::vector< PointType > mdtJoinTreePointType_
bool areSaddlesSwitchables(const TreeType treeType, const int &firstId, const int &secondId) const
MandatoryCriticalPoints()
int enumerateMandatoryExtrema(const PointType pointType, SubLevelSetTree &firstTree, SubLevelSetTree &secondTree, std::vector< int > &mandatoryExtremum, std::vector< std::pair< double, double > > &criticalInterval) const
std::vector< std::vector< double > > vertexPositions_
Position (x,y,z) of each vertex.
std::vector< std::vector< int > > mandatoryMaximumComponentVertices_
List of the vertices forming each of the mandatory maximum components.
int buildMandatoryTree(const TreeType treeType, Graph &mdtTree, std::vector< int > &mdtTreePointComponentId, std::vector< PointType > &mdtTreePointType, std::vector< double > &mdtTreePointLowInterval, std::vector< double > &mdtTreePointUpInterval, std::vector< int > &mdtTreeEdgeSwitchable, const std::vector< int > &mdtExtremumParentSaddle, const std::vector< int > &mdtSaddleParentSaddle, const std::vector< bool > &isExtremumSimplified, const std::vector< bool > &isSaddleSimplified, const std::vector< std::pair< double, double > > &extremumInterval, const std::vector< std::pair< int, int > > &mandatorySaddleVertices, const int extremaNumber, const int saddleNumber, const PointType extremumType, const PointType saddleType, const PointType otherExtremumType, const double globalOtherExtremumValue) const
int * outputMandatoryMinimum_
Void pointer to the output mandatory minima components.
int getVertexSuperArcId(const int &vertexId, const SubLevelSetTree *tree) const
std::vector< int > mdtSplitSaddleParentSaddleId_
void getSubTreeSuperArcIds(const SubLevelSetTree *tree, const int &rootSuperArcId, std::vector< int > &subTreeSuperArcId) const
SubLevelSetTree lowerJoinTree_
Join tree of the lower bound scalar field.
std::vector< std::vector< int > > mergedMinimaId_
std::vector< int > upperMinimumList_
List of vertex id of the minima in the upper bound scalar field.
void * inputUpperBoundField_
Void pointer to the input upper bound scalar field.
int * outputMandatoryJoinSaddle_
Void pointer to the output mandatory join saddles components.
std::vector< int > mdtJoinSaddleParentSaddleId_
SubLevelSetTree lowerSplitTree_
Split tree of the lower bound scalar field.
std::vector< PointType > mdtSplitTreePointType_
void * inputLowerBoundField_
Void pointer to the input lower bound scalar field.
std::vector< int > mdtMaximumParentSaddleId_
std::vector< std::pair< std::pair< int, int >, double > > mdtMinJoinSaddlePair_
Pairs ( (M,S), d(M,S) ) Of minima and join saddles.
int simplify(const double normalizedThreshold, const TreeType treeType, const std::vector< std::pair< std::pair< int, int >, double > > &extremaSaddlePair, const std::vector< std::vector< int > > &mergedExtrema, const int numberOfExtrema, std::vector< bool > &extremumSimplified, std::vector< bool > &saddleSimplified, std::vector< int > &extremumParentSaddle, std::vector< int > &saddleParentSaddle) const
std::vector< int > lowerMaximumList_
List of vertex id of the maxima in the lower bound scalar field.
std::vector< double > upperVertexScalars_
Copy of the input upper scalar field converted in double.
std::vector< double > mdtJoinTreePointUpInterval_
std::vector< bool > isMdtMinimumSimplified_
std::vector< int > mdtJoinTreePointComponentId_
std::vector< double > lowerVertexScalars_
Copy of the input lower scalar field converted in double.
std::vector< bool > isMdtSplitSaddleSimplified_
double normalizedThreshold_
Value of the simplification threshold.
SubLevelSetTree upperJoinTree_
Join tree of the upper bound scalar field.
int findCommonAncestorNodeId(const SubLevelSetTree *tree, const int &vertexId0, const int &vertexId1) const
std::vector< int > mdtSplitTreePointComponentId_
std::vector< std::pair< int, int > > mandatoryJoinSaddleVertex_
Pair of mandatory vertices for each join saddle component.
std::vector< bool > isMdtMaximumSimplified_
SubLevelSetTree upperSplitTree_
Split tree of the upper bound scalar field.
std::vector< int > vertexSoSoffsets_
Offsets.
double globalMaximumValue_
std::vector< std::pair< double, double > > mandatoryMaximumInterval_
Critical interval for each maximum component.
double getGlobalMinimum() const
int getSubTreeRootSuperArcId(const SubLevelSetTree *tree, const int &startingSuperArcId, const double &targetValue) const
std::vector< double > mdtSplitTreePointXCoord_
int computeExtremumComponent(const PointType &pointType, const SubLevelSetTree &tree, const int seedVertexId, const std::vector< double > &vertexScalars, std::vector< int > &componentVertexList) const
int buildPairs(const TreeType treeType, const std::vector< std::pair< int, int > > &saddleList, const std::vector< std::vector< int > > &mergedExtrema, const std::vector< std::pair< double, double > > &extremumInterval, SubLevelSetTree &lowerTree, SubLevelSetTree &upperTree, std::vector< std::pair< std::pair< int, int >, double > > &extremaSaddlePair) const
TODO : Replace SubLevelSetTrees by scalar fields for vertex value.
std::vector< std::vector< int > > mandatorySplitSaddleComponentVertices_
std::vector< int > lowerMinimumList_
List of vertex id of the minima in the lower bound scalar field.
int * outputMandatorySplitSaddle_
Void pointer to the output mandatory split saddles components.
std::vector< std::vector< int > > mandatoryJoinSaddleComponentVertices_
std::vector< int > mdtMinimumParentSaddleId_
std::vector< std::vector< int > > mergedMaximaId_
std::vector< std::pair< double, double > > mandatoryMinimumInterval_
Critical interval for each minimum component.
std::vector< double > mdtSplitTreePointYCoord_
std::vector< double > mdtJoinTreePointYCoord_
double globalMinimumValue_
std::vector< double > mdtSplitTreePointUpInterval_
int * outputMandatoryMaximum_
Void pointer to the output mandatory maxima components.
std::vector< bool > isMdtJoinSaddleSimplified_
std::vector< int > upperMaximumList_
List of vertex id of the maxima in the upper bound scalar field.
std::vector< int > mandatoryMinimumVertex_
Mandatory vertex for each minimum component.
std::vector< int > mandatoryMaximumVertex_
Mandatory vertex for each maximum component.
int getUpSuperArcId(const int &neighborId) const
int getNumberOfDownSuperArcs() const
int getNumberOfUpSuperArcs() const
int getDownSuperArcId(const int &neighborId) const
const std::vector< int > * getExtremumList() const
const Node * getNode(const int &nodeId) const
double getNodeScalar(const int &nodeId) const
int getNumberOfSuperArcs() const
const SuperArc * getSuperArc(const int &superArcId) const
int getVertexScalar(const int &vertexId, double &scalar)
int getNumberOfRegularNodes() const
int getRegularNodeId(const int &arcNodeId) const
Comparison between critical point pairs ( (Extremum,Saddle), dist(M,S) )
Comparison between mandatory saddles (Saddle id, Number of merged extrema)
Comparison of the second member of a std::pair<int,double>