GCC Code Coverage Report

Directory: ./
File: lib/geogram/delaunay/CDT_2d.h
Date: 2026-09-07 02:28:19
Exec Total Coverage
Lines: 156 184 84.8%
Functions: 16 20 80.0%
Branches: 244 396 61.6%
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 void set_delaunay(bool delaunay) {
123
2/3
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 229 times.
✗ Branch 2 not taken.
 261 delaunay_ = delaunay;
124 }
125
126 /**
127 * \brief Gets the number of triangles
128 */
129 index_t nT() const {
130
28/46
✓ Branch 0 taken 475695 times.
✓ Branch 1 taken 8205 times.
✓ Branch 2 taken 158933 times.
✓ Branch 3 taken 23572 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 7 taken 6534 times.
✓ Branch 8 taken 106 times.
✓ Branch 9 taken 6534 times.
✓ Branch 10 taken 106 times.
✓ Branch 11 taken 498 times.
✓ Branch 12 taken 9 times.
✓ Branch 13 taken 6036 times.
✓ Branch 14 taken 97 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✓ Branch 29 taken 150272 times.
✓ Branch 30 taken 8177 times.
✓ Branch 31 taken 150272 times.
✓ Branch 32 taken 8177 times.
✓ Branch 33 taken 70937 times.
✓ Branch 34 taken 8177 times.
✓ Branch 35 taken 103086 times.
✓ Branch 36 taken 8067 times.
✓ Branch 37 taken 103086 times.
✓ Branch 38 taken 8067 times.
✓ Branch 39 taken 40652 times.
✓ Branch 40 taken 4 times.
✓ Branch 41 taken 286772 times.
✓ Branch 42 taken 222686 times.
✗ Branch 43 not taken.
✓ Branch 44 taken 259998 times.
✗ Branch 45 not taken.
✓ Branch 46 taken 1704732 times.
 4106365 return T_.size()/3;
131 }
132
133 /**
134 * \brief Gets the number of vertices
135 */
136 index_t nv() const {
137
10/18
✓ Branch 0 taken 296 times.
✓ Branch 1 taken 46694 times.
✓ Branch 2 taken 151847 times.
✓ Branch 3 taken 96463 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 32 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 32 times.
✓ Branch 14 taken 84 times.
✓ Branch 15 taken 427 times.
✓ Branch 16 taken 286740 times.
✓ Branch 17 taken 222718 times.
 805439 return nv_;
138 }
139
140 /**
141 * \brief Gets the number of constraints
142 */
143 index_t ncnstr() const {
144 5933 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 index_t Tv(index_t t, index_t lv) const {
154 geo_debug_assert(t<nT());
155 geo_debug_assert(lv<3);
156
23/43
✓ Branch 1 taken 475695 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2672460 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✓ Branch 27 taken 942853 times.
✓ Branch 28 taken 650431 times.
✓ Branch 29 taken 334339 times.
✓ Branch 30 taken 238191 times.
✗ Branch 33 not taken.
✓ Branch 34 taken 23837 times.
✓ Branch 35 taken 2463 times.
✓ Branch 36 taken 521 times.
✓ Branch 38 taken 19367 times.
✓ Branch 39 taken 18922 times.
✓ Branch 40 taken 12779 times.
✓ Branch 41 taken 26178 times.
✓ Branch 42 taken 1543 times.
✓ Branch 43 taken 11236 times.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✓ Branch 46 taken 57602 times.
✓ Branch 47 taken 2347 times.
✓ Branch 48 taken 53729 times.
✓ Branch 49 taken 3873 times.
✓ Branch 50 taken 50005 times.
✓ Branch 51 taken 3724 times.
✓ Branch 53 taken 429424 times.
✓ Branch 54 taken 89137 times.
 14042558 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 index_t Tv_find(index_t t, index_t v) const {
166 geo_debug_assert(t<nT());
167 geo_debug_assert(v<nv());
168 604479 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 index_t Tadj(index_t t, index_t le) const {
179 geo_debug_assert(t<nT());
180 geo_debug_assert(le<3);
181
31/48
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 424 times.
✓ Branch 5 taken 18330 times.
✓ Branch 6 taken 19178 times.
✓ Branch 7 taken 424 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 21830 times.
✓ Branch 13 taken 550700 times.
✓ Branch 15 taken 8127 times.
✓ Branch 16 taken 72448 times.
✓ Branch 18 taken 70925 times.
✓ Branch 19 taken 12 times.
✓ Branch 20 taken 70937 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 70937 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 32268 times.
✓ Branch 25 taken 212542 times.
✓ Branch 26 taken 276990 times.
✓ Branch 27 taken 32268 times.
✓ Branch 28 taken 40396 times.
✗ Branch 29 not taken.
✓ Branch 30 taken 40295 times.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✓ Branch 33 taken 40445 times.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✓ Branch 36 taken 2817626 times.
✓ Branch 37 taken 8933502 times.
✗ Branch 38 not taken.
✓ Branch 39 taken 23837 times.
✓ Branch 40 taken 204097 times.
✓ Branch 41 taken 1490845 times.
✓ Branch 42 taken 38289 times.
✓ Branch 43 taken 18960 times.
✓ Branch 44 taken 15776 times.
✓ Branch 45 taken 62400 times.
✗ Branch 46 not taken.
✓ Branch 47 taken 23858 times.
✓ Branch 48 taken 142919 times.
✓ Branch 49 taken 119001 times.
 23820010 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 index_t Tadj_find(index_t t1, index_t t2) const {
192 geo_debug_assert(t1<nT());
193 geo_debug_assert(t2<nT());
194 8007002 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 index_t vT(index_t v) const {
204 geo_debug_assert(v < nv());
205
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5849 times.
 492854 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 index_t Tedge_cnstr_first(index_t t, index_t le) const {
236 geo_debug_assert(t < nT());
237 geo_debug_assert(le < 3);
238
13/28
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 8127 times.
✓ Branch 13 taken 72448 times.
✓ Branch 15 taken 40396 times.
✓ Branch 16 taken 256 times.
✓ Branch 17 taken 40295 times.
✓ Branch 18 taken 357 times.
✓ Branch 19 taken 40445 times.
✓ Branch 20 taken 207 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✓ Branch 24 taken 23837 times.
✓ Branch 25 taken 9790 times.
✓ Branch 26 taken 1694942 times.
✓ Branch 27 taken 5933 times.
✓ Branch 28 taken 42925 times.
 2314914 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 index_t edge_cnstr_next(index_t ecit) const {
249 66586 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 index_t edge_cnstr(index_t ecit) const {
260
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 8323 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 index_t Tedge_cnstr_nb(index_t t, index_t le) const {
275 index_t result = 0;
276 103086 for(
277 index_t ecit = Tedge_cnstr_first(t,le);
278
3/6
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 32268 times.
✓ Branch 4 taken 31014 times.
✓ Branch 5 taken 70818 times.
 134100 ecit != NO_INDEX;
279 ecit = edge_cnstr_next(ecit)
280 ) {
281 31014 ++result;
282 }
283 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 void Tset_flag(index_t t, index_t flag) {
343 geo_debug_assert(t < nT());
344 geo_debug_assert(flag < 8);
345
3/4
✓ Branch 4 taken 32268 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 43272 times.
✓ Branch 7 taken 27546 times.
 109620 Tflags_[t] |= Numeric::uint8(1u << flag);
346 79335 }
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 void Treset_flag(index_t t, index_t flag) {
354 geo_debug_assert(t < nT());
355 geo_debug_assert(flag < 8);
356 109620 Tflags_[t] &= Numeric::uint8(~(1u << flag));
357 }
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 bool Tflag_is_set(index_t t, index_t flag) {
367 geo_debug_assert(t < nT());
368 geo_debug_assert(flag < 8);
369
8/12
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 6110 times.
✓ Branch 3 taken 13068 times.
✓ Branch 4 taken 79335 times.
✓ Branch 5 taken 70937 times.
✓ Branch 6 taken 70818 times.
✓ Branch 7 taken 206172 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 38289 times.
✓ Branch 11 taken 18960 times.
 606775 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 bool Tis_in_list(index_t t) const {
397 return (
398 (Tflags_[t] &
399 Numeric::uint8((1 << DLIST_NB)-1)
400 ) != 0
401 );
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 193344 DList(CDTBase2d& cdt, index_t list_id) :
423 193344 cdt_(cdt), list_id_(list_id),
424
1/4
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 5849 times.
 14026 back_(NO_INDEX), front_(NO_INDEX) {
425 geo_debug_assert(list_id < DLIST_NB);
426 }
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 542809 DList(CDTBase2d& cdt) :
437 542809 cdt_(cdt), list_id_(NO_INDEX),
438
4/4
✓ Branch 0 taken 178891 times.
✓ Branch 1 taken 427 times.
✓ Branch 2 taken 246202 times.
✓ Branch 3 taken 111356 times.
 536876 back_(NO_INDEX), front_(NO_INDEX) {
439 }
440
441 /**
442 * \brief Initializes a list
443 * \param[in] list_id the DList id, in 0..DLIST_NB-1
444 */
445 void initialize(index_t list_id) {
446 geo_debug_assert(!initialized());
447 geo_debug_assert(list_id < DLIST_NB);
448 425093 list_id_ = list_id;
449 425093 }
450
451 /**
452 * \brief Tests whether a DList is initialized
453 */
454 bool initialized() const {
455
8/8
✓ Branch 0 taken 171752 times.
✓ Branch 1 taken 111164 times.
✓ Branch 2 taken 2193 times.
✓ Branch 3 taken 5934 times.
✓ Branch 4 taken 72257 times.
✓ Branch 5 taken 191 times.
✓ Branch 6 taken 11353 times.
✓ Branch 7 taken 16368 times.
 391212 return (list_id_ != NO_INDEX);
456 }
457
458 736153 ~DList() {
459
2/2
✓ Branch 0 taken 618437 times.
✓ Branch 1 taken 117716 times.
 736153 if(initialized()) {
460 clear();
461 }
462 736153 }
463
464 bool empty() const {
465 geo_debug_assert(initialized());
466 geo_debug_assert(
467 (back_==NO_INDEX)==(front_==NO_INDEX)
468 );
469
10/14
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 6534 times.
✓ Branch 3 taken 106 times.
✓ Branch 4 taken 103086 times.
✓ Branch 5 taken 8067 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4 times.
✓ Branch 8 taken 1704732 times.
✓ Branch 9 taken 252051 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 21353 times.
✓ Branch 12 taken 78602 times.
✓ Branch 13 taken 258522 times.
 2433057 return (back_==NO_INDEX);
470 }
471
472 bool contains(index_t t) const {
473 geo_debug_assert(initialized());
474
2/2
✓ Branch 0 taken 38289 times.
✓ Branch 1 taken 18960 times.
 57249 return cdt_.Tflag_is_set(t, list_id_);
475 }
476
477 index_t front() const {
478 geo_debug_assert(initialized());
479 259998 return front_;
480 }
481
482 index_t back() const {
483 geo_debug_assert(initialized());
484 return back_;
485 }
486
487 index_t next(index_t t) const {
488 geo_debug_assert(initialized());
489 geo_debug_assert(contains(t));
490 23837 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 void clear() {
500
5/10
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 11353 times.
✓ Branch 3 taken 258009 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 3968 times.
✓ Branch 7 taken 5765 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 618437 times.
 897532 for(index_t t=front_; t!=NO_INDEX; t = cdt_.Tnext_[t]) {
501 15321 cdt_.Treset_flag(t,list_id_);
502 }
503 263774 back_ = NO_INDEX;
504
1/2
✓ Branch 2 taken 5765 times.
✗ Branch 3 not taken.
 263774 front_ = NO_INDEX;
505 }
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 1868630 void push_back(index_t t) {
517 geo_debug_assert(initialized());
518 geo_debug_assert(!cdt_.Tis_in_list(t));
519
2/2
✓ Branch 0 taken 400377 times.
✓ Branch 1 taken 1468253 times.
 1868630 cdt_.Tset_flag(t,list_id_);
520
2/2
✓ Branch 0 taken 400377 times.
✓ Branch 1 taken 1468253 times.
 1868630 if(empty()) {
521 400377 back_ = t;
522 400377 front_ = t;
523 400377 cdt_.Tnext_[t] = NO_INDEX;
524 400377 cdt_.Tprev_[t] = NO_INDEX;
525 } else {
526 1468253 cdt_.Tnext_[t] = NO_INDEX;
527 1468253 cdt_.Tnext_[back_] = t;
528 1468253 cdt_.Tprev_[t] = back_;
529 1468253 back_ = t;
530 }
531 1868630 }
532
533 1892954 index_t pop_back() {
534 geo_debug_assert(initialized());
535 geo_debug_assert(!empty());
536 1892954 index_t t = back_;
537
2/2
✓ Branch 0 taken 393969 times.
✓ Branch 1 taken 1498985 times.
 1892954 back_ = cdt_.Tprev_[back_];
538
2/2
✓ Branch 0 taken 393969 times.
✓ Branch 1 taken 1498985 times.
 1892954 if(back_ == NO_INDEX) {
539 geo_debug_assert(front_ == t);
540 393969 front_ = NO_INDEX;
541 } else {
542 1498985 cdt_.Tnext_[back_] = NO_INDEX;
543 }
544 geo_debug_assert(contains(t));
545 1892954 cdt_.Treset_flag(t,list_id_);
546 1892954 return t;
547 }
548
549 39645 void push_front(index_t t) {
550 geo_debug_assert(initialized());
551 geo_debug_assert(!cdt_.Tis_in_list(t));
552
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 39645 times.
 39645 cdt_.Tset_flag(t,list_id_);
553
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 39645 times.
 39645 if(empty()) {
554 ✗ back_ = t;
555 ✗ front_ = t;
556 ✗ cdt_.Tnext_[t] = NO_INDEX;
557 ✗ cdt_.Tprev_[t] = NO_INDEX;
558 } else {
559 39645 cdt_.Tprev_[t] = NO_INDEX;
560 39645 cdt_.Tprev_[front_] = t;
561 39645 cdt_.Tnext_[t] = front_;
562 39645 front_ = t;
563 }
564 39645 }
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 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 2362219 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 geo_debug_assert(t < nT());
747 geo_debug_assert(v1 < nv());
748 geo_debug_assert(v2 < nv());
749 geo_debug_assert(v3 < nv());
750 geo_debug_assert(adj1 < nT() || adj1 == NO_INDEX);
751 geo_debug_assert(adj2 < nT() || adj2 == NO_INDEX);
752 geo_debug_assert(adj3 < nT() || adj3 == NO_INDEX);
753 geo_debug_assert(v1 != v2);
754 geo_debug_assert(v2 != v3);
755 geo_debug_assert(v3 != v1);
756 geo_debug_assert(adj1 != adj2 || adj1 == NO_INDEX);
757 geo_debug_assert(adj2 != adj3 || adj2 == NO_INDEX);
758 geo_debug_assert(adj3 != adj1 || adj3 == NO_INDEX);
759 geo_debug_assert(orient2d(v1,v2,v3) != ZERO);
760 2362219 T_[3*t ] = v1;
761 2362219 T_[3*t+1] = v2;
762 2362219 T_[3*t+2] = v3;
763 2362219 Tadj_[3*t ] = adj1;
764 2362219 Tadj_[3*t+1] = adj2;
765 2362219 Tadj_[3*t+2] = adj3;
766 2362219 Tecnstr_first_[3*t] = e1cnstr;
767 2362219 Tecnstr_first_[3*t+1] = e2cnstr;
768 2362219 Tecnstr_first_[3*t+2] = e3cnstr;
769 2362219 v2T_[v1] = t;
770 2362219 v2T_[v2] = t;
771 2362219 v2T_[v3] = t;
772 2362219 }
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 149973 void Trot(index_t t, index_t lv) {
782 geo_debug_assert(t < nT());
783 geo_debug_assert(lv < 3);
784
2/2
✓ Branch 0 taken 88187 times.
✓ Branch 1 taken 61786 times.
 149973 if(lv != 0) {
785 88187 index_t i = 3*t+lv;
786 88187 index_t j = 3*t+((lv+1)%3);
787 88187 index_t k = 3*t+((lv+2)%3);
788 88187 Tset(
789 t,
790 T_[i], T_[j], T_[k],
791 Tadj_[i], Tadj_[j], Tadj_[k],
792 Tecnstr_first_[i], Tecnstr_first_[j], Tecnstr_first_[k]
793 );
794 }
795 149973 }
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 void Tadj_set(index_t t, index_t le, index_t adj) {
819 geo_debug_assert(t < nT());
820 geo_debug_assert(adj < nT());
821 geo_debug_assert(le < 3);
822 3007666 Tadj_[3*t+le] = adj;
823 }
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
2/2
✓ Branch 0 taken 3007666 times.
✓ Branch 1 taken 308606 times.
 3316272 void Tadj_back_connect(
851 index_t t1, index_t le1, index_t prev_t2_adj_e2
852 ) {
853 geo_debug_assert(t1 < nT());
854 geo_debug_assert(le1 < 3);
855 index_t t2 = Tadj(t1,le1);
856
2/2
✓ Branch 0 taken 3007666 times.
✓ Branch 1 taken 308606 times.
 3316272 if(t2 == NO_INDEX) {
857 return;
858 }
859 index_t le2 = Tadj_find(t2,prev_t2_adj_e2);
860 Tadj_set(t2,le2,t1);
861 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
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 686417 times.
 686417 index_t Tnew() {
869 index_t t = nT();
870 686417 index_t nc = (t+1)*3; // new number of corners
871 686417 T_.resize(nc, NO_INDEX);
872 686417 Tadj_.resize(nc, NO_INDEX);
873 686417 Tecnstr_first_.resize(nc, NO_INDEX);
874 686417 Tflags_.resize(t+1,0);
875 686417 Tnext_.resize(t+1,NO_INDEX);
876 686417 Tprev_.resize(t+1,NO_INDEX);
877 686417 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 void Tset_edge_cnstr_first(
888 index_t t, index_t le, index_t ecit
889 ) {
890 geo_debug_assert(t < nT());
891 geo_debug_assert(le < 3);
892
4/4
✓ Branch 0 taken 2193 times.
✓ Branch 1 taken 5934 times.
✓ Branch 2 taken 72257 times.
✓ Branch 3 taken 191 times.
 223494 Tecnstr_first_[3*t+le] = ecit;
893 3409107 }
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 261920 void Tadd_edge_cnstr(
902 index_t t, index_t le, index_t cnstr_id
903 ) {
904 geo_debug_assert(t < nT());
905 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 261920 for(
914 index_t ecit = Tedge_cnstr_first(t,le);
915
2/2
✓ Branch 0 taken 36580 times.
✓ Branch 1 taken 258522 times.
 295102 ecit != NO_INDEX;
916 ecit = edge_cnstr_next(ecit)
917 ) {
918
2/2
✓ Branch 0 taken 33182 times.
✓ Branch 1 taken 3398 times.
 36580 if(edge_cnstr(ecit) == cnstr_id) {
919 return;
920 }
921 }
922
2/2
✓ Branch 0 taken 255605 times.
✓ Branch 1 taken 2917 times.
 258522 ecnstr_val_.push_back(cnstr_id);
923
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 258522 times.
 258522 ecnstr_next_.push_back(Tedge_cnstr_first(t,le));
924 258522 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 261920 void Tadd_edge_cnstr_with_neighbor(
935 index_t t, index_t le, index_t cnstr_id
936 ) {
937 geo_debug_assert(t < nT());
938 geo_debug_assert(le < 3);
939 #ifdef GEO_DEBUG
940 index_t t_e_cnstr_first = Tedge_cnstr_first(t,le);
941 #endif
942 261920 Tadd_edge_cnstr(t, le, cnstr_id);
943 index_t t2 = Tadj(t,le);
944
2/2
✓ Branch 0 taken 142919 times.
✓ Branch 1 taken 119001 times.
 261920 if(t2 != NO_INDEX) {
945 index_t le2 = Tadj_find(t2,t);
946 // Sanity check: make sure the two edges always share the
947 // same constraint list.
948 geo_debug_assert(Tedge_cnstr_first(t2,le2) == t_e_cnstr_first);
949 Tset_edge_cnstr_first(t2,le2,Tedge_cnstr_first(t,le));
950 }
951 261920 }
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 bool Tedge_is_constrained(index_t t, index_t le) const {
961 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 264287 void for_each_T_around_v(
974 index_t v, std::function<bool(index_t t, index_t lv)> doit
975 ) {
976 index_t t = vT(v);
977 index_t lv = NO_INDEX;
978 do {
979 lv = Tv_find(t,v);
980
2/2
✓ Branch 0 taken 216182 times.
✓ Branch 1 taken 202227 times.
 418409 if(doit(t,lv)) {
981 return;
982 }
983
1/2
✓ Branch 0 taken 216182 times.
✗ Branch 1 not taken.
 216182 t = Tadj(t, (lv+1)%3);
984
5/6
✓ Branch 0 taken 216182 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 62060 times.
✓ Branch 3 taken 154122 times.
✓ Branch 4 taken 154122 times.
✓ Branch 5 taken 62060 times.
 278242 } while(t != vT(v) && t != NO_INDEX);
985
986 // We are done, this was an interior vertex
987
1/2
✓ Branch 0 taken 62060 times.
✗ Branch 1 not taken.
 62060 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 t = vT(v);
995 lv = Tv_find(t,v);
996 62060 t = Tadj(t, (lv+2)%3);
997
1/2
✓ Branch 0 taken 100152 times.
✗ Branch 1 not taken.
 100152 while(t != NO_INDEX) {
998 lv = Tv_find(t,v);
999
2/2
✓ Branch 0 taken 38092 times.
✓ Branch 1 taken 62060 times.
 100152 if(doit(t,lv)) {
1000 return;
1001 }
1002 38092 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 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
27/32
✓ Branch 0 taken 1781040 times.
✓ Branch 1 taken 890464 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2671504 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 7 taken 942853 times.
✓ Branch 8 taken 650431 times.
✓ Branch 10 taken 334339 times.
✓ Branch 11 taken 238191 times.
✓ Branch 12 taken 21830 times.
✓ Branch 13 taken 550700 times.
✓ Branch 14 taken 5723 times.
✓ Branch 15 taken 2404 times.
✓ Branch 17 taken 6943 times.
✓ Branch 18 taken 4029 times.
✓ Branch 20 taken 23500 times.
✓ Branch 21 taken 358 times.
✓ Branch 23 taken 252407 times.
✓ Branch 24 taken 166002 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 418409 times.
✓ Branch 27 taken 30171 times.
✓ Branch 28 taken 31889 times.
✓ Branch 29 taken 84393 times.
✓ Branch 30 taken 15759 times.
✗ Branch 31 not taken.
✓ Branch 32 taken 100152 times.
✓ Branch 33 taken 79758 times.
✓ Branch 34 taken 63161 times.
✓ Branch 35 taken 1768863 times.
✓ Branch 36 taken 1238803 times.
 13921471 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 geo_debug_assert(T[result] == v);
1093 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 890820 void Tcheck(index_t t) const {
1116
1/2
✓ Branch 0 taken 890820 times.
✗ Branch 1 not taken.
 890820 if(t == NO_INDEX) {
1117 return;
1118 }
1119
2/2
✓ Branch 0 taken 2672460 times.
✓ Branch 1 taken 890820 times.
 3563280 for(index_t e=0; e<3; ++e) {
1120
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2672460 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2672460 geo_assert(Tv(t,e) != Tv(t,(e+1)%3));
1121
2/2
✓ Branch 0 taken 956 times.
✓ Branch 1 taken 2671504 times.
 2672460 if(Tadj(t,e) == NO_INDEX) {
1122 956 continue;
1123 }
1124
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2671504 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2671504 geo_assert(Tadj(t,e) != Tadj(t,(e+1)%3));
1125 index_t t2 = Tadj(t,e);
1126 index_t e2 = Tadj_find(t2,t);
1127
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2671504 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2671504 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 void debug_Tcheck(index_t t) const {
1138 #ifdef GEO_DEBUG
1139 Tcheck(t);
1140 #else
1141 geo_argused(t);
1142 #endif
1143 }
1144
1145 /**
1146 * \brief Consistency combinatorial check for all the triangles
1147 * \details aborts if inconsistency is detected
1148 */
1149 64 void check_combinatorics() const {
1150
2/2
✓ Branch 0 taken 890820 times.
✓ Branch 1 taken 64 times.
 1781768 for(index_t t=0; t<nT(); ++t) {
1151 890820 Tcheck(t);
1152 }
1153 64 }
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 void debug_check_combinatorics() const {
1161 #ifdef GEO_DEBUG
1162 check_combinatorics();
1163 #endif
1164 }
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 void debug_check_geometry() const {
1178 #ifdef GEO_DEBUG
1179 check_geometry();
1180 #endif
1181 }
1182
1183
1184 public:
1185 /**
1186 * \brief Checks both combinatorics and geometry,
1187 * aborts on unconsistency
1188 */
1189 void check_consistency() const {
1190
2/4
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 32 times.
✗ Branch 5 not taken.
 64 check_combinatorics();
1191
2/4
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 32 times.
✗ Branch 5 not taken.
 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 void debug_check_consistency() const {
1201 debug_check_combinatorics();
1202 debug_check_geometry();
1203 }
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 ✗ bool segment_segment_intersect(
1214 index_t u1, index_t u2, index_t v1, index_t v2
1215 ) const {
1216 ✗ if(orient2d(u1,u2,v1)*orient2d(u1,u2,v2) > 0) {
1217 return false;
1218 }
1219 ✗ 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 ✗ bool segment_edge_intersect(
1233 index_t v1, index_t v2, index_t t, index_t le
1234 ) const {
1235 ✗ index_t u1 = Tv(t,(le + 1)%3);
1236 ✗ index_t u2 = Tv(t,(le + 2)%3);
1237 ✗ 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 const vec2 point(index_t v) const {
1500 geo_debug_assert(v < nv());
1501
1/2
✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
 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 53819 void insert_constraint(index_t v1, index_t v2, index_t operand_bits=0) {
1587 53819 constraints_.push_back(bindex(v1,v2,bindex::KEEP_ORDER));
1588
2/2
✓ Branch 0 taken 52220 times.
✓ Branch 1 taken 1599 times.
 53819 cnstr_operand_bits_.push_back(operand_bits);
1589 53819 CDTBase2d::insert_constraint(v1,v2);
1590 53819 }
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 16346 ExactPoint(vec2(x1,y1)),
1625 16346 ExactPoint(vec2(x2,y1)),
1626 16346 ExactPoint(vec2(x2,y2)),
1627 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 const ExactPoint& point(index_t v) const {
1637 geo_debug_assert(v < nv());
1638 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 index_t vertex_id(index_t v) const {
1647 geo_debug_assert(v < nv());
1648
1/2
✓ Branch 1 taken 30285 times.
✗ Branch 2 not taken.
 57831 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 void set_vertex_id(index_t v, index_t id) {
1657 geo_debug_assert(v < nv());
1658 81 id_[v] = id;
1659 }
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