GCC Code Coverage Report

Directory: ./
File: lib/geogram/voronoi/generic_RVD_vertex.h
Date: 2026-09-07 02:36:43
Exec Total Coverage
Lines: 241 251 96.0%
Functions: 57 65 87.7%
Branches: 83 198 41.9%
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_VORONOI_GENERIC_RVD_VERTEX
41 #define GEOGRAM_VORONOI_GENERIC_RVD_VERTEX
42
43 #include <geogram/basic/common.h>
44 #include <geogram/mesh/mesh.h>
45 #include <geogram/delaunay/delaunay_nn.h>
46 #include <geogram/basic/assert.h>
47 #include <geogram/basic/process.h>
48 #include <geogram/basic/attributes.h>
49
50 /**
51 * \file geogram/voronoi/generic_RVD_vertex.h
52 * \brief Types and utilities for manipulating vertices in geometric
53 * and symbolic forms in restricted Voronoi diagrams.
54 * \note This file contains functions and classes used by the
55 * internal implementation of GEO::GenericVoronoiDiagram.
56 * Except some special uses, e.g. subclassing GEO::IntegrationSimplex,
57 * they are not meant to be used directly by client code.
58 */
59
60 namespace GEOGen {
61
62 using GEO::Delaunay; /**< \brief type for nD Delaunay triangulation */
63 using GEO::index_t; /**< \brief type for indices (vertex and facet id) */
64 using GEO::signed_index_t; /**< \brief type for indices (can be <0) */
65 using GEO::coord_index_t; /**< \brief type for coordinate indices */
66 using GEO::Sign; /**< \brief type for signs (POSITIVE,ZERO or NEGATIVE) */
67
68 using GEO::Mesh;
69
70 /**
71 * \brief Small_set is similar to std::set, but with fixed
72 * maximum size (and no dynamic memory allocation).
73 * \details Used by GenericVoronoiDiagram to store vertices equations
74 * (represented as plane indices triplets).
75 * \note This is an internal implementation class, not meant to be
76 * used by client code.
77 */
78 template <class T, index_t DIM>
79 class small_set {
80
81 /** \brief This class type */
82 typedef small_set<T, DIM> thisclass;
83
84 public:
85 /** \brief A random access iterator to elements */
86 typedef T* iterator;
87
88 /** \brief A random access iterator to const elements */
89 typedef const T* const_iterator;
90
91 /** \brief Reference to element */
92 typedef T& reference;
93
94 /** \brief Type of the elements */
95 typedef T value_type;
96
97 /**
98 * \brief Constructs an empty small_set.
99 */
100 36342926 small_set() :
101 36342926 size_(0) {
102 36342926 }
103
104 /**
105 * \brief Gets the number of element in this small_set.
106 */
107 1982204 index_t size() const {
108 1982204 return size_;
109 }
110
111 /**
112 * \brief Gets the maximum number of elements that can be
113 * stored in this small_set.
114 */
115 index_t capacity() const {
116 return (index_t) DIM;
117 }
118
119 /**
120 * \brief Gets an iterator to the first element.
121 */
122 9153693 iterator begin() {
123 9153693 return data_;
124 }
125
126 /**
127 * \brief Gets an iterator one position past the last element.
128 */
129 32106444 iterator end() {
130 32106444 return data_ + size_;
131 }
132
133 /**
134 * \brief Gets an iterator one position past the last element
135 * that can be stored.
136 */
137 7277160 iterator end_of_storage() {
138 7277160 return data_ + DIM;
139 }
140
141 /**
142 * \brief Gets a const iterator to the first element..
143 */
144 550690361 const_iterator begin() const {
145 550690361 return data_;
146 }
147
148 /**
149 * \brief Gets a const iterator one position past the last element.
150 */
151 432306914 const_iterator end() const {
152 432306914 return data_ + size_;
153 }
154
155 /**
156 * \brief Gets a const iterator one position past the last element
157 * that can be stored
158 */
159 const_iterator end_of_storage() const {
160 return data_ + DIM;
161 }
162
163 /**
164 * \brief Insert a new element.
165 * \param[in] x a const reference to the element to be inserted
166 * \return an iterator to the inserted element
167 * \note Throws an assertion failure if maximum capacity is reached
168 */
169 3312752 iterator insert(const T& x) {
170 3312752 return insert(x, find_i(x));
171 }
172
173 /**
174 * \brief Inserts a new element at a specified location..
175 * \param[in] x a const reference to the element to be inserted
176 * \param[in] where an iterator to the location where \p x should be
177 * inserted
178 * \return an iterator to the inserted element (\p = where)
179 * \note Throws an assertion failure if maximum capacity is reached
180 */
181 3312752 iterator insert(const T& x, iterator where) {
182
2/2
✓ Branch 1 taken 2462423 times.
✓ Branch 2 taken 850329 times.
 3312752 if(where == end()) {
183 2462423 *where = x;
184 2462423 grow();
185 2462423 return where;
186 }
187
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 850329 times.
 850329 if(*where == x) {
188 ✗ return where;
189 }
190 850329 grow();
191
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 850329 times.
 850329 if(where == end() - 1) {
192 ✗ *where = x;
193 ✗ return where;
194 }
195
2/2
✓ Branch 1 taken 1026204 times.
✓ Branch 2 taken 850329 times.
 1876533 for(iterator i = end() - 1; i != where; i--) {
196
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 1026204 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 1026204 geo_debug_assert(i != begin());
197 1026204 *i = *(i - 1);
198 }
199 850329 *where = x;
200 #ifdef GEO_DEBUG
201
2/2
✓ Branch 2 taken 1473021 times.
✓ Branch 3 taken 850329 times.
 2323350 for(iterator i = begin(); i != end() - 1; ++i) {
202
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1473021 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 1473021 geo_debug_assert(*i < *(i + 1));
203 }
204 #endif
205 850329 return where;
206 }
207
208 /**
209 * \brief Clears this small_set.
210 */
211 1982204 void clear() {
212 1982204 size_ = 0;
213 1982204 }
214
215 /**
216 * \brief Finds an element by value.
217 * \param[in] x a const reference to the value of the element
218 * \return an iterator to the element or end() if not found
219 */
220 iterator find(const T& x) {
221 iterator result = find_i(x);
222 if(*result != x) {
223 result = end();
224 }
225 return result;
226 }
227
228 /**
229 * \brief Finds an element by value.
230 * \param[in] x a const reference to the value of the element
231 * \return a const iterator to the element or end() if not found
232 */
233 const_iterator find(const T& x) const {
234 const_iterator result = find_i(x);
235 if(*result != x) {
236 result = end();
237 }
238 return result;
239 }
240
241 /**
242 * \brief Appends an element to the end of the list.
243 * \param[in] x a const reference to the value of the element
244 * \pre \p x is greater than all the stored elements
245 */
246 3964408 void push_back(const T& x) {
247 #ifdef GEO_DEBUG
248
2/2
✓ Branch 2 taken 1982204 times.
✓ Branch 3 taken 3964408 times.
 5946612 for(iterator i = begin(); i != end(); ++i) {
249
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1982204 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 1982204 geo_debug_assert(*i < x);
250 }
251 #endif
252 3964408 *end() = x;
253 3964408 grow();
254 3964408 }
255
256 /**
257 * \brief Displays the stored elements.
258 */
259 void print(std::ostream& out) const {
260 out << "[ ";
261 for(const_iterator it = begin(); it != end(); ++it) {
262 out << *it << " ";
263 }
264 out << "]";
265 }
266
267 /**
268 * \brief Direct access to an element.
269 * \param[in] i index of the element
270 * \return a reference to the element
271 */
272 T& operator[] (signed_index_t i) {
273 geo_debug_assert(i >= 0);
274 geo_debug_assert(begin() + i < end());
275 return begin()[i];
276 }
277
278 /**
279 * \brief Direct access to an element.
280 * \param[in] i index of the element
281 * \return a const reference to the element
282 */
283 ✗ const T& operator[] (signed_index_t i) const {
284 ✗ geo_debug_assert(i >= 0);
285 ✗ geo_debug_assert(begin() + i < end());
286 ✗ return begin()[i];
287 }
288
289 protected:
290 /**
291 * \brief Increases the size of this small_set.
292 * \details Cannot grow past the maximum size.
293 */
294 7277160 void grow() {
295
1/6
✗ Branch 2 not taken.
✓ Branch 3 taken 7277160 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 7277160 geo_debug_assert(end() != end_of_storage());
296 7277160 size_++;
297 7277160 }
298
299 // Note: maybe we should start from end() instead of begin()
300 // since negative indices are inserted first.
301
302 /**
303 * \brief Finds where an element is or where it should be inserted
304 * from its value.
305 * \param[in] x a const reference to the value of the element
306 * \return an iterator to the location where the element should be
307 * found or inserted
308 */
309 3312752 iterator find_i(const T& x) {
310 3312752 iterator result = begin();
311
6/6
✓ Branch 1 taken 5119081 times.
✓ Branch 2 taken 2462423 times.
✓ Branch 3 taken 4268752 times.
✓ Branch 4 taken 850329 times.
✓ Branch 5 taken 4268752 times.
✓ Branch 6 taken 3312752 times.
 7581504 while(result != end() && *result < x) {
312 4268752 result++;
313 }
314 3312752 return result;
315 }
316
317 /**
318 * \brief Finds where an element should be located from its value.
319 * \param[in] x a const reference to the value of the element
320 * \return a const iterator to the location where the element should be
321 * found.
322 */
323 const_iterator find_i(const T& x) const {
324 const_iterator result = begin();
325 while(result != end() && *result < x) {
326 result++;
327 }
328 return result;
329 }
330
331 protected:
332 T data_[DIM];
333 index_t size_;
334 };
335
336 /**
337 * \brief Displays the contents of a small_set to a std::ostream.
338 */
339 template <class T, index_t DIM>
340 inline std::ostream& operator<< (
341 std::ostream& out,
342 const small_set<T, DIM>& S) {
343 S.print(out);
344 return out;
345 }
346
347 /**
348 * \brief Computes the intersection between two small_set%s.
349 * \param[in] S1 the first set
350 * \param[in] S2 the second set
351 * \param[out] I where to store the intersection
352 */
353 template <class T, index_t DIM1, index_t DIM2, index_t DIM3>
354 1982204 inline void sets_intersect(
355 const small_set<T, DIM1>& S1,
356 const small_set<T, DIM2>& S2,
357 small_set<T, DIM3>& I
358 ) {
359 1982204 I.clear();
360 1982204 auto i1 = S1.begin();
361 1982204 auto i2 = S2.begin();
362
6/6
✓ Branch 1 taken 7343615 times.
✓ Branch 2 taken 1401236 times.
✓ Branch 4 taken 6762647 times.
✓ Branch 5 taken 580968 times.
✓ Branch 6 taken 6762647 times.
✓ Branch 7 taken 1982204 times.
 8744851 while(i1 < S1.end() && i2 < S2.end()) {
363
2/2
✓ Branch 0 taken 1401236 times.
✓ Branch 1 taken 5361411 times.
 6762647 if(*i1 < *i2) {
364 1401236 ++i1;
365 }
366
2/2
✓ Branch 0 taken 1397003 times.
✓ Branch 1 taken 3964408 times.
 5361411 else if(*i2 < *i1) {
367 1397003 ++i2;
368 }
369 else {
370 3964408 I.push_back(*i1);
371 3964408 ++i1;
372 3964408 ++i2;
373 }
374 }
375 1982204 }
376
377 /**
378 * \brief A set of three integers that encodes the
379 * equation of a vertex in GenericVoronoiDiagram.
380 *
381 * \details
382 * - Each positive entry i denotes the bisector of the segment that connects
383 * the center vertex to the i-th vertex (note that the center vertex
384 * needs to be stored elsewhere, but is known when a RVD is used,
385 * since we know which dual cell we are processing).
386 *
387 * - Each negative entry i denotes the i-th face in the boundary TriMesh.
388 * Note: indexing starts with 1 (resp. -1), 0 is kept for error codes.
389 *
390 * - There is some additional information for the following
391 * two configurations:
392 * - boundary vertex: (nb_boundary_facets = 3)
393 * the index of the boundary vertex is returned
394 * by get_boundary_vertex()
395 * - intersection between boundary edge and bisector:
396 * (nb_boundary_facets = 2)
397 * the indices v1,v2 of the extremities of the boundary edges
398 * are obtained by get_boundary_edge(v1,v2)
399 *
400 * Doing so avoids recomputing vertices that we already know
401 * (and avoids numerical problems when the boundary surface has
402 * coplanar (or nearly coplanar) facets).
403 * It also allows using exact predicates (not implemented yet).
404 *
405 * \note This is an internal implementation class, not meant to be
406 * used by client code.
407 */
408 class SymbolicVertex : public small_set<GEO::signed_index_t, 3> {
409
410 /** \brief This class type */
411 typedef SymbolicVertex thisclass;
412
413 /** \brief The base class of this class */
414 typedef small_set<GEO::signed_index_t, 3> baseclass;
415
416 public:
417 /**
418 * \brief Creates an uninitialized SymbolicVertex.
419 */
420 36342926 SymbolicVertex() :
421 36342926 v1_(0),
422 36342926 v2_(0) {
423 36342926 }
424
425 /**
426 * \brief Adds a bisector to the symbolic representation.
427 */
428 1982204 void add_bisector(index_t i) {
429
1/2
✓ Branch 1 taken 1982204 times.
✗ Branch 2 not taken.
 1982204 baseclass::insert(signed_index_t(i) + 1);
430 1982204 }
431
432 /**
433 * \brief Adds a boundary facet to the symbolic representation.
434 */
435 1330548 void add_boundary_facet(index_t i) {
436
1/2
✓ Branch 1 taken 1330548 times.
✗ Branch 2 not taken.
 1330548 baseclass::insert(-signed_index_t(i) - 1);
437 1330548 }
438
439 /**
440 * \brief Gets the number of boundary facets in the
441 * symbolic representation.
442 */
443 81451770 index_t nb_boundary_facets() const {
444 81451770 index_t result = 0;
445 81451770 for(auto it = baseclass::begin();
446
6/6
✓ Branch 1 taken 159560460 times.
✓ Branch 2 taken 11870829 times.
✓ Branch 3 taken 89979519 times.
✓ Branch 4 taken 69580941 times.
✓ Branch 5 taken 89979519 times.
✓ Branch 6 taken 81451770 times.
 171431289 it != baseclass::end() && *it < 0; ++it) {
447 89979519 result++;
448 }
449 81451770 return result;
450 }
451
452 /**
453 * \brief Gets the number of bisectors in the symbolic representation.
454 */
455 123396469 index_t nb_bisectors() const {
456 123396469 index_t result = 0;
457 123396469 for(auto it = baseclass::end() - 1;
458
6/6
✓ Branch 1 taken 362198702 times.
✓ Branch 2 taken 91065075 times.
✓ Branch 3 taken 329867308 times.
✓ Branch 4 taken 32331394 times.
✓ Branch 5 taken 329867308 times.
✓ Branch 6 taken 123396469 times.
 453263777 it != baseclass::begin() - 1 && *it > 0; --it) {
459 329867308 result++;
460 }
461 123396469 return result;
462 }
463
464 /**
465 * \brief Casts a signed_index_t into an (unsigned) index_t.
466 * \details In debug mode, throws an assertion failure
467 * exception whenever \p x is negative.
468 */
469 133401096 static index_t to_unsigned_int(signed_index_t x) {
470
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 133401096 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 133401096 geo_debug_assert(x >= 0);
471 133401096 return (index_t) (x);
472 }
473
474 /**
475 * \brief Gets a bisector
476 * \param[in] i local index of the bisector
477 * \return the index of the Delaunay vertex that corresponds to
478 * the second extremity of the bisector
479 * \pre i < nb_bisectors()
480 */
481 121390690 index_t bisector(signed_index_t i) const {
482
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 121390690 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 121390690 geo_debug_assert(i < (signed_index_t) nb_bisectors());
483 121390690 return to_unsigned_int((baseclass::end()[-1 - i]) - 1);
484 }
485
486 /**
487 * \brief Gets a boundary facet
488 * \param[in] i local index of the boundary facet
489 * \return the index of the mesh facet
490 * \pre i < nb_boundary_facets()
491 */
492 12010406 index_t boundary_facet(signed_index_t i) const {
493
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 12010406 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 12010406 geo_debug_assert(i < (signed_index_t) nb_boundary_facets());
494 12010406 return to_unsigned_int(-(baseclass::begin()[i]) - 1);
495 }
496
497 /**
498 * \brief Tests whether a bisector is present in the
499 * symbolic representation this vertex.
500 * \param[in] i global index of the bisector
501 */
502 bool has_bisector(index_t i) const {
503 return baseclass::find(signed_index_t(i) + 1) != baseclass::end();
504 }
505
506 /**
507 * \brief Tests whether a boundary facet is present in the
508 * symbolic representation of this vertex.
509 * \param[in] i global index of the boundary facet
510 */
511 bool has_boundary_facet(index_t i) const {
512 return baseclass::find(-signed_index_t(i) - 1) != baseclass::end();
513 }
514
515 /**
516 * \brief Gets the global index of the boundary vertex that corresponds
517 * to this vertex.
518 * \pre nb_boundary_facets() == 3
519 */
520 5838957 index_t get_boundary_vertex() const {
521
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 5838957 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 5838957 geo_debug_assert(nb_boundary_facets() == 3);
522
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 5838957 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 5838957 geo_debug_assert(v1_ != 0);
523 5838957 return v1_ - 1;
524 }
525
526 /**
527 * \brief Gets the global indices of the boundary vertices that
528 * define the boundary edge on which this vertex is located.
529 * \param[out] v1 index of the first extremity of the boundary edge
530 * \param[out] v2 index of the second extremity of the boundary edge
531 * \pre nb_boundary_facets() == 2
532 */
533 6357867 void get_boundary_edge(index_t& v1, index_t& v2) const {
534
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 6357867 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 6357867 geo_debug_assert(nb_boundary_facets() == 2);
535
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 6357867 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 6357867 geo_debug_assert(v1_ != 0);
536
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 6357867 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 6357867 geo_debug_assert(v2_ != 0);
537 6357867 v1 = v1_ - 1;
538 6357867 v2 = v2_ - 1;
539 6357867 }
540
541 /**
542 * \brief Sets the boundary vertex on which this vertex is located.
543 * \param[in] v global index of the boundary vertex
544 */
545 443516 void set_boundary_vertex(index_t v) {
546 443516 v1_ = v + 1;
547 443516 v2_ = 0;
548 443516 }
549
550 /**
551 * \brief Sets the boundary edge on which this vertex is located.
552 * \param[in] v1 global index of the first boundary vertex
553 * \param[in] v2 global index of the second boundary vertex
554 */
555 786428 void set_boundary_edge(index_t v1, index_t v2) {
556 786428 v1_ = v1 + 1;
557 786428 v2_ = v2 + 1;
558 786428 }
559
560 /**
561 * \brief Copies a boundary edge from the symbolic representation
562 * of another vertex.
563 */
564 510170 void copy_boundary_edge_from(const thisclass& rhs) {
565
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 510170 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 510170 geo_debug_assert(rhs.nb_boundary_facets() == 2);
566
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 510170 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 510170 geo_debug_assert(rhs.nb_bisectors() == 1);
567
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 510170 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 510170 geo_debug_assert(rhs.v1_ > 0);
568
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 510170 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 510170 geo_debug_assert(rhs.v2_ > 0);
569 510170 v1_ = rhs.v1_;
570 510170 v2_ = rhs.v2_;
571 510170 }
572
573 /**
574 * \brief Computes the symbolic representation of the intersection
575 * between a segment and a bisector.
576 * \details Computes the intersection between
577 * the segment [\p v1, \p v2] and the bisector \p E
578 *
579 * \return false if there was a problem
580 * (happens sometimes in finite precision mode)
581 */
582 1982204 bool intersect_symbolic(
583 const thisclass& v1,
584 const thisclass& v2,
585 index_t E
586 ) {
587
588 // Compute the symbolic representation as the intersection
589 // of three planes.
590 1982204 sets_intersect(v1, v2, *this);
591 // this computes the set of planes that contain
592 // the edge [v1,v2]
593
594 1982204 add_bisector(E); // the intersection is on E.
595
596 // Compute the symbolic representation as intersection between
597 // bisector and boundary edge
598 // (it's redundant and less elegant than the representation
599 // as planes interactions,
600 // but we need this to handle degenerate configurations properly,
601 // and to use exact predicates with original boundary vertices
602 // coordinates).
603
604
2/2
✓ Branch 1 taken 1296598 times.
✓ Branch 2 taken 685606 times.
 1982204 if(nb_boundary_facets() == 2) {
605 // If *this is on the intersection of two boundary facets,
606 // then *this is on
607 // a boundary edge, and we need to retrieve the indices of the
608 // two extremities of this boundary edge.
609
610 1296598 index_t nb1 = v1.nb_boundary_facets();
611 1296598 index_t nb2 = v2.nb_boundary_facets();
612
4/4
✓ Branch 0 taken 984574 times.
✓ Branch 1 taken 312024 times.
✓ Branch 2 taken 786428 times.
✓ Branch 3 taken 198146 times.
 1296598 if(nb1 == 3 && nb2 == 3) {
613 // If v1 and v2 are boundary vertices,
614 // then I is on the boundary
615 // edge that connects v1 and v2
616 786428 set_boundary_edge(
617 v1.get_boundary_vertex(),
618 v2.get_boundary_vertex()
619 );
620
2/2
✓ Branch 0 taken 312024 times.
✓ Branch 1 taken 198146 times.
 510170 } else if(nb1 == 2) {
621
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 312024 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 312024 geo_debug_assert(nb_boundary_facets() == 2);
622 // If v1 is on a boundary edge,
623 // then I is on the same boundary edge as v1
624 312024 copy_boundary_edge_from(v1);
625
1/2
✓ Branch 0 taken 198146 times.
✗ Branch 1 not taken.
 198146 } else if(nb2 == 2) {
626
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 198146 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 198146 geo_debug_assert(nb_boundary_facets() == 2);
627 // If v2 is on a boundary edge,
628 // then I is on the same boundary edge as v2
629 198146 copy_boundary_edge_from(v2);
630 }
631 }
632
633 // Sanity check: problem detected here, we
634 // notify the caller that will use a workaround
635 // (see clip_by_plane())
636
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1982204 times.
 1982204 if(baseclass::size() != 3) {
637 ✗ return false;
638 }
639 1982204 return true;
640 }
641
642 private:
643 index_t v1_;
644 index_t v2_;
645 };
646
647
648 /**
649 * \brief An allocator for points that are created
650 * from intersections in GenericVoronoiDiagram.
651 *
652 * \details Implementation is an array of chunk. We do not use
653 * std::deque since we want to control the chunk size,
654 * and we want to clear it without deallocating
655 * memory to avoid many calls to memory allocator (which
656 * would probably slow down the Windows version a lot, there
657 * seems to be a global multithreading lock on malloc()).
658 *
659 * In most cases, only the first chunk is used
660 * (but some degenerate cases may use more). There seems
661 * to be no measurable overhead as compared to a contiguous
662 * array in our scenario.
663 *
664 * \note This is an internal implementation class, not meant to be
665 * used by client code.
666 */
667
668 class PointAllocator {
669 public:
670 /**
671 * \brief Creates a new empty PointAllocator.
672 * \param[in] dim dimension of the points to be allocated
673 */
674 748 PointAllocator(coord_index_t dim) :
675 748 size_(0),
676 748 capacity_(0),
677 748 dimension_(dim) {
678 748 }
679
680 /**
681 * \brief Allocates a new point.
682 * \return a pointer to the coordinates of the new point. Memory
683 * ownership remains to this PointAllocator.
684 */
685 39607406 double* new_item() {
686
2/2
✓ Branch 0 taken 863 times.
✓ Branch 1 taken 39606543 times.
 39607406 if(size_ == capacity_) {
687 863 grow();
688 }
689 39607406 size_++;
690 39607406 return item(size_ - 1);
691 }
692
693 /**
694 * \brief Clears this PointAllocator.
695 */
696 4113392 void clear() {
697 4113392 size_ = 0;
698 4113392 }
699
700 /**
701 * \brief PointAllocator destructor
702 * \details This releases all allocated chunks.
703 */
704 748 ~PointAllocator() {
705
2/2
✓ Branch 1 taken 863 times.
✓ Branch 2 taken 748 times.
 1611 for(index_t c = 0; c < chunks_.size(); c++) {
706 863 GEO::Memory::aligned_free(chunks_[c]);
707 }
708 748 }
709
710 /**
711 * \brief Gets the dimension of the points stored in
712 * this PointAllocator.
713 */
714 ✗ coord_index_t dimension() const {
715 ✗ return dimension_;
716 }
717
718 protected:
719 /**
720 * \brief Constants that determine the size of a chunk.
721 */
722 enum {
723 CHUNK_SHIFT = 8,
724 CHUNK_SIZE = 1 << CHUNK_SHIFT,
725 CHUNK_MASK = CHUNK_SIZE - 1
726 };
727
728 /**
729 * \brief Allocates a new chunk of memory.
730 */
731 863 void grow() {
732 1726 chunks_.push_back(
733
1/2
✓ Branch 1 taken 863 times.
✗ Branch 2 not taken.
 863 reinterpret_cast<double*>(
734 863 GEO::Memory::aligned_malloc(
735 863 index_t(CHUNK_SIZE) * dimension_ * sizeof(double)
736 )
737 )
738 );
739 863 capacity_ += CHUNK_SIZE;
740 863 }
741
742 /**
743 * \brief Gets a pointer to one of the allocated points from its index.
744 * \param[in] i the index of the point in this PointAllocator
745 * \return A pointer to the coordinates of the point
746 */
747 39607406 double* item(index_t i) {
748
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 39607406 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 39607406 geo_debug_assert(i < size_);
749 39607406 return &(chunks_[i >> CHUNK_SHIFT][(i & CHUNK_MASK) * dimension_]);
750 }
751
752 private:
753 index_t size_;
754 index_t capacity_;
755 coord_index_t dimension_;
756 std::vector<double*> chunks_;
757 };
758
759 /**
760 * \brief Flags associated with edges.
761 */
762 enum {
763 ORIGINAL = 1, /**< Edge belongs to the input surface */
764 INTERSECT = 2 /**< Edge was generated by an intersection */
765 };
766
767 /**
768 * \brief A set of EdgeFlags
769 * \details EdgeFlag%s are combined with bitewise or.
770 */
771 typedef index_t EdgeFlags;
772
773 /**
774 * \brief An individual edge flag
775 */
776 typedef index_t EdgeFlag;
777
778 /**
779 * \brief Internal representation of vertices
780 * in GenericVoronoiDiagram.
781 * \details Vertex has both
782 * geometrical and symbolic representations.
783 * \note This is an internal implementation class, not meant to be
784 * used by client code (except in some particular case, such as
785 * subclassing GEO::IntegrationSimplex).
786 */
787 class Vertex {
788
789 /** \brief This class type */
790 typedef Vertex thisclass;
791
792 public:
793 /**
794 * \brief Creates a new Vertex
795 * \param[in] p geometric location at the vertex, shared with caller
796 * \param[in] w weight
797 * \param[in] f facet of the input mesh this Vertex comes from
798 * \param[in] sym symbolic representation
799 */
800 Vertex(
801 const double* p, double w, signed_index_t f,
802 const SymbolicVertex& sym
803 ) :
804 point_(p),
805 weight_(w),
806 f_(f),
807 seed_(-1),
808 sym_(sym),
809 flags_(ORIGINAL) {
810 }
811
812 /**
813 * \brief Creates a new Vertex
814 * \param[in] p geometric location at the vertex, shared with caller
815 * \param[in] w weight
816 * \param[in] f facet of the input mesh this Vertex comes from
817 */
818 1300944 Vertex(const double* p, double w, signed_index_t f) :
819 1300944 point_(p),
820 1300944 weight_(w),
821 1300944 f_(f),
822 1300944 seed_(-1),
823 1300944 flags_(ORIGINAL) {
824 1300944 }
825
826 /**
827 * \brief Creates an uninitialized Vertex.
828 */
829 35041982 Vertex() :
830 35041982 point_(nullptr),
831 35041982 weight_(1.0),
832 35041982 f_(-1),
833 35041982 seed_(-1),
834 35041982 flags_(0) {
835 35041982 }
836
837 /**
838 * \brief Gets the geometric location at this Vertex.
839 * \return a const pointer to the coordinates
840 */
841 605752069 const double* point() const {
842 605752069 return point_;
843 }
844
845 /**
846 * \brief Sets the geometric location at this vertex.
847 * \param[in] p the geometric location, shared with caller
848 */
849 42845250 void set_point(const double* p) {
850 42845250 point_ = p;
851 42845250 }
852
853 /**
854 * \brief Gets Vertex weight.
855 * \details Used by non-uniform centroidal
856 * Voronoi tesselation.
857 */
858 79214812 double weight() const {
859 79214812 return weight_;
860 }
861
862 /**
863 * \brief Sets the vertex weight.
864 * \details Used by non-uniform centroidal
865 * Voronoi tesselation..
866 */
867 42845250 void set_weight(double w) {
868 42845250 weight_ = w;
869 42845250 }
870
871 /**
872 * \brief Gets the adjacent seed.
873 * \return the global index of the adjacent seed
874 */
875 27938169 signed_index_t adjacent_seed() const {
876 27938169 return seed_;
877 }
878
879 /**
880 * \brief Sets the adjacent seed.
881 * \param[in] s the global index of the adjacent seed
882 */
883 24842791 void set_adjacent_seed(signed_index_t s) {
884 24842791 seed_ = s;
885 24842791 }
886
887 /** Symbolic representation */
888
889 /**
890 * \brief Gets the symbolic representation.
891 */
892 198087228 const SymbolicVertex& sym() const {
893 198087228 return sym_;
894 }
895
896 /**
897 * \brief Gets the symbolic representation.
898 */
899 4518230 SymbolicVertex& sym() {
900 4518230 return sym_;
901 }
902
903 /**
904 * \brief Gets the adjacent facet.
905 * \return the global index of the adjacent facet
906 */
907 29470806 signed_index_t adjacent_facet() const {
908 29470806 return f_;
909 }
910
911 /**
912 * \brief Sets the adjacent facet.
913 * \param[in] f the global index of the adjacent facet
914 */
915 12419977 void set_adjacent_facet(signed_index_t f) {
916 12419977 f_ = f;
917 12419977 }
918
919 /**
920 * \brief Implicit conversion that accesses the geometric location.
921 * \details With this implicit conversions, we can have template
922 * arguments for RestrictedVoronoiDiagram that take
923 * const double* as arguments instead of Vertices.
924 * \return a const pointer to the coordinates
925 */
926 47670152 operator const double* () const {
927 47670152 return point_;
928 }
929
930 /**
931 * \brief Clears this Vertex.
932 */
933 void clear() {
934 flags_ = 0;
935 f_ = -1;
936 }
937
938 /**
939 * \brief Sets an EdgeFlag in this Vertex.
940 */
941 12422814 void set_flag(EdgeFlag f) {
942 12422814 flags_ |= f;
943 12422814 }
944
945 /**
946 * \brief Resets an EdgeFlag in this Vertex.
947 */
948 void unset_flag(EdgeFlag f) {
949 flags_ &= ~f;
950 }
951
952 /**
953 * \brief Tests an EdgeFlag in this Vertex.
954 */
955 bool check_flag(EdgeFlag f) const {
956 return (flags_ & f) != 0;
957 }
958
959 /**
960 * \brief Copies adjacent facet and edge flags from another Vertex.
961 */
962 12419977 void copy_edge_from(const Vertex& rhs) {
963 12419977 set_adjacent_facet(rhs.adjacent_facet());
964 12419977 flags_ = rhs.flags_;
965 12419977 }
966
967 /**
968 * \brief Computes the intersection between
969 * a segment and a bisector.
970 * \details Computes the intersection between
971 * the segment [vq1, vq2] and the bisector
972 * of [p1,p2]..
973 * \tparam DIM dimension, specified as a template
974 * argument for efficiency considerations
975 */
976 template <index_t DIM>
977 31520634 void intersect_geom(
978 PointAllocator& target_intersections,
979 const Vertex& vq1, const Vertex& vq2,
980 const double* p1, const double* p2
981 ) {
982 31520634 const double* q1 = vq1.point();
983 31520634 const double* q2 = vq2.point();
984 31520634 double* Ipoint = target_intersections.new_item();
985 31520634 set_point(Ipoint);
986 31520634 double d = 0.0, l1 = 0.0, l2 = 0.0;
987
2/2
✓ Branch 0 taken 85054517 times.
✓ Branch 1 taken 15778491 times.
 201520624 for(coord_index_t c = 0; c < DIM; ++c) {
988 169999990 double n = p1[c] - p2[c];
989 169999990 d -= n * (p2[c] + p1[c]);
990 169999990 l1 += q2[c] * n;
991 169999990 l2 += q1[c] * n;
992 }
993 31520634 d = 0.5 * d;
994 31520634 l1 = ::fabs(l1 + d);
995 31520634 l2 = ::fabs(l2 + d);
996 31520634 double l12 = l1 + l2;
997
2/2
✓ Branch 0 taken 15736831 times.
✓ Branch 1 taken 41660 times.
 31520634 if(l12 > 1e-30) {
998 31437314 l1 /= l12;
999 31437314 l2 /= l12;
1000 } else {
1001 83320 l1 = 0.5;
1002 83320 l2 = 0.5;
1003 }
1004
2/2
✓ Branch 0 taken 85054517 times.
✓ Branch 1 taken 15778491 times.
 201520624 for(coord_index_t c = 0; c < DIM; ++c) {
1005 169999990 Ipoint[c] = l1 * q1[c] + l2 * q2[c];
1006 }
1007 31520634 set_weight(l1 * vq1.weight() + l2 * vq2.weight());
1008 31520634 }
1009
1010 /**
1011 * \brief Computes the side of this vertex relative
1012 * to a bisector.
1013 * \details This version is not exact.
1014 * \param[in] p1 first extremity of the bisector
1015 * \param[in] p2 second extremity of the bisector
1016 * \return POSITIVE if this vertex is on p1's side,
1017 * NEGATIVE if this vertex is on p2's side, and ZERO
1018 * if this vertex is on the bisector of [p1,p2].
1019 */
1020 template <index_t DIM>
1021 41891537 Sign side_fast(
1022 const double* p1, const double* p2
1023 ) const {
1024 41891537 double r = 0.0;
1025
2/2
✓ Branch 0 taken 116516526 times.
✓ Branch 1 taken 20954589 times.
 274871666 for(index_t c = 0; c < DIM; ++c) {
1026 232980129 r += GEO::geo_sqr(p2[c] - point()[c]);
1027 232980129 r -= GEO::geo_sqr(p1[c] - point()[c]);
1028 }
1029
1/2
✓ Branch 1 taken 20954589 times.
✗ Branch 2 not taken.
 83783074 return GEO::geo_sgn(r);
1030 }
1031
1032 private:
1033 const double* point_;
1034 double weight_;
1035
1036 /**
1037 * The facet adjacent to the edge
1038 * incident to this vertex.
1039 */
1040 signed_index_t f_;
1041
1042 /**
1043 * indicates the seed of the bisector that generated the
1044 * edge that has this vertex and the previous one as
1045 * extremities (or -1 if border).
1046 */
1047 signed_index_t seed_;
1048
1049 /** The symbolic representation of this vertex. */
1050 SymbolicVertex sym_;
1051
1052 /**
1053 * Indicates the type of edge
1054 * (virtual, original or intersection).
1055 */
1056 EdgeFlags flags_;
1057 };
1058 }
1059
1060 #endif
1061