GCC Code Coverage Report


Directory: ./
File: lib/geogram/delaunay/CDT_2d.cpp
Date: 2026-09-07 02:36:43
Exec Total Coverage
Lines: 776 1042 74.5%
Functions: 50 64 78.1%
Branches: 575 1678 34.3%

Line Branch Exec Source
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 // Reference: S. W. Sloan, a fast algorithm for generating
41 // constrained Delaunay triangulation, 1992,
42 // Computers and Structures
43 //
44 // Specificities of this implementation:
45 //
46 // - Edges are systematically manipulated through triangles,
47 // and these triangles are rotated in-place in the mesh,
48 // in such a way that the edge we are talking about is
49 // systematically edge 0 (with vertices 1 and 2).
50 //
51 // - The constraint-enforcing step manipulates a queue Q
52 // of edges encoded this way. It examines pairs of triangles
53 // t1,t2=Tadj(t1,0), decides whether to swap their common
54 // edge (based on convexity test and intersection of
55 // t1's edge 0 with the constraint). In fact, this intersection
56 // test only depends on the combinatorics of (t1,t2) (two cases)
57 // and the position of t1' vertex 0 relative to the constraint
58 // (two cases), that makes 4 cases in total. In these cases,
59 // - t1 can either leave Q or be enqueued again
60 // - t2 was always in Q already (because it has an edge
61 // that has an intersection with the constraint), but
62 // there is one case where it leaves Q
63 // DList has an O(1) function to test whether an element is in the list (using
64 // flags associated with the elements). It is used in one case: when t2 is
65 // is not in Q, it means there is no intersection.
66
67 #include <geogram/delaunay/CDT_2d.h>
68 #include <geogram/mesh/mesh_reorder.h>
69 #include <geogram/basic/numeric.h>
70 #include <geogram/basic/boolean_expression.h>
71
72 #ifndef GEOGRAM_PSM
73 #include <geogram/mesh/mesh.h>
74 #include <geogram/mesh/mesh_io.h>
75 #endif
76
77 // Used by debugging functions and statistics
78 #include <geogram/mesh/index.h>
79 #include <geogram/basic/debug_stream.h>
80 #include <set>
81 #include <deque>
82 #include <stack>
83 //#define CDT_NAIVE // use naive per-edge method (kept for reference/debugging)
84
85 #ifdef GEO_DEBUG
86 //#define CDT_DEBUG // display *lots* of messages and activates costly checks
87 #endif
88
89 #ifdef CDT_DEBUG
90 #define CDT_LOG(X) std::cerr << X << std::endl
91 #else
92 #define CDT_LOG(X)
93 #endif
94
95 namespace GEO {
96
97 651 CDTBase2d::CDTBase2d() :
98 651 nv_(0),
99 651 ncnstr_(0),
100 651 delaunay_(true),
101 651 exact_incircle_(true),
102 651 exact_intersections_(true) {
103 651 }
104
105 1302 CDTBase2d::~CDTBase2d() {
106 1302 }
107
108 23572 void CDTBase2d::clear() {
109 23572 nv_ = 0;
110 23572 ncnstr_ = 0;
111 23572 T_.resize(0);
112 23572 Tadj_.resize(0);
113 23572 v2T_.resize(0);
114 23572 Tflags_.resize(0);
115 23572 Tecnstr_first_.resize(0);
116 23572 ecnstr_val_.resize(0);
117 23572 ecnstr_next_.resize(0);
118 23572 Tnext_.resize(0);
119 23572 Tprev_.resize(0);
120 23572 }
121
122 15537 void CDTBase2d::create_enclosing_triangle(
123 index_t v0, index_t v1, index_t v2
124 ) {
125 15537 nv_ = 3;
126 15537 v2T_.resize(3);
127
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 15537 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
15537 geo_debug_assert(v0 <= 3);
128
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 15537 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
15537 geo_debug_assert(v1 <= 3);
129
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 15537 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
15537 geo_debug_assert(v2 <= 3);
130 15537 index_t t0 = Tnew();
131 15537 Tset(t0, v0, v1, v2, NO_INDEX, NO_INDEX, NO_INDEX);
132 15537 orient_012_ = orient2d(0,1,2);
133
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 15537 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
15537 geo_assert(orient_012_ != ZERO);
134 15537 }
135
136 8173 void CDTBase2d::create_enclosing_quad(
137 index_t v0, index_t v1, index_t v2, index_t v3
138 ) {
139 8173 nv_ = 4;
140 8173 v2T_.resize(4);
141
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 8173 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
8173 geo_debug_assert(v0 <= 4);
142
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 8173 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
8173 geo_debug_assert(v1 <= 4);
143
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 8173 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
8173 geo_debug_assert(v2 <= 4);
144
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 8173 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
8173 geo_debug_assert(v3 <= 4);
145 8173 index_t t0 = Tnew();
146 8173 index_t t1 = Tnew();
147 8173 Tset(t0, v0, v1, v3, t1, NO_INDEX, NO_INDEX);
148 8173 Tset(t1, v3, v1, v2, NO_INDEX, NO_INDEX, t0);
149 8173 orient_012_ = orient2d(0,1,2);
150
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 8173 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
8173 geo_debug_assert(is_convex_quad(t0));
151
2/2
✓ Branch 1 taken 6130 times.
✓ Branch 2 taken 2043 times.
8173 if(Sign(incircle(v0,v1,v2,v3)*orient_012_) == POSITIVE) {
152 6130 swap_edge(t0);
153 }
154 8173 }
155
156 222718 void CDTBase2d::begin_insert_transaction() {
157 222718 }
158
159 222712 void CDTBase2d::commit_insert_transaction() {
160 222712 }
161
162 6 void CDTBase2d::rollback_insert_transaction() {
163 6 }
164
165 515221 index_t CDTBase2d::insert(index_t v, index_t hint) {
166 515221 bool keep_duplicates = false;
167
2/2
✓ Branch 1 taken 286738 times.
✓ Branch 2 taken 228483 times.
515221 if(v == nv()) {
168
1/2
✓ Branch 1 taken 286738 times.
✗ Branch 2 not taken.
286738 v2T_.push_back(NO_INDEX);
169 286738 ++nv_;
170 } else {
171 // We are inserting a vertex in the middle of the
172 // list, which means we are doing batch-insertion,
173 // then we will not discard duplicates.
174 228483 keep_duplicates = true;
175
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 228483 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
228483 geo_debug_assert(v < nv_);
176 }
177
178 // Phase 1: find triangle that contains vertex i
179
1/2
✓ Branch 1 taken 515221 times.
✗ Branch 2 not taken.
515221 begin_insert_transaction();
180 Sign o[3];
181
1/2
✓ Branch 1 taken 515221 times.
✗ Branch 2 not taken.
515221 index_t t = locate(v,hint,o);
182 515221 int nb_z = (o[0] == ZERO) + (o[1] == ZERO) + (o[2] == ZERO);
183
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 515221 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
515221 geo_debug_assert(nb_z != 3);
184
185 // Duplicated vertex
186
2/2
✓ Branch 0 taken 157665 times.
✓ Branch 1 taken 357556 times.
515221 if(nb_z == 2) {
187 CDT_LOG("duplicated vertex");
188
3/4
✓ Branch 0 taken 69208 times.
✓ Branch 1 taken 88457 times.
✓ Branch 3 taken 69208 times.
✗ Branch 4 not taken.
157665 v = (o[0] != ZERO) ? Tv(t,0) :
189
3/4
✓ Branch 0 taken 44406 times.
✓ Branch 1 taken 44051 times.
✓ Branch 3 taken 44406 times.
✗ Branch 4 not taken.
88457 (o[1] != ZERO) ? Tv(t,1) :
190
1/2
✓ Branch 1 taken 44051 times.
✗ Branch 2 not taken.
44051 Tv(t,2) ;
191
2/2
✓ Branch 0 taken 151894 times.
✓ Branch 1 taken 5771 times.
157665 if(!keep_duplicates) {
192 151894 v2T_.pop_back();
193 151894 --nv_;
194 }
195 // Used by optional predicate cache management in derived classes.
196 // locate() computed some orient_2d() predicates, that may be
197 // stored in a temporary buffer, discard it.
198
1/2
✓ Branch 1 taken 157665 times.
✗ Branch 2 not taken.
157665 rollback_insert_transaction();
199 157665 return v;
200 }
201
202 // Used by optional predicate cache management in derived classes.
203 // Copy the computed orient_2d() values to predicate cache.
204
1/2
✓ Branch 1 taken 357556 times.
✗ Branch 2 not taken.
357556 commit_insert_transaction();
205
206 // Stack of triangle edges to examine for flipping. Ignored in
207 // non-Delaunay mode (ignored if !delaunay_)
208 // Note: it is always edge 0 that we examine, since new
209 // triangles are always created with v as vertex 0.
210 357556 DList S(*this);
211
2/2
✓ Branch 0 taken 246200 times.
✓ Branch 1 taken 111356 times.
357556 if(delaunay_) {
212
1/2
✓ Branch 1 taken 246200 times.
✗ Branch 2 not taken.
246200 S.initialize(DLIST_S_ID);
213 }
214
215 // Phase 2: split triangle
216 // Particular case: v is on edge
217
2/2
✓ Branch 0 taken 74627 times.
✓ Branch 1 taken 282929 times.
357556 if(nb_z == 1) {
218 CDT_LOG("insert vertex on edge");
219
2/2
✓ Branch 0 taken 34305 times.
✓ Branch 1 taken 40322 times.
108932 index_t le = (o[0] == ZERO) ? 0 :
220
2/2
✓ Branch 0 taken 18156 times.
✓ Branch 1 taken 16149 times.
34305 (o[1] == ZERO) ? 1 :
221 2 ;
222
1/2
✓ Branch 1 taken 74627 times.
✗ Branch 2 not taken.
74627 insert_vertex_in_edge(v,t,le,S);
223 } else {
224 CDT_LOG("insert vertex in triangle");
225
1/2
✓ Branch 1 taken 282929 times.
✗ Branch 2 not taken.
282929 insert_vertex_in_triangle(v,t,S);
226 }
227
228 // Phase 3: recursively restore Delaunay conditions for the neighbors
229 // of the new vertex
230
2/2
✓ Branch 0 taken 246200 times.
✓ Branch 1 taken 111356 times.
357556 if(delaunay_) {
231
1/2
✓ Branch 1 taken 246200 times.
✗ Branch 2 not taken.
246200 Delaunayize_vertex_neighbors(v,S);
232 }
233
234 #ifdef CDT_DEBUG
235 debug_check_consistency();
236 #endif
237 357556 return v;
238 357556 }
239
240
241 179300 void CDTBase2d::insert_constraint(index_t i, index_t j) {
242
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 179300 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
179300 geo_debug_assert(i < nv());
243
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 179300 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
179300 geo_debug_assert(j < nv());
244 CDT_LOG("insert constraint: " << i << "-" << j);
245 #ifdef CDT_DEBUG
246 debug_check_consistency();
247 #endif
248 179300 ++ncnstr_;
249
250 // Index of first vertex coming from constraints intersection
251 // (keep track of it to re-Delaunayize their neighborhoods).
252 179300 index_t first_v_isect = nv_;
253
254 #ifndef CDT_NAIVE
255
1/2
✓ Branch 1 taken 179300 times.
✗ Branch 2 not taken.
179300 DList Q(*this, DLIST_Q_ID); // Queue of edges to constrain
256 179300 DList N(*this); // New edges to re-Delaunayize (ignored if !delaunay_)
257
2/2
✓ Branch 0 taken 178873 times.
✓ Branch 1 taken 427 times.
179300 if(delaunay_) {
258
1/2
✓ Branch 1 taken 178873 times.
✗ Branch 2 not taken.
178873 N.initialize(DLIST_N_ID);
259 }
260
2/2
✓ Branch 0 taken 258470 times.
✓ Branch 1 taken 179300 times.
437770 while(i != j) {
261
262 // Step 1: find all the edges that have an intersection
263 // with the constraint [i,j], enqueue them in Q.
264 // Stop at vertex on constraint or constraint intersection
265 // if any (returned in k)
266
1/2
✓ Branch 1 taken 258470 times.
✗ Branch 2 not taken.
258470 index_t k = find_intersected_edges(i,j,Q);
267
268 // If we found a constraint intersection,
269 // we need to Delaunayize the neigborhood
270 // of the newly created vertex. Then we
271 // need to find the intersected edges again,
272 // since they may have changed.
273
6/6
✓ Branch 0 taken 257957 times.
✓ Branch 1 taken 513 times.
✓ Branch 2 taken 257444 times.
✓ Branch 3 taken 513 times.
✓ Branch 4 taken 5765 times.
✓ Branch 5 taken 251679 times.
258470 if(delaunay_ && exact_intersections_ && k >= first_v_isect) {
274
2/8
✓ Branch 1 taken 5765 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 5765 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
5765 geo_debug_assert(insert(k) == k);
275
1/2
✓ Branch 1 taken 5765 times.
✗ Branch 2 not taken.
5765 Q.clear();
276
1/2
✓ Branch 1 taken 5765 times.
✗ Branch 2 not taken.
5765 Delaunayize_vertex_neighbors(k);
277 #ifdef CDT_DEBUG
278 debug_check_geometry();
279 #endif
280
1/2
✓ Branch 1 taken 5765 times.
✗ Branch 2 not taken.
5765 index_t new_k = find_intersected_edges(i,j,Q);
281
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 5765 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
5765 geo_assert(new_k == k);
282 }
283
284 // Step 2: constrain edges
285
1/2
✓ Branch 1 taken 258470 times.
✗ Branch 2 not taken.
258470 constrain_edges(i,k,Q,N);
286
287 // Step 3: restore Delaunay condition
288
2/2
✓ Branch 0 taken 257957 times.
✓ Branch 1 taken 513 times.
258470 if(delaunay_) {
289
1/2
✓ Branch 1 taken 257957 times.
✗ Branch 2 not taken.
257957 Delaunayize_new_edges(N);
290 #ifdef CDT_DEBUG
291 debug_check_geometry();
292 #endif
293 }
294
295 258470 i = k;
296 }
297 #else
298 DList Q(*this, DLIST_Q_ID); // Queue of edges to constrain
299 vector<Edge> N; // New edges to re-Delaunayize
300 while(i != j) {
301 index_t k = find_intersected_edges(i,j,Q);
302
303 // If we found a constraint intersection,
304 // we need to Delaunayize the neigborhood
305 // of the newly created vertex. Then we
306 // need to find the intersected edges again,
307 // since they may have changed.
308 if(delaunay_ && exact_intersections_ && k >= first_v_isect) {
309 geo_debug_assert(insert(k) == k);
310 Q.clear();
311 Delaunayize_vertex_neighbors(k);
312 #ifdef CDT_DEBUG
313 debug_check_geometry();
314 #endif
315 index_t new_k = find_intersected_edges(i,j,Q);
316 geo_assert(new_k == k);
317 }
318
319 constrain_edges_naive(i,k,Q,N);
320 debug_check_combinatorics();
321 if(delaunay_) {
322 Delaunayize_new_edges_naive(N);
323 #ifdef CDT_DEBUG
324 debug_check_geometry();
325 #endif
326 }
327 debug_check_combinatorics();
328 i = k;
329 }
330 #endif
331 // Delaunayize neighborhood of vertices yielded by constraint
332 // intersections now if not done before (that is, if intersections
333 // are not exact, like in the default CDT2d class). If intersections
334 // are exact, it was done before (right after creating the intersection)
335
4/4
✓ Branch 0 taken 178873 times.
✓ Branch 1 taken 427 times.
✓ Branch 2 taken 427 times.
✓ Branch 3 taken 178446 times.
179300 if(delaunay_ && !exact_intersections_) {
336
2/2
✓ Branch 1 taken 84 times.
✓ Branch 2 taken 427 times.
511 for(index_t v=first_v_isect; v<nv(); ++v) {
337
1/2
✓ Branch 1 taken 84 times.
✗ Branch 2 not taken.
84 Delaunayize_vertex_neighbors(v);
338 }
339 }
340
341 #ifdef CDT_DEBUG
342 debug_check_consistency();
343 #endif
344 179300 }
345
346 5849 void CDTBase2d::Delaunayize_vertex_neighbors(index_t v) {
347 CDT_LOG("Delaunayize_vertex_neighbors " << v);
348
349 // Delaunayize triangles around vertices coming from
350 // constraint intersections
351
1/2
✓ Branch 1 taken 5849 times.
✗ Branch 2 not taken.
5849 DList S(*this, DLIST_S_ID);
352
353
2/8
✓ Branch 1 taken 5849 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 5849 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
5849 geo_assert(vT(v) != NO_INDEX);
354
355 // We cannot use for_each_triangle_around_vertex()
356 // because we need to Trot() t during traveral,
357 // to have v has t's vertex 0
358 // But the good news is that v is never on the border,
359 // (because it comes from an edge *intersection*),
360 // hence traversal is easier.
361
1/2
✓ Branch 1 taken 5849 times.
✗ Branch 2 not taken.
5849 index_t t0 = vT(v); // Need to store it, because we Trot()
362 5849 index_t t = t0;
363 do {
364
1/2
✓ Branch 1 taken 23858 times.
✗ Branch 2 not taken.
23858 index_t lv = Tv_find(t,v);
365
1/2
✓ Branch 1 taken 23858 times.
✗ Branch 2 not taken.
23858 Trot(t,lv);
366
2/8
✓ Branch 1 taken 23858 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 23858 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
23858 geo_debug_assert(Tv(t,0) == v);
367
1/2
✓ Branch 1 taken 23858 times.
✗ Branch 2 not taken.
23858 S.push_back(t);
368
1/2
✓ Branch 1 taken 23858 times.
✗ Branch 2 not taken.
23858 t = Tadj(t, 1);
369
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 23858 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
23858 geo_assert(t != NO_INDEX);
370
2/2
✓ Branch 0 taken 18009 times.
✓ Branch 1 taken 5849 times.
23858 } while(t != t0);
371
1/2
✓ Branch 1 taken 5849 times.
✗ Branch 2 not taken.
5849 Delaunayize_vertex_neighbors(v,S);
372 5849 }
373
374 /**
375 * \brief Used by the implementation of find_intersected_edges()
376 * \details During traversal of a constrained edge [i,j], we can be on
377 * a vertex (then v != NO_INDEX) or
378 * on a triangle (then t != NO_INDEX).
379 * We also keep track of the previous vertex (prev_v) and previous
380 * triangle(prev_t) in order to make sure we do not go backwards.
381 */
382 struct CDT2d_ConstraintWalker {
383 /**
384 * \brief ConstraintWalker constructor
385 * \param[in] i_in , j_in extremities of the constrained edge
386 */
387 264235 CDT2d_ConstraintWalker(index_t i_in, index_t j_in) :
388 264235 i(i_in), j(j_in),
389 264235 t_prev(NO_INDEX), v_prev(NO_INDEX),
390 264235 t(NO_INDEX), v(i_in),
391 264235 v_cnstr(NO_INDEX)
392 {
393 264235 }
394 index_t i, j;
395 index_t t_prev, v_prev;
396 index_t t, v;
397 index_t v_cnstr;
398 };
399
400 264235 index_t CDTBase2d::find_intersected_edges(index_t i, index_t j, DList& Q) {
401 CDT_LOG("Find intersected edges: " << i << "-" << j);
402 264235 CDT2d_ConstraintWalker W(i,j);
403 // Stop at the first encountered vertex or constraint intersection.
404
4/4
✓ Branch 0 taken 264235 times.
✓ Branch 1 taken 324107 times.
✓ Branch 2 taken 59872 times.
✓ Branch 3 taken 264235 times.
588342 while(W.v == i || W.v == NO_INDEX) {
405 CDT_LOG(
406 " t=" << int(W.t) << " v=" << int(W.v) << " "
407 "t_prev=" << int(W.t_prev) << " v_prev=" << int(W.v_prev)
408 << " "
409 );
410
2/2
✓ Branch 0 taken 264235 times.
✓ Branch 1 taken 59872 times.
324107 if(W.v != NO_INDEX) {
411
1/2
✓ Branch 1 taken 264235 times.
✗ Branch 2 not taken.
264235 walk_constraint_v(W);
412 } else {
413
1/2
✓ Branch 1 taken 59872 times.
✗ Branch 2 not taken.
59872 walk_constraint_t(W,Q);
414 }
415 }
416 264235 return W.v;
417 }
418
419 // The two functions below are more complicated than I wished, but are
420 // simpler than it looks like. There are two main different cases:
421 //
422 // - walk_constraint_v(): we are on a vertex.
423 // traverse all the triangles around v and find the one that
424 // has an intersection. For instance, when we start from vertex i,
425 // and also when the previous step encountered a vertex exactly on
426 // the constrained segment. It is the annoying case where one has
427 // to traverse the triangles incident to v (using the function
428 // for_each_T_around_v() that takes a lambda).
429 //
430 // - walk_constraint_t(): we are on an edge intersection.
431 // propagate to the neighbor of t accross the intersected edge.
432 // It is the "generic" case, simpler (the next triangle is
433 // determined by the edge of t that is intersected).
434 //
435 // There are three things that makes things slightly
436 // more complicated:
437 // - each case has two sub-cases, depending on whether the next
438 // intersection is an existing vertex.
439 // - if an existing edge is embedded in the constraint, one needs
440 // to flag that edge as a constraint.
441 // - we need to test whether we are arrived at vertex j
442
443 264235 void CDTBase2d::walk_constraint_v(CDT2d_ConstraintWalker& W) {
444
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 264235 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
264235 geo_debug_assert(W.v != NO_INDEX);
445
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 264235 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
264235 geo_debug_assert(W.t == NO_INDEX);
446
447 264235 index_t t_next = NO_INDEX;
448 264235 index_t v_next = NO_INDEX;
449
450
1/2
✓ Branch 1 taken 264235 times.
✗ Branch 2 not taken.
264235 for_each_T_around_v(
451
1/2
✓ Branch 1 taken 264235 times.
✗ Branch 2 not taken.
528470 W.v, [&](index_t t_around_v, index_t le) {
452
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 518992 times.
518992 if(t_around_v == W.t_prev) { // Don't go backwards !
453 return false;
454 }
455 518992 index_t v1 = Tv(t_around_v, (le + 1)%3);
456 518992 index_t v2 = Tv(t_around_v, (le + 2)%3);
457
4/4
✓ Branch 0 taken 429772 times.
✓ Branch 1 taken 89220 times.
✓ Branch 2 taken 72955 times.
✓ Branch 3 taken 356817 times.
518992 if(v1 == W.j || v2 == W.j) { // Are we arrived at j ?
458 162175 v_next = W.j;
459 // Edge is flagged as constraint here, because
460 // it will not be seen by constraint enforcement.
461
2/2
✓ Branch 0 taken 89220 times.
✓ Branch 1 taken 72955 times.
162175 index_t le_cnstr_edge = (v1 == W.j) ? (le+2)%3 : (le+1)%3;
462 162175 Tadd_edge_cnstr_with_neighbor(
463 162175 t_around_v, le_cnstr_edge, ncnstr_-1
464 );
465 CDT_LOG(
466 " During cnstr " << W.i << "-" << W.j << ": "
467 << " Constrained edge "
468 << Tv(t_around_v, (le_cnstr_edge+1)%3) << "-"
469 << Tv(t_around_v, (le_cnstr_edge+2)%3)
470 );
471 162175 return true;
472 }
473 356817 Sign o1 = orient2d(W.i,W.j,v1);
474 356817 Sign o2 = orient2d(W.i,W.j,v2);
475 356817 Sign o3 = orient2d(v1,v2,W.j);
476 356817 Sign o4 = orient_012_; // equivalent to orient2d(v1,v2,i)
477
4/4
✓ Branch 0 taken 40739 times.
✓ Branch 1 taken 316078 times.
✓ Branch 2 taken 17056 times.
✓ Branch 3 taken 23683 times.
356817 if(o1*o2 < 0 && o3*o4 < 0) {
478 17056 Trot(t_around_v,le); // so that le becomes edge 0
479 17056 t_next = t_around_v; // will be added to Q during next round
480 17056 return true;
481 } else {
482 // Special case: v1 or v2 is exactly on [i,j]
483 // Edge is flagged as constraint here, because
484 // it will not be seen by constraint enforcement.
485
3/8
✓ Branch 0 taken 91514 times.
✓ Branch 1 taken 248247 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 91514 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
339761 geo_debug_assert(o1 != ZERO || o2 != ZERO);
486
5/6
✓ Branch 0 taken 91514 times.
✓ Branch 1 taken 248247 times.
✓ Branch 2 taken 78559 times.
✓ Branch 3 taken 12955 times.
✓ Branch 4 taken 78559 times.
✗ Branch 5 not taken.
339761 if(o1 == ZERO && o3*o4 < 0 && v1 != W.v_prev) {
487 78559 v_next = v1;
488 78559 Tadd_edge_cnstr_with_neighbor(
489 78559 t_around_v, (le + 2)%3, ncnstr_-1
490 );
491 78559 return true;
492
5/6
✓ Branch 0 taken 54875 times.
✓ Branch 1 taken 206327 times.
✓ Branch 2 taken 6445 times.
✓ Branch 3 taken 48430 times.
✓ Branch 4 taken 6445 times.
✗ Branch 5 not taken.
261202 } else if(o2 == ZERO && o3*o4 < 0 && v2 != W.v_prev) {
493 6445 v_next = v2;
494 6445 Tadd_edge_cnstr_with_neighbor(
495 6445 t_around_v, (le + 1)%3, ncnstr_-1
496 );
497 6445 return true;
498 }
499 }
500 254757 return false;
501 }
502 ); // End of for_each_T_around_v() loop
503 264235 W.t_prev = W.t;
504 264235 W.v_prev = W.v;
505 264235 W.t = t_next;
506 264235 W.v = v_next;
507 264235 }
508
509 59872 void CDTBase2d::walk_constraint_t(CDT2d_ConstraintWalker& W, DList& Q) {
510
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 59872 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
59872 geo_debug_assert(W.v == NO_INDEX);
511
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 59872 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
59872 geo_debug_assert(W.t != NO_INDEX);
512
513 59872 index_t v_next = NO_INDEX;
514 59872 index_t t_next = NO_INDEX;
515
516
8/8
✓ Branch 1 taken 57492 times.
✓ Branch 2 taken 2380 times.
✓ Branch 4 taken 53689 times.
✓ Branch 5 taken 3803 times.
✓ Branch 7 taken 3744 times.
✓ Branch 8 taken 49945 times.
✓ Branch 9 taken 9927 times.
✓ Branch 10 taken 49945 times.
59872 if(Tv(W.t,0) == W.j || Tv(W.t,1) == W.j || Tv(W.t,2) == W.j) {
517 9927 v_next = W.j; // Are we arrived at j ?
518 } else {
519 // Test the three edges of the triangle
520
1/2
✓ Branch 0 taken 78104 times.
✗ Branch 1 not taken.
78104 for(index_t le = 0; le<3; ++le) {
521
2/2
✓ Branch 1 taken 15829 times.
✓ Branch 2 taken 62275 times.
78104 if(Tadj(W.t,le) == W.t_prev) { // Do not go backwards !
522 15829 continue;
523 }
524 // Test whether [v1,v2] intersects the support line of (i,j).
525 // No need to test the *segment* [i,j]: we know the line enters
526 // the triangle (it is how we came here), and we know it leaves
527 // it, else j would have been one of the triangle's vertices.
528 62275 index_t v1 = Tv(W.t, (le + 1)%3);
529 62275 index_t v2 = Tv(W.t, (le + 2)%3);
530 62275 Sign o1 = orient2d(W.i,W.j,v1);
531 62275 Sign o2 = orient2d(W.i,W.j,v2);
532
2/2
✓ Branch 0 taken 48749 times.
✓ Branch 1 taken 13526 times.
62275 if(o1*o2 < 0) {
533 // [v1,v2] has a frank intersection with [i,j]
534 48749 Trot(W.t,le); // So that edge 0 is intersected edge
535
2/2
✓ Branch 1 taken 5933 times.
✓ Branch 2 taken 42816 times.
48749 if(Tedge_is_constrained(W.t,0)) {
536 CDT_LOG(" Cnstr isect with:" << v1 << "-" << v2);
537 11866 v_next = create_intersection(
538 5933 ncnstr()-1, W.i, W.j,
539 edge_cnstr(Tedge_cnstr_first(W.t,0)), v1, v2
540 );
541 5933 insert_vertex_in_edge(v_next,W.t,0);
542 // Mark new edge as constraint if walker was previously
543 // on a vertex.
544
2/2
✓ Branch 0 taken 3373 times.
✓ Branch 1 taken 2560 times.
5933 if(W.v_prev != NO_INDEX) {
545 3373 Tadd_edge_cnstr_with_neighbor(W.t,2,ncnstr_-1);
546 }
547 } else {
548 CDT_LOG(" Isect: t=" << W.t <<" E=" << v1 <<"-"<< v2);
549 42816 Q.push_back(W.t);
550 42816 t_next = Tadj(W.t,0);
551 }
552 48749 break;
553 } else { // Special case: v1 or v2 is exactly on [i,j]
554
3/8
✓ Branch 0 taken 17 times.
✓ Branch 1 taken 13509 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 17 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
13526 geo_debug_assert(o1 != ZERO || o2 != ZERO);
555
2/2
✓ Branch 0 taken 17 times.
✓ Branch 1 taken 13509 times.
13526 if(o1 == ZERO) {
556 17 v_next = v1;
557 17 break;
558
2/2
✓ Branch 0 taken 1179 times.
✓ Branch 1 taken 12330 times.
13509 } else if(o2 == ZERO) {
559 1179 v_next = v2;
560 1179 break;
561 }
562 }
563 }
564 }
565 59872 W.t_prev = W.t;
566 59872 W.v_prev = W.v;
567 59872 W.t = t_next;
568 59872 W.v = v_next;
569 59872 }
570
571 258470 void CDTBase2d::constrain_edges(index_t i, index_t j, DList& Q, DList& N) {
572
573 #ifdef CDT_DEBUG
574 // The function find_edge_intersections() is super complicated,
575 // so in debug mode I make sure it did its job correctly (by testing
576 // *all* edge intersections).
577 check_edge_intersections(i,j,Q);
578 #endif
579 // Called each time edge le of triangle t has no isect with cnstr,
580 // (then it is a "new edge")
581 38846 auto new_edge = [&](index_t t,index_t le) {
582 38846 Trot(t,le);
583 38846 if(
584
6/6
✓ Branch 1 taken 12732 times.
✓ Branch 2 taken 26114 times.
✓ Branch 4 taken 1464 times.
✓ Branch 5 taken 11268 times.
✓ Branch 6 taken 11268 times.
✓ Branch 7 taken 27578 times.
66424 (Tv(t,1) == i && Tv(t,2) == j) ||
585
1/4
✗ Branch 1 not taken.
✓ Branch 2 taken 27578 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
27578 (Tv(t,1) == j && Tv(t,2) == i)
586 ) {
587 // Set constraint flag if the new edge is the constrained edge
588 11268 Tadd_edge_cnstr_with_neighbor(t,0,ncnstr_-1);
589 } else {
590 // Memorize new edge as "to be Delaunayized"
591
2/2
✓ Branch 1 taken 11210 times.
✓ Branch 2 taken 16368 times.
27578 if(N.initialized()) {
592 11210 N.push_back(t);
593 }
594 }
595 297316 };
596
597 // Called each time edge le of triangle t still has an isect with cnstr
598 // (then it is queued again)
599 18154 auto isect_edge = [&](index_t t, index_t le) {
600 18154 Trot(t,le);
601 18154 Q.push_front(t);
602 276624 };
603
604
3/4
✓ Branch 1 taken 336782 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 78312 times.
✓ Branch 4 taken 258470 times.
336782 while(!Q.empty()) {
605
1/2
✓ Branch 1 taken 78312 times.
✗ Branch 2 not taken.
78312 index_t t1 = Q.pop_back();
606
3/4
✓ Branch 1 taken 78312 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 21312 times.
✓ Branch 4 taken 57000 times.
78312 if(!is_convex_quad(t1)) {
607 // Sanity check: if the only remaining edge to flip does
608 // not form a convex quad, it means we are going to
609 // flip forever ! (shoud not happen)
610
2/8
✓ Branch 1 taken 21312 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 21312 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
21312 geo_assert(!Q.empty());
611
1/2
✓ Branch 1 taken 21312 times.
✗ Branch 2 not taken.
21312 Q.push_front(t1);
612 } else {
613
1/2
✓ Branch 1 taken 57000 times.
✗ Branch 2 not taken.
57000 index_t t2 = Tadj(t1,0);
614
1/2
✓ Branch 1 taken 57000 times.
✗ Branch 2 not taken.
57000 bool no_isect = !Q.contains(t2);
615
1/2
✓ Branch 1 taken 57000 times.
✗ Branch 2 not taken.
57000 index_t v0 = Tv(t1,0);
616
7/10
✓ Branch 1 taken 57000 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 38129 times.
✓ Branch 4 taken 18871 times.
✓ Branch 6 taken 38129 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 38129 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 18771 times.
✓ Branch 12 taken 19358 times.
57000 bool t2v0_t1v2 = (Q.contains(t2) && Tv(t2,0) == Tv(t1,2));
617
7/10
✓ Branch 1 taken 57000 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 38129 times.
✓ Branch 4 taken 18871 times.
✓ Branch 6 taken 38129 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 38129 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 19358 times.
✓ Branch 12 taken 18771 times.
57000 bool t2v0_t1v1 = (Q.contains(t2) && Tv(t2,0) == Tv(t1,1));
618 57000 geo_argused(t2v0_t1v1);
619
620
2/2
✓ Branch 0 taken 18871 times.
✓ Branch 1 taken 38129 times.
57000 if(no_isect) {
621
1/2
✓ Branch 1 taken 18871 times.
✗ Branch 2 not taken.
18871 swap_edge(t1);
622
2/8
✓ Branch 1 taken 18871 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 18871 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
18871 geo_debug_assert(!segment_edge_intersect(i,j,t1,2));
623
1/2
✓ Branch 1 taken 18871 times.
✗ Branch 2 not taken.
18871 new_edge(t1,2);
624 } else {
625 // See comment at beginning of file
626 // (a small variation in Sloan's
627 // method that makes better use of the combinatorics)
628
1/2
✓ Branch 1 taken 38129 times.
✗ Branch 2 not taken.
38129 Sign o = Sign(orient2d(i,j,v0) * orient_012_);
629
2/2
✓ Branch 0 taken 18771 times.
✓ Branch 1 taken 19358 times.
38129 if(t2v0_t1v2) {
630
1/2
✓ Branch 1 taken 18771 times.
✗ Branch 2 not taken.
18771 swap_edge(t1,false); // "new t1 on top"
631
2/2
✓ Branch 0 taken 9573 times.
✓ Branch 1 taken 9198 times.
18771 if(o >= 0) {
632
2/8
✓ Branch 1 taken 9573 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9573 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
9573 geo_debug_assert(!segment_edge_intersect(i,j,t1,2));
633
2/8
✓ Branch 1 taken 9573 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9573 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
9573 geo_debug_assert( segment_edge_intersect(i,j,t2,0));
634
1/2
✓ Branch 1 taken 9573 times.
✗ Branch 2 not taken.
9573 new_edge(t1,2);
635 } else {
636
2/8
✓ Branch 1 taken 9198 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9198 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
9198 geo_debug_assert( segment_edge_intersect(i,j,t1,2));
637
2/8
✓ Branch 1 taken 9198 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9198 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
9198 geo_debug_assert( segment_edge_intersect(i,j,t2,0));
638
1/2
✓ Branch 1 taken 9198 times.
✗ Branch 2 not taken.
9198 isect_edge(t1,2);
639 }
640 } else {
641
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 19358 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
19358 geo_debug_assert(t2v0_t1v1);
642
1/2
✓ Branch 1 taken 19358 times.
✗ Branch 2 not taken.
19358 swap_edge(t1,true); // "new t1 on bottom"
643
2/2
✓ Branch 0 taken 8956 times.
✓ Branch 1 taken 10402 times.
19358 if(o > 0) {
644
2/8
✓ Branch 1 taken 8956 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 8956 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
8956 geo_debug_assert( segment_edge_intersect(i,j,t1,1));
645
2/8
✓ Branch 1 taken 8956 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 8956 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
8956 geo_debug_assert( segment_edge_intersect(i,j,t2,0));
646
1/2
✓ Branch 1 taken 8956 times.
✗ Branch 2 not taken.
8956 isect_edge(t1,1);
647 } else {
648
2/8
✓ Branch 1 taken 10402 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 10402 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
10402 geo_debug_assert(!segment_edge_intersect(i,j,t1,1));
649
2/8
✓ Branch 1 taken 10402 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 10402 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
10402 geo_debug_assert( segment_edge_intersect(i,j,t2,0));
650
1/2
✓ Branch 1 taken 10402 times.
✗ Branch 2 not taken.
10402 new_edge(t1,1);
651 }
652 }
653 }
654 }
655 }
656 258470 }
657
658 252049 void CDTBase2d::Delaunayize_vertex_neighbors(index_t v, DList& S) {
659 CDT_LOG("Delaunayize_vertex_neighbors");
660 252049 index_t count = 0;
661
2/2
✓ Branch 1 taken 1702083 times.
✓ Branch 2 taken 252049 times.
1954132 while(!S.empty()) {
662 // NASA programming style: all loops have
663 // a maximum number of iterations
664 1702083 ++count;
665
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1702083 times.
1702083 if(count > 10*nT()) {
666 Logger::warn("CDT2d")
667 << "Emergency exit in Delaunayize_vertex_neighbors()"
668 << std::endl;
669 S.clear();
670 geo_assert_not_reached; // For now, assert fail.
671 }
672 1702083 index_t t1 = S.pop_back();
673
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 1702083 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
1702083 geo_debug_assert(Tv(t1,0) == v);
674
2/2
✓ Branch 1 taken 9636 times.
✓ Branch 2 taken 1692447 times.
1702083 if(Tedge_is_constrained(t1,0)) {
675 9636 continue;
676 }
677 1692447 index_t t2 = Tadj(t1,0);
678
2/2
✓ Branch 0 taken 203876 times.
✓ Branch 1 taken 1488571 times.
1692447 if(t2 == NO_INDEX) {
679 203876 continue;
680 }
681
5/6
✓ Branch 0 taken 1488571 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 229919 times.
✓ Branch 4 taken 1258652 times.
✓ Branch 5 taken 229919 times.
✓ Branch 6 taken 1258652 times.
1488571 if(!exact_incircle_ && !is_convex_quad(t1)) {
682 229919 continue;
683 }
684
685 1258652 index_t v1 = Tv(t2,0);
686 1258652 index_t v2 = Tv(t2,1);
687 1258652 index_t v3 = Tv(t2,2);
688
2/2
✓ Branch 1 taken 504852 times.
✓ Branch 2 taken 753800 times.
1258652 if(Sign(incircle(v1,v2,v3,v)*orient_012_) == POSITIVE) {
689 504852 swap_edge(t1);
690 504852 S.push_back(t1);
691 504852 S.push_back(t2);
692 }
693 }
694 CDT_LOG("/Delaunayize_vertex_neighbors");
695 252049 }
696
697 257957 void CDTBase2d::Delaunayize_new_edges(DList& N) {
698 257957 index_t count = 0;
699 257957 bool swap_occured = true;
700
2/2
✓ Branch 0 taken 259966 times.
✓ Branch 1 taken 257957 times.
517923 while(swap_occured) {
701 259966 swap_occured = false;
702 // NASA programming style: all loops have
703 // a maximum number of iterations
704 259966 ++count;
705
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 259966 times.
259966 if(count > 10*nT()) {
706 Logger::warn("CDT2d")
707 << "Emergency exit in Delaunayize_new_edges()"
708 << std::endl;
709 break;
710 }
711
2/2
✓ Branch 2 taken 23453 times.
✓ Branch 3 taken 259966 times.
283419 for(index_t t1 = N.front(); t1 != NO_INDEX; t1 = N.next(t1)) {
712
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 23453 times.
23453 if(Tedge_is_constrained(t1,0)) {
713 continue;
714 }
715 23453 index_t v1 = Tv(t1,1);
716 23453 index_t v2 = Tv(t1,2);
717 23453 index_t v0 = Tv(t1,0);
718 23453 index_t t2 = Tadj(t1,0);
719
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 23453 times.
23453 if(t2 == NO_INDEX) {
720 continue;
721 }
722
5/6
✓ Branch 0 taken 23453 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 12741 times.
✓ Branch 4 taken 10712 times.
✓ Branch 5 taken 12741 times.
✓ Branch 6 taken 10712 times.
23453 if(!exact_incircle_ && !is_convex_quad(t1)) {
723 12741 continue;
724 }
725 10712 index_t e2 = Tadj_find(t2,t1);
726 10712 index_t v3 = Tv(t2,e2);
727
2/2
✓ Branch 1 taken 2949 times.
✓ Branch 2 taken 7763 times.
10712 if(Sign(incircle(v0,v1,v2,v3)*orient_012_) == POSITIVE) {
728 // t2 may also encode a new edge, we need to preserve it,
729 // by chosing the right swap:
730
2/2
✓ Branch 2 taken 2431 times.
✓ Branch 3 taken 518 times.
2949 if(Tv(t2,0) == Tv(t1,1)) {
731 2431 swap_edge(t1, true); // t2 on top
732 2431 Trot(t1,1);
733 } else {
734 518 swap_edge(t1, false); // t1 on top
735 518 Trot(t1,2);
736 }
737 2949 swap_occured = true;
738 }
739 }
740 }
741 257957 N.clear();
742 257957 }
743
744 515221 index_t CDTBase2d::locate(index_t v, index_t hint, Sign* o) const {
745 Sign o_local[3];
746
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 515221 times.
515221 if(o == nullptr) {
747 o = o_local;
748 }
749
750 // Efficient locate, "walking the triangulation"
751
1/2
✓ Branch 1 taken 515221 times.
✗ Branch 2 not taken.
515221 index_t t_pred = nT()+1; // Needs to be different from NO_INDEX
752
2/2
✓ Branch 0 taken 292535 times.
✓ Branch 1 taken 222686 times.
515221 index_t t = (hint == NO_INDEX) ?
753
2/4
✓ Branch 1 taken 292535 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 292535 times.
✗ Branch 5 not taken.
292535 index_t(Numeric::random_int32()) % nT() :
754 515221 hint ;
755 #ifdef GEO_DEBUG
756 515221 index_t nb_traversed_t = 0;
757 #endif
758
759 6120976 still_walking:
760 {
761 #ifdef GEO_DEBUG
762 6120976 ++nb_traversed_t;
763 #endif
764
765 // Infinite loop are not supposed to happen, but
766 // let us detect them, just in case...
767
2/8
✓ Branch 1 taken 6120976 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 6120976 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
6120976 geo_debug_assert(nb_traversed_t <= 2*nT());
768
769 // You will land here if we try to locate a point outside
770 // the boundary
771 6120976 bool point_outside_boundary = (t == NO_INDEX);
772
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 6120976 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
6120976 geo_assert(!point_outside_boundary);
773
774 index_t tv[3];
775
1/2
✓ Branch 1 taken 6120976 times.
✗ Branch 2 not taken.
6120976 tv[0] = Tv(t,0);
776
1/2
✓ Branch 1 taken 6120976 times.
✗ Branch 2 not taken.
6120976 tv[1] = Tv(t,1);
777
1/2
✓ Branch 1 taken 6120976 times.
✗ Branch 2 not taken.
6120976 tv[2] = Tv(t,2);
778
779 // Start from a random edge
780
1/2
✓ Branch 1 taken 6120976 times.
✗ Branch 2 not taken.
6120976 index_t e0 = index_t(Numeric::random_int32()) % 3;
781
2/2
✓ Branch 0 taken 11877849 times.
✓ Branch 1 taken 515221 times.
12393070 for(index_t de = 0; de < 3; ++de) {
782 11877849 index_t le = (e0 + de) % 3;
783
784
1/2
✓ Branch 1 taken 11877849 times.
✗ Branch 2 not taken.
11877849 index_t t_next = Tadj(t,le);
785
786 // If the candidate next triangle is the
787 // one we came from, then we know already that
788 // the orientation is positive, thus we examine
789 // the next candidate (or exit the loop if they
790 // are exhausted).
791 //
792 // (here is why intial value of t_pred needs to be
793 // different from NO_INDEX)
794
2/2
✓ Branch 0 taken 2847448 times.
✓ Branch 1 taken 9030401 times.
11877849 if(t_next == t_pred) {
795 2847448 o[le] = POSITIVE;
796 2847448 continue ;
797 }
798
799 // To test the orientation of p w.r.t. the facet f of
800 // t, we replace vertex number f with p in t (same
801 // convention as in CGAL).
802 9030401 index_t v_bkp = tv[le];
803 9030401 tv[le] = v;
804
1/2
✓ Branch 1 taken 9030401 times.
✗ Branch 2 not taken.
9030401 o[le] = Sign(orient_012_ * orient2d(tv[0], tv[1], tv[2]));
805
806 // If the orientation is not negative, then we cannot
807 // walk towards t_next, and examine the next candidate
808 // (or exit the loop if they are exhausted).
809
2/2
✓ Branch 0 taken 3424646 times.
✓ Branch 1 taken 5605755 times.
9030401 if(o[le] != NEGATIVE) {
810 3424646 tv[le] = v_bkp;
811 3424646 continue;
812 }
813
814 // If we reach this point, then t_next is a valid
815 // successor, thus we are still walking.
816 5605755 t_pred = t;
817 5605755 t = t_next;
818 5605755 goto still_walking;
819 }
820 }
821
822 515221 return t;
823 }
824
825 8071 void CDTBase2d::remove_external_triangles(bool remove_internal_holes) {
826
827
2/2
✓ Branch 0 taken 8067 times.
✓ Branch 1 taken 4 times.
8071 if(remove_internal_holes) {
828
1/2
✓ Branch 1 taken 8067 times.
✗ Branch 2 not taken.
8067 DList S(*this, DLIST_S_ID);
829
830 // Step 1: get triangles adjacent to the border,
831 // mark them as visited, classify them
832
3/4
✓ Branch 1 taken 111149 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 103082 times.
✓ Branch 4 taken 8067 times.
111149 for(index_t t=0; t<nT(); ++t) {
833
2/2
✓ Branch 0 taken 244798 times.
✓ Branch 1 taken 70814 times.
315612 for(index_t le=0; le<3; ++le) {
834
3/4
✓ Branch 1 taken 244798 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 32268 times.
✓ Branch 4 taken 212530 times.
244798 if(Tadj(t,le) == NO_INDEX) {
835
1/2
✓ Branch 1 taken 32268 times.
✗ Branch 2 not taken.
32268 bool outside = ((Tedge_cnstr_nb(t,le)%2) == 0);
836
1/2
✓ Branch 1 taken 32268 times.
✗ Branch 2 not taken.
32268 Tset_flag(t, T_VISITED_FLAG);
837
1/2
✓ Branch 0 taken 32268 times.
✗ Branch 1 not taken.
32268 if(outside) {
838
1/2
✓ Branch 1 taken 32268 times.
✗ Branch 2 not taken.
32268 Tset_flag(t, T_MARKED_FLAG);
839 }
840
1/2
✓ Branch 1 taken 32268 times.
✗ Branch 2 not taken.
32268 S.push_back(t);
841 32268 break;
842 }
843 }
844 }
845
846 // Step 2: recursive traversal
847
3/4
✓ Branch 1 taken 111149 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 103082 times.
✓ Branch 4 taken 8067 times.
111149 while(!S.empty()) {
848
1/2
✓ Branch 1 taken 103082 times.
✗ Branch 2 not taken.
103082 index_t t1 = S.pop_back();
849
1/2
✓ Branch 1 taken 103082 times.
✗ Branch 2 not taken.
103082 bool t1_outside = Tflag_is_set(t1, T_MARKED_FLAG);
850
2/2
✓ Branch 0 taken 309246 times.
✓ Branch 1 taken 103082 times.
412328 for(index_t le=0; le<3; ++le) {
851
1/2
✓ Branch 1 taken 309246 times.
✗ Branch 2 not taken.
309246 index_t t2 = Tadj(t1,le);
852 309246 if(
853
4/4
✓ Branch 0 taken 276978 times.
✓ Branch 1 taken 32268 times.
✓ Branch 2 taken 70814 times.
✓ Branch 3 taken 238432 times.
586224 t2 != NO_INDEX &&
854
3/4
✓ Branch 1 taken 276978 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 70814 times.
✓ Branch 4 taken 206164 times.
276978 !Tflag_is_set(t2,T_VISITED_FLAG)
855 ) {
856 bool t2_outside =
857
1/2
✓ Branch 1 taken 70814 times.
✗ Branch 2 not taken.
70814 t1_outside ^ ((Tedge_cnstr_nb(t1,le)%2) != 0);
858
1/2
✓ Branch 1 taken 70814 times.
✗ Branch 2 not taken.
70814 Tset_flag(t2, T_VISITED_FLAG);
859
2/2
✓ Branch 0 taken 43270 times.
✓ Branch 1 taken 27544 times.
70814 if(t2_outside) {
860
1/2
✓ Branch 1 taken 43270 times.
✗ Branch 2 not taken.
43270 Tset_flag(t2, T_MARKED_FLAG);
861 }
862
1/2
✓ Branch 1 taken 70814 times.
✗ Branch 2 not taken.
70814 S.push_back(t2);
863 }
864 }
865 }
866
867 // Step 3: reset visited flag
868
3/4
✓ Branch 1 taken 111149 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 103082 times.
✓ Branch 4 taken 8067 times.
111149 for(index_t t=0; t<nT(); ++t) {
869
1/2
✓ Branch 1 taken 103082 times.
✗ Branch 2 not taken.
103082 Treset_flag(t, T_VISITED_FLAG);
870 }
871
872 8067 } else {
873
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
4 DList S(*this, DLIST_S_ID);
874
875 // Step 1: get triangles adjacent to the border
876
3/4
✓ Branch 1 taken 40656 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 40652 times.
✓ Branch 4 taken 4 times.
40656 for(index_t t=0; t<nT(); ++t) {
877 40652 if(
878 // TODO: replace with parity check ?
879
3/6
✓ Branch 1 taken 40652 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 40396 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 40396 times.
✗ Branch 7 not taken.
81048 (!Tedge_is_constrained(t,0) && Tadj(t,0) == NO_INDEX) ||
880
8/12
✓ Branch 0 taken 40396 times.
✓ Branch 1 taken 256 times.
✓ Branch 3 taken 40652 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 40295 times.
✓ Branch 6 taken 357 times.
✓ Branch 8 taken 40295 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 40295 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 40652 times.
121700 (!Tedge_is_constrained(t,1) && Tadj(t,1) == NO_INDEX) ||
881
5/8
✓ Branch 1 taken 40652 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 40445 times.
✓ Branch 4 taken 207 times.
✓ Branch 6 taken 40445 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 40445 times.
40652 (!Tedge_is_constrained(t,2) && Tadj(t,2) == NO_INDEX)
882 ) {
883 Tset_flag(t, T_MARKED_FLAG);
884 S.push_back(t);
885 }
886 }
887
888 // Step 2: recursive traversal
889
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
4 while(!S.empty()) {
890 index_t t1 = S.pop_back();
891 for(index_t le=0; le<3; ++le) {
892 index_t t2 = Tadj(t1,le);
893 if(
894 t2 != NO_INDEX &&
895 !Tedge_is_constrained(t1,le) &&
896 !Tflag_is_set(t2,T_MARKED_FLAG)
897 ) {
898 Tset_flag(t2, T_MARKED_FLAG);
899 S.push_back(t2);
900 }
901 }
902 }
903 4 }
904
905 // Step 3: remove marked triangles
906 8071 remove_marked_triangles();
907 8071 }
908
909
910 8177 void CDTBase2d::remove_marked_triangles() {
911 // Step 1: compute old2new map
912 // (use Tnext_'s storage, that we do not need now)
913 8177 vector<index_t>& old2new = Tnext_;
914 8177 index_t cur_t_new = 0;
915
2/2
✓ Branch 1 taken 150268 times.
✓ Branch 2 taken 8177 times.
158445 for(index_t t=0; t<nT(); ++t) {
916
2/2
✓ Branch 1 taken 79333 times.
✓ Branch 2 taken 70935 times.
150268 if(Tflag_is_set(t,T_MARKED_FLAG)) {
917 79333 old2new[t] = NO_INDEX;
918 } else {
919 70935 old2new[t] = cur_t_new;
920 70935 ++cur_t_new;
921 }
922 }
923 8177 index_t nT_new = cur_t_new;
924
925 // Step 2: translate adjacency and move triangles
926
2/2
✓ Branch 1 taken 150268 times.
✓ Branch 2 taken 8177 times.
158445 for(index_t t=0; t<nT(); ++t) {
927 150268 index_t t_new = old2new[t];
928
2/2
✓ Branch 0 taken 79333 times.
✓ Branch 1 taken 70935 times.
150268 if(t_new == NO_INDEX) {
929 79333 continue;
930 }
931 70935 index_t adj0 = Tadj(t,0);
932
2/2
✓ Branch 0 taken 70923 times.
✓ Branch 1 taken 12 times.
70935 if(adj0 != NO_INDEX) {
933 70923 adj0 = old2new[adj0];
934 }
935 70935 index_t adj1 = Tadj(t,1);
936
1/2
✓ Branch 0 taken 70935 times.
✗ Branch 1 not taken.
70935 if(adj1 != NO_INDEX) {
937 70935 adj1 = old2new[adj1];
938 }
939 70935 index_t adj2 = Tadj(t,2);
940
1/2
✓ Branch 0 taken 70935 times.
✗ Branch 1 not taken.
70935 if(adj2 != NO_INDEX) {
941 70935 adj2 = old2new[adj2];
942 }
943 70935 Tset(
944 t_new,
945 Tv(t,0), Tv(t,1), Tv(t,2),
946 adj0, adj1, adj2,
947 Tedge_cnstr_first(t,0),
948 Tedge_cnstr_first(t,1),
949 Tedge_cnstr_first(t,2)
950 );
951 70935 Tflags_[t_new] = 0;
952 }
953
954 // Step 3: resize arrays
955 8177 T_.resize(3*nT_new);
956 8177 Tadj_.resize(3*nT_new);
957 8177 Tflags_.resize(nT_new);
958 8177 Tecnstr_first_.resize(3*nT_new);
959 8177 Tnext_.resize(nT_new);
960 8177 Tprev_.resize(nT_new);
961
962 // Step 4: fix v2T_
963
2/2
✓ Branch 1 taken 70935 times.
✓ Branch 2 taken 8177 times.
79112 for(index_t t=0; t<nT(); ++t) {
964 70935 v2T_[Tv(t,0)] = t;
965 70935 v2T_[Tv(t,1)] = t;
966 70935 v2T_[Tv(t,2)] = t;
967 }
968 8177 }
969
970
971 /***************** Triangulation surgery (boring code ahead) *********/
972
973 80560 void CDTBase2d::insert_vertex_in_edge(
974 index_t v, index_t t, index_t le1, DList& S
975 ) {
976 80560 index_t cnstr_first = Tedge_cnstr_first(t,le1);
977 80560 index_t t1 = t;
978 80560 index_t t2 = Tadj(t1,le1);
979 80560 index_t v1 = Tv(t1,le1);
980 80560 index_t v2 = Tv(t1,(le1+1)%3);
981 80560 index_t v3 = Tv(t1,(le1+2)%3);
982 80560 index_t t1_adj2 = Tadj(t1,(le1+1)%3);
983 80560 index_t t1_adj3 = Tadj(t1,(le1+2)%3);
984
2/2
✓ Branch 0 taken 8112 times.
✓ Branch 1 taken 72448 times.
80560 if(t2 != NO_INDEX) {
985 CDT_LOG(" insert vertex in internal edge");
986 // New vertex is on an edge of t1 and t1 has a neighbor
987 // accross that edge. Discard the two triangles t1 and t2
988 // adjacent to the edge, and create four new triangles
989 // (t1 and t2 are recycled).
990 8112 index_t le2 = Tadj_find(t2,t1);
991
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 8112 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
8112 geo_debug_assert(Tv(t2, (le2+1)%3) == v3);
992
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 8112 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
8112 geo_debug_assert(Tv(t2, (le2+2)%3) == v2);
993 8112 index_t v4 = Tv(t2,le2);
994 8112 index_t t2_adj2 = Tadj(t2,(le2+1)%3);
995 8112 index_t t2_adj3 = Tadj(t2,(le2+2)%3);
996 8112 index_t t3 = Tnew();
997 8112 index_t t4 = Tnew();
998 8112 Tset(t1,v,v1,v2,t1_adj3,t2,t4);
999 8112 Tset(t2,v,v2,v4,t2_adj2,t3,t1);
1000 8112 Tset(t3,v,v4,v3,t2_adj3,t4,t2);
1001 8112 Tset(t4,v,v3,v1,t1_adj2,t1,t3);
1002 8112 Tadj_back_connect(t1,0,t1);
1003 8112 Tadj_back_connect(t2,0,t2);
1004 8112 Tadj_back_connect(t3,0,t2);
1005 8112 Tadj_back_connect(t4,0,t1);
1006 8112 Tset_edge_cnstr_first(t1,1,cnstr_first);
1007 8112 Tset_edge_cnstr_first(t2,2,cnstr_first);
1008 8112 Tset_edge_cnstr_first(t3,1,cnstr_first);
1009 8112 Tset_edge_cnstr_first(t4,2,cnstr_first);
1010
2/2
✓ Branch 1 taken 2178 times.
✓ Branch 2 taken 5934 times.
8112 if(S.initialized()) {
1011 2178 S.push_back(t1);
1012 2178 S.push_back(t2);
1013 2178 S.push_back(t3);
1014 2178 S.push_back(t4);
1015 }
1016 } else {
1017 CDT_LOG(" insert vertex in border edge");
1018 // New vertex is on an edge of t1 and t1 has no neighbor
1019 // accross that edge. Discard t1 and replace it with two
1020 // new triangles (recycle t1).
1021 72448 t2 = Tnew();
1022 72448 Tset(t1,v,v1,v2,t1_adj3,NO_INDEX,t2);
1023 72448 Tset(t2,v,v3,v1,t1_adj2,t1,NO_INDEX);
1024 72448 Tadj_back_connect(t1,0,t1);
1025 72448 Tadj_back_connect(t2,0,t1);
1026 72448 Tset_edge_cnstr_first(t1,1,cnstr_first);
1027 72448 Tset_edge_cnstr_first(t2,2,cnstr_first);
1028
2/2
✓ Branch 1 taken 72257 times.
✓ Branch 2 taken 191 times.
72448 if(S.initialized()) {
1029 72257 S.push_back(t1);
1030 72257 S.push_back(t2);
1031 }
1032 }
1033 80560 }
1034
1035 282929 void CDTBase2d::insert_vertex_in_triangle(index_t v, index_t t, DList& S) {
1036 // New vertex is in t1. Discard t1 and replace it with three
1037 // new triangles (recycle t1).
1038 282929 index_t t1 = t;
1039 282929 index_t v1 = Tv(t1,0);
1040 282929 index_t v2 = Tv(t1,1);
1041 282929 index_t v3 = Tv(t1,2);
1042 282929 index_t adj1 = Tadj(t1,0);
1043 282929 index_t adj2 = Tadj(t1,1);
1044 282929 index_t adj3 = Tadj(t1,2);
1045 282929 index_t t2 = Tnew();
1046 282929 index_t t3 = Tnew();
1047 282929 Tset(t1,v,v2,v3,adj1,t2,t3);
1048 282929 Tset(t2,v,v3,v1,adj2,t3,t1);
1049 282929 Tset(t3,v,v1,v2,adj3,t1,t2);
1050 282929 Tadj_back_connect(t1,0,t1);
1051 282929 Tadj_back_connect(t2,0,t1);
1052 282929 Tadj_back_connect(t3,0,t1);
1053
2/2
✓ Branch 1 taken 171765 times.
✓ Branch 2 taken 111164 times.
282929 if(S.initialized()) {
1054 171765 S.push_back(t1);
1055 171765 S.push_back(t2);
1056 171765 S.push_back(t3);
1057 }
1058 282929 }
1059
1060 570931 void CDTBase2d::swap_edge(index_t t1, bool swap_t1_t2) {
1061
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 570931 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
570931 geo_debug_assert(!Tedge_is_constrained(t1,0));
1062 570931 index_t v1 = Tv(t1,0);
1063 570931 index_t v2 = Tv(t1,1);
1064 570931 index_t v3 = Tv(t1,2);
1065 570931 index_t t1_adj2 = Tadj(t1,1);
1066 570931 index_t t1_adj3 = Tadj(t1,2);
1067 570931 index_t t2 = Tadj(t1,0);
1068 570931 index_t le2 = Tadj_find(t2,t1);
1069 570931 index_t v4 = Tv(t2,le2);
1070
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 570931 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
570931 geo_debug_assert(Tv(t2,(le2+1)%3) == v3);
1071
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 570931 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
570931 geo_debug_assert(Tv(t2,(le2+2)%3) == v2);
1072
1073 570931 debug_Tcheck(t1);
1074 570931 debug_Tcheck(t2);
1075
1076 570931 index_t t2_adj2 = Tadj(t2,(le2+1)%3);
1077 570931 index_t t2_adj3 = Tadj(t2,(le2+2)%3);
1078
2/2
✓ Branch 0 taken 21789 times.
✓ Branch 1 taken 549142 times.
570931 if(swap_t1_t2) {
1079 21789 Tset(t2,v1,v4,v3,t2_adj3,t1_adj2,t1);
1080 21789 Tset(t1,v1,v2,v4,t2_adj2,t2,t1_adj3);
1081 21789 Tadj_back_connect(t2,0,t2);
1082 21789 Tadj_back_connect(t2,1,t1);
1083 21789 Tadj_back_connect(t1,0,t2);
1084 21789 Tadj_back_connect(t1,2,t1);
1085 } else {
1086 549142 Tset(t1,v1,v4,v3,t2_adj3,t1_adj2,t2);
1087 549142 Tset(t2,v1,v2,v4,t2_adj2,t1,t1_adj3);
1088 549142 Tadj_back_connect(t1,0,t2);
1089 549142 Tadj_back_connect(t1,1,t1);
1090 549142 Tadj_back_connect(t2,0,t2);
1091 549142 Tadj_back_connect(t2,2,t1);
1092 }
1093
1094 570931 debug_Tcheck(t1);
1095 570931 debug_Tcheck(t2);
1096 570931 }
1097
1098 /***************** Geometry ***********************/
1099
1100 1598509 bool CDTBase2d::is_convex_quad(index_t t) const {
1101 1598509 index_t v1 = Tv(t,0);
1102 1598509 index_t v2 = Tv(t,1);
1103 1598509 index_t v3 = Tv(t,2);
1104 1598509 index_t t2 = Tadj(t,0);
1105 1598509 index_t le2 = Tadj_find(t2,t);
1106 1598509 index_t v4 = Tv(t2,le2);
1107 // t and Tadj(t,0) have the correct orientation,
1108 // so one just needs to check the orientation of
1109 // the two triangles that would be generated by
1110 // an edge flip.
1111 return
1112
2/2
✓ Branch 1 taken 1470901 times.
✓ Branch 2 taken 127608 times.
3069410 orient2d(v1,v4,v3) == orient_012_ &&
1113
2/2
✓ Branch 1 taken 1334537 times.
✓ Branch 2 taken 136364 times.
3069410 orient2d(v4,v1,v2) == orient_012_ ;
1114 }
1115
1116 /******* Debugging ******************************************************/
1117
1118 93330 void CDTBase2d::check_geometry() const {
1119
3/4
✓ Branch 0 taken 93266 times.
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 93266 times.
93330 if(delaunay_ && exact_incircle_) {
1120 for(index_t t=0; t<nT(); ++t) {
1121 for(index_t le=0; le<3; ++le) {
1122 // geo_assert(Tedge_is_Delaunay(t,le));
1123 }
1124 }
1125 }
1126 93330 }
1127
1128
1129 bool CDTBase2d::Tedge_is_Delaunay(index_t t1, index_t le1) const {
1130 if(Tedge_is_constrained(t1,le1)) {
1131 return true;
1132 }
1133 index_t t2 = Tadj(t1,le1);
1134 if(t2 == NO_INDEX) {
1135 return true;
1136 }
1137 index_t le2 = Tadj_find(t2,t1);
1138 index_t v1 = Tv(t1,le1);
1139 index_t v2 = Tv(t1,(le1+1)%3);
1140 index_t v3 = Tv(t1,(le1+2)%3);
1141 index_t v4 = Tv(t2,le2);
1142
1143 // If we do not check that, we assert fail
1144 // whenever there is a vertex inserted in
1145 // the macroborder of the triangle
1146 if(
1147 orient2d(v1,v4,v3) != orient_012_ ||
1148 orient2d(v4,v1,v2) != orient_012_
1149 ) {
1150 return true;
1151 }
1152
1153 return Sign(incircle(v1,v2,v3,v4)*orient_012_) <= 0;
1154 }
1155
1156 void CDTBase2d::check_edge_intersections(
1157 index_t v1, index_t v2, const DList& Q
1158 ) {
1159 std::set<Edge> I;
1160 auto make_edge = [](index_t w1, index_t w2)->Edge {
1161 return std::make_pair(std::min(w1,w2), std::max(w1,w2));
1162 };
1163 for(index_t t=Q.front(); t!=NO_INDEX; t = Q.next(t)) {
1164 geo_assert(segment_edge_intersect(v1,v2,t,0));
1165 I.insert(make_edge(Tv(t,1), Tv(t,2)));
1166 }
1167 for(index_t t=0; t<nT(); ++t) {
1168 for(index_t le=0; le<3; ++le) {
1169 if(segment_edge_intersect(v1,v2,t,le)) {
1170 index_t w1 = Tv(t,(le+1)%3);
1171 index_t w2 = Tv(t,(le+2)%3);
1172 geo_assert(I.find(make_edge(w1,w2)) != I.end());
1173 }
1174 }
1175 }
1176 }
1177
1178 /*** Naive versions of algorithms, for reference / debugging if need be ***/
1179
1180 index_t CDTBase2d::locate_naive(index_t v, index_t hint, Sign* o) const {
1181 geo_argused(hint);
1182 Sign o_local[3];
1183 if(o == nullptr) {
1184 o = o_local;
1185 }
1186
1187 for(index_t t=0; t<nT(); ++t) {
1188 index_t i = Tv(t,0);
1189 index_t j = Tv(t,1);
1190 index_t k = Tv(t,2);
1191 o[0] = orient2d(v,j,k);
1192 o[1] = orient2d(v,k,i);
1193 o[2] = orient2d(v,i,j);
1194 if(o[0]*o[1] >= 0 && o[1]*o[2] >= 0 && o[2]*o[0] >= 0) {
1195 return t;
1196 }
1197 }
1198 geo_assert_not_reached;
1199 }
1200
1201 void CDTBase2d::Delaunayize_new_edges_naive(vector<Edge>& N) {
1202 CDT_LOG("Delaunayize_new_edges_naive()");
1203 for(Edge E: N) {
1204 index_t v1 = std::min(E.first,E.second);
1205 index_t v2 = std::max(E.first,E.second);
1206 if(v2 < v1) {
1207 std::swap(v1,v2);
1208 }
1209 CDT_LOG("new edge: " << v1 << " " << v2);
1210 }
1211 index_t count = 0;
1212 bool swap_occured = true;
1213 while(swap_occured) {
1214 // NASA programming style: all loops have
1215 // a maximum number of iterations
1216 ++count;
1217 if(count > 10*nT()) {
1218 Logger::warn("CDT2d")
1219 << "Emergency exit in Delaunayize_new_edges_naive()"
1220 << std::endl;
1221 return;
1222 }
1223 swap_occured = false;
1224 for(Edge& E: N) {
1225 index_t t1 = eT(E);
1226 if(Tedge_is_constrained(t1,0)) {
1227 continue;
1228 }
1229 index_t v1 = Tv(t1,1);
1230 index_t v2 = Tv(t1,2);
1231 index_t v0 = Tv(t1,0);
1232 index_t t2 = Tadj(t1,0);
1233 if(t2 == NO_INDEX) {
1234 continue;
1235 }
1236 index_t e2 = Tadj_find(t2,t1);
1237 index_t v3 = Tv(t2,e2);
1238 if(!exact_incircle_ && !is_convex_quad(t1)) {
1239 continue;
1240 }
1241 if(Sign(incircle(v0,v1,v2,v3)*orient_012_) == POSITIVE) {
1242 CDT_LOG("swap " << v1 << " " << v2
1243 << " ---> "
1244 << v0 << " " << v3 );
1245 swap_edge(t1);
1246 E = std::make_pair(Tv(t1,0), Tv(t1,1));
1247 swap_occured = true;
1248 }
1249 }
1250 }
1251 N.resize(0);
1252 CDT_LOG("/Delaunayize_new_edges_naive()");
1253 }
1254
1255 void CDTBase2d::constrain_edges_naive(
1256 index_t i, index_t j, DList& Q_in, vector<Edge>& N
1257 ) {
1258 CDT_LOG("Q size=" << Q_in.size());
1259
1260 std::deque<Edge> Q;
1261 for(index_t t=Q_in.front(); t != NO_INDEX; t = Q_in.next(t)) {
1262 Q.push_back(std::make_pair(Tv(t,1), Tv(t,2)));
1263 }
1264 Q_in.clear();
1265
1266 for(index_t t=0; t<nT(); ++t) {
1267 geo_debug_assert(!Tis_in_list(t));
1268 }
1269
1270 while(Q.size() != 0) {
1271 Edge E = Q.back();
1272 Q.pop_back();
1273 if(!is_convex_quad(eT(E))) {
1274 if(Q.size() == 0) {
1275 CDT_LOG("... infinite iteration");
1276 abort();
1277 }
1278 Q.push_front(E);
1279 } else {
1280 index_t t = eT(E);
1281 swap_edge(t);
1282 E = std::make_pair(Tv(t,0), Tv(t,1));
1283 if(segment_segment_intersect(i,j,E.first,E.second)) {
1284 Q.push_front(E);
1285 } else {
1286 if(
1287 (E.first == i && E.second == j) ||
1288 (E.first == j && E.second == i)
1289 ) {
1290 Tadd_edge_cnstr_with_neighbor(eT(E),0,ncnstr_-1);
1291 } else {
1292 N.push_back(E);
1293 }
1294 }
1295 }
1296 }
1297 }
1298
1299 /********************************************************************/
1300
1301 32 CDT2d::CDT2d() {
1302 32 exact_intersections_ = false;
1303 32 exact_incircle_ = false;
1304 32 }
1305
1306 64 CDT2d::~CDT2d() {
1307 64 }
1308
1309 void CDT2d::clear() {
1310 CDTBase2d::clear();
1311 point_.resize(0);
1312 }
1313
1314 32 void CDT2d::create_enclosing_triangle(
1315 const vec2& p1, const vec2& p2, const vec2& p3
1316 ) {
1317
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 32 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
32 geo_assert(nv() == 0);
1318
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 32 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
32 geo_assert(nT() == 0);
1319 32 point_.push_back(p1);
1320 32 point_.push_back(p2);
1321 32 point_.push_back(p3);
1322 32 CDTBase2d::create_enclosing_triangle(0,1,2);
1323 32 }
1324
1325 void CDT2d::create_enclosing_quad(
1326 const vec2& p1, const vec2& p2, const vec2& p3, const vec2& p4
1327 ) {
1328 geo_assert(nv() == 0);
1329 geo_assert(nT() == 0);
1330 point_.push_back(p1);
1331 point_.push_back(p2);
1332 point_.push_back(p3);
1333 point_.push_back(p4);
1334 CDTBase2d::create_enclosing_quad(0,1,2,3);
1335 }
1336
1337 9081215 Sign CDT2d::orient2d(index_t i, index_t j, index_t k) const {
1338
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 9081215 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
9081215 geo_debug_assert(i < nv());
1339
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 9081215 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
9081215 geo_debug_assert(j < nv());
1340
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 9081215 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
9081215 geo_debug_assert(k < nv());
1341 9081215 return PCK::orient_2d(
1342 9081215 point_[i].data(), point_[j].data(), point_[k].data()
1343 9081215 );
1344 }
1345
1346 935458 Sign CDT2d::incircle(index_t i, index_t j, index_t k, index_t l) const {
1347
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 935458 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
935458 geo_debug_assert(i < nv());
1348
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 935458 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
935458 geo_debug_assert(j < nv());
1349
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 935458 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
935458 geo_debug_assert(k < nv());
1350
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 935458 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
935458 geo_debug_assert(l < nv());
1351 935458 return PCK::in_circle_2d_SOS(
1352 935458 point_[i].data(), point_[j].data(), point_[k].data(),
1353 935458 point_[l].data()
1354 935458 );
1355 }
1356
1357 168 index_t CDT2d::create_intersection(
1358 index_t E1, index_t i, index_t j,
1359 index_t E2, index_t k, index_t l
1360 ) {
1361 168 geo_argused(E1);
1362 168 geo_argused(E2);
1363
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 168 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
168 geo_debug_assert(i < nv());
1364
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 168 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
168 geo_debug_assert(j < nv());
1365
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 168 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
168 geo_debug_assert(k < nv());
1366
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 168 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
168 geo_debug_assert(l < nv());
1367
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 168 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
168 geo_debug_assert(E1 < ncnstr());
1368
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 168 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
168 geo_debug_assert(E2 < ncnstr());
1369
2/4
✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 168 times.
✗ Branch 5 not taken.
168 vec2 U = point_[j] - point_[i];
1370
2/4
✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 168 times.
✗ Branch 5 not taken.
168 vec2 V = point_[l] - point_[k];
1371
2/4
✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 168 times.
✗ Branch 5 not taken.
168 vec2 D = point_[k] - point_[i];
1372 168 double delta = det(U,V);
1373 168 double t = det(D,V)/delta;
1374
1/2
✓ Branch 2 taken 168 times.
✗ Branch 3 not taken.
168 vec2 P = point_[i] + t*U;
1375
1/2
✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
168 point_.push_back(P);
1376
1/2
✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
168 v2T_.push_back(NO_INDEX);
1377 168 index_t v = nv_;
1378 168 ++nv_;
1379 168 return v;
1380 }
1381
1382 32 void CDT2d::insert(
1383 index_t nb_points, const double* points,
1384 index_t* indices, bool remove_unreferenced_vertices
1385 ) {
1386 CDT_LOG("Inserting " << nb_points << " points");
1387
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
32 debug_check_consistency();
1388
1389 // Compute spatial sort
1390 32 vector<index_t> sorted_indices;
1391
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
32 compute_BRIO_order(nb_points, points, sorted_indices, 2, 2);
1392
1393 // Insert vertices one by one, following the order given
1394 // by spatial sort.
1395
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
32 if(remove_unreferenced_vertices) {
1396
1397 // Pre-allocate memory
1398 point_.reserve(point_.size()+nb_points);
1399 v2T_.reserve(v2T_.size()+nb_points);
1400
1401 // Insert the points and vertices one by one, following
1402 // spatial sort order.
1403 index_t hint = NO_INDEX;
1404 for(index_t i=0; i<nb_points; ++i) {
1405 point_.push_back(vec2(points+2*sorted_indices[i]));
1406 index_t v = CDTBase2d::insert(point_.size()-1, hint);
1407
1408 // If it was a duplicated point, then remove the point
1409 if(point_.size() > nv()) {
1410 point_.pop_back();
1411 }
1412
1413 indices[sorted_indices[i]] = v;
1414 hint = vT(v);
1415 }
1416
1417 } else {
1418
1419 // Insert all the points in the point_ vector
1420 32 index_t v_offset = nv();
1421
1/2
✓ Branch 2 taken 32 times.
✗ Branch 3 not taken.
32 point_.reserve(point_.size()+nb_points);
1422
2/2
✓ Branch 0 taken 222718 times.
✓ Branch 1 taken 32 times.
222750 for(index_t i=0; i<nb_points; ++i) {
1423
1/2
✓ Branch 2 taken 222718 times.
✗ Branch 3 not taken.
222718 point_.push_back(vec2(points+2*i));
1424 }
1425
1426 // Resize vertex-to-triangle array accordingly,
1427 // and update number of points
1428
1/2
✓ Branch 2 taken 32 times.
✗ Branch 3 not taken.
32 v2T_.resize(v2T_.size()+nb_points, NO_INDEX);
1429 32 nv_+=nb_points;
1430
1431 // Now insert the vertices in the triangulation,
1432 // following the order of spatial search (but
1433 // this will not change the order of the points,
1434 // in contrast with the "remove_unreferenced_vertices"
1435 // alternative). In the end, each duplicated point
1436 // v has vT(v) == NO_INDEX (no incident triangle).
1437 32 index_t hint = NO_INDEX;
1438
2/2
✓ Branch 0 taken 222718 times.
✓ Branch 1 taken 32 times.
222750 for(index_t i=0; i<nb_points; ++i) {
1439 668154 index_t v = CDTBase2d::insert(
1440
2/4
✓ Branch 1 taken 222718 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 222718 times.
✗ Branch 5 not taken.
222718 v_offset+sorted_indices[i], hint
1441 );
1442
1/2
✓ Branch 0 taken 222718 times.
✗ Branch 1 not taken.
222718 if(indices != nullptr) {
1443
1/2
✓ Branch 1 taken 222718 times.
✗ Branch 2 not taken.
222718 indices[sorted_indices[i]] = v;
1444 }
1445
1/2
✓ Branch 1 taken 222718 times.
✗ Branch 2 not taken.
222718 hint = vT(v);
1446 }
1447 }
1448 CDT_LOG("Inserted.");
1449
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
32 debug_check_consistency();
1450 32 }
1451
1452 void CDT2d::save(const std::string& filename) const {
1453 #ifndef GEOGRAM_PSM
1454 Mesh M;
1455 M.vertices.set_dimension(2);
1456 for(const vec2& P: point_) {
1457 M.vertices.create_vertex(P.data());
1458 }
1459 for(index_t t=0; t<nT(); ++t) {
1460 index_t i = Tv(t,0);
1461 index_t j = Tv(t,1);
1462 index_t k = Tv(t,2);
1463 M.facets.create_triangle(i,j,k);
1464 }
1465
1466
1467 Attribute<double> tex_coord;
1468 tex_coord.create_vector_attribute(
1469 M.facet_corners.attributes(), "tex_coord", 2
1470 );
1471 static double triangle_tex[3][2] = {
1472 {0.0, 0.0},
1473 {1.0, 0.0},
1474 {0.0, 1.0}
1475 };
1476 for(index_t c: M.facet_corners) {
1477 tex_coord[2*c] = triangle_tex[c%3][0];
1478 tex_coord[2*c+1] = triangle_tex[c%3][1];
1479 }
1480
1481 Attribute<bool> constraint(M.facet_corners.attributes(), "constraint");
1482 for(index_t c: M.facet_corners) {
1483 index_t t = c/3;
1484 index_t lv = c%3;
1485 constraint[c] =
1486 Tedge_is_constrained(t, (lv+1)%3) ||
1487 Tedge_is_constrained(t, (lv+2)%3) ;
1488 }
1489
1490 for(index_t t=0; t<nT(); ++t) {
1491 for(index_t le=0; le<3; ++le) {
1492 if(Tedge_is_constrained(t,le)) {
1493 index_t v1 = Tv(t, (le+1)%3);
1494 index_t v2 = Tv(t, (le+2)%3);
1495 M.edges.create_edge(v1,v2);
1496 }
1497 }
1498 }
1499
1500
1501 M.facets.connect();
1502
1503 mesh_save(M, filename);
1504 #else
1505 if(!String::string_ends_with(filename,".obj")) {
1506 Logger::err("CDT_2d")
1507 << "save() only supports .obj file format in PSM"
1508 << std::endl;
1509 return;
1510 }
1511 std::ofstream out(filename);
1512 for(index_t v=0; v<nv(); ++v) {
1513 out << "v " << point(v) << " " << 0.0 << std::endl;
1514 }
1515 for(index_t t=0; t<nT(); ++t) {
1516 out << "f " << Tv(t,0)+1 << " " << Tv(t,1)+1 << " " << Tv(t,2)+1
1517 << std::endl;
1518 }
1519
1520 for(index_t t=0; t<nT(); ++t) {
1521 for(index_t le=0; le<3; ++le) {
1522 if(Tedge_is_constrained(t,le)) {
1523 index_t v1 = Tv(t,(le+1)%3);
1524 index_t v2 = Tv(t,(le+2)%3);
1525 out << "l " << v1+1 << " " << v2+1 << std::endl;
1526 }
1527 }
1528 }
1529 #endif
1530 }
1531
1532 /**************************************************************************/
1533
1534 391 ExactCDT2d::ExactCDT2d():
1535 391 use_pred_cache_insert_buffer_(false) {
1536 #ifdef GEOGRAM_USE_EXACT_NT
1537 CDTBase2d::exact_incircle_ = true;
1538 #else
1539 // Since incircle() with expansions computes approximated
1540 // lifted coordinate, we need to activate additional
1541 // checks for Delaunayization.
1542 391 CDTBase2d::exact_incircle_ = false;
1543 #endif
1544 391 }
1545
1546 782 ExactCDT2d::~ExactCDT2d() {
1547 782 }
1548
1549 8067 void ExactCDT2d::clear() {
1550 8067 point_.resize(0);
1551 8067 id_.resize(0);
1552 8067 pred_cache_.clear();
1553 8067 pred_cache_insert_buffer_.resize(0);
1554 8067 use_pred_cache_insert_buffer_ = false;
1555 8067 cnstr_operand_bits_.resize(0);
1556 8067 constraints_.resize(0);
1557 8067 CDTBase2d::clear();
1558 #ifndef GEOGRAM_USE_EXACT_NT
1559 8067 length_.resize(0);
1560 #endif
1561 8067 }
1562
1563 void ExactCDT2d::create_enclosing_triangle(
1564 const ExactPoint& p1, const ExactPoint& p2, const ExactPoint& p3
1565 ) {
1566 geo_assert(nv() == 0);
1567 geo_assert(nT() == 0);
1568 #ifndef GEOGRAM_USE_EXACT_NT
1569 geo_debug_assert(length_.size() == 0);
1570 #endif
1571 add_point(p1);
1572 add_point(p2);
1573 add_point(p3);
1574 CDTBase2d::create_enclosing_triangle(0,1,2);
1575 }
1576
1577 8173 void ExactCDT2d::create_enclosing_quad(
1578 const ExactPoint& p1, const ExactPoint& p2,
1579 const ExactPoint& p3, const ExactPoint& p4
1580 ) {
1581
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 8173 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
8173 geo_assert(nv() == 0);
1582
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 8173 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
8173 geo_assert(nT() == 0);
1583 #ifndef GEOGRAM_USE_EXACT_NT
1584
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 8173 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
8173 geo_debug_assert(length_.size() == 0);
1585 #endif
1586 8173 add_point(p1);
1587 8173 add_point(p2);
1588 8173 add_point(p3);
1589 8173 add_point(p4);
1590 8173 CDTBase2d::create_enclosing_quad(0,1,2,3);
1591 8173 }
1592
1593 46601 index_t ExactCDT2d::insert(const ExactPoint& p, index_t id, index_t hint) {
1594 #ifndef GEOGRAM_USE_EXACT_NT
1595
1/6
✗ Branch 2 not taken.
✓ Branch 3 taken 46601 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
46601 geo_debug_assert(nv() == length_.size());
1596 #endif
1597 46601 debug_check_consistency();
1598 46601 add_point(p,id);
1599 46601 index_t v = CDTBase2d::insert(point_.size()-1, hint);
1600 // If inserted point already existed in
1601 // triangulation, then nv() did not increase
1602
2/2
✓ Branch 2 taken 47 times.
✓ Branch 3 taken 46554 times.
46601 if(point_.size() > nv()) {
1603 47 point_.pop_back();
1604 47 id_.pop_back();
1605 #ifndef GEOGRAM_USE_EXACT_NT
1606 47 length_.pop_back();
1607 #endif
1608 }
1609 46601 debug_check_consistency();
1610 #ifndef GEOGRAM_USE_EXACT_NT
1611
1/6
✗ Branch 2 not taken.
✓ Branch 3 taken 46601 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
46601 geo_debug_assert(nv() == length_.size());
1612 #endif
1613 46601 return v;
1614 }
1615
1616 79293 void ExactCDT2d::add_point(const ExactPoint& p, index_t id) {
1617 79293 point_.push_back(p);
1618 79293 id_.push_back(id);
1619 #ifndef GEOGRAM_USE_EXACT_NT
1620 158586 length_.push_back(
1621
6/12
✓ Branch 1 taken 79293 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 79293 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 79293 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 79293 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 79293 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 79293 times.
✗ Branch 17 not taken.
158586 (geo_sqr(p.x) + geo_sqr(p.y)).estimate() /
1622
4/8
✓ Branch 1 taken 79293 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 79293 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 79293 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 79293 times.
✗ Branch 11 not taken.
79293 geo_sqr(p.w).estimate()
1623 );
1624 #endif
1625 79293 }
1626
1627 46682 void ExactCDT2d::begin_insert_transaction() {
1628 46682 use_pred_cache_insert_buffer_ = true;
1629 46682 }
1630
1631 46554 void ExactCDT2d::commit_insert_transaction() {
1632
2/2
✓ Branch 2 taken 321789 times.
✓ Branch 3 taken 46554 times.
414897 for(const auto& it: pred_cache_insert_buffer_) {
1633
1/2
✓ Branch 1 taken 321789 times.
✗ Branch 2 not taken.
321789 pred_cache_[it.first] = it.second;
1634 }
1635 46554 pred_cache_insert_buffer_.resize(0);
1636 46554 use_pred_cache_insert_buffer_ = false;
1637 46554 }
1638
1639 128 void ExactCDT2d::rollback_insert_transaction() {
1640 128 pred_cache_insert_buffer_.resize(0);
1641 128 use_pred_cache_insert_buffer_ = false;
1642 128 }
1643
1644 /**
1645 * \brief Tests the parity of the permutation of a list of
1646 * three distinct indices with respect to the canonical order.
1647 */
1648 1289032 static bool odd_order(index_t i, index_t j, index_t k) {
1649 // Implementation: sort the elements (bubble sort is OK for
1650 // such a small number), and invert parity each time
1651 // two elements are swapped.
1652 1289032 index_t tab[3] = { i, j, k};
1653 1289032 const int N = 3;
1654 1289032 bool result = false;
1655
2/2
✓ Branch 0 taken 2578064 times.
✓ Branch 1 taken 1289032 times.
3867096 for (int I = 0; I < N - 1; ++I) {
1656
2/2
✓ Branch 0 taken 3867096 times.
✓ Branch 1 taken 2578064 times.
6445160 for (int J = 0; J < N - I - 1; ++J) {
1657
2/2
✓ Branch 0 taken 2447239 times.
✓ Branch 1 taken 1419857 times.
3867096 if (tab[J] > tab[J + 1]) {
1658 2447239 std::swap(tab[J], tab[J + 1]);
1659 2447239 result = !result;
1660 }
1661 }
1662 }
1663 1289032 return result;
1664 }
1665
1666 1289032 Sign ExactCDT2d::orient2d(index_t i, index_t j, index_t k) const {
1667
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 1289032 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
1289032 geo_debug_assert(i < nv());
1668
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 1289032 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
1289032 geo_debug_assert(j < nv());
1669
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 1289032 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
1289032 geo_debug_assert(k < nv());
1670
1671
1/2
✓ Branch 1 taken 1289032 times.
✗ Branch 2 not taken.
1289032 trindex K(i, j, k);
1672
1673
2/2
✓ Branch 0 taken 323322 times.
✓ Branch 1 taken 965710 times.
1289032 if(use_pred_cache_insert_buffer_) {
1674 646644 Sign result = PCK::orient_2d(
1675
1/2
✓ Branch 1 taken 323322 times.
✗ Branch 2 not taken.
323322 point_[K.indices[0]],
1676
1/2
✓ Branch 1 taken 323322 times.
✗ Branch 2 not taken.
323322 point_[K.indices[1]],
1677
2/4
✓ Branch 1 taken 323322 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 323322 times.
✗ Branch 5 not taken.
323322 point_[K.indices[2]]
1678 323322 );
1679
2/4
✓ Branch 1 taken 323322 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 323322 times.
✗ Branch 5 not taken.
323322 pred_cache_insert_buffer_.push_back(std::make_pair(K, result));
1680
2/2
✓ Branch 1 taken 138297 times.
✓ Branch 2 taken 185025 times.
323322 if(odd_order(i,j,k)) {
1681 138297 result = Sign(-result);
1682 }
1683 323322 return result;
1684 }
1685
1686 bool inserted;
1687 965710 std::map<trindex, Sign>::iterator it;
1688
2/4
✓ Branch 1 taken 965710 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 965710 times.
✗ Branch 5 not taken.
965710 std::tie(it,inserted) = pred_cache_.insert(std::make_pair(K,ZERO));
1689 Sign result;
1690
1691
2/2
✓ Branch 0 taken 347616 times.
✓ Branch 1 taken 618094 times.
965710 if(inserted) {
1692 695232 result = PCK::orient_2d(
1693
1/2
✓ Branch 1 taken 347616 times.
✗ Branch 2 not taken.
347616 point_[K.indices[0]],
1694
1/2
✓ Branch 1 taken 347616 times.
✗ Branch 2 not taken.
347616 point_[K.indices[1]],
1695
2/4
✓ Branch 1 taken 347616 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 347616 times.
✗ Branch 5 not taken.
347616 point_[K.indices[2]]
1696 );
1697 347616 it->second = result;
1698 } else {
1699 618094 result = it->second;
1700 }
1701
1702
2/2
✓ Branch 1 taken 414350 times.
✓ Branch 2 taken 551360 times.
965710 if(odd_order(i,j,k)) {
1703 414350 result = Sign(-result);
1704 }
1705
1706 965710 return result;
1707 }
1708
1709 156583 Sign ExactCDT2d::incircle(index_t i,index_t j,index_t k,index_t l) const {
1710 #ifdef GEOGRAM_USE_EXACT_NT
1711 return PCK::incircle_2d_SOS(point_[i], point_[j], point_[k], point_[l]);
1712 #else
1713 939498 return PCK::incircle_2d_SOS_with_lengths(
1714 156583 point_[i], point_[j], point_[k], point_[l],
1715 156583 length_[i], length_[j], length_[k], length_[l]
1716 156583 );
1717 #endif
1718 }
1719
1720 81 index_t ExactCDT2d::create_intersection(
1721 index_t E1, index_t i, index_t j,
1722 index_t E2, index_t k, index_t l
1723 ) {
1724
1725 81 geo_argused(i);
1726 81 geo_argused(j);
1727 81 geo_argused(k);
1728 81 geo_argused(l);
1729
1730 // Here we could use i,j,k,l directly, but it is *much better* to take
1731 // the original extremities of the constrained segments, since they have
1732 // simpler coordinates ! (i,j,k,l might be themselves vertices created
1733 // from constraints intersections, whereas constraint extremities can
1734 // only be initial vertices).
1735
1/2
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
81 i = constraints_[E1].indices[0];
1736
1/2
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
81 j = constraints_[E1].indices[1];
1737
1/2
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
81 k = constraints_[E2].indices[0];
1738
1/2
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
81 l = constraints_[E2].indices[1];
1739
1740
3/6
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 81 times.
✗ Branch 8 not taken.
81 exact::vec2h U = point_[j] - point_[i];
1741
3/6
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 81 times.
✗ Branch 8 not taken.
81 exact::vec2h V = point_[l] - point_[k];
1742
3/6
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 81 times.
✗ Branch 8 not taken.
81 exact::vec2h D = point_[k] - point_[i];
1743
1744 exact::rational t(
1745
2/4
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
162 det2x2(D.x, D.y, V.x, V.y) * U.w,
1746
2/4
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
162 det2x2(U.x, U.y, V.x, V.y) * D.w
1747
1/2
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
81 );
1748
1749
4/8
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 81 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 81 times.
✗ Branch 11 not taken.
81 point_.push_back(mix(t, point_[i], point_[j]));
1750
1/2
✓ Branch 3 taken 81 times.
✗ Branch 4 not taken.
81 Numeric::optimize_number_representation(*point_.rbegin());
1751
1752 #ifndef GEOGRAM_USE_EXACT_NT
1753 {
1754 81 const ExactPoint& p = *point_.rbegin();
1755 162 length_.push_back(
1756
6/12
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 81 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 81 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 81 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 81 times.
✗ Branch 17 not taken.
162 (geo_sqr(p.x) + geo_sqr(p.y)).estimate() /
1757
4/8
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 81 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 81 times.
✗ Branch 11 not taken.
81 geo_sqr(p.w).estimate()
1758 );
1759 }
1760 #endif
1761
1762
1763
1/2
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
81 id_.push_back(NO_INDEX);
1764 81 index_t x = point_.size()-1;
1765
1766
1/2
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
81 CDTBase2d::v2T_.push_back(NO_INDEX);
1767
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
81 geo_debug_assert(x == CDTBase2d::nv_);
1768 81 ++CDTBase2d::nv_;
1769
1770 81 return x;
1771 81 }
1772
1773 106 void ExactCDT2d::classify_triangles(
1774 const std::string& expr, bool mark_only
1775 ) {
1776
2/4
✓ Branch 1 taken 106 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 106 times.
106 if(expr == "union_cnstr_operand_bits_is_operand_id") {
1777 classify_triangles_union_cnstr_operand_bits_is_operand_id(
1778 mark_only
1779 );
1780 return;
1781 }
1782
1783
2/4
✓ Branch 1 taken 106 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 106 times.
✗ Branch 5 not taken.
106 facet_inclusion_bits_.assign(nT(), 0);
1784
1785
1/2
✓ Branch 1 taken 106 times.
✗ Branch 2 not taken.
106 DList S(*this, DLIST_S_ID);
1786
1787 // Step 1: get triangles adjacent to the border,
1788 // mark them as visited, classify them
1789
3/4
✓ Branch 1 taken 6640 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6534 times.
✓ Branch 4 taken 106 times.
6640 for(index_t t=0; t<nT(); ++t) {
1790
2/2
✓ Branch 0 taken 18754 times.
✓ Branch 1 taken 6110 times.
24864 for(index_t le=0; le<3; ++le) {
1791
3/4
✓ Branch 1 taken 18754 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 424 times.
✓ Branch 4 taken 18330 times.
18754 if(Tadj(t,le) == NO_INDEX) {
1792
1/2
✓ Branch 1 taken 424 times.
✗ Branch 2 not taken.
424 Tset_flag(t, T_VISITED_FLAG);
1793
1/2
✓ Branch 1 taken 424 times.
✗ Branch 2 not taken.
424 S.push_back(t);
1794 424 break;
1795 }
1796 }
1797 }
1798
1799 // Step 2: recursive traversal
1800
3/4
✓ Branch 1 taken 6640 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6534 times.
✓ Branch 4 taken 106 times.
6640 while(!S.empty()) {
1801
1/2
✓ Branch 1 taken 6534 times.
✗ Branch 2 not taken.
6534 index_t t1 = S.pop_back();
1802
1/2
✓ Branch 1 taken 6534 times.
✗ Branch 2 not taken.
6534 index_t t1_bits = facet_inclusion_bits_[t1];
1803
2/2
✓ Branch 0 taken 19602 times.
✓ Branch 1 taken 6534 times.
26136 for(index_t le=0; le<3; ++le) {
1804
1/2
✓ Branch 1 taken 19602 times.
✗ Branch 2 not taken.
19602 index_t t2 = Tadj(t1,le);
1805 19602 if(
1806
4/4
✓ Branch 0 taken 19178 times.
✓ Branch 1 taken 424 times.
✓ Branch 2 taken 6110 times.
✓ Branch 3 taken 13492 times.
38780 t2 != NO_INDEX &&
1807
3/4
✓ Branch 1 taken 19178 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6110 times.
✓ Branch 4 taken 13068 times.
19178 !Tflag_is_set(t2,T_VISITED_FLAG)
1808 ) {
1809 // t2 is included in the same operands as t1,
1810 // except for the operands that touch the boundary
1811 // between t1 and t2, for which inclusion changes
1812 6110 index_t t2_bits = t1_bits;
1813 6110 for(
1814
1/2
✓ Branch 1 taken 6110 times.
✗ Branch 2 not taken.
6110 index_t ecit = Tedge_cnstr_first(t1,le);
1815
2/2
✓ Branch 0 taken 2391 times.
✓ Branch 1 taken 6110 times.
8501 ecit != NO_INDEX;
1816 2391 ecit = edge_cnstr_next(ecit)
1817 ) {
1818
1/2
✓ Branch 1 taken 2391 times.
✗ Branch 2 not taken.
2391 index_t cnstr = edge_cnstr(ecit);
1819
2/4
✓ Branch 1 taken 2391 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2391 times.
✗ Branch 5 not taken.
2391 t2_bits ^= cnstr_operand_bits_[cnstr];
1820 }
1821
1/2
✓ Branch 1 taken 6110 times.
✗ Branch 2 not taken.
6110 facet_inclusion_bits_[t2] = t2_bits;
1822
1/2
✓ Branch 1 taken 6110 times.
✗ Branch 2 not taken.
6110 Tset_flag(t2, T_VISITED_FLAG);
1823
1/2
✓ Branch 1 taken 6110 times.
✗ Branch 2 not taken.
6110 S.push_back(t2);
1824 }
1825 }
1826 }
1827
1828 // Step 3: reset visited flag
1829
3/4
✓ Branch 1 taken 6640 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6534 times.
✓ Branch 4 taken 106 times.
6640 for(index_t t=0; t<nT(); ++t) {
1830
1/2
✓ Branch 1 taken 6534 times.
✗ Branch 2 not taken.
6534 Treset_flag(t, T_VISITED_FLAG);
1831 }
1832
1833 // Step 4: mark triangles to be deleted
1834
3/4
✓ Branch 1 taken 106 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 9 times.
✓ Branch 4 taken 97 times.
106 if(expr == "intersection") {
1835 9 index_t all_bits_set = 0;
1836
2/2
✓ Branch 2 taken 222 times.
✓ Branch 3 taken 9 times.
240 for(index_t e_operand_bits: cnstr_operand_bits_) {
1837 222 all_bits_set |= e_operand_bits;
1838 }
1839
3/4
✓ Branch 1 taken 507 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 498 times.
✓ Branch 4 taken 9 times.
507 for(index_t t=0; t<nT(); ++t) {
1840
3/4
✓ Branch 1 taken 498 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 431 times.
✓ Branch 4 taken 67 times.
498 if(facet_inclusion_bits_[t] != all_bits_set) {
1841
1/2
✓ Branch 1 taken 431 times.
✗ Branch 2 not taken.
431 Tset_flag(t, T_MARKED_FLAG);
1842 }
1843 }
1844 } else {
1845
8/14
✓ Branch 1 taken 97 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 88 times.
✓ Branch 4 taken 9 times.
✓ Branch 6 taken 88 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 9 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 97 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 88 times.
✓ Branch 16 taken 9 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
185 BooleanExpression E(expr == "union" ? "*" : expr);
1846
3/4
✓ Branch 1 taken 6133 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6036 times.
✓ Branch 4 taken 97 times.
6133 for(index_t t=0; t<nT(); ++t) {
1847
4/6
✓ Branch 1 taken 6036 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6036 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 3364 times.
✓ Branch 7 taken 2672 times.
6036 if(!E(facet_inclusion_bits_[t])) {
1848
1/2
✓ Branch 1 taken 3364 times.
✗ Branch 2 not taken.
3364 Tset_flag(t, T_MARKED_FLAG);
1849 }
1850 }
1851 97 }
1852
1/2
✓ Branch 0 taken 106 times.
✗ Branch 1 not taken.
106 if(!mark_only) {
1853
1/2
✓ Branch 1 taken 106 times.
✗ Branch 2 not taken.
106 remove_marked_triangles();
1854 }
1855 106 }
1856
1857 typedef std::set<index_t> SparseBits;
1858
1859 inline void sparse_bits_flip_bit(SparseBits& bits, index_t bit) {
1860 auto it = bits.find(bit);
1861 if(it == bits.end()) {
1862 bits.insert(bit);
1863 } else {
1864 bits.erase(it);
1865 }
1866 }
1867
1868 inline bool sparse_bits_is_zero(const SparseBits& bits) {
1869 return (bits.size() == 0);
1870 }
1871
1872 void ExactCDT2d::classify_triangles_union_cnstr_operand_bits_is_operand_id(
1873 bool mark_only
1874 ) {
1875
1876 DList S(*this, DLIST_S_ID);
1877 std::stack<SparseBits> Sbits;
1878
1879 // Step 1: get triangles adjacent to the border,
1880 // mark them as visited, classify them as to-delete
1881 for(index_t t=0; t<nT(); ++t) {
1882 for(index_t le=0; le<3; ++le) {
1883 if(Tadj(t,le) == NO_INDEX) {
1884 Tset_flag(t, T_VISITED_FLAG);
1885 Tset_flag(t, T_MARKED_FLAG);
1886 S.push_back(t);
1887 Sbits.push(SparseBits());
1888 break;
1889 }
1890 }
1891 }
1892
1893 // Step 2: recursive traversal
1894 while(!S.empty()) {
1895 index_t t1 = S.pop_back();
1896 std::set<index_t> t1_bits = Sbits.top();
1897 Sbits.pop();
1898 for(index_t le=0; le<3; ++le) {
1899 index_t t2 = Tadj(t1,le);
1900 if(
1901 t2 != NO_INDEX &&
1902 !Tflag_is_set(t2,T_VISITED_FLAG)
1903 ) {
1904 // t2 is included in the same operands as t1,
1905 // except for the operands that touch the boundary
1906 // between t1 and t2, for which inclusion changes
1907 SparseBits t2_bits = t1_bits;
1908 for(
1909 index_t ecit = Tedge_cnstr_first(t1,le);
1910 ecit != NO_INDEX;
1911 ecit = edge_cnstr_next(ecit)
1912 ) {
1913 index_t cnstr = edge_cnstr(ecit);
1914 sparse_bits_flip_bit(
1915 t2_bits, cnstr_operand_bits_[cnstr]
1916 );
1917 }
1918 if(sparse_bits_is_zero(t2_bits)) {
1919 Tset_flag(t2, T_MARKED_FLAG);
1920 }
1921 Tset_flag(t2, T_VISITED_FLAG);
1922 S.push_back(t2);
1923 Sbits.push(t2_bits);
1924 }
1925 }
1926 }
1927
1928 // Step 3: reset visited flag
1929 for(index_t t=0; t<nT(); ++t) {
1930 Treset_flag(t, T_VISITED_FLAG);
1931 }
1932
1933 if(!mark_only) {
1934 remove_marked_triangles();
1935 }
1936 }
1937
1938 void ExactCDT2d::save(const std::string& filename) const {
1939 #ifndef GEOGRAM_PSM
1940 Mesh M;
1941 Attribute<index_t> nb_cnstr(M.edges.attributes(),"nb_cnstr");
1942 M.vertices.set_dimension(2);
1943 for(const ExactPoint& P: point_) {
1944 double w = P.w.estimate();
1945 vec2 p(P.x.estimate() / w, P.y.estimate() / w);
1946 M.vertices.create_vertex(p.data());
1947 }
1948 for(index_t t=0; t<nT(); ++t) {
1949 index_t i = Tv(t,0);
1950 index_t j = Tv(t,1);
1951 index_t k = Tv(t,2);
1952 M.facets.create_triangle(i,j,k);
1953
1954 for(index_t le=0; le<3; ++le) {
1955 if(Tedge_is_constrained(t,le)) {
1956 index_t e = M.edges.create_edge(
1957 Tv(t,(le+1)%3), Tv(t,(le+2)%3)
1958 );
1959 nb_cnstr[e] = Tedge_cnstr_nb(t,le);
1960 }
1961 }
1962 }
1963 M.facets.connect();
1964 M.vertices.remove_isolated();
1965 mesh_save(M, filename);
1966 #else
1967 if(!String::string_ends_with(filename,".obj")) {
1968 Logger::err("CDT_2d")
1969 << "save() only supports .obj file format in PSM"
1970 << std::endl;
1971 return;
1972 }
1973 std::ofstream out(filename);
1974 for(const ExactPoint& P: point_) {
1975 double w = P.w.estimate();
1976 vec2 p(P.x.estimate() / w, P.y.estimate() / w);
1977 out << "v " << p << " " << 0.0 << std::endl;
1978 }
1979 for(index_t t=0; t<nT(); ++t) {
1980 out << "f " << Tv(t,0)+1 << " " << Tv(t,1)+1 << " " << Tv(t,2)+1
1981 << std::endl;
1982 }
1983
1984 for(index_t t=0; t<nT(); ++t) {
1985 for(index_t le=0; le<3; ++le) {
1986 if(Tedge_is_constrained(t,le)) {
1987 index_t v1 = Tv(t,(le+1)%3);
1988 index_t v2 = Tv(t,(le+2)%3);
1989 out << "l " << v1+1 << " " << v2+1 << std::endl;
1990 }
1991 }
1992 }
1993 #endif
1994 }
1995
1996 /***************************************************************************/
1997 }
1998