GCC Code Coverage Report

Directory: ./
File: lib/geogram/voronoi/generic_RVD_polygon.h
Date: 2026-09-07 02:28:19
Exec Total Coverage
Lines: 53 56 94.6%
Functions: 5 21 23.8%
Branches: 39 46 84.8%
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_POLYGON
41 #define GEOGRAM_VORONOI_GENERIC_RVD_POLYGON
42
43 #include <geogram/basic/common.h>
44 #include <geogram/voronoi/generic_RVD_vertex.h>
45 #include <geogram/basic/attributes.h>
46
47 /**
48 * \file geogram/voronoi/generic_RVD_polygon.h
49 * \brief Internal representation of polygons for GenericVoronoiDiagram.
50 * \note This file contains functions and classes used by the internal
51 * implementation of GEO::GenericVoronoiDiagram.
52 * They are not meant to be used directly by client code.
53 */
54
55 namespace GEOGen {
56
57 /**
58 * \brief Internal representation of polygons for GenericVoronoiDiagram.
59 * \details Stores both geometrical and symbolic representations.
60 * \note This is an internal implementation class used by
61 * GEO::RestrictedVoronoiDiagram. It is not meant to be
62 * used directly by client code.
63 */
64
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 954 class Polygon {
65 public:
66
67 /**
68 * \brief Gets the number of vertices.
69 */
70 index_t nb_vertices() const {
71 return index_t(vertex_.size());
72 }
73
74 /**
75 * \brief Gets a vertex by index.
76 * \param[in] i index of the Vertex in this Polygon
77 * \return a const reference to the Vertex at index \p i
78 * \pre \p i < nb_vertices()
79 */
80 const Vertex& vertex(index_t i) const {
81 geo_debug_assert(i < nb_vertices());
82 return vertex_[i];
83 }
84
85 /**
86 * \brief Gets a vertex by index.
87 * \param[in] i index of the Vertex in this Polygon
88 * \return a reference to the Vertex at index \p i
89 * \pre \p i < nb_vertices()
90 */
91 Vertex& vertex(index_t i) {
92 geo_debug_assert(i < nb_vertices());
93 return vertex_[i];
94 }
95
96 /**
97 * \brief Gets the index of the successor of a Vertex.
98 * \param[in] i index of the Vertex in this Polygon
99 * \return the index of the successor of Vertex \p i
100 * \pre \p i < nb_vertices()
101 */
102 index_t next_vertex(index_t i) const {
103 geo_debug_assert(i < nb_vertices());
104 return
105
2/2
✓ Branch 0 taken 243906 times.
✓ Branch 1 taken 121953 times.
 365859 (i == nb_vertices() - 1) ? 0 : (i + 1)
106 ;
107 }
108
109 /**
110 * \brief Gets the index of the predecessor of a Vertex.
111 * \param[in] i index of the Vertex in this Polygon
112 * \return the index of the predecessor of Vertex \p i
113 * \pre \p ii < nb_vertices()
114 */
115 index_t prev_vertex(index_t i) const {
116 geo_debug_assert(i < nb_vertices());
117 return (i == 0) ? (nb_vertices() - 1) : (i - 1);
118 }
119
120 /**
121 * \brief Adds a Vertex to this Polygon.
122 * \param[in] v the vertex to be added. It is copied.
123 * \return the address of the stored vertex.
124 */
125 Vertex* add_vertex(const Vertex& v) {
126 79892214 vertex_.push_back(v);
127 54890292 return &*(vertex_.rbegin());
128 }
129
130 /**
131 * \brief Clears this Polygon.
132 */
133 void clear() {
134 19890855 vertex_.resize(0);
135 }
136
137 /**
138 * \brief Resizes this Polygon.
139 * \param[in] sz new size
140 */
141 void resize(index_t sz) {
142 vertex_.resize(sz);
143 }
144
145 /**
146 * \brief Assigns a mesh facet to this Polygon.
147 * \details The facet from the initial mesh is converted into
148 * the internal geometric/symbolic representation.
149 * \param[in] mesh the mesh from which the facet is copied
150 * \param[in] f the index of the facet in \p mesh
151 * \param[in] symbolic if true, symbolic information is copied
152 * \param[in] vertex_weight a reference to a vertex attribute
153 * that stores weights. If not bound, then 1.0 is used for
154 * the weights.
155 */
156 void initialize_from_mesh_facet(
157 const Mesh* mesh, index_t f, bool symbolic,
158 const GEO::Attribute<double>& vertex_weight
159 );
160
161 /**
162 * \brief Clips a polygon with a plane.
163 * \details Computes the intersection between this Polygon
164 * and the half-space determined by the positive side
165 * of the bisector of segment [i,j] (on the same side as vertex i).
166 *
167 * \param[out] target where to store the intersection
168 * \param[out] target_intersections
169 * where to allocate the generated vertices
170 * \param[in] mesh the input mesh, used by the symbolic information
171 * \param[in] delaunay the Delaunay triangulation
172 * \param[in] i index of one extremity of bisector in \p delaunay
173 * \param[in] j index of the other extremity of the bisector
174 * in \p delaunay
175 * \param[in] exact if true, exact predicates are used.
176 * Implies symbolic.
177 * \param[in] symbolic if true, symbolic representation
178 * of vertices is computed
179 */
180 template <index_t DIM>
181 38943156 void clip_by_plane(
182 Polygon& target, PointAllocator& target_intersections,
183 const Mesh* mesh, const Delaunay* delaunay,
184 index_t i, index_t j,
185 bool exact, bool symbolic
186 ) {
187
2/2
✓ Branch 0 taken 2041999 times.
✓ Branch 1 taken 17429579 times.
 38943156 if(exact) {
188 4083998 clip_by_plane_exact<DIM>(
189 target, target_intersections, mesh, delaunay, i, j
190 );
191 } else {
192 34859158 clip_by_plane_fast<DIM>(
193 target, target_intersections, delaunay, i, j, symbolic
194 );
195 }
196 38943156 }
197
198 /**
199 * \brief Overwrites this Polygon with the contents of another
200 * polygon.
201 * \param[in] rhs a const reference to the polygon to be copied.
202 */
203 void copy(const Polygon& rhs) {
204 vertex_ = rhs.vertex_;
205 ✗ }
206
207 /**
208 * \brief Swaps the contents of this Polygon and another polygon.
209 * \param[in,out] rhs a reference to the Polygon to be swapped with
210 * this one.
211 */
212 void swap(Polygon& rhs) {
213 vertex_.swap(rhs.vertex_);
214 }
215
216 protected:
217 /**
218 * \brief Clips a Polygon with a plane (fast inexact version).
219 * \details Computes the intersection between this Polygon
220 * and the half-space determined by the positive side
221 * of the bisector of segment [i,j] (the side of i).
222 * This version uses a "fused" predicates-constructions
223 * strategy (and reuses the computations from the predicates
224 * to accelerate the constructions).
225 *
226 * \param[out] target where to store the intersection
227 * \param[out] target_intersections
228 * where to allocate the generated vertices
229 * \param[in] delaunay the Delaunay triangulation
230 * \param[in] i index of one extremity of bisector in \p delaunay
231 * \param[in] j index of the other extremity
232 * of the bisector in \p delaunay
233 * \param[in] symbolic if true, symbolic representation
234 * of vertices is computed
235 *
236 * \internal
237 * \note Profiling revealed that this routine is where
238 * the system spends the largest amount of time
239 * (no big surprise...).
240 */
241 template <index_t DIM>
242 34859158 void clip_by_plane_fast(
243 Polygon& target, PointAllocator& target_intersections,
244 const Delaunay* delaunay, index_t i, index_t j,
245 bool symbolic
246 ) const {
247 target.clear();
248
1/2
✓ Branch 0 taken 17429579 times.
✗ Branch 1 not taken.
 34859158 if(nb_vertices() == 0) {
249 return;
250 }
251
252 const double* geo_restrict pi = delaunay->vertex_ptr(i);
253 geo_assume_aligned(pi, geo_dim_alignment(DIM));
254 const double* geo_restrict pj = delaunay->vertex_ptr(j);
255 geo_assume_aligned(pj, geo_dim_alignment(DIM));
256
257 // Compute d = n . m, where n is the
258 // normal vector of the bisector [pi,pj]
259 // and m the middle point of the bisector.
260 geo_decl_aligned(double d);
261 d = 0;
262
2/2
✓ Branch 0 taken 91195854 times.
✓ Branch 1 taken 17429579 times.
 217250866 for(coord_index_t c = 0; c < DIM; ++c) {
263 182391708 d += (pi[c] + pj[c]) * (pi[c] - pj[c]);
264 }
265
266 // The predecessor of the first vertex is the last vertex
267 34859158 index_t prev_k = nb_vertices() - 1;
268 const Vertex* prev_vk = &(vertex(prev_k));
269 const double* geo_restrict prev_pk = prev_vk->point();
270 geo_assume_aligned(prev_pk, geo_dim_alignment(DIM));
271
272 // We compute:
273 // prev_l = prev_vk . n
274 geo_decl_aligned(double prev_l);
275 prev_l = 0.0;
276
2/2
✓ Branch 0 taken 91195854 times.
✓ Branch 1 taken 17429579 times.
 217250866 for(coord_index_t c = 0; c < DIM; ++c) {
277 182391708 prev_l += prev_pk[c] * (pi[c] - pj[c]);
278 }
279
280 // We compute:
281 // side1(pi,pj,q) = sign(2*q.n - n.m) = sign(2*l - d)
282 34859158 GEO::Sign prev_status = GEO::geo_sgn(2.0 * prev_l - d);
283
284
2/2
✓ Branch 0 taken 65912997 times.
✓ Branch 1 taken 17429579 times.
 333370304 for(index_t k = 0; k < nb_vertices(); k++) {
285 const Vertex* vk = &(vertex(k));
286 const double* pk = vk->point();
287
288 // We compute: l = vk . n
289 geo_decl_aligned(double l);
290 l = 0.0;
291
2/2
✓ Branch 0 taken 343952850 times.
✓ Branch 1 taken 65912997 times.
 819731694 for(coord_index_t c = 0; c < DIM; ++c) {
292 687905700 l += pk[c] * (pi[c] - pj[c]);
293 }
294
295 // We compute:
296 // side1(pi,pj,q) = sign(2*q.n - n.m) = sign(2*l - d)
297
2/2
✓ Branch 0 taken 22730215 times.
✓ Branch 1 taken 43182782 times.
 131825994 GEO::Sign status = GEO::geo_sgn(2.0 * l - d);
298
299 // If status of edge extremities differ,
300 // then there is an intersection.
301
2/2
✓ Branch 0 taken 22730215 times.
✓ Branch 1 taken 43182782 times.
 131825994 if(status != prev_status && (prev_status != 0)) {
302 Vertex I;
303 45460430 double* Ipoint = target_intersections.new_item();
304 I.set_point(Ipoint);
305
2/2
✓ Branch 0 taken 588025 times.
✓ Branch 1 taken 22142190 times.
 45460430 if(symbolic) {
306
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 588025 times.
 1176050 if(
307 1176050 !I.sym().intersect_symbolic(
308 prev_vk->sym(), vk->sym(), j
309 )
310 ) {
311 // We encountered a problem. As a workaround,
312 // we copy prev_vk into the result.
313 ✗ I = *prev_vk;
314 }
315 }
316
317 // Compute lambda1 and lambda2, the
318 // barycentric coordinates of the intersection I
319 // in the segment [prev_vk vk]
320 // Note that d and l (used for the predicates)
321 // are reused here.
322 45460430 double denom = 2.0 * (prev_l - l);
323 double lambda1, lambda2;
324
325 // Shit happens ! [Forrest Gump]
326
1/2
✓ Branch 0 taken 22730215 times.
✗ Branch 1 not taken.
 45460430 if(::fabs(denom) < 1e-20) {
327 lambda1 = 0.5;
328 lambda2 = 0.5;
329 } else {
330 45460430 lambda1 = (d - 2.0 * l) / denom;
331 // Note: lambda2 is also given
332 // by (2.0*l2-d)/denom
333 // (but 1.0 - lambda1 is a bit
334 // faster to compute...)
335 45460430 lambda2 = 1.0 - lambda1;
336 }
337 // Compute intersection I by weighting
338 // the edge extremities with the barycentric
339 // coordinates lambda1 and lambda2
340
2/2
✓ Branch 0 taken 116089455 times.
✓ Branch 1 taken 22730215 times.
 277639340 for(coord_index_t c = 0; c < DIM; ++c) {
341 232178910 Ipoint[c] =
342 232178910 lambda1 * prev_pk[c] +
343 232178910 lambda2 * pk[c];
344 }
345 45460430 I.set_weight(
346
2/2
✓ Branch 0 taken 11363689 times.
✓ Branch 1 taken 11366526 times.
 45460430 lambda1 * prev_vk->weight() + lambda2 * vk->weight()
347 );
348
2/2
✓ Branch 0 taken 11363689 times.
✓ Branch 1 taken 11366526 times.
 45460430 if(status > 0) {
349 I.copy_edge_from(*prev_vk);
350 22727378 I.set_adjacent_seed(signed_index_t(j));
351 } else {
352 I.set_flag(INTERSECT);
353 I.set_adjacent_seed(vk->adjacent_seed());
354 }
355 target.add_vertex(I);
356 }
357
2/2
✓ Branch 0 taken 48145567 times.
✓ Branch 1 taken 17767430 times.
 131825994 if(status > 0) {
358 target.add_vertex(*vk);
359 }
360 prev_vk = vk;
361 prev_pk = pk;
362 prev_status = status;
363 prev_k = k;
364 prev_l = l;
365 }
366 }
367
368 /**
369 * \brief Clips a Polygon with a plane (exact version).
370 * \details Computes the intersection between this Polygon
371 * and the half-space determined by the positive side
372 * of the bisector of segment [i,j] (the side of i).
373 * This version uses symbolically perturbed exact predicates.
374 *
375 * \param[out] target where to store the intersection
376 * \param[out] target_intersections
377 * where to allocate the generated vertices
378 * \param[in] mesh the input mesh (used by exact predicates)
379 * \param[in] delaunay the Delaunay triangulation
380 * \param[in] i index of one extremity of bisector in \p delaunay
381 * \param[in] j index of the other extremity of
382 * the bisector in \p delaunay
383 */
384 template <index_t DIM>
385 4083998 void clip_by_plane_exact(
386 Polygon& target, PointAllocator& target_intersections,
387 const Mesh* mesh, const Delaunay* delaunay,
388 index_t i, index_t j
389 ) {
390 target.clear();
391
1/2
✓ Branch 0 taken 2041999 times.
✗ Branch 1 not taken.
 4083998 if(nb_vertices() == 0) {
392 return;
393 }
394
395 const double* pi = delaunay->vertex_ptr(i);
396 const double* pj = delaunay->vertex_ptr(j);
397
398 // The predecessor of the first vertex is the last vertex
399 4083998 index_t prev_k = nb_vertices() - 1;
400 const Vertex* prev_vk = &(vertex(prev_k));
401 4083998 Sign prev_status = side_exact(
402 mesh, delaunay, *prev_vk, pi, pj, DIM
403 );
404
405
2/2
✓ Branch 0 taken 7548471 times.
✓ Branch 1 taken 2041999 times.
 38361880 for(index_t k = 0; k < nb_vertices(); ++k) {
406 const Vertex* vk = &(vertex(k));
407 15096942 Sign status = side_exact(mesh, delaunay, *vk, pi, pj, DIM);
408
409 // If status of edge extremities differ,
410 // there is an intersection.
411
2/2
✓ Branch 0 taken 1013876 times.
✓ Branch 1 taken 6534595 times.
 15096942 if(status != prev_status && (prev_status != 0)) {
412
413 Vertex I;
414
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1013876 times.
 2027752 if(
415 2027752 !I.sym().intersect_symbolic(
416 prev_vk->sym(), vk->sym(), j
417 )
418 ) {
419 // We encountered a problem. As a workaround,
420 // we copy prev_vk into the result.
421 ✗ I = *prev_vk;
422 // geo_assert_not_reached ;
423 // not supposed to happen in exact mode
424 }
425 2027752 I.intersect_geom<DIM>(
426 target_intersections, *prev_vk, *vk, pi, pj
427 );
428
2/2
✓ Branch 0 taken 506938 times.
✓ Branch 1 taken 506938 times.
 2027752 if(status > 0) {
429 I.copy_edge_from(*prev_vk);
430 1013876 I.set_adjacent_seed(signed_index_t(j));
431 } else {
432 I.set_flag(INTERSECT);
433 I.set_adjacent_seed(vk->adjacent_seed());
434 }
435 target.add_vertex(I);
436 }
437
2/2
✓ Branch 0 taken 6744725 times.
✓ Branch 1 taken 803746 times.
 15096942 if(status > 0) {
438 target.add_vertex(*vk);
439 }
440 prev_vk = vk;
441 prev_status = status;
442 prev_k = k;
443 }
444 }
445
446 /**
447 * \brief Returns the position of a point
448 * relative to a bisector (exact version).
449 * \details Position of q relative to the bisector Pi(i,j).
450 * The symbolic representation of q is used. Symbolic
451 * perturbation is applied to degenerate configurations,
452 * therefore ZERO is never returned.
453 * \param[in] mesh the input mesh
454 * \param[in] delaunay the Delaunay triangulation
455 * \param[in] q query point
456 * \param[in] pi one extremity of the bisector
457 * \param[in] pj the other extremity of the bisector
458 * \param[in] dim dimension of the points
459 * \return POSITIVE if q is on pi's side, NEGATIVE otherwise
460 * (ZERO is never encountered thanks to globally coherent
461 * symbolic perturbations).
462 */
463 static Sign side_exact(
464 const Mesh* mesh, const Delaunay* delaunay,
465 const Vertex& q, const double* pi, const double* pj,
466 coord_index_t dim
467 );
468
469 private:
470 GEO::vector<Vertex> vertex_;
471 };
472 }
473
474 #endif
475