TTK
Loading...
Searching...
No Matches
DiscreteGradient_Template.h
Go to the documentation of this file.
1
14
15#pragma once
16
17#include <DiscreteGradient.h>
18#include <random>
19
20using ttk::SimplexId;
21using ttk::dcg::Cell;
22using ttk::dcg::CellExt;
23using ttk::dcg::DiscreteGradient;
24
25template <typename dataType, typename triangulationType>
26dataType DiscreteGradient::getPersistence(
27 const Cell &up,
28 const Cell &down,
29 const dataType *const scalars,
30 const triangulationType &triangulation) const {
31
32 return scalars[getCellGreaterVertex(up, triangulation)]
33 - scalars[getCellGreaterVertex(down, triangulation)];
34}
35
36template <typename triangulationType>
37int DiscreteGradient::buildGradient(const triangulationType &triangulation,
38 bool bypassCache,
39 const std::vector<bool> *updateMask) {
40
41 return buildGradient<double, triangulationType>(
42 triangulation, bypassCache, updateMask);
43}
44
45template <typename dataType, typename triangulationType>
46int DiscreteGradient::buildGradient(const triangulationType &triangulation,
47 bool bypassCache,
48 const std::vector<bool> *updateMask) {
49
50#ifdef TTK_ENABLE_MPI_TIME
51 ttk::Timer t_mpi;
52 ttk::startMPITimer(t_mpi, ttk::MPIrank_, ttk::MPIsize_);
53#endif
54
55 auto &cacheHandler = *triangulation.getGradientCacheHandler();
56 const auto findGradient
57 = [this, &cacheHandler]() -> AbstractTriangulation::gradientType * {
58 if(this->inputScalarField_.first == nullptr) {
59 return {};
60 }
61 return cacheHandler.get(this->inputScalarField_);
62 };
63
64#ifdef TTK_ENABLE_OPENMP
65 if(!bypassCache && omp_in_parallel()) {
66 this->printWrn(
67 "buildGradient() called inside a parallel region, disabling cache...");
68 bypassCache = true;
69 }
70#endif // TTK_ENABLE_OPENMP
71
72 // set member variables at each buildGradient() call
73 this->dimensionality_ = triangulation.getCellVertexNumber(0) - 1;
74 this->numberOfVertices_ = triangulation.getNumberOfVertices();
75
76 this->gradient_ = !bypassCache ? findGradient() : &this->localGradient_;
77 this->setReturnSaddleConnectors(bypassCache);
78 if(this->gradient_ == nullptr || bypassCache || this->newParameters()) {
79 if(this->gradient_ == nullptr && !bypassCache) {
80 // add new cache entry
81 cacheHandler.insert(this->inputScalarField_, {});
82 this->gradient_ = cacheHandler.get(this->inputScalarField_);
83 }
84
85 // allocate gradient memory
86 this->initMemory(triangulation);
87 Timer tm{};
88 if(updateMask) {
89 this->processLowerStarsWithMask(
90 this->inputOffsets_, triangulation, updateMask);
91 this->printMsg("Update cached discrete gradient", 1.0,
92 tm.getElapsedTime(), this->threadNumber_);
93 } else if(this->BackEnd == BACKEND::STOCHASTIC_BACKEND) {
94 this->processLowerStarsStochastic<dataType, triangulationType>(
95 this->inputOffsets_, triangulation);
96 this->printMsg("Built gradient (stochastic algorithm)", 1.0,
97 tm.getElapsedTime(), this->threadNumber_);
98
99 } else if(this->BackEnd == BACKEND::CLASSIC_BACKEND) {
100 this->processLowerStars(this->inputOffsets_, triangulation);
101 this->printMsg("Built gradient (homotopic expansion)", 1.0,
102 tm.getElapsedTime(), this->threadNumber_);
103#ifdef TTK_ENABLE_MPI_TIME
104 double elapsedTime
105 = ttk::endMPITimer(t_mpi, ttk::MPIrank_, ttk::MPIsize_);
106 if(ttk::MPIrank_ == 0) {
107 printMsg("Computation performed using " + std::to_string(ttk::MPIsize_)
108 + " MPI processes lasted: " + std::to_string(elapsedTime));
109 }
110#endif
111 }
112 } else {
113 this->printMsg("Fetched cached discrete gradient");
114 if(updateMask) {
115 Timer tm{};
116 this->processLowerStarsWithMask(
117 this->inputOffsets_, triangulation, updateMask);
118 this->printMsg("Update cached discrete gradient", 1.0,
119 tm.getElapsedTime(), this->threadNumber_);
120 }
121 }
122 return 0;
123}
124
125template <typename triangulationType>
126int DiscreteGradient::setCriticalPoints(
127 const std::array<std::vector<SimplexId>, 4> &criticalCellsByDim,
128 std::vector<std::array<float, 3>> &points,
129 std::vector<char> &cellDimensions,
130 std::vector<SimplexId> &cellIds,
131 std::vector<char> &isOnBoundary,
132 std::vector<SimplexId> &PLVertexIdentifiers,
133 const triangulationType &triangulation) const {
134
135 std::array<size_t, 5> partSums{};
136 for(size_t i = 0; i < criticalCellsByDim.size(); ++i) {
137 partSums[i + 1] = partSums[i] + criticalCellsByDim[i].size();
138 }
139
140 const auto nCritPoints = partSums.back();
141
142 points.resize(nCritPoints);
143 cellDimensions.resize(nCritPoints);
144 cellIds.resize(nCritPoints);
145 isOnBoundary.resize(nCritPoints);
146 PLVertexIdentifiers.resize(nCritPoints);
147
148 for(size_t i = 0; i < criticalCellsByDim.size(); ++i) {
149#ifdef TTK_ENABLE_OPENMP
150#pragma omp parallel for num_threads(threadNumber_)
151#endif // TTK_ENABLE_OPENMP
152 for(size_t j = 0; j < criticalCellsByDim[i].size(); ++j) {
153 const SimplexId cellId = criticalCellsByDim[i][j];
154 const int cellDim = i;
155 const auto o{partSums[i] + j};
156
157 triangulation.getCellIncenter(cellId, i, points[o].data());
158 cellDimensions[o] = cellDim;
159#ifdef TTK_ENABLE_MPI
160 ttk::SimplexId globalId{-1};
161 triangulation.getDistributedGlobalCellId(cellId, cellDim, globalId);
162 cellIds[o] = globalId;
163#else
164 cellIds[o] = cellId;
165#endif // TTK_ENABLE_MPI
166 const Cell cell{static_cast<int>(i), cellId};
167 isOnBoundary[o] = this->isBoundary(cell, triangulation);
168 PLVertexIdentifiers[o] = this->getCellGreaterVertex(cell, triangulation);
169 }
170 }
171
172 std::vector<std::vector<std::string>> rows(this->dimensionality_ + 1);
173 for(int i = 0; i < this->dimensionality_ + 1; ++i) {
174 rows[i]
175 = std::vector<std::string>{"#" + std::to_string(i) + "-cell(s)",
176 std::to_string(criticalCellsByDim[i].size())};
177 }
178 this->printMsg(rows);
179
180 return 0;
181}
182
183template <typename triangulationType>
184int DiscreteGradient::setCriticalPoints(
185 std::vector<std::array<float, 3>> &points,
186 std::vector<char> &cellDimensions,
187 std::vector<SimplexId> &cellIds,
188 std::vector<char> &isOnBoundary,
189 std::vector<SimplexId> &PLVertexIdentifiers,
190 const triangulationType &triangulation) const {
191
192 std::array<std::vector<SimplexId>, 4> criticalCellsByDim;
193 getCriticalPoints(criticalCellsByDim, triangulation);
194 setCriticalPoints(criticalCellsByDim, points, cellDimensions, cellIds,
195 isOnBoundary, PLVertexIdentifiers, triangulation);
196
197 return 0;
198}
199
200template <typename triangulationType>
201int DiscreteGradient::getCriticalPoints(
202 std::array<std::vector<SimplexId>, 4> &criticalCellsByDim,
203 const triangulationType &triangulation) const {
204
205 const auto dims{this->getNumberOfDimensions()};
206 for(int i = 0; i < dims; ++i) {
207
208 // map: store critical cell per dimension per thread
209 std::vector<std::vector<SimplexId>> critCellsPerThread(this->threadNumber_);
210 const auto numberOfCells{this->getNumberOfCells(i, triangulation)};
211
212 // use static scheduling to ensure that critical cells
213 // are sorted by id
214
215#ifdef TTK_ENABLE_OPENMP
216#pragma omp parallel for num_threads(this->threadNumber_) schedule(static)
217#endif // TTK_ENABLE_OPENMP
218 for(SimplexId j = 0; j < numberOfCells; ++j) {
219#ifdef TTK_ENABLE_OPENMP
220 const auto tid = omp_get_thread_num();
221#else
222 const auto tid = 0;
223#endif // TTK_ENABLE_OPENMP
224 if(this->isCellCritical(i, j)) {
225 // Only non-ghost critical simplices are taken into consideration
226#ifdef TTK_ENABLE_MPI
227 if(triangulation.getSimplexRank(j, i) == ttk::MPIrank_)
228#endif
229 critCellsPerThread[tid].emplace_back(j);
230 }
231 }
232
233 // reduce: aggregate critical cells per thread
234 criticalCellsByDim[i] = std::move(critCellsPerThread[0]);
235 for(size_t j = 1; j < critCellsPerThread.size(); ++j) {
236 const auto &vec{critCellsPerThread[j]};
237 criticalCellsByDim[i].insert(
238 criticalCellsByDim[i].end(), vec.begin(), vec.end());
239 }
240 }
241
242 return 0;
243}
244
245template <typename triangulationType>
246int DiscreteGradient::getCriticalPoints(
247 std::vector<Cell> &criticalPoints,
248 const triangulationType &triangulation) const {
249
250 // foreach dimension
251 const int numberOfDimensions = getNumberOfDimensions();
252 for(int i = 0; i < numberOfDimensions; ++i) {
253
254 // foreach cell of that dimension
255 const SimplexId numberOfCells = getNumberOfCells(i, triangulation);
256 for(SimplexId j = 0; j < numberOfCells; ++j) {
257 const Cell cell(i, j);
258
259 if(isCellCritical(cell)) {
260 criticalPoints.push_back(cell);
261 }
262 }
263 }
264
265 return 0;
266}
267
268template <typename triangulationType>
269SimplexId DiscreteGradient::getNumberOfCells(
270 const int dimension, const triangulationType &triangulation) const {
271
272 if(dimension > this->dimensionality_ || dimension < 0) {
273 return -1;
274 }
275
276 switch(dimension) {
277 case 0:
278 return triangulation.getNumberOfVertices();
279 break;
280
281 case 1:
282 return triangulation.getNumberOfEdges();
283 break;
284
285 case 2:
286 return triangulation.getNumberOfTriangles();
287 break;
288
289 case 3:
290 return triangulation.getNumberOfCells();
291 break;
292 }
293
294 return -1;
295}
296
297template <typename triangulationType>
298inline void DiscreteGradient::lowerStarWithMask(
299 lowerStarType &ls,
300 const SimplexId a,
301 const SimplexId *const offsets,
302 const triangulationType &triangulation,
303 const std::vector<bool> *updateMask) const {
304
305 // make sure that ls is cleared
306 for(auto &vec : ls) {
307 vec.clear();
308 }
309
310 // a belongs to its lower star
311 ls[0].emplace_back(CellExt{0, a});
312
313 if(updateMask != nullptr) {
314 for(auto &vec : ls[0]) {
315 int vertexId = vec.id_;
316 int edgeId = (*gradient_)[0][vertexId];
317
318 (*gradient_)[0][vertexId] = -1;
319 if(edgeId != -1) {
320 (*gradient_)[1][edgeId] = -1;
321 }
322 }
323 }
324
325 // store lower edges
326 const auto nedges = triangulation.getVertexEdgeNumber(a);
327 ls[1].reserve(nedges);
328 for(SimplexId i = 0; i < nedges; i++) {
329 SimplexId edgeId;
330 triangulation.getVertexEdge(a, i, edgeId);
331 SimplexId vertexId;
332 triangulation.getEdgeVertex(edgeId, 0, vertexId);
333 if(vertexId == a) {
334 triangulation.getEdgeVertex(edgeId, 1, vertexId);
335 }
336 if(offsets[vertexId] < offsets[a]) {
337 ls[1].emplace_back(CellExt{1, edgeId, {offsets[vertexId], -1, -1}, {}});
338 }
339 }
340
341 if(updateMask != nullptr) {
342 for(auto &vec : ls[1]) {
343 int edgeId = vec.id_;
344 int vertexId = (*gradient_)[1][edgeId];
345 int triangleId = (*gradient_)[2][edgeId];
346
347 (*gradient_)[1][edgeId] = -1;
348 if(vertexId != -1) {
349 (*gradient_)[0][vertexId] = -1;
350 }
351
352 (*gradient_)[2][edgeId] = -1;
353 if(triangleId != -1) {
354 (*gradient_)[3][triangleId] = -1;
355 }
356 }
357 }
358
359 if(ls[1].size() < 2) {
360 // at least two edges in the lower star for one triangle
361 return;
362 }
363
364 const auto processTriangle
365 = [&](const SimplexId triangleId, const SimplexId v0, const SimplexId v1,
366 const SimplexId v2) {
367 std::array<SimplexId, 3> lowVerts{-1, -1, -1};
368 if(v0 == a) {
369 lowVerts[0] = offsets[v1];
370 lowVerts[1] = offsets[v2];
371 } else if(v1 == a) {
372 lowVerts[0] = offsets[v0];
373 lowVerts[1] = offsets[v2];
374 } else if(v2 == a) {
375 lowVerts[0] = offsets[v0];
376 lowVerts[1] = offsets[v1];
377 }
378 // higher order vertex first
379 if(lowVerts[0] < lowVerts[1]) {
380 std::swap(lowVerts[0], lowVerts[1]);
381 }
382 if(offsets[a] > lowVerts[0]) { // triangle in lowerStar
383 uint8_t j{}, k{};
384 // store edges indices of current triangle
385 std::array<uint8_t, 3> faces{};
386 for(const auto &e : ls[1]) {
387 if(e.lowVerts_[0] == lowVerts[0] || e.lowVerts_[0] == lowVerts[1]) {
388 faces[k++] = j;
389 }
390 j++;
391 }
392 ls[2].emplace_back(CellExt{2, triangleId, lowVerts, faces});
393 }
394 };
395
396 if(dimensionality_ == 2) {
397 // store lower triangles
398
399 // use optimised triangulation methods:
400 // getVertexStar instead of getVertexTriangle
401 // getCellVertex instead of getTriangleVertex
402 const auto ncells = triangulation.getVertexStarNumber(a);
403 ls[2].reserve(ncells);
404 for(SimplexId i = 0; i < ncells; ++i) {
405 SimplexId cellId;
406 triangulation.getVertexStar(a, i, cellId);
407 SimplexId v0{}, v1{}, v2{};
408 triangulation.getCellVertex(cellId, 0, v0);
409 triangulation.getCellVertex(cellId, 1, v1);
410 triangulation.getCellVertex(cellId, 2, v2);
411 processTriangle(cellId, v0, v1, v2);
412 }
413
414 if(updateMask != nullptr) {
415 for(auto &vec : ls[2]) {
416 int triangleId = vec.id_;
417 int edgeId = (*gradient_)[3][triangleId];
418
419 (*gradient_)[3][triangleId] = -1;
420 if(edgeId != -1) {
421 (*gradient_)[2][edgeId] = -1;
422 }
423 }
424 }
425
426 } else if(dimensionality_ == 3) {
427 // store lower triangles
428 const auto ntri = triangulation.getVertexTriangleNumber(a);
429 ls[2].reserve(ntri);
430 for(SimplexId i = 0; i < ntri; i++) {
431 SimplexId triangleId;
432 triangulation.getVertexTriangle(a, i, triangleId);
433 SimplexId v0{}, v1{}, v2{};
434 triangulation.getTriangleVertex(triangleId, 0, v0);
435 triangulation.getTriangleVertex(triangleId, 1, v1);
436 triangulation.getTriangleVertex(triangleId, 2, v2);
437 processTriangle(triangleId, v0, v1, v2);
438 }
439
440 if(updateMask != nullptr) {
441 for(auto &vec : ls[2]) {
442 int triangleId = vec.id_;
443 int edgeId = (*gradient_)[3][triangleId];
444 int tetraId = (*gradient_)[4][triangleId];
445
446 (*gradient_)[3][triangleId] = -1;
447 if(edgeId != -1) {
448 (*gradient_)[2][edgeId] = -1;
449 }
450
451 (*gradient_)[4][triangleId] = -1;
452 if(tetraId != -1) {
453 (*gradient_)[5][tetraId] = -1;
454 }
455 }
456 }
457
458 // at least three triangles in the lower star for one tetra
459 if(ls[2].size() >= 3) {
460 // store lower tetra
461 const auto ncells = triangulation.getVertexStarNumber(a);
462 ls[3].reserve(ncells);
463 for(SimplexId i = 0; i < ncells; ++i) {
464 SimplexId cellId;
465 triangulation.getVertexStar(a, i, cellId);
466 std::array<SimplexId, 3> lowVerts{-1, -1, -1};
467 SimplexId v0{}, v1{}, v2{}, v3{};
468 triangulation.getCellVertex(cellId, 0, v0);
469 triangulation.getCellVertex(cellId, 1, v1);
470 triangulation.getCellVertex(cellId, 2, v2);
471 triangulation.getCellVertex(cellId, 3, v3);
472 if(v0 == a) {
473 lowVerts[0] = offsets[v1];
474 lowVerts[1] = offsets[v2];
475 lowVerts[2] = offsets[v3];
476 } else if(v1 == a) {
477 lowVerts[0] = offsets[v0];
478 lowVerts[1] = offsets[v2];
479 lowVerts[2] = offsets[v3];
480 } else if(v2 == a) {
481 lowVerts[0] = offsets[v0];
482 lowVerts[1] = offsets[v1];
483 lowVerts[2] = offsets[v3];
484 } else if(v3 == a) {
485 lowVerts[0] = offsets[v0];
486 lowVerts[1] = offsets[v1];
487 lowVerts[2] = offsets[v2];
488 }
489 if(offsets[a] > *std::max_element(
490 lowVerts.begin(), lowVerts.end())) { // tetra in lowerStar
491
492 // higher order vertex first
493 std::sort(lowVerts.rbegin(), lowVerts.rend());
494
495 uint8_t j{}, k{};
496 // store triangles indices of current tetra
497 std::array<uint8_t, 3> faces{};
498 for(const auto &t : ls[2]) {
499 // lowVerts & t.lowVerts are ordered, no need to check if
500 // t.lowVerts[0] == lowVerts[2] or t.lowVerts[1] == lowVerts[0]
501 if((t.lowVerts_[0] == lowVerts[0]
502 && (t.lowVerts_[1] == lowVerts[1]
503 || t.lowVerts_[1] == lowVerts[2]))
504 || (t.lowVerts_[0] == lowVerts[1]
505 && t.lowVerts_[1] == lowVerts[2])) {
506 faces[k++] = j;
507 }
508 j++;
509 }
510
511 ls[3].emplace_back(CellExt{3, cellId, lowVerts, faces});
512 }
513 }
514
515 if(updateMask != nullptr) {
516 for(auto &vec : ls[3]) {
517 int tetraId = vec.id_;
518 int triangleId = (*gradient_)[5][tetraId];
519 (*gradient_)[5][tetraId] = -1;
520 if(triangleId != -1) {
521 (*gradient_)[4][triangleId] = -1;
522 }
523 }
524 }
525 }
526 }
527}
528
529template <typename triangulationType>
530inline void
531 DiscreteGradient::lowerStar(lowerStarType &ls,
532 const SimplexId a,
533 const SimplexId *const offsets,
534 const triangulationType &triangulation) const {
535
536 // WARNING
537 // If you modify this function, please make sure to also report your edit to
538 // the other implementation of this function, `lowerStarWithMask`.
539
540 // make sure that ls is cleared
541 for(auto &vec : ls) {
542 vec.clear();
543 }
544
545 // a belongs to its lower star
546 CellExt const localCellExt{0, a};
547 ls[0].emplace_back(localCellExt);
548
549 // store lower edges
550 const auto nedges = triangulation.getVertexEdgeNumber(a);
551 ls[1].reserve(nedges);
552 for(SimplexId i = 0; i < nedges; i++) {
553 SimplexId edgeId;
554 triangulation.getVertexEdge(a, i, edgeId);
555 SimplexId vertexId;
556 triangulation.getEdgeVertex(edgeId, 0, vertexId);
557 if(vertexId == a) {
558 triangulation.getEdgeVertex(edgeId, 1, vertexId);
559 }
560 if(offsets[vertexId] < offsets[a]) {
561 ls[1].emplace_back(CellExt{1, edgeId, {offsets[vertexId], -1, -1}, {}});
562 }
563 }
564
565 if(ls[1].size() < 2) {
566 // at least two edges in the lower star for one triangle
567 return;
568 }
569
570 const auto processTriangle
571 = [&](const SimplexId triangleId, const SimplexId v0, const SimplexId v1,
572 const SimplexId v2) {
573 std::array<SimplexId, 3> lowVerts{-1, -1, -1};
574 if(v0 == a) {
575 lowVerts[0] = offsets[v1];
576 lowVerts[1] = offsets[v2];
577 } else if(v1 == a) {
578 lowVerts[0] = offsets[v0];
579 lowVerts[1] = offsets[v2];
580 } else if(v2 == a) {
581 lowVerts[0] = offsets[v0];
582 lowVerts[1] = offsets[v1];
583 }
584 // higher order vertex first
585 if(lowVerts[0] < lowVerts[1]) {
586 std::swap(lowVerts[0], lowVerts[1]);
587 }
588 if(offsets[a] > lowVerts[0]) { // triangle in lowerStar
589 uint8_t j{}, k{};
590 // store edges indices of current triangle
591 std::array<uint8_t, 3> faces{};
592 for(const auto &e : ls[1]) {
593 if(e.lowVerts_[0] == lowVerts[0] || e.lowVerts_[0] == lowVerts[1]) {
594 faces[k++] = j;
595 }
596 j++;
597 }
598 CellExt const localCellExt2{2, triangleId, lowVerts, faces};
599 ls[2].emplace_back(localCellExt2);
600 }
601 };
602
603 if(dimensionality_ == 2) {
604 // store lower triangles
605
606 // use optimised triangulation methods:
607 // getVertexStar instead of getVertexTriangle
608 // getCellVertex instead of getTriangleVertex
609 const auto ncells = triangulation.getVertexStarNumber(a);
610 ls[2].reserve(ncells);
611 for(SimplexId i = 0; i < ncells; ++i) {
612 SimplexId cellId;
613 triangulation.getVertexStar(a, i, cellId);
614 SimplexId v0{}, v1{}, v2{};
615 triangulation.getCellVertex(cellId, 0, v0);
616 triangulation.getCellVertex(cellId, 1, v1);
617 triangulation.getCellVertex(cellId, 2, v2);
618 processTriangle(cellId, v0, v1, v2);
619 }
620 } else if(dimensionality_ == 3) {
621 // store lower triangles
622 const auto ntri = triangulation.getVertexTriangleNumber(a);
623 ls[2].reserve(ntri);
624 for(SimplexId i = 0; i < ntri; i++) {
625 SimplexId triangleId;
626 triangulation.getVertexTriangle(a, i, triangleId);
627 SimplexId v0{}, v1{}, v2{};
628 triangulation.getTriangleVertex(triangleId, 0, v0);
629 triangulation.getTriangleVertex(triangleId, 1, v1);
630 triangulation.getTriangleVertex(triangleId, 2, v2);
631 processTriangle(triangleId, v0, v1, v2);
632 }
633
634 // at least three triangles in the lower star for one tetra
635 if(ls[2].size() >= 3) {
636 // store lower tetra
637 const auto ncells = triangulation.getVertexStarNumber(a);
638 ls[3].reserve(ncells);
639 for(SimplexId i = 0; i < ncells; ++i) {
640 SimplexId cellId;
641 triangulation.getVertexStar(a, i, cellId);
642 std::array<SimplexId, 3> lowVerts{-1, -1, -1};
643 SimplexId v0{}, v1{}, v2{}, v3{};
644 triangulation.getCellVertex(cellId, 0, v0);
645 triangulation.getCellVertex(cellId, 1, v1);
646 triangulation.getCellVertex(cellId, 2, v2);
647 triangulation.getCellVertex(cellId, 3, v3);
648 if(v0 == a) {
649 lowVerts[0] = offsets[v1];
650 lowVerts[1] = offsets[v2];
651 lowVerts[2] = offsets[v3];
652 } else if(v1 == a) {
653 lowVerts[0] = offsets[v0];
654 lowVerts[1] = offsets[v2];
655 lowVerts[2] = offsets[v3];
656 } else if(v2 == a) {
657 lowVerts[0] = offsets[v0];
658 lowVerts[1] = offsets[v1];
659 lowVerts[2] = offsets[v3];
660 } else if(v3 == a) {
661 lowVerts[0] = offsets[v0];
662 lowVerts[1] = offsets[v1];
663 lowVerts[2] = offsets[v2];
664 }
665 if(offsets[a] > *std::max_element(
666 lowVerts.begin(), lowVerts.end())) { // tetra in lowerStar
667
668 // higher order vertex first
669 std::sort(lowVerts.rbegin(), lowVerts.rend());
670
671 uint8_t j{}, k{};
672 // store triangles indices of current tetra
673 std::array<uint8_t, 3> faces{};
674 for(const auto &t : ls[2]) {
675 // lowVerts & t.lowVerts are ordered, no need to check if
676 // t.lowVerts[0] == lowVerts[2] or t.lowVerts[1] == lowVerts[0]
677 if((t.lowVerts_[0] == lowVerts[0]
678 && (t.lowVerts_[1] == lowVerts[1]
679 || t.lowVerts_[1] == lowVerts[2]))
680 || (t.lowVerts_[0] == lowVerts[1]
681 && t.lowVerts_[1] == lowVerts[2])) {
682 faces[k++] = j;
683 }
684 j++;
685 }
686
687 CellExt const localCellExt3{3, cellId, lowVerts, faces};
688 ls[3].emplace_back(localCellExt3);
689 }
690 }
691 }
692 }
693}
694
695template <typename triangulationType>
696inline void DiscreteGradient::pairCells(
697 CellExt &alpha, CellExt &beta, const triangulationType &triangulation) {
698#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
699 char localBId{0}, localAId{0};
700 SimplexId a{}, b{};
701
702 if(beta.dim_ == 1) {
703
704 for(SimplexId i = 0; i < 2; ++i) {
705 triangulation.getEdgeVertex(beta.id_, i, a);
706 if(a == alpha.id_) {
707 localAId = i;
708 break;
709 }
710 }
711 const auto nedges = triangulation.getVertexEdgeNumber(alpha.id_);
712 for(SimplexId i = 0; i < nedges; ++i) {
713 triangulation.getVertexEdge(alpha.id_, i, b);
714 if(b == beta.id_) {
715 localBId = i;
716 break;
717 }
718 }
719 } else if(beta.dim_ == 2) {
720 for(SimplexId i = 0; i < 3; ++i) {
721 triangulation.getTriangleEdge(beta.id_, i, a);
722 if(a == alpha.id_) {
723 localAId = i;
724 break;
725 }
726 }
727 const auto ntri = triangulation.getEdgeTriangleNumber(alpha.id_);
728 for(SimplexId i = 0; i < ntri; ++i) {
729 triangulation.getEdgeTriangle(alpha.id_, i, b);
730 if(b == beta.id_) {
731 localBId = i;
732 break;
733 }
734 }
735 } else {
736 for(SimplexId i = 0; i < 4; ++i) {
737 triangulation.getCellTriangle(beta.id_, i, a);
738 if(a == alpha.id_) {
739 localAId = i;
740 break;
741 }
742 }
743 const auto ntetra = triangulation.getTriangleStarNumber(alpha.id_);
744 for(SimplexId i = 0; i < ntetra; ++i) {
745 triangulation.getTriangleStar(alpha.id_, i, b);
746 if(b == beta.id_) {
747 localBId = i;
748 break;
749 }
750 }
751 }
752 (*gradient_)[2 * alpha.dim_][alpha.id_] = localBId;
753 (*gradient_)[2 * alpha.dim_ + 1][beta.id_] = localAId;
754#else
755 TTK_FORCE_USE(triangulation);
756 (*gradient_)[2 * alpha.dim_][alpha.id_] = beta.id_;
757 (*gradient_)[2 * alpha.dim_ + 1][beta.id_] = alpha.id_;
758#endif // TTK_ENABLE_DCG_OPTIMIZE_MEMORY
759 alpha.paired_ = true;
760 beta.paired_ = true;
761}
762
763template <typename triangulationType>
764int DiscreteGradient::processLowerStars(
765 const SimplexId *const offsets, const triangulationType &triangulation) {
766
767 // WARNING
768 // If you modify this function, please make sure to also report your edit to
769 // the other implementation of this function, `processLowerStarsWithMask`.
770
771 /* Compute gradient */
772
773 auto nverts = triangulation.getNumberOfVertices();
774
775 // Comparison function for Cells inside priority queues
776 const auto orderCells = [&](const CellExt &a, const CellExt &b) -> bool {
777 return a.lowVerts_ > b.lowVerts_;
778 };
779
780 // Type alias for priority queues
781 using pqType
782 = std::priority_queue<std::reference_wrapper<CellExt>,
783 std::vector<std::reference_wrapper<CellExt>>,
784 decltype(orderCells)>;
785
786 // To reduce allocations, priority queues and lowerStar objects are
787 // cleaned & reused between iterations.
788
789 // Priority queues are pushed at the beginning and popped at the
790 // end. To pop the minimum, elements should be sorted in a
791 // decreasing order.
792 pqType pqZero{orderCells}, pqOne{orderCells};
793
794 // store lower star structure
795 lowerStarType Lx;
796
797 // In case the vertex is a ghost, the gradient is computed but may result in a
798 // false pairing if the ghost vertex belongs to the second layer of ghosts. To
799 // only produce correct pairing, one can add the following test: If the vertex
800 // is a ghost and all of its neighbor vertices are also ghosts, then the
801 // computation should not be done.
802
803#ifdef TTK_ENABLE_OPENMP
804#pragma omp parallel for num_threads(threadNumber_) \
805 firstprivate(Lx, pqZero, pqOne)
806#endif // TTK_ENABLE_OPENMP
807 for(SimplexId x = 0; x < nverts; x++) {
808
809 // clear priority queues (they should be empty at the end of the
810 // previous iteration)
811 while(!pqZero.empty()) {
812 pqZero.pop();
813 }
814 while(!pqOne.empty()) {
815 pqOne.pop();
816 }
817
818 // Insert into pqOne cofacets of cell c_alpha such as numUnpairedFaces == 1
819 const auto insertCofacets = [&](const CellExt &ca, lowerStarType &ls) {
820 if(ca.dim_ == 1) {
821 for(auto &beta : ls[2]) {
822 if(ls[1][beta.faces_[0]].id_ == ca.id_
823 || ls[1][beta.faces_[1]].id_ == ca.id_) {
824 // edge ca belongs to triangle beta
825 if(numUnpairedFacesTriangle(beta, ls).first == 1) {
826 pqOne.push(beta);
827 }
828 }
829 }
830
831 } else if(ca.dim_ == 2) {
832 for(auto &beta : ls[3]) {
833 if(ls[2][beta.faces_[0]].id_ == ca.id_
834 || ls[2][beta.faces_[1]].id_ == ca.id_
835 || ls[2][beta.faces_[2]].id_ == ca.id_) {
836 // triangle ca belongs to tetra beta
837 if(numUnpairedFacesTetra(beta, ls).first == 1) {
838 pqOne.push(beta);
839 }
840 }
841 }
842 }
843 };
844
845 lowerStar(Lx, x, offsets, triangulation);
846
847 // Lx[1] empty => x is a local minimum
848 if(!Lx[1].empty()) {
849 // get delta: 1-cell (edge) with minimal G value (steeper gradient)
850 size_t minId = 0;
851 for(size_t i = 1; i < Lx[1].size(); ++i) {
852 const auto &a = Lx[1][minId].lowVerts_[0];
853 const auto &b = Lx[1][i].lowVerts_[0];
854 if(a > b) {
855 // edge[i] < edge[0]
856 minId = i;
857 }
858 }
859
860 auto &c_delta = Lx[1][minId];
861
862 // store x (0-cell) -> delta (1-cell) V-path
863 pairCells(Lx[0][0], c_delta, triangulation);
864
865 // push every 1-cell in Lx that is not delta into pqZero
866 for(auto &alpha : Lx[1]) {
867 if(alpha.id_ != c_delta.id_) {
868 pqZero.push(alpha);
869 }
870 }
871
872 // push into pqOne every coface of delta in Lx (2-cells only,
873 // 3-cells have not any facet paired yet) such that
874 // numUnpairedFaces == 1
875 insertCofacets(c_delta, Lx);
876
877 while(!pqOne.empty() || !pqZero.empty()) {
878 while(!pqOne.empty()) {
879 auto &c_alpha = pqOne.top().get();
880 pqOne.pop();
881 auto unpairedFaces = numUnpairedFaces(c_alpha, Lx);
882 if(unpairedFaces.first == 0) {
883 pqZero.push(c_alpha);
884 } else {
885 auto &c_pair_alpha = Lx[c_alpha.dim_ - 1][unpairedFaces.second];
886
887 // store (pair_alpha) -> (alpha) V-path
888 pairCells(c_pair_alpha, c_alpha, triangulation);
889
890 // add cofaces of c_alpha and c_pair_alpha to pqOne
891 insertCofacets(c_alpha, Lx);
892 insertCofacets(c_pair_alpha, Lx);
893 }
894 }
895
896 // skip pair_alpha from pqZero:
897 // cells in pqZero are not critical if already paired
898 while(!pqZero.empty() && pqZero.top().get().paired_) {
899 pqZero.pop();
900 }
901
902 if(!pqZero.empty()) {
903 auto &c_gamma = pqZero.top().get();
904 pqZero.pop();
905
906 // gamma is a critical cell
907 // mark gamma as paired
908 c_gamma.paired_ = true;
909
910 // add cofacets of c_gamma to pqOne
911 insertCofacets(c_gamma, Lx);
912 }
913 }
914 }
915 }
916
917 return 0;
918}
919
920template <typename triangulationType>
921int DiscreteGradient::processLowerStarsWithMask(
922 const SimplexId *const offsets,
923 const triangulationType &triangulation,
924 const std::vector<bool> *updateMask) {
925
926 auto nverts = triangulation.getNumberOfVertices();
927
928 /* Compute gradient */
929
930 // Comparison function for Cells inside priority queues
931 const auto orderCells = [&](const CellExt &a, const CellExt &b) -> bool {
932 return a.lowVerts_ > b.lowVerts_;
933 };
934
935 // Type alias for priority queues
936 using pqType
937 = std::priority_queue<std::reference_wrapper<CellExt>,
938 std::vector<std::reference_wrapper<CellExt>>,
939 decltype(orderCells)>;
940
941 // To reduce allocations, priority queues and lowerStar objects are
942 // cleaned & reused between iterations.
943
944 // Priority queues are pushed at the beginning and popped at the
945 // end. To pop the minimum, elements should be sorted in a
946 // decreasing order.
947 pqType pqZero{orderCells}, pqOne{orderCells};
948
949 // store lower star structure
950 lowerStarType Lx;
951
952#ifdef TTK_ENABLE_OPENMP
953#pragma omp parallel for num_threads(threadNumber_) \
954 firstprivate(Lx, pqZero, pqOne)
955#endif // TTK_ENABLE_OPENMP
956 for(SimplexId x = 0; x < nverts; x++) {
957
958 if((updateMask != nullptr) && ((*updateMask)[x] == false)) {
959 continue;
960 }
961
962 // clear priority queues (they should be empty at the end of the
963 // previous iteration)
964 while(!pqZero.empty()) {
965 pqZero.pop();
966 }
967 while(!pqOne.empty()) {
968 pqOne.pop();
969 }
970
971 // Insert into pqOne cofacets of cell c_alpha such as numUnpairedFaces == 1
972 const auto insertCofacets = [&](const CellExt &ca, lowerStarType &ls) {
973 if(ca.dim_ == 1) {
974 for(auto &beta : ls[2]) {
975 if(ls[1][beta.faces_[0]].id_ == ca.id_
976 || ls[1][beta.faces_[1]].id_ == ca.id_) {
977 // edge ca belongs to triangle beta
978 if(numUnpairedFacesTriangle(beta, ls).first == 1) {
979 pqOne.push(beta);
980 }
981 }
982 }
983
984 } else if(ca.dim_ == 2) {
985 for(auto &beta : ls[3]) {
986 if(ls[2][beta.faces_[0]].id_ == ca.id_
987 || ls[2][beta.faces_[1]].id_ == ca.id_
988 || ls[2][beta.faces_[2]].id_ == ca.id_) {
989 // triangle ca belongs to tetra beta
990 if(numUnpairedFacesTetra(beta, ls).first == 1) {
991 pqOne.push(beta);
992 }
993 }
994 }
995 }
996 };
997
998 lowerStarWithMask(Lx, x, offsets, triangulation, updateMask);
999 // In case the vertex is a ghost, the gradient of the
1000 // simplices of its star is set to GHOST_GRADIENT
1001#ifdef TTK_ENABLE_MPI
1002 if(ttk::isRunningWithMPI()
1003 && triangulation.getVertexRank(x) != ttk::MPIrank_) {
1004 int sizeDim = Lx.size();
1005 for(int i = 0; i < sizeDim; i++) {
1006 int nCells = Lx[i].size();
1007 for(int j = 0; j < nCells; j++) {
1008 setCellToGhost(Lx[i][j].dim_, Lx[i][j].id_);
1009 }
1010 }
1011 } else
1012#endif // TTK_ENABLE_MPI
1013
1014 {
1015 // Lx[1] empty => x is a local minimum
1016 if(!Lx[1].empty()) {
1017 // get delta: 1-cell (edge) with minimal G value (steeper gradient)
1018 size_t minId = 0;
1019 for(size_t i = 1; i < Lx[1].size(); ++i) {
1020 const auto &a = Lx[1][minId].lowVerts_[0];
1021 const auto &b = Lx[1][i].lowVerts_[0];
1022 if(a > b) {
1023 // edge[i] < edge[0]
1024 minId = i;
1025 }
1026 }
1027
1028 auto &c_delta = Lx[1][minId];
1029
1030 // store x (0-cell) -> delta (1-cell) V-path
1031 pairCells(Lx[0][0], c_delta, triangulation);
1032
1033 // push every 1-cell in Lx that is not delta into pqZero
1034 for(auto &alpha : Lx[1]) {
1035 if(alpha.id_ != c_delta.id_) {
1036 pqZero.push(alpha);
1037 }
1038 }
1039
1040 // push into pqOne every coface of delta in Lx (2-cells only,
1041 // 3-cells have not any facet paired yet) such that
1042 // numUnpairedFaces == 1
1043 insertCofacets(c_delta, Lx);
1044
1045 while(!pqOne.empty() || !pqZero.empty()) {
1046 while(!pqOne.empty()) {
1047 auto &c_alpha = pqOne.top().get();
1048 pqOne.pop();
1049 auto unpairedFaces = numUnpairedFaces(c_alpha, Lx);
1050 if(unpairedFaces.first == 0) {
1051 pqZero.push(c_alpha);
1052 } else {
1053 auto &c_pair_alpha = Lx[c_alpha.dim_ - 1][unpairedFaces.second];
1054
1055 // store (pair_alpha) -> (alpha) V-path
1056 pairCells(c_pair_alpha, c_alpha, triangulation);
1057
1058 // add cofaces of c_alpha and c_pair_alpha to pqOne
1059 insertCofacets(c_alpha, Lx);
1060 insertCofacets(c_pair_alpha, Lx);
1061 }
1062 }
1063
1064 // skip pair_alpha from pqZero:
1065 // cells in pqZero are not critical if already paired
1066 while(!pqZero.empty() && pqZero.top().get().paired_) {
1067 pqZero.pop();
1068 }
1069
1070 if(!pqZero.empty()) {
1071 auto &c_gamma = pqZero.top().get();
1072 pqZero.pop();
1073
1074 // gamma is a critical cell
1075 // mark gamma as paired
1076 c_gamma.paired_ = true;
1077
1078 // add cofacets of c_gamma to pqOne
1079 insertCofacets(c_gamma, Lx);
1080 }
1081 }
1082 }
1083 }
1084 }
1085
1086 return 0;
1087}
1088
1089template <typename dataType, typename triangulationType>
1090int DiscreteGradient::processLowerStarsStochastic(
1091 const SimplexId *const offsets, const triangulationType &triangulation) {
1092
1093 // WARNING
1094 // If you modify this function, please make sure to also report your edit to
1095 // the other implementation of this function, `processLowerStarsWithMask`.
1096
1097 /* Compute gradient */
1098
1099 auto nverts = triangulation.getNumberOfVertices();
1100
1101 // Comparison function for Cells inside priority queues
1102 const auto orderCells = [&](const CellExt &a, const CellExt &b) -> bool {
1103 return a.lowVerts_ > b.lowVerts_;
1104 };
1105
1106 // Type alias for priority queues
1107 using pqType
1108 = std::priority_queue<std::reference_wrapper<CellExt>,
1109 std::vector<std::reference_wrapper<CellExt>>,
1110 decltype(orderCells)>;
1111
1112 // To reduce allocations, priority queues and lowerStar objects are
1113 // cleaned & reused between iterations.
1114
1115 // Priority queues are pushed at the beginning and popped at the
1116 // end. To pop the minimum, elements should be sorted in a
1117 // decreasing order.
1118 pqType pqZero{orderCells}, pqOne{orderCells};
1119
1120 // store lower star structure
1121 lowerStarType Lx;
1122
1123 // Compute edge length in directions dx, dy and dz
1124 float threshold = 10e-12;
1125 const auto nedges = triangulation.getVertexEdgeNumber(0);
1126 float x1, x2, x3;
1127 float y1, y2, y3;
1128 triangulation.getVertexPoint(0, x1, x2, x3);
1129 std::vector<float> edgeLengths(3);
1130 for(SimplexId i = 0; i < nedges; i++) {
1131 SimplexId edgeId;
1132 triangulation.getVertexEdge(0, i, edgeId);
1133 SimplexId vertexId;
1134 triangulation.getEdgeVertex(edgeId, 0, vertexId);
1135 if(vertexId == 0)
1136 triangulation.getEdgeVertex(edgeId, 1, vertexId);
1137 triangulation.getVertexPoint(vertexId, y1, y2, y3);
1138 if(std::abs(y1 - x1) < threshold && std::abs(y2 - x2) < threshold) {
1139 edgeLengths[2] = std::abs(x3 - y3);
1140 }
1141 if(std::abs(y2 - x2) < threshold && std::abs(y3 - x3) < threshold) {
1142 edgeLengths[0] = std::abs(x1 - y1);
1143 }
1144 if(std::abs(y1 - x1) < threshold && std::abs(y3 - x3) < threshold) {
1145 edgeLengths[1] = std::abs(x2 - y2);
1146 }
1147 }
1148
1149#ifdef TTK_ENABLE_OPENMP
1150#pragma omp parallel for num_threads(threadNumber_) \
1151 firstprivate(Lx, pqZero, pqOne)
1152#endif // TTK_ENABLE_OPENMP
1153 for(SimplexId x = 0; x < nverts; x++) {
1154
1155 // clear priority queues (they should be empty at the end of the
1156 // previous iteration)
1157 while(!pqZero.empty()) {
1158 pqZero.pop();
1159 }
1160 while(!pqOne.empty()) {
1161 pqOne.pop();
1162 }
1163
1164 // Insert into pqOne cofacets of cell c_alpha such as numUnpairedFaces == 1
1165 const auto insertCofacets = [&](const CellExt &ca, lowerStarType &ls) {
1166 if(ca.dim_ == 1) {
1167 for(auto &beta : ls[2]) {
1168 if(ls[1][beta.faces_[0]].id_ == ca.id_
1169 || ls[1][beta.faces_[1]].id_ == ca.id_) {
1170 // edge ca belongs to triangle beta
1171 if(numUnpairedFacesTriangle(beta, ls).first == 1) {
1172 pqOne.push(beta);
1173 }
1174 }
1175 }
1176
1177 } else if(ca.dim_ == 2) {
1178 for(auto &beta : ls[3]) {
1179 if(ls[2][beta.faces_[0]].id_ == ca.id_
1180 || ls[2][beta.faces_[1]].id_ == ca.id_
1181 || ls[2][beta.faces_[2]].id_ == ca.id_) {
1182 // triangle ca belongs to tetra beta
1183 if(numUnpairedFacesTetra(beta, ls).first == 1) {
1184 pqOne.push(beta);
1185 }
1186 }
1187 }
1188 }
1189 };
1190
1191 lowerStar(Lx, x, offsets, triangulation);
1192 // In case the vertex is a ghost, the gradient of the
1193 // simplices of its star is set to GHOST_GRADIENT
1194#ifdef TTK_ENABLE_MPI
1195 if(ttk::isRunningWithMPI()
1196 && triangulation.getVertexRank(x) != ttk::MPIrank_) {
1197 int sizeDim = Lx.size();
1198 for(int i = 0; i < sizeDim; i++) {
1199 int nCells = Lx[i].size();
1200 for(int j = 0; j < nCells; j++) {
1201 setCellToGhost(Lx[i][j].dim_, Lx[i][j].id_);
1202 }
1203 }
1204 } else
1205#endif // TTK_ENABLE_MPI
1206
1207 {
1208 // Lx[1] empty => x is a local minimum
1209 if(!Lx[1].empty()) {
1210
1211 size_t minId = 0;
1212 std::array<float, 3> xCoords;
1213 triangulation.getVertexPoint(x, xCoords[0], xCoords[1], xCoords[2]);
1214 // build stencil
1215 std::vector<SimplexId>
1216 stencilIds; // in order +dx, -dx, +dy, -dy, +dz, -dz
1217 std::vector<float> stencilLengths = edgeLengths;
1218
1219 buildStencil(x, xCoords, triangulation, stencilIds, stencilLengths);
1220
1221 using gradType
1222 = std::conditional_t<std::is_same_v<dataType, double>, double, float>;
1223
1224 const dataType *scalars
1225 = static_cast<dataType const *>(this->inputScalarField_.first);
1226 gradType grad[3];
1227 grad[0] = (scalars[stencilIds[0]] - scalars[stencilIds[1]])
1228 / stencilLengths[0];
1229 grad[1] = (scalars[stencilIds[2]] - scalars[stencilIds[3]])
1230 / stencilLengths[1];
1231 grad[2] = triangulation.getDimensionality() == 3
1232 ? (scalars[stencilIds[4]] - scalars[stencilIds[5]])
1233 / stencilLengths[2]
1234 : 0;
1235
1236 std::vector<double> repartitionBounds;
1237 repartitionBounds.push_back(0);
1238 std::vector<int> indexInLowerStar;
1239 indexInLowerStar.push_back(0);
1240 double totalWeight = 0;
1241 for(size_t i = 0; i < Lx[1].size(); ++i) {
1242 SimplexId vertexId;
1243 triangulation.getEdgeVertex(Lx[1][i].id_, 0, vertexId);
1244 if(vertexId == x)
1245 triangulation.getEdgeVertex(Lx[1][i].id_, 1, vertexId);
1246 if(std::find(stencilIds.begin(), stencilIds.end(), vertexId)
1247 == stencilIds.end())
1248 continue;
1249 std::array<float, 3> newCoords;
1250 triangulation.getVertexPoint(
1251 vertexId, newCoords[0], newCoords[1], newCoords[2]);
1252
1253 float scalarProduct = -(newCoords[0] - xCoords[0]) * grad[0]
1254 - (newCoords[1] - xCoords[1]) * grad[1]
1255 - (newCoords[2] - xCoords[2]) * grad[2];
1256 if(scalarProduct > 0) {
1257 double previousBound
1258 = repartitionBounds[repartitionBounds.size() - 1];
1259 repartitionBounds.push_back(scalarProduct + previousBound);
1260 totalWeight += scalarProduct;
1261 indexInLowerStar.push_back(i);
1262 }
1263 }
1264
1265 for(size_t i = 1; i < repartitionBounds.size(); ++i) {
1266 repartitionBounds[i] /= totalWeight;
1267 }
1268
1269 std::mt19937 gen;
1270 gen.seed(Seed + x);
1271 std::uniform_real_distribution<float> dis(0.0, 1.0);
1272 float random_number = dis(gen);
1273 size_t it = 0;
1274
1275 while(repartitionBounds[it] < random_number
1276 && repartitionBounds.size() > 1) {
1277 it++;
1278 }
1279
1280 minId = indexInLowerStar[it];
1281 auto &c_delta = Lx[1][minId];
1282
1283 // store x (0-cell) -> delta (1-cell) V-path
1284 pairCells(Lx[0][0], c_delta, triangulation);
1285
1286 // push every 1-cell in Lx that is not delta into pqZero
1287 for(auto &alpha : Lx[1]) {
1288 if(alpha.id_ != c_delta.id_) {
1289 pqZero.push(alpha);
1290 }
1291 }
1292
1293 // push into pqOne every coface of delta in Lx (2-cells only,
1294 // 3-cells have not any facet paired yet) such that
1295 // numUnpairedFaces == 1
1296 insertCofacets(c_delta, Lx);
1297
1298 while(!pqOne.empty() || !pqZero.empty()) {
1299 while(!pqOne.empty()) {
1300 auto &c_alpha = pqOne.top().get();
1301 pqOne.pop();
1302 auto unpairedFaces = numUnpairedFaces(c_alpha, Lx);
1303 if(unpairedFaces.first == 0) {
1304 pqZero.push(c_alpha);
1305 } else {
1306 auto &c_pair_alpha = Lx[c_alpha.dim_ - 1][unpairedFaces.second];
1307
1308 // store (pair_alpha) -> (alpha) V-path
1309 pairCells(c_pair_alpha, c_alpha, triangulation);
1310
1311 // add cofaces of c_alpha and c_pair_alpha to pqOne
1312 insertCofacets(c_alpha, Lx);
1313 insertCofacets(c_pair_alpha, Lx);
1314 }
1315 }
1316
1317 // skip pair_alpha from pqZero:
1318 // cells in pqZero are not critical if already paired
1319 while(!pqZero.empty() && pqZero.top().get().paired_) {
1320 pqZero.pop();
1321 }
1322
1323 if(!pqZero.empty()) {
1324 auto &c_gamma = pqZero.top().get();
1325 pqZero.pop();
1326
1327 // gamma is a critical cell
1328 // mark gamma as paired
1329 c_gamma.paired_ = true;
1330
1331 // add cofacets of c_gamma to pqOne
1332 insertCofacets(c_gamma, Lx);
1333 }
1334 }
1335 }
1336 }
1337 }
1338 return 0;
1339}
1340
1341template <typename triangulationType>
1342void DiscreteGradient::buildStencil(const SimplexId &x,
1343 const std::array<float, 3> xCoords,
1344 const triangulationType &triangulation,
1345 std::vector<SimplexId> &stencilIds,
1346 std::vector<float> &stencilLength) {
1347 float threshold = 10e-12;
1348 const auto nedges = triangulation.getVertexEdgeNumber(x);
1349 int dimension = triangulation.getDimensionality();
1350 stencilIds.resize(dimension * 2, -1);
1351 for(SimplexId i = 0; i < nedges; i++) {
1352 SimplexId edgeId;
1353 triangulation.getVertexEdge(x, i, edgeId);
1354 SimplexId vertexId;
1355 triangulation.getEdgeVertex(edgeId, 0, vertexId);
1356 if(vertexId == x) {
1357 triangulation.getEdgeVertex(edgeId, 1, vertexId);
1358 }
1359 std::array<float, 3> currentCoords;
1360 triangulation.getVertexPoint(
1361 vertexId, currentCoords[0], currentCoords[1], currentCoords[2]);
1362 if(std::abs(currentCoords[0] - xCoords[0]) < threshold
1363 && std::abs(currentCoords[1] - xCoords[1]) < threshold
1364 && dimension == 3) {
1365 if(currentCoords[2] - xCoords[2] > threshold) {
1366 stencilIds[4] = vertexId;
1367 } else if(currentCoords[2] - xCoords[2] < -threshold) {
1368 stencilIds[5] = vertexId;
1369 }
1370 } else if(std::abs(currentCoords[0] - xCoords[0]) < threshold
1371 && std::abs(currentCoords[2] - xCoords[2]) < threshold) {
1372 if(currentCoords[1] - xCoords[1] > threshold) {
1373 stencilIds[2] = vertexId;
1374 } else if(currentCoords[1] - xCoords[1] < -threshold) {
1375 stencilIds[3] = vertexId;
1376 }
1377 } else if(std::abs(currentCoords[1] - xCoords[1]) < threshold
1378 && std::abs(currentCoords[2] - xCoords[2]) < threshold) {
1379 if(currentCoords[0] - xCoords[0] > threshold) {
1380 stencilIds[0] = vertexId;
1381 } else if(currentCoords[0] - xCoords[0] < -threshold) {
1382 stencilIds[1] = vertexId;
1383 }
1384 }
1385 }
1386 for(int i = 0; i < dimension; i++) {
1387 if(stencilIds[2 * i] != -1 && stencilIds[2 * i + 1] != -1) {
1388 stencilLength[i] *= 2;
1389 }
1390 }
1391 for(size_t i = 0; i < stencilIds.size(); i++) {
1392 if(stencilIds[i] == -1)
1393 stencilIds[i] = x;
1394 }
1395}
1396
1397template <typename triangulationType>
1398bool DiscreteGradient::isBoundary(
1399 const Cell &cell, const triangulationType &triangulation) const {
1400
1401 if(cell.dim_ > this->dimensionality_ || cell.dim_ < 0) {
1402 return false;
1403 }
1404
1405 const auto vert{this->getCellGreaterVertex(cell, triangulation)};
1406 return triangulation.isVertexOnBoundary(vert);
1407}
1408
1409template <typename triangulationType>
1410SimplexId
1411 DiscreteGradient::getPairedCell(const Cell &cell,
1412 const triangulationType &triangulation,
1413 bool isReverse) const {
1414
1415 // ensure that getPairedCell(Cell, boolean) calls are rejected
1416 static_assert(
1417 std::is_base_of<AbstractTriangulation, triangulationType>(),
1418 "triangulationType should be an AbstractTriangulation derivative");
1419
1420#ifndef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
1421 TTK_FORCE_USE(triangulation);
1422#endif // !TTK_ENABLE_DCG_OPTIMIZE_MEMORY
1423
1424 if((cell.dim_ > this->dimensionality_ - 1 && !isReverse)
1425 || (cell.dim_ > this->dimensionality_ && isReverse) || cell.dim_ < 0) {
1426 return -1;
1427 }
1428
1429 SimplexId id{-1};
1430
1431 if(cell.dim_ == 0) {
1432 if(!isReverse) {
1433#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
1434 const auto locId{(*gradient_)[0][cell.id_]};
1435 if(locId != -1) {
1436 triangulation.getVertexEdge(cell.id_, locId, id);
1437 }
1438#else
1439 id = (*gradient_)[0][cell.id_];
1440#endif
1441 }
1442 }
1443
1444 else if(cell.dim_ == 1) {
1445 if(isReverse) {
1446#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
1447 const auto locId{(*gradient_)[1][cell.id_]};
1448 if(locId != -1) {
1449 triangulation.getEdgeVertex(cell.id_, locId, id);
1450 }
1451#else
1452 id = (*gradient_)[1][cell.id_];
1453#endif
1454 } else {
1455#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
1456 const auto locId{(*gradient_)[2][cell.id_]};
1457 if(locId != -1) {
1458 triangulation.getEdgeTriangle(cell.id_, locId, id);
1459 }
1460#else
1461 id = (*gradient_)[2][cell.id_];
1462#endif
1463 }
1464 }
1465
1466 else if(cell.dim_ == 2) {
1467 if(isReverse) {
1468#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
1469 const auto locId{(*gradient_)[3][cell.id_]};
1470 if(locId != -1) {
1471 triangulation.getTriangleEdge(cell.id_, locId, id);
1472 }
1473#else
1474 id = (*gradient_)[3][cell.id_];
1475#endif
1476 } else {
1477#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
1478 const auto locId{(*gradient_)[4][cell.id_]};
1479 if(locId != -1) {
1480 triangulation.getTriangleStar(cell.id_, locId, id);
1481 }
1482#else
1483 id = (*gradient_)[4][cell.id_];
1484#endif
1485 }
1486 }
1487
1488 else if(cell.dim_ == 3) {
1489 if(isReverse) {
1490#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
1491 const auto locId{(*gradient_)[5][cell.id_]};
1492 if(locId != -1) {
1493 triangulation.getCellTriangle(cell.id_, locId, id);
1494 }
1495#else
1496 id = (*gradient_)[5][cell.id_];
1497#endif
1498 }
1499 }
1500
1501 return id;
1502}
1503
1504template <typename triangulationType>
1505int DiscreteGradient::getDescendingPath(
1506 const Cell &cell,
1507 std::vector<Cell> &vpath,
1508 const triangulationType &triangulation) const {
1509
1510 if(cell.dim_ == 0) {
1511 // assume that cellId is a vertex
1512 SimplexId currentId = cell.id_;
1513 SimplexId connectedEdgeId;
1514 do {
1515 // add a vertex
1516 const Cell vertex(0, currentId);
1517 vpath.push_back(vertex);
1518
1519 if(isCellCritical(vertex)
1520#ifdef TTK_ENABLE_MPI
1521 || triangulation.getVertexRank(currentId) != ttk::MPIrank_
1522#endif
1523 ) {
1524 break;
1525 }
1526
1527 connectedEdgeId = getPairedCell(vertex, triangulation);
1528 if(connectedEdgeId == -1) {
1529 break;
1530 }
1531
1532 // add an edge
1533 const Cell edge(1, connectedEdgeId);
1534 vpath.push_back(edge);
1535
1536 if(isCellCritical(edge)) {
1537 break;
1538 }
1539
1540 for(int i = 0; i < 2; ++i) {
1541 SimplexId vertexId;
1542 triangulation.getEdgeVertex(connectedEdgeId, i, vertexId);
1543
1544 if(vertexId != currentId) {
1545 currentId = vertexId;
1546 break;
1547 }
1548 }
1549
1550 } while(connectedEdgeId != -1);
1551 }
1552
1553 return 0;
1554}
1555
1556template <typename triangulationType>
1557bool DiscreteGradient::getDescendingPathThroughWall(
1558 const Cell &saddle2,
1559 const Cell &saddle1,
1560 const std::vector<bool> &isVisited,
1561 std::vector<Cell> *const vpath,
1562 const triangulationType &triangulation,
1563 const bool stopIfMultiConnected,
1564 const bool enableCycleDetector) const {
1565
1566 // debug
1567 const SimplexId numberOfEdges = triangulation.getNumberOfEdges();
1568 std::vector<bool> isCycle;
1569 if(enableCycleDetector) {
1570 isCycle.resize(numberOfEdges, false);
1571 }
1572
1573 if(dimensionality_ == 3) {
1574 // add the 2-saddle to the path
1575 if(vpath != nullptr) {
1576 vpath->push_back(saddle2);
1577 }
1578
1579 SimplexId currentId = -1;
1580 {
1581 int nconnections = 0;
1582 for(int i = 0; i < 3; ++i) {
1583 SimplexId edgeId;
1584 triangulation.getTriangleEdge(saddle2.id_, i, edgeId);
1585 if(isVisited[edgeId]) {
1586 // saddle2 can be adjacent to saddle1 on the wall
1587 if(isSaddle1(Cell(1, edgeId))) {
1588 if(vpath != nullptr) {
1589 vpath->push_back(Cell(1, edgeId));
1590 }
1591 return false;
1592 }
1593
1594 currentId = edgeId;
1595 ++nconnections;
1596 }
1597 }
1598 if(stopIfMultiConnected && nconnections > 1) {
1599 return true;
1600 }
1601 }
1602
1603 int oldId;
1604 do {
1605
1606 // debug
1607 if(enableCycleDetector) {
1608 if(not isCycle[currentId]) {
1609 isCycle[currentId] = true;
1610 } else {
1611 this->printErr("Cycle detected on the wall of 1-saddle "
1612 + std::to_string(saddle1.id_));
1613 break;
1614 }
1615 }
1616
1617 oldId = currentId;
1618
1619 // add an edge
1620 const Cell edge(1, currentId);
1621 if(vpath != nullptr) {
1622 vpath->push_back(edge);
1623 }
1624
1625 if(isCellCritical(edge)) {
1626 break;
1627 }
1628
1629 const SimplexId connectedTriangleId = getPairedCell(edge, triangulation);
1630
1631 // add a triangle
1632 const Cell triangle(2, connectedTriangleId);
1633 if(vpath != nullptr) {
1634 vpath->push_back(triangle);
1635 }
1636
1637 if(isCellCritical(triangle)) {
1638 break;
1639 }
1640
1641 int nconnections = 0;
1642 for(int i = 0; i < 3; ++i) {
1643 SimplexId edgeId;
1644 triangulation.getTriangleEdge(connectedTriangleId, i, edgeId);
1645
1646 if(isVisited[edgeId] and edgeId != oldId) {
1647 currentId = edgeId;
1648 ++nconnections;
1649 }
1650 }
1651 if(stopIfMultiConnected && nconnections > 1) {
1652 return true;
1653 }
1654
1655 // stop at convergence caused by boundary effect
1656 } while(currentId != oldId);
1657 }
1658
1659 return false;
1660}
1661
1662template <typename triangulationType>
1663int DiscreteGradient::getAscendingPath(const Cell &cell,
1664 std::vector<Cell> &vpath,
1665 const triangulationType &triangulation,
1666 const bool enableCycleDetector) const {
1667
1668 const SimplexId numberOfCells = triangulation.getNumberOfCells();
1669 std::vector<bool> isCycle;
1670 if(enableCycleDetector) {
1671 isCycle.resize(numberOfCells, false);
1672 }
1673
1674 if(dimensionality_ == 2) {
1675 if(cell.dim_ == 2) {
1676 // assume that cellId is a triangle
1677 SimplexId currentId = cell.id_;
1678 SimplexId oldId;
1679 do {
1680 oldId = currentId;
1681
1682 // add a triangle
1683 const Cell triangle(2, currentId);
1684 vpath.push_back(triangle);
1685
1686 if(isCellCritical(triangle)
1687#ifdef TTK_ENABLE_MPI
1688 || triangulation.getTriangleRank(currentId) != ttk::MPIrank_
1689#endif
1690 ) {
1691 break;
1692 }
1693
1694 const SimplexId connectedEdgeId
1695 = getPairedCell(triangle, triangulation, true);
1696 if(connectedEdgeId == -1) {
1697 break;
1698 }
1699
1700 // add an edge
1701 const Cell edge(1, connectedEdgeId);
1702 vpath.push_back(edge);
1703
1704 if(isCellCritical(edge)) {
1705 break;
1706 }
1707
1708 const SimplexId starNumber
1709 = triangulation.getEdgeStarNumber(connectedEdgeId);
1710 for(SimplexId i = 0; i < starNumber; ++i) {
1711 SimplexId starId;
1712 triangulation.getEdgeStar(connectedEdgeId, i, starId);
1713
1714 if(starId != currentId) {
1715 currentId = starId;
1716 break;
1717 }
1718 }
1719
1720 // stop at convergence caused by boundary effect
1721 } while(currentId != oldId);
1722 }
1723 } else if(dimensionality_ == 3) {
1724 if(cell.dim_ == 3) {
1725 // assume that cellId is a tetra
1726 SimplexId currentId = cell.id_;
1727 SimplexId oldId;
1728 do {
1729
1730 // debug
1731 if(enableCycleDetector) {
1732 if(not isCycle[currentId]) {
1733 isCycle[currentId] = true;
1734 } else {
1735 this->printErr("cycle detected in the path from tetra "
1736 + std::to_string(cell.id_));
1737 break;
1738 }
1739 }
1740
1741 oldId = currentId;
1742
1743 // add a tetra
1744 const Cell tetra(3, currentId);
1745 vpath.push_back(tetra);
1746
1747 if(isCellCritical(tetra)
1748#ifdef TTK_ENABLE_MPI
1749 || triangulation.getCellRank(currentId) != ttk::MPIrank_
1750#endif
1751 ) {
1752 break;
1753 }
1754
1755 const SimplexId connectedTriangleId
1756 = getPairedCell(tetra, triangulation, true);
1757 if(connectedTriangleId == -1) {
1758 break;
1759 }
1760
1761 // add a triangle
1762 const Cell triangle(2, connectedTriangleId);
1763 vpath.push_back(triangle);
1764
1765 if(isCellCritical(triangle)) {
1766 break;
1767 }
1768
1769 const SimplexId starNumber
1770 = triangulation.getTriangleStarNumber(connectedTriangleId);
1771 for(SimplexId i = 0; i < starNumber; ++i) {
1772 SimplexId starId;
1773 triangulation.getTriangleStar(connectedTriangleId, i, starId);
1774
1775 if(starId != currentId) {
1776 currentId = starId;
1777 break;
1778 }
1779 }
1780
1781 // stop at convergence caused by boundary effect
1782 } while(currentId != oldId);
1783 }
1784 }
1785
1786 return 0;
1787}
1788
1789template <typename triangulationType>
1790bool DiscreteGradient::getAscendingPathThroughWall(
1791 const Cell &saddle1,
1792 const Cell &saddle2,
1793 const std::vector<bool> &isVisited,
1794 std::vector<Cell> *const vpath,
1795 const triangulationType &triangulation,
1796 const bool stopIfMultiConnected,
1797 const bool enableCycleDetector,
1798 bool *const cycleFound) const {
1799
1800 // debug
1801 const SimplexId numberOfTriangles = triangulation.getNumberOfTriangles();
1802 std::vector<bool> isCycle;
1803 if(enableCycleDetector) {
1804 isCycle.resize(numberOfTriangles, false);
1805 }
1806
1807 if(dimensionality_ == 3) {
1808 // add the 1-saddle to the path
1809 if(vpath != nullptr) {
1810 vpath->push_back(saddle1);
1811 }
1812
1813 SimplexId currentId = -1;
1814 {
1815 int nconnections = 0;
1816 const SimplexId triangleNumber
1817 = triangulation.getEdgeTriangleNumber(saddle1.id_);
1818 for(SimplexId i = 0; i < triangleNumber; ++i) {
1819 SimplexId triangleId;
1820 triangulation.getEdgeTriangle(saddle1.id_, i, triangleId);
1821 if(isVisited[triangleId]) {
1822 // saddle1 can be adjacent to saddle2 on the wall
1823 if(isSaddle2(Cell(2, triangleId))) {
1824 if(vpath != nullptr) {
1825 vpath->push_back(Cell(2, triangleId));
1826 }
1827 return false;
1828 }
1829
1830 currentId = triangleId;
1831 ++nconnections;
1832 }
1833 }
1834 if(stopIfMultiConnected && nconnections > 1) {
1835 return true;
1836 }
1837 }
1838
1839 if(currentId == -1) {
1840 return true;
1841 }
1842
1843 SimplexId oldId;
1844 do {
1845
1846 // debug
1847 if(enableCycleDetector) {
1848 if(not isCycle[currentId]) {
1849 isCycle[currentId] = true;
1850 } else {
1851 if(cycleFound)
1852 *cycleFound = true;
1853 else
1854 this->printErr("Cycle detected on the wall of 2-saddle "
1855 + std::to_string(saddle2.id_));
1856 break;
1857 }
1858 }
1859
1860 oldId = currentId;
1861
1862 // add a triangle
1863 const Cell triangle(2, currentId);
1864 if(vpath != nullptr) {
1865 vpath->push_back(triangle);
1866 }
1867
1868 if(isCellCritical(triangle)) {
1869 break;
1870 }
1871
1872 const SimplexId connectedEdgeId
1873 = getPairedCell(triangle, triangulation, true);
1874
1875 // add an edge
1876 const Cell edge(1, connectedEdgeId);
1877 if(vpath != nullptr) {
1878 vpath->push_back(edge);
1879 }
1880
1881 if(isCellCritical(edge)) {
1882 break;
1883 }
1884
1885 int nconnections = 0;
1886 const SimplexId triangleNumber
1887 = triangulation.getEdgeTriangleNumber(connectedEdgeId);
1888 for(SimplexId i = 0; i < triangleNumber; ++i) {
1889 SimplexId triangleId;
1890 triangulation.getEdgeTriangle(connectedEdgeId, i, triangleId);
1891
1892 if(isVisited[triangleId] and triangleId != oldId) {
1893 currentId = triangleId;
1894 ++nconnections;
1895 }
1896 }
1897 if(stopIfMultiConnected && nconnections > 1) {
1898 return true;
1899 }
1900
1901 // stop at convergence caused by boundary effect
1902 } while(currentId != oldId);
1903 }
1904 return false;
1905}
1906
1907template <typename triangulationType>
1908bool DiscreteGradient::detectGradientCycle(
1909 const Cell &cell, const triangulationType &triangulation) const {
1910 if(dimensionality_ == 3) {
1911 if(cell.dim_ == 1) {
1912 const SimplexId originId = getPairedCell(cell, triangulation);
1913 if(originId == -1)
1914 return false;
1915
1916 std::queue<SimplexId> bfs;
1917 bfs.push(originId);
1918 std::vector<bool> isVisited(triangulation.getNumberOfTriangles(), false);
1919
1920 // BFS traversal
1921 while(!bfs.empty()) {
1922 const SimplexId triangleId = bfs.front();
1923 bfs.pop();
1924
1925 isVisited[triangleId] = true;
1926
1927 for(int j = 0; j < 3; ++j) {
1928 SimplexId edgeId;
1929 triangulation.getTriangleEdge(triangleId, j, edgeId);
1930
1931 const SimplexId pairedCellId
1932 = getPairedCell(Cell(1, edgeId), triangulation);
1933
1934 if(triangleId == pairedCellId or pairedCellId == -1)
1935 continue;
1936
1937 if(isVisited[pairedCellId])
1938 return true;
1939 else
1940 bfs.push(pairedCellId);
1941 }
1942 }
1943 }
1944 }
1945 return false;
1946}
1947
1948template <typename triangulationType>
1949int DiscreteGradient::getDescendingWall(
1950 const Cell &cell,
1951 VisitedMask &mask,
1952 const triangulationType &triangulation,
1953 std::vector<Cell> *const wall,
1954 std::vector<SimplexId> *const saddles) const {
1955
1956 if(saddles != nullptr) {
1957 saddles->clear();
1958 }
1959
1960 if(dimensionality_ == 3) {
1961 if(cell.dim_ == 2) {
1962 // assume that cellId is a triangle
1963 const SimplexId originId = cell.id_;
1964
1965 std::queue<SimplexId> bfs;
1966 bfs.push(originId);
1967
1968 // BFS traversal
1969 while(!bfs.empty()) {
1970 const SimplexId triangleId = bfs.front();
1971 bfs.pop();
1972
1973 if(!mask.isVisited_[triangleId]) {
1974 mask.isVisited_[triangleId] = true;
1975 mask.visitedIds_.emplace_back(triangleId);
1976
1977 // add the triangle
1978 if(wall != nullptr) {
1979 wall->push_back(Cell(2, triangleId));
1980 }
1981
1982 for(int j = 0; j < 3; ++j) {
1983 SimplexId edgeId;
1984 triangulation.getTriangleEdge(triangleId, j, edgeId);
1985
1986 if((saddles != nullptr) and isSaddle1(Cell(1, edgeId))) {
1987 saddles->emplace_back(edgeId);
1988 }
1989
1990 const SimplexId pairedCellId
1991 = getPairedCell(Cell(1, edgeId), triangulation);
1992
1993 if(pairedCellId != -1 and pairedCellId != triangleId) {
1994 bfs.push(pairedCellId);
1995 }
1996 }
1997 }
1998 }
1999
2000 if(saddles != nullptr && saddles->size() > 1) {
2001 std::sort(saddles->begin(), saddles->end());
2002 const auto last = std::unique(saddles->begin(), saddles->end());
2003 saddles->erase(last, saddles->end());
2004 }
2005 }
2006 }
2007
2008 return 0;
2009}
2010
2011template <typename triangulationType>
2012int DiscreteGradient::getAscendingWall(
2013 const Cell &cell,
2014 VisitedMask &mask,
2015 const triangulationType &triangulation,
2016 std::vector<Cell> *const wall,
2017 std::vector<SimplexId> *const saddles) const {
2018
2019 if(saddles != nullptr) {
2020 saddles->clear();
2021 }
2022
2023 if(dimensionality_ == 3) {
2024 if(cell.dim_ == 1) {
2025 // assume that cellId is an edge
2026 const SimplexId originId = cell.id_;
2027
2028 std::queue<SimplexId> bfs;
2029 bfs.push(originId);
2030
2031 // BFS traversal
2032 while(!bfs.empty()) {
2033 const SimplexId edgeId = bfs.front();
2034 bfs.pop();
2035
2036 if(!mask.isVisited_[edgeId]) {
2037 mask.isVisited_[edgeId] = true;
2038 mask.visitedIds_.emplace_back(edgeId);
2039
2040 // add the edge
2041 if(wall != nullptr) {
2042 wall->push_back(Cell(1, edgeId));
2043 }
2044
2045 const SimplexId triangleNumber
2046 = triangulation.getEdgeTriangleNumber(edgeId);
2047 for(SimplexId j = 0; j < triangleNumber; ++j) {
2048 SimplexId triangleId;
2049 triangulation.getEdgeTriangle(edgeId, j, triangleId);
2050
2051 if((saddles != nullptr) and isSaddle2(Cell(2, triangleId))) {
2052 saddles->emplace_back(triangleId);
2053 }
2054
2055 const SimplexId pairedCellId
2056 = getPairedCell(Cell(2, triangleId), triangulation, true);
2057
2058 if(pairedCellId != -1 and pairedCellId != edgeId) {
2059 bfs.push(pairedCellId);
2060 }
2061 }
2062 }
2063 }
2064
2065 if(saddles != nullptr && saddles->size() > 1) {
2066 std::sort(saddles->begin(), saddles->end());
2067 const auto last = std::unique(saddles->begin(), saddles->end());
2068 saddles->erase(last, saddles->end());
2069 }
2070 }
2071 }
2072
2073 return 0;
2074}
2075
2076template <typename triangulationType>
2077int DiscreteGradient::reverseAscendingPath(
2078 const std::vector<Cell> &vpath,
2079 const triangulationType &triangulation) const {
2080
2081 if(dimensionality_ == 2) {
2082 // assume that the first cell is an edge
2083 const SimplexId numberOfCellsInPath = vpath.size();
2084 for(SimplexId i = 0; i < numberOfCellsInPath; i += 2) {
2085 const SimplexId edgeId = vpath[i].id_;
2086 const SimplexId triangleId = vpath[i + 1].id_;
2087
2088#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
2089 for(int k = 0; k < 3; ++k) {
2090 SimplexId tmp;
2091 triangulation.getCellEdge(triangleId, k, tmp);
2092 if(tmp == edgeId) {
2093 (*gradient_)[3][triangleId] = k;
2094 break;
2095 }
2096 }
2097 for(int k = 0; k < triangulation.getEdgeStarNumber(edgeId); ++k) {
2098 SimplexId tmp;
2099 triangulation.getEdgeStar(edgeId, k, tmp);
2100 if(tmp == triangleId) {
2101 (*gradient_)[2][edgeId] = k;
2102 break;
2103 }
2104 }
2105#else
2106 TTK_FORCE_USE(triangulation);
2107 (*gradient_)[3][triangleId] = edgeId;
2108 (*gradient_)[2][edgeId] = triangleId;
2109#endif
2110 }
2111 } else if(dimensionality_ == 3) {
2112 // assume that the first cell is a triangle
2113 const SimplexId numberOfCellsInPath = vpath.size();
2114 for(SimplexId i = 0; i < numberOfCellsInPath; i += 2) {
2115 const SimplexId triangleId = vpath[i].id_;
2116 const SimplexId tetraId = vpath[i + 1].id_;
2117
2118#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
2119 for(int k = 0; k < 4; ++k) {
2120 SimplexId tmp;
2121 triangulation.getCellTriangle(tetraId, k, tmp);
2122 if(tmp == triangleId) {
2123 (*gradient_)[5][tetraId] = k;
2124 break;
2125 }
2126 }
2127 for(int k = 0; k < triangulation.getTriangleStarNumber(triangleId); ++k) {
2128 SimplexId tmp;
2129 triangulation.getTriangleStar(triangleId, k, tmp);
2130 if(tmp == tetraId) {
2131 (*gradient_)[4][triangleId] = k;
2132 break;
2133 }
2134 }
2135#else
2136 (*gradient_)[5][tetraId] = triangleId;
2137 (*gradient_)[4][triangleId] = tetraId;
2138#endif
2139 }
2140 }
2141
2142 return 0;
2143}
2144
2145template <typename triangulationType>
2146int DiscreteGradient::reverseDescendingPath(
2147 const std::vector<Cell> &vpath,
2148 const triangulationType &triangulation) const {
2149
2150 // assume that the first cell is an edge
2151 for(size_t i = 0; i < vpath.size(); i += 2) {
2152 const SimplexId edgeId = vpath[i].id_;
2153 const SimplexId vertId = vpath[i + 1].id_;
2154
2155#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
2156 const auto nneighs = triangulation.getVertexEdgeNumber();
2157 for(int k = 0; k < nneighs; ++k) {
2158 SimplexId tmp;
2159 triangulation.getVertexEdge(vertId, k, tmp);
2160 if(tmp == edgeId) {
2161 (*gradient_)[0][vertId] = k;
2162 break;
2163 }
2164 }
2165 const auto nverts = triangulation.getEdgeStarNumber(edgeId);
2166 for(int k = 0; k < nverts; ++k) {
2167 SimplexId tmp;
2168 triangulation.getEdgeVertex(edgeId, k, tmp);
2169 if(tmp == vertId) {
2170 (*gradient_)[1][edgeId] = k;
2171 break;
2172 }
2173 }
2174#else
2175 TTK_FORCE_USE(triangulation);
2176 (*gradient_)[0][vertId] = edgeId;
2177 (*gradient_)[1][edgeId] = vertId;
2178#endif
2179 }
2180
2181 return 0;
2182}
2183
2184template <typename triangulationType>
2185int DiscreteGradient::reverseAscendingPathOnWall(
2186 const std::vector<Cell> &vpath,
2187 const triangulationType &triangulation,
2188 bool cancelReversal) const {
2189
2190 if(dimensionality_ == 3) {
2191 // assume that the first cell is an edge
2192 if(cancelReversal) {
2193 if(vpath.empty())
2194 return 0;
2195 (*gradient_)[2][vpath[0].id_] = NULL_GRADIENT;
2196 // assume that the last cell is a triangle
2197 if(vpath.size() <= 1)
2198 return 0;
2199 (*gradient_)[3][vpath[vpath.size() - 1].id_] = NULL_GRADIENT;
2200 }
2201 const SimplexId numberOfCellsInPath = vpath.size();
2202 const SimplexId startIndex = (cancelReversal ? 2 : 0);
2203 for(SimplexId i = startIndex; i < numberOfCellsInPath; i += 2) {
2204 const SimplexId vpathEdgeIndex = i;
2205 const SimplexId vpathTriangleIndex = (cancelReversal ? i - 1 : i + 1);
2206 const SimplexId edgeId = vpath[vpathEdgeIndex].id_;
2207 const SimplexId triangleId = vpath[vpathTriangleIndex].id_;
2208
2209#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
2210 for(int k = 0; k < 3; ++k) {
2211 SimplexId tmp;
2212 triangulation.getTriangleEdge(triangleId, k, tmp);
2213 if(tmp == edgeId) {
2214 (*gradient_)[3][triangleId] = k;
2215 break;
2216 }
2217 }
2218 for(int k = 0; k < triangulation.getEdgeTriangleNumber(edgeId); ++k) {
2219 SimplexId tmp;
2220 triangulation.getEdgeTriangle(edgeId, k, tmp);
2221 if(tmp == triangleId) {
2222 (*gradient_)[2][edgeId] = k;
2223 break;
2224 }
2225 }
2226#else
2227 TTK_FORCE_USE(triangulation);
2228 (*gradient_)[3][triangleId] = edgeId;
2229 (*gradient_)[2][edgeId] = triangleId;
2230#endif
2231 }
2232 }
2233
2234 return 0;
2235}
2236
2237template <typename triangulationType>
2238int DiscreteGradient::reverseDescendingPathOnWall(
2239 const std::vector<Cell> &vpath,
2240 const triangulationType &triangulation) const {
2241
2242 if(dimensionality_ == 3) {
2243 // assume that the first cell is a triangle
2244 const SimplexId numberOfCellsInPath = vpath.size();
2245 for(SimplexId i = 0; i < numberOfCellsInPath; i += 2) {
2246 const SimplexId triangleId = vpath[i].id_;
2247 const SimplexId edgeId = vpath[i + 1].id_;
2248
2249#ifdef TTK_ENABLE_DCG_OPTIMIZE_MEMORY
2250 for(int k = 0; k < 3; ++k) {
2251 SimplexId tmp;
2252 triangulation.getTriangleEdge(triangleId, k, tmp);
2253 if(tmp == edgeId) {
2254 (*gradient_)[3][triangleId] = k;
2255 break;
2256 }
2257 }
2258 for(int k = 0; k < triangulation.getEdgeTriangleNumber(edgeId); ++k) {
2259 SimplexId tmp;
2260 triangulation.getEdgeTriangle(edgeId, k, tmp);
2261 if(tmp == triangleId) {
2262 (*gradient_)[2][edgeId] = k;
2263 break;
2264 }
2265 }
2266#else
2267 TTK_FORCE_USE(triangulation);
2268 (*gradient_)[2][edgeId] = triangleId;
2269 (*gradient_)[3][triangleId] = edgeId;
2270#endif
2271 }
2272 }
2273
2274 return 0;
2275}
2276
2277template <typename triangulationType>
2278ttk::SimplexId DiscreteGradient::getCellGreaterVertex(
2279 const Cell c, const triangulationType &triangulation) const {
2280
2281 const auto offsets = this->inputOffsets_;
2282
2283 auto cellDim = c.dim_;
2284 auto cellId = c.id_;
2285
2286 SimplexId vertexId = -1;
2287 if(cellDim == 0) {
2288 vertexId = cellId;
2289 }
2290
2291 else if(cellDim == 1) {
2292 SimplexId v0;
2293 SimplexId v1;
2294 triangulation.getEdgeVertex(cellId, 0, v0);
2295 triangulation.getEdgeVertex(cellId, 1, v1);
2296
2297 if(offsets[v0] > offsets[v1]) {
2298 vertexId = v0;
2299 } else {
2300 vertexId = v1;
2301 }
2302 }
2303
2304 else if(cellDim == 2) {
2305 SimplexId v0{}, v1{}, v2{};
2306 triangulation.getTriangleVertex(cellId, 0, v0);
2307 triangulation.getTriangleVertex(cellId, 1, v1);
2308 triangulation.getTriangleVertex(cellId, 2, v2);
2309 if(offsets[v0] > offsets[v1] && offsets[v0] > offsets[v2]) {
2310 vertexId = v0;
2311 } else if(offsets[v1] > offsets[v0] && offsets[v1] > offsets[v2]) {
2312 vertexId = v1;
2313 } else {
2314 vertexId = v2;
2315 }
2316 }
2317
2318 else if(cellDim == 3) {
2319 SimplexId v0{}, v1{}, v2{}, v3{};
2320 triangulation.getCellVertex(cellId, 0, v0);
2321 triangulation.getCellVertex(cellId, 1, v1);
2322 triangulation.getCellVertex(cellId, 2, v2);
2323 triangulation.getCellVertex(cellId, 3, v3);
2324 if(offsets[v0] > offsets[v1] && offsets[v0] > offsets[v2]
2325 && offsets[v0] > offsets[v3]) {
2326 vertexId = v0;
2327 } else if(offsets[v1] > offsets[v0] && offsets[v1] > offsets[v2]
2328 && offsets[v1] > offsets[v3]) {
2329 vertexId = v1;
2330 } else if(offsets[v2] > offsets[v0] && offsets[v2] > offsets[v1]
2331 && offsets[v2] > offsets[v3]) {
2332 vertexId = v2;
2333 } else {
2334 vertexId = v3;
2335 }
2336 }
2337 return vertexId;
2338}
2339
2340template <typename triangulationType>
2341ttk::SimplexId DiscreteGradient::getCellLowerVertex(
2342 const Cell c, const triangulationType &triangulation) const {
2343
2344 const auto offsets = this->inputOffsets_;
2345
2346 auto cellDim = c.dim_;
2347 auto cellId = c.id_;
2348
2349 SimplexId vertexId = -1;
2350 if(cellDim == 0) {
2351 vertexId = cellId;
2352 }
2353
2354 else if(cellDim == 1) {
2355 SimplexId v0;
2356 SimplexId v1;
2357 triangulation.getEdgeVertex(cellId, 0, v0);
2358 triangulation.getEdgeVertex(cellId, 1, v1);
2359
2360 if(offsets[v0] < offsets[v1]) {
2361 vertexId = v0;
2362 } else {
2363 vertexId = v1;
2364 }
2365 }
2366
2367 else if(cellDim == 2) {
2368 SimplexId v0{}, v1{}, v2{};
2369 triangulation.getTriangleVertex(cellId, 0, v0);
2370 triangulation.getTriangleVertex(cellId, 1, v1);
2371 triangulation.getTriangleVertex(cellId, 2, v2);
2372 if(offsets[v0] < offsets[v1] && offsets[v0] < offsets[v2]) {
2373 vertexId = v0;
2374 } else if(offsets[v1] < offsets[v0] && offsets[v1] < offsets[v2]) {
2375 vertexId = v1;
2376 } else {
2377 vertexId = v2;
2378 }
2379 }
2380
2381 else if(cellDim == 3) {
2382 SimplexId v0{}, v1{}, v2{}, v3{};
2383 triangulation.getCellVertex(cellId, 0, v0);
2384 triangulation.getCellVertex(cellId, 1, v1);
2385 triangulation.getCellVertex(cellId, 2, v2);
2386 triangulation.getCellVertex(cellId, 3, v3);
2387 if(offsets[v0] < offsets[v1] && offsets[v0] < offsets[v2]
2388 && offsets[v0] < offsets[v3]) {
2389 vertexId = v0;
2390 } else if(offsets[v1] < offsets[v0] && offsets[v1] < offsets[v2]
2391 && offsets[v1] < offsets[v3]) {
2392 vertexId = v1;
2393 } else if(offsets[v2] < offsets[v0] && offsets[v2] < offsets[v1]
2394 && offsets[v2] < offsets[v3]) {
2395 vertexId = v2;
2396 } else {
2397 vertexId = v3;
2398 }
2399 }
2400 return vertexId;
2401}
2402
2403template <typename triangulationType>
2404int DiscreteGradient::setGradientGlyphs(
2405 std::vector<std::array<float, 3>> &points,
2406 std::vector<char> &points_pairOrigins,
2407 std::vector<char> &cells_pairTypes,
2408 std::vector<SimplexId> &cellIds,
2409 std::vector<char> &cellDimensions,
2410 const triangulationType &triangulation) const {
2411
2412 const auto nDims = this->getNumberOfDimensions();
2413
2414 // number of glyphs per dimension
2415 std::vector<size_t> nGlyphsPerDim(nDims);
2416
2417#ifdef TTK_ENABLE_OPENMP
2418#pragma omp parallel for num_threads(threadNumber_)
2419#endif // TTK_ENABLE_OPENMP
2420 for(int i = 0; i < nDims - 1; ++i) {
2421 const auto nCells = this->getNumberOfCells(i, triangulation);
2422 for(SimplexId j = 0; j < nCells; ++j) {
2423 if(this->getPairedCell(Cell{i, j}, triangulation) > -1) {
2424 nGlyphsPerDim[i]++;
2425 }
2426 }
2427 }
2428
2429 // partial sum of number of gradient glyphs
2430 std::vector<size_t> offsets(nDims + 1);
2431 for(SimplexId i = 0; i < nDims; ++i) {
2432 offsets[i + 1] = offsets[i] + nGlyphsPerDim[i];
2433 }
2434
2435 // total number of glyphs
2436 const auto nGlyphs = offsets.back();
2437
2438 // resize arrays accordingly
2439 points.resize(2 * nGlyphs);
2440 points_pairOrigins.resize(2 * nGlyphs);
2441 cells_pairTypes.resize(nGlyphs);
2442 cellIds.resize(2 * nGlyphs);
2443 cellDimensions.resize(2 * nGlyphs);
2444
2445#ifdef TTK_ENABLE_OPENMP
2446#pragma omp parallel for num_threads(threadNumber_)
2447#endif // TTK_ENABLE_OPENMP
2448 for(int i = 0; i < nDims - 1; ++i) {
2449 const SimplexId nCells = getNumberOfCells(i, triangulation);
2450 size_t nProcessedGlyphs{offsets[i]};
2451 for(SimplexId j = 0; j < nCells; ++j) {
2452 const Cell c{i, j};
2453 const auto pcid = this->getPairedCell(c, triangulation);
2454 if(pcid > -1) {
2455 const Cell pc{i + 1, pcid};
2456 triangulation.getCellIncenter(
2457 c.id_, c.dim_, points[2 * nProcessedGlyphs].data());
2458 triangulation.getCellIncenter(
2459 pc.id_, pc.dim_, points[2 * nProcessedGlyphs + 1].data());
2460 points_pairOrigins[2 * nProcessedGlyphs] = 0;
2461 points_pairOrigins[2 * nProcessedGlyphs + 1] = 1;
2462 cells_pairTypes[nProcessedGlyphs] = i;
2463#ifdef TTK_ENABLE_MPI
2464 ttk::SimplexId globalId{-1};
2465 triangulation.getDistributedGlobalCellId(j, i, globalId);
2466 cellIds[2 * nProcessedGlyphs + 0] = globalId;
2467 triangulation.getDistributedGlobalCellId(pcid, i + 1, globalId);
2468 cellIds[2 * nProcessedGlyphs + 1] = globalId;
2469#else
2470 cellIds[2 * nProcessedGlyphs + 0] = j;
2471 cellIds[2 * nProcessedGlyphs + 1] = pcid;
2472#endif // TTK_ENABLE_MPI
2473 cellDimensions[2 * nProcessedGlyphs + 0] = i;
2474 cellDimensions[2 * nProcessedGlyphs + 1] = i + 1;
2475 nProcessedGlyphs++;
2476 }
2477 }
2478 }
2479
2480 return 0;
2481}
2482
2483template <typename dataType, typename triangulationType>
2484int DiscreteGradient::simplifySaddleSaddleConnections1(
2485 const std::vector<std::pair<SimplexId, char>> &criticalPoints,
2486 const std::vector<char> &isPL,
2487 const int iterationThreshold,
2488 const bool allowBoundary,
2489 const bool allowBruteForce,
2490 const bool returnSaddleConnectors,
2491 const triangulationType &triangulation) {
2492 Timer t;
2493
2494 // Part 0 : get removable cells
2495 std::vector<char> isRemovableSaddle1;
2496 std::vector<SimplexId> pl2dmt_saddle1(numberOfVertices_, -1);
2497 getRemovableSaddles1<dataType>(criticalPoints, allowBoundary,
2498 isRemovableSaddle1, pl2dmt_saddle1,
2499 triangulation);
2500
2501 std::vector<char> isRemovableSaddle2;
2502 std::vector<SimplexId> pl2dmt_saddle2(numberOfVertices_, -1);
2503 getRemovableSaddles2<dataType>(criticalPoints, allowBoundary,
2504 isRemovableSaddle2, pl2dmt_saddle2,
2505 triangulation);
2506
2507 // Part 1 : initialization
2508 std::vector<VPath> vpaths;
2509 std::vector<CriticalPoint> dmt_criticalPoints;
2510 std::vector<SimplexId> saddle1Index;
2511 std::vector<SimplexId> saddle2Index;
2512 initializeSaddleSaddleConnections1<dataType>(
2513 isRemovableSaddle1, isRemovableSaddle2, allowBruteForce, vpaths,
2514 dmt_criticalPoints, saddle1Index, saddle2Index, triangulation);
2515
2516 // Part 2 : push the vpaths and order by persistence
2517 SaddleSaddleVPathComparator<dataType> cmp_f;
2518 std::set<std::tuple<dataType, SimplexId, SimplexId>,
2519 SaddleSaddleVPathComparator<dataType>>
2520 S(cmp_f);
2521 orderSaddleSaddleConnections1<dataType>(vpaths, dmt_criticalPoints, S);
2522
2523 // Part 3 : process the vpaths
2524 processSaddleSaddleConnections1<dataType>(
2525 iterationThreshold, isPL, allowBoundary, allowBruteForce,
2526 returnSaddleConnectors, S, pl2dmt_saddle1, pl2dmt_saddle2,
2527 isRemovableSaddle1, isRemovableSaddle2, vpaths, dmt_criticalPoints,
2528 saddle1Index, saddle2Index, triangulation);
2529
2530 this->printMsg("Saddle-Saddle pairs simplified #1", 1.0, t.getElapsedTime(),
2531 this->threadNumber_);
2532
2533 return 0;
2534}
2535
2536template <typename dataType, typename triangulationType>
2537int DiscreteGradient::simplifySaddleSaddleConnections2(
2538 const std::vector<std::pair<SimplexId, char>> &criticalPoints,
2539 const std::vector<char> &isPL,
2540 const int iterationThreshold,
2541 const bool allowBoundary,
2542 const bool allowBruteForce,
2543 const bool returnSaddleConnectors,
2544 const triangulationType &triangulation) {
2545 Timer t;
2546
2547 // Part 0 : get removable cells
2548 std::vector<char> isRemovableSaddle1;
2549 std::vector<SimplexId> pl2dmt_saddle1(numberOfVertices_, -1);
2550 getRemovableSaddles1<dataType>(criticalPoints, allowBoundary,
2551 isRemovableSaddle1, pl2dmt_saddle1,
2552 triangulation);
2553
2554 std::vector<char> isRemovableSaddle2;
2555 std::vector<SimplexId> pl2dmt_saddle2(numberOfVertices_, -1);
2556 getRemovableSaddles2<dataType>(criticalPoints, allowBoundary,
2557 isRemovableSaddle2, pl2dmt_saddle2,
2558 triangulation);
2559
2560 // Part 1 : initialization
2561 std::vector<VPath> vpaths;
2562 std::vector<CriticalPoint> dmt_criticalPoints;
2563 std::vector<SimplexId> saddle1Index;
2564 std::vector<SimplexId> saddle2Index;
2565 initializeSaddleSaddleConnections2<dataType>(
2566 isRemovableSaddle1, isRemovableSaddle2, allowBruteForce, vpaths,
2567 dmt_criticalPoints, saddle1Index, saddle2Index, triangulation);
2568
2569 // Part 2 : push the vpaths and order by persistence
2570 SaddleSaddleVPathComparator<dataType> cmp_f;
2571 std::set<std::tuple<dataType, SimplexId, SimplexId>,
2572 SaddleSaddleVPathComparator<dataType>>
2573 S(cmp_f);
2574 orderSaddleSaddleConnections2<dataType>(vpaths, dmt_criticalPoints, S);
2575
2576 // Part 3 : process the vpaths
2577 processSaddleSaddleConnections2<dataType>(
2578 iterationThreshold, isPL, allowBoundary, allowBruteForce,
2579 returnSaddleConnectors, S, pl2dmt_saddle1, pl2dmt_saddle2,
2580 isRemovableSaddle1, isRemovableSaddle2, vpaths, dmt_criticalPoints,
2581 saddle1Index, saddle2Index, triangulation);
2582
2583 this->printMsg("Saddle-Saddle pairs simplified #2", 1.0, t.getElapsedTime(),
2584 this->threadNumber_);
2585
2586 return 0;
2587}
2588
2589template <typename dataType, typename triangulationType>
2590int DiscreteGradient::filterSaddleConnectors(
2591 const bool allowBoundary, const triangulationType &triangulation) {
2592
2593 const bool allowBruteForce = false;
2594 const bool returnSaddleConnectors = true;
2595
2596 // get the node type of a contour tree node (for compatibility with
2597 // ScalarFieldCriticalPoints)
2598 auto getNodeType = [&](const ftm::Node *node) {
2599 const int upDegree = node->getNumberOfUpSuperArcs();
2600 const int downDegree = node->getNumberOfDownSuperArcs();
2601 const int degree = upDegree + downDegree;
2602
2603 // saddle point
2604 if(degree > 1) {
2605 if(upDegree == 2 and downDegree == 1) {
2606 return 2;
2607 } else if(upDegree == 1 and downDegree == 2) {
2608 return 1;
2609 }
2610 }
2611 // local extremum
2612 else {
2613 if(upDegree) {
2614 return 0;
2615 } else {
2616 return 3;
2617 }
2618 }
2619
2620 return -1;
2621 };
2622
2623 std::vector<std::pair<SimplexId, char>> cpset;
2624
2625 const auto *const scalars
2626 = static_cast<const dataType *>(inputScalarField_.first);
2627 const auto *const offsets = inputOffsets_;
2628
2629 ftm::FTMTree contourTree_{};
2630 contourTree_.setDebugLevel(debugLevel_);
2631 contourTree_.setVertexScalars(scalars);
2632 contourTree_.setTreeType(ftm::TreeType::Contour);
2633 contourTree_.setVertexSoSoffsets(offsets);
2634 contourTree_.setThreadNumber(threadNumber_);
2635 contourTree_.setSegmentation(false);
2636 contourTree_.build<dataType>(&triangulation);
2637 ftm::FTMTree_MT *tree = contourTree_.getTree(ftm::TreeType::Contour);
2638
2639 const SimplexId numberOfNodes = tree->getNumberOfNodes();
2640 for(SimplexId nodeId = 0; nodeId < numberOfNodes; ++nodeId) {
2641 const ftm::Node *node = tree->getNode(nodeId);
2642 const SimplexId vertexId = node->getVertexId();
2643
2644 cpset.push_back(std::make_pair(vertexId, getNodeType(node)));
2645 }
2646
2647 std::vector<char> isPL;
2648 getCriticalPointMap(cpset, isPL);
2649
2650 simplifySaddleSaddleConnections1<dataType>(
2651 cpset, isPL, IterationThreshold, allowBoundary, allowBruteForce,
2652 returnSaddleConnectors, triangulation);
2653 simplifySaddleSaddleConnections2<dataType>(
2654 cpset, isPL, IterationThreshold, allowBoundary, allowBruteForce,
2655 returnSaddleConnectors, triangulation);
2656
2657 return 0;
2658}
2659
2660template <typename dataType, typename triangulationType>
2661int DiscreteGradient::getRemovableSaddles1(
2662 const std::vector<std::pair<SimplexId, char>> &criticalPoints,
2663 const bool allowBoundary,
2664 std::vector<char> &isRemovableSaddle,
2665 std::vector<SimplexId> &pl2dmt_saddle,
2666 const triangulationType &triangulation) {
2667
2668 const SimplexId numberOfEdges = triangulation.getNumberOfEdges();
2669 isRemovableSaddle.resize(numberOfEdges);
2670
2671 std::vector<int> dmt1Saddle2PL_(numberOfEdges, -1);
2672
2673 // by default : 1-saddle is removable
2674#ifdef TTK_ENABLE_OPENMP
2675#pragma omp parallel for num_threads(threadNumber_)
2676#endif
2677 for(SimplexId i = 0; i < numberOfEdges; ++i) {
2678 const Cell saddleCandidate(1, i);
2679 isRemovableSaddle[i] = isSaddle1(saddleCandidate);
2680 }
2681
2682 // is [edgeId] in star of PL-1saddle?
2683 for(auto &criticalPoint : criticalPoints) {
2684 const SimplexId criticalPointId = criticalPoint.first;
2685 const char criticalPointType = criticalPoint.second;
2686
2687 if(criticalPointType == static_cast<char>(CriticalType::Saddle1)) {
2688 if(!allowBoundary and triangulation.isVertexOnBoundary(criticalPointId)) {
2689 continue;
2690 }
2691
2692 SimplexId numberOfSaddles = 0;
2693 SimplexId saddleId = -1;
2694 const SimplexId edgeNumber
2695 = triangulation.getVertexEdgeNumber(criticalPointId);
2696 for(SimplexId i = 0; i < edgeNumber; ++i) {
2697 SimplexId edgeId;
2698 triangulation.getVertexEdge(criticalPointId, i, edgeId);
2699 const Cell saddleCandidate(1, edgeId);
2700
2701 if(isSaddle1(saddleCandidate) and dmt1Saddle2PL_[edgeId] == -1) {
2702 saddleId = edgeId;
2703 ++numberOfSaddles;
2704 }
2705 }
2706
2707 // only one DMT-1saddle in the star so this one is non-removable
2708 if(numberOfSaddles == 1) {
2709 if(dmt1Saddle2PL_[saddleId] == -1
2710 and pl2dmt_saddle[criticalPointId] == -1) {
2711 dmt1Saddle2PL_[saddleId] = criticalPointId;
2712 pl2dmt_saddle[criticalPointId] = saddleId;
2713 isRemovableSaddle[saddleId] = false;
2714 }
2715 }
2716 }
2717 }
2718
2719 return 0;
2720}
2721
2722template <typename dataType, typename triangulationType>
2723int DiscreteGradient::getRemovableSaddles2(
2724 const std::vector<std::pair<SimplexId, char>> &criticalPoints,
2725 const bool allowBoundary,
2726 std::vector<char> &isRemovableSaddle,
2727 std::vector<SimplexId> &pl2dmt_saddle,
2728 const triangulationType &triangulation) {
2729
2730 const SimplexId numberOfTriangles = triangulation.getNumberOfTriangles();
2731 isRemovableSaddle.resize(numberOfTriangles);
2732
2733 std::vector<int> dmt2Saddle2PL_(numberOfTriangles, -1);
2734
2735 // by default : 2-saddle is removable
2736#ifdef TTK_ENABLE_OPENMP
2737#pragma omp parallel for num_threads(threadNumber_)
2738#endif
2739 for(SimplexId i = 0; i < numberOfTriangles; ++i) {
2740 const Cell saddleCandidate(2, i);
2741 isRemovableSaddle[i] = isSaddle2(saddleCandidate);
2742 }
2743
2744 // is [triangleId] in star of PL-2saddle?
2745 for(auto &criticalPoint : criticalPoints) {
2746 const SimplexId criticalPointId = criticalPoint.first;
2747 const char criticalPointType = criticalPoint.second;
2748
2749 if(criticalPointType == static_cast<char>(CriticalType::Saddle2)) {
2750 if(!allowBoundary and triangulation.isVertexOnBoundary(criticalPointId)) {
2751 continue;
2752 }
2753
2754 SimplexId numberOfSaddles = 0;
2755 SimplexId saddleId = -1;
2756 const SimplexId triangleNumber
2757 = triangulation.getVertexTriangleNumber(criticalPointId);
2758 for(SimplexId i = 0; i < triangleNumber; ++i) {
2759 SimplexId triangleId;
2760 triangulation.getVertexTriangle(criticalPointId, i, triangleId);
2761 const Cell saddleCandidate(2, triangleId);
2762
2763 if(isSaddle2(saddleCandidate) and dmt2Saddle2PL_[triangleId] == -1) {
2764 saddleId = triangleId;
2765 ++numberOfSaddles;
2766 }
2767 }
2768
2769 // only one DMT-2saddle in the star so this one is non-removable
2770 if(numberOfSaddles == 1) {
2771 if(dmt2Saddle2PL_[saddleId] == -1
2772 and pl2dmt_saddle[criticalPointId] == -1) {
2773 dmt2Saddle2PL_[saddleId] = criticalPointId;
2774 pl2dmt_saddle[criticalPointId] = saddleId;
2775 isRemovableSaddle[saddleId] = false;
2776 }
2777 }
2778 }
2779 }
2780
2781 return 0;
2782}
2783
2784template <typename dataType, typename triangulationType>
2785int DiscreteGradient::initializeSaddleSaddleConnections1(
2786 const std::vector<char> &isRemovableSaddle1,
2787 const std::vector<char> &isRemovableSaddle2,
2788 const bool allowBruteForce,
2789 std::vector<VPath> &vpaths,
2790 std::vector<CriticalPoint> &criticalPoints,
2791 std::vector<SimplexId> &saddle1Index,
2792 std::vector<SimplexId> &saddle2Index,
2793 const triangulationType &triangulation) const {
2794 Timer t;
2795
2796 const auto *const scalars
2797 = static_cast<const dataType *>(inputScalarField_.first);
2798
2799 const int maximumDim = dimensionality_;
2800 const int saddle2Dim = maximumDim - 1;
2801 const int saddle1Dim = saddle2Dim - 1;
2802
2803 // Part 1 : build initial structures
2804 // add the 2-saddles to CriticalPointList
2805 const SimplexId numberOfSaddle2Candidates
2806 = getNumberOfCells(saddle2Dim, triangulation);
2807 saddle2Index.resize(numberOfSaddle2Candidates, -1);
2808 for(SimplexId i = 0; i < numberOfSaddle2Candidates; ++i) {
2809 if(allowBruteForce or isRemovableSaddle2[i]) {
2810 const Cell saddle2Candidate(saddle2Dim, i);
2811
2812 if(isSaddle2(saddle2Candidate)) {
2813 const SimplexId index = criticalPoints.size();
2814 saddle2Index[i] = index;
2815 criticalPoints.emplace_back(saddle2Candidate);
2816 }
2817 }
2818 }
2819 const SimplexId numberOf2Saddles = criticalPoints.size();
2820
2821 // add the 1-saddles to CriticalPointList
2822 const SimplexId numberOfSaddle1Candidates
2823 = getNumberOfCells(saddle1Dim, triangulation);
2824 saddle1Index.resize(numberOfSaddle1Candidates, -1);
2825 for(SimplexId i = 0; i < numberOfSaddle1Candidates; ++i) {
2826 if(isRemovableSaddle1[i]) {
2827 const Cell saddle1Candidate(saddle1Dim, i);
2828
2829 const SimplexId index = criticalPoints.size();
2830 saddle1Index[i] = index;
2831 criticalPoints.emplace_back(saddle1Candidate);
2832 }
2833 }
2834
2835 // Part 2 : update the structures
2836 // apriori: by default construction, the vpaths and segments are not valid
2837 std::vector<bool> isVisited(numberOfSaddle2Candidates, false);
2838 std::vector<SimplexId> visitedTriangles{};
2839
2840 for(SimplexId i = 0; i < numberOf2Saddles; ++i) {
2841 const SimplexId destinationIndex = i;
2842 CriticalPoint &destination = criticalPoints[destinationIndex];
2843 const Cell &saddle2 = destination.cell_;
2844 VisitedMask mask{isVisited, visitedTriangles};
2845
2846 std::vector<SimplexId> saddles1;
2847 getDescendingWall(saddle2, mask, triangulation, nullptr, &saddles1);
2848
2849 for(auto &saddle1Id : saddles1) {
2850 if(!isRemovableSaddle1[saddle1Id]) {
2851 continue;
2852 }
2853
2854 const Cell &saddle1 = Cell(1, saddle1Id);
2855
2856 std::vector<Cell> path;
2857 const bool isMultiConnected = getAscendingPathThroughWall(
2858 saddle1, saddle2, isVisited, &path, triangulation, true);
2859
2860 if(!isMultiConnected) {
2861 const SimplexId sourceIndex = saddle1Index[saddle1Id];
2862 CriticalPoint &source = criticalPoints[sourceIndex];
2863
2864 // update source and destination
2865 const SimplexId sourceSlot = source.addSlot();
2866 const SimplexId destinationSlot = destination.addSlot();
2867
2868 // update vpath
2869 const auto persistence
2870 = getPersistence(saddle2, saddle1, scalars, triangulation);
2871
2872 vpaths.push_back(VPath(true, -1, sourceIndex, destinationIndex,
2873 sourceSlot, destinationSlot, persistence));
2874 }
2875 }
2876 }
2877
2878 // Part 3 : initialize the last structures
2879 const SimplexId numberOfCriticalPoints = criticalPoints.size();
2880 for(SimplexId i = 0; i < numberOfCriticalPoints; ++i) {
2881 CriticalPoint &cp = criticalPoints[i];
2882
2883 const SimplexId numberOfSlots = cp.numberOfSlots_;
2884 cp.vpaths_.resize(numberOfSlots);
2885 cp.numberOfSlots_ = 0;
2886 }
2887
2888 const SimplexId numberOfVPaths = vpaths.size();
2889#ifdef TTK_ENABLE_OPENMP
2890#pragma omp parallel for num_threads(threadNumber_)
2891#endif
2892 for(SimplexId i = 0; i < numberOfVPaths; ++i) {
2893 const VPath &vpath = vpaths[i];
2894
2895 if(vpath.isValid_) {
2896 const SimplexId sourceIndex = vpath.source_;
2897 const SimplexId destinationIndex = vpath.destination_;
2898
2899 const SimplexId sourceSlot = vpath.sourceSlot_;
2900 const SimplexId destinationSlot = vpath.destinationSlot_;
2901
2902 CriticalPoint &source = criticalPoints[sourceIndex];
2903 CriticalPoint &destination = criticalPoints[destinationIndex];
2904
2905 source.vpaths_[sourceSlot] = i;
2906 destination.vpaths_[destinationSlot] = i;
2907 }
2908 }
2909
2910 this->printMsg(
2911 " Initialization step #1", 1.0, t.getElapsedTime(), this->threadNumber_);
2912
2913 return 0;
2914}
2915
2916template <typename dataType, typename triangulationType>
2917int DiscreteGradient::initializeSaddleSaddleConnections2(
2918 const std::vector<char> &isRemovableSaddle1,
2919 const std::vector<char> &isRemovableSaddle2,
2920 const bool allowBruteForce,
2921 std::vector<VPath> &vpaths,
2922 std::vector<CriticalPoint> &criticalPoints,
2923 std::vector<SimplexId> &saddle1Index,
2924 std::vector<SimplexId> &saddle2Index,
2925 const triangulationType &triangulation) const {
2926 Timer t;
2927
2928 const auto *const scalars
2929 = static_cast<const dataType *>(inputScalarField_.first);
2930
2931 const int maximumDim = dimensionality_;
2932 const int saddle2Dim = maximumDim - 1;
2933 const int saddle1Dim = saddle2Dim - 1;
2934
2935 // Part 1 : build initial structures
2936 // add the 1-saddles to CriticalPointList
2937 const SimplexId numberOfSaddle1Candidates
2938 = getNumberOfCells(saddle1Dim, triangulation);
2939 saddle1Index.resize(numberOfSaddle1Candidates, -1);
2940 for(SimplexId i = 0; i < numberOfSaddle1Candidates; ++i) {
2941 if(isRemovableSaddle1[i]) {
2942 const Cell saddle1Candidate(saddle1Dim, i);
2943
2944 const SimplexId index = criticalPoints.size();
2945 saddle1Index[i] = index;
2946 criticalPoints.emplace_back(saddle1Candidate);
2947 }
2948 }
2949 const SimplexId numberOf1Saddles = criticalPoints.size();
2950
2951 // add the 2-saddles to CriticalPointList
2952 const SimplexId numberOfSaddle2Candidates
2953 = getNumberOfCells(saddle2Dim, triangulation);
2954 saddle2Index.resize(numberOfSaddle2Candidates, -1);
2955 for(SimplexId i = 0; i < numberOfSaddle2Candidates; ++i) {
2956 if(allowBruteForce or isRemovableSaddle2[i]) {
2957 const Cell saddle2Candidate(saddle2Dim, i);
2958
2959 if(isSaddle2(saddle2Candidate)) {
2960 const SimplexId index = criticalPoints.size();
2961 saddle2Index[i] = index;
2962 criticalPoints.emplace_back(saddle2Candidate);
2963 }
2964 }
2965 }
2966
2967 // Part 2 : update the structures
2968 // apriori: by default construction, the vpaths and segments are not valid
2969 std::vector<bool> isVisited(numberOfSaddle1Candidates, false);
2970 std::vector<SimplexId> visitedEdges{};
2971 for(SimplexId i = 0; i < numberOf1Saddles; ++i) {
2972 const SimplexId sourceIndex = i;
2973 CriticalPoint &source = criticalPoints[sourceIndex];
2974 const Cell &saddle1 = source.cell_;
2975
2976 std::vector<SimplexId> saddles2;
2977 VisitedMask mask{isVisited, visitedEdges};
2978 getAscendingWall(saddle1, mask, triangulation, nullptr, &saddles2);
2979
2980 for(auto &saddle2Id : saddles2) {
2981 if(!isRemovableSaddle2[saddle2Id]) {
2982 continue;
2983 }
2984
2985 const Cell &saddle2 = Cell(2, saddle2Id);
2986
2987 std::vector<Cell> path;
2988 const bool isMultiConnected = getDescendingPathThroughWall(
2989 saddle2, saddle1, isVisited, &path, triangulation, true);
2990
2991 if(!isMultiConnected) {
2992 const SimplexId destinationIndex = saddle2Index[saddle2Id];
2993 CriticalPoint &destination = criticalPoints[destinationIndex];
2994
2995 // update source and destination
2996 const SimplexId sourceSlot = source.addSlot();
2997 const SimplexId destinationSlot = destination.addSlot();
2998
2999 // update vpath
3000 const auto persistence
3001 = getPersistence(saddle2, saddle1, scalars, triangulation);
3002
3003 vpaths.push_back(VPath(true, -1, sourceIndex, destinationIndex,
3004 sourceSlot, destinationSlot, persistence));
3005 }
3006 }
3007 }
3008
3009 // Part 3 : initialize the last structures
3010 const SimplexId numberOfCriticalPoints = criticalPoints.size();
3011 for(SimplexId i = 0; i < numberOfCriticalPoints; ++i) {
3012 CriticalPoint &cp = criticalPoints[i];
3013
3014 const SimplexId numberOfSlots = cp.numberOfSlots_;
3015 cp.vpaths_.resize(numberOfSlots);
3016 cp.numberOfSlots_ = 0;
3017 }
3018
3019 const SimplexId numberOfVPaths = vpaths.size();
3020#ifdef TTK_ENABLE_OPENMP
3021#pragma omp parallel for num_threads(threadNumber_)
3022#endif
3023 for(SimplexId i = 0; i < numberOfVPaths; ++i) {
3024 const VPath &vpath = vpaths[i];
3025
3026 if(vpath.isValid_) {
3027 const SimplexId sourceIndex = vpath.source_;
3028 const SimplexId destinationIndex = vpath.destination_;
3029
3030 const SimplexId sourceSlot = vpath.sourceSlot_;
3031 const SimplexId destinationSlot = vpath.destinationSlot_;
3032
3033 CriticalPoint &source = criticalPoints[sourceIndex];
3034 CriticalPoint &destination = criticalPoints[destinationIndex];
3035
3036 source.vpaths_[sourceSlot] = i;
3037 destination.vpaths_[destinationSlot] = i;
3038 }
3039 }
3040
3041 this->printMsg(
3042 " Initialization step #2", 1.0, t.getElapsedTime(), this->threadNumber_);
3043
3044 return 0;
3045}
3046
3047template <typename dataType>
3048int DiscreteGradient::orderSaddleSaddleConnections1(
3049 const std::vector<VPath> &vpaths,
3050 std::vector<CriticalPoint> &criticalPoints,
3051 std::set<std::tuple<dataType, SimplexId, SimplexId>,
3052 SaddleSaddleVPathComparator<dataType>> &S) {
3053 Timer t;
3054
3055 const SimplexId numberOfVPaths = vpaths.size();
3056 for(SimplexId i = 0; i < numberOfVPaths; ++i) {
3057 const VPath &vpath = vpaths[i];
3058
3059 if(vpath.isValid_) {
3060 const SimplexId saddleId = criticalPoints[vpath.destination_].cell_.id_;
3061 S.insert(std::make_tuple(vpath.persistence_, i, saddleId));
3062 }
3063 }
3064
3065 this->printMsg(
3066 " Ordering of the vpaths #1", 1.0, t.getElapsedTime(), this->threadNumber_);
3067
3068 return 0;
3069}
3070
3071template <typename dataType>
3072int DiscreteGradient::orderSaddleSaddleConnections2(
3073 const std::vector<VPath> &vpaths,
3074 std::vector<CriticalPoint> &criticalPoints,
3075 std::set<std::tuple<dataType, SimplexId, SimplexId>,
3076 SaddleSaddleVPathComparator<dataType>> &S) {
3077 Timer t;
3078
3079 const SimplexId numberOfVPaths = vpaths.size();
3080 for(SimplexId i = 0; i < numberOfVPaths; ++i) {
3081 const VPath &vpath = vpaths[i];
3082
3083 if(vpath.isValid_) {
3084 const SimplexId saddleId = criticalPoints[vpath.source_].cell_.id_;
3085 S.insert(std::make_tuple(vpath.persistence_, i, saddleId));
3086 }
3087 }
3088
3089 this->printMsg(
3090 " Ordering of the vpaths #2", 1.0, t.getElapsedTime(), this->threadNumber_);
3091
3092 return 0;
3093}
3094
3095template <typename dataType, typename triangulationType>
3096int DiscreteGradient::processSaddleSaddleConnections1(
3097 const int iterationThreshold,
3098 const std::vector<char> &isPL,
3099 const bool allowBoundary,
3100 const bool allowBruteForce,
3101 const bool returnSaddleConnectors,
3102 std::set<std::tuple<dataType, SimplexId, SimplexId>,
3103 SaddleSaddleVPathComparator<dataType>> &S,
3104 std::vector<SimplexId> &pl2dmt_saddle1,
3105 std::vector<SimplexId> &pl2dmt_saddle2,
3106 std::vector<char> &isRemovableSaddle1,
3107 std::vector<char> &isRemovableSaddle2,
3108 std::vector<VPath> &vpaths,
3109 std::vector<CriticalPoint> &criticalPoints,
3110 std::vector<SimplexId> &saddle1Index,
3111 std::vector<SimplexId> &saddle2Index,
3112 const triangulationType &triangulation) {
3113 Timer t;
3114
3115 const auto *const scalars
3116 = static_cast<const dataType *>(inputScalarField_.first);
3117
3118 const SimplexId numberOfEdges = triangulation.getNumberOfEdges();
3119 const SimplexId numberOfTriangles = triangulation.getNumberOfTriangles();
3120 const SimplexId optimizedSize = std::max(numberOfEdges, numberOfTriangles);
3121 std::vector<bool> isVisited(optimizedSize, false);
3122 std::vector<SimplexId> visitedCells{};
3123
3124 int numberOfIterations{};
3125 while(!S.empty()) {
3126 if(iterationThreshold >= 0 and numberOfIterations >= iterationThreshold) {
3127 break;
3128 }
3129
3130 auto ptr = S.begin();
3131 const SimplexId vpathId = std::get<1>(*ptr);
3132 S.erase(ptr);
3133 VPath &vpath = vpaths[vpathId];
3134
3135 if(vpath.isValid_) {
3136 if(returnSaddleConnectors) {
3137 const dataType persistence = vpath.persistence_;
3138 if(persistence > SaddleConnectorsPersistenceThreshold) {
3139 break;
3140 }
3141 }
3142
3143 const Cell &minSaddle1 = criticalPoints[vpath.source_].cell_;
3144 const Cell &minSaddle2 = criticalPoints[vpath.destination_].cell_;
3145
3146 std::vector<SimplexId> saddles1;
3147 VisitedMask mask{isVisited, visitedCells};
3148 getDescendingWall(minSaddle2, mask, triangulation, nullptr, &saddles1);
3149
3150 // check if at least one connection exists
3151 auto isFound
3152 = std::find_if(saddles1.begin(), saddles1.end(),
3153 [&](const auto &s) { return s == minSaddle1.id_; });
3154 if(isFound == saddles1.end()) {
3155 ++numberOfIterations;
3156 continue;
3157 }
3158
3159 // check if there is multiple connections
3160 std::vector<Cell> path;
3161 const bool isMultiConnected = getAscendingPathThroughWall(
3162 minSaddle1, minSaddle2, isVisited, &path, triangulation, true);
3163
3164 if(isMultiConnected) {
3165 ++numberOfIterations;
3166 continue;
3167 }
3168
3169 // filter by 1-saddle condition
3170 if(vpath.isValid_) {
3171 const Cell &dmt_saddle1 = criticalPoints[vpath.source_].cell_;
3172 const SimplexId dmt_saddle1Id = dmt_saddle1.id_;
3173
3174 if(isSaddle1(dmt_saddle1)) {
3175 for(int i = 0; i < 2; ++i) {
3176 SimplexId vertexId;
3177 triangulation.getEdgeVertex(dmt_saddle1Id, i, vertexId);
3178
3179 if(isPL[vertexId] != static_cast<char>(CriticalType::Saddle1)) {
3180 continue;
3181 }
3182
3183 if(!allowBoundary and triangulation.isVertexOnBoundary(vertexId)) {
3184 continue;
3185 }
3186
3187 if(pl2dmt_saddle1[vertexId] == -1) {
3188 const SimplexId pl_saddle1Id = vertexId;
3189
3190 SimplexId numberOfRemainingSaddles1 = 0;
3191
3192 SimplexId savedId = -1;
3193 const SimplexId edgeNumber
3194 = triangulation.getVertexEdgeNumber(pl_saddle1Id);
3195 for(SimplexId j = 0; j < edgeNumber; ++j) {
3196 SimplexId edgeId;
3197 triangulation.getVertexEdge(pl_saddle1Id, j, edgeId);
3198
3199 if(edgeId != dmt_saddle1Id and isSaddle1(Cell(1, edgeId))
3200 and isRemovableSaddle1[edgeId]) {
3201 ++numberOfRemainingSaddles1;
3202 savedId = edgeId;
3203 }
3204 }
3205
3206 if(numberOfRemainingSaddles1 == 0) {
3207 isRemovableSaddle1[dmt_saddle1Id] = false;
3208 pl2dmt_saddle1[vertexId] = dmt_saddle1Id;
3209 // dmt1Saddle2PL_[dmt_saddle1Id] = vertexId;
3210 vpath.invalidate();
3211 break;
3212 }
3213 if(numberOfRemainingSaddles1 == 1) {
3214 isRemovableSaddle1[dmt_saddle1Id] = false;
3215 isRemovableSaddle1[savedId] = false;
3216 pl2dmt_saddle1[vertexId] = savedId;
3217 // dmt1Saddle2PL_[savedId] = vertexId;
3218 break;
3219 }
3220 } else if(pl2dmt_saddle1[vertexId] == dmt_saddle1Id) {
3221 vpath.invalidate();
3222 break;
3223 }
3224 }
3225 } else {
3226 vpath.invalidate();
3227 }
3228 }
3229
3230 // filter by 2-saddle condition
3231 if(!allowBruteForce and vpath.isValid_) {
3232 const Cell &dmt_saddle2 = criticalPoints[vpath.destination_].cell_;
3233 const SimplexId dmt_saddle2Id = dmt_saddle2.id_;
3234
3235 if(isSaddle2(dmt_saddle2)) {
3236 for(int i = 0; i < 3; ++i) {
3237 SimplexId vertexId;
3238 triangulation.getTriangleVertex(dmt_saddle2Id, i, vertexId);
3239
3240 if(isPL[vertexId] != static_cast<char>(CriticalType::Saddle2)) {
3241 continue;
3242 }
3243
3244 if(!allowBoundary and triangulation.isVertexOnBoundary(vertexId)) {
3245 continue;
3246 }
3247
3248 if(pl2dmt_saddle2[vertexId] == -1) {
3249 const SimplexId pl_saddle2Id = vertexId;
3250
3251 SimplexId numberOfRemainingSaddles2 = 0;
3252
3253 SimplexId savedId = -1;
3254 const SimplexId triangleNumber
3255 = triangulation.getVertexTriangleNumber(pl_saddle2Id);
3256 for(SimplexId j = 0; j < triangleNumber; ++j) {
3257 SimplexId triangleId;
3258 triangulation.getVertexTriangle(pl_saddle2Id, j, triangleId);
3259
3260 if(triangleId != dmt_saddle2Id
3261 and isSaddle2(Cell(2, triangleId))
3262 and isRemovableSaddle2[triangleId]) {
3263 ++numberOfRemainingSaddles2;
3264 savedId = triangleId;
3265 }
3266 }
3267
3268 if(numberOfRemainingSaddles2 == 0) {
3269 isRemovableSaddle2[dmt_saddle2Id] = false;
3270 pl2dmt_saddle2[vertexId] = dmt_saddle2Id;
3271 // dmt2Saddle2PL_[dmt_saddle2Id] = vertexId;
3272 vpath.invalidate();
3273 break;
3274 }
3275 if(numberOfRemainingSaddles2 == 1) {
3276 isRemovableSaddle2[dmt_saddle2Id] = false;
3277 isRemovableSaddle2[savedId] = false;
3278 pl2dmt_saddle2[vertexId] = savedId;
3279 // dmt2Saddle2PL_[savedId] = vertexId;
3280 break;
3281 }
3282 } else if(pl2dmt_saddle2[vertexId] == dmt_saddle2Id) {
3283 vpath.invalidate();
3284 break;
3285 }
3286 }
3287 } else {
3288 vpath.invalidate();
3289 }
3290 }
3291
3292 if(vpath.isValid_) {
3293 reverseAscendingPathOnWall(path, triangulation);
3294 }
3295 }
3296
3297 if(vpath.isValid_) {
3298
3299 const SimplexId sourceId = vpath.source_;
3300 const SimplexId destinationId = vpath.destination_;
3301
3302 // invalidate vpaths connected to destination
3303 std::vector<SimplexId> newSourceIds;
3304 CriticalPoint &destination = criticalPoints[destinationId];
3305 for(auto &destinationVPathId : destination.vpaths_) {
3306 VPath &destinationVPath = vpaths[destinationVPathId];
3307
3308 if(destinationVPath.isValid_ and destinationVPath.source_ != sourceId) {
3309 // save critical point
3310 const SimplexId newSourceId = destinationVPath.source_;
3311 newSourceIds.push_back(newSourceId);
3312
3313 // clear vpath
3314 destinationVPath.invalidate();
3315 }
3316 }
3317
3318 // invalidate vpaths connected to source and save the critical points to
3319 // update
3320 std::vector<SimplexId> newDestinationIds;
3321 CriticalPoint &source = criticalPoints[sourceId];
3322 for(auto &sourceVPathId : source.vpaths_) {
3323 VPath &sourceVPath = vpaths[sourceVPathId];
3324
3325 if(sourceVPath.isValid_ and sourceVPath.destination_ != destinationId) {
3326 // save critical point
3327 const SimplexId newDestinationId = sourceVPath.destination_;
3328 newDestinationIds.push_back(newDestinationId);
3329
3330 CriticalPoint &newDestination = criticalPoints[newDestinationId];
3331 for(auto &newDestinationVPathId : newDestination.vpaths_) {
3332 VPath &newDestinationVPath = vpaths[newDestinationVPathId];
3333 if(newDestinationVPath.isValid_
3334 and newDestinationVPath.source_ != sourceId) {
3335
3336 // clear vpath
3337 newDestinationVPath.invalidate();
3338 }
3339 }
3340
3341 // clear vpath
3342 sourceVPath.invalidate();
3343 }
3344 }
3345
3346 // finally invalidate current vpath and critical points
3347 vpath.invalidate();
3348 source.clear();
3349 destination.clear();
3350
3351 // look at the gradient : reconnect locally the critical points
3352 for(auto &newDestinationId : newDestinationIds) {
3353 CriticalPoint &newDestination = criticalPoints[newDestinationId];
3354 const Cell &saddle2 = newDestination.cell_;
3355
3356 std::vector<SimplexId> saddles1;
3357 VisitedMask mask{isVisited, visitedCells};
3358 getDescendingWall(saddle2, mask, triangulation, nullptr, &saddles1);
3359
3360 for(auto &saddle1Id : saddles1) {
3361 const Cell saddle1(1, saddle1Id);
3362
3363 std::vector<Cell> path;
3364 const bool isMultiConnected = getAscendingPathThroughWall(
3365 saddle1, saddle2, isVisited, &path, triangulation, true);
3366
3367 if(isMultiConnected) {
3368 continue;
3369 }
3370
3371 SimplexId newSourceId = saddle1Index[saddle1Id];
3372
3373 // connection to a new saddle1 (not present in the graph before)
3374 if(newSourceId == -1) {
3375 if(!isRemovableSaddle1[saddle1Id]) {
3376 continue;
3377 }
3378
3379 const SimplexId newCriticalPointId = criticalPoints.size();
3380 saddle1Index[saddle1Id] = newCriticalPointId;
3381 criticalPoints.emplace_back(saddle1);
3382
3383 newSourceId = newCriticalPointId;
3384 }
3385 CriticalPoint &newSource = criticalPoints[newSourceId];
3386
3387 // update vpaths
3388 const SimplexId newVPathId = vpaths.size();
3389 const auto persistence
3390 = getPersistence(saddle2, saddle1, scalars, triangulation);
3391
3392 vpaths.push_back(VPath(
3393 true, -1, newSourceId, newDestinationId, -1, -1, persistence));
3394
3395 // update criticalPoints
3396 newDestination.vpaths_.push_back(newVPathId);
3397 newSource.vpaths_.push_back(newVPathId);
3398
3399 // update set
3400 S.insert(
3401 std::make_tuple(persistence, newVPathId, newDestination.cell_.id_));
3402 }
3403 }
3404
3405 // look at the gradient : get the links not predicted by the graph
3406 for(auto &newSourceId : newSourceIds) {
3407 CriticalPoint &newSource = criticalPoints[newSourceId];
3408 const Cell &saddle1 = newSource.cell_;
3409
3410 std::vector<SimplexId> saddles2;
3411 VisitedMask mask{isVisited, visitedCells};
3412 getAscendingWall(saddle1, mask, triangulation, nullptr, &saddles2);
3413
3414 for(auto &saddle2Id : saddles2) {
3415 const Cell saddle2(2, saddle2Id);
3416
3417 std::vector<Cell> path;
3418 const bool isMultiConnected = getDescendingPathThroughWall(
3419 saddle2, saddle1, isVisited, &path, triangulation, true);
3420
3421 if(isMultiConnected) {
3422 continue;
3423 }
3424
3425 const SimplexId newDestinationId = saddle2Index[saddle2Id];
3426
3427 // connection to a new saddle2 (not present in the graph before)
3428 if(newDestinationId == -1) {
3429 continue;
3430 }
3431
3432 CriticalPoint &newDestination = criticalPoints[newDestinationId];
3433
3434 // check existence of the possibly newVPath in the graph
3435 bool alreadyExists = false;
3436 for(auto &newDestinationVPathId : newDestination.vpaths_) {
3437 const VPath &newDestinationVPath = vpaths[newDestinationVPathId];
3438
3439 if(newDestinationVPath.isValid_
3440 and newDestinationVPath.source_ == newSourceId) {
3441 alreadyExists = true;
3442 break;
3443 }
3444 }
3445
3446 if(alreadyExists) {
3447 continue;
3448 }
3449
3450 // update vpaths
3451 const SimplexId newVPathId = vpaths.size();
3452 const auto persistence
3453 = getPersistence(saddle2, saddle1, scalars, triangulation);
3454
3455 vpaths.push_back(VPath(
3456 true, -1, newSourceId, newDestinationId, -1, -1, persistence));
3457
3458 // update criticalPoints
3459 newDestination.vpaths_.push_back(newVPathId);
3460 newSource.vpaths_.push_back(newVPathId);
3461
3462 // update set
3463 S.insert(
3464 std::make_tuple(persistence, newVPathId, newDestination.cell_.id_));
3465 }
3466 }
3467 }
3468
3469 ++numberOfIterations;
3470 }
3471
3472 this->printMsg(" Processing of the vpaths #1", 1.0, t.getElapsedTime(),
3473 this->threadNumber_);
3474
3475 return 0;
3476}
3477
3478template <typename dataType, typename triangulationType>
3479int DiscreteGradient::processSaddleSaddleConnections2(
3480 const int iterationThreshold,
3481 const std::vector<char> &isPL,
3482 const bool allowBoundary,
3483 const bool allowBruteForce,
3484 const bool returnSaddleConnectors,
3485 std::set<std::tuple<dataType, SimplexId, SimplexId>,
3486 SaddleSaddleVPathComparator<dataType>> &S,
3487 std::vector<SimplexId> &pl2dmt_saddle1,
3488 std::vector<SimplexId> &pl2dmt_saddle2,
3489 std::vector<char> &isRemovableSaddle1,
3490 std::vector<char> &isRemovableSaddle2,
3491 std::vector<VPath> &vpaths,
3492 std::vector<CriticalPoint> &criticalPoints,
3493 std::vector<SimplexId> &saddle1Index,
3494 std::vector<SimplexId> &saddle2Index,
3495 const triangulationType &triangulation) {
3496 Timer t;
3497
3498 this->printMsg("Saddle connector persistence threshold: "
3499 + std::to_string(this->SaddleConnectorsPersistenceThreshold));
3500
3501 const auto *const scalars
3502 = static_cast<const dataType *>(inputScalarField_.first);
3503
3504 const SimplexId numberOfEdges = triangulation.getNumberOfEdges();
3505 const SimplexId numberOfTriangles = triangulation.getNumberOfTriangles();
3506 const SimplexId optimizedSize = std::max(numberOfEdges, numberOfTriangles);
3507 std::vector<bool> isVisited(optimizedSize, false);
3508 std::vector<SimplexId> visitedIds{};
3509
3510 int numberOfIterations{};
3511 while(!S.empty()) {
3512 if(iterationThreshold >= 0 and numberOfIterations >= iterationThreshold) {
3513 break;
3514 }
3515
3516 auto ptr = S.begin();
3517 const SimplexId vpathId = std::get<1>(*ptr);
3518 S.erase(ptr);
3519 VPath &vpath = vpaths[vpathId];
3520
3521 if(vpath.isValid_) {
3522 if(returnSaddleConnectors) {
3523 const dataType persistence = vpath.persistence_;
3524 if(persistence > SaddleConnectorsPersistenceThreshold) {
3525 break;
3526 }
3527 }
3528
3529 const Cell &minSaddle1 = criticalPoints[vpath.source_].cell_;
3530 const Cell &minSaddle2 = criticalPoints[vpath.destination_].cell_;
3531
3532 std::vector<SimplexId> saddles2;
3533 VisitedMask mask{isVisited, visitedIds};
3534 getAscendingWall(minSaddle1, mask, triangulation, nullptr, &saddles2);
3535
3536 // check if at least one connection exists
3537 auto isFound
3538 = std::find_if(saddles2.begin(), saddles2.end(),
3539 [&](const auto &s) { return s == minSaddle2.id_; });
3540
3541 if(isFound == saddles2.end()) {
3542 ++numberOfIterations;
3543 continue;
3544 }
3545
3546 // check if there is multiple connections
3547 std::vector<Cell> path;
3548 const bool isMultiConnected = getDescendingPathThroughWall(
3549 minSaddle2, minSaddle1, isVisited, &path, triangulation, true);
3550
3551 if(isMultiConnected) {
3552 ++numberOfIterations;
3553 continue;
3554 }
3555
3556 // filter by 1-saddle condition
3557 if(vpath.isValid_) {
3558 const Cell &dmt_saddle1 = criticalPoints[vpath.source_].cell_;
3559 const SimplexId dmt_saddle1Id = dmt_saddle1.id_;
3560
3561 if(isSaddle1(dmt_saddle1)) {
3562 for(int i = 0; i < 2; ++i) {
3563 SimplexId vertexId;
3564 triangulation.getEdgeVertex(dmt_saddle1Id, i, vertexId);
3565
3566 if(isPL[vertexId] != static_cast<char>(CriticalType::Saddle1)) {
3567 continue;
3568 }
3569
3570 if(!allowBoundary and triangulation.isVertexOnBoundary(vertexId)) {
3571 continue;
3572 }
3573
3574 if(pl2dmt_saddle1[vertexId] == -1) {
3575 const SimplexId pl_saddle1Id = vertexId;
3576
3577 SimplexId numberOfRemainingSaddles1 = 0;
3578
3579 SimplexId savedId = -1;
3580 const SimplexId edgeNumber
3581 = triangulation.getVertexEdgeNumber(pl_saddle1Id);
3582 for(SimplexId j = 0; j < edgeNumber; ++j) {
3583 SimplexId edgeId;
3584 triangulation.getVertexEdge(pl_saddle1Id, j, edgeId);
3585
3586 if(edgeId != dmt_saddle1Id and isSaddle1(Cell(1, edgeId))
3587 and isRemovableSaddle1[edgeId]) {
3588 ++numberOfRemainingSaddles1;
3589 savedId = edgeId;
3590 }
3591 }
3592
3593 if(numberOfRemainingSaddles1 == 0) {
3594 isRemovableSaddle1[dmt_saddle1Id] = false;
3595 pl2dmt_saddle1[vertexId] = dmt_saddle1Id;
3596 // dmt1Saddle2PL_[dmt_saddle1Id] = vertexId;
3597 vpath.invalidate();
3598 break;
3599 }
3600 if(numberOfRemainingSaddles1 == 1) {
3601 isRemovableSaddle1[dmt_saddle1Id] = false;
3602 isRemovableSaddle1[savedId] = false;
3603 pl2dmt_saddle1[vertexId] = savedId;
3604 // dmt1Saddle2PL_[savedId] = vertexId;
3605 break;
3606 }
3607 } else if(pl2dmt_saddle1[vertexId] == dmt_saddle1Id) {
3608 vpath.invalidate();
3609 break;
3610 }
3611 }
3612 } else {
3613 vpath.invalidate();
3614 }
3615 }
3616
3617 // filter by 2-saddle condition
3618 if(!allowBruteForce and vpath.isValid_) {
3619 const Cell &dmt_saddle2 = criticalPoints[vpath.destination_].cell_;
3620 const SimplexId dmt_saddle2Id = dmt_saddle2.id_;
3621
3622 if(isSaddle2(dmt_saddle2)) {
3623 for(int i = 0; i < 3; ++i) {
3624 SimplexId vertexId;
3625 triangulation.getTriangleVertex(dmt_saddle2Id, i, vertexId);
3626
3627 if(isPL[vertexId] != static_cast<char>(CriticalType::Saddle2)) {
3628 continue;
3629 }
3630
3631 if(!allowBoundary and triangulation.isVertexOnBoundary(vertexId)) {
3632 continue;
3633 }
3634
3635 if(pl2dmt_saddle2[vertexId] == -1) {
3636 const SimplexId pl_saddle2Id = vertexId;
3637
3638 SimplexId numberOfRemainingSaddles2 = 0;
3639
3640 SimplexId savedId = -1;
3641 const SimplexId triangleNumber
3642 = triangulation.getVertexTriangleNumber(pl_saddle2Id);
3643 for(SimplexId j = 0; j < triangleNumber; ++j) {
3644 SimplexId triangleId;
3645 triangulation.getVertexTriangle(pl_saddle2Id, j, triangleId);
3646
3647 if(triangleId != dmt_saddle2Id
3648 and isSaddle2(Cell(2, triangleId))
3649 and isRemovableSaddle2[triangleId]) {
3650 ++numberOfRemainingSaddles2;
3651 savedId = triangleId;
3652 }
3653 }
3654
3655 if(!numberOfRemainingSaddles2) {
3656 isRemovableSaddle2[dmt_saddle2Id] = false;
3657 pl2dmt_saddle2[vertexId] = dmt_saddle2Id;
3658 vpath.invalidate();
3659 break;
3660 }
3661 if(numberOfRemainingSaddles2 == 1) {
3662 isRemovableSaddle2[dmt_saddle2Id] = false;
3663 isRemovableSaddle2[savedId] = false;
3664 pl2dmt_saddle2[vertexId] = savedId;
3665 break;
3666 }
3667 } else if(pl2dmt_saddle2[vertexId] == dmt_saddle2Id) {
3668 vpath.invalidate();
3669 break;
3670 }
3671 }
3672 } else {
3673 vpath.invalidate();
3674 }
3675 }
3676
3677 if(vpath.isValid_) {
3678 reverseDescendingPathOnWall(path, triangulation);
3679 }
3680 }
3681
3682 if(vpath.isValid_) {
3683 // add persistence pair to collection if necessary
3684 // if(CollectPersistencePairs and outputPersistencePairs_) {
3685 // const Cell &minSaddle1 = criticalPoints[vpath.source_].cell_;
3686 // const Cell &minSaddle2 = criticalPoints[vpath.destination_].cell_;
3687 // outputPersistencePairs_->push_back({minSaddle1, minSaddle2});
3688 //}
3689
3690 const SimplexId sourceId = vpath.source_;
3691 const SimplexId destinationId = vpath.destination_;
3692
3693 // invalidate vpaths connected to source
3694 std::vector<SimplexId> newDestinationIds;
3695 CriticalPoint &source = criticalPoints[sourceId];
3696 for(auto &sourceVPathId : source.vpaths_) {
3697 VPath &sourceVPath = vpaths[sourceVPathId];
3698
3699 if(sourceVPath.isValid_ and sourceVPath.destination_ != destinationId) {
3700 // save critical point
3701 const SimplexId newDestinationId = sourceVPath.destination_;
3702 newDestinationIds.push_back(newDestinationId);
3703
3704 // clear vpath
3705 sourceVPath.invalidate();
3706 }
3707 }
3708
3709 // invalidate vpaths connected to destination and save the critical
3710 // points to update
3711 std::vector<SimplexId> newSourceIds;
3712 CriticalPoint &destination = criticalPoints[destinationId];
3713 for(auto &destinationVPathId : destination.vpaths_) {
3714 VPath &destinationVPath = vpaths[destinationVPathId];
3715
3716 if(destinationVPath.isValid_ and destinationVPath.source_ != sourceId) {
3717 // save critical point
3718 const SimplexId newSourceId = destinationVPath.source_;
3719 newSourceIds.push_back(newSourceId);
3720
3721 CriticalPoint &newSource = criticalPoints[newSourceId];
3722 for(auto &newSourceVPathId : newSource.vpaths_) {
3723 VPath &newSourceVPath = vpaths[newSourceVPathId];
3724 if(newSourceVPath.isValid_
3725 and newSourceVPath.destination_ != destinationId) {
3726
3727 // clear vpath
3728 newSourceVPath.invalidate();
3729 }
3730 }
3731
3732 // clear vpath
3733 destinationVPath.invalidate();
3734 }
3735 }
3736
3737 // finally invalidate current vpath and critical points
3738 vpath.invalidate();
3739 source.clear();
3740 destination.clear();
3741
3742 // look at the gradient : reconnect locally the critical points
3743 for(auto &newSourceId : newSourceIds) {
3744 CriticalPoint &newSource = criticalPoints[newSourceId];
3745 const Cell &saddle1 = newSource.cell_;
3746
3747 std::vector<SimplexId> saddles2;
3748 VisitedMask mask{isVisited, visitedIds};
3749 getAscendingWall(saddle1, mask, triangulation, nullptr, &saddles2);
3750
3751 for(auto &saddle2Id : saddles2) {
3752 const Cell saddle2(2, saddle2Id);
3753
3754 const bool isMultiConnected = getDescendingPathThroughWall(
3755 saddle2, saddle1, isVisited, nullptr, triangulation, true);
3756
3757 if(isMultiConnected) {
3758 continue;
3759 }
3760
3761 SimplexId newDestinationId = saddle2Index[saddle2Id];
3762
3763 // connection to a new saddle2 (not present in the graph before)
3764 if(newDestinationId == -1) {
3765 if(!isRemovableSaddle2[saddle2Id]) {
3766 continue;
3767 }
3768
3769 const SimplexId newCriticalPointId = criticalPoints.size();
3770 saddle2Index[saddle2Id] = newCriticalPointId;
3771 criticalPoints.emplace_back(saddle2);
3772
3773 newDestinationId = newCriticalPointId;
3774 }
3775
3776 CriticalPoint &newDestination = criticalPoints[newDestinationId];
3777
3778 // update vpaths
3779 const SimplexId newVPathId = vpaths.size();
3780 const auto persistence
3781 = getPersistence(saddle2, saddle1, scalars, triangulation);
3782
3783 vpaths.push_back(VPath(
3784 true, -1, newSourceId, newDestinationId, -1, -1, persistence));
3785
3786 // update criticalPoints
3787 newDestination.vpaths_.push_back(newVPathId);
3788 newSource.vpaths_.push_back(newVPathId);
3789
3790 // update set
3791 S.insert(
3792 std::make_tuple(persistence, newVPathId, newSource.cell_.id_));
3793 }
3794 }
3795
3796 // look at the gradient : get the links not predicted by the graph
3797 for(auto &newDestinationId : newDestinationIds) {
3798 CriticalPoint &newDestination = criticalPoints[newDestinationId];
3799 const Cell &saddle2 = newDestination.cell_;
3800
3801 std::vector<SimplexId> saddles1;
3802 VisitedMask mask{isVisited, visitedIds};
3803 getDescendingWall(saddle2, mask, triangulation, nullptr, &saddles1);
3804
3805 for(auto &saddle1Id : saddles1) {
3806 const Cell saddle1(1, saddle1Id);
3807
3808 std::vector<Cell> path;
3809 const bool isMultiConnected = getAscendingPathThroughWall(
3810 saddle1, saddle2, isVisited, &path, triangulation, true);
3811
3812 if(isMultiConnected) {
3813 continue;
3814 }
3815
3816 const SimplexId newSourceId = saddle1Index[saddle1Id];
3817
3818 if(newSourceId == -1) {
3819 continue;
3820 }
3821
3822 CriticalPoint &newSource = criticalPoints[newSourceId];
3823
3824 // check existence of the possibly newVPath in the graph
3825 bool alreadyExists = false;
3826 for(auto &newSourceVPathId : newSource.vpaths_) {
3827 const VPath &newSourceVPath = vpaths[newSourceVPathId];
3828
3829 if(newSourceVPath.isValid_
3830 and newSourceVPath.destination_ == newDestinationId) {
3831 alreadyExists = true;
3832 break;
3833 }
3834 }
3835
3836 if(alreadyExists) {
3837 continue;
3838 }
3839
3840 // update vpaths
3841 const SimplexId newVPathId = vpaths.size();
3842 const auto persistence
3843 = getPersistence(saddle2, saddle1, scalars, triangulation);
3844
3845 vpaths.push_back(VPath(
3846 true, -1, newSourceId, newDestinationId, -1, -1, persistence));
3847
3848 // update criticalPoints
3849 newDestination.vpaths_.push_back(newVPathId);
3850 newSource.vpaths_.push_back(newVPathId);
3851
3852 // update set
3853 S.insert(
3854 std::make_tuple(persistence, newVPathId, newSource.cell_.id_));
3855 }
3856 }
3857 }
3858
3859 ++numberOfIterations;
3860 }
3861
3862 this->printMsg(" Processing of the vpaths #2", 1.0, t.getElapsedTime(),
3863 this->threadNumber_);
3864
3865 return 0;
3866}
TTK_FORCE_USE
#define TTK_FORCE_USE(x)
Force the compiler to use the function/method parameter.
Definition BaseClass.h:57
DiscreteGradient::getPersistence
dataType getPersistence(const Cell &up, const Cell &down, const dataType *const scalars, const triangulationType &triangulation) const
Definition DiscreteGradient_Template.h:26
ttk::AbstractTriangulation::gradientType
std::array< std::vector< gradIdType >, 6 > gradientType
Discrete gradient internal struct.
Definition AbstractTriangulation.h:4002
ttk::Debug::debugLevel_
int debugLevel_
Definition Debug.h:379
ttk::Debug::printWrn
int printWrn(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition Debug.h:159
ttk::Debug::printErr
int printErr(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition Debug.h:149
ttk::Timer::getElapsedTime
double getElapsedTime()
Definition Timer.h:15
ttk::dcg::DiscreteGradient
TTK discreteGradient processing package.
Definition DiscreteGradient.h:251
ttk::dcg::DiscreteGradient::getCriticalPointMap
int getCriticalPointMap(const std::vector< std::pair< SimplexId, char > > &criticalPoints, std::vector< char > &isPL)
Definition DiscreteGradient.cpp:290
ttk::dcg::DiscreteGradient::getCellLowerVertex
SimplexId getCellLowerVertex(const Cell c, const triangulationType &triangulation) const
ttk::dcg::DiscreteGradient::reverseAscendingPath
int reverseAscendingPath(const std::vector< Cell > &vpath, const triangulationType &triangulation) const
Definition DiscreteGradient_Template.h:2077
ttk::dcg::DiscreteGradient::reverseDescendingPath
int reverseDescendingPath(const std::vector< Cell > &vpath, const triangulationType &triangulation) const
Definition DiscreteGradient_Template.h:2146
ttk::dcg::DiscreteGradient::setGradientGlyphs
int setGradientGlyphs(std::vector< std::array< float, 3 > > &points, std::vector< char > &points_pairOrigins, std::vector< char > &cells_pairTypes, std::vector< SimplexId > &cellsIds, std::vector< char > &cellsDimensions, const triangulationType &triangulation) const
Definition DiscreteGradient_Template.h:2404
ttk::dcg::DiscreteGradient::reverseAscendingPathOnWall
int reverseAscendingPathOnWall(const std::vector< Cell > &vpath, const triangulationType &triangulation, bool cancelReversal=false) const
Definition DiscreteGradient_Template.h:2185
ttk::dcg::DiscreteGradient::getDescendingPathThroughWall
bool getDescendingPathThroughWall(const Cell &saddle2, const Cell &saddle1, const std::vector< bool > &isVisited, std::vector< Cell > *const vpath, const triangulationType &triangulation, const bool stopIfMultiConnected=false, const bool enableCycleDetector=false) const
Definition DiscreteGradient_Template.h:1557
ttk::dcg::DiscreteGradient::isSaddle1
bool isSaddle1(const Cell &cell) const
Definition DiscreteGradient.cpp:140
ttk::dcg::DiscreteGradient::isSaddle2
bool isSaddle2(const Cell &cell) const
Definition DiscreteGradient.cpp:149
ttk::dcg::DiscreteGradient::inputScalarField_
AbstractTriangulation::gradientKeyType inputScalarField_
Definition DiscreteGradient.h:951
ttk::dcg::DiscreteGradient::simplifySaddleSaddleConnections1
int simplifySaddleSaddleConnections1(const std::vector< std::pair< SimplexId, char > > &criticalPoints, const std::vector< char > &isPL, const int iterationThreshold, const bool allowBoundary, const bool allowBruteForce, const bool returnSaddleConnectors, const triangulationType &triangulation)
Definition DiscreteGradient_Template.h:2484
ttk::dcg::DiscreteGradient::getNumberOfCells
SimplexId getNumberOfCells(const int dimension, const triangulationType &triangulation) const
Definition DiscreteGradient_Template.h:269
ttk::dcg::DiscreteGradient::isBoundary
bool isBoundary(const Cell &cell, const triangulationType &triangulation) const
Definition DiscreteGradient_Template.h:1398
ttk::dcg::DiscreteGradient::isCellCritical
bool isCellCritical(const int cellDim, const SimplexId cellId) const
Definition DiscreteGradient.cpp:174
ttk::dcg::DiscreteGradient::getPairedCell
SimplexId getPairedCell(const Cell &cell, const triangulationType &triangulation, bool isReverse=false) const
Definition DiscreteGradient_Template.h:1411
ttk::dcg::DiscreteGradient::getAscendingWall
int getAscendingWall(const Cell &cell, VisitedMask &mask, const triangulationType &triangulation, std::vector< Cell > *const wall=nullptr, std::vector< SimplexId > *const saddles=nullptr) const
Definition DiscreteGradient_Template.h:2012
ttk::dcg::DiscreteGradient::reverseDescendingPathOnWall
int reverseDescendingPathOnWall(const std::vector< Cell > &vpath, const triangulationType &triangulation) const
Definition DiscreteGradient_Template.h:2238
ttk::dcg::DiscreteGradient::filterSaddleConnectors
int filterSaddleConnectors(const bool allowBoundary, const triangulationType &triangulation)
Definition DiscreteGradient_Template.h:2590
ttk::dcg::DiscreteGradient::getCellGreaterVertex
SimplexId getCellGreaterVertex(const Cell c, const triangulationType &triangulation) const
ttk::dcg::DiscreteGradient::getCriticalPoints
int getCriticalPoints(std::array< std::vector< SimplexId >, 4 > &criticalCellsByDim, const triangulationType &triangulation) const
Definition DiscreteGradient_Template.h:201
ttk::dcg::DiscreteGradient::getAscendingPathThroughWall
bool getAscendingPathThroughWall(const Cell &saddle1, const Cell &saddle2, const std::vector< bool > &isVisited, std::vector< Cell > *const vpath, const triangulationType &triangulation, const bool stopIfMultiConnected=false, const bool enableCycleDetector=false, bool *const cycleFound=nullptr) const
Definition DiscreteGradient_Template.h:1790
ttk::dcg::DiscreteGradient::getPersistence
dataType getPersistence(const Cell &up, const Cell &down, const dataType *const scalars, const triangulationType &triangulation) const
Definition DiscreteGradient_Template.h:26
ttk::dcg::DiscreteGradient::setCriticalPoints
int setCriticalPoints(const std::array< std::vector< SimplexId >, 4 > &criticalCellsByDim, std::vector< std::array< float, 3 > > &points, std::vector< char > &cellDimensions, std::vector< SimplexId > &cellIds, std::vector< char > &isOnBoundary, std::vector< SimplexId > &PLVertexIdentifiers, const triangulationType &triangulation) const
Definition DiscreteGradient_Template.h:126
ttk::dcg::DiscreteGradient::getDescendingPath
int getDescendingPath(const Cell &cell, std::vector< Cell > &vpath, const triangulationType &triangulation) const
Definition DiscreteGradient_Template.h:1505
ttk::dcg::DiscreteGradient::simplifySaddleSaddleConnections2
int simplifySaddleSaddleConnections2(const std::vector< std::pair< SimplexId, char > > &criticalPoints, const std::vector< char > &isPL, const int iterationThreshold, const bool allowBoundary, const bool allowBruteForce, const bool returnSaddleConnectors, const triangulationType &triangulation)
Definition DiscreteGradient_Template.h:2537
ttk::dcg::DiscreteGradient::buildGradient
int buildGradient(const triangulationType &triangulation, bool bypassCache=false, const std::vector< bool > *updateMask=nullptr)
Definition DiscreteGradient_Template.h:46
ttk::dcg::DiscreteGradient::detectGradientCycle
bool detectGradientCycle(const Cell &cell, const triangulationType &triangulation) const
Definition DiscreteGradient_Template.h:1908
ttk::dcg::DiscreteGradient::getAscendingPath
int getAscendingPath(const Cell &cell, std::vector< Cell > &vpath, const triangulationType &triangulation, const bool enableCycleDetector=false) const
Definition DiscreteGradient_Template.h:1663
ttk::dcg::DiscreteGradient::localGradient_
AbstractTriangulation::gradientType localGradient_
Definition DiscreteGradient.h:949
ttk::dcg::DiscreteGradient::getDescendingWall
int getDescendingWall(const Cell &cell, VisitedMask &mask, const triangulationType &triangulation, std::vector< Cell > *const wall=nullptr, std::vector< SimplexId > *const saddles=nullptr) const
Definition DiscreteGradient_Template.h:1949
ttk::dcg::DiscreteGradient::gradient_
AbstractTriangulation::gradientType * gradient_
Definition DiscreteGradient.h:954
ttk::dcg::DiscreteGradient::inputOffsets_
const SimplexId * inputOffsets_
Definition DiscreteGradient.h:955
ttk::dcg::DiscreteGradient::setReturnSaddleConnectors
void setReturnSaddleConnectors(const bool &returnSaddleConnectors)
Definition DiscreteGradient.h:310
ttk::ftm::FTMTree_CT::setDebugLevel
int setDebugLevel(const int &d) override
Definition FTMTree_CT.h:79
ttk::ftm::FTMTree_MT
Definition FTMTree_MT.h:83
ttk::ftm::Node::getVertexId
SimplexId getVertexId() const
Definition FTMNode.h:54
ttk::MPIsize_
COMMON_EXPORTS int MPIsize_
Definition BaseClass.cpp:10
ttk::SimplexId
int SimplexId
Identifier type for simplices of any dimension.
Definition DataTypes.h:22
ttk::MPIrank_
COMMON_EXPORTS int MPIrank_
Definition BaseClass.cpp:9
end
T end(std::pair< T, T > &p)
Definition ripser.cpp:503
ttk::VisitedMask
Auto-cleaning re-usable graph propagations data structure.
Definition VisitedMask.h:27
ttk::VisitedMask::isVisited_
std::vector< bool > & isVisited_
Definition VisitedMask.h:28
ttk::VisitedMask::visitedIds_
std::vector< SimplexId > & visitedIds_
Definition VisitedMask.h:29
ttk::dcg::CellExt
Extended Cell structure for processLowerStars.
Definition DiscreteGradient.h:62
ttk::dcg::CellExt::lowVerts_
const std::array< SimplexId, 3 > lowVerts_
Definition DiscreteGradient.h:74
ttk::dcg::CellExt::faces_
const std::array< uint8_t, 3 > faces_
Definition DiscreteGradient.h:78
printMsg
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/| (_) |"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)
Generated on Fri Apr 10 2026 10:02:07 for TTK by doxygen 1.13.2