Geogram Version 1.8.5
A programming library of geometric algorithms
Loading...
Searching...
No Matches
generic_RVD_cell.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2000-2022 Inria
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 * * Redistributions in binary form must reproduce the above copyright notice,
11 * this list of conditions and the following disclaimer in the documentation
12 * and/or other materials provided with the distribution.
13 * * Neither the name of the ALICE Project-Team nor the names of its
14 * contributors may be used to endorse or promote products derived from this
15 * software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 * POSSIBILITY OF SUCH DAMAGE.
28 *
29 * Contact: Bruno Levy
30 *
31 * https://www.inria.fr/fr/bruno-levy
32 *
33 * Inria,
34 * Domaine de Voluceau,
35 * 78150 Le Chesnay - Rocquencourt
36 * FRANCE
37 *
38 */
39
40#ifndef GEOGRAM_VORONOI_GENERIC_RVD_CELL
41#define GEOGRAM_VORONOI_GENERIC_RVD_CELL
42
47#include <iosfwd>
48#include <stack>
49
59namespace GEOGen {
60
61 using GEO::Mesh;
62
69 class GEOGRAM_API ConvexCell {
70
72 typedef ConvexCell thisclass;
73
74 public:
75
76 static const index_t NO_TRIANGLE = index_t(-1);
77 static const index_t NO_VERTEX = index_t(-1);
78 static const index_t END_OF_LIST = index_t(-1);
79
84 TRI_IS_FREE = 0,
85 TRI_IS_CONFLICT = 1,
86 TRI_IS_USED = 2
87 };
88
99 struct Triangle {
100
111 index_t v0, index_t v1, index_t v2,
112 index_t f0, index_t f1, index_t f2
113 ) :
114 next_(END_OF_LIST),
115 status_(TRI_IS_FREE),
116 id_(-1) {
117 v[0] = v0;
118 v[1] = v1;
119 v[2] = v2;
120 t[0] = f0;
121 t[1] = f1;
122 t[2] = f2;
123 }
124
129 next_(END_OF_LIST),
130 status_(TRI_IS_FREE),
131 id_(-1) {
132 v[0] = NO_VERTEX;
133 v[1] = NO_VERTEX;
134 v[2] = NO_VERTEX;
135 t[0] = NO_TRIANGLE;
136 t[1] = NO_TRIANGLE;
137 t[2] = NO_TRIANGLE;
138 }
139
140 GEOGen::Vertex dual_;
141 index_t v[3]; // The 3 vertices of this triangle
142 index_t t[3]; // The 3 triangles adjacent to this triangle
143 index_t next_; // Linked list management.
144 TriangleStatus status_;
145 // TRI_IS_FREE,TRI_IS_USED or TRI_IS_CONFLICT.
146 signed_index_t id_;
147 };
148
155 class Vertex {
156 public:
161 t(-1),
162 id_(-1) {
163 }
164
165 signed_index_t t; // One triangle incident to this vertex
166 signed_index_t id_;
167 };
168
173 ConvexCell(coord_index_t dim) :
174 first_free_(END_OF_LIST),
175 v_to_t_dirty_(false),
176 intersections_(dim),
177 symbolic_is_surface_(false),
178 cell_id_(-1) {
179 }
180
188 void copy(const ConvexCell& rhs);
189
194 coord_index_t dimension() const {
195 return intersections_.dimension();
196 }
197
201 void clear() {
202 first_free_ = END_OF_LIST;
203 triangles_.resize(0);
204 vertices_.resize(0);
205 v_to_t_dirty_ = false;
206 intersections_.clear();
207 }
208
218 symbolic_is_surface_ = x;
219 }
220
232 const Mesh* mesh, index_t t, bool symbolic,
233 const GEO::Attribute<double>& vertex_weight
234 );
235
236
245 Mesh* mesh, bool symbolic
246 );
247
248
262 void convert_to_mesh(Mesh* mesh, bool copy_symbolic_info = false);
263
281 template <index_t DIM>
282 signed_index_t clip_by_plane(
283 const Mesh* mesh, const Delaunay* delaunay,
284 index_t i, index_t j,
285 bool exact, bool symbolic
286 ) {
287 index_t new_v = create_vertex();
288 set_vertex_id(index_t(new_v), signed_index_t(j)+1);
289 index_t conflict_begin, conflict_end;
290
291 // Phase I: Determine the conflict zone and chain the triangles
292 // Note: they are not immediately deleted, since we need the
293 // geometric information in the triangles to compute the new
294 // intersections.
295 get_conflict_list<DIM>(
296 mesh, delaunay, i, j, exact, conflict_begin, conflict_end
297 );
298
299 // Special case: the clipping plane did not clip anything.
300 if(conflict_begin == END_OF_LIST) {
301 return signed_index_t(new_v);
302 }
303
304 // Phase II: Find a triangle on the border of the conflict zone
305 // (by traversing the conflict list).
306 index_t first_conflict_t;
307 index_t first_conflict_e;
308 bool found_h = find_triangle_on_border(
309 conflict_begin, conflict_end,
310 first_conflict_t, first_conflict_e
311 );
312
313 // The clipping plane removed everything !
314 // (note: cannot be empty conflict list, since this
315 // case was detected by previous test at the end of
316 // phase I)
317 if(!found_h) {
318 clear();
319 return -1;
320 }
321
322 // Phase III: Triangulate hole.
323 triangulate_hole<DIM>(
324 delaunay, i, j, symbolic,
325 first_conflict_t, first_conflict_e,
326 new_v
327 );
328
329 // Phase IV: Merge the conflict zone into the free list.
330 merge_into_free_list(conflict_begin, conflict_end);
331 return signed_index_t(new_v);
332 }
333
339 index_t max_t() const {
340 return triangles_.size();
341 }
342
346 index_t nb_t() const {
347 index_t result = 0;
348 for(index_t t = 0; t < max_t(); t++) {
349 if(triangle_is_used(t)) {
350 result++;
351 }
352 }
353 return result;
354 }
355
359 index_t max_v() const {
360 return vertices_.size();
361 }
362
366 bool triangle_is_free(index_t t) const {
367 geo_debug_assert(t != NO_TRIANGLE);
368 geo_debug_assert(t < max_t());
369 return triangles_[t].status_ == TRI_IS_FREE;
370 }
371
380 bool triangle_is_valid(index_t t) const {
381 return !triangle_is_free(t);
382 }
383
389 bool triangle_is_used(index_t t) const {
390 geo_debug_assert(t != NO_TRIANGLE);
391 geo_debug_assert(t < max_t());
392 return triangles_[t].status_ == TRI_IS_USED;
393 }
394
401 bool triangle_is_conflict(index_t t) const {
402 geo_debug_assert(t != NO_TRIANGLE);
403 geo_debug_assert(t < max_t());
404 return triangles_[t].status_ == TRI_IS_CONFLICT;
405 }
406
414 index_t triangle_vertex(index_t t, index_t iv) const {
415 geo_debug_assert(iv < 3);
416 geo_debug_assert(triangle_is_valid(t));
417 return triangles_[t].v[iv];
418 }
419
426 index_t triangle_adjacent(index_t t, index_t e) const {
427 geo_debug_assert(e < 3);
428 geo_debug_assert(triangle_is_valid(t));
429 return triangles_[t].t[e];
430 }
431
438 void set_triangle_vertex(index_t t, index_t iv, index_t v) {
439 geo_debug_assert(t != NO_TRIANGLE);
440 geo_debug_assert(t < max_t());
441 geo_debug_assert(iv < 3);
442 geo_debug_assert(v < max_v());
443 triangles_[t].v[iv] = v;
444 }
445
453 index_t t, index_t e, index_t t2
454 ) {
455 geo_debug_assert(t != NO_TRIANGLE);
456 geo_debug_assert(t < max_t());
457 geo_debug_assert(e < 3);
458 geo_debug_assert(t2 < max_t());
459 triangles_[t].t[e] = t2;
460 }
461
469 index_t find_triangle_vertex(index_t t, index_t v) const {
470 geo_debug_assert(t != NO_TRIANGLE);
471 geo_debug_assert(t < max_t());
472 geo_debug_assert(v < max_v());
473
474 // The following expression is 10% faster than using
475 // if() statements (multiply by boolean result of test).
476 // Thank to Laurent Alonso for this idea.
477 index_t result = index_t(
478 (triangles_[t].v[1] == v) | ((triangles_[t].v[2] == v) * 2)
479 );
480 geo_debug_assert(triangles_[t].v[result] == v);
481 return result;
482 }
483
493 index_t t1, index_t t2
494 ) const {
495 geo_debug_assert(t1 != NO_TRIANGLE);
496 geo_debug_assert(t1 < max_t());
497 geo_debug_assert(t2 != NO_TRIANGLE);
498 geo_debug_assert(t2 < max_t());
499
500 // The following expression is 10% faster than using
501 // if() statements (multiply by boolean result of test).
502 // Thank to Laurent Alonso for this idea.
503 index_t result = index_t(
504 (triangles_[t1].t[1] == t2) | ((triangles_[t1].t[2] == t2) * 2)
505 );
506
507 geo_debug_assert(triangles_[t1].t[result] == t2);
508
509 return result;
510 }
511
519 signed_index_t vertex_triangle(index_t v) const {
520 if(v_to_t_dirty_) {
521 const_cast<ConvexCell*>(this)->init_v_to_t();
522 }
523 geo_debug_assert(v != NO_VERTEX);
524 geo_debug_assert(v < max_v());
525 return vertices_[v].t;
526 }
527
534 void set_vertex_triangle(index_t v, index_t t) {
535 geo_debug_assert(v != NO_VERTEX);
536 geo_debug_assert(index_t(v) < max_v());
537 vertices_[v].t = signed_index_t(t);
538 }
539
549 const GEOGen::Vertex& triangle_dual(index_t t) const {
550 geo_debug_assert(triangle_is_valid(t));
551 return triangles_[t].dual_;
552 }
553
564 geo_debug_assert(triangle_is_valid(t));
565 return triangles_[t].dual_;
566 }
567
575 signed_index_t cell_id() const {
576 return cell_id_;
577 }
578
586 void set_cell_id(signed_index_t i) {
587 cell_id_ = i;
588 }
589
596 signed_index_t triangle_id(index_t t) const {
597 geo_debug_assert(triangle_is_valid(t));
598 return triangles_[t].id_;
599 }
600
606 void set_triangle_id(index_t t, signed_index_t id) {
607 geo_debug_assert(triangle_is_valid(t));
608 triangles_[t].id_ = id;
609 }
610
616 signed_index_t vertex_id(index_t v) const {
617 geo_debug_assert(v != NO_VERTEX);
618 geo_debug_assert(v < max_v());
619 return vertices_[v].id_;
620 }
621
627 void set_vertex_id(index_t v, signed_index_t id) {
628 geo_debug_assert(v != NO_VERTEX);
629 geo_debug_assert(v < max_v());
630 vertices_[v].id_ = id;
631 }
632
638 class Corner {
639 public:
644 t(NO_TRIANGLE),
645 v(3) {
646 }
647
654 Corner(index_t t_in, index_t v_in) :
655 t(t_in),
656 v(v_in) {
657 }
658
666 bool operator== (const Corner& rhs) const {
667 return t == rhs.t && v == rhs.v;
668 }
669
677 bool operator!= (const Corner& rhs) const {
678 return t != rhs.t || v != rhs.v;
679 }
680
681 index_t t;
682 index_t v;
683 };
684
691 index_t t2 = triangle_adjacent(c.t, plus1mod3(c.v));
692 index_t v = triangle_vertex(c.t, c.v);
693 c.v = find_triangle_vertex(t2, v);
694 c.t = t2;
695 }
696
701 void init_v_to_t() {
702 v_to_t_dirty_ = false;
703 for(index_t v = 0; v < max_v(); v++) {
704 set_vertex_triangle(v, NO_TRIANGLE);
705 }
706 for(index_t t = 0; t < max_t(); t++) {
707 if(triangle_is_used(t)) {
708 for(index_t iv = 0; iv < 3; iv++) {
709 set_vertex_triangle(triangle_vertex(t, iv), t);
710 }
711 }
712 }
713 }
714
720 index_t create_triangle() {
721 if(first_free_ == END_OF_LIST) {
722 grow();
723 }
724 index_t result = first_free_;
725 first_free_ = next_triangle(first_free_);
726 mark_as_used(result);
727 set_triangle_id(result, -1);
728 return result;
729 }
730
740 index_t create_triangle(index_t v0, index_t v1, index_t v2) {
741 index_t t = create_triangle();
742 triangles_[t].v[0] = v0;
743 triangles_[t].v[1] = v1;
744 triangles_[t].v[2] = v2;
745 return t;
746 }
747
761 index_t v0, index_t v1, index_t v2,
762 index_t t0, index_t t1, index_t t2
763 ) {
764 index_t t = create_triangle();
765 triangles_[t].v[0] = v0;
766 triangles_[t].v[1] = v1;
767 triangles_[t].v[2] = v2;
768 triangles_[t].t[0] = t0;
769 triangles_[t].t[1] = t1;
770 triangles_[t].t[2] = t2;
771 return t;
772 }
773
785 index_t t,
786 index_t v0, index_t v1, index_t v2,
787 index_t t0, index_t t1, index_t t2
788 ) {
789 geo_debug_assert(t != NO_TRIANGLE);
790 geo_debug_assert(t < max_t());
791 triangles_[t].v[0] = v0;
792 triangles_[t].v[1] = v1;
793 triangles_[t].v[2] = v2;
794 triangles_[t].t[0] = t0;
795 triangles_[t].t[1] = t1;
796 triangles_[t].t[2] = t2;
797 }
798
817 const double* p,
818 double w,
819 index_t v0, index_t v1, index_t v2,
820 index_t t0, index_t t1, index_t t2
821 ) {
822 index_t t = create_triangle(v0, v1, v2, t0, t1, t2);
823 triangle_dual(t).set_point(p);
824 triangle_dual(t).set_weight(w);
825 return t;
826 }
827
845 const double* p,
846 index_t v0, index_t v1, index_t v2,
847 index_t t0, index_t t1, index_t t2
848 ) {
849 index_t t = create_triangle(v0, v1, v2, t0, t1, t2);
850 double* np = intersections_.new_item();
851 for(coord_index_t c = 0; c < dimension(); ++c) {
852 np[c] = p[c];
853 }
854 triangle_dual(t).set_point(np);
855 return t;
856 }
857
862 index_t create_vertex() {
863 v_to_t_dirty_ = true;
864 vertices_.push_back(Vertex());
865 return vertices_.size() - 1;
866 }
867
886 template <index_t DIM>
888 const Delaunay* delaunay,
889 index_t i, index_t j, bool symbolic,
890 index_t t1, index_t t1ebord,
891 index_t v_in
892 ) {
893 index_t t = t1;
894 index_t e = t1ebord;
895 index_t t_adj = triangle_adjacent(t,e);
896 geo_debug_assert(t_adj != index_t(-1));
897
898 geo_debug_assert(triangle_is_conflict(t));
899 geo_debug_assert(!triangle_is_conflict(t_adj));
900
901 index_t new_t_first = index_t(-1);
902 index_t new_t_prev = index_t(-1);
903
904 do {
905
906 index_t v1 = triangle_vertex(t, plus1mod3(e));
907 index_t v2 = triangle_vertex(t, minus1mod3(e));
908
909 // Create new triangle
910 index_t new_t = create_triangle(v_in, v1, v2);
911
912 triangle_dual(new_t).intersect_geom<DIM>(
913 intersections_,
914 triangle_dual(t),
915 triangle_dual(triangle_adjacent(t, e)),
916 delaunay->vertex_ptr(i), delaunay->vertex_ptr(j)
917 );
918
919 if(symbolic) {
920 triangle_dual(new_t).sym().intersect_symbolic(
921 triangle_dual(t).sym(),
922 triangle_dual(triangle_adjacent(t, e)).sym(),
923 j
924 );
925 }
926
927 // Connect new triangle to triangle on the other
928 // side of the conflict zone.
929 set_triangle_adjacent(new_t, 0, t_adj);
930 index_t adj_e = triangle_adjacent_index(t_adj, t);
931 set_triangle_adjacent(t_adj, adj_e, new_t);
932
933
934 // Move to next triangle
935 e = plus1mod3(e);
936 t_adj = index_t(triangle_adjacent(t,e));
937 while(triangle_is_conflict(t_adj)) {
938 t = t_adj;
939 e = minus1mod3(find_triangle_vertex(t,v2));
940 t_adj = index_t(triangle_adjacent(t,e));
941 geo_debug_assert(t_adj != index_t(-1));
942 }
943
944 if(new_t_prev == index_t(-1)) {
945 new_t_first = new_t;
946 } else {
947 set_triangle_adjacent(new_t_prev, 1, new_t);
948 set_triangle_adjacent(new_t, 2, new_t_prev);
949 }
950
951 new_t_prev = new_t;
952
953 } while((t != t1) || (e != t1ebord));
954
955 // Connect last triangle to first triangle
956 set_triangle_adjacent(new_t_prev, 1, new_t_first);
957 set_triangle_adjacent(new_t_first, 2, new_t_prev);
958
959 return new_t_prev;
960 }
961
979 template <index_t DIM>
981 const Mesh* mesh, const Delaunay* delaunay,
982 index_t i, index_t j, bool exact,
983 index_t& conflict_begin, index_t& conflict_end
984 ) {
985 conflict_begin = END_OF_LIST;
986 conflict_end = END_OF_LIST;
987 if(exact) {
988 // In exact mode, we classify each vertex
989 // using the exact predicate. Note that
990 // "climbing/walking" from a random vertex
991 // would be more efficient, but it would require a
992 // "comparison" exact predicate (with 8 different
993 // versions according to the configuration of the
994 // two vertices to be compared)
995 for(index_t t = 0; t < max_t(); t++) {
996 if(triangle_is_used(t)) {
997 Sign s = side<DIM>(
998 mesh, delaunay,
999 triangle_dual(t),
1000 i, j,
1001 exact
1002 );
1003 if(s == GEO::NEGATIVE) {
1004 append_triangle_to_conflict_list(
1005 t, conflict_begin, conflict_end
1006 );
1007 }
1008 }
1009 }
1010 } else {
1011 // In non-exact mode, we first detect the
1012 // vertex that is furthest away along the
1013 // normal vector of the clipping bisector,
1014 // and then we propagate using a flood-fill
1015 // algorithm. This strategy ensures that
1016 // the conflict zone remains connected, even
1017 // in the presence of numerical errors. Clearly
1018 // it does not ensure validity, but in practice
1019 // it improves resistance to degeneracies.
1020 index_t furthest_t =
1021 find_furthest_point_linear_scan<DIM>(
1022 delaunay, i, j
1023 );
1024 propagate_conflict_list<DIM>(
1025 mesh, delaunay, furthest_t,
1026 i, j, exact,
1027 conflict_begin, conflict_end
1028 );
1029 }
1030 }
1031
1043 template <index_t DIM>
1045 const Delaunay* delaunay, index_t i, index_t j
1046 ) const {
1047 index_t result = NO_TRIANGLE;
1048 double furthest_dist = 0.0;
1049 for(index_t t = 0; t < max_t(); ++t) {
1050 if(triangle_is_used(t)) {
1051 double d = signed_bisector_distance<DIM>(
1052 delaunay, i, j, triangle_dual(t).point()
1053 );
1054 if(d < furthest_dist) {
1055 result = t;
1056 furthest_dist = d;
1057 }
1058 }
1059 }
1060 return (furthest_dist < 0) ? result : NO_TRIANGLE;
1061 }
1062
1074 template <index_t DIM>
1076 const Delaunay* delaunay, index_t i, index_t j, const double* q
1077 ) {
1078 const double* pi = delaunay->vertex_ptr(i);
1079 const double* pj = delaunay->vertex_ptr(j);
1080 double result = 0;
1081 for(coord_index_t c = 0; c < DIM; ++c) {
1082 result += GEO::geo_sqr(q[c] - pj[c]);
1083 result -= GEO::geo_sqr(q[c] - pi[c]);
1084 }
1085 return result;
1086 }
1087
1104 template <index_t DIM>
1106 const Mesh* mesh, const Delaunay* delaunay,
1107 index_t first_t,
1108 index_t i, index_t j, bool exact,
1109 index_t& conflict_begin, index_t& conflict_end
1110 ) {
1111 conflict_begin = END_OF_LIST;
1112 conflict_end = END_OF_LIST;
1113
1114 // Special case, clipping plane does not clip anything
1115 if(first_t == NO_TRIANGLE) {
1116 return;
1117 }
1118
1119 std::stack<index_t> S;
1120 S.push(first_t);
1121 append_triangle_to_conflict_list(
1122 first_t, conflict_begin, conflict_end
1123 );
1124 while(!S.empty()) {
1125 index_t t = S.top();
1126 S.pop();
1127 for(unsigned int e = 0; e < 3; e++) {
1128 index_t neigh = index_t(triangle_adjacent(t, e));
1129 if(!triangle_is_conflict(neigh)) {
1130 if(
1131 side<DIM>(
1132 mesh, delaunay, triangle_dual(neigh),
1133 i, j, exact
1134 ) == GEO::NEGATIVE
1135 ) {
1136 S.push(neigh);
1137 append_triangle_to_conflict_list(
1138 neigh, conflict_begin, conflict_end
1139 );
1140 }
1141 }
1142 }
1143 }
1144 }
1145
1164 template <index_t DIM>
1165 Sign side(
1166 const Mesh* mesh, const Delaunay* delaunay,
1167 const GEOGen::Vertex& v,
1168 index_t i, index_t j, bool exact
1169 ) const {
1170 Sign result = GEO::ZERO;
1171 if(exact) {
1172 result = side_exact(
1173 mesh, delaunay, v,
1174 delaunay->vertex_ptr(i),
1175 delaunay->vertex_ptr(j),
1176 DIM,
1177 symbolic_is_surface_
1178 );
1179 } else {
1180 result = v.side_fast<DIM>(
1181 delaunay->vertex_ptr(i),
1182 delaunay->vertex_ptr(j)
1183 );
1184 }
1185 return result;
1186 }
1187
1207 const Mesh* mesh, const Delaunay* delaunay,
1208 const GEOGen::Vertex& v,
1209 const double* pi, const double* pj,
1210 coord_index_t dim,
1211 bool symbolic_is_surface = false
1212 ) const;
1213
1227 index_t conflict_begin, index_t conflict_end,
1228 index_t& t, index_t& e
1229 ) const {
1230 GEO::geo_argused(conflict_end);
1231 t = conflict_begin;
1232 do {
1233 for(e = 0; e < 3; ++e) {
1234 index_t nt = triangle_adjacent(t, e);
1235 if(triangle_is_used(nt)) {
1236 return true;
1237 }
1238 }
1239 t = next_triangle(t);
1240 } while(t != END_OF_LIST);
1241 return false;
1242 }
1243
1249 index_t next_triangle(index_t t) const {
1250 geo_debug_assert(t != NO_TRIANGLE);
1251 geo_debug_assert(t < max_t());
1252 return triangles_[t].next_;
1253 }
1254
1262 void set_next_triangle(index_t t, index_t t2) {
1263 geo_debug_assert(t != NO_TRIANGLE);
1264 geo_debug_assert(t < max_t());
1265 triangles_[t].next_ = t2;
1266 }
1267
1273 void mark_as_free(index_t t) {
1274 geo_debug_assert(t != NO_TRIANGLE);
1275 geo_debug_assert(t < max_t());
1276 triangles_[t].status_ = TRI_IS_FREE;
1277 }
1278
1282 void mark_as_conflict(index_t t) {
1283 geo_debug_assert(t != NO_TRIANGLE);
1284 geo_debug_assert(t < max_t());
1285 triangles_[t].status_ = TRI_IS_CONFLICT;
1286 }
1287
1291 void mark_as_used(index_t t) {
1292 geo_debug_assert(t != NO_TRIANGLE);
1293 geo_debug_assert(t < max_t());
1294 triangles_[t].status_ = TRI_IS_USED;
1295 }
1296
1307 index_t t, index_t& conflict_begin, index_t& conflict_end
1308 ) {
1309 geo_debug_assert(triangle_is_used(t));
1310 set_next_triangle(t, conflict_begin);
1311 mark_as_conflict(t);
1312 conflict_begin = t;
1313 if(conflict_end == END_OF_LIST) {
1314 conflict_end = t;
1315 }
1316 }
1317
1324 void merge_into_free_list(index_t list_begin, index_t list_end) {
1325 if(list_begin != END_OF_LIST) {
1326 geo_debug_assert(list_end != END_OF_LIST);
1327
1328 index_t cur = list_begin;
1329 while(cur != list_end) {
1330 mark_as_free(cur);
1331 cur = next_triangle(cur);
1332 }
1333 mark_as_free(list_end);
1334 set_next_triangle(list_end, first_free_);
1335 first_free_ = list_begin;
1336 }
1337 }
1338
1344 void grow() {
1345 geo_debug_assert(first_free_ == END_OF_LIST);
1346 first_free_ = triangles_.size();
1347 triangles_.push_back(Triangle());
1348 }
1349
1354 std::ostream& show_stats(std::ostream& os) const;
1355
1367 static index_t global_facet_id(
1368 const Mesh* mesh, index_t t, index_t lf
1369 ) {
1370 index_t t2 = mesh->cells.tet_adjacent(t, lf);
1371 if(t2 != GEO::NO_CELL && t2 > t) {
1372 index_t lf2 = mesh->cells.find_tet_adjacent(
1373 t2, t
1374 );
1375 geo_debug_assert(lf2 != GEO::NO_FACET);
1376 return index_t(4 * t2 + lf2);
1377 }
1378 return 4 * t + lf;
1379 }
1380
1381 private:
1382 GEO::vector<Triangle> triangles_;
1383 GEO::vector<Vertex> vertices_;
1384 index_t first_free_;
1385 bool v_to_t_dirty_;
1386 PointAllocator intersections_;
1387 bool symbolic_is_surface_;
1388 signed_index_t cell_id_;
1389
1390 static index_t plus1mod3_[3];
1391 static index_t minus1mod3_[3];
1392
1398 static index_t plus1mod3(index_t i) {
1399 geo_debug_assert(i < 3);
1400 return plus1mod3_[i];
1401 }
1402
1408 static index_t minus1mod3(index_t i) {
1409 geo_debug_assert(i < 3);
1410 return minus1mod3_[i];
1411 }
1412 };
1413}
1414
1415#endif
1416
A function to suppress unused parameters compilation warnings.
#define geo_debug_assert(x)
Verifies that a condition is met.
Definition assert.h:195
Generic mechanism for attributes.
Common include file, providing basic definitions. Should be included before anything else by all head...
A Corner corresponds to a vertex seen from a triangle.
Corner(index_t t_in, index_t v_in)
Creates a Corner from a triangle index and local vertex index.
Corner()
Creates an uninitialized Corner.
Represents a facet of this ConvexCell in dual form.
Vertex()
Creates a new uninitialized Vertex.
Computes the intersection between a set of halfspaces.
std::ostream & show_stats(std::ostream &os) const
Displays the number of free,used,conflict triangles.
void initialize_from_surface_mesh(Mesh *mesh, bool symbolic)
Copies a Mesh into a ConvexCell.
void mark_as_free(index_t t)
Specify that a triangle is free.
void grow()
Allocates a new triangle.
index_t create_triangle(index_t v0, index_t v1, index_t v2)
Creates a new triangle with specified vertices.
index_t create_vertex()
Creates a new vertex.
bool find_triangle_on_border(index_t conflict_begin, index_t conflict_end, index_t &t, index_t &e) const
Gets a triangle and an edge on the internal border of the conflict zone.
const GEOGen::Vertex & triangle_dual(index_t t) const
Gets the dual vertex that corresponds to a triangle.
static index_t global_facet_id(const Mesh *mesh, index_t t, index_t lf)
Computes a unique global facet id from a mesh tetrahedron and local facet index.
void set_vertex_id(index_t v, signed_index_t id)
Sets the id of a vertex.
void set_triangle_id(index_t t, signed_index_t id)
Sets the id of a triangle.
Sign side(const Mesh *mesh, const Delaunay *delaunay, const GEOGen::Vertex &v, index_t i, index_t j, bool exact) const
Tests on which side a vertex is relative to a bisector.
index_t create_triangle()
Creates a new uninitialized triangle.
index_t max_v() const
Gets the maximum valid vertex index plus one.
TriangleStatus
Represents the current state of a triangle.
void get_conflict_list(const Mesh *mesh, const Delaunay *delaunay, index_t i, index_t j, bool exact, index_t &conflict_begin, index_t &conflict_end)
Determines the conflict zone.
static double signed_bisector_distance(const Delaunay *delaunay, index_t i, index_t j, const double *q)
Evaluates the equation of a bisector at a given point.
void set_vertex_triangle(index_t v, index_t t)
Stores in a vertex the index of one the triangles incident to it.
void mark_as_used(index_t t)
Specify that a triangle is used.
bool triangle_is_free(index_t t) const
Tests whether a given triangle is free.
index_t find_furthest_point_linear_scan(const Delaunay *delaunay, index_t i, index_t j) const
Finds the index of the vertex furthest away on the negative side of a bisector.
void initialize_from_mesh_tetrahedron(const Mesh *mesh, index_t t, bool symbolic, const GEO::Attribute< double > &vertex_weight)
Assigns a mesh tetrahedron to this ConvexCell.
Sign side_exact(const Mesh *mesh, const Delaunay *delaunay, const GEOGen::Vertex &v, const double *pi, const double *pj, coord_index_t dim, bool symbolic_is_surface=false) const
Tests on which side a vertex is relative to a bisector using exact predicates.
void propagate_conflict_list(const Mesh *mesh, const Delaunay *delaunay, index_t first_t, index_t i, index_t j, bool exact, index_t &conflict_begin, index_t &conflict_end)
Computes the conflict list by propagation from a conflict triangle.
signed_index_t cell_id() const
Gets the id of this ConvexCell.
void set_cell_id(signed_index_t i)
Sets the id of this ConvexCell.
void set_symbolic_is_surface(bool x)
Specifies that symbolic information is relative to surfacic mesh (rather than volumetric mesh).
signed_index_t clip_by_plane(const Mesh *mesh, const Delaunay *delaunay, index_t i, index_t j, bool exact, bool symbolic)
Clips this ConvexCell with a plane.
index_t nb_t() const
Gets the number of used triangles.
signed_index_t vertex_id(index_t v) const
Gets the id of a vertex.
void clear()
Clears this ConvexCell.
void set_next_triangle(index_t t, index_t t2)
Sets the successor of a triangle.
void merge_into_free_list(index_t list_begin, index_t list_end)
Merges a list of triangles into the free list.
index_t triangle_adjacent_index(index_t t1, index_t t2) const
Finds the edge along which two triangles are adjacent.
void copy(const ConvexCell &rhs)
Copies a ConvexCell.
index_t triangle_adjacent(index_t t, index_t e) const
Gets the index of a triangle adjacent to another one.
bool triangle_is_valid(index_t t) const
Tests whether a given triangle is valid.
void convert_to_mesh(Mesh *mesh, bool copy_symbolic_info=false)
Copies a ConvexCell into a Mesh.
index_t find_triangle_vertex(index_t t, index_t v) const
Finds the local index of a triangle vertex.
void mark_as_conflict(index_t t)
Specify that a triangle belongs to the conflict zone.
index_t create_triangle_copy(const double *p, index_t v0, index_t v1, index_t v2, index_t t0, index_t t1, index_t t2)
Creates a new triangles with specified vertices, adjacent triangles and geometric location at the dua...
index_t create_triangle(index_t v0, index_t v1, index_t v2, index_t t0, index_t t1, index_t t2)
Creates a new triangles with specified vertices and adjacent triangles.
void set_triangle_adjacent(index_t t, index_t e, index_t t2)
Sets a triangle adjacency.
void append_triangle_to_conflict_list(index_t t, index_t &conflict_begin, index_t &conflict_end)
Appends a triangle to the conflict list.
coord_index_t dimension() const
Gets the dimension of this ConvexCell.
void set_triangle(index_t t, index_t v0, index_t v1, index_t v2, index_t t0, index_t t1, index_t t2)
Sets the vertices and adjacent triangles of a triangle.
void set_triangle_vertex(index_t t, index_t iv, index_t v)
Sets a vertex of a triangle.
ConvexCell(coord_index_t dim)
ConvexCell constructor.
void init_v_to_t()
Updates the cache that stores for each vertex a triangle incident to it.
signed_index_t triangle_id(index_t t) const
Gets the id of a triangle.
index_t triangulate_hole(const Delaunay *delaunay, index_t i, index_t j, bool symbolic, index_t t1, index_t t1ebord, index_t v_in)
Triangulates the conflict zone.
bool triangle_is_used(index_t t) const
Tests whether a given triangle is used.
void move_to_next_around_vertex(Corner &c) const
Replaces a corner by the next corner obtained by turing around the vertex.
index_t create_triangle(const double *p, double w, index_t v0, index_t v1, index_t v2, index_t t0, index_t t1, index_t t2)
Creates a new triangles with specified vertices, adjacent triangles and geometric location at the dua...
index_t next_triangle(index_t t) const
Gets the successor of a triangle.
GEOGen::Vertex & triangle_dual(index_t t)
Gets the dual vertex that corresponds to a triangle.
bool triangle_is_conflict(index_t t) const
Tests whether a given triangle belongs to the conflict zone.
signed_index_t vertex_triangle(index_t v) const
Gets one of the triangles incident to a vertex.
index_t max_t() const
Gets the maximum valid triangle index plus one.
index_t triangle_vertex(index_t t, index_t iv) const
Gets the index of a triangle vertex.
An allocator for points that are created from intersections in GenericVoronoiDiagram.
Internal representation of vertices in GenericVoronoiDiagram.
Sign side_fast(const double *p1, const double *p2) const
Computes the side of this vertex relative to a bisector.
Manages an attribute attached to a set of object.
Abstract interface for Delaunay triangulation in Nd.
Definition delaunay.h:71
const double * vertex_ptr(index_t i) const
Gets a pointer to a vertex by its global index.
Definition delaunay.h:233
Represents a mesh.
Definition mesh.h:2648
Vector with aligned memory allocation.
Definition memory.h:623
Types and utilities for manipulating vertices in geometric and symbolic forms in restricted Voronoi d...
T geo_sqr(T x)
Gets the square value of a value.
Definition numeric.h:259
void geo_argused(const T &)
Suppresses compiler warnings about unused parameters.
Definition argused.h:60
@ ZERO
Definition numeric.h:228
@ NEGATIVE
Definition numeric.h:226
Represents a vertex of this ConvexCell in dual form.
Triangle()
Creates a new uninitialized Triangle.
Triangle(index_t v0, index_t v1, index_t v2, index_t f0, index_t f1, index_t f2)
Creates a new triangle.