TTK
Loading...
Searching...
No Matches
MarchingTetrahedra.h
Go to the documentation of this file.
1
26
27#pragma once
28
29// base code includes
30#include <Triangulation.h>
31
32#include <climits>
33#include <queue>
34#include <type_traits>
35
36#include <MarchingTetrahedraLookupTables.inl>
37
38using ttk::SimplexId;
39
40namespace ttk {
41 namespace mth {
50 constexpr unsigned long long int getHash(const unsigned long long int a,
51 const unsigned long long int b) {
52 return (a * b + (a * a) + (b * b) + (a * a * a) * (b * b * b))
53 % ULLONG_MAX;
54 }
55
64 inline void getCenter(const std::array<float, 3> pos0,
65 const std::array<float, 3> pos1,
66 std::array<float, 3> &incenter) {
67 incenter[0] = 0.5 * (pos0[0] + pos1[0]);
68 incenter[1] = 0.5 * (pos0[1] + pos1[1]);
69 incenter[2] = 0.5 * (pos0[2] + pos1[2]);
70 }
71
81 inline void getCenter(const std::array<float, 3> pos0,
82 const std::array<float, 3> pos1,
83 const std::array<float, 3> pos2,
84 std::array<float, 3> &incenter) {
85 incenter[0] = 0.3333 * (pos0[0] + pos1[0] + pos2[0]);
86 incenter[1] = 0.3333 * (pos0[1] + pos1[1] + pos2[1]);
87 incenter[2] = 0.3333 * (pos0[2] + pos1[2] + pos2[2]);
88 }
89
100 inline void getCenter(const std::array<float, 3> pos0,
101 const std::array<float, 3> pos1,
102 const std::array<float, 3> pos2,
103 const std::array<float, 3> pos3,
104 std::array<float, 3> &incenter) {
105 incenter[0] = 0.25 * (pos0[0] + pos1[0] + pos2[0] + pos3[0]);
106 incenter[1] = 0.25 * (pos0[1] + pos1[1] + pos2[1] + pos3[1]);
107 incenter[2] = 0.25 * (pos0[2] + pos1[2] + pos2[2] + pos3[2]);
108 }
109
119 inline void interpolatePoints(const std::array<float, 3> pos0,
120 const std::array<float, 3> pos1,
121 const float lambda,
122 std::array<float, 3> &result) {
123
124 result[0] = lambda * pos0[0] + (1 - lambda) * pos1[0];
125 result[1] = lambda * pos0[1] + (1 - lambda) * pos1[1];
126 result[2] = lambda * pos0[2] + (1 - lambda) * pos1[2];
127 }
128 } // namespace mth
129
130 class MarchingTetrahedra : public virtual Debug {
131 public:
132 MarchingTetrahedra();
133
135 enum class SURFACE_MODE {
136 SM_SEPARATORS = 0,
137 SM_BOUNDARIES = 1,
138 SM_BOUNDARIES_DETAILED = 2
139 };
140
147 template <typename dataType, typename triangulationType>
148 inline int execute(const dataType *const scalars,
149 const triangulationType &triangulation);
150
163 template <typename triangulationType>
164 int computeMarchingCases_2D(unsigned char *const tetCases,
165 size_t *const numEdges,
166 const unsigned long long *const scalars,
167 const size_t *const triangleCounter,
168 const triangulationType &triangulation) const;
169
181 template <typename triangulationType>
182 int writeSeparators_2D(const unsigned char *const tetCases,
183 const size_t *numEdges,
184 const unsigned long long *const scalars,
185 const triangulationType &triangulation);
186
198 template <typename triangulationType>
199 int writeBoundaries_2D(const unsigned char *const tetCases,
200 const size_t *numEdges,
201 const unsigned long long *const scalars,
202 const triangulationType &triangulation);
203
215 template <typename triangulationType>
216 int writeBoundariesDetailed_2D(const unsigned char *const tetCases,
217 const size_t *numEdges,
218 const unsigned long long *const scalars,
219 const triangulationType &triangulation);
220
221 /* 3D Datasets */
222
235 template <typename triangulationType>
236 int computeMarchingCases_3D(unsigned char *const tetCases,
237 size_t *const numTriangles,
238 const unsigned long long *const scalars,
239 const size_t *const triangleCounter,
240 const triangulationType &triangulation) const;
241
253 template <typename triangulationType>
254 int writeSeparators_3D(const unsigned char *const tetCases,
255 const size_t *numTriangles,
256 const unsigned long long *const scalars,
257 const triangulationType &triangulation);
258
270 template <typename triangulationType>
271 int writeBoundaries_3D(const unsigned char *const tetCases,
272 const size_t *numTriangles,
273 const unsigned long long *const scalars,
274 const triangulationType &triangulation);
275
287 template <typename triangulationType>
288 int writeBoundariesDetailed_3D(const unsigned char *const tetCases,
289 const size_t *numTriangles,
290 const unsigned long long *const scalars,
291 const triangulationType &triangulation);
292
293 protected:
294 // Output options
295 SURFACE_MODE SurfaceMode{SURFACE_MODE::SM_SEPARATORS};
296
297 // Output data
298 SimplexId output_numberOfPoints_{};
299 SimplexId output_numberOfCells_{};
300 std::vector<float> output_points_;
301 std::vector<unsigned long long> output_cells_labels_;
302 std::vector<SimplexId> output_cells_connectivity_;
303 };
304} // namespace ttk
305
306template <typename dataType, typename triangulationType>
307int ttk::MarchingTetrahedra::execute(const dataType *const scalars,
308 const triangulationType &triangulation) {
309
310 Timer t;
311
312 if(scalars == nullptr)
313 return this->printErr("Input scalar field pointer is null.");
314
315 const SimplexId nV = triangulation.getNumberOfVertices();
316 const SimplexId nC = triangulation.getNumberOfCells();
317 const int dim = triangulation.getDimensionality();
318
319 std::vector<unsigned char> tetCases;
320 std::vector<size_t> numberOfTetCases;
321 std::vector<unsigned long long> cScalars;
322
323 cScalars.resize(nV);
324 tetCases.resize(nC);
325
326#ifdef TTK_ENABLE_OPENMP
327 numberOfTetCases.resize(this->threadNumber_);
328#else
329 numberOfTetCases.resize(1);
330#endif // TTK_ENABLE_OPENMP
331
332 for(SimplexId vert = 0; vert < nV; vert++)
333 std::memcpy(&cScalars[vert], &scalars[vert], sizeof(dataType));
334
335 if(dim == 2) {
336 if(SurfaceMode == SURFACE_MODE::SM_SEPARATORS) {
337 computeMarchingCases_2D(&tetCases[0], &numberOfTetCases[0], &cScalars[0],
338 triangleNumberLookup, triangulation);
339 writeSeparators_2D(
340 &tetCases[0], &numberOfTetCases[0], &cScalars[0], triangulation);
341 } else if(SurfaceMode == SURFACE_MODE::SM_BOUNDARIES) {
342 computeMarchingCases_2D(&tetCases[0], &numberOfTetCases[0], &cScalars[0],
343 triangleNumberLookupBoundary, triangulation);
344 writeBoundaries_2D(
345 &tetCases[0], &numberOfTetCases[0], &cScalars[0], triangulation);
346 } else if(SurfaceMode == SURFACE_MODE::SM_BOUNDARIES_DETAILED) {
347 computeMarchingCases_2D(&tetCases[0], &numberOfTetCases[0], &cScalars[0],
348 triangleNumberLookupBoundaryDetailed,
349 triangulation);
350 writeBoundariesDetailed_2D(
351 &tetCases[0], &numberOfTetCases[0], &cScalars[0], triangulation);
352 }
353 } else if(dim == 3) {
354 if(SurfaceMode == SURFACE_MODE::SM_SEPARATORS) {
355 computeMarchingCases_3D(&tetCases[0], &numberOfTetCases[0], &cScalars[0],
356 tetLookupNumWallTriangles, triangulation);
357 writeSeparators_3D(
358 &tetCases[0], &numberOfTetCases[0], &cScalars[0], triangulation);
359 } else if(SurfaceMode == SURFACE_MODE::SM_BOUNDARIES) {
360 computeMarchingCases_3D(&tetCases[0], &numberOfTetCases[0], &cScalars[0],
361 tetLookupNumTrianglesBoundaries, triangulation);
362 writeBoundaries_3D(
363 &tetCases[0], &numberOfTetCases[0], &cScalars[0], triangulation);
364 } else if(SurfaceMode == SURFACE_MODE::SM_BOUNDARIES_DETAILED) {
365 computeMarchingCases_3D(&tetCases[0], &numberOfTetCases[0], &cScalars[0],
366 tetLookupNumTrianglesDetailedBoundary,
367 triangulation);
368 writeBoundariesDetailed_3D(
369 &tetCases[0], &numberOfTetCases[0], &cScalars[0], triangulation);
370 }
371 } else {
372 return this->printErr("Data of dimension " + std::to_string(dim)
373 + "not suported");
374 }
375
376 this->printMsg("Data-set ("
377 + std::to_string(triangulation.getNumberOfVertices())
378 + " points) processed",
379 1.0, t.getElapsedTime(), this->threadNumber_);
380
381 return 0;
382}
383
384template <typename triangulationType>
385int ttk::MarchingTetrahedra::computeMarchingCases_2D(
386 unsigned char *const tetCases,
387 size_t *const numEdges,
388 const unsigned long long *const scalars,
389 const size_t *const triangleCounter,
390 const triangulationType &triangulation) const {
391
392 ttk::Timer localTimer;
393
394 // print the progress of the current subprocedure (currently 0%)
395 this->printMsg("Computing separator cases", 0, 0, this->threadNumber_,
396 ttk::debug::LineMode::REPLACE);
397
398 const SimplexId numTetra = triangulation.getNumberOfCells();
399
400#ifdef TTK_ENABLE_OPENMP
401#pragma omp parallel num_threads(this->threadNumber_)
402 {
403 SimplexId threadEdges = 0;
404 const int tid = omp_get_thread_num();
405
406#pragma omp for schedule(static)
407#else // TTK_ENABLE_OPENMP
408 SimplexId threadEdges = 0;
409 const int tid = 0;
410#endif // TTK_ENABLE_OPENMP
411 for(SimplexId tet = 0; tet < numTetra; ++tet) {
412 std::array<SimplexId, 3> vertices{};
413 triangulation.getCellVertex(tet, 0, vertices[0]);
414 triangulation.getCellVertex(tet, 1, vertices[1]);
415 triangulation.getCellVertex(tet, 2, vertices[2]);
416
417 const std::array<unsigned long long, 3> label
418 = {scalars[vertices[0]], scalars[vertices[1]], scalars[vertices[2]]};
419
420 // Set the third bit to 0 or 1
421 const unsigned char index0 = (label[0] == label[1]) ? 0x00 : 0x04;
422
423 // Set the first and second bit to 0, 1 or 2
424 const unsigned char index1 = (label[0] == label[2]) ? 0x00
425 : (label[1] == label[2]) ? 0x01
426 : 0x02;
427
428 tetCases[tet] = index0 | index1;
429 threadEdges += triangleCounter[tetCases[tet]];
430 }
431 numEdges[tid] = threadEdges;
432#ifdef TTK_ENABLE_OPENMP
433 }
434#endif // TTK_ENABLE_OPENMP
435
436 this->printMsg("Computed separator cases", 1, localTimer.getElapsedTime(),
437 this->threadNumber_);
438
439 return 0;
440}
441
442template <typename triangulationType>
443int ttk::MarchingTetrahedra::writeSeparators_2D(
444 const unsigned char *const tetCases,
445 const size_t *numEdges,
446 const unsigned long long *const scalars,
447 const triangulationType &triangulation) {
448
449 ttk::Timer localTimer;
450
451 this->printMsg("Writing Boundaries", 0, 0, this->threadNumber_,
452 ttk::debug::LineMode::REPLACE);
453
454#ifdef TTK_ENABLE_OPENMP
455 std::vector<size_t> edgeStartIndex(this->threadNumber_ + 1);
456
457 edgeStartIndex[0] = 0;
458
459 // Count triangle number and create iterator start indices
460 for(int t = 1; t <= this->threadNumber_; ++t) {
461 edgeStartIndex[t] = numEdges[t - 1] + edgeStartIndex[t - 1];
462 }
463
464 const size_t numTotalEdges = edgeStartIndex[this->threadNumber_];
465#else // TTK_ENABLE_OPENMP
466 const size_t numTotalEdges = numEdges[0];
467#endif // TTK_ENABLE_OPENMP
468
469 output_points_.resize(6 * numTotalEdges);
470 output_cells_connectivity_.resize(2 * numTotalEdges);
471 output_cells_labels_.resize(numTotalEdges);
472 output_numberOfPoints_ = 2 * numTotalEdges;
473 output_numberOfCells_ = numTotalEdges;
474
475 const SimplexId numTets = triangulation.getNumberOfCells();
476 float *p = output_points_.data();
477 SimplexId *c = output_cells_connectivity_.data();
478 unsigned long long *m = output_cells_labels_.data();
479
480#ifdef TTK_ENABLE_OPENMP
481#pragma omp parallel num_threads(this->threadNumber_) firstprivate(p, c, m)
482 {
483 const int tid = omp_get_thread_num();
484 size_t numThreadIndex = edgeStartIndex[tid];
485
486 p += (numThreadIndex * 6);
487 c += (numThreadIndex * 2);
488 m += numThreadIndex;
489 numThreadIndex = 2 * numThreadIndex;
490
491#pragma omp for schedule(static)
492#else // TTK_ENABLE_OPENMP
493 size_t numThreadIndex = 0;
494#endif // TTK_ENABLE_OPENMP
495 for(SimplexId tet = 0; tet < numTets; ++tet) {
496 if(!triangleLookupIsMultiLabel[tetCases[tet]]) {
497 continue;
498 }
499
500 const int *edgeVerts = triangleLookupEdgeVerts[tetCases[tet]];
501
502 std::array<SimplexId, 3> vertices{};
503 triangulation.getCellVertex(tet, 0, vertices[0]);
504 triangulation.getCellVertex(tet, 1, vertices[1]);
505 triangulation.getCellVertex(tet, 2, vertices[2]);
506
507 std::array<std::array<float, 3>, 3> vertPos{};
508 triangulation.getVertexPoint(
509 vertices[0], vertPos[0][0], vertPos[0][1], vertPos[0][2]);
510 triangulation.getVertexPoint(
511 vertices[1], vertPos[1][0], vertPos[1][1], vertPos[1][2]);
512 triangulation.getVertexPoint(
513 vertices[2], vertPos[2][0], vertPos[2][1], vertPos[2][2]);
514
515 const std::array<unsigned long long, 3> label
516 = {scalars[vertices[0]], scalars[vertices[1]], scalars[vertices[2]]};
517
518 if(triangleLookupIs2Label[tetCases[tet]]) {
519 std::array<std::array<float, 3>, 2> eC{};
520 mth::getCenter(vertPos[edgeVerts[0]], vertPos[edgeVerts[1]], eC[0]);
521 mth::getCenter(vertPos[edgeVerts[2]], vertPos[edgeVerts[3]], eC[1]);
522
523 // Create a hash from vertex label 0 and 1
524 const unsigned long long sparseID
525 = mth::getHash(label[edgeVerts[0]], label[edgeVerts[1]]);
526
527 // Write the edge endpoints, cell indices, and label
528 p[0] = eC[0][0];
529 p[1] = eC[0][1];
530 p[2] = eC[0][2];
531 p[3] = eC[1][0];
532 p[4] = eC[1][1];
533 p[5] = eC[1][2];
534 p += 6;
535
536 c[0] = numThreadIndex;
537 c[1] = numThreadIndex + 1;
538 c += 2;
539 numThreadIndex += 2;
540
541 m[0] = sparseID;
542 m += 1;
543
544 } else {
545 std::array<std::array<float, 3>, 4> eC{};
546 mth::getCenter(vertPos[0], vertPos[1], eC[0]);
547 mth::getCenter(vertPos[0], vertPos[2], eC[1]);
548 mth::getCenter(vertPos[1], vertPos[2], eC[2]);
549 mth::getCenter(vertPos[0], vertPos[1], vertPos[2], eC[3]);
550
551 // Create a hash from all vertex label combinations
552 const std::array<unsigned long long, 3> sparseID
553 = {mth::getHash(label[0], label[1]), mth::getHash(label[0], label[2]),
554 mth::getHash(label[1], label[2])};
555
556 // Write all three lines, each connected to the triangle center
557 p[0] = eC[0][0];
558 p[1] = eC[0][1];
559 p[2] = eC[0][2];
560 p[3] = eC[3][0];
561 p[4] = eC[3][1];
562 p[5] = eC[3][2];
563 p[6] = eC[1][0];
564 p[7] = eC[1][1];
565 p[8] = eC[1][2];
566 p[9] = eC[3][0];
567 p[10] = eC[3][1];
568 p[11] = eC[3][2];
569 p[12] = eC[2][0];
570 p[13] = eC[2][1];
571 p[14] = eC[2][2];
572 p[15] = eC[3][0];
573 p[16] = eC[3][1];
574 p[17] = eC[3][2];
575 p += 18;
576
577 c[0] = numThreadIndex + 0;
578 c[1] = numThreadIndex + 1;
579 c[2] = numThreadIndex + 2;
580 c[3] = numThreadIndex + 3;
581 c[4] = numThreadIndex + 4;
582 c[5] = numThreadIndex + 5;
583 c += 6;
584 numThreadIndex += 6;
585
586 m[0] = sparseID[0];
587 m[1] = sparseID[1];
588 m[2] = sparseID[2];
589 m += 3;
590 }
591 }
592#ifdef TTK_ENABLE_OPENMP
593 }
594#endif // TTK_ENABLE_OPENMP
595
596 this->printMsg(
597 "Wrote Boundaries", 1, localTimer.getElapsedTime(), this->threadNumber_);
598
599 return 0;
600}
601
602template <typename triangulationType>
603int ttk::MarchingTetrahedra::writeBoundaries_2D(
604 const unsigned char *const tetCases,
605 const size_t *numEdges,
606 const unsigned long long *const scalars,
607 const triangulationType &triangulation) {
608
609 ttk::Timer localTimer;
610
611 this->printMsg("Writing Boundaries", 0, 0, this->threadNumber_,
612 ttk::debug::LineMode::REPLACE);
613
614#ifdef TTK_ENABLE_OPENMP
615 std::vector<size_t> edgeStartIndex(this->threadNumber_ + 1);
616
617 edgeStartIndex[0] = 0;
618
619 // Count triangle number and create iterator start indices
620 for(int t = 1; t <= this->threadNumber_; ++t) {
621 edgeStartIndex[t] = numEdges[t - 1] + edgeStartIndex[t - 1];
622 }
623
624 size_t const numTotalEdges = edgeStartIndex[this->threadNumber_];
625#else // TTK_ENABLE_OPENMP
626 size_t numTotalEdges = numEdges[0];
627#endif // TTK_ENABLE_OPENMP
628
629 output_points_.resize(6 * numTotalEdges);
630 output_cells_connectivity_.resize(2 * numTotalEdges);
631 output_cells_labels_.resize(numTotalEdges);
632 output_numberOfPoints_ = 2 * numTotalEdges;
633 output_numberOfCells_ = numTotalEdges;
634
635 float *p = output_points_.data();
636 SimplexId *c = output_cells_connectivity_.data();
637 unsigned long long *m = output_cells_labels_.data();
638
639 const SimplexId numTets = triangulation.getNumberOfCells();
640
641#ifdef TTK_ENABLE_OPENMP
642#pragma omp parallel num_threads(this->threadNumber_) firstprivate(p, c, m)
643 {
644 const int tid = omp_get_thread_num();
645 size_t numThreadIndex = edgeStartIndex[tid];
646
647 p += (numThreadIndex * 6);
648 c += (numThreadIndex * 2);
649 m += numThreadIndex;
650
651 numThreadIndex = 2 * numThreadIndex;
652
653#pragma omp for schedule(static)
654#else // TTK_ENABLE_OPENMP
655 size_t numThreadIndex = 0;
656#endif // TTK_ENABLE_OPENMP
657 for(SimplexId tet = 0; tet < numTets; ++tet) {
658 if(!triangleLookupIsMultiLabel[tetCases[tet]]) {
659 continue;
660 }
661
662 std::array<SimplexId, 3> vertices{};
663 triangulation.getCellVertex(tet, 0, vertices[0]);
664 triangulation.getCellVertex(tet, 1, vertices[1]);
665 triangulation.getCellVertex(tet, 2, vertices[2]);
666
667 std::array<std::array<float, 3>, 3> vertPos{};
668 triangulation.getVertexPoint(
669 vertices[0], vertPos[0][0], vertPos[0][1], vertPos[0][2]);
670 triangulation.getVertexPoint(
671 vertices[1], vertPos[1][0], vertPos[1][1], vertPos[1][2]);
672 triangulation.getVertexPoint(
673 vertices[2], vertPos[2][0], vertPos[2][1], vertPos[2][2]);
674
675 const std::array<unsigned long long, 3> label
676 = {scalars[vertices[0]], scalars[vertices[1]], scalars[vertices[2]]};
677
678 if(triangleLookupIs2Label[tetCases[tet]]) {
679 const int *edgeVerts = triangleLookupEdgeVerts[tetCases[tet]];
680
681 // Write the edge endpoints, cell indices, and label
682 p[0] = vertPos[edgeVerts[0]][0];
683 p[1] = vertPos[edgeVerts[0]][1];
684 p[2] = vertPos[edgeVerts[0]][2];
685 p[3] = vertPos[edgeVerts[2]][0];
686 p[4] = vertPos[edgeVerts[2]][1];
687 p[5] = vertPos[edgeVerts[2]][2];
688 p += 6;
689
690 c[0] = numThreadIndex;
691 c[1] = numThreadIndex + 1;
692 c += 2;
693 numThreadIndex += 2;
694
695 m[0] = label[edgeVerts[0]];
696 m += 1;
697 }
698 }
699#ifdef TTK_ENABLE_OPENMP
700 }
701#endif // TTK_ENABLE_OPENMP
702
703 this->printMsg(
704 "Wrote Boundaries", 1, localTimer.getElapsedTime(), this->threadNumber_);
705
706 return 0;
707}
708
709template <typename triangulationType>
710int ttk::MarchingTetrahedra::writeBoundariesDetailed_2D(
711 const unsigned char *const tetCases,
712 const size_t *numEdges,
713 const unsigned long long *const scalars,
714 const triangulationType &triangulation) {
715
716 ttk::Timer localTimer;
717
718 this->printMsg("Writing Boundaries", 0, 0, this->threadNumber_,
719 ttk::debug::LineMode::REPLACE);
720
721 constexpr float diff = 0.02;
722 constexpr float d0 = 0.5 + diff;
723 constexpr float d1 = 0.5 - diff;
724 constexpr float dc = 2 * diff;
725
726#ifdef TTK_ENABLE_OPENMP
727 std::vector<size_t> edgeStartIndex(this->threadNumber_ + 1);
728
729 edgeStartIndex[0] = 0;
730
731 // Count triangle numbers and create iterator start indices
732 for(int t = 1; t <= this->threadNumber_; ++t) {
733 edgeStartIndex[t] = numEdges[t - 1] + edgeStartIndex[t - 1];
734 }
735
736 size_t const numTotalEdges = edgeStartIndex[this->threadNumber_];
737#else // TTK_ENABLE_OPENMP
738 size_t numTotalEdges = numEdges[0];
739#endif // TTK_ENABLE_OPENMP
740
741 output_points_.resize(6 * numTotalEdges);
742 output_cells_connectivity_.resize(2 * numTotalEdges);
743 output_cells_labels_.resize(numTotalEdges);
744 output_numberOfPoints_ = 2 * numTotalEdges;
745 output_numberOfCells_ = numTotalEdges;
746
747 float *p = output_points_.data();
748 SimplexId *c = output_cells_connectivity_.data();
749 unsigned long long *m = output_cells_labels_.data();
750
751 const SimplexId numTets = triangulation.getNumberOfCells();
752
753#ifdef TTK_ENABLE_OPENMP
754#pragma omp parallel num_threads(this->threadNumber_) firstprivate(p, c, m)
755 {
756 const int tid = omp_get_thread_num();
757 size_t numThreadIndex = edgeStartIndex[tid];
758
759 p += (numThreadIndex * 6);
760 c += (numThreadIndex * 2);
761 m += numThreadIndex;
762
763 numThreadIndex = 2 * numThreadIndex;
764
765#pragma omp for schedule(static)
766#else // TTK_ENABLE_OPENMP
767 size_t numThreadIndex = 0;
768#endif // TTK_ENABLE_OPENMP
769 for(SimplexId tet = 0; tet < numTets; ++tet) {
770 if(!triangleLookupIsMultiLabel[tetCases[tet]]) {
771 continue;
772 }
773
774 const int *vIds = triangleLookupEdgeVerts[tetCases[tet]];
775
776 std::array<SimplexId, 3> vertices{};
777 triangulation.getCellVertex(tet, 0, vertices[0]);
778 triangulation.getCellVertex(tet, 1, vertices[1]);
779 triangulation.getCellVertex(tet, 2, vertices[2]);
780
781 std::array<std::array<float, 3>, 3> vPos{};
782 triangulation.getVertexPoint(
783 vertices[0], vPos[0][0], vPos[0][1], vPos[0][2]);
784 triangulation.getVertexPoint(
785 vertices[1], vPos[1][0], vPos[1][1], vPos[1][2]);
786 triangulation.getVertexPoint(
787 vertices[2], vPos[2][0], vPos[2][1], vPos[2][2]);
788
789 const std::array<unsigned long long, 3> label
790 = {scalars[vertices[0]], scalars[vertices[1]], scalars[vertices[2]]};
791
792 if(triangleLookupIs2Label[tetCases[tet]]) {
793 std::array<float, 3> vert00{}, vert01{}, vert10{}, vert11{};
794
795 mth::interpolatePoints(vPos[vIds[0]], vPos[vIds[1]], d0, vert00);
796 mth::interpolatePoints(vPos[vIds[2]], vPos[vIds[3]], d0, vert01);
797 mth::interpolatePoints(vPos[vIds[0]], vPos[vIds[1]], d1, vert10);
798 mth::interpolatePoints(vPos[vIds[2]], vPos[vIds[3]], d1, vert11);
799
800 // Write the edge endpoints, cell indices, and labels
801 p[0] = vert00[0];
802 p[1] = vert00[1];
803 p[2] = vert00[2];
804 p[3] = vert01[0];
805 p[4] = vert01[1];
806 p[5] = vert01[2];
807 p[6] = vert10[0];
808 p[7] = vert10[1];
809 p[8] = vert10[2];
810 p[9] = vert11[0];
811 p[10] = vert11[1];
812 p[11] = vert11[2];
813 p += 12;
814
815 c[0] = numThreadIndex;
816 c[1] = numThreadIndex + 1;
817 c[2] = numThreadIndex + 2;
818 c[3] = numThreadIndex + 3;
819 c += 4;
820 numThreadIndex += 4;
821
822 m[0] = label[vIds[0]];
823 m[1] = label[vIds[1]];
824 m += 2;
825
826 } else {
827 std::array<float, 3> vert00{}, vert01{}, vert10{}, vert11{}, vert20{},
828 vert21{};
829 mth::interpolatePoints(vPos[0], vPos[1], d0, vert00);
830 mth::interpolatePoints(vPos[0], vPos[2], d0, vert01);
831 mth::interpolatePoints(vPos[1], vPos[0], d0, vert10);
832 mth::interpolatePoints(vPos[1], vPos[2], d0, vert11);
833 mth::interpolatePoints(vPos[2], vPos[0], d0, vert20);
834 mth::interpolatePoints(vPos[2], vPos[1], d0, vert21);
835
836 std::array<float, 3> triCenter{};
837 mth::getCenter(vPos[0], vPos[1], vPos[2], triCenter);
838
839 std::array<float, 3> triS0{}, triS1{}, triS2{};
840 mth::interpolatePoints(vPos[0], triCenter, dc, triS0);
841 mth::interpolatePoints(vPos[1], triCenter, dc, triS1);
842 mth::interpolatePoints(vPos[2], triCenter, dc, triS2);
843
844 // Write the edge endpoints, cell indices, and labels
845 p[0] = vert00[0];
846 p[1] = vert00[1];
847 p[2] = vert00[2];
848 p[3] = triS0[0];
849 p[4] = triS0[1];
850 p[5] = triS0[2];
851 p[6] = triS0[0];
852 p[7] = triS0[1];
853 p[8] = triS0[2];
854 p[9] = vert01[0];
855 p[10] = vert01[1];
856 p[11] = vert01[2];
857
858 p[12] = vert10[0];
859 p[13] = vert10[1];
860 p[14] = vert10[2];
861 p[15] = triS1[0];
862 p[16] = triS1[1];
863 p[17] = triS1[2];
864 p[18] = triS1[0];
865 p[19] = triS1[1];
866 p[20] = triS1[2];
867 p[21] = vert11[0];
868 p[22] = vert11[1];
869 p[23] = vert11[2];
870
871 p[24] = vert20[0];
872 p[25] = vert20[1];
873 p[26] = vert20[2];
874 p[27] = triS2[0];
875 p[28] = triS2[1];
876 p[29] = triS2[2];
877 p[30] = triS2[0];
878 p[31] = triS2[1];
879 p[32] = triS2[2];
880 p[33] = vert21[0];
881 p[34] = vert21[1];
882 p[35] = vert21[2];
883 p += 36;
884
885 c[0] = numThreadIndex + 0;
886 c[1] = numThreadIndex + 1;
887 c[2] = numThreadIndex + 2;
888 c[3] = numThreadIndex + 3;
889
890 c[4] = numThreadIndex + 4;
891 c[5] = numThreadIndex + 5;
892 c[6] = numThreadIndex + 6;
893 c[7] = numThreadIndex + 7;
894
895 c[8] = numThreadIndex + 8;
896 c[9] = numThreadIndex + 9;
897 c[10] = numThreadIndex + 10;
898 c[11] = numThreadIndex + 11;
899 c += 12;
900 numThreadIndex += 12;
901
902 m[0] = label[0];
903 m[1] = label[0];
904 m[2] = label[1];
905 m[3] = label[1];
906 m[4] = label[2];
907 m[5] = label[2];
908 m += 6;
909 }
910 }
911#ifdef TTK_ENABLE_OPENMP
912 }
913#endif // TTK_ENABLE_OPENMP
914
915 this->printMsg(
916 "Wrote Boundaries", 1, localTimer.getElapsedTime(), this->threadNumber_);
917
918 return 0;
919}
920
921template <typename triangulationType>
922int ttk::MarchingTetrahedra::computeMarchingCases_3D(
923 unsigned char *const tetCases,
924 size_t *const numTriangles,
925 const unsigned long long *const scalars,
926 const size_t *const triangleCounter,
927 const triangulationType &triangulation) const {
928
929 ttk::Timer localTimer;
930
931 // print the progress of the current subprocedure (currently 0%)
932 this->printMsg("Computing separator cases", 0, 0, this->threadNumber_,
933 ttk::debug::LineMode::REPLACE);
934
935 const SimplexId numTetra = triangulation.getNumberOfCells();
936
937#ifdef TTK_ENABLE_OPENMP
938#pragma omp parallel num_threads(this->threadNumber_)
939 {
940 SimplexId threadTriangles = 0;
941 const int tid = omp_get_thread_num();
942
943#pragma omp for schedule(static)
944#else
945 SimplexId threadTriangles = 0;
946 const int tid = 0;
947#endif
948 for(SimplexId tet = 0; tet < numTetra; ++tet) {
949
950 std::array<SimplexId, 4> vertices{};
951 triangulation.getCellVertex(tet, 0, vertices[0]);
952 triangulation.getCellVertex(tet, 1, vertices[1]);
953 triangulation.getCellVertex(tet, 2, vertices[2]);
954 triangulation.getCellVertex(tet, 3, vertices[3]);
955
956 const std::array<unsigned long long, 4> label
957 = {scalars[vertices[0]], scalars[vertices[1]], scalars[vertices[2]],
958 scalars[vertices[3]]};
959
960 // Set the fifth bit to 0 or 1
961 const unsigned char index1 = (label[0] == label[1]) ? 0x00 : 0x10;
962
963 // Set the fourth and third bit to 0, 1 or 2
964 const unsigned char index2 = (label[0] == label[2]) ? 0x00
965 : (label[1] == label[2]) ? 0x04
966 : 0x08;
967
968 // Set the first and second bit to 0, 1, 2 or 3
969 const unsigned char index3 = (label[0] == label[3]) ? 0x00
970 : (label[1] == label[3]) ? 0x01
971 : (label[2] == label[3]) ? 0x02
972 : 0x03;
973
974 tetCases[tet] = index1 | index2 | index3;
975 threadTriangles += triangleCounter[tetCases[tet]];
976 }
977 numTriangles[tid] = threadTriangles;
978#ifdef TTK_ENABLE_OPENMP
979 }
980#endif // TTK_ENABLE_OPENMP
981
982 this->printMsg("Computed separator cases", 1, localTimer.getElapsedTime(),
983 this->threadNumber_);
984
985 return 0;
986}
987
988template <typename triangulationType>
989int ttk::MarchingTetrahedra::writeSeparators_3D(
990 const unsigned char *const tetCases,
991 const size_t *numTriangles,
992 const unsigned long long *const scalars,
993 const triangulationType &triangulation) {
994
995 ttk::Timer localTimer;
996
997 this->printMsg("Writing separators", 0, 0, this->threadNumber_,
998 ttk::debug::LineMode::REPLACE);
999
1000#ifdef TTK_ENABLE_OPENMP
1001 std::vector<size_t> triangleStartIndex(this->threadNumber_ + 1);
1002
1003 triangleStartIndex[0] = 0;
1004
1005 // Count triangle number and create iterator start indices
1006 for(int t = 1; t <= this->threadNumber_; ++t) {
1007 triangleStartIndex[t] = numTriangles[t - 1] + triangleStartIndex[t - 1];
1008 }
1009
1010 size_t const numTotalTriangles = triangleStartIndex[this->threadNumber_];
1011#else // TTK_ENABLE_OPENMP
1012 size_t numTotalTriangles = numTriangles[0];
1013#endif // TTK_ENABLE_OPENMP
1014
1015 output_points_.resize(9 * numTotalTriangles);
1016 output_cells_connectivity_.resize(3 * numTotalTriangles);
1017 output_cells_labels_.resize(numTotalTriangles);
1018 output_numberOfPoints_ = 3 * numTotalTriangles;
1019 output_numberOfCells_ = numTotalTriangles;
1020
1021 float *p = output_points_.data();
1022 SimplexId *c = output_cells_connectivity_.data();
1023 unsigned long long *m = output_cells_labels_.data();
1024
1025 const SimplexId numTets = triangulation.getNumberOfCells();
1026
1027#ifdef TTK_ENABLE_OPENMP
1028#pragma omp parallel num_threads(this->threadNumber_) firstprivate(p, c, m)
1029 {
1030 const int tid = omp_get_thread_num();
1031 size_t numThreadIndex = triangleStartIndex[tid];
1032
1033 p += (numThreadIndex * 9);
1034 c += (numThreadIndex * 3);
1035 m += numThreadIndex;
1036
1037 numThreadIndex = 3 * numThreadIndex;
1038
1039#pragma omp for schedule(static)
1040#else // TTK_ENABLE_OPENMP
1041 size_t numThreadIndex = 0;
1042#endif // TTK_ENABLE_OPENMP
1043 for(SimplexId tet = 0; tet < numTets; ++tet) {
1044 if(!tetLookupIsMultiLabel[tetCases[tet]]) {
1045 continue;
1046 }
1047
1048 const int *tetEdgeIndices = tetLookupWall[tetCases[tet]];
1049 const int *tetVertLabel = tetLookupWallLabel[tetCases[tet]];
1050
1051 std::array<SimplexId, 4> vertices{};
1052 triangulation.getCellVertex(tet, 0, vertices[0]);
1053 triangulation.getCellVertex(tet, 1, vertices[1]);
1054 triangulation.getCellVertex(tet, 2, vertices[2]);
1055 triangulation.getCellVertex(tet, 3, vertices[3]);
1056
1057 const std::array<unsigned long long, 4> label
1058 = {scalars[vertices[0]], scalars[vertices[1]], scalars[vertices[2]],
1059 scalars[vertices[3]]};
1060
1061 std::array<std::array<float, 3>, 4> vertPos{};
1062 triangulation.getVertexPoint(
1063 vertices[0], vertPos[0][0], vertPos[0][1], vertPos[0][2]);
1064 triangulation.getVertexPoint(
1065 vertices[1], vertPos[1][0], vertPos[1][1], vertPos[1][2]);
1066 triangulation.getVertexPoint(
1067 vertices[2], vertPos[2][0], vertPos[2][1], vertPos[2][2]);
1068 triangulation.getVertexPoint(
1069 vertices[3], vertPos[3][0], vertPos[3][1], vertPos[3][2]);
1070
1071 std::array<std::array<float, 3>, 10> eC{};
1072 // 6 edge centers
1073 mth::getCenter(vertPos[0], vertPos[1], eC[0]);
1074 mth::getCenter(vertPos[0], vertPos[2], eC[1]);
1075 mth::getCenter(vertPos[0], vertPos[3], eC[2]);
1076 mth::getCenter(vertPos[1], vertPos[2], eC[3]);
1077 mth::getCenter(vertPos[1], vertPos[3], eC[4]);
1078 mth::getCenter(vertPos[2], vertPos[3], eC[5]);
1079
1080 // 4 triangle centers
1081 mth::getCenter(vertPos[0], vertPos[1], vertPos[2], eC[6]);
1082 mth::getCenter(vertPos[0], vertPos[1], vertPos[3], eC[7]);
1083 mth::getCenter(vertPos[0], vertPos[2], vertPos[3], eC[8]);
1084 mth::getCenter(vertPos[1], vertPos[2], vertPos[3], eC[9]);
1085
1086 if(tetEdgeIndices[0] == 10) { // 4 labels on tetraeder
1087 std::array<float, 3> tetCenter{};
1088 mth::getCenter(
1089 vertPos[0], vertPos[1], vertPos[2], vertPos[3], tetCenter);
1090
1091 // Create a hashes from all four label combinations
1092 unsigned long long const sparseMSIds[6] = {
1093 mth::getHash(label[0], label[1]), mth::getHash(label[0], label[2]),
1094 mth::getHash(label[0], label[3]), mth::getHash(label[1], label[2]),
1095 mth::getHash(label[1], label[3]), mth::getHash(label[2], label[3])};
1096
1097 // Write the triangle endpoints, cell indices, and label
1098 p[0] = eC[7][0];
1099 p[1] = eC[7][1];
1100 p[2] = eC[7][2];
1101 p[3] = eC[0][0];
1102 p[4] = eC[0][1];
1103 p[5] = eC[0][2];
1104 p[6] = tetCenter[0];
1105 p[7] = tetCenter[1];
1106 p[8] = tetCenter[2];
1107 p[9] = eC[0][0];
1108 p[10] = eC[0][1];
1109 p[11] = eC[0][2];
1110 p[12] = eC[6][0];
1111 p[13] = eC[6][1];
1112 p[14] = eC[6][2];
1113 p[15] = tetCenter[0];
1114 p[16] = tetCenter[1];
1115 p[17] = tetCenter[2];
1116 p[18] = eC[8][0];
1117 p[19] = eC[8][1];
1118 p[20] = eC[8][2];
1119 p[21] = eC[1][0];
1120 p[22] = eC[1][1];
1121 p[23] = eC[1][2];
1122 p[24] = tetCenter[0];
1123 p[25] = tetCenter[1];
1124 p[26] = tetCenter[2];
1125 p[27] = eC[1][0];
1126 p[28] = eC[1][1];
1127 p[29] = eC[1][2];
1128 p[30] = eC[6][0];
1129 p[31] = eC[6][1];
1130 p[32] = eC[6][2];
1131 p[33] = tetCenter[0];
1132 p[34] = tetCenter[1];
1133 p[35] = tetCenter[2];
1134 p[36] = eC[8][0];
1135 p[37] = eC[8][1];
1136 p[38] = eC[8][2];
1137 p[39] = eC[2][0];
1138 p[40] = eC[2][1];
1139 p[41] = eC[2][2];
1140 p[42] = tetCenter[0];
1141 p[43] = tetCenter[1];
1142 p[44] = tetCenter[2];
1143 p[45] = eC[2][0];
1144 p[46] = eC[2][1];
1145 p[47] = eC[2][2];
1146 p[48] = eC[7][0];
1147 p[49] = eC[7][1];
1148 p[50] = eC[7][2];
1149 p[51] = tetCenter[0];
1150 p[52] = tetCenter[1];
1151 p[53] = tetCenter[2];
1152 p[54] = eC[6][0];
1153 p[55] = eC[6][1];
1154 p[56] = eC[6][2];
1155 p[57] = eC[3][0];
1156 p[58] = eC[3][1];
1157 p[59] = eC[3][2];
1158 p[60] = tetCenter[0];
1159 p[61] = tetCenter[1];
1160 p[62] = tetCenter[2];
1161 p[63] = eC[3][0];
1162 p[64] = eC[3][1];
1163 p[65] = eC[3][2];
1164 p[66] = eC[9][0];
1165 p[67] = eC[9][1];
1166 p[68] = eC[9][2];
1167 p[69] = tetCenter[0];
1168 p[70] = tetCenter[1];
1169 p[71] = tetCenter[2];
1170 p[72] = eC[7][0];
1171 p[73] = eC[7][1];
1172 p[74] = eC[7][2];
1173 p[75] = eC[4][0];
1174 p[76] = eC[4][1];
1175 p[77] = eC[4][2];
1176 p[78] = tetCenter[0];
1177 p[79] = tetCenter[1];
1178 p[80] = tetCenter[2];
1179 p[81] = eC[4][0];
1180 p[82] = eC[4][1];
1181 p[83] = eC[4][2];
1182 p[84] = eC[9][0];
1183 p[85] = eC[9][1];
1184 p[86] = eC[9][2];
1185 p[87] = tetCenter[0];
1186 p[88] = tetCenter[1];
1187 p[89] = tetCenter[2];
1188 p[90] = eC[9][0];
1189 p[91] = eC[9][1];
1190 p[92] = eC[9][2];
1191 p[93] = eC[5][0];
1192 p[94] = eC[5][1];
1193 p[95] = eC[5][2];
1194 p[96] = tetCenter[0];
1195 p[97] = tetCenter[1];
1196 p[98] = tetCenter[2];
1197 p[99] = eC[5][0];
1198 p[100] = eC[5][1];
1199 p[101] = eC[5][2];
1200 p[102] = eC[8][0];
1201 p[103] = eC[8][1];
1202 p[104] = eC[8][2];
1203 p[105] = tetCenter[0];
1204 p[106] = tetCenter[1];
1205 p[107] = tetCenter[2];
1206 p += 108;
1207
1208 c[0] = numThreadIndex + 0;
1209 c[1] = numThreadIndex + 1;
1210 c[2] = numThreadIndex + 2;
1211 c[3] = numThreadIndex + 3;
1212 c[4] = numThreadIndex + 4;
1213 c[5] = numThreadIndex + 5;
1214 c[6] = numThreadIndex + 6;
1215 c[7] = numThreadIndex + 7;
1216 c[8] = numThreadIndex + 8;
1217 c[9] = numThreadIndex + 9;
1218 c[10] = numThreadIndex + 10;
1219 c[11] = numThreadIndex + 11;
1220 c[12] = numThreadIndex + 12;
1221 c[13] = numThreadIndex + 13;
1222 c[14] = numThreadIndex + 14;
1223 c[15] = numThreadIndex + 15;
1224 c[16] = numThreadIndex + 16;
1225 c[17] = numThreadIndex + 17;
1226 c[18] = numThreadIndex + 18;
1227 c[19] = numThreadIndex + 19;
1228 c[20] = numThreadIndex + 20;
1229 c[21] = numThreadIndex + 21;
1230 c[22] = numThreadIndex + 22;
1231 c[23] = numThreadIndex + 23;
1232 c[24] = numThreadIndex + 24;
1233 c[25] = numThreadIndex + 25;
1234 c[26] = numThreadIndex + 26;
1235 c[27] = numThreadIndex + 27;
1236 c[28] = numThreadIndex + 28;
1237 c[29] = numThreadIndex + 29;
1238 c[30] = numThreadIndex + 30;
1239 c[31] = numThreadIndex + 31;
1240 c[32] = numThreadIndex + 32;
1241 c[33] = numThreadIndex + 33;
1242 c[34] = numThreadIndex + 34;
1243 c[35] = numThreadIndex + 35;
1244 c += 36;
1245 numThreadIndex += 36;
1246
1247 m[0] = sparseMSIds[0];
1248 m[1] = sparseMSIds[0];
1249 m[2] = sparseMSIds[1];
1250 m[3] = sparseMSIds[1];
1251 m[4] = sparseMSIds[2];
1252 m[5] = sparseMSIds[2];
1253 m[6] = sparseMSIds[3];
1254 m[7] = sparseMSIds[3];
1255 m[8] = sparseMSIds[4];
1256 m[9] = sparseMSIds[4];
1257 m[10] = sparseMSIds[5];
1258 m[11] = sparseMSIds[5];
1259 m += 12;
1260
1261 } else { // 2 or 3 labels on tetraeder
1262 const size_t numTris = tetLookupNumWallTriangles[tetCases[tet]];
1263
1264 for(size_t t = 0; t < numTris; ++t) {
1265 // Write the triangle endpoints, cell indices, and label
1266 p[0] = eC[tetEdgeIndices[(t * 3)]][0];
1267 p[1] = eC[tetEdgeIndices[(t * 3)]][1];
1268 p[2] = eC[tetEdgeIndices[(t * 3)]][2];
1269 p[3] = eC[tetEdgeIndices[(t * 3) + 1]][0];
1270 p[4] = eC[tetEdgeIndices[(t * 3) + 1]][1];
1271 p[5] = eC[tetEdgeIndices[(t * 3) + 1]][2];
1272 p[6] = eC[tetEdgeIndices[(t * 3) + 2]][0];
1273 p[7] = eC[tetEdgeIndices[(t * 3) + 2]][1];
1274 p[8] = eC[tetEdgeIndices[(t * 3) + 2]][2];
1275 p += 9;
1276
1277 c[0] = numThreadIndex + 0;
1278 c[1] = numThreadIndex + 1;
1279 c[2] = numThreadIndex + 2;
1280 c += 3;
1281 numThreadIndex += 3;
1282
1283 // Create hash from the corresponsing vertex labels
1284 m[0] = mth::getHash(
1285 label[tetVertLabel[t * 2]], label[tetVertLabel[(t * 2) + 1]]);
1286 m += 1;
1287 }
1288 }
1289 }
1290#ifdef TTK_ENABLE_OPENMP
1291 }
1292#endif // TTK_ENABLE_OPENMP
1293
1294 this->printMsg(
1295 "Wrote separators", 1, localTimer.getElapsedTime(), this->threadNumber_);
1296
1297 return 0;
1298}
1299
1300template <typename triangulationType>
1301int ttk::MarchingTetrahedra::writeBoundaries_3D(
1302 const unsigned char *const tetCases,
1303 const size_t *numTriangles,
1304 const unsigned long long *const scalars,
1305 const triangulationType &triangulation) {
1306
1307 ttk::Timer localTimer;
1308
1309 this->printMsg("Writing Boundaries", 0, 0, this->threadNumber_,
1310 ttk::debug::LineMode::REPLACE);
1311
1312#ifdef TTK_ENABLE_OPENMP
1313 std::vector<size_t> triangleStartIndex(this->threadNumber_ + 1);
1314
1315 triangleStartIndex[0] = 0;
1316
1317 // Count triangle number and create iterator start indices
1318 for(int t = 1; t <= this->threadNumber_; ++t) {
1319 triangleStartIndex[t] = numTriangles[t - 1] + triangleStartIndex[t - 1];
1320 }
1321
1322 size_t const numTotalTriangles = triangleStartIndex[this->threadNumber_];
1323#else // TTK_ENABLE_OPENMP
1324 size_t numTotalTriangles = numTriangles[0];
1325#endif // TTK_ENABLE_OPENMP
1326
1327 output_points_.resize(9 * numTotalTriangles);
1328 output_cells_connectivity_.resize(3 * numTotalTriangles);
1329 output_cells_labels_.resize(numTotalTriangles);
1330 output_numberOfPoints_ = 3 * numTotalTriangles;
1331 output_numberOfCells_ = numTotalTriangles;
1332
1333 float *p = output_points_.data();
1334 SimplexId *c = output_cells_connectivity_.data();
1335 unsigned long long *m = output_cells_labels_.data();
1336
1337 const SimplexId numTets = triangulation.getNumberOfCells();
1338
1339#ifdef TTK_ENABLE_OPENMP
1340#pragma omp parallel num_threads(this->threadNumber_) firstprivate(p, c, m)
1341 {
1342 const int tid = omp_get_thread_num();
1343 size_t numThreadIndex = triangleStartIndex[tid];
1344
1345 p += (numThreadIndex * 9);
1346 c += (numThreadIndex * 3);
1347 m += numThreadIndex;
1348
1349 numThreadIndex = 3 * numThreadIndex;
1350
1351#pragma omp for schedule(static)
1352#else // TTK_ENABLE_OPENMP
1353 size_t numThreadIndex = 0;
1354#endif // TTK_ENABLE_OPENMP
1355 for(SimplexId tet = 0; tet < numTets; ++tet) {
1356 if(!tetLookupFast[tetCases[tet]]) {
1357 continue;
1358 }
1359
1360 std::array<SimplexId, 4> vertices{};
1361 triangulation.getCellVertex(tet, 0, vertices[0]);
1362 triangulation.getCellVertex(tet, 1, vertices[1]);
1363 triangulation.getCellVertex(tet, 2, vertices[2]);
1364 triangulation.getCellVertex(tet, 3, vertices[3]);
1365
1366 const std::array<unsigned long long, 4> label
1367 = {scalars[vertices[0]], scalars[vertices[1]], scalars[vertices[2]],
1368 scalars[vertices[3]]};
1369
1370 // Get the tetrahedron local vertex ids that all three have the same
1371 // label
1372 const int id0 = tetLookupFastTri[tetLookupFastCase[tetCases[tet]]][0];
1373 const int id1 = tetLookupFastTri[tetLookupFastCase[tetCases[tet]]][1];
1374 const int id2 = tetLookupFastTri[tetLookupFastCase[tetCases[tet]]][2];
1375
1376 std::array<std::array<float, 3>, 4> vertPos{};
1377 triangulation.getVertexPoint(
1378 vertices[0], vertPos[0][0], vertPos[0][1], vertPos[0][2]);
1379 triangulation.getVertexPoint(
1380 vertices[1], vertPos[1][0], vertPos[1][1], vertPos[1][2]);
1381 triangulation.getVertexPoint(
1382 vertices[2], vertPos[2][0], vertPos[2][1], vertPos[2][2]);
1383 triangulation.getVertexPoint(
1384 vertices[3], vertPos[3][0], vertPos[3][1], vertPos[3][2]);
1385
1386 // Write the triangle endpoints, cell indices, and label
1387 p[0] = vertPos[id0][0];
1388 p[1] = vertPos[id0][1];
1389 p[2] = vertPos[id0][2];
1390 p[3] = vertPos[id1][0];
1391 p[4] = vertPos[id1][1];
1392 p[5] = vertPos[id1][2];
1393 p[6] = vertPos[id2][0];
1394 p[7] = vertPos[id2][1];
1395 p[8] = vertPos[id2][2];
1396 p += 9;
1397
1398 c[0] = numThreadIndex + 0;
1399 c[1] = numThreadIndex + 1;
1400 c[2] = numThreadIndex + 2;
1401 c += 3;
1402 numThreadIndex += 3;
1403
1404 m[0] = label[id0];
1405 m += 1;
1406 }
1407#ifdef TTK_ENABLE_OPENMP
1408 }
1409#endif // TTK_ENABLE_OPENMP
1410
1411 this->printMsg(
1412 "Wrote Boundaries", 1, localTimer.getElapsedTime(), this->threadNumber_);
1413
1414 return 0;
1415}
1416
1417template <typename triangulationType>
1418int ttk::MarchingTetrahedra::writeBoundariesDetailed_3D(
1419 const unsigned char *const tetCases,
1420 const size_t *numTriangles,
1421 const unsigned long long *const scalars,
1422 const triangulationType &triangulation) {
1423
1424 ttk::Timer localTimer;
1425
1426 this->printMsg("Writing detailed boundaries", 0, 0, this->threadNumber_,
1427 ttk::debug::LineMode::REPLACE);
1428
1429 constexpr float diff = 0.02;
1430 constexpr float d0 = 0.5 + diff;
1431 constexpr float d1 = 0.5 - diff;
1432 constexpr float dc = diff * 2;
1433
1434#ifdef TTK_ENABLE_OPENMP
1435 std::vector<size_t> triangleStartIndex(this->threadNumber_ + 1);
1436
1437 triangleStartIndex[0] = 0;
1438
1439 // Count triangle number and create iterator start indices
1440 for(int t = 1; t <= this->threadNumber_; ++t) {
1441 triangleStartIndex[t] = numTriangles[t - 1] + triangleStartIndex[t - 1];
1442 }
1443
1444 size_t const numTotalTriangles = triangleStartIndex[this->threadNumber_];
1445#else // TTK_ENABLE_OPENMP
1446 size_t numTotalTriangles = numTriangles[0];
1447#endif // TTK_ENABLE_OPENMP
1448
1449 output_points_.resize(9 * numTotalTriangles);
1450 output_cells_connectivity_.resize(3 * numTotalTriangles);
1451 output_cells_labels_.resize(numTotalTriangles);
1452 output_numberOfPoints_ = 3 * numTotalTriangles;
1453 output_numberOfCells_ = numTotalTriangles;
1454
1455 float *p = output_points_.data();
1456 SimplexId *c = output_cells_connectivity_.data();
1457 unsigned long long *m = output_cells_labels_.data();
1458
1459 const SimplexId numTets = triangulation.getNumberOfCells();
1460
1461#ifdef TTK_ENABLE_OPENMP
1462#pragma omp parallel num_threads(this->threadNumber_) firstprivate(p, c, m)
1463 {
1464 const int tid = omp_get_thread_num();
1465 size_t numThreadIndex = triangleStartIndex[tid];
1466
1467 p += (numThreadIndex * 9);
1468 c += (numThreadIndex * 3);
1469 m += numThreadIndex;
1470
1471 numThreadIndex = 3 * numThreadIndex;
1472
1473#pragma omp for schedule(static)
1474#else // TTK_ENABLE_OPENMP
1475 size_t numThreadIndex = 0;
1476#endif // TTK_ENABLE_OPENMP
1477 for(SimplexId tet = 0; tet < numTets; ++tet) {
1478 if(!tetLookupIsMultiLabel[tetCases[tet]]) {
1479 continue;
1480 }
1481
1482 std::array<SimplexId, 4> vertices{};
1483 triangulation.getCellVertex(tet, 0, vertices[0]);
1484 triangulation.getCellVertex(tet, 1, vertices[1]);
1485 triangulation.getCellVertex(tet, 2, vertices[2]);
1486 triangulation.getCellVertex(tet, 3, vertices[3]);
1487
1488 const std::array<unsigned long long, 4> label
1489 = {scalars[vertices[0]], scalars[vertices[1]], scalars[vertices[2]],
1490 scalars[vertices[3]]};
1491
1492 std::array<std::array<float, 3>, 4> vPos{};
1493 triangulation.getVertexPoint(
1494 vertices[0], vPos[0][0], vPos[0][1], vPos[0][2]);
1495 triangulation.getVertexPoint(
1496 vertices[1], vPos[1][0], vPos[1][1], vPos[1][2]);
1497 triangulation.getVertexPoint(
1498 vertices[2], vPos[2][0], vPos[2][1], vPos[2][2]);
1499 triangulation.getVertexPoint(
1500 vertices[3], vPos[3][0], vPos[3][1], vPos[3][2]);
1501
1502 if(tetLookupIs2Label[tetCases[tet]]) { // 2 labels (eg. AAAB / AABB)
1503 const int *vIds = tetLookupSplitBasins2Label[tetCases[tet]];
1504
1505 std::array<float, 3> vert00{}, vert01{}, vert02{}, vert10{}, vert11{},
1506 vert12{};
1507
1508 mth::interpolatePoints(vPos[vIds[0]], vPos[vIds[1]], d0, vert00);
1509 mth::interpolatePoints(vPos[vIds[2]], vPos[vIds[3]], d0, vert01);
1510 mth::interpolatePoints(vPos[vIds[4]], vPos[vIds[5]], d0, vert02);
1511 mth::interpolatePoints(vPos[vIds[0]], vPos[vIds[1]], d1, vert10);
1512 mth::interpolatePoints(vPos[vIds[2]], vPos[vIds[3]], d1, vert11);
1513 mth::interpolatePoints(vPos[vIds[4]], vPos[vIds[5]], d1, vert12);
1514
1515 // Write the triangle endpoints, cell indices, and label
1516 p[0] = vert00[0];
1517 p[1] = vert00[1];
1518 p[2] = vert00[2];
1519 p[3] = vert01[0];
1520 p[4] = vert01[1];
1521 p[5] = vert01[2];
1522 p[6] = vert02[0];
1523 p[7] = vert02[1];
1524 p[8] = vert02[2];
1525
1526 p[9] = vert10[0];
1527 p[10] = vert10[1];
1528 p[11] = vert10[2];
1529 p[12] = vert11[0];
1530 p[13] = vert11[1];
1531 p[14] = vert11[2];
1532 p[15] = vert12[0];
1533 p[16] = vert12[1];
1534 p[17] = vert12[2];
1535 p += 18;
1536
1537 c[0] = numThreadIndex + 0;
1538 c[1] = numThreadIndex + 1;
1539 c[2] = numThreadIndex + 2;
1540 c[3] = numThreadIndex + 3;
1541 c[4] = numThreadIndex + 4;
1542 c[5] = numThreadIndex + 5;
1543 c += 6;
1544 numThreadIndex += 6;
1545
1546 m[0] = label[vIds[0]];
1547 m[1] = label[vIds[1]];
1548 m += 2;
1549
1550 if(vIds[6] != -1) { // 2 vertices per label (e.g. AABB)
1551 std::array<float, 3> vert03{}, vert13{};
1552
1553 mth::interpolatePoints(vPos[vIds[6]], vPos[vIds[7]], d0, vert03);
1554 mth::interpolatePoints(vPos[vIds[6]], vPos[vIds[7]], d1, vert13);
1555
1556 // Write the triangle endpoints, cell indices, and label
1557 p[0] = vert00[0];
1558 p[1] = vert00[1];
1559 p[2] = vert00[2];
1560 p[3] = vert01[0];
1561 p[4] = vert01[1];
1562 p[5] = vert01[2];
1563 p[6] = vert03[0];
1564 p[7] = vert03[1];
1565 p[8] = vert03[2];
1566 p[9] = vert10[0];
1567 p[10] = vert10[1];
1568 p[11] = vert10[2];
1569 p[12] = vert11[0];
1570 p[13] = vert11[1];
1571 p[14] = vert11[2];
1572 p[15] = vert13[0];
1573 p[16] = vert13[1];
1574 p[17] = vert13[2];
1575 p += 18;
1576
1577 c[0] = numThreadIndex + 0;
1578 c[1] = numThreadIndex + 1;
1579 c[2] = numThreadIndex + 2;
1580 c[3] = numThreadIndex + 3;
1581 c[4] = numThreadIndex + 4;
1582 c[5] = numThreadIndex + 5;
1583 c += 6;
1584 numThreadIndex += 6;
1585
1586 m[0] = label[vIds[0]];
1587 m[1] = label[vIds[1]];
1588 m += 2;
1589 }
1590 } else if(tetLookupIs3Label[tetCases[tet]]) {
1591 const int *vIds = tetLookupSplitBasisns3Label[tetCases[tet]];
1592
1593 std::array<std::array<float, 3>, 4> triCenter{};
1594 mth::getCenter(vPos[0], vPos[1], vPos[2], triCenter[0]);
1595 mth::getCenter(vPos[0], vPos[1], vPos[3], triCenter[1]);
1596 mth::getCenter(vPos[0], vPos[2], vPos[3], triCenter[2]);
1597 mth::getCenter(vPos[1], vPos[2], vPos[3], triCenter[3]);
1598
1599 std::array<std::array<float, 3>, 10> edgeCenters{};
1600
1601 mth::getCenter(vPos[0], vPos[1], edgeCenters[0]);
1602 mth::getCenter(vPos[0], vPos[2], edgeCenters[1]);
1603 mth::getCenter(vPos[0], vPos[3], edgeCenters[2]);
1604 mth::getCenter(vPos[1], vPos[2], edgeCenters[3]);
1605 mth::getCenter(vPos[1], vPos[3], edgeCenters[4]);
1606 mth::getCenter(vPos[2], vPos[3], edgeCenters[5]);
1607
1608 std::array<float, 3> edge00{}, edge01{}, edge02{}, edge03{}, tri00{},
1609 tri01{}, edge10{}, edge11{}, edge12{}, tri10{}, tri11{}, edge20{},
1610 edge21{}, edge22{}, tri20{}, tri21{};
1611
1612 mth::interpolatePoints(vPos[vIds[0]], edgeCenters[vIds[1]], dc, edge00);
1613 mth::interpolatePoints(vPos[vIds[0]], edgeCenters[vIds[2]], dc, edge01);
1614 mth::interpolatePoints(vPos[vIds[0]], triCenter[vIds[3]], dc, tri00);
1615 mth::interpolatePoints(vPos[vIds[4]], edgeCenters[vIds[5]], dc, edge02);
1616 mth::interpolatePoints(vPos[vIds[4]], edgeCenters[vIds[6]], dc, edge03);
1617 mth::interpolatePoints(vPos[vIds[4]], triCenter[vIds[7]], dc, tri01);
1618
1619 mth::interpolatePoints(vPos[vIds[8]], edgeCenters[vIds[1]], dc, edge10);
1620 mth::interpolatePoints(vPos[vIds[8]], edgeCenters[vIds[5]], dc, edge11);
1621 mth::interpolatePoints(
1622 vPos[vIds[8]], edgeCenters[vIds[10]], dc, edge12);
1623 mth::interpolatePoints(vPos[vIds[8]], triCenter[vIds[3]], dc, tri10);
1624 mth::interpolatePoints(vPos[vIds[8]], triCenter[vIds[7]], dc, tri11);
1625
1626 mth::interpolatePoints(vPos[vIds[9]], edgeCenters[vIds[2]], dc, edge20);
1627 mth::interpolatePoints(vPos[vIds[9]], edgeCenters[vIds[6]], dc, edge21);
1628 mth::interpolatePoints(
1629 vPos[vIds[9]], edgeCenters[vIds[10]], dc, edge22);
1630 mth::interpolatePoints(vPos[vIds[9]], triCenter[vIds[3]], dc, tri20);
1631 mth::interpolatePoints(vPos[vIds[9]], triCenter[vIds[7]], dc, tri21);
1632
1633 // Write the triangle endpoints, cell indices, and labels
1634 // Label 0
1635 p[0] = edge00[0];
1636 p[1] = edge00[1];
1637 p[2] = edge00[2];
1638 p[3] = edge02[0];
1639 p[4] = edge02[1];
1640 p[5] = edge02[2];
1641 p[6] = tri00[0];
1642 p[7] = tri00[1];
1643 p[8] = tri00[2];
1644 p[9] = edge02[0];
1645 p[10] = edge02[1];
1646 p[11] = edge02[2];
1647 p[12] = tri00[0];
1648 p[13] = tri00[1];
1649 p[14] = tri00[2];
1650 p[15] = tri01[0];
1651 p[16] = tri01[1];
1652 p[17] = tri01[2];
1653 p[18] = edge01[0];
1654 p[19] = edge01[1];
1655 p[20] = edge01[2];
1656 p[21] = edge03[0];
1657 p[22] = edge03[1];
1658 p[23] = edge03[2];
1659 p[24] = tri00[0];
1660 p[25] = tri00[1];
1661 p[26] = tri00[2];
1662 p[27] = edge03[0];
1663 p[28] = edge03[1];
1664 p[29] = edge03[2];
1665 p[30] = tri00[0];
1666 p[31] = tri00[1];
1667 p[32] = tri00[2];
1668 p[33] = tri01[0];
1669 p[34] = tri01[1];
1670 p[35] = tri01[2];
1671
1672 // Label 1
1673 p[36] = edge10[0];
1674 p[37] = edge10[1];
1675 p[38] = edge10[2];
1676 p[39] = edge11[0];
1677 p[40] = edge11[1];
1678 p[41] = edge11[2];
1679 p[42] = tri10[0];
1680 p[43] = tri10[1];
1681 p[44] = tri10[2];
1682 p[45] = edge11[0];
1683 p[46] = edge11[1];
1684 p[47] = edge11[2];
1685 p[48] = tri10[0];
1686 p[49] = tri10[1];
1687 p[50] = tri10[2];
1688 p[51] = tri11[0];
1689 p[52] = tri11[1];
1690 p[53] = tri11[2];
1691 p[54] = edge12[0];
1692 p[55] = edge12[1];
1693 p[56] = edge12[2];
1694 p[57] = tri10[0];
1695 p[58] = tri10[1];
1696 p[59] = tri10[2];
1697 p[60] = tri11[0];
1698 p[61] = tri11[1];
1699 p[62] = tri11[2];
1700
1701 // Label 2
1702 p[63] = edge20[0];
1703 p[64] = edge20[1];
1704 p[65] = edge20[2];
1705 p[66] = edge21[0];
1706 p[67] = edge21[1];
1707 p[68] = edge21[2];
1708 p[69] = tri20[0];
1709 p[70] = tri20[1];
1710 p[71] = tri20[2];
1711 p[72] = edge21[0];
1712 p[73] = edge21[1];
1713 p[74] = edge21[2];
1714 p[75] = tri20[0];
1715 p[76] = tri20[1];
1716 p[77] = tri20[2];
1717 p[78] = tri21[0];
1718 p[79] = tri21[1];
1719 p[80] = tri21[2];
1720 p[81] = edge22[0];
1721 p[82] = edge22[1];
1722 p[83] = edge22[2];
1723 p[84] = tri20[0];
1724 p[85] = tri20[1];
1725 p[86] = tri20[2];
1726 p[87] = tri21[0];
1727 p[88] = tri21[1];
1728 p[89] = tri21[2];
1729 p += 90;
1730
1731 c[0] = numThreadIndex + 0;
1732 c[1] = numThreadIndex + 1;
1733 c[2] = numThreadIndex + 2;
1734 c[3] = numThreadIndex + 3;
1735 c[4] = numThreadIndex + 4;
1736 c[5] = numThreadIndex + 5;
1737 c[6] = numThreadIndex + 6;
1738 c[7] = numThreadIndex + 7;
1739 c[8] = numThreadIndex + 8;
1740 c[9] = numThreadIndex + 9;
1741 c[10] = numThreadIndex + 10;
1742 c[11] = numThreadIndex + 11;
1743 c[12] = numThreadIndex + 12;
1744 c[13] = numThreadIndex + 13;
1745 c[14] = numThreadIndex + 14;
1746 c[15] = numThreadIndex + 15;
1747 c[16] = numThreadIndex + 16;
1748 c[17] = numThreadIndex + 17;
1749 c[18] = numThreadIndex + 18;
1750 c[19] = numThreadIndex + 19;
1751 c[20] = numThreadIndex + 20;
1752 c[21] = numThreadIndex + 21;
1753 c[22] = numThreadIndex + 22;
1754 c[23] = numThreadIndex + 23;
1755 c[24] = numThreadIndex + 24;
1756 c[25] = numThreadIndex + 25;
1757 c[26] = numThreadIndex + 26;
1758 c[27] = numThreadIndex + 27;
1759 c[28] = numThreadIndex + 28;
1760 c[29] = numThreadIndex + 29;
1761 c += 30;
1762 numThreadIndex += 30;
1763
1764 m[0] = label[vIds[0]];
1765 m[1] = label[vIds[0]];
1766 m[2] = label[vIds[0]];
1767 m[3] = label[vIds[0]];
1768 m[4] = label[vIds[8]];
1769 m[5] = label[vIds[8]];
1770 m[6] = label[vIds[8]];
1771 m[7] = label[vIds[9]];
1772 m[8] = label[vIds[9]];
1773 m[9] = label[vIds[9]];
1774 m += 10;
1775
1776 } else { // 4 labels
1777 std::array<float, 3> tetCenter{};
1778 mth::getCenter(vPos[0], vPos[1], vPos[2], vPos[3], tetCenter);
1779
1780 // the 4 triangle centers
1781 std::array<std::array<float, 3>, 4> triCenter{};
1782 mth::getCenter(vPos[0], vPos[1], vPos[2], triCenter[0]);
1783 mth::getCenter(vPos[0], vPos[1], vPos[3], triCenter[1]);
1784 mth::getCenter(vPos[0], vPos[2], vPos[3], triCenter[2]);
1785 mth::getCenter(vPos[1], vPos[2], vPos[3], triCenter[3]);
1786
1787 std::array<float, 3> vert00{}, vert01{}, vert02{}, vert0tet{},
1788 vert0t0{}, vert0t1{}, vert0t2{}, vert10{}, vert11{}, vert12{},
1789 vert1tet{}, vert1t0{}, vert1t1{}, vert1t2{}, vert20{}, vert21{},
1790 vert22{}, vert2tet{}, vert2t0{}, vert2t1{}, vert2t2{}, vert30{},
1791 vert31{}, vert32{}, vert3tet{}, vert3t0{}, vert3t1{}, vert3t2{};
1792
1793 mth::interpolatePoints(vPos[0], vPos[1], d0, vert00);
1794 mth::interpolatePoints(vPos[0], vPos[2], d0, vert01);
1795 mth::interpolatePoints(vPos[0], vPos[3], d0, vert02);
1796 mth::interpolatePoints(vPos[0], tetCenter, dc, vert0tet);
1797 mth::interpolatePoints(vPos[0], triCenter[0], dc, vert0t0);
1798 mth::interpolatePoints(vPos[0], triCenter[1], dc, vert0t1);
1799 mth::interpolatePoints(vPos[0], triCenter[2], dc, vert0t2);
1800
1801 mth::interpolatePoints(vPos[1], vPos[0], d0, vert10);
1802 mth::interpolatePoints(vPos[1], vPos[2], d0, vert11);
1803 mth::interpolatePoints(vPos[1], vPos[3], d0, vert12);
1804 mth::interpolatePoints(vPos[1], tetCenter, dc, vert1tet);
1805 mth::interpolatePoints(vPos[1], triCenter[0], dc, vert1t0);
1806 mth::interpolatePoints(vPos[1], triCenter[1], dc, vert1t1);
1807 mth::interpolatePoints(vPos[1], triCenter[3], dc, vert1t2);
1808
1809 mth::interpolatePoints(vPos[2], vPos[0], d0, vert20);
1810 mth::interpolatePoints(vPos[2], vPos[1], d0, vert21);
1811 mth::interpolatePoints(vPos[2], vPos[3], d0, vert22);
1812 mth::interpolatePoints(vPos[2], tetCenter, dc, vert2tet);
1813 mth::interpolatePoints(vPos[2], triCenter[0], dc, vert2t0);
1814 mth::interpolatePoints(vPos[2], triCenter[2], dc, vert2t1);
1815 mth::interpolatePoints(vPos[2], triCenter[3], dc, vert2t2);
1816
1817 mth::interpolatePoints(vPos[3], vPos[0], d0, vert30);
1818 mth::interpolatePoints(vPos[3], vPos[1], d0, vert31);
1819 mth::interpolatePoints(vPos[3], vPos[2], d0, vert32);
1820 mth::interpolatePoints(vPos[3], tetCenter, dc, vert3tet);
1821 mth::interpolatePoints(vPos[3], triCenter[1], dc, vert3t0);
1822 mth::interpolatePoints(vPos[3], triCenter[2], dc, vert3t1);
1823 mth::interpolatePoints(vPos[3], triCenter[3], dc, vert3t2);
1824
1825 // Write the triangle endpoints, cell indices, and labels
1826 // Label Vert 0
1827 p[0] = vert00[0];
1828 p[1] = vert00[1];
1829 p[2] = vert00[2];
1830 p[3] = vert0t0[0];
1831 p[4] = vert0t0[1];
1832 p[5] = vert0t0[2];
1833 p[6] = vert0tet[0];
1834 p[7] = vert0tet[1];
1835 p[8] = vert0tet[2];
1836 p[9] = vert00[0];
1837 p[10] = vert00[1];
1838 p[11] = vert00[2];
1839 p[12] = vert0t1[0];
1840 p[13] = vert0t1[1];
1841 p[14] = vert0t1[2];
1842 p[15] = vert0tet[0];
1843 p[16] = vert0tet[1];
1844 p[17] = vert0tet[2];
1845 p[18] = vert01[0];
1846 p[19] = vert01[1];
1847 p[20] = vert01[2];
1848 p[21] = vert0t0[0];
1849 p[22] = vert0t0[1];
1850 p[23] = vert0t0[2];
1851 p[24] = vert0tet[0];
1852 p[25] = vert0tet[1];
1853 p[26] = vert0tet[2];
1854 p[27] = vert01[0];
1855 p[28] = vert01[1];
1856 p[29] = vert01[2];
1857 p[30] = vert0t2[0];
1858 p[31] = vert0t2[1];
1859 p[32] = vert0t2[2];
1860 p[33] = vert0tet[0];
1861 p[34] = vert0tet[1];
1862 p[35] = vert0tet[2];
1863 p[36] = vert02[0];
1864 p[37] = vert02[1];
1865 p[38] = vert02[2];
1866 p[39] = vert0t2[0];
1867 p[40] = vert0t2[1];
1868 p[41] = vert0t2[2];
1869 p[42] = vert0tet[0];
1870 p[43] = vert0tet[1];
1871 p[44] = vert0tet[2];
1872 p[45] = vert02[0];
1873 p[46] = vert02[1];
1874 p[47] = vert02[2];
1875 p[48] = vert0t1[0];
1876 p[49] = vert0t1[1];
1877 p[50] = vert0t1[2];
1878 p[51] = vert0tet[0];
1879 p[52] = vert0tet[1];
1880 p[53] = vert0tet[2];
1881
1882 // Label Vert 1
1883 p[54] = vert10[0];
1884 p[55] = vert10[1];
1885 p[56] = vert10[2];
1886 p[57] = vert1t0[0];
1887 p[58] = vert1t0[1];
1888 p[59] = vert1t0[2];
1889 p[60] = vert1tet[0];
1890 p[61] = vert1tet[1];
1891 p[62] = vert1tet[2];
1892 p[63] = vert10[0];
1893 p[64] = vert10[1];
1894 p[65] = vert10[2];
1895 p[66] = vert1t1[0];
1896 p[67] = vert1t1[1];
1897 p[68] = vert1t1[2];
1898 p[69] = vert1tet[0];
1899 p[70] = vert1tet[1];
1900 p[71] = vert1tet[2];
1901 p[72] = vert11[0];
1902 p[73] = vert11[1];
1903 p[74] = vert11[2];
1904 p[75] = vert1t0[0];
1905 p[76] = vert1t0[1];
1906 p[77] = vert1t0[2];
1907 p[78] = vert1tet[0];
1908 p[79] = vert1tet[1];
1909 p[80] = vert1tet[2];
1910 p[81] = vert11[0];
1911 p[82] = vert11[1];
1912 p[83] = vert11[2];
1913 p[84] = vert1t2[0];
1914 p[85] = vert1t2[1];
1915 p[86] = vert1t2[2];
1916 p[87] = vert1tet[0];
1917 p[88] = vert1tet[1];
1918 p[89] = vert1tet[2];
1919 p[90] = vert12[0];
1920 p[91] = vert12[1];
1921 p[92] = vert12[2];
1922 p[93] = vert1t2[0];
1923 p[94] = vert1t2[1];
1924 p[95] = vert1t2[2];
1925 p[96] = vert1tet[0];
1926 p[97] = vert1tet[1];
1927 p[98] = vert1tet[2];
1928 p[99] = vert12[0];
1929 p[100] = vert12[1];
1930 p[101] = vert12[2];
1931 p[102] = vert1t1[0];
1932 p[103] = vert1t1[1];
1933 p[104] = vert1t1[2];
1934 p[105] = vert1tet[0];
1935 p[106] = vert1tet[1];
1936 p[107] = vert1tet[2];
1937
1938 // Label Vert 2
1939 p[108] = vert20[0];
1940 p[109] = vert20[1];
1941 p[110] = vert20[2];
1942 p[111] = vert2t0[0];
1943 p[112] = vert2t0[1];
1944 p[113] = vert2t0[2];
1945 p[114] = vert2tet[0];
1946 p[115] = vert2tet[1];
1947 p[116] = vert2tet[2];
1948 p[117] = vert20[0];
1949 p[118] = vert20[1];
1950 p[119] = vert20[2];
1951 p[120] = vert2t1[0];
1952 p[121] = vert2t1[1];
1953 p[122] = vert2t1[2];
1954 p[123] = vert2tet[0];
1955 p[124] = vert2tet[1];
1956 p[125] = vert2tet[2];
1957 p[126] = vert21[0];
1958 p[127] = vert21[1];
1959 p[128] = vert21[2];
1960 p[129] = vert2t0[0];
1961 p[130] = vert2t0[1];
1962 p[131] = vert2t0[2];
1963 p[132] = vert2tet[0];
1964 p[133] = vert2tet[1];
1965 p[134] = vert2tet[2];
1966 p[135] = vert21[0];
1967 p[136] = vert21[1];
1968 p[137] = vert21[2];
1969 p[138] = vert2t2[0];
1970 p[139] = vert2t2[1];
1971 p[140] = vert2t2[2];
1972 p[141] = vert2tet[0];
1973 p[142] = vert2tet[1];
1974 p[143] = vert2tet[2];
1975 p[144] = vert22[0];
1976 p[145] = vert22[1];
1977 p[146] = vert22[2];
1978 p[147] = vert2t2[0];
1979 p[148] = vert2t2[1];
1980 p[149] = vert2t2[2];
1981 p[150] = vert2tet[0];
1982 p[151] = vert2tet[1];
1983 p[152] = vert2tet[2];
1984 p[153] = vert22[0];
1985 p[154] = vert22[1];
1986 p[155] = vert22[2];
1987 p[156] = vert2t1[0];
1988 p[157] = vert2t1[1];
1989 p[158] = vert2t1[2];
1990 p[159] = vert2tet[0];
1991 p[160] = vert2tet[1];
1992 p[161] = vert2tet[2];
1993
1994 // Label Vert 3
1995 p[162] = vert30[0];
1996 p[163] = vert30[1];
1997 p[164] = vert30[2];
1998 p[165] = vert3t0[0];
1999 p[166] = vert3t0[1];
2000 p[167] = vert3t0[2];
2001 p[168] = vert3tet[0];
2002 p[169] = vert3tet[1];
2003 p[170] = vert3tet[2];
2004 p[171] = vert30[0];
2005 p[172] = vert30[1];
2006 p[173] = vert30[2];
2007 p[174] = vert3t1[0];
2008 p[175] = vert3t1[1];
2009 p[176] = vert3t1[2];
2010 p[177] = vert3tet[0];
2011 p[178] = vert3tet[1];
2012 p[179] = vert3tet[2];
2013 p[180] = vert31[0];
2014 p[181] = vert31[1];
2015 p[182] = vert31[2];
2016 p[183] = vert3t0[0];
2017 p[184] = vert3t0[1];
2018 p[185] = vert3t0[2];
2019 p[186] = vert3tet[0];
2020 p[187] = vert3tet[1];
2021 p[188] = vert3tet[2];
2022 p[189] = vert31[0];
2023 p[190] = vert31[1];
2024 p[191] = vert31[2];
2025 p[192] = vert3t2[0];
2026 p[193] = vert3t2[1];
2027 p[194] = vert3t2[2];
2028 p[195] = vert3tet[0];
2029 p[196] = vert3tet[1];
2030 p[197] = vert3tet[2];
2031 p[198] = vert32[0];
2032 p[199] = vert32[1];
2033 p[200] = vert32[2];
2034 p[201] = vert3t2[0];
2035 p[202] = vert3t2[1];
2036 p[203] = vert3t2[2];
2037 p[204] = vert3tet[0];
2038 p[205] = vert3tet[1];
2039 p[206] = vert3tet[2];
2040 p[207] = vert32[0];
2041 p[208] = vert32[1];
2042 p[209] = vert32[2];
2043 p[210] = vert3t1[0];
2044 p[211] = vert3t1[1];
2045 p[212] = vert3t1[2];
2046 p[213] = vert3tet[0];
2047 p[214] = vert3tet[1];
2048 p[215] = vert3tet[2];
2049 p += 216;
2050
2051 c[0] = numThreadIndex + 0;
2052 c[1] = numThreadIndex + 1;
2053 c[2] = numThreadIndex + 2;
2054 c[3] = numThreadIndex + 3;
2055 c[4] = numThreadIndex + 4;
2056 c[5] = numThreadIndex + 5;
2057 c[6] = numThreadIndex + 6;
2058 c[7] = numThreadIndex + 7;
2059 c[8] = numThreadIndex + 8;
2060 c[9] = numThreadIndex + 9;
2061 c[10] = numThreadIndex + 10;
2062 c[11] = numThreadIndex + 11;
2063 c[12] = numThreadIndex + 12;
2064 c[13] = numThreadIndex + 13;
2065 c[14] = numThreadIndex + 14;
2066 c[15] = numThreadIndex + 15;
2067 c[16] = numThreadIndex + 16;
2068 c[17] = numThreadIndex + 17;
2069 c[18] = numThreadIndex + 18;
2070 c[19] = numThreadIndex + 19;
2071 c[20] = numThreadIndex + 20;
2072 c[21] = numThreadIndex + 21;
2073 c[22] = numThreadIndex + 22;
2074 c[23] = numThreadIndex + 23;
2075 c[24] = numThreadIndex + 24;
2076 c[25] = numThreadIndex + 25;
2077 c[26] = numThreadIndex + 26;
2078 c[27] = numThreadIndex + 27;
2079 c[28] = numThreadIndex + 28;
2080 c[29] = numThreadIndex + 29;
2081 c[30] = numThreadIndex + 30;
2082 c[31] = numThreadIndex + 31;
2083 c[32] = numThreadIndex + 32;
2084 c[33] = numThreadIndex + 33;
2085 c[34] = numThreadIndex + 34;
2086 c[35] = numThreadIndex + 35;
2087 c[36] = numThreadIndex + 36;
2088 c[37] = numThreadIndex + 37;
2089 c[38] = numThreadIndex + 38;
2090 c[39] = numThreadIndex + 39;
2091 c[40] = numThreadIndex + 40;
2092 c[41] = numThreadIndex + 41;
2093 c[42] = numThreadIndex + 42;
2094 c[43] = numThreadIndex + 43;
2095 c[44] = numThreadIndex + 44;
2096 c[45] = numThreadIndex + 45;
2097 c[46] = numThreadIndex + 46;
2098 c[47] = numThreadIndex + 47;
2099 c[48] = numThreadIndex + 48;
2100 c[49] = numThreadIndex + 49;
2101 c[50] = numThreadIndex + 50;
2102 c[51] = numThreadIndex + 51;
2103 c[52] = numThreadIndex + 52;
2104 c[53] = numThreadIndex + 53;
2105 c[54] = numThreadIndex + 54;
2106 c[55] = numThreadIndex + 55;
2107 c[56] = numThreadIndex + 56;
2108 c[57] = numThreadIndex + 57;
2109 c[58] = numThreadIndex + 58;
2110 c[59] = numThreadIndex + 59;
2111 c[60] = numThreadIndex + 60;
2112 c[61] = numThreadIndex + 61;
2113 c[62] = numThreadIndex + 62;
2114 c[63] = numThreadIndex + 63;
2115 c[64] = numThreadIndex + 64;
2116 c[65] = numThreadIndex + 65;
2117 c[66] = numThreadIndex + 66;
2118 c[67] = numThreadIndex + 67;
2119 c[68] = numThreadIndex + 68;
2120 c[69] = numThreadIndex + 69;
2121 c[70] = numThreadIndex + 70;
2122 c[71] = numThreadIndex + 71;
2123 c += 72;
2124 numThreadIndex += 72;
2125
2126 m[0] = label[0];
2127 m[1] = label[0];
2128 m[2] = label[0];
2129 m[3] = label[0];
2130 m[4] = label[0];
2131 m[5] = label[0];
2132 m[6] = label[1];
2133 m[7] = label[1];
2134 m[8] = label[1];
2135 m[9] = label[1];
2136 m[10] = label[1];
2137 m[11] = label[1];
2138 m[12] = label[2];
2139 m[13] = label[2];
2140 m[14] = label[2];
2141 m[15] = label[2];
2142 m[16] = label[2];
2143 m[17] = label[2];
2144 m[18] = label[3];
2145 m[19] = label[3];
2146 m[20] = label[3];
2147 m[21] = label[3];
2148 m[22] = label[3];
2149 m[23] = label[3];
2150 m += 24;
2151 }
2152 }
2153#ifdef TTK_ENABLE_OPENMP
2154 }
2155#endif // TTK_ENABLE_OPENMP
2156
2157 this->printMsg("Wrote detailed boundaries", 1, localTimer.getElapsedTime(),
2158 this->threadNumber_);
2159
2160 return 0;
2161}
triangleLookupIs2Label
constexpr bool triangleLookupIs2Label[7]
Definition MarchingTetrahedraLookupTables.inl:450
tetLookupSplitBasisns3Label
constexpr int tetLookupSplitBasisns3Label[27][11]
Definition MarchingTetrahedraLookupTables.inl:408
triangleLookupIsMultiLabel
constexpr bool triangleLookupIsMultiLabel[7]
Definition MarchingTetrahedraLookupTables.inl:439
tetLookupIsMultiLabel
constexpr bool tetLookupIsMultiLabel[28]
Definition MarchingTetrahedraLookupTables.inl:194
tetLookupIs2Label
constexpr bool tetLookupIs2Label[28]
Definition MarchingTetrahedraLookupTables.inl:226
tetLookupFastTri
constexpr int tetLookupFastTri[4][3]
Definition MarchingTetrahedraLookupTables.inl:363
tetLookupWallLabel
constexpr int tetLookupWallLabel[27][10]
Definition MarchingTetrahedraLookupTables.inl:56
triangleNumberLookup
constexpr size_t triangleNumberLookup[7]
Definition MarchingTetrahedraLookupTables.inl:461
triangleLookupEdgeVerts
constexpr int triangleLookupEdgeVerts[7][4]
Definition MarchingTetrahedraLookupTables.inl:497
tetLookupNumTrianglesDetailedBoundary
constexpr size_t tetLookupNumTrianglesDetailedBoundary[28]
Definition MarchingTetrahedraLookupTables.inl:125
tetLookupNumTrianglesBoundaries
constexpr size_t tetLookupNumTrianglesBoundaries[28]
Definition MarchingTetrahedraLookupTables.inl:162
triangleNumberLookupBoundary
constexpr size_t triangleNumberLookupBoundary[7]
Definition MarchingTetrahedraLookupTables.inl:472
tetLookupSplitBasins2Label
constexpr int tetLookupSplitBasins2Label[22][8]
Definition MarchingTetrahedraLookupTables.inl:374
tetLookupIs3Label
constexpr bool tetLookupIs3Label[28]
Definition MarchingTetrahedraLookupTables.inl:258
triangleNumberLookupBoundaryDetailed
constexpr size_t triangleNumberLookupBoundaryDetailed[7]
Definition MarchingTetrahedraLookupTables.inl:483
tetLookupNumWallTriangles
constexpr size_t tetLookupNumWallTriangles[28]
Definition MarchingTetrahedraLookupTables.inl:89
tetLookupFastCase
constexpr int tetLookupFastCase[28]
Definition MarchingTetrahedraLookupTables.inl:293
tetLookupFast
constexpr bool tetLookupFast[28]
Definition MarchingTetrahedraLookupTables.inl:327
ttk::Debug
Minimalist debugging class.
Definition Debug.h:88
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::MarchingTetrahedra
TTK processing package for Marching Tetra/Triangles computations.
Definition MarchingTetrahedra.h:130
ttk::MarchingTetrahedra::writeBoundaries_3D
int writeBoundaries_3D(const unsigned char *const tetCases, const size_t *numTriangles, const unsigned long long *const scalars, const triangulationType &triangulation)
Definition MarchingTetrahedra.h:1301
ttk::MarchingTetrahedra::writeBoundariesDetailed_3D
int writeBoundariesDetailed_3D(const unsigned char *const tetCases, const size_t *numTriangles, const unsigned long long *const scalars, const triangulationType &triangulation)
Definition MarchingTetrahedra.h:1418
ttk::MarchingTetrahedra::computeMarchingCases_2D
int computeMarchingCases_2D(unsigned char *const tetCases, size_t *const numEdges, const unsigned long long *const scalars, const size_t *const triangleCounter, const triangulationType &triangulation) const
Definition MarchingTetrahedra.h:385
ttk::MarchingTetrahedra::writeBoundaries_2D
int writeBoundaries_2D(const unsigned char *const tetCases, const size_t *numEdges, const unsigned long long *const scalars, const triangulationType &triangulation)
Definition MarchingTetrahedra.h:603
ttk::MarchingTetrahedra::output_points_
std::vector< float > output_points_
Definition MarchingTetrahedra.h:300
ttk::MarchingTetrahedra::SURFACE_MODE
SURFACE_MODE
Type of 2-separatrix output.
Definition MarchingTetrahedra.h:135
ttk::MarchingTetrahedra::output_cells_labels_
std::vector< unsigned long long > output_cells_labels_
Definition MarchingTetrahedra.h:301
ttk::MarchingTetrahedra::computeMarchingCases_3D
int computeMarchingCases_3D(unsigned char *const tetCases, size_t *const numTriangles, const unsigned long long *const scalars, const size_t *const triangleCounter, const triangulationType &triangulation) const
Definition MarchingTetrahedra.h:922
ttk::MarchingTetrahedra::output_cells_connectivity_
std::vector< SimplexId > output_cells_connectivity_
Definition MarchingTetrahedra.h:302
ttk::MarchingTetrahedra::writeSeparators_2D
int writeSeparators_2D(const unsigned char *const tetCases, const size_t *numEdges, const unsigned long long *const scalars, const triangulationType &triangulation)
Definition MarchingTetrahedra.h:443
ttk::MarchingTetrahedra::writeSeparators_3D
int writeSeparators_3D(const unsigned char *const tetCases, const size_t *numTriangles, const unsigned long long *const scalars, const triangulationType &triangulation)
Definition MarchingTetrahedra.h:989
ttk::MarchingTetrahedra::writeBoundariesDetailed_2D
int writeBoundariesDetailed_2D(const unsigned char *const tetCases, const size_t *numEdges, const unsigned long long *const scalars, const triangulationType &triangulation)
Definition MarchingTetrahedra.h:710
ttk::MarchingTetrahedra::execute
int execute(const dataType *const scalars, const triangulationType &triangulation)
Definition MarchingTetrahedra.h:307
ttk::Timer::getElapsedTime
double getElapsedTime()
Definition Timer.h:15
tetLookupWall
constexpr int tetLookupWall[28][15]
Lookuptables used in MarchingTetrahedra algorithm.
Definition MarchingTetrahedraLookupTables.inl:19
std::to_string
std::string to_string(__int128)
Definition ripserpy.cpp:99
ttk::mth::getCenter
void getCenter(const std::array< float, 3 > pos0, const std::array< float, 3 > pos1, std::array< float, 3 > &incenter)
Get the center of two points.
Definition MarchingTetrahedra.h:64
ttk::mth::getHash
constexpr unsigned long long int getHash(const unsigned long long int a, const unsigned long long int b)
Get a hash value from two keys.
Definition MarchingTetrahedra.h:50
ttk::mth::interpolatePoints
void interpolatePoints(const std::array< float, 3 > pos0, const std::array< float, 3 > pos1, const float lambda, std::array< float, 3 > &result)
Interpolate between two points (lambda = 0 -> pos1 / 1 -> pos0)
Definition MarchingTetrahedra.h:119
ttk
The Topology ToolKit.
Definition AbstractTriangulation.h:51
ttk::SimplexId
int SimplexId
Identifier type for simplices of any dimension.
Definition DataTypes.h:22
printMsg
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)
Generated on Mon May 6 2024 07:44:29 for TTK by doxygen 1.9.8