TTK
Loading...
Searching...
No Matches
MorseSmaleQuadrangulation.h
Go to the documentation of this file.
1
18
19#pragma once
20
21// base code includes
23#include <Dijkstra.h>
24#include <Geometry.h>
25#include <Triangulation.h>
26
27#include <array>
28#include <cmath>
29#include <map>
30#include <numeric>
31#include <queue>
32#include <set>
33
34namespace ttk {
35
36 class MorseSmaleQuadrangulation : virtual public Debug {
37 public:
39 this->setDebugMsgPrefix("MorseSmaleQuadrangulation");
40 }
41
42 inline void setCriticalPoints(const unsigned int number,
43 void *const points,
44 void *const ids,
45 void *const cellIds,
46 void *const type) {
47 criticalPointsNumber_ = number;
48 criticalPoints_ = static_cast<float *>(points);
49 criticalPointsIdentifier_ = static_cast<SimplexId *>(ids);
50 criticalPointsCellIds_ = static_cast<SimplexId *>(cellIds);
51 criticalPointsType_ = static_cast<unsigned char *>(type);
52 }
53
54 inline void setSeparatrices(const unsigned int number,
55 void *const cellIds,
56 void *const cellDims,
57 void *const mask,
58 void *const points) {
59 separatriceNumber_ = number;
60 sepCellIds_ = static_cast<SimplexId *>(cellIds);
61 sepCellDims_ = static_cast<unsigned char *>(cellDims);
62 sepMask_ = static_cast<unsigned char *>(mask);
63 sepPoints_ = static_cast<float *>(points);
64 }
65
66 inline void setDualQuadrangulation(const bool input) {
67 DualQuadrangulation = input;
68 }
69 inline void setShowResError(const bool value) {
70 ShowResError = value;
71 }
72 inline void
74 if(triangl != nullptr) {
78 verticesNumber_ = triangl->getNumberOfVertices();
80 }
81 }
82
83 template <typename triangulationType>
84 int execute(const triangulationType &triangulation);
85
86 private:
95 template <typename triangulationType>
96 size_t findSeparatrixMiddle(const size_t a,
97 const size_t b,
98 const triangulationType &triangulation);
99
109 int findSepsVertices(const std::vector<size_t> &seps,
110 std::vector<LongSimplexId> &srcs,
111 std::vector<LongSimplexId> &dsts) const;
112
122 template <typename triangulationType>
123 int quadrangulate(size_t &ndegen, const triangulationType &triangulation);
124
133 int dualQuadrangulate();
134
143 template <typename triangulationType>
144 int subdivise(const triangulationType &triangulation);
145
155 template <typename triangulationType>
156 int detectCellSeps(const triangulationType &triangulation);
157
164 template <typename triangulationType>
165 bool checkSurfaceCloseness(const triangulationType &triangulation) const;
166
170 void clearData();
171
175 using Quad = std::array<LongSimplexId, 4>;
176
184 template <typename triangulationType>
185 int subdiviseDegenerateQuads(std::vector<Quad> &outputSubd,
186 const triangulationType &triangulation);
187
188 // number of vertices in triangulation
189 SimplexId verticesNumber_{};
190
191 // number of critical points from the Morse-Smale complex
192 SimplexId criticalPointsNumber_{};
193 // critical points 3d coordinates
194 float *criticalPoints_{};
195 // mapping points id -> cells id
196 SimplexId *criticalPointsCellIds_{};
197 // mapping point id -> TTK identifier
198 SimplexId *criticalPointsIdentifier_{};
199 // critical point type: 0 minimum, 1 saddle point, 2 maximum
200 unsigned char *criticalPointsType_{};
201
202 // number of separatrices data
203 SimplexId separatriceNumber_{};
204 // separatrices points cellIds (to be linked to critical points cellIds)
205 SimplexId *sepCellIds_{};
206 // separatrices mask scalar field (0 for critical points, 1 otherwise)
207 unsigned char *sepMask_{};
208 // separatrices cell dimension: 0 for vertices, 1 for edges, 2 for triangles
209 unsigned char *sepCellDims_{};
210 // separatrices points
211 float *sepPoints_{};
212
213 // index of separatrices beginnings in separatrices arrays
214 std::vector<size_t> sepBegs_{};
215 // index of separatrices endings in separatrices arrays
216 std::vector<size_t> sepEnds_{};
217 // separatrices middles index in output points array
218 std::vector<SimplexId> sepMids_{};
219 // sub-segmentation of Morse-Smale cells
220 std::vector<SimplexId> morseSeg_{};
221 // morseSeg_ id -> separatrices that border quads id
222 std::vector<std::pair<SimplexId, std::vector<size_t>>> quadSeps_{};
223
224 protected:
225 // array of output quads
226 std::vector<Quad> outputCells_{};
227 // array of output vertices (generated middles of duplicated separatrices)
228 std::vector<float> outputPoints_{};
229 // array of output vertices identifiers
230 std::vector<SimplexId> outputPointsIds_{};
231 // 0: critical points, 1: edge middle, 2: quad barycenter
232 std::vector<SimplexId> outputPointsTypes_{};
233 // for critical points, their id, for sep middles, the sep id and
234 // for barycenters the parent quad id
235 std::vector<SimplexId> outputPointsCells_{};
236 // triangulation subdivision to detect cells
238
239 // if dual quadrangulation
241 // display result despite error
242 bool ShowResError{false};
243 };
244} // namespace ttk
245
246template <typename triangulationType>
247int ttk::MorseSmaleQuadrangulation::detectCellSeps(
248 const triangulationType &triangulation) {
249
250 ExplicitTriangulation newT{};
251 bs.execute(triangulation, newT);
252
253 newT.setDebugLevel(this->debugLevel_);
254 newT.setThreadNumber(this->threadNumber_);
255 newT.preconditionEdges();
256 newT.preconditionVertexNeighbors();
257 newT.preconditionVertexEdges();
258 newT.preconditionVertexTriangles();
259 newT.preconditionEdgeTriangles();
260 newT.preconditionTriangleEdges();
261
262 auto nEdges = triangulation.getNumberOfEdges();
263
264 // store separatrix index on subdivised triangulation edges
265 std::vector<SimplexId> edgeOnSep(newT.getNumberOfEdges(), -1);
266 // count the number of critical points encountered
267 size_t critPoints{0};
268 if(sepMask_[0] == 0) {
269 critPoints++;
270 }
271
272 // id of critical point in new triangulation
273 auto critPointId = [&](const SimplexId a) -> SimplexId {
274 if(sepCellDims_[a] == 0) {
275 return sepCellIds_[a];
276 }
277 if(sepCellDims_[a] == 1) {
278 return verticesNumber_ + sepCellIds_[a];
279 }
280 if(sepCellDims_[a] == 2) {
281 return verticesNumber_ + nEdges + sepCellIds_[a];
282 }
283 return -1;
284 };
285
286 // init first vertex id
287 SimplexId prev{critPointId(0)};
288
289 // put separatrices indices onto edges
290 // begin loop at one to get the first edge
291 for(SimplexId i = 1; i < separatriceNumber_; ++i) {
292
293 // for computing the separatrix id
294 if(sepMask_[i] == 0) {
295 critPoints++;
296 // beginning of a new separatrix
297 if(critPoints % 2 == 1) {
298 prev = critPointId(i);
299 continue;
300 }
301 }
302
303 // current separatrix id is critPoints // 2
304 auto currSepId
305 = (critPoints % 2 == 0) ? critPoints / 2 - 1 : critPoints / 2;
306
307 // current point
308 SimplexId curr{critPointId(i)};
309
310 // get edge id
311 auto ne = newT.getVertexEdgeNumber(curr);
312 for(SimplexId j = 0; j < ne; ++j) {
313 SimplexId e{};
314 newT.getVertexEdge(curr, j, e);
315 SimplexId e0{}, e1{};
316 newT.getEdgeVertex(e, 0, e0);
317 newT.getEdgeVertex(e, 1, e1);
318
319 if((e0 == prev && e1 == curr) || (e1 == prev && e0 == curr)) {
320 edgeOnSep[e] = currSepId;
321 break;
322 }
323 }
324
325 prev = curr;
326 }
327
328 // if vertex is saddle in the subdivised triangulation
329 std::vector<bool> isSaddle(newT.getNumberOfVertices(), false);
330
331 for(SimplexId i = 0; i < criticalPointsNumber_; ++i) {
332 // keep only saddle points
333 if(criticalPointsType_[i] != 1) {
334 continue;
335 }
336 isSaddle[verticesNumber_ + criticalPointsCellIds_[i]] = true;
337 }
338
339 auto getTriangleEdges
340 = [&](const SimplexId tr, SimplexId &e0, SimplexId &e1, SimplexId &e2) {
341 newT.getTriangleEdge(tr, 0, e0);
342 newT.getTriangleEdge(tr, 1, e1);
343 newT.getTriangleEdge(tr, 2, e2);
344 };
345
346 // get the indices of the two separatrices around a triangle
347 // (having a saddle point as vertex)
348 auto sepIdAroundTriangle = [&](const SimplexId tr) {
349 SimplexId e0{}, e1{}, e2{};
350 getTriangleEdges(tr, e0, e1, e2);
351
352 std::set<SimplexId> sepId{};
353 if(edgeOnSep[e0] != -1 && edgeOnSep[e1] != -1) {
354 sepId.emplace(edgeOnSep[e0]);
355 sepId.emplace(edgeOnSep[e1]);
356 } else if(edgeOnSep[e1] != -1 && edgeOnSep[e2] != -1) {
357 sepId.emplace(edgeOnSep[e1]);
358 sepId.emplace(edgeOnSep[e2]);
359 } else if(edgeOnSep[e0] != -1 && edgeOnSep[e2] != -1) {
360 sepId.emplace(edgeOnSep[e0]);
361 sepId.emplace(edgeOnSep[e2]);
362 }
363 return sepId;
364 };
365
366 auto hasTriangleSaddle = [&](const SimplexId tr) {
367 SimplexId a, b, c;
368 newT.getTriangleVertex(tr, 0, a);
369 newT.getTriangleVertex(tr, 1, b);
370 newT.getTriangleVertex(tr, 2, c);
371 return (isSaddle[a] || isSaddle[b] || isSaddle[c]);
372 };
373
374 // propagate from triangles around saddle
375 std::vector<SimplexId> processed(newT.getNumberOfTriangles(), -1);
376
377 // store the saddles id
378 std::vector<SimplexId> saddles{};
379
380 for(SimplexId i = 0; i < criticalPointsNumber_; ++i) {
381 // keep only saddle points
382 if(criticalPointsType_[i] == 1) {
383 saddles.emplace_back(verticesNumber_ + criticalPointsCellIds_[i]);
384 }
385 }
386
387 // look around the saddle points
388 for(size_t i = 0; i < saddles.size(); ++i) {
389 SimplexId saddle = saddles[i];
390
391 auto sadtri = newT.getVertexTriangleNumber(saddle);
392
393#ifdef TTK_ENABLE_OPENMP
394#pragma omp parallel for num_threads(threadNumber_)
395#endif // TTK_ENABLE_OPENMP
396 for(SimplexId j = 0; j < sadtri; ++j) {
397 std::queue<SimplexId> toProcess{};
398 SimplexId tr;
399 newT.getVertexTriangle(saddle, j, tr);
400
401 std::set<SimplexId> sepIdBeg = sepIdAroundTriangle(tr);
402 // current iteration id
403 SimplexId iter = i * sadtri + j;
404
405 toProcess.push(tr);
406
407 while(!toProcess.empty()) {
408 auto curr = toProcess.front();
409 toProcess.pop();
410
411 // mark current triangle, skip already processed
412#ifdef TTK_ENABLE_OPENMP
413 if(processed[curr] == -1 || processed[curr] > iter) {
414#pragma omp atomic write
415 processed[curr] = iter;
416 }
417 // skip if curr marked by thread
418 else if(processed[curr] == iter) {
419 continue;
420 }
421 // stop BFS if thread with higher iteration id has reached curr
422 else if(processed[curr] < iter) {
423 break;
424 }
425#else
426 if(processed[curr] != -1) {
427 continue;
428 }
429 processed[curr] = 0;
430#endif // TTK_ENABLE_OPENMP
431
432 // check for saddle at vertices
433 bool hasSaddle = hasTriangleSaddle(curr);
434
435 // check for separatrices on edges
436 auto sepIdEnd = sepIdAroundTriangle(curr);
437
438 if(hasSaddle && sepIdEnd.size() == 2) {
439 // check that seps are different from beginning
440 std::vector<size_t> cellSeps{};
441 std::set_union(sepIdBeg.begin(), sepIdBeg.end(), sepIdEnd.begin(),
442 sepIdEnd.end(), std::back_inserter(cellSeps));
443 if(cellSeps.size() > 2) {
444 // found it
445#ifdef TTK_ENABLE_OPENMP
446#pragma omp critical
447#endif // TTK_ENABLE_OPENMP
448 {
449 // keep indices in sync
450 quadSeps_.emplace_back(iter, cellSeps);
451 }
452 }
453 }
454
455 // mark vertices from the original mesh with cell id
456 for(SimplexId k = 0; k < 3; ++k) {
457 SimplexId vert{};
458 newT.getTriangleVertex(curr, k, vert);
459 if(vert >= verticesNumber_) {
460 continue;
461 }
462 bool vertOnSep = false;
463 for(SimplexId l = 0; l < newT.getVertexEdgeNumber(vert); ++l) {
464 SimplexId e{};
465 newT.getVertexEdge(vert, l, e);
466 if(edgeOnSep[e] != -1) {
467 vertOnSep = true;
468 break;
469 }
470 }
471 if(!vertOnSep) {
472#ifdef TTK_ENABLE_OPENMP
473#pragma omp atomic write
474#endif // TTK_ENABLE_OPENMP
475 morseSeg_[vert] = iter;
476 }
477 }
478
479 // look for neighboring triangles
480 std::array<SimplexId, 3> edges{};
481 getTriangleEdges(curr, edges[0], edges[1], edges[2]);
482 for(const auto e : edges) {
483 // do not cross separatrices
484 if(edgeOnSep[e] != -1) {
485 continue;
486 }
487 for(SimplexId k = 0; k < newT.getEdgeTriangleNumber(e); ++k) {
488 SimplexId neigh{};
489 newT.getEdgeTriangle(e, k, neigh);
490 // push only non processed triangles
491 if(processed[neigh] == -1) {
492 toProcess.push(neigh);
493 }
494#ifdef TTK_ENABLE_OPENMP
495 // push if neigh marked during a previous iteration
496 else if(processed[neigh] > iter) {
497 toProcess.push(neigh);
498 }
499 // stop pushing neighbors if curr marked in a newer iteration
500 else if(processed[neigh] < iter) {
501 break;
502 }
503#endif // TTK_ENABLE_OPENMP
504 }
505 }
506 }
507 }
508 }
509
510 if(this->threadNumber_ > 1) {
511 // sort quadSeps_ according to cellIds_ to get a deterministic
512 // output when filled in parallel
513 TTK_PSORT(this->threadNumber_, quadSeps_.begin(), quadSeps_.end());
514 // (by default, pairs are sorted by their first element)
515 }
516
517 return 0;
518}
519
520template <typename triangulationType>
521int ttk::MorseSmaleQuadrangulation::quadrangulate(
522 size_t &ndegen, const triangulationType &triangulation) {
523 // quadrangle vertices are either extrema or saddle points
524
525 // separatrices bounds indices and cell ids
526 std::vector<size_t> sepFlatEdges{};
527
528 for(SimplexId i = 0; i < separatriceNumber_; ++i) {
529 if(sepMask_[i] == 1) {
530 continue;
531 }
532 sepFlatEdges.emplace_back(i);
533 }
534
535 // number of separatrices
536 auto numSeps = sepFlatEdges.size() / 2;
537
538 // clear data members
539 sepBegs_.resize(numSeps);
540 sepEnds_.resize(numSeps);
541 morseSeg_.resize(verticesNumber_);
542 std::fill(morseSeg_.begin(), morseSeg_.end(), -1);
543 quadSeps_.clear();
544
545 // fill in data arrays
546 for(size_t i = 0; i < numSeps; ++i) {
547 // separatrices bounds
548 sepBegs_[i] = sepFlatEdges[2 * i];
549 sepEnds_[i] = sepFlatEdges[2 * i + 1];
550 }
551
552 detectCellSeps(triangulation);
553
554 outputCells_.reserve(quadSeps_.size());
555
556 for(const auto &qsp : quadSeps_) {
557
558 const auto &qs{qsp.second};
559
560 std::vector<LongSimplexId> srcs{};
561 std::vector<LongSimplexId> dsts{};
562
563 findSepsVertices(qs, srcs, dsts);
564
565 // remove duplicates
566 std::set<LongSimplexId> srcs_set(srcs.begin(), srcs.end());
567 std::set<LongSimplexId> dsts_set(dsts.begin(), dsts.end());
568 srcs.assign(srcs_set.begin(), srcs_set.end());
569 dsts.assign(dsts_set.begin(), dsts_set.end());
570
571 // filter out separatrices whose sources are not in contact with
572 // the current cell
573
574 bool found = true;
575
576 if(dsts.size() != 2) {
577 found = false;
578 }
579
580 if(srcs.size() == 2) {
581 outputCells_.emplace_back(Quad{dsts[0], srcs[0], dsts[1], srcs[1]});
582 } else if(srcs.size() == 1) {
583 outputCells_.emplace_back(Quad{dsts[0], srcs[0], dsts[1], srcs[0]});
584 ndegen++;
585 } else {
586 found = false;
587 }
588
589 if(!found) {
590 this->printMsg("Missing quadrangle", ttk::debug::Priority::DETAIL);
591 }
592 }
593
594 return 0;
595}
596
597template <typename triangulationType>
598size_t ttk::MorseSmaleQuadrangulation::findSeparatrixMiddle(
599 const size_t a, const size_t b, const triangulationType &triangulation) {
600
601 const int dim = 3;
602
603 std::vector<float> distFromA(b - a + 1);
604 std::array<float, dim> prev{}, curr{};
605
606 if(distFromA.empty()) {
607 return 0;
608 }
609
610 curr[0] = sepPoints_[dim * a];
611 curr[1] = sepPoints_[dim * a + 1];
612 curr[2] = sepPoints_[dim * a + 2];
613
614 // integrate distances at every point of this separatrix
615 for(size_t i = 1; i < b - a + 1; ++i) {
616 std::swap(curr, prev);
617 curr[0] = sepPoints_[dim * (a + i)];
618 curr[1] = sepPoints_[dim * (a + i) + 1];
619 curr[2] = sepPoints_[dim * (a + i) + 2];
620 distFromA[i]
621 = distFromA[i - 1] + ttk::Geometry::distance(&curr[0], &prev[0]);
622 }
623
624 auto distAB = distFromA.back();
625 for(auto &el : distFromA) {
626 el = std::abs(el - distAB / 2.0);
627 }
628
629 // index in separatrices point data array of separatrix middle
630 auto pos = a + std::min_element(distFromA.begin(), distFromA.end())
631 - distFromA.begin();
632
633 // new point!
634 outputPoints_.emplace_back(sepPoints_[dim * pos]);
635 outputPoints_.emplace_back(sepPoints_[dim * pos + 1]);
636 outputPoints_.emplace_back(sepPoints_[dim * pos + 2]);
637
638 SimplexId id = pos;
639
640 // new point identifier (on the triangular mesh)
641 switch(sepCellDims_[pos]) {
642 case 0:
643 outputPointsIds_.emplace_back(sepCellIds_[pos]);
644 break;
645 case 1: {
646 // take the first vertex of the edge
647 triangulation.getEdgeVertex(sepCellIds_[pos], 0, id);
648 outputPointsIds_.emplace_back(id);
649 break;
650 }
651 case 2: {
652 // take the first vertex of the triangle
653 triangulation.getTriangleVertex(sepCellIds_[pos], 0, id);
654 outputPointsIds_.emplace_back(id);
655 break;
656 }
657 default:
658 break;
659 }
660
661 outputPointsTypes_.emplace_back(1);
662
663 return id;
664}
665
666template <typename triangulationType>
667int ttk::MorseSmaleQuadrangulation::subdiviseDegenerateQuads(
668 std::vector<Quad> &outputSubd, const triangulationType &triangulation) {
669
670 for(size_t i = 0; i < outputCells_.size(); ++i) {
671 auto q = outputCells_[i];
672 auto seps = quadSeps_[i].second;
673
674 // don't deal with normal quadrangles
675 if(q[1] != q[3]) {
676 continue;
677 }
678
679 std::vector<LongSimplexId> srcs{};
680 std::vector<LongSimplexId> dsts{};
681
682 findSepsVertices(seps, srcs, dsts);
683
684 // identify the extremum that is twice dest
685 int count_vi = 0, count_vk = 0;
686 for(const auto &s : dsts) {
687 if(s == q[0]) {
688 count_vi++;
689 }
690 if(s == q[2]) {
691 count_vk++;
692 }
693 }
694 // extremum index
695 LongSimplexId vert2Seps = count_vi > count_vk ? q[0] : q[2];
696 LongSimplexId vert1Sep = count_vi > count_vk ? q[2] : q[0];
697 // the two seps from j to vert2Seps
698 std::vector<size_t> borderseps{};
699 for(size_t j = 0; j < seps.size(); ++j) {
700 if(dsts[j] == vert2Seps && srcs[j] == q[1]) {
701 borderseps.emplace_back(seps[j]);
702 }
703 }
704
705 if(borderseps.size() < 2) {
706 continue;
707 }
708
709 // find a midpoint between the two extrema on the triangulation
710
711 std::vector<SimplexId> boundi{criticalPointsIdentifier_[q[0]]};
712 std::vector<SimplexId> boundk{criticalPointsIdentifier_[q[2]]};
713 std::array<std::vector<float>, 6> outputDists{};
714
716 criticalPointsIdentifier_[q[0]], triangulation, outputDists[0], boundk);
718 criticalPointsIdentifier_[q[2]], triangulation, outputDists[1], boundi);
719
720 auto inf = std::numeric_limits<float>::infinity();
721 std::vector<float> sum(outputDists[0].size(), inf);
722
723 for(size_t j = 0; j < sum.size(); ++j) {
724 auto m = outputDists[0][j];
725 auto n = outputDists[1][j];
726 if(m == inf || n == inf) {
727 continue;
728 }
729 // cost to minimize
730 sum[j] = m + n + std::abs(m - n);
731 }
732
733 auto insertNewPoint
734 = [&](const SimplexId a, const size_t idx, const SimplexId type) {
735 float x, y, z;
736 triangulation.getVertexPoint(a, x, y, z);
737 outputPoints_.emplace_back(x);
738 outputPoints_.emplace_back(y);
739 outputPoints_.emplace_back(z);
740 outputPointsIds_.emplace_back(a);
741 outputPointsTypes_.emplace_back(type);
742 outputPointsCells_.emplace_back(idx);
743 return outputPointsIds_.size() - 1;
744 };
745
746 auto v0 = std::min_element(sum.begin(), sum.end()) - sum.begin();
747 auto v0Pos = static_cast<LongSimplexId>(insertNewPoint(v0, i, 3));
748
749 // find two other points
750
751 std::vector<SimplexId> bounds{criticalPointsIdentifier_[q[0]],
752 criticalPointsIdentifier_[q[1]],
753 criticalPointsIdentifier_[q[2]]};
754
755 auto m0Pos = sepMids_[borderseps[0]];
756 auto m1Pos = sepMids_[borderseps[1]];
757 auto m0 = outputPointsIds_[m0Pos];
758 auto m1 = outputPointsIds_[m1Pos];
759
760 Dijkstra::shortestPath(criticalPointsIdentifier_[vert1Sep], triangulation,
761 outputDists[2], bounds);
762 Dijkstra::shortestPath(v0, triangulation, outputDists[3], bounds);
763 Dijkstra::shortestPath(m0, triangulation, outputDists[4], bounds);
764 Dijkstra::shortestPath(m1, triangulation, outputDists[5], bounds);
765
766 std::fill(sum.begin(), sum.end(), inf);
767
768 for(size_t j = 0; j < sum.size(); ++j) {
769 auto m = outputDists[2][j];
770 auto n = outputDists[3][j];
771 auto o = outputDists[4][j];
772 if(m == inf || n == inf || o == inf) {
773 continue;
774 }
775 // cost to minimize
776 sum[j] = m + n + o + std::abs(m - n) + std::abs(m - o) + std::abs(n - o);
777 }
778
779 auto v1 = std::min_element(sum.begin(), sum.end()) - sum.begin();
780 auto v1Pos = static_cast<LongSimplexId>(insertNewPoint(v1, i, 4));
781 std::fill(sum.begin(), sum.end(), inf);
782
783 for(size_t j = 0; j < sum.size(); ++j) {
784 auto m = outputDists[2][j];
785 auto n = outputDists[3][j];
786 auto o = outputDists[5][j];
787 if(m == inf || n == inf || o == inf) {
788 continue;
789 }
790 // cost to minimize
791 sum[j] = m + n + o + std::abs(m - n) + std::abs(m - o) + std::abs(n - o);
792 }
793
794 auto v2 = std::min_element(sum.begin(), sum.end()) - sum.begin();
795 auto v2Pos = static_cast<LongSimplexId>(insertNewPoint(v2, i, 4));
796
797 outputSubd.emplace_back(Quad{vert2Seps, m0Pos, v1Pos, v0Pos});
798 outputSubd.emplace_back(Quad{vert2Seps, m1Pos, v2Pos, v0Pos});
799 outputSubd.emplace_back(Quad{q[1], m0Pos, v1Pos, vert1Sep});
800 outputSubd.emplace_back(Quad{q[1], m1Pos, v2Pos, vert1Sep});
801 outputSubd.emplace_back(Quad{vert1Sep, v1Pos, v0Pos, v2Pos});
802 }
803 return 0;
804}
805
806template <typename triangulationType>
807int ttk::MorseSmaleQuadrangulation::subdivise(
808 const triangulationType &triangulation) {
809
810 // separatrices middles index in output points array
811 sepMids_.resize(sepBegs_.size());
812
813 for(size_t i = 0; i < sepMids_.size(); ++i) {
814 // separatrices middles
815 sepMids_[i] = outputPoints_.size() / 3; // before insertion at next line
816 findSeparatrixMiddle(sepBegs_[i], sepEnds_[i], triangulation);
817 outputPointsCells_.emplace_back(i);
818 }
819
820 // for each output quad, its barycenter position in outputPoints_
821 std::vector<size_t> cellBary(outputCells_.size());
822
823 std::array<std::vector<float>, 4> outputDists{};
824
825 // hold quad subdivision
826 decltype(outputCells_) outputSubd{};
827 outputSubd.reserve(4 * outputCells_.size());
828
829 for(size_t i = 0; i < outputCells_.size(); ++i) {
830 auto q = outputCells_[i];
831 auto seps = quadSeps_[i].second;
832
833 // skip degenerate case here
834 if(q[1] == q[3]) {
835 continue;
836 }
837
838 std::vector<LongSimplexId> sepMids(seps.size());
839 std::vector<SimplexId> midsNearestVertex(seps.size());
840 for(size_t j = 0; j < seps.size(); ++j) {
841 sepMids[j] = sepMids_[seps[j]];
842 midsNearestVertex[j] = outputPointsIds_[sepMids_[seps[j]]];
843 }
844
845 // find barycenter of current cell (c.f. QuadrangulationSubdivision.cpp)
846
847 // bound Dijkstra by parent quad vertices
848 std::vector<SimplexId> bounds{
849 criticalPointsIdentifier_[q[0]], criticalPointsIdentifier_[q[1]],
850 criticalPointsIdentifier_[q[2]], criticalPointsIdentifier_[q[3]]};
851
852 // Dijkstra propagation mask
853 std::vector<bool> mask(morseSeg_.size(), false);
854 // restrict Dijkstra propagation to current cell
855 for(size_t j = 0; j < morseSeg_.size(); ++j) {
856 if(morseSeg_[j] == quadSeps_[i].first) {
857 mask[j] = true;
858 }
859 }
860 // also allow to propagate on separatrices
861 for(const auto s : seps) {
862 for(size_t j = sepBegs_[s]; j <= sepEnds_[s]; ++j) {
863 if(sepCellDims_[j] == 1) {
864 auto e = sepCellIds_[j];
865 SimplexId e0{}, e1{};
866 triangulation.getEdgeVertex(e, 0, e0);
867 triangulation.getEdgeVertex(e, 1, e1);
868 mask[e0] = true;
869 mask[e1] = true;
870 }
871 }
872 }
873
874#ifdef TTK_ENABLE_OPENMP
875#pragma omp parallel for num_threads(threadNumber_)
876#endif // TTK_ENABLE_OPENMP
877 for(size_t j = 0; j < outputDists.size(); ++j) {
878 Dijkstra::shortestPath(midsNearestVertex[j], triangulation,
879 outputDists.at(j), std::vector<SimplexId>(), mask);
880 }
881
882 auto inf = std::numeric_limits<float>::infinity();
883 std::vector<float> sum(outputDists[0].size(), inf);
884
885#ifdef TTK_ENABLE_OPENMP
886#pragma omp parallel for num_threads(threadNumber_)
887#endif // TTK_ENABLE_OPENMP
888 for(size_t j = 0; j < sum.size(); ++j) {
889 // skip if vertex j not in cell i
890 if(morseSeg_[j] != quadSeps_[i].first) {
891 continue;
892 }
893 auto m = outputDists[0][j];
894 auto n = outputDists[1][j];
895 auto o = outputDists[2][j];
896 auto p = outputDists[3][j];
897 if(m == inf || n == inf || o == inf || p == inf) {
898 continue;
899 }
900 // cost to minimize
901 sum[j] = m + n + o + p + std::abs(m - o) + std::abs(n - p);
902 }
903
904 size_t verticesInCell
905 = verticesNumber_
906 - std::count(
907 sum.begin(), sum.end(), std::numeric_limits<float>::infinity());
908
909 const size_t thresholdVertsInCell{50};
910 if(verticesInCell <= thresholdVertsInCell) {
911 this->printMsg("Small cell detected");
912 }
913
914 SimplexId baryId{};
915 if(verticesInCell == 0) {
916 this->printMsg("Barycenter of cell " + std::to_string(i) + " not found");
917 // snap bary on sepMids[0]
918 baryId = outputPointsIds_[sepMids[0]];
919 } else {
920 baryId = std::min_element(sum.begin(), sum.end()) - sum.begin();
921 }
922
923 LongSimplexId baryPos = outputPointsIds_.size();
924 {
925 float x, y, z;
926 triangulation.getVertexPoint(baryId, x, y, z);
927 outputPoints_.emplace_back(x);
928 outputPoints_.emplace_back(y);
929 outputPoints_.emplace_back(z);
930 outputPointsIds_.emplace_back(baryId);
931 outputPointsTypes_.emplace_back(2);
932 outputPointsCells_.emplace_back(i);
933 }
934
935 std::vector<LongSimplexId> srcs{};
936 std::vector<LongSimplexId> dsts{};
937
938 findSepsVertices(seps, srcs, dsts);
939
940 auto sepsQuadVertex = [&](const size_t a, const size_t b) {
941 if(srcs[a] == srcs[b]) {
942 outputSubd.emplace_back(Quad{srcs[a], sepMids[a], baryPos, sepMids[b]});
943 }
944 if(dsts[a] == dsts[b]) {
945 outputSubd.emplace_back(Quad{dsts[a], sepMids[a], baryPos, sepMids[b]});
946 }
947 };
948
949 // test every pair of current quad seps for a common vertex, if yes, use
950 // their middles and the quad barycenter
951 sepsQuadVertex(0, 1);
952 sepsQuadVertex(0, 2);
953 sepsQuadVertex(0, 3);
954 sepsQuadVertex(1, 2);
955 sepsQuadVertex(1, 3);
956 sepsQuadVertex(2, 3);
957 }
958
959 subdiviseDegenerateQuads(outputSubd, triangulation);
960
961 // overwrite old quads
962 outputCells_ = std::move(outputSubd);
963
964 return 0;
965}
966
967template <typename triangulationType>
968bool ttk::MorseSmaleQuadrangulation::checkSurfaceCloseness(
969 const triangulationType &triangulation) const {
970
971 bool triangulationClosed{true};
972 // sweep over all vertices to check if one is on a boundary
973 for(SimplexId i = 0; i < verticesNumber_; ++i) {
974 if(triangulation.isVertexOnBoundary(i)) {
975 triangulationClosed = false;
976 break;
977 }
978 }
979 // quadrangles edges -> quadrangles
980 std::map<std::pair<LongSimplexId, LongSimplexId>, std::set<size_t>>
981 quadEdges{};
982
983 // sweep over quadrangulation edges
984 for(size_t i = 0; i < outputCells_.size(); ++i) {
985 auto q = outputCells_[i];
986 // store edges in order
987 if(q[0] < q[1]) {
988 quadEdges[std::make_pair(q[0], q[1])].emplace(i);
989 } else {
990 quadEdges[std::make_pair(q[1], q[0])].emplace(i);
991 }
992 if(q[1] < q[2]) {
993 quadEdges[std::make_pair(q[1], q[2])].emplace(i);
994 } else {
995 quadEdges[std::make_pair(q[2], q[1])].emplace(i);
996 }
997 if(q[2] < q[3]) {
998 quadEdges[std::make_pair(q[2], q[3])].emplace(i);
999 } else {
1000 quadEdges[std::make_pair(q[3], q[2])].emplace(i);
1001 }
1002 if(q[3] < q[0]) {
1003 quadEdges[std::make_pair(q[3], q[0])].emplace(i);
1004 } else {
1005 quadEdges[std::make_pair(q[0], q[3])].emplace(i);
1006 }
1007 }
1008
1009 bool quadrangulationClosed{true};
1010 for(const auto &e : quadEdges) {
1011 if(e.second.size() < 2) {
1012 quadrangulationClosed = false;
1013 break;
1014 }
1015 }
1016
1017 // each edge should be shared by two different quadrangles
1018 return triangulationClosed == quadrangulationClosed;
1019}
1020
1021// main routine
1022template <typename triangulationType>
1024 const triangulationType &triangulation) {
1025
1026 Timer tm;
1027
1028 // sanity check
1029 if(separatriceNumber_ == 0) {
1030 this->printErr("Unable to perform quadrangulation without separatrices");
1031 return 1;
1032 }
1033
1034 // clear output
1035 clearData();
1036 outputPoints_.resize(3 * criticalPointsNumber_);
1037 outputPointsIds_.resize(criticalPointsNumber_);
1038 outputPointsTypes_.resize(criticalPointsNumber_);
1039 outputPointsCells_.resize(criticalPointsNumber_);
1040
1041 // fill in critical points 3d coordinates and identifiers
1042 for(SimplexId i = 0; i < criticalPointsNumber_; ++i) {
1043 outputPoints_[3 * i] = criticalPoints_[3 * i];
1044 outputPoints_[3 * i + 1] = criticalPoints_[3 * i + 1];
1045 outputPoints_[3 * i + 2] = criticalPoints_[3 * i + 2];
1046 outputPointsIds_[i] = criticalPointsIdentifier_[i];
1047 outputPointsTypes_[i] = 0;
1048 outputPointsCells_[i] = i;
1049 }
1050
1051 // number of degenerate quadrangles
1052 size_t ndegen = 0;
1053
1054 // direct quadrangulation with saddle points
1055 int ret = quadrangulate(ndegen, triangulation);
1056
1057 if(ret == 0) {
1058 subdivise(triangulation);
1059 } else {
1060 // clean, log & early return
1061 clearData();
1062 this->printErr("Unable to generate a quadrangulation from the given "
1063 "Morse-Smale complex");
1064 return 1;
1065 }
1066
1067 if(DualQuadrangulation) {
1068 dualQuadrangulate();
1069 }
1070
1071 if(!checkSurfaceCloseness(triangulation)) {
1072 // log, clean & early return
1073 this->printErr("Output surface does not match input surface closeness");
1074 if(!ShowResError) {
1075 clearData();
1076 return 1;
1077 }
1078 }
1079
1080 std::string s_degen{
1081 ndegen > 0 ? "(" + std::to_string(ndegen) + " degenerated) " : ""};
1082
1083 this->printMsg("Produced " + std::to_string(this->outputCells_.size())
1084 + " quads " + s_degen + "("
1085 + std::to_string(this->outputPointsIds_.size()) + " points)",
1086 1.0, tm.getElapsedTime(), this->threadNumber_);
1087
1088 return 0;
1089}
#define TTK_PSORT(NTHREADS,...)
Parallel sort macro.
Definition: OpenMP.h:46
AbstractTriangulation is an interface class that defines an interface for efficient traversal methods...
virtual SimplexId getNumberOfVertices() const
Subdivise a triangulation according to triangle barycenter.
void preconditionTriangulation(AbstractTriangulation *const triangulation)
int execute(const triangulationType &inputTriangl, ExplicitTriangulation &outputTriangl)
int threadNumber_
Definition: BaseClass.h:95
Minimalist debugging class.
Definition: Debug.h:88
int debugLevel_
Definition: Debug.h:379
void setDebugMsgPrefix(const std::string &prefix)
Definition: Debug.h:364
TTK processing package for the topological simplification of scalar data.
std::vector< SimplexId > outputPointsCells_
std::vector< SimplexId > outputPointsIds_
void setDualQuadrangulation(const bool input)
std::vector< SimplexId > outputPointsTypes_
void setSeparatrices(const unsigned int number, void *const cellIds, void *const cellDims, void *const mask, void *const points)
void preconditionTriangulation(AbstractTriangulation *const triangl)
void setCriticalPoints(const unsigned int number, void *const points, void *const ids, void *const cellIds, void *const type)
int execute(const triangulationType &triangulation)
double getElapsedTime()
Definition: Timer.h:15
int shortestPath(const SimplexId source, const triangulationType &triangulation, std::vector< T > &outputDists, const std::vector< SimplexId > &bounds=std::vector< SimplexId >(), const std::vector< bool > &mask=std::vector< bool >())
Compute the Dijkstra shortest path from source.
Definition: Dijkstra.h:26
T distance(const T *p0, const T *p1, const int &dimension=3)
Definition: Geometry.cpp:344
The Topology ToolKit.
long long int LongSimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:15
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22
printMsg(debug::output::GREEN+"                           "+debug::output::ENDCOLOR+debug::output::GREEN+"▒"+debug::output::ENDCOLOR+debug::output::GREEN+"▒▒▒▒▒▒▒▒▒▒▒▒▒░"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)