GCC Code Coverage Report

Directory: ./
File: lib/geogram/voronoi/generic_RVD_cell.cpp
Date: 2026-09-07 02:25:23
Exec Total Coverage
Lines: 85 111 76.6%
Functions: 4 6 66.7%
Branches: 58 105 55.2%
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 #include <geogram/voronoi/generic_RVD_cell.h>
41 #include <geogram/mesh/mesh_halfedges.h>
42 #include <geogram/numerics/predicates.h>
43
44 namespace GEOGen {
45
46 index_t ConvexCell::plus1mod3_[3] = {1, 2, 0};
47 index_t ConvexCell::minus1mod3_[3] = {2, 0, 1};
48
49 ✗ std::ostream& ConvexCell::show_stats(std::ostream& os) const {
50 unsigned int nb_free = 0;
51 unsigned int nb_conflict = 0;
52 unsigned int nb_used = 0;
53 ✗ for(unsigned int t = 0; t < max_t(); t++) {
54 ✗ switch(triangles_[t].status_) {
55 ✗ case TRI_IS_FREE:
56 ✗ nb_free++;
57 ✗ break;
58 ✗ case TRI_IS_USED:
59 ✗ nb_used++;
60 ✗ break;
61 ✗ case TRI_IS_CONFLICT:
62 ✗ nb_conflict++;
63 ✗ break;
64 }
65 }
66 return os << "Nb tot = " << max_t()
67 << " free=" << nb_free
68 << " used=" << nb_used
69 << " conflict=" << nb_conflict
70 ✗ << std::endl;
71 }
72
73 42058330 Sign ConvexCell::side_exact(
74 const Mesh* mesh, const Delaunay* delaunay,
75 const GEOGen::Vertex& q,
76 const double* pi, const double* pj,
77 coord_index_t dim,
78 bool symbolic_is_surface
79 ) const {
80
81
4/5
✓ Branch 0 taken 30284273 times.
✓ Branch 1 taken 10619305 times.
✓ Branch 2 taken 891685 times.
✓ Branch 3 taken 263067 times.
✗ Branch 4 not taken.
 42058330 switch(q.sym().nb_boundary_facets()) {
82 case 0:
83 {
84 // The point q is the intersection between
85 // three bisectors [pi b0], [pi b1] and [pi b2]
86 // (and a tet [q0 q1 q2 q3])
87
88 index_t b0 = q.sym().bisector(0);
89 index_t b1 = q.sym().bisector(1);
90 index_t b2 = q.sym().bisector(2);
91
92
2/2
✓ Branch 0 taken 30218722 times.
✓ Branch 1 taken 65551 times.
 30284273 if(dim == 3) {
93 // 3d is a special case for side4()
94 // (intrinsic dim == ambient dim)
95 // therefore embedding tet q0,q1,q2,q3 is not needed.
96 30218722 return GEO::PCK::side4_3d_SOS(
97 pi,
98 delaunay->vertex_ptr(b0),
99 delaunay->vertex_ptr(b1),
100 delaunay->vertex_ptr(b2),
101 pj
102 30218722 );
103 } else {
104 geo_debug_assert(cell_id() >= 0);
105 65551 index_t t = index_t(cell_id());
106 65551 return GEO::PCK::side4_SOS(
107 pi,
108 delaunay->vertex_ptr(b0),
109 delaunay->vertex_ptr(b1),
110 delaunay->vertex_ptr(b2),
111 pj,
112 mesh->vertices.point_ptr(mesh->cells.tet_vertex(t, 0)),
113 mesh->vertices.point_ptr(mesh->cells.tet_vertex(t, 1)),
114 mesh->vertices.point_ptr(mesh->cells.tet_vertex(t, 2)),
115 mesh->vertices.point_ptr(mesh->cells.tet_vertex(t, 3)),
116 dim
117 65551 );
118 }
119 }
120
121 case 1:
122 {
123 // The point q is the intersection between
124 // a facet (f0,f1,f2) of the surface and two
125 // bisectors [pi b0] and [pi b1].
126
127 index_t b0 = q.sym().bisector(0);
128 index_t b1 = q.sym().bisector(1);
129 index_t f = q.sym().boundary_facet(0);
130
131
2/2
✓ Branch 0 taken 10040147 times.
✓ Branch 1 taken 579158 times.
 10619305 if(symbolic_is_surface) {
132 index_t c = mesh->facets.corners_begin(f);
133 const double* q0 = mesh->vertices.point_ptr(
134 mesh->facet_corners.vertex(c)
135 );
136 10040147 const double* q1 = mesh->vertices.point_ptr(
137 mesh->facet_corners.vertex(c+1)
138 );
139 10040147 const double* q2 = mesh->vertices.point_ptr(
140 mesh->facet_corners.vertex(c+2)
141 );
142
143 10040147 return GEO::PCK::side3_SOS(
144 pi,
145 delaunay->vertex_ptr(b0),
146 delaunay->vertex_ptr(b1),
147 pj,
148 q0, q1, q2, dim
149 10040147 );
150
151 } else {
152 579158 index_t t = f / 4;
153 579158 index_t lf = f % 4;
154 index_t j0 = mesh->cells.tet_vertex(
155 t, GEO::MeshCells::local_tet_facet_vertex_index(lf, 0)
156 );
157 index_t j1 = mesh->cells.tet_vertex(
158 t, GEO::MeshCells::local_tet_facet_vertex_index(lf, 1)
159 );
160 index_t j2 = mesh->cells.tet_vertex(
161 t, GEO::MeshCells::local_tet_facet_vertex_index(lf, 2)
162 );
163
164 579158 return GEO::PCK::side3_SOS(
165 pi,
166 delaunay->vertex_ptr(b0),
167 delaunay->vertex_ptr(b1),
168 pj,
169 mesh->vertices.point_ptr(j0),
170 mesh->vertices.point_ptr(j1),
171 mesh->vertices.point_ptr(j2),
172 dim
173 579158 );
174 }
175 }
176
177 case 2:
178 {
179 // The point q is the intersection between
180 // two facets of the surface (i.e. an edge [e0 e1])
181 // and one bisector [pi b0].
182 // i.e. it's a vertex of the surface.
183 index_t b0 = q.sym().bisector(0);
184 index_t e0, e1;
185 q.sym().get_boundary_edge(e0, e1);
186 891685 return GEO::PCK::side2_SOS(
187 pi, delaunay->vertex_ptr(b0), pj,
188 mesh->vertices.point_ptr(e0),
189 mesh->vertices.point_ptr(e1),
190 dim
191 );
192 }
193
194 case 3:
195 {
196 // The point q is the intersection between
197 // three facets of the surface
198 // (i.e. a vertex v0 of the surface).
199 index_t v0 = q.sym().get_boundary_vertex();
200 263067 return GEO::PCK::side1_SOS(
201 pi, pj, mesh->vertices.point_ptr(v0), dim
202 263067 );
203 }
204 }
205 ✗ geo_assert_not_reached;
206 }
207
208 809425 void ConvexCell::initialize_from_mesh_tetrahedron(
209 const Mesh* mesh, index_t t, bool symbolic,
210 const GEO::Attribute<double>& vertex_weight
211 ) {
212 809425 clear();
213
214 index_t v0 = mesh->cells.tet_vertex(t, 0);
215 index_t v1 = mesh->cells.tet_vertex(t, 1);
216 index_t v2 = mesh->cells.tet_vertex(t, 2);
217 index_t v3 = mesh->cells.tet_vertex(t, 3);
218
219 809425 signed_index_t t0 = signed_index_t(mesh->cells.tet_adjacent(t, 0));
220 809425 signed_index_t t1 = signed_index_t(mesh->cells.tet_adjacent(t, 1));
221 809425 signed_index_t t2 = signed_index_t(mesh->cells.tet_adjacent(t, 2));
222 809425 signed_index_t t3 = signed_index_t(mesh->cells.tet_adjacent(t, 3));
223
224 809425 create_vertex();
225 809425 create_vertex();
226 809425 create_vertex();
227 809425 create_vertex();
228
229
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 809425 times.
 809425 set_cell_id(signed_index_t(t));
230
231
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 809425 times.
 809425 set_vertex_id(0, (t0 == signed_index_t(GEO::NO_CELL)) ? 0 : -t0 - 1);
232 809425 set_vertex_id(1, (t1 == signed_index_t(GEO::NO_CELL)) ? 0 : -t1 - 1);
233 809425 set_vertex_id(2, (t2 == signed_index_t(GEO::NO_CELL)) ? 0 : -t2 - 1);
234
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 809425 times.
 809425 set_vertex_id(3, (t3 == signed_index_t(GEO::NO_CELL)) ? 0 : -t3 - 1);
235
236 double w0 = 1.0;
237 double w1 = 1.0;
238 double w2 = 1.0;
239 double w3 = 1.0;
240
241 if(vertex_weight.is_bound()) {
242 ✗ w0 = vertex_weight[v0];
243 ✗ w1 = vertex_weight[v1];
244 ✗ w2 = vertex_weight[v2];
245 ✗ w3 = vertex_weight[v3];
246 }
247
248 create_triangle(mesh->vertices.point_ptr(v0), w0, 2, 1, 3, 2, 1, 3);
249 create_triangle(mesh->vertices.point_ptr(v1), w1, 3, 0, 2, 3, 0, 2);
250 create_triangle(mesh->vertices.point_ptr(v2), w2, 0, 3, 1, 0, 3, 1);
251 create_triangle(mesh->vertices.point_ptr(v3), w3, 2, 0, 1, 2, 0, 1);
252
253
2/2
✓ Branch 0 taken 19547 times.
✓ Branch 1 taken 789878 times.
 809425 if(symbolic) {
254
255 19547 index_t f0 = global_facet_id(mesh, t, 0);
256 19547 index_t f1 = global_facet_id(mesh, t, 1);
257 19547 index_t f2 = global_facet_id(mesh, t, 2);
258 19547 index_t f3 = global_facet_id(mesh, t, 3);
259
260 triangle_dual(0).sym().set_boundary_vertex(v0);
261 triangle_dual(0).sym().add_boundary_facet(f1);
262 triangle_dual(0).sym().add_boundary_facet(f2);
263 triangle_dual(0).sym().add_boundary_facet(f3);
264
265 triangle_dual(1).sym().set_boundary_vertex(v1);
266 triangle_dual(1).sym().add_boundary_facet(f2);
267 triangle_dual(1).sym().add_boundary_facet(f3);
268 triangle_dual(1).sym().add_boundary_facet(f0);
269
270 triangle_dual(2).sym().set_boundary_vertex(v2);
271 triangle_dual(2).sym().add_boundary_facet(f3);
272 triangle_dual(2).sym().add_boundary_facet(f0);
273 triangle_dual(2).sym().add_boundary_facet(f1);
274
275 triangle_dual(3).sym().set_boundary_vertex(v3);
276 triangle_dual(3).sym().add_boundary_facet(f0);
277 triangle_dual(3).sym().add_boundary_facet(f1);
278 triangle_dual(3).sym().add_boundary_facet(f2);
279 }
280 809425 }
281
282 18 void ConvexCell::initialize_from_surface_mesh(
283 Mesh* mesh, bool symbolic
284 ) {
285 18 clear();
286
287
2/2
✓ Branch 0 taken 108 times.
✓ Branch 1 taken 18 times.
 126 for(index_t f = 0; f < mesh->facets.nb(); ++f) {
288 index_t v = create_vertex();
289 108 set_vertex_id(v,-1-signed_index_t(f));
290 }
291 GEO::vector<GEO::MeshHalfedges::Halfedge> v2h(mesh->vertices.nb());
292
293 GEO::MeshHalfedges MH(*mesh);
294
2/2
✓ Branch 0 taken 108 times.
✓ Branch 1 taken 18 times.
 126 for(index_t f = 0; f < mesh->facets.nb(); ++f) {
295 432 for(index_t c = mesh->facets.corners_begin(f);
296
2/2
✓ Branch 0 taken 432 times.
✓ Branch 1 taken 108 times.
 540 c < mesh->facets.corners_end(f); ++c
297 ) {
298 index_t v = mesh->facet_corners.vertex(c);
299 432 v2h[v] = GEO::MeshHalfedges::Halfedge(f, c);
300 }
301 }
302
303
2/2
✓ Branch 0 taken 144 times.
✓ Branch 1 taken 18 times.
 162 for(index_t v = 0; v < mesh->vertices.nb(); ++v) {
304 index_t fi[3];
305 index_t va[3];
306 index_t cur = 0;
307 144 GEO::MeshHalfedges::Halfedge H = v2h[v];
308 do {
309 // All the vertices of the input mesh should be
310 // incident to three facets exactly (this is because
311 // the ConvexCell is represented in dual form).
312
1/8
✗ Branch 0 not taken.
✓ Branch 1 taken 432 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 432 geo_assert(cur < 3);
313 432 fi[cur] = H.facet;
314
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 432 times.
 432 index_t ca = mesh->facets.next_corner_around_facet(
315 H.facet, H.corner
316 );
317 432 va[cur] = mesh->facet_corners.vertex(ca);
318
1/2
✓ Branch 1 taken 432 times.
✗ Branch 2 not taken.
 432 bool ok = MH.move_to_prev_around_vertex(H);
319
1/8
✗ Branch 0 not taken.
✓ Branch 1 taken 432 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 432 geo_assert(ok);
320
2/2
✓ Branch 0 taken 144 times.
✓ Branch 1 taken 288 times.
 432 ++cur;
321
2/2
✓ Branch 0 taken 288 times.
✓ Branch 1 taken 144 times.
 432 } while(H != v2h[v]);
322
323 // Note: va[] order is different, because of
324 // Mesh numbering -> Triangulation numbering
325 // conversion !
326
1/2
✓ Branch 1 taken 144 times.
✗ Branch 2 not taken.
 144 create_triangle(
327 mesh->vertices.point_ptr(v), 1.0,
328 fi[0], fi[1], fi[2], va[2], va[0], va[1]
329 );
330
1/2
✓ Branch 0 taken 144 times.
✗ Branch 1 not taken.
 144 if(symbolic) {
331 triangle_dual(v).sym().add_boundary_facet(fi[0]);
332 triangle_dual(v).sym().add_boundary_facet(fi[1]);
333 triangle_dual(v).sym().add_boundary_facet(fi[2]);
334 triangle_dual(v).sym().set_boundary_vertex(v);
335 }
336 }
337
1/2
✓ Branch 0 taken 18 times.
✗ Branch 1 not taken.
 18 if(symbolic) {
338 set_symbolic_is_surface(true);
339 }
340 18 }
341
342
343
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
 18 void ConvexCell::convert_to_mesh(Mesh* mesh, bool copy_symbolic_info) {
344 GEO::vector<index_t> tri_to_v(max_t());
345
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 18 mesh->clear();
346
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 18 mesh->vertices.set_dimension(3);
347
348 index_t cur_v = 0;
349
2/2
✓ Branch 0 taken 13377 times.
✓ Branch 1 taken 18 times.
 26772 for(index_t t = 0; t < max_t(); ++t) {
350
2/2
✓ Branch 0 taken 13266 times.
✓ Branch 1 taken 111 times.
 13377 if(triangle_is_valid(t)) {
351
1/2
✓ Branch 1 taken 13266 times.
✗ Branch 2 not taken.
 13266 mesh->vertices.create_vertex(triangle_dual(t).point());
352 13266 tri_to_v[t] = cur_v;
353 13266 ++cur_v;
354 }
355 }
356 GEO::Attribute<signed_index_t> facet_id;
357
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
 18 if(copy_symbolic_info) {
358 ✗ facet_id.bind(mesh->facets.attributes(), "id");
359 }
360 GEO::vector<index_t> facet_vertices;
361
2/2
✓ Branch 0 taken 26760 times.
✓ Branch 1 taken 18 times.
 53556 for(index_t v = 0; v < max_v(); v++) {
362
1/2
✓ Branch 1 taken 26760 times.
✗ Branch 2 not taken.
 26760 facet_vertices.resize(0);
363 signed_index_t t = vertex_triangle(v);
364
2/2
✓ Branch 0 taken 6669 times.
✓ Branch 1 taken 20091 times.
 26760 if(t != -1) {
365 Corner first_c(
366 index_t(t), find_triangle_vertex(index_t(t), v)
367
2/2
✓ Branch 0 taken 4704 times.
✓ Branch 1 taken 1965 times.
 6669 );
368 6669 Corner c = first_c;
369 do {
370
2/2
✓ Branch 0 taken 39696 times.
✓ Branch 1 taken 102 times.
 39798 facet_vertices.push_back(tri_to_v[c.t]);
371 39798 move_to_next_around_vertex(c);
372
2/2
✓ Branch 0 taken 33129 times.
✓ Branch 1 taken 6669 times.
 39798 } while(c != first_c);
373
374
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 6669 times.
✓ Branch 3 taken 6669 times.
✗ Branch 4 not taken.
 6669 index_t f = mesh->facets.create_polygon(facet_vertices.size());
375
2/2
✓ Branch 0 taken 39798 times.
✓ Branch 1 taken 6669 times.
 46467 for(index_t lv=0; lv<facet_vertices.size(); ++lv) {
376
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 39798 times.
 39798 mesh->facets.set_vertex(f,lv,facet_vertices[lv]);
377 }
378 if(facet_id.is_bound()) {
379 ✗ facet_id[f] = vertex_id(v);
380 }
381 }
382 }
383
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 18 mesh->facets.connect();
384 18 }
385
386 ✗ void ConvexCell::copy(const ConvexCell& rhs) {
387 geo_debug_assert(
388 intersections_.dimension() == rhs.intersections_.dimension()
389 );
390 triangles_ = rhs.triangles_;
391 vertices_ = rhs.vertices_;
392 ✗ first_free_ = rhs.first_free_;
393 ✗ v_to_t_dirty_ = rhs.v_to_t_dirty_;
394 ✗ symbolic_is_surface_ = rhs.symbolic_is_surface_;
395 ✗ cell_id_ = rhs.cell_id_;
396 ✗ }
397
398 }
399