GCC Code Coverage Report

Directory: ./
File: lib/geogram/delaunay/CDT_2d.h
Date: 2026-09-07 02:36:43
Exec Total Coverage
Lines: 321 347 92.5%
Functions: 57 59 96.6%
Branches: 139 548 25.4%
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 #ifndef GEOGRAM_DELAUNAY_CDT_2D
41 #define GEOGRAM_DELAUNAY_CDT_2D
42
43 #include <geogram/basic/common.h>
44 #include <geogram/basic/geometry.h>
45 #include <geogram/numerics/predicates.h>
46 #include <geogram/numerics/exact_geometry.h>
47 #include <geogram/mesh/index.h>
48 #include <functional>
49
50 /**
51 * \file geogram/delaunay/CDT_2d.h
52 * \brief Constained Delaunay triangulation in 2D
53 * \details See documentation of CDT at the end of this file.
54 */
55
56 namespace GEO {
57
58 /**
59 * \brief Forward declaration of a small data structure
60 * used internally by CDTBase2d::find_intersected_edges()
61 */
62 struct CDT2d_ConstraintWalker;
63
64 /**
65 * \brief Base class for constrained Delaunay triangulation
66 * \details Manages the combinatorics of the constrained Delaunay
67 * triangulation. The points need to be stored elsewhere, and manipulated
68 * through indices, with two predicates:
69 * - orient2d(i,j,k)
70 * - incircle(i,j,k,l)
71 * and one construction:
72 * - create_intersection(i,j,k,l)
73 * See \ref CDT2d for an example of implementation
74 */
75 class GEOGRAM_API CDTBase2d {
76 public:
77 /**
78 * \brief CDTBase2d constructor
79 */
80 CDTBase2d();
81
82 /**
83 * \brief CDTBase2d destructor
84 */
85 virtual ~CDTBase2d();
86
87 /**
88 * \brief Removes everything from this triangulation
89 */
90 virtual void clear();
91
92 /**
93 * \brief Inserts a constraint
94 * \param[in] i , j the indices of the two vertices
95 * of the constrained segment
96 */
97 void insert_constraint(index_t i, index_t j);
98
99 /**
100 * \brief Recursively removes all the triangles adjacent to
101 * the border, and keeps what's surrounded by constraints
102 * \param[in] remove_internal_holes if set, triangles inside
103 * the internal closed loops of constrained edges are removed
104 * as well.
105 * \details If \p remove_internal_holes is set,
106 * closed loops inside holes are considered as
107 * "matter" (and kept), and so on and so forth. This also works
108 * if there are overlapping constraints (what counts is the
109 * number of constraints associated with each triangle edge).
110 * Note that this does not work if there is a
111 * chain of constrained internal edges (as opposed to a loop).
112 */
113 void remove_external_triangles(bool remove_internal_holes=false);
114
115 /**
116 * \brief Specifies whether a constrained Delaunay
117 * triangulation should be constructed, or just a
118 * plain constrained triangulation
119 * \param[in] delaunay true if a Delaunay triangulation
120 * should be constructed, false otherwise.
121 */
122 260 void set_delaunay(bool delaunay) {
123 260 delaunay_ = delaunay;
124 260 }
125
126 /**
127 * \brief Gets the number of triangles
128 */
129 289212785 index_t nT() const {
130 289212785 return T_.size()/3;
131 }
132
133 /**
134 * \brief Gets the number of vertices
135 */
136 44748718 index_t nv() const {
137 44748718 return nv_;
138 }
139
140 /**
141 * \brief Gets the number of constraints
142 */
143 6269 index_t ncnstr() const {
144 6269 return ncnstr_;
145 }
146
147 /**
148 * \brief Gets a vertex of a triangle
149 * \param[in] t the triangle
150 * \param[in] lv the local index of the vertex, in 0,1,2
151 * \return the global index of the vertex
152 */
153 76449087 index_t Tv(index_t t, index_t lv) const {
154
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 76449087 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 76449087 geo_debug_assert(t<nT());
155
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 76449087 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 76449087 geo_debug_assert(lv<3);
156 76449087 return T_[3*t+lv];
157 }
158
159 /**
160 * \brief Finds the local index of a vertex in a triangle
161 * \param[in] t the triangle
162 * \param[in] v the vertex
163 * \return lv such that Tv(t,lv) = v
164 */
165 605052 index_t Tv_find(index_t t, index_t v) const {
166
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 605052 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 605052 geo_debug_assert(t<nT());
167
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 605052 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 605052 geo_debug_assert(v<nv());
168 605052 return find_3(T_.data()+3*t, v);
169 }
170
171 /**
172 * \brief Gets a triangle adjacent to a triangle
173 * \param[in] t the triangle
174 * \param[in] le the local edge index, in 0,1,2
175 * \return the triangle adjacent to \p t accross \p le,
176 * or NO_INDEX if there is no such triangle
177 */
178 114167911 index_t Tadj(index_t t, index_t le) const {
179
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 114167911 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 114167911 geo_debug_assert(t<nT());
180
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 114167911 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 114167911 geo_debug_assert(le<3);
181 114167911 return Tadj_[3*t+le];
182 }
183
184 /**
185 * \brief Finds the edge accross which a triangle is
186 * adjacent to another one
187 * \param[in] t1 , t2 the two triangles
188 * \return the local edge index le such that
189 * Tadj(t1,le) = t2
190 */
191 23174627 index_t Tadj_find(index_t t1, index_t t2) const {
192
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 23174627 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 23174627 geo_debug_assert(t1<nT());
193
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 23174627 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 23174627 geo_debug_assert(t2<nT());
194 23174627 return find_3(Tadj_.data()+3*t1, t2);
195 }
196
197 /**
198 * \brief Gets a triangle incident to a given vertex
199 * \param[in] v a vertex
200 * \return a triangle t such that there exists lv in
201 * 0,1,2 such that Tv(t,lv) = v
202 */
203 777907 index_t vT(index_t v) const {
204
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 777907 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 777907 geo_debug_assert(v < nv());
205 777907 return v2T_[v];
206 }
207
208
209 /**
210 * \brief Gets the constraint associated with an edge
211 * \details
212 * When constraining segments on a CDT2d by calling
213 * insert_constraint(), some segments in the triangulation may
214 * be included in several constraints (it the constraints are
215 * co-linear and overlapping).
216 * One iterates on the constraints associated with an edge as follows:
217 * \code
218 * for(
219 * index_t ecit = Tedge_cnstr_first(t,le);
220 * ecit != NO_INDEX;
221 * ecit = edge_cnstr_next(ecit)
222 * ) {
223 * index_t cnstr = edge_cnstr(ecit);
224 * ... // do something with cnstr
225 * }
226 * \endcode
227 * where 'cnstr' corresponds to the value of ncnstr() when
228 * insert_constraint() was called for that constraint.
229 * \param[in] t a triangle
230 * \param[in] le local edge index, in 0,1,2
231 * \return the edge constraints iterator associated with this edge or
232 * NO_INDEX if the edge is not constrained.
233 * \see edge_cnstr_next(), edge_cnstr()
234 */
235 6945466 index_t Tedge_cnstr_first(index_t t, index_t le) const {
236
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 6945466 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 6945466 geo_debug_assert(t < nT());
237
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 6945466 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 6945466 geo_debug_assert(le < 3);
238 6945466 return Tecnstr_first_[3*t+le];
239 }
240
241 /**
242 * \brief Gets the successor of an edge constraint iterator
243 * \param[in] ecit the edge constraint iterator
244 * \return the edge constraint iterator to the successor of \p ecit
245 * or NO_INDEX if \p ecit is the last of the list.
246 * \see Tedge_cnstr_first(), edge_cnstr()
247 */
248 66520 index_t edge_cnstr_next(index_t ecit) const {
249 66520 return ecnstr_next_[ecit];
250 }
251
252 /**
253 * \brief Gets an edge constraint from an edge constraint iterator
254 * \param[in] ecit the edge constraint iterator. Needs to be a valid
255 * iterator, different from NO_INDEX.
256 * \return the edge constraint associated with the iterator.
257 * \see Tedge_cnstr_first(), edge_cnstr_next()
258 */
259 44806 index_t edge_cnstr(index_t ecit) const {
260 44806 return ecnstr_val_[ecit];
261 }
262
263 /**
264 * \brief Gets the number of constraints associated with a triange edge
265 * \details There can be several constraints associated with the same
266 * edge, whenever there are overlapping constraints. For instance,
267 * this function is useful to test the parity of the number of
268 * constraints when classifying inside/outside triangles
269 * in a CSG operation.
270 * \param[in] t the triangle
271 * \param[in] le the local index of the edge (0,1,2) in the triangle
272 * \return the number of constraints associated with the edge
273 */
274 103082 index_t Tedge_cnstr_nb(index_t t, index_t le) const {
275 103082 index_t result = 0;
276 103082 for(
277 103082 index_t ecit = Tedge_cnstr_first(t,le);
278
2/2
✓ Branch 0 taken 30997 times.
✓ Branch 1 taken 103082 times.
 134079 ecit != NO_INDEX;
279 30997 ecit = edge_cnstr_next(ecit)
280 ) {
281 30997 ++result;
282 }
283 103082 return result;
284 }
285
286
287 /**
288 * \brief Saves this CDT to a geogram mesh file.
289 * \param[in] filename where to save this CDT
290 */
291 virtual void save(const std::string& filename) const = 0;
292
293 /**
294 * \brief Tests whether a triangle edge is Delaunay
295 * \details returns true also for constrained edges and edges on borders
296 */
297 bool Tedge_is_Delaunay(index_t t, index_t le) const;
298
299 protected:
300 virtual void begin_insert_transaction();
301 virtual void commit_insert_transaction();
302 virtual void rollback_insert_transaction();
303
304 /**
305 * \brief Inserts a new point
306 * \param[in] v the index of the new point, supposed to be
307 * equal to nv()
308 * \param[in] hint an optional triangle, not too far away
309 * from the point to be inserted
310 * \return the index of the created point. May be different
311 * from v if the point already existed in the triangulation
312 */
313 index_t insert(index_t v, index_t hint = NO_INDEX);
314
315 /**
316 * \brief Creates the combinatorics for a first large enclosing
317 * triangle
318 * \param[in] v1 , v2 , v3 the three vertices of the first triangle,
319 * in 0,1,2
320 * \details create_enclosing_triangle() or create_enclosing_quad()
321 * need to be called before anything else
322 */
323 void create_enclosing_triangle(index_t v1, index_t v2, index_t v3);
324
325 /**
326 * \brief Creates the combinatorics for a first large enclosing
327 * quad
328 * \param[in] v1 , v2 , v3 , v4 the four vertices of the quad,
329 * in 0,1,2,3
330 * \details create_enclosing_triangle() or create_enclosing_quad()
331 * need to be called before anything else
332 */
333 void create_enclosing_quad(
334 index_t v1, index_t v2, index_t v3, index_t v4
335 );
336
337 /**
338 * \brief Sets a triangle flag
339 * \param[in] t the triangle
340 * \param[in] flag the flag, in 0..7
341 */
342 2094140 void Tset_flag(index_t t, index_t flag) {
343
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 2094140 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2094140 geo_debug_assert(t < nT());
344
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2094140 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 2094140 geo_debug_assert(flag < 8);
345 2094140 Tflags_[t] |= Numeric::uint8(1u << flag);
346 2094140 }
347
348 /**
349 * \brief Resets a triangle flag
350 * \param[in] t the triangle
351 * \param[in] flag the flag, in 0..7
352 */
353 2014807 void Treset_flag(index_t t, index_t flag) {
354
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 2014807 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2014807 geo_debug_assert(t < nT());
355
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2014807 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 2014807 geo_debug_assert(flag < 8);
356 2014807 Tflags_[t] &= Numeric::uint8(~(1u << flag));
357 2014807 }
358
359 /**
360 * \brief Tests a triangle flag
361 * \param[in] t the triangle
362 * \param[in] flag the flag, in 0..7
363 * \retval true if the flag is set
364 * \retval false otherwise
365 */
366 2633970 bool Tflag_is_set(index_t t, index_t flag) {
367
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 2633970 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2633970 geo_debug_assert(t < nT());
368
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2633970 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 2633970 geo_debug_assert(flag < 8);
369 2633970 return ((Tflags_[t] & (1u << flag)) != 0);
370 }
371
372 /**
373 * \brief Constants for list_id
374 */
375 enum {
376 DLIST_S_ID=0,
377 DLIST_Q_ID=1,
378 DLIST_N_ID=2,
379 DLIST_NB=3
380 };
381
382 /**
383 * \brief Constants for triangle flags
384 */
385 enum {
386 T_MARKED_FLAG = DLIST_NB,
387 T_VISITED_FLAG = DLIST_NB+1
388 };
389
390 /**
391 * \brief Tests whether a triangle is in a DList
392 * \param[in] t the triangle
393 * \retval true if the triangle is in a list
394 * \retval false otherwise
395 */
396 1905191 bool Tis_in_list(index_t t) const {
397 return (
398 1905191 (Tflags_[t] &
399 Numeric::uint8((1 << DLIST_NB)-1)
400 ) != 0
401 1905191 );
402 }
403
404 /**
405 * \brief Doubly connected triangle list
406 * \details DList is used to implement:
407 * - the stack S of triangles to flip in insert()
408 * - the queue Q of intersected edges in
409 * detect_intersected_edges() and constrain_edges()
410 * - the list N of new edges in constrain_edges()
411 * Everything is stored in CDBase
412 * vectors Tnext_, Tprev_ and Tflags_. As
413 * a consequence, the same triangle can be only
414 * in a single DList at the same time.
415 */
416 struct DList {
417 /**
418 * \brief Constructs an empty DList
419 * \param[in] cdt a reference to the CDTBase2d
420 * \param[in] list_id the DList id, in 0..DLIST_NB-1
421 */
422 193326 DList(CDTBase2d& cdt, index_t list_id) :
423 193326 cdt_(cdt), list_id_(list_id),
424 193326 back_(NO_INDEX), front_(NO_INDEX) {
425
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 193326 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 193326 geo_debug_assert(list_id < DLIST_NB);
426 193326 }
427
428 /**
429 * \brief Creates an uninitialized DList
430 * \details One cannot do anything with an
431 * uninitialized Dlist, except:
432 * - initializing it with DList::initialize()
433 * - testing its status with DList::initialized()
434 * - display it with DList::show()
435 */
436 542789 DList(CDTBase2d& cdt) :
437 542789 cdt_(cdt), list_id_(NO_INDEX),
438 542789 back_(NO_INDEX), front_(NO_INDEX) {
439 542789 }
440
441 /**
442 * \brief Initializes a list
443 * \param[in] list_id the DList id, in 0..DLIST_NB-1
444 */
445 425073 void initialize(index_t list_id) {
446
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 425073 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 425073 geo_debug_assert(!initialized());
447
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 425073 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 425073 geo_debug_assert(list_id < DLIST_NB);
448 425073 list_id_ = list_id;
449 425073 }
450
451 /**
452 * \brief Tests whether a DList is initialized
453 */
454 13940561 bool initialized() const {
455 13940561 return (list_id_ != NO_INDEX);
456 }
457
458 736115 ~DList() {
459
2/2
✓ Branch 1 taken 618399 times.
✓ Branch 2 taken 117716 times.
 736115 if(initialized()) {
460 618399 clear();
461 }
462 736115 }
463
464 6225221 bool empty() const {
465
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 6225221 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 6225221 geo_debug_assert(initialized());
466
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 6225221 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 6225221 geo_debug_assert(
467 (back_==NO_INDEX)==(front_==NO_INDEX)
468 );
469 6225221 return (back_==NO_INDEX);
470 }
471
472 2084464 bool contains(index_t t) const {
473
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 2084464 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2084464 geo_debug_assert(initialized());
474 2084464 return cdt_.Tflag_is_set(t, list_id_);
475 }
476
477 259966 index_t front() const {
478
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 259966 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 259966 geo_debug_assert(initialized());
479 259966 return front_;
480 }
481
482 index_t back() const {
483 geo_debug_assert(initialized());
484 return back_;
485 }
486
487 23453 index_t next(index_t t) const {
488
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 23453 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 23453 geo_debug_assert(initialized());
489
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 23453 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 23453 geo_debug_assert(contains(t));
490 23453 return cdt_.Tnext_[t];
491 }
492
493 index_t prev(index_t t) const {
494 geo_debug_assert(initialized());
495 geo_debug_assert(contains(t));
496 return cdt_.Tprev_[t];
497 }
498
499 882121 void clear() {
500
2/2
✓ Branch 1 taken 15180 times.
✓ Branch 2 taken 882121 times.
 897301 for(index_t t=front_; t!=NO_INDEX; t = cdt_.Tnext_[t]) {
501 15180 cdt_.Treset_flag(t,list_id_);
502 }
503 882121 back_ = NO_INDEX;
504 882121 front_ = NO_INDEX;
505 882121 }
506
507 index_t size() const {
508 geo_debug_assert(initialized());
509 index_t result = 0;
510 for(index_t t=front(); t!=NO_INDEX; t = next(t)) {
511 ++result;
512 }
513 return result;
514 }
515
516 1865725 void push_back(index_t t) {
517
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 1865725 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 1865725 geo_debug_assert(initialized());
518
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 1865725 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 1865725 geo_debug_assert(!cdt_.Tis_in_list(t));
519 1865725 cdt_.Tset_flag(t,list_id_);
520
2/2
✓ Branch 1 taken 400564 times.
✓ Branch 2 taken 1465161 times.
 1865725 if(empty()) {
521 400564 back_ = t;
522 400564 front_ = t;
523 400564 cdt_.Tnext_[t] = NO_INDEX;
524 400564 cdt_.Tprev_[t] = NO_INDEX;
525 } else {
526 1465161 cdt_.Tnext_[t] = NO_INDEX;
527 1465161 cdt_.Tnext_[back_] = t;
528 1465161 cdt_.Tprev_[t] = back_;
529 1465161 back_ = t;
530 }
531 1865725 }
532
533 1890011 index_t pop_back() {
534
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 1890011 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 1890011 geo_debug_assert(initialized());
535
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 1890011 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 1890011 geo_debug_assert(!empty());
536 1890011 index_t t = back_;
537 1890011 back_ = cdt_.Tprev_[back_];
538
2/2
✓ Branch 0 taken 394159 times.
✓ Branch 1 taken 1495852 times.
 1890011 if(back_ == NO_INDEX) {
539
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 394159 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 394159 geo_debug_assert(front_ == t);
540 394159 front_ = NO_INDEX;
541 } else {
542 1495852 cdt_.Tnext_[back_] = NO_INDEX;
543 }
544
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 1890011 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 1890011 geo_debug_assert(contains(t));
545 1890011 cdt_.Treset_flag(t,list_id_);
546 1890011 return t;
547 }
548
549 39466 void push_front(index_t t) {
550
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 39466 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 39466 geo_debug_assert(initialized());
551
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 39466 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 39466 geo_debug_assert(!cdt_.Tis_in_list(t));
552 39466 cdt_.Tset_flag(t,list_id_);
553
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 39466 times.
 39466 if(empty()) {
554 ✗ back_ = t;
555 ✗ front_ = t;
556 ✗ cdt_.Tnext_[t] = NO_INDEX;
557 ✗ cdt_.Tprev_[t] = NO_INDEX;
558 } else {
559 39466 cdt_.Tprev_[t] = NO_INDEX;
560 39466 cdt_.Tprev_[front_] = t;
561 39466 cdt_.Tnext_[t] = front_;
562 39466 front_ = t;
563 }
564 39466 }
565
566 index_t pop_front() {
567 geo_debug_assert(initialized());
568 geo_debug_assert(!empty());
569 index_t t = front_;
570 front_ = cdt_.Tnext_[front_];
571 if(front_ == NO_INDEX) {
572 geo_debug_assert(back_ == t);
573 back_ = NO_INDEX;
574 } else {
575 cdt_.Tprev_[front_] = NO_INDEX;
576 }
577 geo_debug_assert(contains(t));
578 cdt_.Treset_flag(t,list_id_);
579 return t;
580 }
581
582 void remove(index_t t) {
583 geo_debug_assert(initialized());
584 if(t == front_) {
585 pop_front();
586 } else if(t == back_) {
587 pop_back();
588 } else {
589 geo_debug_assert(contains(t));
590 index_t t_prev = cdt_.Tprev_[t];
591 index_t t_next = cdt_.Tnext_[t];
592 cdt_.Tprev_[t_next] = t_prev;
593 cdt_.Tnext_[t_prev] = t_next;
594 cdt_.Treset_flag(t,list_id_);
595 }
596 }
597
598 void show(std::ostream& out = std::cerr) const {
599 switch(list_id_) {
600 case DLIST_S_ID:
601 out << "S";
602 break;
603 case DLIST_Q_ID:
604 out << "Q";
605 break;
606 case DLIST_N_ID:
607 out << "N";
608 break;
609 case NO_INDEX:
610 out << "<uninitialized list>";
611 break;
612 default:
613 out << "<unknown list id:" << list_id_ << ">";
614 break;
615 }
616 out << "=";
617 for(index_t t=front(); t!=NO_INDEX; t = next(t)) {
618 out << t << ";";
619 }
620 out << std::endl;
621 }
622
623 private:
624 CDTBase2d& cdt_;
625 index_t list_id_;
626 index_t back_;
627 index_t front_;
628 };
629
630 /**
631 * \brief Inserts a vertex in an edge
632 * \param[in] v the vertex to be inserted
633 * \param[in] t a triangle incident to the edge
634 * \param[in] le the local index of the edge in \p t
635 * \param[out] S DList of created triangles, ignored if uninitialized
636 */
637 void insert_vertex_in_edge(index_t v, index_t t, index_t le, DList& S);
638
639 /**
640 * \brief Inserts a vertex in an edge
641 * \param[in] v the vertex to be inserted
642 * \param[in] t a triangle incident to the edge
643 * \param[in] le the local index of the edge in \p t
644 */
645 5933 void insert_vertex_in_edge(index_t v, index_t t, index_t le) {
646 5933 DList S(*this);
647
1/2
✓ Branch 1 taken 5933 times.
✗ Branch 2 not taken.
 5933 insert_vertex_in_edge(v,t,le,S);
648 5933 }
649
650 /**
651 * \brief Inserts a vertex in a triangle
652 * \param[in] v the vertex to be inserted
653 * \param[in] t the triangle
654 * \param[out] S optional DList of created triangles
655 */
656 void insert_vertex_in_triangle(index_t v, index_t t, DList& S);
657
658 /**
659 * \brief Finds the edges intersected by a constraint
660 * \param[in] i , j the two vertices of the constraint
661 * \param[out] Q for each intersected edge, a triangle t
662 * will be pushed-back to Q, such that vT(t,1) and
663 * vT(t,2) are the extremities of the intersected edge.
664 * In addition, each triangle t is marked.
665 * \details If a vertex k that is exactly on the constraint
666 * is found, then traversal stops there and k is returned.
667 * One can find the remaining intersections by continuing
668 * to call the function with (k,j) until \p j is returned.
669 * \return the first vertex on [i,j] encountered when
670 * traversing the segment [i,j].
671 */
672 index_t find_intersected_edges(index_t i, index_t j, DList& Q);
673
674 /**
675 * \brief Used by find_intersected_edges()
676 */
677 void walk_constraint_v(CDT2d_ConstraintWalker& W);
678
679 /**
680 * \brief Used by find_intersected_edges()
681 */
682 void walk_constraint_t(CDT2d_ConstraintWalker& W, DList& Q);
683
684 /**
685 * \brief Constrains an edge by iteratively flipping
686 * the intersected edges.
687 * \param[in] i , j the extremities of the edge
688 * \param[in] Q the list of intersected edges, computed by
689 * find_intersected_edges()
690 * \param[out] N optional DList with the new edges
691 * that need to be re-Delaunized by find_intersected_edges(),
692 * ignored if uninitialized
693 */
694 void constrain_edges(index_t i, index_t j, DList& Q, DList& N);
695
696 /**
697 * \brief Restores Delaunay condition starting from the
698 * triangles incident to a given vertex.
699 * \details This version uses internally a stack, initialized
700 * with the triangles incident to the vertex.
701 * \param[in] from_v the vertex. Cannot be a vertex incident to the
702 * border.
703 */
704 void Delaunayize_vertex_neighbors(index_t from_v);
705
706 /**
707 * \brief Restores Delaunay condition starting from the
708 * triangles incident to a given vertex.
709 * \param[in] v the vertex
710 * \param[in] S a stack of triangles, initialized with
711 * the triangles incident to the vertex. Each triangle t
712 * is Trot()-ed in such a way that the vertex v
713 * corresponds to Vt(t,0)
714 * \details Each time a triangle edge is swapped, the
715 * two new neighbors are recursively examined.
716 */
717 void Delaunayize_vertex_neighbors(index_t v, DList& S);
718
719 /**
720 * \brief Restores Delaunay condition for a set of
721 * edges after inserting a constrained edge
722 * \param[in] N the edges for which Delaunay condition
723 * should be restored.
724 */
725 void Delaunayize_new_edges(DList& N);
726
727
728 /**
729 * \brief Sets all the combinatorial information
730 * of a triangle and edge flags
731 * \param[in] t the triangle
732 * \param[in] v1 , v2 , v3 the three vertices
733 * \param[in] adj1 , adj2 , adj3 the three triangles
734 * adjacent to \p t
735 * \param[in] e1cnstr , e2cnstr , e3cnstr optional
736 * edge constraints
737 */
738 2358567 void Tset(
739 index_t t,
740 index_t v1, index_t v2, index_t v3,
741 index_t adj1, index_t adj2, index_t adj3,
742 index_t e1cnstr = NO_INDEX,
743 index_t e2cnstr = NO_INDEX,
744 index_t e3cnstr = NO_INDEX
745 ) {
746
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 2358567 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2358567 geo_debug_assert(t < nT());
747
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 2358567 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2358567 geo_debug_assert(v1 < nv());
748
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 2358567 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2358567 geo_debug_assert(v2 < nv());
749
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 2358567 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2358567 geo_debug_assert(v3 < nv());
750
4/10
✓ Branch 1 taken 265118 times.
✓ Branch 2 taken 2093449 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 265118 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2358567 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 2358567 geo_debug_assert(adj1 < nT() || adj1 == NO_INDEX);
751
4/10
✓ Branch 1 taken 173404 times.
✓ Branch 2 taken 2185163 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 173404 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2358567 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 2358567 geo_debug_assert(adj2 < nT() || adj2 == NO_INDEX);
752
4/10
✓ Branch 1 taken 167592 times.
✓ Branch 2 taken 2190975 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 167592 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2358567 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 2358567 geo_debug_assert(adj3 < nT() || adj3 == NO_INDEX);
753
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2358567 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 2358567 geo_debug_assert(v1 != v2);
754
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2358567 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 2358567 geo_debug_assert(v2 != v3);
755
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2358567 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 2358567 geo_debug_assert(v3 != v1);
756
3/8
✓ Branch 0 taken 60672 times.
✓ Branch 1 taken 2297895 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 60672 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 2358567 geo_debug_assert(adj1 != adj2 || adj1 == NO_INDEX);
757
3/8
✓ Branch 0 taken 33622 times.
✓ Branch 1 taken 2324945 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 33622 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 2358567 geo_debug_assert(adj2 != adj3 || adj2 == NO_INDEX);
758
3/8
✓ Branch 0 taken 52436 times.
✓ Branch 1 taken 2306131 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 52436 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 2358567 geo_debug_assert(adj3 != adj1 || adj3 == NO_INDEX);
759
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 2358567 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2358567 geo_debug_assert(orient2d(v1,v2,v3) != ZERO);
760 2358567 T_[3*t ] = v1;
761 2358567 T_[3*t+1] = v2;
762 2358567 T_[3*t+2] = v3;
763 2358567 Tadj_[3*t ] = adj1;
764 2358567 Tadj_[3*t+1] = adj2;
765 2358567 Tadj_[3*t+2] = adj3;
766 2358567 Tecnstr_first_[3*t] = e1cnstr;
767 2358567 Tecnstr_first_[3*t+1] = e2cnstr;
768 2358567 Tecnstr_first_[3*t+2] = e3cnstr;
769 2358567 v2T_[v1] = t;
770 2358567 v2T_[v2] = t;
771 2358567 v2T_[v3] = t;
772 2358567 }
773
774 /**
775 * \brief Rotates indices in triangle t in such a way
776 * that a given vertex becomes vertex 0
777 * \details On exit, vertex \p lv of \p t becomes vertex 0
778 * \param[in] t a triangle index
779 * \param[in] lv local vertex index in 0,1,2
780 */
781 149612 void Trot(index_t t, index_t lv) {
782
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 149612 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 149612 geo_debug_assert(t < nT());
783
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 149612 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 149612 geo_debug_assert(lv < 3);
784
2/2
✓ Branch 0 taken 87756 times.
✓ Branch 1 taken 61856 times.
 149612 if(lv != 0) {
785 87756 index_t i = 3*t+lv;
786 87756 index_t j = 3*t+((lv+1)%3);
787 87756 index_t k = 3*t+((lv+2)%3);
788 789804 Tset(
789 t,
790 87756 T_[i], T_[j], T_[k],
791 87756 Tadj_[i], Tadj_[j], Tadj_[k],
792 87756 Tecnstr_first_[i], Tecnstr_first_[j], Tecnstr_first_[k]
793 );
794 }
795 149612 }
796
797 /**
798 * \brief Swaps an edge.
799 * \details Swaps edge 0 of \p t1.
800 * Vertex 0 of \p t1 is vertex 0 of
801 * the two new triangles.
802 * \param[in] t1 a triangle index. Its edge
803 * opposite to vertex 0 is swapped
804 * \param[in] swap_t1_t2 if set, swap which triangle will be
805 * t1 and which triangle will be Tadj(t1,0) in the
806 * new pair of triange (needed for two configurations
807 * of the optimized constraint enforcement algorithm).
808 */
809 void swap_edge(index_t t1, bool swap_t1_t2=false);
810
811 /**
812 * \brief Sets a triangle adjacency relation
813 * \param[in] t a triangle
814 * \param[in] le local edge index, in 0,1,2
815 * \param[in] adj the triangle adjacent to \p t
816 * accross \p le
817 */
818 3002056 void Tadj_set(index_t t, index_t le, index_t adj) {
819
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 3002056 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 3002056 geo_debug_assert(t < nT());
820
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 3002056 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 3002056 geo_debug_assert(adj < nT());
821
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 3002056 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 3002056 geo_debug_assert(le < 3);
822 3002056 Tadj_[3*t+le] = adj;
823 3002056 }
824
825 /**
826 * \brief Gets the neighboring triangle vertex
827 * opposite to a given vertex
828 */
829 index_t Topp(index_t t, index_t e=0) const {
830 index_t t2 = Tadj(t,e);
831 if(t2 == NO_INDEX) {
832 return NO_INDEX;
833 }
834 index_t e2 = Tadj_find(t2,t);
835 return Tv(t2,e2);
836 }
837
838 /**
839 * \brief After having changed connections from triangle
840 * to a neighbor, creates connections from neighbor
841 * to triangle.
842 * \details edge flags are copied from the neighbor to \p t1.
843 * If there is no triangle accross \p le1, then
844 * nothing is done
845 * \param[in] t1 a triangle
846 * \param[in] le1 a local edge of \p t1, in 0,1,2
847 * \param[in] prev_t2_adj_e2 the triangle adjacent to t2 that
848 * \p t1 will replace, where t2 = Tadj(t1,le1)
849 */
850 3309855 void Tadj_back_connect(
851 index_t t1, index_t le1, index_t prev_t2_adj_e2
852 ) {
853
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 3309855 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 3309855 geo_debug_assert(t1 < nT());
854
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 3309855 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 3309855 geo_debug_assert(le1 < 3);
855 3309855 index_t t2 = Tadj(t1,le1);
856
2/2
✓ Branch 0 taken 307799 times.
✓ Branch 1 taken 3002056 times.
 3309855 if(t2 == NO_INDEX) {
857 307799 return;
858 }
859 3002056 index_t le2 = Tadj_find(t2,prev_t2_adj_e2);
860 3002056 Tadj_set(t2,le2,t1);
861 3002056 Tset_edge_cnstr_first(t1,le1,Tedge_cnstr_first(t2,le2));
862 }
863
864 /**
865 * \brief Creates a new triangle
866 * \return the index of the new triange
867 */
868 686413 index_t Tnew() {
869 686413 index_t t = nT();
870 686413 index_t nc = (t+1)*3; // new number of corners
871 686413 T_.resize(nc, NO_INDEX);
872 686413 Tadj_.resize(nc, NO_INDEX);
873 686413 Tecnstr_first_.resize(nc, NO_INDEX);
874
1/2
✓ Branch 1 taken 686413 times.
✗ Branch 2 not taken.
 686413 Tflags_.resize(t+1,0);
875 686413 Tnext_.resize(t+1,NO_INDEX);
876 686413 Tprev_.resize(t+1,NO_INDEX);
877 686413 return t;
878 }
879
880 /**
881 * \brief Sets the constraints list associated with an edge
882 * \param[in] t a triangle
883 * \param[in] le local edge index, in 0,1,2
884 * \param[in] ecit the edge constraint iterator that points to
885 * the first constraint associated with the edge
886 */
887 3580689 void Tset_edge_cnstr_first(
888 index_t t, index_t le, index_t ecit
889 ) {
890
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 3580689 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 3580689 geo_debug_assert(t < nT());
891
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 3580689 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 3580689 geo_debug_assert(le < 3);
892 3580689 Tecnstr_first_[3*t+le] = ecit;
893 3580689 }
894
895 /**
896 * \brief Adds a constraint to a triangle edge
897 * \param[in] t a triangle
898 * \param[in] le local edge index, in 0,1,2
899 * \param[in] cnstr_id the constraint
900 */
901 261820 void Tadd_edge_cnstr(
902 index_t t, index_t le, index_t cnstr_id
903 ) {
904
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 261820 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 261820 geo_debug_assert(t < nT());
905
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 261820 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 261820 geo_debug_assert(le < 3);
906 // Check whether the edge is already constrained with the
907 // same constraint.
908 // TODO (if possible): understand how this can happen and
909 // remove this bloc of code that is not super elegant
910 // (it seems to be when we arrive at j and coming from a vertex
911 // traversed by the edge, both conditions make the constraint
912 // added to the traversed edge).
913 261820 for(
914 261820 index_t ecit = Tedge_cnstr_first(t,le);
915
2/2
✓ Branch 0 taken 36482 times.
✓ Branch 1 taken 258470 times.
 294952 ecit != NO_INDEX;
916 33132 ecit = edge_cnstr_next(ecit)
917 ) {
918
2/2
✓ Branch 1 taken 3350 times.
✓ Branch 2 taken 33132 times.
 36482 if(edge_cnstr(ecit) == cnstr_id) {
919 3350 return;
920 }
921 }
922 258470 ecnstr_val_.push_back(cnstr_id);
923
2/4
✓ Branch 1 taken 258470 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 258470 times.
✗ Branch 5 not taken.
 258470 ecnstr_next_.push_back(Tedge_cnstr_first(t,le));
924 258470 Tset_edge_cnstr_first(t,le, ecnstr_val_.size()-1);
925 }
926
927 /**
928 * \brief Adds a constraint to a triangle edge and to the neighboring
929 * edge if it exists
930 * \param[in] t a triangle
931 * \param[in] le local edge index, in 0,1,2
932 * \param[in] cnstr_id the constraint
933 */
934 261820 void Tadd_edge_cnstr_with_neighbor(
935 index_t t, index_t le, index_t cnstr_id
936 ) {
937
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 261820 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 261820 geo_debug_assert(t < nT());
938
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 261820 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 261820 geo_debug_assert(le < 3);
939 #ifdef GEO_DEBUG
940 261820 index_t t_e_cnstr_first = Tedge_cnstr_first(t,le);
941 #endif
942 261820 Tadd_edge_cnstr(t, le, cnstr_id);
943 261820 index_t t2 = Tadj(t,le);
944
2/2
✓ Branch 0 taken 142819 times.
✓ Branch 1 taken 119001 times.
 261820 if(t2 != NO_INDEX) {
945 142819 index_t le2 = Tadj_find(t2,t);
946 // Sanity check: make sure the two edges always share the
947 // same constraint list.
948
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 142819 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 142819 geo_debug_assert(Tedge_cnstr_first(t2,le2) == t_e_cnstr_first);
949 142819 Tset_edge_cnstr_first(t2,le2,Tedge_cnstr_first(t,le));
950 }
951 261820 }
952
953 /**
954 * \brief Tests whether an edge is constrained
955 * \param[in] t a triangle
956 * \param[in] le local edge index, in 0,1,2
957 * \retval true if the edge is constrained
958 * \retval false otherwise
959 */
960 2467172 bool Tedge_is_constrained(index_t t, index_t le) const {
961 2467172 return (Tedge_cnstr_first(t,le) != NO_INDEX);
962 }
963
964 /**
965 * \brief Calls a user-defined function for each triangle
966 * around a vertex
967 * \param[in] v the vertex
968 * \param[in] doit the function, that takes as argument the
969 * current triangle t and the local index lv of \p v in t.
970 * The function returns true if iteration is finished and can be
971 * exited, false otherwise.
972 */
973 264235 void for_each_T_around_v(
974 index_t v, std::function<bool(index_t t, index_t lv)> doit
975 ) {
976 264235 index_t t = vT(v);
977 264235 index_t lv = NO_INDEX;
978 do {
979 419087 lv = Tv_find(t,v);
980
2/2
✓ Branch 1 taken 202033 times.
✓ Branch 2 taken 217054 times.
 419087 if(doit(t,lv)) {
981 202033 return;
982 }
983 217054 t = Tadj(t, (lv+1)%3);
984
5/6
✓ Branch 1 taken 217054 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 154852 times.
✓ Branch 4 taken 62202 times.
✓ Branch 5 taken 154852 times.
✓ Branch 6 taken 62202 times.
 217054 } while(t != vT(v) && t != NO_INDEX);
985
986 // We are done, this was an interior vertex
987
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 62202 times.
 62202 if(t != NO_INDEX) {
988 ✗ return;
989 }
990
991 // It was a vertex on the border, so we need
992 // to traverse the triangle fan in the other
993 // direction until we reach the border again
994 62202 t = vT(v);
995 62202 lv = Tv_find(t,v);
996 62202 t = Tadj(t, (lv+2)%3);
997
1/2
✓ Branch 0 taken 99905 times.
✗ Branch 1 not taken.
 99905 while(t != NO_INDEX) {
998 99905 lv = Tv_find(t,v);
999
2/2
✓ Branch 1 taken 62202 times.
✓ Branch 2 taken 37703 times.
 99905 if(doit(t,lv)) {
1000 62202 return;
1001 }
1002 37703 t = Tadj(t, (lv+2)%3);
1003 }
1004 }
1005
1006
1007 /**
1008 * \brief Locates a vertex
1009 * \param[in] v the vertex index
1010 * \param[in] hint an optional triangle, not too far away from the
1011 * point to be inserted
1012 * \param[out] orient a pointer to the three orientations in the
1013 * triangle. If one of them is zero, the point is on an edge, and
1014 * if two of them are zero, it is on a vertex.
1015 * \return a triangle that contains \p v
1016 */
1017 index_t locate(
1018 index_t v, index_t hint = NO_INDEX, Sign* orient = nullptr
1019 ) const;
1020
1021 /**
1022 * \brief Tests whether triange t and its neighbor accross edge 0 form
1023 * a strictly convex quad
1024 * \retval true if triange \p t and its neighbor accross edge 0 form
1025 * a strictly convex quad
1026 * \retval false otherwise
1027 */
1028 bool is_convex_quad(index_t t) const;
1029
1030 /**
1031 * \brief Orientation predicate
1032 * \param[in] i , j , k three vertices
1033 * \return the sign of det(pj-pi,pk-pi)
1034 */
1035 virtual Sign orient2d(index_t i,index_t j,index_t k) const=0;
1036
1037 /**
1038 * \brief Incircle predicate
1039 * \param[in] i , j , k the three vertices of a triangle
1040 * \param[in] l another vertex
1041 * \retval POSITIVE if \p l is inside the circumscribed circle of
1042 * the triangle
1043 * \retval ZERO if \p l is on the circumscribed circle of
1044 * the triangle
1045 * \retval NEGATIVE if \p l is outside the circumscribed circle of
1046 * the triangle
1047 */
1048 virtual Sign incircle(index_t i,index_t j,index_t k,index_t l) const=0;
1049
1050 /**
1051 * \brief Given two segments that have an intersection, create the
1052 * intersection
1053 * \details The intersection is given both as the indices of segment
1054 * extremities (i,j) and (k,l), that one can use to retreive the
1055 * points in derived classes, and constraint indices E1 and E2, that
1056 * derived classes may use to retreive symbolic information attached
1057 * to the constraint
1058 * \param[in] E1 the index of the first edge, corresponding to the
1059 * value of ncnstr() when insert_constraint() was called for
1060 * that edge
1061 * \param[in] i , j the vertices of the first segment
1062 * \param[in] E2 the index of the second edge, corresponding to the
1063 * value of ncnstr() when insert_constraint() was called for
1064 * that edge
1065 * \param[in] k , l the vertices of the second segment
1066 * \return the index of a newly created vertex that corresponds to
1067 * the intersection between [\p i , \p j] and [\p k , \p l]
1068 */
1069 virtual index_t create_intersection(
1070 index_t E1, index_t i, index_t j,
1071 index_t E2, index_t k, index_t l
1072 ) = 0;
1073
1074 /**
1075 * \brief Finds the index of an integer in an array of three integers.
1076 * \param[in] T a const pointer to an array of three integers
1077 * \param[in] v the integer to retrieve in \p T
1078 * \return the index (0,1 or 2) of \p v in \p T
1079 * \pre The three entries of \p T are different and one of them is
1080 * equal to \p v.
1081 */
1082 23779679 static inline index_t find_3(const index_t* T, index_t v) {
1083 // The following expression is 10% faster than using
1084 // if() statements. This uses the C++ norm, that
1085 // ensures that the 'true' boolean value converted to
1086 // an int is always 1. With most compilers, this avoids
1087 // generating branching instructions.
1088 // Thank to Laurent Alonso for this idea.
1089
2/2
✓ Branch 0 taken 8923859 times.
✓ Branch 1 taken 14855820 times.
 23779679 index_t result = index_t( (T[1] == v) | ((T[2] == v) * 2) );
1090 // Sanity check, important if it was T[0], not explicitly
1091 // tested (detects input that does not meet the precondition).
1092
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 23779679 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 23779679 geo_debug_assert(T[result] == v);
1093 23779679 return result;
1094 }
1095
1096 /*******************************************************************/
1097
1098 /**
1099 * \brief Removes all the triangles that have the flag T_MARKED_FLAG
1100 * set.
1101 * \details This compresses the triangle array, and updates triangle
1102 * adjacencies, as well as the vertex to triangle array. Note that
1103 * it uses Tnext_'s storage for internal bookkeeping, so this function
1104 * should not be used if there exists a non-empty DList.
1105 */
1106 void remove_marked_triangles();
1107
1108 /******************** Debugging ************************************/
1109
1110 /**
1111 * \brief Consistency check for a triangle
1112 * \details aborts if inconsistency is detected
1113 * \param[in] t the triangle to be tested
1114 */
1115 6210534 void Tcheck(index_t t) const {
1116
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6210534 times.
 6210534 if(t == NO_INDEX) {
1117 ✗ return;
1118 }
1119
2/2
✓ Branch 0 taken 18631602 times.
✓ Branch 1 taken 6210534 times.
 24842136 for(index_t e=0; e<3; ++e) {
1120
1/6
✗ Branch 2 not taken.
✓ Branch 3 taken 18631602 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 18631602 geo_assert(Tv(t,e) != Tv(t,(e+1)%3));
1121
2/2
✓ Branch 1 taken 790114 times.
✓ Branch 2 taken 17841488 times.
 18631602 if(Tadj(t,e) == NO_INDEX) {
1122 790114 continue;
1123 }
1124
1/6
✗ Branch 2 not taken.
✓ Branch 3 taken 17841488 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 17841488 geo_assert(Tadj(t,e) != Tadj(t,(e+1)%3));
1125 17841488 index_t t2 = Tadj(t,e);
1126 17841488 index_t e2 = Tadj_find(t2,t);
1127
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 17841488 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 17841488 geo_assert(Tadj(t2,e2) == t);
1128 }
1129 }
1130
1131 /**
1132 * \brief Consistency check for a triangle in debug mode,
1133 * ignored in release mode
1134 * \details aborts if inconsistency is detected
1135 * \param[in] t the triangle to be tested
1136 */
1137 2283724 void debug_Tcheck(index_t t) const {
1138 #ifdef GEO_DEBUG
1139 2283724 Tcheck(t);
1140 #else
1141 geo_argused(t);
1142 #endif
1143 2283724 }
1144
1145 /**
1146 * \brief Consistency combinatorial check for all the triangles
1147 * \details aborts if inconsistency is detected
1148 */
1149 93330 void check_combinatorics() const {
1150
2/2
✓ Branch 1 taken 3926810 times.
✓ Branch 2 taken 93330 times.
 4020140 for(index_t t=0; t<nT(); ++t) {
1151 3926810 Tcheck(t);
1152 }
1153 93330 }
1154
1155 /**
1156 * \brief Consistency combinatorial check for all the triangles
1157 * in debug mode, ignored in release mode
1158 * \details aborts if inconsistency is detected
1159 */
1160 93266 void debug_check_combinatorics() const {
1161 #ifdef GEO_DEBUG
1162 93266 check_combinatorics();
1163 #endif
1164 93266 }
1165
1166 /**
1167 * \brief Consistency geometrical check for all the triangles
1168 * \details aborts if inconsistency is detected
1169 */
1170 virtual void check_geometry() const;
1171
1172 /**
1173 * \brief Consistency geometrical check for all the triangles
1174 * in debug mode, ignored in release mode
1175 * \details aborts if inconsistency is detected
1176 */
1177 93266 void debug_check_geometry() const {
1178 #ifdef GEO_DEBUG
1179 93266 check_geometry();
1180 #endif
1181 93266 }
1182
1183
1184 public:
1185 /**
1186 * \brief Checks both combinatorics and geometry,
1187 * aborts on unconsistency
1188 */
1189 64 void check_consistency() const {
1190 64 check_combinatorics();
1191 64 check_geometry();
1192 64 }
1193
1194 protected:
1195 /**
1196 * \brief Checks both combinatorics and geometry
1197 * in debug mode, ignored in release mode,
1198 * aborts on unconsistency
1199 */
1200 93266 void debug_check_consistency() const {
1201 93266 debug_check_combinatorics();
1202 93266 debug_check_geometry();
1203 93266 }
1204
1205 /**
1206 * \brief Tests whether two segments have a frank intersection
1207 * \param[in] u1 , u2 the two extremities of the first segment
1208 * \param[in] v1 , v2 the two extremities of the second segment
1209 * \retval true if \p u1 , \p u2 has a frank intersection with
1210 * \p v1 , \p v2
1211 * \retval false otherwise
1212 */
1213 95129 bool segment_segment_intersect(
1214 index_t u1, index_t u2, index_t v1, index_t v2
1215 ) const {
1216
2/2
✓ Branch 2 taken 1154 times.
✓ Branch 3 taken 93975 times.
 95129 if(orient2d(u1,u2,v1)*orient2d(u1,u2,v2) > 0) {
1217 1154 return false;
1218 }
1219 93975 return (orient2d(v1,v2,u1)*orient2d(v1,v2,u2) < 0);
1220 }
1221
1222 /**
1223 * \brief Tests whether an edge triangle and a segment have a frank
1224 * intersection
1225 * \param[in] v1 , v2 the two extremities of the segment
1226 * \param[in] t a triangle
1227 * \param[in] le local edge index in 0,1,2
1228 * \retval true if edge \p le of \p t has a frank intersection with edge
1229 * \p v1 , \p v2
1230 * \retval false otherwise
1231 */
1232 95129 bool segment_edge_intersect(
1233 index_t v1, index_t v2, index_t t, index_t le
1234 ) const {
1235 95129 index_t u1 = Tv(t,(le + 1)%3);
1236 95129 index_t u2 = Tv(t,(le + 2)%3);
1237 95129 return segment_segment_intersect(u1,u2,v1,v2);
1238 }
1239
1240 /**
1241 * \brief Checks that the edges stored in a DList exactly correspond
1242 * to all edge intersections between a segment and the triangle edges
1243 * \param[in] v1 , v2 the two vertices of the constrained segment
1244 * \param[in] Q a list of triangle. For each triangle in Q, edge 0
1245 * is supposed to have an intersection with \p v1 , \p v2
1246 */
1247 void check_edge_intersections(
1248 index_t v1, index_t v2, const DList& Q
1249 );
1250
1251 typedef std::pair<index_t, index_t> Edge;
1252
1253 /**
1254 * \brief Gets a triangle incident a a given edge
1255 * \param[in] E the edge
1256 * \return a triangle with E as its edge 0
1257 */
1258 ✗ index_t eT(Edge E) {
1259 ✗ index_t v1 = E.first;
1260 ✗ index_t v2 = E.second;
1261 ✗ index_t result = NO_INDEX;
1262 ✗ for_each_T_around_v(
1263 ✗ v1, [&](index_t t, index_t lv)->bool {
1264 ✗ if(Tv(t, (lv+1)%3) == v2) {
1265 ✗ if(Tv(t, (lv+2)%3) != v1) {
1266 ✗ Trot(t, (lv+2)%3);
1267 }
1268 ✗ result = t;
1269 ✗ return true;
1270 ✗ } else if(Tv(t, (lv+1)%3) == v1) {
1271 ✗ if(Tv(t, (lv+2)%3) != v2) {
1272 ✗ Trot(t, (lv+2)%3);
1273 }
1274 ✗ result = t;
1275 ✗ return true;
1276 }
1277 ✗ return false;
1278 }
1279 );
1280 ✗ geo_debug_assert(result != NO_INDEX);
1281 ✗ geo_debug_assert(
1282 (Tv(result,1) == v1 && Tv(result,2) == v2) ||
1283 (Tv(result,1) == v2 && Tv(result,2) == v1)
1284 );
1285 ✗ return result;
1286 }
1287
1288 /**
1289 * \brief Simpler version of locate() kept for reference
1290 * \see locate()
1291 */
1292 index_t locate_naive(
1293 index_t v, index_t hint = NO_INDEX, Sign* orient = nullptr
1294 ) const;
1295
1296 /**
1297 * \brief Simpler version of constrain_edges() kept for reference
1298 * \see constrain_edges()
1299 */
1300 void constrain_edges_naive(
1301 index_t i, index_t j, DList& Q, vector<Edge>& N
1302 );
1303
1304 /**
1305 * \brief Simpler version of Delaunayize_new_edges() that uses a vector
1306 * instead of a DList, kept for reference
1307 */
1308 void Delaunayize_new_edges_naive(vector<Edge>& N);
1309
1310 protected:
1311 index_t nv_;
1312 index_t ncnstr_;
1313 vector<index_t> T_; /**< triangles vertices array */
1314 vector<index_t> Tadj_; /**< triangles adjacency array */
1315 vector<index_t> v2T_; /**< vertex to triangle back pointer */
1316 vector<uint8_t> Tflags_; /**< triangle flags */
1317 vector<index_t> Tecnstr_first_; /**< index in edge constraints list */
1318 vector<index_t> ecnstr_val_; /**< edge constraints list, indices */
1319 vector<index_t> ecnstr_next_; /**< edge constraints list, links */
1320 vector<index_t> Tnext_; /**< doubly connected triangle list */
1321 vector<index_t> Tprev_; /**< doubly connected triangle list */
1322 bool delaunay_; /**< if set, compute a CDT, else just a CT */
1323 Sign orient_012_; /**< global triangles orientation */
1324 bool exact_incircle_; /**< true if incircle() is 100% exact */
1325 bool exact_intersections_; /**< true if intersections are 100% exact */
1326 };
1327
1328 /*****************************************************************/
1329
1330 /**
1331 * \brief Constrained Delaunay triangulation
1332 * \details
1333 * Example:
1334 * \code
1335 * CDT cdt;
1336 * vec2 p1(.,.), p2(.,.), p3(.,.);
1337 * cdt.create_enclosing_triangle(p1,p2,p3);
1338 * // or create_enclosing_quad() or create_enclosing_rect()
1339 * // insert points
1340 * for(...) {
1341 * vec2 p(.,.);
1342 * index_t v = cdt.insert(p);
1343 * ...
1344 * }
1345 * // insert constraints
1346 * for(...) {
1347 * index_t v1=..., v2=...;
1348 * cdt.insert_constraint(v1,v2);
1349 * }
1350 * // get triangles
1351 * for(index_t t=0; t<cdt.nT(); ++t) {
1352 * index_t v1 = cdt.Tv(t,0);
1353 * index_t v2 = cdt.Tv(t,1);
1354 * index_t v3 = cdt.Tv(t,2);
1355 * ... do something with v1,v2,v3
1356 * }
1357 * \endcode
1358 * If some constraints are intersecting, new vertices are generated. They
1359 * can be accessed using the function vec2 CDT::point(index_t v). Vertices
1360 * coming from an intersection are between indices nv1 and CDT::nv(),
1361 * where nv1 is the value of CDT::nv() before inserting the constraints
1362 * (nv1 corresponds to the number of times CDT::insert() was called plus
1363 * the number of points in the enclosing polygon). Note that like input
1364 * points, constraint intersections are represented using double-precision
1365 * floating point numbers, which is not always sufficient to ensure
1366 * robustness. If the input has intersecting constraints and bullet-proof
1367 * guarantees are needed, one can use ExactCDT2d instead.
1368 *
1369 * If you want only a constrained triangulation (not Delaunay),
1370 * you can call CDT::set_Delaunay(false) before inserting the points.
1371 *
1372 * If you have many points to insert, you can use the function:
1373 * \code
1374 * void CDT::insert(index_t nb_points, const double* points);
1375 * \endcode
1376 * It is much much faster than inserting the points one by one. Internally
1377 * it uses Amenta et.al's BRIO method (multi-resolution spatial sort).
1378 */
1379 class GEOGRAM_API CDT2d: public CDTBase2d {
1380 public:
1381
1382 CDT2d();
1383
1384 ~CDT2d() override;
1385
1386 /**
1387 * \copydoc CDTBase2d::clear()
1388 */
1389 void clear() override;
1390
1391 /**
1392 * \brief Creates a first large enclosing triangle
1393 * \param[in] p1 , p2 , p3 the three vertices of the first triangle
1394 * \details create_enclosing_triangle(), create_enclosing_rectangle()
1395 * or create_enclosing_quad() need to be called before anything else.
1396 * Note that create_enclosing_triangle() creates three vertices, with
1397 * indices 0,1,2. Then, the first call to insert() creates vertex 3.
1398 * It is important to know when calling insert_constraint().
1399 */
1400 void create_enclosing_triangle(
1401 const vec2& p1, const vec2& p2, const vec2& p3
1402 );
1403
1404 /**
1405 * \brief Creates a first large enclosing quad
1406 * \param[in] p1 , p2 , p3 , p4 the four vertices of the quad
1407 * \details The quad needs to be convex.
1408 * create_enclosing_triangle(), create_enclosing_rectangle()
1409 * or create_enclosing_quad() need to be called before anything else.
1410 * Note that create_enclosing_quad() creates four vertices, with
1411 * indices 0,1,2,3. Then, the first call to insert() creates vertex 4.
1412 * It is important to know when calling insert_constraint().
1413 */
1414 void create_enclosing_quad(
1415 const vec2& p1, const vec2& p2, const vec2& p3, const vec2& p4
1416 );
1417
1418 /**
1419 * \brief Creates a first large enclosing rectangle
1420 * \param[in] x1 , y1 , x2 , y2 rectangle bounds
1421 * \details create_enclosing_triangle(), create_enclosing_rectangle()
1422 * or create_enclosing_quad() need to be called before anything else.
1423 * Note that create_enclosing_rectangle() creates four vertices, with
1424 * indices 0,1,2,3. Then, the first call to insert() creates vertex 4.
1425 * It is important to know when calling insert_constraint().
1426 */
1427 void create_enclosing_rectangle(
1428 double x1, double y1, double x2, double y2
1429 ) {
1430 create_enclosing_quad(
1431 vec2(x1,y1),
1432 vec2(x2,y1),
1433 vec2(x2,y2),
1434 vec2(x1,y2)
1435 );
1436 }
1437
1438 /**
1439 * \brief Inserts a point
1440 * \param[in] p the point to be inserted
1441 * \param[in] hint a triangle not too far away from the point to
1442 * be inserted
1443 * \return the index of the created point. Duplicated points are
1444 * detected (and then the index of the existing point is returned)
1445 */
1446 index_t insert(const vec2& p, index_t hint = NO_INDEX) {
1447 debug_check_consistency();
1448 point_.push_back(p);
1449 index_t v = CDTBase2d::insert(point_.size()-1, hint);
1450 // If inserted point already existed in
1451 // triangulation, then nv() did not increase
1452 if(point_.size() > nv()) {
1453 point_.pop_back();
1454 }
1455 debug_check_consistency();
1456 return v;
1457 }
1458
1459 /**
1460 * \brief Batch-inserts a set of point
1461 * \details In general, it is much faster than calling
1462 * insert() multiple times. Internally it uses a spatial
1463 * sort (Amenta et.al's BRIO method).
1464 * On exit, the optional \p indices array contains the index mapping.
1465 * indices[i] may be different from i if there were duplicated points.
1466 * If there may be duplicated points and if one wants to insert
1467 * constraint using CDTBase2d::insert_constraint(), one needs to call
1468 * insert_constraint(indices[i],indices[j]) to get the correct
1469 * translation of the indices
1470 * \param[in] points a contiguous array of all point coordinates
1471 * \param[in] nb_points number of points
1472 * \param[out] indices an optional pointer to an array of size \p
1473 * nb_points of indices. On exit, indices[i] contains the index
1474 * of the mesh vertex that corresponds to the i-th point. If there
1475 * are duplicated points or if \p remove_unreferenced_vertices is set,
1476 * then indices[i] may be different from i
1477 * \param[in] remove_unreferenced_vertices if set, then duplicated
1478 * vertices are not stored in the vertices array. Internally, this
1479 * systematically changes the order of the points. For this reason,
1480 * if this flag is set, then one needs to do index mapping with
1481 * \p indices, even when there is no duplicated point
1482 */
1483 void insert(
1484 index_t nb_points, const double* points,
1485 index_t* indices = nullptr,
1486 bool remove_unreferenced_vertices = false
1487 );
1488
1489 /**
1490 * \copydoc CDTBase2d::save()
1491 */
1492 void save(const std::string& filename) const override;
1493
1494 /**
1495 * \brief Gets a point by index
1496 * \param[in] v point index
1497 * \return the point at index \p v
1498 */
1499 168 const vec2 point(index_t v) const {
1500
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(v < nv());
1501 168 return point_[v];
1502 }
1503
1504 protected:
1505 /**
1506 * \copydoc CDTBase2d::orient_2d()
1507 */
1508 Sign orient2d(index_t i, index_t j, index_t k) const override;
1509
1510 /**
1511 * \copydoc CDTBase2d::incircle()
1512 */
1513 Sign incircle(index_t i,index_t j,index_t k,index_t l) const override;
1514
1515 /**
1516 * \copydoc CDTBase2d::create_intersection()
1517 */
1518 index_t create_intersection(
1519 index_t E1, index_t i, index_t j,
1520 index_t E2, index_t k, index_t l
1521 ) override;
1522
1523 protected:
1524 vector<vec2> point_;
1525 };
1526
1527 /*****************************************************************/
1528
1529 /**
1530 * \brief Constrained Delaunay Triangulation with vertices that are
1531 * exact points. Can be used to implement 2D CSG.
1532 * \details Points are represented using exact 2d homogeneous coordinates.
1533 * Unlike CDT2d, this ensures exact representation of constraints
1534 * intersections with guaranteed behavior.
1535 * Under the hood, it inherits CDTBase2d (constrained Delaunay
1536 * triangulation), and redefines orient2d(), incircle2d() and
1537 * create_intersection() using vectors with homogeneous coordinates
1538 * stored as arithmetic expansions (vec2HE) or arbitrary-precision
1539 * floating point numbers (vec2HEx) if compiled with Tessael's
1540 * geogramplus extension package.
1541 * \see CDT2d
1542 */
1543 class GEOGRAM_API ExactCDT2d : public CDTBase2d {
1544 public:
1545 typedef exact::vec2h ExactPoint;
1546
1547 /**
1548 * \brief ExactCDT2d constructor
1549 */
1550 ExactCDT2d();
1551
1552 /**
1553 * \brief ExactCDT2d destructor
1554 */
1555 ~ExactCDT2d() override;
1556
1557 /**
1558 * \copydoc CDTBase2d::clear()
1559 */
1560 void clear() override;
1561
1562 /**
1563 * \brief Inserts a point
1564 * \param[in] p the point to be inserted
1565 * \param[in] hint a triangle not too far away from the point to
1566 * be inserted
1567 * \param[in] id an opaque identifier attached to the vertex that
1568 * can be used by client code for instance to keep relations
1569 * with a mesh. It can be queried using the vertex_id() function.
1570 * \return the index of the created point. Duplicated points are
1571 * detected (and then the index of the existing point is returned)
1572 */
1573 index_t insert(
1574 const ExactPoint& p, index_t id=0, index_t hint = NO_INDEX
1575 );
1576
1577 /**
1578 * \brief Inserts a constraint
1579 * \param[in] v1 , v2 the two extremities of the constraint, as returned
1580 * by insert()
1581 * \param[in] operand_bits optional bitfield used by 2D CSG, indicating
1582 * on which primitive boundaries the constraint is. Each bit set
1583 * corresponds to a primitive. It is used by the classify() function.
1584 * \see CDTBase::insert_constraint() and classify()
1585 */
1586 53801 void insert_constraint(index_t v1, index_t v2, index_t operand_bits=0) {
1587
2/4
✓ Branch 1 taken 53801 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 53801 times.
✗ Branch 5 not taken.
 53801 constraints_.push_back(bindex(v1,v2,bindex::KEEP_ORDER));
1588 53801 cnstr_operand_bits_.push_back(operand_bits);
1589 53801 CDTBase2d::insert_constraint(v1,v2);
1590 53801 }
1591
1592 /**
1593 * \brief Creates a first large enclosing triangle
1594 * \param[in] p1 , p2 , p3 the three vertices of the triangle
1595 * \details create_enclosing_triangle(), create_enclosing_rectangle() or
1596 * create_enclosing_quad() need to be called before anything else
1597 */
1598 void create_enclosing_triangle(
1599 const ExactPoint& p1, const ExactPoint& p2, const ExactPoint& p3
1600 );
1601
1602 /**
1603 * \brief Creates a first large enclosing quad
1604 * \param[in] p1 , p2 , p3 , p4 the four vertices of the quad
1605 * \details The quad needs to be convex.
1606 * create_enclosing_triangle(), create_enclosing_rectangle() or
1607 * create_enclosing_quad() need to be called before anything else
1608 */
1609 void create_enclosing_quad(
1610 const ExactPoint& p1, const ExactPoint& p2,
1611 const ExactPoint& p3, const ExactPoint& p4
1612 );
1613
1614 /**
1615 * \brief Creates a first large enclosing rectangle
1616 * \param[in] x1 , y1 , x2 , y2 rectangle bounds
1617 * \details create_enclosing_triangle(), create_enclosing_rectangle() or
1618 * create_enclosing_quad() need to be called before anything else
1619 */
1620 8173 void create_enclosing_rectangle(
1621 double x1, double y1, double x2, double y2
1622 ) {
1623
1/2
✓ Branch 1 taken 8173 times.
✗ Branch 2 not taken.
 8173 create_enclosing_quad(
1624
1/2
✓ Branch 2 taken 8173 times.
✗ Branch 3 not taken.
 16346 ExactPoint(vec2(x1,y1)),
1625
1/2
✓ Branch 2 taken 8173 times.
✗ Branch 3 not taken.
 16346 ExactPoint(vec2(x2,y1)),
1626
1/2
✓ Branch 2 taken 8173 times.
✗ Branch 3 not taken.
 16346 ExactPoint(vec2(x2,y2)),
1627
1/2
✓ Branch 2 taken 8173 times.
✗ Branch 3 not taken.
 16346 ExactPoint(vec2(x1,y2))
1628 );
1629 8173 }
1630
1631 /**
1632 * \brief Gets a point by vertex index
1633 * \param[in] v vertex index
1634 * \return the point at index \p v
1635 */
1636 81 const ExactPoint& point(index_t v) const {
1637
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 81 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 81 geo_debug_assert(v < nv());
1638 81 return point_[v];
1639 }
1640
1641 /**
1642 * \brief Gets a vertex id by vertex index
1643 * \param[in] v vertex index
1644 * \return the point at index \p v
1645 */
1646 173481 index_t vertex_id(index_t v) const {
1647
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 173481 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 173481 geo_debug_assert(v < nv());
1648 173481 return id_[v];
1649 }
1650
1651 /**
1652 * \brief Sets a vertex id by vertex index
1653 * \param[in] v vertex index
1654 * \param[in] id vertex id
1655 */
1656 81 void set_vertex_id(index_t v, index_t id) {
1657
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 81 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 81 geo_debug_assert(v < nv());
1658 81 id_[v] = id;
1659 81 }
1660
1661 /**
1662 * \brief Used by 2D CSG operations, discards triangles according
1663 * to a boolean operation.
1664 * \details Discards all the triangles that are outside the object
1665 * defined by the boolean expression. It uses the operand bits
1666 * associated with the constraints.
1667 * \param[in] boolean_expression a string with the boolean expression,
1668 * as defined by BooleanExpression constructor. Each variable
1669 * corresponds to an operand bit associated with the constraints. There
1670 * can be up to 32 operands.
1671 * \param mark_only if set, triangles to be discarded are marked (but
1672 * not discarded)
1673 * \see BooleanExpression, insert_constraint()
1674 */
1675 void classify_triangles(
1676 const std::string& boolean_expression, bool mark_only=false
1677 );
1678
1679 /**
1680 * \see CDTBase2d::save()
1681 */
1682 void save(const std::string& filename) const override;
1683
1684 protected:
1685
1686 /**
1687 * \details Does the same thing as classify("union", mark_only) with
1688 * the difference that cnstr_operand_bits_ contains operand ids
1689 * rather than operand bits. Internally, facet inclusion bits are
1690 * represented by std::set<index_t> (there is a stack of them to
1691 * keep track of inclusion status during traversal). This lets
1692 * compute 2D unions with an arbitrary number of primitives.
1693 * It is used by the projection() primitive of CSGBuilder(), that
1694 * computes a 2D union of all the projected triangles.
1695 */
1696 void classify_triangles_union_cnstr_operand_bits_is_operand_id(
1697 bool mark_only=false
1698 );
1699
1700 void add_point(const ExactPoint& p, index_t id = NO_INDEX);
1701 void begin_insert_transaction() override;
1702 void commit_insert_transaction() override;
1703 void rollback_insert_transaction() override;
1704
1705 /**
1706 * \copydoc CDTBase2d::orient_2d()
1707 */
1708 Sign orient2d(index_t i, index_t j, index_t k) const override;
1709
1710 /**
1711 * \copydoc CDTBase2d::incircle()
1712 */
1713 Sign incircle(index_t i,index_t j,index_t k,index_t l) const override;
1714
1715 /**
1716 * \copydoc CDTBase2d::create_intersection()
1717 */
1718 index_t create_intersection(
1719 index_t E1, index_t i, index_t j,
1720 index_t E2, index_t k, index_t l
1721 ) override;
1722
1723 protected:
1724 vector<ExactPoint> point_;
1725 #ifndef GEOGRAM_USE_EXACT_NT
1726 vector<double> length_;
1727 #endif
1728 vector<index_t> id_;
1729 vector<index_t> cnstr_operand_bits_;
1730 vector<index_t> facet_inclusion_bits_;
1731 mutable std::map<trindex, Sign> pred_cache_;
1732 bool use_pred_cache_insert_buffer_;
1733 mutable std::vector<std::pair<trindex, Sign>> pred_cache_insert_buffer_;
1734 vector<bindex> constraints_;
1735 };
1736
1737 /*****************************************************************/
1738
1739 }
1740
1741 #endif
1742