TTK
Loading...
Searching...
No Matches
PersistenceDiagramCompression.h
Go to the documentation of this file.
1//
2// Created by max on 24/05/18.
3//
4
5#pragma once
6
9
10template <typename triangulationType>
11int ttk::TopologicalCompression::ReadPersistenceGeometry(
12 FILE *fm, const triangulationType &triangulation) {
13
14 std::vector<std::tuple<double, int>> mappingsSortedPerValue;
15
16 double min = 0;
17 double max = 0;
18 int nbConstraints = 0;
19
20 int numberOfBytesRead = 0;
21 if(!ZFPOnly) {
22 numberOfBytesRead
23 += ReadPersistenceIndex(fm, mapping_, mappingsSortedPerValue,
24 criticalConstraints_, min, max, nbConstraints);
25 this->printMsg("Successfully read geomap.");
26 }
27
28 // Prepare array reconstruction.
29 int nx = 1 + dataExtent_[1] - dataExtent_[0];
30 int ny = 1 + dataExtent_[3] - dataExtent_[2];
31 int nz = 1 + dataExtent_[5] - dataExtent_[4];
32 int vertexNumber = nx * ny * nz;
33
34 decompressedData_.resize(vertexNumber);
35 if(ZFPTolerance < 0.0) {
36
37 // 2.a. (2.) Assign values to points thanks to topology indices.
38 for(int i = 0; i < vertexNumber; ++i) {
39 int seg = segmentation_[i];
40 auto end = mapping_.end();
41 auto it = std::lower_bound(
42 mapping_.begin(), mapping_.end(), std::make_tuple(0, seg), cmp);
43 if(it != end) {
44 std::tuple<double, int> tt = *it;
45 double value = std::get<0>(tt);
46 int sseg = std::get<1>(tt);
47 if(seg != sseg) {
48 this->printErr("Decompression mismatch (" + std::to_string(seg) + ", "
49 + std::to_string(sseg) + ")");
50 }
51 decompressedData_[i] = value;
52 } else {
53 this->printErr("Could not find " + std::to_string(seg) + " index.");
54 std::tuple<double, int> tt = *it;
55 double value = std::get<0>(tt);
56 decompressedData_[i] = value;
57 }
58 }
59
60 this->printMsg("Successfully assigned geomap.");
61
62 } else {
63#ifdef TTK_ENABLE_ZFP
64 // 2.b. (2.) Read with ZFP.
65 numberOfBytesRead
66 += CompressWithZFP(fm, true, decompressedData_, nx, ny, nz, ZFPTolerance);
67 this->printMsg("Successfully read with ZFP.");
68#else
69 this->printErr(
70 "Attempted to read with ZFP a ZFP block but ZFP is not installed.");
71 return -5;
72#endif
73 }
74
75 // No SQ.
76 if(SQMethodInt == 0 || SQMethodInt == 3) {
77 for(int i = 0; i < (int)criticalConstraints_.size(); ++i) {
78 std::tuple<int, double, int> t = criticalConstraints_[i];
79 int id = std::get<0>(t);
80 double val = std::get<1>(t);
81 decompressedData_[id] = val;
82 }
83 }
84
85 // double tolerance = Tolerance;
86 if(min == max) {
87 this->printWrn("Empty scalar field range.");
88 } else {
89 // tolerance *= (max - min);
90 }
91
92 if(SQMethodInt == 1 || SQMethodInt == 2)
93 return 0;
94
95 if(ZFPOnly)
96 return 0;
97
98 // 2.b. (3.) Crop whatever doesn't fit in topological intervals.
99 CropIntervals(mapping_, mappingsSortedPerValue, min, max, vertexNumber,
101 this->printMsg("Successfully cropped bad intervals.");
102
103 // 2.b. (4.) Apply topological simplification with min/max constraints
104 PerformSimplification<double>(criticalConstraints_, nbConstraints,
105 vertexNumber, decompressedData_.data(),
106 triangulation);
107 this->printMsg("Successfully performed simplification.");
108
109 this->rawFileLength += numberOfBytesRead;
110
111 return 0;
112}
113
114template <typename dataType, typename triangulationType>
115int ttk::TopologicalCompression::PerformSimplification(
116 const std::vector<std::tuple<int, double, int>> &constraints,
117 int nbConstraints,
118 int vertexNumber,
119 double *array,
120 const triangulationType &triangulation) {
121
122 std::vector<int> inputOffsets(vertexNumber);
123 std::vector<SimplexId> critConstraints(nbConstraints);
124 std::vector<double> inArray(vertexNumber);
125 decompressedOffsets_.resize(vertexNumber);
126 int status = 0;
127
128 // Offsets
129 for(int i = 0; i < vertexNumber; ++i)
130 inputOffsets[i] = i;
131
132 // Preprocess simplification.
133 for(int i = 0; i < nbConstraints; ++i) {
134 std::tuple<int, double, int> t = constraints[i];
135 int id = std::get<0>(t);
136 double val = std::get<1>(t);
137 int type = std::get<2>(t);
138
139 array[id] = val;
140
141 // Smoothe neighborhood (along with offsets).
142 SimplexId neighborNumber = triangulation.getVertexNeighborNumber(id);
143 for(SimplexId j = 0; j < neighborNumber; ++j) {
144 SimplexId neighbor{-1};
145 triangulation.getVertexNeighbor(id, j, neighbor);
146
147 if(type == 1) { // Local_maximum.
148 if(array[neighbor] > val)
149 array[neighbor] = val;
150 if(array[neighbor] == val
151 && inputOffsets[neighbor] > inputOffsets[id]) {
152 std::swap(inputOffsets[id], inputOffsets[neighbor]);
153 }
154 } else if(type == -1) { // Local_minimum.
155 if(array[neighbor] < val)
156 array[neighbor] = val;
157 if(array[neighbor] == val
158 && inputOffsets[neighbor] < inputOffsets[id]) {
159 std::swap(inputOffsets[id], inputOffsets[neighbor]);
160 }
161 } else if(type == 0) { // Saddle
162 }
163 }
164
165 critConstraints[i] = id;
166 }
167
168 for(int i = 0; i < vertexNumber; ++i) {
169 inArray[i] = array[i];
170 }
171 for(int i = 0; i < vertexNumber; ++i)
172 decompressedOffsets_[i] = 0;
173
174 std::vector<SimplexId> vertsOrder(vertexNumber);
175 sortVertices(vertexNumber, array, inputOffsets.data(), vertsOrder.data(),
176 this->threadNumber_);
177
178 status = topologicalSimplification.execute<double>(
179 inArray.data(), array, critConstraints.data(), vertsOrder.data(),
180 decompressedOffsets_.data(), nbConstraints, triangulation);
181
182 return status;
183}
184
185template <typename dataType>
187 std::vector<std::tuple<dataType, int>> &mappings,
188 std::vector<std::tuple<dataType, int>> &mappingsSortedPerValue,
189 double min,
190 double max,
191 int vertexNumber,
192 double *array,
193 std::vector<int> &segmentation) const {
194
195 int numberOfMisses = 0;
196 for(int i = 0; i < vertexNumber; ++i) {
197 int seg = segmentation[i];
198 auto end = mappings.end();
199 auto it = lower_bound(
200 mappings.begin(), mappings.end(), std::make_tuple(0, seg), cmp);
201 if(it != end) {
202 std::tuple<dataType, int> tt = *it;
203 double value = std::get<0>(tt);
204 int sseg = std::get<1>(tt);
205 if(seg != sseg) {
206 this->printErr("Decompression mismatch.");
207 }
208
209 auto it2 = lower_bound(mappingsSortedPerValue.begin(),
210 mappingsSortedPerValue.end(),
211 std::make_tuple(value, 0), cmp2);
212
213 if(it2 != mappingsSortedPerValue.end()
214 && it2 != mappingsSortedPerValue.begin()
215 && (it2 + 1) != mappingsSortedPerValue.end()) {
216 dataType vv0 = std::get<0>(*(it2 - 1));
217 // double vv1 = std::get<0>(*(it2));
218 dataType vv2 = std::get<0>(*(++it2));
219 // double m0 = (vv0 + vv1) / 2;
220 // double m1 = (vv1 + vv2) / 2;
221
222 if(array[i] < vv0) {
223 numberOfMisses++;
224 array[i] = vv0;
225 } else if(array[i] > vv2) {
226 numberOfMisses++;
227 array[i] = vv2;
228 }
229 } else {
230 if((it2 == mappingsSortedPerValue.end()
231 || (it2 + 1) == mappingsSortedPerValue.end())
232 && array[i] > max) {
233 numberOfMisses++;
234 array[i] = max;
235 } else if((it2 == mappingsSortedPerValue.begin()
236 || (it2 - 1) == mappingsSortedPerValue.begin())
237 && array[i] < min) {
238 numberOfMisses++;
239 array[i] = min;
240 }
241 }
242 } else {
243 this->printErr("Error looking for topo index.");
244 }
245 }
246
247 if(numberOfMisses > 0) {
248 this->printWrn("Missed " + std::to_string(numberOfMisses) + " values.");
249 }
250}
251
252template <typename dataType, typename triangulationType>
254 std::vector<std::tuple<SimplexId, SimplexId, dataType>> &JTPairs,
255 std::vector<std::tuple<SimplexId, SimplexId, dataType>> &STPairs,
256 const dataType *const inputScalars_,
257 const SimplexId *const inputOffsets,
258 const triangulationType &triangulation) {
259
260 // Compute offsets
261 const SimplexId numberOfVertices = triangulation.getNumberOfVertices();
262 std::vector<SimplexId> voffsets((unsigned long)numberOfVertices);
263 std::copy(inputOffsets, inputOffsets + numberOfVertices, voffsets.begin());
264
265 // Get contour tree
266 ftmTreePP.setVertexScalars(inputScalars_);
267 ftmTreePP.setTreeType(ftm::TreeType::Join_Split);
268 ftmTreePP.setVertexSoSoffsets(voffsets.data());
269 ftmTreePP.setThreadNumber(threadNumber_);
270 ftmTreePP.build<dataType>(&triangulation);
271 ftmTreePP.setSegmentation(false);
272 ftmTreePP.computePersistencePairs<dataType>(JTPairs, true);
273 ftmTreePP.computePersistencePairs<dataType>(STPairs, false);
274
275 return 0;
276}
277
278template <typename dataType, typename triangulationType>
280 int vertexNumber,
281 const dataType *const inputData,
282 const SimplexId *const inputOffsets,
283 dataType *outputData,
284 const double &tol,
285 const triangulationType &triangulation) {
286
287 ttk::Timer t;
288 ttk::Timer t1;
289
290 // 1. Compute persistence pairs.
291 std::vector<std::tuple<dataType, int>> topoIndices;
292
293 // Compute global min & max
294 dataType iter;
295 dataType maxValue = inputData[0];
296 int maxIndex = 0;
297 for(int i = 1; i < vertexNumber; ++i) {
298 iter = inputData[i];
299 if(iter > maxValue) {
300 maxIndex = i;
301 maxValue = iter;
302 }
303 }
304
305 dataType minValue = inputData[0];
306 int minIndex = 0;
307 for(int i = 1; i < vertexNumber; ++i) {
308 iter = inputData[i];
309 if(iter < minValue) {
310 minIndex = i;
311 minValue = iter;
312 }
313 }
314
315 topoIndices.push_back(std::make_tuple(maxValue, maxIndex));
316 topoIndices.push_back(std::make_tuple(minValue, minIndex));
317 double tolerance = 0.01 * tol * (maxValue - minValue);
318 double maxError = 0.01 * MaximumError * (maxValue - minValue);
319
320 this->printMsg(
321 "Computed min/max", 1.0, t.getElapsedTime(), this->threadNumber_);
322 t.reStart();
323
324 bool sqDomain = false;
325 bool sqRange = false;
326
327 const char *sq = SQMethod.c_str();
328 int nbCrit = 0;
329 std::vector<SimplexId> simplifiedConstraints;
330 if(strcmp(sq, "") == 0 && !ZFPOnly) {
331 // No SQ: perform topological control
332
333 std::vector<std::tuple<SimplexId, SimplexId, dataType>> JTPairs;
334 std::vector<std::tuple<SimplexId, SimplexId, dataType>> STPairs;
335 computePersistencePairs<dataType>(
336 JTPairs, STPairs, inputData, inputOffsets, triangulation);
337
338 this->printMsg("Computed persistence pairs", 1.0, t.getElapsedTime(),
339 this->threadNumber_);
340 t.reStart();
341
342 int nbJ = JTPairs.size();
343 int nbS = STPairs.size();
344 std::vector<int> critConstraints(2 * nbJ + 2 * nbS);
345
346 dataType maxEpsilon = 0;
347 dataType epsilonSum2 = 0;
348 dataType epsilonSum1 = 0;
349 dataType persistentSum2 = 0;
350 dataType persistentSum1 = 0;
351
352 // Join
353 for(int i = 0; i < nbJ; ++i) {
354 SimplexId cp1 = std::get<0>(JTPairs[i]);
355 SimplexId cp2 = std::get<1>(JTPairs[i]);
356 dataType idt1 = inputData[cp1];
357 dataType idt2 = inputData[cp2];
358 dataType p1 = std::max(idt2, idt1) - std::min(idt2, idt1);
359 if(p1 > tolerance) {
360 persistentSum2 += (p1 * p1);
361 persistentSum1 += abs<dataType>(p1);
362 int type1 = topologicalSimplification.getCriticalType(
363 cp1, inputOffsets, triangulation);
364 int type2 = topologicalSimplification.getCriticalType(
365 cp2, inputOffsets, triangulation);
366 if(type1 == 0) {
367 // authorizedSaddles->push_back(cp1);
368 topoIndices.push_back(std::make_tuple(idt1, cp1));
369 }
370 if(type2 == 0) {
371 // authorizedSaddles->push_back(cp2);
372 topoIndices.push_back(std::make_tuple(idt2, cp2));
373 }
374
375 nbCrit += 2;
376 critConstraints[2 * i] = cp1;
377 critConstraints[2 * i + 1] = cp2;
378 } else {
379 if(maxEpsilon < abs<dataType>(p1))
380 maxEpsilon = abs<dataType>(p1);
381 epsilonSum2 += (p1 * p1);
382 epsilonSum1 += abs<dataType>(p1);
383 critConstraints[2 * i] = -1;
384 critConstraints[2 * i + 1] = -1;
385 }
386 }
387
388 this->printMsg("Post-processed join pairs", 1.0, t.getElapsedTime(),
389 this->threadNumber_);
390 t.reStart();
391
392 // Split
393 for(int i = nbJ; i < nbJ + nbS; ++i) {
394 int si = i - nbJ;
395 SimplexId cp1 = std::get<0>(STPairs[si]);
396 SimplexId cp2 = std::get<1>(STPairs[si]);
397 dataType idt1 = inputData[cp1];
398 dataType idt2 = inputData[cp2];
399 dataType p1 = std::max(idt2, idt1) - std::min(idt2, idt1);
400 if(p1 > tolerance) {
401 persistentSum2 += (p1 * p1);
402 persistentSum1 += abs<dataType>(p1);
403 // Saddle selection.
404 int type1 = topologicalSimplification.getCriticalType(
405 cp1, inputOffsets, triangulation);
406 int type2 = topologicalSimplification.getCriticalType(
407 cp2, inputOffsets, triangulation);
408 if(type1 == 0) {
409 // authorizedSaddles->push_back(cp1);
410 topoIndices.push_back(std::make_tuple(idt1, cp1));
411 }
412 if(type2 == 0) {
413 // authorizedSaddles->push_back(cp2);
414 topoIndices.push_back(std::make_tuple(idt2, cp2));
415 }
416
417 nbCrit += 2;
418 critConstraints[2 * i] = cp1;
419 critConstraints[2 * i + 1] = cp2;
420 } else {
421 if(maxEpsilon < p1)
422 maxEpsilon = abs<dataType>(p1);
423 epsilonSum2 += (p1 * p1);
424 epsilonSum1 += abs<dataType>(p1);
425 critConstraints[2 * i] = -1;
426 critConstraints[2 * i + 1] = -1;
427 }
428 }
429
430 this->printMsg("Post-processed split pairs", 1.0, t.getElapsedTime(),
431 this->threadNumber_);
432 t.reStart();
433
434 simplifiedConstraints.resize(nbCrit);
435 {
436 int j = 0;
437 for(int i = 0; i < 2 * nbJ + 2 * nbS; ++i) {
438 int c = critConstraints[i];
439 if(c != -1)
440 simplifiedConstraints[j++] = c;
441 }
442 }
443
444 this->printMsg("Got pairs", 1.0, t.getElapsedTime(), this->threadNumber_);
445 t.reStart();
446
447 // 2. Perform topological simplification with constraints.
448 if(UseTopologicalSimplification) {
449 compressedOffsets_.resize(vertexNumber);
450 int status = 0;
451 status = topologicalSimplification.execute<dataType>(
452 inputData, outputData, simplifiedConstraints.data(), inputOffsets,
453 compressedOffsets_.data(), nbCrit, triangulation);
454 if(status != 0) {
455 return status;
456 }
457 }
458
459 this->printMsg(
460 "Performed simplification", 1.0, t.getElapsedTime(), this->threadNumber_);
461 t.reStart();
462
463 } else if(strcmp(sq, "r") == 0 || strcmp(sq, "R") == 0) {
464 // Range-based SQ: simple range quantization (automatic)
465 sqRange = true;
466 } else if(strcmp(sq, "d") == 0 || strcmp(sq, "D") == 0) {
467 // Domain-based SQ: range quantization + later domain control
468 sqDomain = true;
469 } else if(ZFPOnly) {
470 return 0;
471 } else {
472 this->printErr("Unrecognized SQ option (" + SQMethod + ").");
473 return -3;
474 }
475
476 // 3. Subdivision and attribution of a topological index
477
478 auto cmp = [](const std::tuple<dataType, int> &a,
479 const std::tuple<dataType, int> &b) {
480 return std::get<0>(a) < std::get<0>(b);
481 };
482 std::sort(topoIndices.begin(), topoIndices.end(), cmp);
483 std::vector<std::tuple<dataType, int>> segments;
484 // ith rank = ith segment (from bottom)
485 // 0: value
486 // 1: critical point index or -1 if !critical
487 auto l = (int)topoIndices.size();
488 if(l < 1) {
489 this->printMsg("Trivial subdivision performed.");
490 segments.push_back(topoIndices[l - 1]);
491 } else {
492 for(int i = 0; i < l; ++i) {
493 dataType v0 = std::get<0>(topoIndices[i]);
494 if(i == l - 1) {
495 segments.push_back(topoIndices[i]);
496 break;
497 }
498
499 dataType v1 = std::get<0>(topoIndices[i + 1]);
500 int i1 = std::get<1>(topoIndices[i + 1]);
501 auto diff = (double)(v1 - v0);
502
503 if(diff == 0)
504 continue;
505 if(!Subdivide || (diff < (dataType)maxError)) {
506 segments.push_back(std::make_tuple(v1, i1));
507 } else {
508 // Subdivide.
509 double nSegments = std::ceil(diff / maxError);
510 for(int j = 0, nbs = (int)nSegments; j < nbs; ++j) {
511 dataType sample = v0 + j * maxError;
512 int int1 = i1;
513 segments.push_back(std::make_tuple(sample, int1));
514 }
515 }
516 }
517 }
518
519 this->printMsg("Subdivision/topological index attribution", 1.0,
520 t.getElapsedTime(), this->threadNumber_);
521 t.reStart();
522
523 // 4. Affect segment value to all points.
524
525 // (pad buffer to avoid an out-of-bounds access/buffer overflow in
526 // WriteCompactSegmentation while loop)
527 segmentation_.resize(vertexNumber + 32);
528
529 std::sort(segments.begin(), segments.end(), cmp);
530 for(int i = 0; i < vertexNumber; ++i) {
531 dataType scalar = inputData[i];
532
533 auto begin = segments.begin();
534 auto end = segments.end();
535 auto it = std::lower_bound(
536 segments.begin(), segments.end(), std::make_tuple(scalar, -1), cmp);
537
538 if(it != end) {
539 std::tuple<dataType, int> tt = *it;
540 int j = it - begin;
541 int last = (int)segments.size() - 1;
542 if(j < last) {
543 dataType dtv = std::get<0>(tt);
544 if(j > 0) {
545 dtv = 0.5 * (dtv + std::get<0>(*(it - 1)));
546 }
547
548 outputData[i] = dtv;
549 int seg = j;
550 segmentation_[i] = seg;
551 } else {
552 segmentation_[i] = last;
553 outputData[i] = std::get<0>(tt); // maxValue;
554 }
555 }
556 }
557
558 this->printMsg(
559 "Assigned values", 1.0, t.getElapsedTime(), this->threadNumber_);
560 t.reStart();
561
562 // 4.1. Simplify mapping
563 auto segmentsSize = (int)segments.size();
564 std::vector<bool> affectedSegments;
565 affectedSegments.resize(segmentsSize);
566 std::vector<int> oob;
567 for(int i = 0; i < vertexNumber; ++i) {
568 int seg = segmentation_[i];
569 if(seg >= segmentsSize
570 && std::find(oob.begin(), oob.end(), seg) == oob.end())
571 oob.push_back(seg);
572 else if(seg >= 0)
573 affectedSegments[seg] = true;
574 else {
575 this->printErr("Negative segment encountered (" + std::to_string(i)
576 + ")");
577 }
578 }
579 std::vector<int> empty;
580 for(int i = 0; i < segmentsSize; ++i) {
581 if(!affectedSegments[i])
582 empty.push_back(i);
583 }
584
585 std::sort(oob.begin(), oob.end());
586 std::sort(empty.begin(), empty.end());
587 int indexLast = -1;
588 if(oob.size() > empty.size()) {
589 this->printWrn("oob size > empty size: " + std::to_string(oob.size()) + ", "
590 + std::to_string(empty.size()));
591 } else {
592 // Replace
593 for(int i = 0; i < vertexNumber; ++i) {
594 int seg = segmentation_[i];
595 if(seg >= segmentsSize) {
596 auto begin = oob.begin();
597 auto end = oob.end();
598 auto it = std::lower_bound(begin, end, seg);
599 if(it != end) {
600 int j = (int)(it - begin);
601 segmentation_[i] = empty[j];
602 affectedSegments[j] = true;
603 }
604 } else if(!affectedSegments[seg]) {
605 this->printErr("Something impossible happened");
606 }
607 }
608
609 if(empty.size() > oob.size()) {
610 std::vector<int> map2(segmentsSize, -1);
611 for(int i = 0; i < segmentsSize; ++i) {
612 if(affectedSegments[i])
613 continue;
614 for(int j = (segmentsSize - 1); j > i; --j) {
615 if(!affectedSegments[j] || map2[j] > 0)
616 continue;
617 map2[j] = i;
618 affectedSegments[j] = false;
619 affectedSegments[i] = true;
620 break;
621 }
622 }
623
624 bool doneAffecting = false;
625 for(int i = 0; i < segmentsSize; ++i) {
626 if(!affectedSegments[i]) {
627 doneAffecting = true;
628 } else {
629 if(doneAffecting) {
630 this->printWrn("Hole detected at " + std::to_string(i)
631 + "th segment.");
632 } else {
633 indexLast = i;
634 }
635 }
636 }
637
638 for(int i = 0; i < vertexNumber; ++i) {
639 int seg = segmentation_[i];
640 if(map2[seg] > 0)
641 segmentation_[i] = map2[seg];
642 }
643 }
644 }
645
646 this->printMsg(
647 "Simplified mapping", 1.0, t.getElapsedTime(), this->threadNumber_);
648 t.reStart();
649
650 // 5. Expose mapping.
651 std::vector<bool> already(vertexNumber);
652 for(int i = 0; i < vertexNumber; ++i)
653 already[i] = false; // init
654
655 for(int i = 0; i < vertexNumber; ++i) {
656 int vert = segmentation_[i];
657 if(!already[vert]) {
658 already[vert] = true;
659 dataType dttt = outputData[i];
660 mapping_.emplace_back((double)dttt, vert);
661 }
662 }
663
664 this->printMsg(
665 "Exposed mapping", 1.0, t.getElapsedTime(), this->threadNumber_);
666 t.reStart();
667
668 // 6. SQ-D correction step.
669 // Indices stay compact.
670 if(sqDomain && !sqRange) {
671 std::vector<bool> markedVertices(vertexNumber); // false by default
672 std::vector<bool> markedSegments(segmentsSize); // false by default
673 int newIndex = segmentsSize;
674
675 for(int i = 0; i < vertexNumber; ++i) {
676 // Skip processed vertices.
677 if(markedVertices[i])
678 continue;
679
680 int seg = segmentation_[i];
681 bool newSegment = markedSegments[seg];
682 dataType minNewSegment = inputData[i];
683 dataType maxNewSegment = minNewSegment;
684
685 if(!newSegment) {
686 markedSegments[seg] = true;
687 }
688
689 // Neighborhood search with seed = i
690 std::stack<int> s;
691 s.push(i);
692
693 while(!s.empty()) {
694 // Get next element.
695 int vertex = s.top();
696 s.pop();
697
698 // Mark vertex as processed.
699 markedVertices[vertex] = true;
700
701 // New region was detected.
702 if(newSegment) {
703 // Affect new segmentation id to all vertices in region.
704 segmentation_[vertex] = newIndex;
705 // Progressively compute min & max on local region.
706 dataType currentValue = inputData[vertex];
707 if(currentValue > maxNewSegment)
708 maxNewSegment = currentValue;
709 if(currentValue < minNewSegment)
710 minNewSegment = currentValue;
711 }
712
713 // Get neighbors.
714 SimplexId neighborNumber
715 = triangulation.getVertexNeighborNumber(vertex);
716 for(SimplexId j = 0; j < neighborNumber; ++j) {
717 SimplexId neighbor{-1};
718 triangulation.getVertexNeighbor(vertex, j, neighbor);
719
720 // Add current neighbor to processing stack.
721 if(!markedVertices[neighbor] && segmentation_[neighbor] == seg) {
722 s.push(neighbor);
723 }
724 }
725 }
726
727 // Affect a value to the new segmentation index.
728 if(newSegment) {
729 dataType result = (maxNewSegment + minNewSegment) / 2;
730 segments.push_back(std::make_tuple(result, newIndex));
731 mapping_.emplace_back((double)result, newIndex);
732 newIndex++; // Set index for next potential segment.
733 }
734 }
735
736 this->printMsg(
737 "SQ-D correction step", 1.0, t.getElapsedTime(), this->threadNumber_);
738 t.reStart();
739 }
740
741 // 7. [ZFP]: max constraints, min constraints
742 if(!sqDomain && !sqRange && !ZFPOnly) {
743 for(int i = 0; i < nbCrit; ++i) {
744 SimplexId id = simplifiedConstraints[i];
745 dataType val = inputData[id];
746 int type = topologicalSimplification.getCriticalType(
747 id, inputOffsets, triangulation);
748 if(type == -1 // Local_minimum
749 || type == 1 // Local_maximum
750 || type == 0) {
751 criticalConstraints_.emplace_back(id, (double)val, type);
752 } else { // 0 -> Saddle
753 }
754 }
755
756 this->printMsg(
757 "Exposed constraints", 1.0, t.getElapsedTime(), this->threadNumber_);
758 t.reStart();
759 }
760
761 {
762 if(indexLast > -1 && indexLast + 1 != (int)mapping_.size()) {
763 this->printWrn("Possible affectation mismatch ("
764 + std::to_string(indexLast + 1) + ", "
765 + std::to_string(mapping_.size()) + ")");
766 }
767
768 int nSegments = indexLast > -1 ? indexLast + 1 : (int)segments.size() - 1;
769 this->NbSegments = nSegments;
770 this->NbVertices = vertexNumber;
771
772 this->printMsg("Assigned " + std::to_string(nSegments) + " segment(s).",
773 1.0, t1.getElapsedTime(), this->threadNumber_);
774 }
775
776 return 0;
777}
int printWrn(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition: Debug.h:159
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
void reStart()
Definition: Timer.h:21
double getElapsedTime()
Definition: Timer.h:15
int ReadPersistenceIndex(FILE *fm, std::vector< std::tuple< double, int > > &mappings, std::vector< std::tuple< double, int > > &mappingsSortedPerValue, std::vector< std::tuple< int, double, int > > &constraints, double &min, double &max, int &nbConstraints) const
int compressForPersistenceDiagram(int vertexNumber, const dataType *const inputData, const SimplexId *const inputOffset, dataType *outputData, const double &tol, const triangulationType &triangulation)
int computePersistencePairs(std::vector< std::tuple< SimplexId, SimplexId, dataType > > &JTPairs, std::vector< std::tuple< SimplexId, SimplexId, dataType > > &STPairs, const dataType *const inputScalars_, const SimplexId *const inputOffsets, const triangulationType &triangulation)
static bool cmp(const std::tuple< double, int > &a, const std::tuple< double, int > &b)
std::vector< double > decompressedData_
std::vector< std::tuple< int, double, int > > criticalConstraints_
void CropIntervals(std::vector< std::tuple< dataType, int > > &mappings, std::vector< std::tuple< dataType, int > > &mappingsSortedPerValue, double min, double max, int vertexNumber, double *array, std::vector< int > &Seg) const
std::vector< std::tuple< double, int > > mapping_
void sortVertices(const size_t nVerts, const scalarType *const scalars, const idType *const offsets, SimplexId *const order, const int nThreads)
Sort vertices according to scalars disambiguated by offsets.
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22
printMsg(debug::output::GREEN+"                           "+debug::output::ENDCOLOR+debug::output::GREEN+"▒"+debug::output::ENDCOLOR+debug::output::GREEN+"▒▒▒▒▒▒▒▒▒▒▒▒▒░"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)