GCC Code Coverage Report

Directory: ./
File: voronoi/generic_RVD.h
Date: 2026-09-27 03:22:43
Exec Total Coverage
Lines: 576 729 79.0%
Functions: 183 1274 14.4%
Branches: 343 714 48.0%
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
41 #define GEOGRAM_VORONOI_GENERIC_RVD
42
43 #include <geogram/basic/common.h>
44 #include <geogram/basic/numeric.h>
45 #include <geogram/voronoi/generic_RVD_utils.h>
46 #include <geogram/voronoi/RVD_callback.h>
47 #include <geogram/numerics/predicates.h>
48 #include <geogram/mesh/index.h>
49 #include <geogram/basic/geometry_nd.h>
50 #include <geogram/basic/process.h>
51 #include <geogram/basic/attributes.h>
52 #include <geogram/basic/argused.h>
53
54 #include <deque>
55 #include <algorithm>
56 #include <iostream>
57
58 /**
59 * \file geogram/voronoi/generic_RVD.h
60 * \brief Generic implementation of restricted Voronoi diagrams.
61 * \note This file contains functions and classes used by the
62 * internal implementation of GEO::GenericVoronoiDiagram.
63 * They are not meant to be used directly by client
64 * code.
65 */
66
67 namespace GEOGen {
68
69 /**
70 * \brief Computes the intersection between a surface (Mesh) and a
71 * Voronoi diagram (dual of a Delaunay).
72 * \details The surface may be embedded in nD
73 * (the Voronoi diagram is then of dimension n).
74 * \note This is an internal implementation class, not meant to
75 * be used directly, use GEO::RestrictedVoronoiDiagram instead.
76 */
77 template <index_t DIM>
78 class RestrictedVoronoiDiagram {
79
80 /** \brief This class type */
81 typedef RestrictedVoronoiDiagram<DIM> thisclass;
82
83 public:
84 /**
85 * \brief Gets the dimension
86 */
87 706033982 static coord_index_t dimension() {
88 706033982 return DIM;
89 }
90
91 /**
92 * \brief Used to allocate the generated points.
93 */
94 typedef GEOGen::PointAllocator PointAllocator;
95
96 /**
97 * \brief Internal representation of vertices.
98 */
99 typedef GEOGen::Vertex Vertex;
100
101 /**
102 * \brief Internal representation of polygons.
103 */
104 typedef GEOGen::Polygon Polygon;
105
106 /**
107 * \brief Internal representation of volumetric cells.
108 */
109 typedef GEOGen::ConvexCell Polyhedron;
110
111 /********************************************************************/
112
113 /**
114 * \brief Constructs a new RestrictedVoronoiDiagram.
115 * \param[in] delaunay the Delaunay triangulation
116 * \param[in] mesh the input mesh
117 */
118 192 RestrictedVoronoiDiagram(
119 Delaunay* delaunay,
120 GEO::Mesh* mesh
121 ) :
122 192 mesh_(mesh),
123 192 delaunay_(delaunay),
124 192 intersections_(DIM),
125 192 symbolic_(false),
126 192 check_SR_(true),
127
2/4
✓ Branch 1 taken 96 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 96 times.
✗ Branch 5 not taken.
 192 exact_(false)
128 {
129 192 delaunay_nn_ = dynamic_cast<GEO::Delaunay_NearestNeighbors*>(
130
2/2
✓ Branch 0 taken 48 times.
✓ Branch 1 taken 48 times.
 192 delaunay_
131 );
132 192 dimension_ = DIM;
133 192 facets_begin_ = UNSPECIFIED_RANGE;
134 192 facets_end_ = UNSPECIFIED_RANGE;
135 192 tets_begin_ = UNSPECIFIED_RANGE;
136 192 tets_end_ = UNSPECIFIED_RANGE;
137 192 connected_components_priority_ = false;
138 192 facet_seed_marking_ = nullptr;
139 192 connected_component_changed_ = false;
140 192 current_connected_component_ = 0;
141 192 cur_stamp_ = NO_INDEX;
142 192 current_facet_ = GEO::max_index_t();
143 192 current_seed_ = GEO::max_index_t();
144 192 current_polygon_ = nullptr;
145 192 current_tet_ = GEO::max_index_t();
146 192 current_polyhedron_ = nullptr;
147 192 }
148
149 /**
150 * \brief Sets traveral priority.
151 * \details If connected_components_priority is set,
152 * then the connected components of the
153 * restricted Voronoi cells will be traversed
154 * one by one.
155 */
156 44 void set_connected_components_priority(bool x) {
157 44 connected_components_priority_ = x;
158 44 }
159
160
161 /**
162 * \brief Tests whether connected components priority is
163 * set.
164 * \details If connected_components_priority is set,
165 * then the connected components of the
166 * restricted Voronoi cells will be traversed
167 * one by one.
168 * \retval true if connected components priority is used.
169 * \retval false otherwise.
170 */
171 ✗ bool connected_components_priority() const {
172 ✗ return connected_components_priority_;
173 }
174
175 /**
176 * \brief Gets the input mesh.
177 */
178 ✗ const GEO::Mesh* mesh() const {
179 ✗ return mesh_;
180 }
181
182 /**
183 * \brief Gets the input mesh.
184 */
185 8 GEO::Mesh* mesh() {
186 8 return mesh_;
187 }
188
189 /**
190 * \brief Gets the Delaunay triangulation.
191 */
192 70006 const Delaunay* delaunay() const {
193 70006 return delaunay_;
194 }
195
196 /**
197 * \brief Gets the Delaunay triangulation.
198 */
199 1280 Delaunay* delaunay() {
200 1280 return delaunay_;
201 }
202
203 /**
204 * \brief Sets the Delaunay triangulation.
205 */
206 96 void set_delaunay(Delaunay* delaunay) {
207 96 delaunay_ = delaunay;
208 96 delaunay_nn_ = dynamic_cast<GEO::Delaunay_NearestNeighbors*>(
209
1/2
✓ Branch 0 taken 48 times.
✗ Branch 1 not taken.
 96 delaunay_
210 );
211 96 }
212
213 /**
214 * \brief Sets the input mesh.
215 */
216 96 void set_mesh(GEO::Mesh* mesh) {
217 96 mesh_ = mesh;
218 96 }
219
220 /**
221 * \brief Sets the facets range.
222 * \details Computations can be restricted to a contiguous facet range.
223 * \param[in] facets_begin first facet in the range.
224 * \param[in] facets_end one position past the last facet in the range.
225 */
226 96 void set_facets_range(index_t facets_begin, index_t facets_end) {
227
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 48 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 96 geo_debug_assert(facets_end >= facets_begin);
228 96 facets_begin_ = facets_begin;
229 96 facets_end_ = facets_end;
230 96 }
231
232 /**
233 * \brief Sets the tetrahedra range.
234 * \details Computations can be restricted to a contiguous
235 * tetrahedra range.
236 * \param[in] tets_begin first tetrahedron in the range.
237 * \param[in] tets_end one position past the last
238 * tetrahedron in the range.
239 */
240 32 void set_tetrahedra_range(index_t tets_begin, index_t tets_end) {
241
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 16 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 32 geo_debug_assert(tets_end >= tets_begin);
242 32 tets_begin_ = tets_begin;
243 32 tets_end_ = tets_end;
244 32 }
245
246 /**
247 * \brief Gets the number of facets in the current range.
248 * \see set_facets_range()
249 */
250 index_t nb_facets_in_range() const {
251 return facets_end_ - facets_begin_;
252 }
253
254 /**
255 * \brief Gets the number of tetrahedra in the current range.
256 * \see set_tetrahedra_range()
257 */
258 index_t nb_tetrahedra_in_range() const {
259 return tets_end_ - tets_begin_;
260 }
261
262 /**
263 * \brief Gets the index of the mesh facet currently processed.
264 * \details Can be used in surfacic traversals (and not volumetric
265 * traversals).
266 */
267 420818 index_t current_facet() const {
268 420818 return current_facet_;
269 }
270
271 /**
272 * \brief Gets the index of the Delaunay vertex currently processed.
273 * \details Can be used in both surfacic traversals and volumetric
274 * traversals.
275 */
276 index_t current_seed() const {
277 return current_seed_;
278 }
279
280 /**
281 * \brief Gets the current polygon.
282 * \details The current polygon corresponds to the
283 * intersection between the current facet
284 * and the Voronoi cell of the current seed. Can be used
285 * in surfacic traversals (and not volumetric traversals).
286 */
287 192516944 const Polygon& current_polygon() const {
288 192516944 return *current_polygon_;
289 }
290
291 /**
292 * \brief Gets the undex of the mesh tetrahedron currently processed.
293 * \details Can be used in volumetric traversals (and not in surfacic
294 * traversals).
295 */
296 index_t current_tet() const {
297 return current_tet_;
298 }
299
300 /**
301 * \brief Gets the current cell.
302 * \details The current cell corresponds to the
303 * intersection between the current tetrahedron
304 * and the Voronoi cell of the current seed.
305 * Can be used in volumetric traversals (and not in
306 * surfacic traversals).
307 */
308 66613104 const Polyhedron& current_polyhedron() const {
309 66613104 return *current_polyhedron_;
310 }
311
312 /**
313 * \brief Sets symbolic mode.
314 * \details If exact mode is active, symbolic mode is enforced.
315 * \param[in] x if set, the symbolic representation of the intersections
316 * are computed.
317 */
318 188 void set_symbolic(bool x) {
319 188 symbolic_ = x;
320 // exact mode requires symbolic mode.
321
2/2
✓ Branch 0 taken 78 times.
✓ Branch 1 taken 16 times.
 188 if(exact_) {
322 156 symbolic_ = true;
323 }
324 188 }
325
326 /**
327 * \brief Tests whether symbolic mode is active.
328 */
329 94 bool symbolic() const {
330 94 return symbolic_;
331 }
332
333 /**
334 * \brief Specifies whether exact predicates should be used.
335 * \details If exact predicates are used, symbolic mode is ensured.
336 * \param[in] x if set, exact predicates are used.
337 */
338 206 void set_exact_predicates(bool x) {
339 206 exact_ = x;
340 // exact mode requires symbolic mode.
341
2/2
✓ Branch 0 taken 55 times.
✓ Branch 1 taken 48 times.
 206 if(exact_) {
342 110 symbolic_ = true;
343 }
344 206 }
345
346 /**
347 * \brief Tests whether exact predicates are used.
348 */
349 96 bool exact_predicates() const {
350 96 return exact_;
351 }
352
353 /**
354 * \brief Specifies whether radius of security should be enforced.
355 */
356 280 void set_check_SR(bool x) {
357 280 check_SR_ = x;
358 280 }
359
360 /**
361 * \brief Tests whether radius of security is enforced.
362 * \retval true if radius of security test is used.
363 * \retval false otherwise.
364 */
365 96 bool check_SR() const {
366 96 return check_SR_;
367 }
368
369 /**
370 * \brief Gets the PointAllocator.
371 * \return a pointer to the PointAllocator, used
372 * to create the new vertices generated by
373 * intersections.
374 */
375 ✗ PointAllocator* point_allocator() {
376 ✗ return &intersections_;
377 }
378
379 protected:
380 /**
381 * \name Adapter classes for surfacic computation
382 * @{
383 */
384
385 /**
386 * \brief Adapter class used internally to implement for_each_polygon()
387 * \details Overrides constness checks, to allow using temporaries as
388 * argument of for_each_xxx().
389 * \tparam ACTION the user action class.
390 */
391 template <class ACTION>
392 class PolygonAction {
393 public:
394 /**
395 * \brief Creates a new PolygonAction around a user ACTION instance.
396 * \param[in] do_it the user ACTION instance
397 */
398 86 PolygonAction(const ACTION& do_it) :
399 86 do_it_(do_it) {
400 86 }
401
402 /**
403 * \brief Callback called for each polygon.
404 * \details Routes the callback to the wrapped user action class.
405 * \param[in] v index of current Delaunay seed
406 * \param[in] f index of current mesh facet
407 * \param[in] P intersection between current mesh facet
408 * and the Voronoi cell of \p v
409 */
410 779598 void operator() (
411 index_t v,
412 index_t f,
413 const Polygon& P
414 ) const {
415 779598 GEO::geo_argused(f);
416 779598 const_cast<ACTION&> ( do_it_)(v, P);
417 779598 }
418
419 protected:
420 const ACTION& do_it_;
421 };
422
423 /**
424 * \brief Adapter class used internally to implement
425 * for_each_triangle().
426 * \details Overrides constness checks, to allow using temporaries as
427 * argument of for_each_xxx().
428 * \tparam ACTION the user action class
429 */
430 template <class ACTION>
431 class TriangleAction {
432 public:
433 /**
434 * \brief Creates a new TriangleAction that wraps a
435 * user ACTION instance.
436 * \param[in] do_it the user ACTION instance
437 */
438 2360 TriangleAction(const ACTION& do_it) :
439 2360 do_it_(do_it) {
440 2360 }
441
442 /**
443 * \brief Callback called for each integration simplex.
444 * \details Decomposes the polygon \p P into triangles and
445 * calls the callback of the wrapped user action class
446 * for each triangle.
447 * \param[in] v index of current Delaunay seed
448 * \param[in] f index of current mesh facet
449 * \param[in] P intersection between current mesh facet and
450 * the Voronoi cell of \p v
451 */
452 18741602 void operator() (
453 index_t v,
454 index_t f,
455 const Polygon& P
456 ) const {
457 18741602 GEO::geo_argused(f);
458
2/2
✓ Branch 1 taken 18072767 times.
✓ Branch 2 taken 9370801 times.
 54887136 for(index_t i = 1; i + 1 < P.nb_vertices(); i++) {
459 48680692 const_cast<ACTION&> (do_it_)(
460 36145534 v, P.vertex(0), P.vertex(i), P.vertex(i + 1)
461 );
462 }
463 18741602 }
464
465 protected:
466 const ACTION& do_it_;
467 };
468
469 /**
470 * \brief Adapter class used internally
471 * to implement for_each_halfedge().
472 * \details Overrides constness checks, to allow using temporaries as
473 * argument of for_each_xxx().
474 * \tparam ACTION the user action class..
475 */
476 template <class ACTION>
477 class HalfedgeAction {
478 public:
479 /**
480 * \brief Creates a new HalfedgeAction that wraps
481 * a user ACTION instance.
482 * \param[in] do_it the user ACTION instance
483 */
484 HalfedgeAction(const ACTION& do_it) :
485 do_it_(do_it) {
486 }
487
488 /**
489 * \brief Callback called for each integration simplex.
490 * \details Calls the callback of the wrapped
491 * user action class for each edge that has the INTERSECT flag.
492 * \param[in] v index of current Delaunay seed
493 * \param[in] f index of current mesh facet
494 * \param[in] P intersection between current mesh facet
495 * and the Voronoi cell of \p v
496 */
497 void operator() (
498 index_t v,
499 index_t f,
500 const Polygon& P
501 ) const {
502 GEO::geo_argused(v);
503 GEO::geo_argused(f);
504 for(index_t i = 0; i < P.nb_vertices(); i++) {
505 if(P.vertex(i).check_flag(INTERSECT)) {
506 index_t j = P.next_vertex(i);
507 const_cast<ACTION&> (do_it_)(
508 P.vertex(i), P.vertex(j)
509 );
510 }
511 }
512 }
513
514 protected:
515 const ACTION& do_it_;
516 };
517
518 /**
519 * \brief Adapter class used internally to implement
520 * for_each_border_halfedge().
521 * \details Overrides constness checks, to allow using temporaries as
522 * argument of for_each_xxx().
523 * \tparam ACTION the user action class..
524 */
525 template <class ACTION>
526 class BorderHalfedgeAction {
527 public:
528 /**
529 * \brief Creates a new BorderHalfedgeAction that wraps
530 * a user ACTION instance.
531 * \param[in] do_it the user ACTION instance
532 */
533 BorderHalfedgeAction(const ACTION& do_it) :
534 do_it_(do_it) {
535 }
536
537 /**
538 * \brief Callback called for each integration simplex.
539 * \details Calls the callback of the wrapped
540 * user action class for each edge that is on the
541 * border of the input surface.
542 * \param[in] v index of current Delaunay seed
543 * \param[in] f index of current mesh facet
544 * \param[in] P intersection between current mesh facet and
545 * the Voronoi cell of \p v
546 */
547 void operator() (
548 index_t v,
549 index_t f,
550 const Polygon& P
551 ) const {
552 GEO::geo_argused(f);
553 for(index_t i = 0; i < P.nb_vertices(); i++) {
554 if(P.vertex(i).check_flag(ORIGINAL)) {
555 if(P.vertex(i).adjacent_facet() == -1) {
556 index_t j = P.next_vertex(i);
557 const_cast<ACTION&> (do_it_)(
558 v, P.vertex(i), P.vertex(j)
559 );
560 }
561 }
562 }
563 }
564
565 private:
566 const ACTION& do_it_;
567 };
568
569 /**
570 * \brief Adapter class used internally to implement
571 * for_each_primal_triangle()
572 * \details Overrides constness checks, to allow using temporaries as
573 * argument of for_each_xxx()
574 */
575 template <class ACTION>
576 class PrimalTriangleAction {
577 public:
578 /**
579 * \brief Creates a new PrimalTriangleAction that wraps
580 * a user ACTION instance.
581 * \param[in] do_it the user ACTION instance
582 */
583 ✗ PrimalTriangleAction(const ACTION& do_it) :
584 ✗ do_it_(do_it) {
585 ✗ }
586
587 /**
588 * \brief Callback called for each primal triangle.
589 * \param[in] iv1 index of current Delaunay seed
590 * \param[in] f index of current mesh facet
591 * \param[in] P intersection between current mesh facet and
592 * the Voronoi cell of \p v
593 */
594 ✗ void operator() (
595 index_t iv1,
596 index_t f,
597 const Polygon& P
598 ) const {
599 ✗ GEO::geo_argused(f);
600 ✗ for(index_t i = 0; i < P.nb_vertices(); i++) {
601 ✗ const Vertex& ve = P.vertex(i);
602 // Primal triangles correspond to vertices of
603 // the RVD that are on two bisectors.
604 ✗ if(ve.sym().nb_bisectors() == 2) {
605 ✗ index_t iv2 = ve.sym().bisector(0);
606 ✗ index_t iv3 = ve.sym().bisector(1);
607 // This test generates triangle (iv1,iv2,iv3)
608 // only once (i.e. if iv1 is the vertex with
609 // the smallest index).
610 ✗ if(iv1 < iv2 && iv1 < iv3) {
611 ✗ const_cast<ACTION&> (do_it_)(iv1, iv2, iv3);
612 }
613 }
614 }
615 ✗ }
616
617 protected:
618 const ACTION& do_it_;
619 };
620
621 /**
622 * @}
623 * \name Adapter classes for volumetric computation
624 * @{
625 */
626
627 /**
628 * \brief Adapter class used internally to implement
629 * for_each_polyhedron()
630 * \details Overrides constness checks, to allow using temporaries as
631 * argument of for_each_xxx()
632 * \tparam ACTION the user action class
633 */
634 template <class ACTION>
635 class PolyhedronAction {
636 public:
637 /**
638 * \brief Creates a new PolyhedronAction that wraps
639 * a user ACTION instance.
640 * \param[in] do_it the user ACTION instance
641 */
642 8 PolyhedronAction(const ACTION& do_it) :
643 8 do_it_(do_it) {
644 8 }
645
646 /**
647 * \brief Callback called for each polyhedron
648 * \details Routes the callback to the wrapped user action class.
649 * \param[in] v index of current Delaunay seed
650 * \param[in] t index of current mesh tetrahedron
651 * \param[in] C intersection between current mesh tetrahedron
652 * and the Voronoi cell of \p v
653 */
654 39094 void operator() (
655 index_t v,
656 index_t t,
657 const Polyhedron& C
658 ) const {
659 39094 const_cast<ACTION&> ( do_it_)(v, t, C);
660 39094 }
661
662 protected:
663 const ACTION& do_it_;
664 };
665
666 /**
667 * \brief Adapter class used internally to implement
668 * for_each_volumetric_integration_simplex()
669 * \details Overrides constness checks, to allow using temporaries as
670 * argument of for_each_xxx().
671 * \tparam ACTION the user action class. It needs to implement:
672 * operator()(index_t v, signed_index_t v_adj,
673 * index_t t, signed_index_t t_adj,
674 * const Vertex& v1, const Vertex& v2, const Vertex& v3
675 * )
676 * where the parameters are as follows:
677 * - v is the index of the current Voronoi cell
678 * (or Delaunay vertex)
679 * - v_adj is the index of the Voronoi cell adjacent to t accros
680 * facet (\p v1, \p v2, \p v3) or -1 if it does not exists
681 * adjacent to v or -1 if current face is a tetrahedron facet
682 * - t is the index of the current tetrahedron
683 * - t_adj is the index of the tetrahedron adjacent to t accros
684 * facet (\p v1, \p v2, \p v3) or -1 if it does not exists
685 * - v1,v2 and v3 are the three vertices of the facet on the
686 * border of the restricted Voronoi cell.
687 */
688 template <class ACTION>
689 class VolumetricIntegrationSimplexAction {
690 public:
691 /**
692 * \brief Creates a new VolumetricIntegrationSimplexAction
693 * that wraps a user ACTION instance.
694 * \param[in] do_it the user ACTION instance
695 * \param[in] visit_inner_tets if set, all the tetrahedron-cell
696 * intersections are visited, else only tetrahedra on the border
697 * of the restricted Voronoi cell are visited. Since all the
698 * visited triangles are connected to the current Voronoi seed
699 * by a tetrahedron, the computed volume is the same
700 * in both cases.
701 * \param[in] coherent_triangles if set, this ensures that
702 * the polygonal facets of the cells are always triangulated
703 * in a coherent manner when seen from two different cells.
704 * For instance, it is required if a tetrahedral mesh is
705 * reconstructed.
706 */
707 ✗ VolumetricIntegrationSimplexAction(
708 const ACTION& do_it,
709 bool visit_inner_tets = false,
710 bool coherent_triangles = false
711 ) :
712 ✗ do_it_(do_it),
713 ✗ visit_inner_tets_(visit_inner_tets),
714 ✗ coherent_triangles_(coherent_triangles)
715 {
716 ✗ }
717
718 /**
719 * \brief Callback called for each polyhedron
720 * \details Routes the callback to the wrapped user action class.
721 * \param[in] v index of current Delaunay seed
722 * \param[in] t index of current mesh tetrahedron
723 * \param[in] C intersection between current mesh tetrahedron
724 * and the Voronoi cell of \p v
725 */
726 ✗ void operator() (index_t v, index_t t, const Polyhedron& C) const {
727 ✗ for(index_t cv = 0; cv < C.max_v(); ++cv) {
728 ✗ signed_index_t ct = C.vertex_triangle(cv);
729 ✗ if(ct == -1) {
730 ✗ continue;
731 }
732 ✗ geo_debug_assert(C.triangle_is_used(index_t(ct)));
733
734 ✗ signed_index_t adjacent = C.vertex_id(cv);
735 ✗ signed_index_t v_adj = -1;
736 ✗ signed_index_t t_adj = -1;
737
738 ✗ if(adjacent < 0) {
739 // Negative adjacent indices correspond to
740 // tet-tet links (ignored when we want to triangulate
741 // the border of the restricted Voronoi cell while
742 // ignoring internal structures).
743 ✗ if(!visit_inner_tets_) {
744 ✗ continue;
745 }
746 ✗ t_adj = -adjacent - 1;
747 ✗ } else if(adjacent > 0) {
748 // Positive adjacent indices correspond to
749 // Voronoi seed - Voronoi seed link
750 ✗ v_adj = adjacent - 1;
751 }
752 // and adjacent indicex equal to zero corresponds
753 // to tet on border.
754
755 ✗ Polyhedron::Corner c1(
756 index_t(ct),
757 index_t(C.find_triangle_vertex(index_t(ct), cv))
758 );
759
760 // If required, ensure that two polygonal facets
761 // seen from two different volumetric cells will
762 // be triangulated coherently.
763 ✗ if(coherent_triangles_) {
764 ✗ move_to_first_corner_of_facet(C, c1, v);
765 }
766
767 ✗ const Vertex& v1 = C.triangle_dual(c1.t);
768
769 ✗ Polyhedron::Corner c2 = c1;
770 ✗ C.move_to_next_around_vertex(c2);
771 ✗ geo_debug_assert(c2 != c1);
772
773 ✗ Polyhedron::Corner c3 = c2;
774 ✗ C.move_to_next_around_vertex(c3);
775 ✗ geo_debug_assert(c3 != c1);
776 do {
777 ✗ const Vertex& v2 = C.triangle_dual(c2.t);
778 ✗ const Vertex& v3 = C.triangle_dual(c3.t);
779 ✗ const_cast<ACTION&> (do_it_)(
780 v, index_t(v_adj), t, index_t(t_adj), v1, v2, v3
781 );
782 ✗ c2 = c3;
783 ✗ C.move_to_next_around_vertex(c3);
784 ✗ } while(c3 != c1);
785 }
786 ✗ }
787
788 /**
789 * \brief Finds the first corner of a facet in a Polyhedron.
790 * \details This function is used to ensure that a facet is
791 * triangulated coherently when seen from two different
792 * volumetric cells, by generating a fan of triangles
793 * that radiates from the first corner. The global order
794 * used to find the first
795 * corner is defined by the function symbolic_compare().
796 *
797 * \param[in] C the Polyhedron
798 * \param[in,out] c a corner of the facet, replaced by the
799 * first corner of the facet on exit.
800 * \param[in] center_vertex_id index of the current Voronoi seed
801 * (needed to determine the full symbolic information in the
802 * vertices).
803 */
804 ✗ void move_to_first_corner_of_facet(
805 const Polyhedron& C, Polyhedron::Corner& c,
806 index_t center_vertex_id
807 ) const {
808 ✗ Polyhedron::Corner first = c;
809 ✗ Polyhedron::Corner cur = c;
810 do {
811 ✗ if(symbolic_compare(
812 ✗ C.triangle_dual(cur.t),
813 ✗ C.triangle_dual(c.t),
814 center_vertex_id
815 )) {
816 ✗ c = cur;
817 }
818 ✗ C.move_to_next_around_vertex(cur);
819 ✗ } while(cur != first);
820 ✗ }
821
822 /**
823 * \brief Compares the symbolic information of two vertices
824 * in such a way that a global order is defined.
825 * \details This function is used to ensure that a facet is
826 * triangulated coherently when seen from two different
827 * volumetric cells (it uniquely determines the "first" vertex).
828 * \param[in] p1 first vertex to compare
829 * \param[in] p2 second vertex to compare
830 * \param[in] center_vertex_id index of the current Voronoi seed
831 * (needed to determine the full symbolic information in
832 * \p p1 and \p p2).
833 * \return true if p1 is before p2 in the global order,
834 * false otherwise.
835 */
836 ✗ static bool symbolic_compare(
837 const Vertex& p1, const Vertex& p2, index_t center_vertex_id
838 ) {
839 ✗ GEO::signed_quadindex K1(
840 signed_index_t(center_vertex_id),
841 ✗ p1.sym()[0], p1.sym()[1], p1.sym()[2]
842 );
843 ✗ GEO::signed_quadindex K2(
844 signed_index_t(center_vertex_id),
845 ✗ p2.sym()[0], p2.sym()[1], p2.sym()[2]
846 );
847 ✗ return K1 < K2;
848 }
849
850 protected:
851 const ACTION& do_it_;
852 bool visit_inner_tets_;
853 bool coherent_triangles_;
854 };
855
856 /**
857 * \brief Adapter class used internally to implement
858 * for_each_tetrahedron()
859 * \details Overrides constness checks, to allow using temporaries as
860 * argument of for_each_xxx()
861 * \tparam ACTION the user action class. It needs to implement:
862 * operator()(index_t v, signed_index_t v_adj,
863 * index_t t, index_t t_adj,
864 * const Vertex& v0, const Vertex& v1,
865 * const Vertex& v2, const Vertex& v3
866 * )
867 * where the parameters are as follows:
868 * - v is the index of the current Voronoi cell
869 * (or Delaunay vertex)
870 * - v_adj is the index of the Voronoi cell adjacent to t accros
871 * facet (\p v1, \p v2, \p v3) or -1 if it does not exists
872 * adjacent to v or -1 if current face is a tetrahedron facet
873 * - t is the index of the current tetrahedron
874 * - t_adj is the index of the tetrahedron adjacent to t accros
875 * facet (\p v1, \p v2, \p v3) or -1 if it does not exists
876 * - v0,v1,v2 and v3 are the four vertices of tetrahedron.
877 */
878 template <class ACTION>
879 class TetrahedronAction {
880 public:
881 /**
882 * \brief Creates a new TetrahedronAction that wraps
883 * a user ACTION instance.
884 * \param[in] do_it the user ACTION instance
885 */
886 1280 TetrahedronAction(
887 const ACTION& do_it
888 ) :
889 1280 do_it_(do_it)
890 {
891 1280 }
892
893 /**
894 * \brief Callback called for each polyhedron
895 * \details Routes the callback to the wrapped user action class.
896 * \param[in] v index of current Delaunay seed
897 * \param[in] t index of current mesh tetrahedron
898 * \param[in] C intersection between current mesh tetrahedron
899 * and the Voronoi cell of \p v
900 */
901 1579756 void operator() (
902 index_t v,
903 index_t t,
904 const Polyhedron& C
905 ) const {
906
907 // Find a vertex of the current cell,
908 // that will be used as the 'origin'
909 // vertex
910 1579756 const Vertex* v0 = nullptr;
911 index_t t0;
912
2/2
✓ Branch 1 taken 1189307 times.
✓ Branch 2 taken 1793 times.
 2382200 for(t0 = 0; t0 < C.max_t(); ++t0) {
913
2/2
✓ Branch 1 taken 788085 times.
✓ Branch 2 taken 401222 times.
 2378614 if(C.triangle_is_used(t0)) {
914 1576170 v0 = &C.triangle_dual(t0);
915 1576170 break;
916 }
917 }
918
919 // If current cell is empty, return
920
2/2
✓ Branch 0 taken 1793 times.
✓ Branch 1 taken 788085 times.
 1579756 if(v0 == nullptr) {
921 3586 return;
922 }
923
924
2/2
✓ Branch 1 taken 12923894 times.
✓ Branch 2 taken 788085 times.
 27423958 for(index_t cv = 0; cv < C.max_v(); ++cv) {
925
1/2
✓ Branch 1 taken 12923894 times.
✗ Branch 2 not taken.
 25847788 signed_index_t ct = C.vertex_triangle(cv);
926
2/2
✓ Branch 0 taken 7828173 times.
✓ Branch 1 taken 5095721 times.
 25847788 if(ct == -1) {
927 20384856 continue;
928 }
929
2/8
✓ Branch 1 taken 5095721 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 5095721 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 10191442 geo_debug_assert(C.triangle_is_used(index_t(ct)));
930
931
1/2
✓ Branch 1 taken 5095721 times.
✗ Branch 2 not taken.
 10191442 signed_index_t adjacent = C.vertex_id(cv);
932 10191442 signed_index_t v_adj = -1;
933 10191442 signed_index_t t_adj = -1;
934
935
2/2
✓ Branch 0 taken 2017011 times.
✓ Branch 1 taken 3078710 times.
 10191442 if(adjacent < 0) {
936 // Negative adjacent indices correspond to
937 // tet-tet links
938 4034022 t_adj = -adjacent - 1;
939
2/2
✓ Branch 0 taken 2883354 times.
✓ Branch 1 taken 195356 times.
 6157420 } else if(adjacent > 0) {
940 // Positive adjacent indices correspond to
941 // Voronoi seed - Voroni seed link
942 5766708 v_adj = adjacent - 1;
943 }
944 // and adjacent indicex equal to zero corresponds
945 // to tet on border.
946
947
1/2
✓ Branch 1 taken 5095721 times.
✗ Branch 2 not taken.
 10191442 Polyhedron::Corner c1(
948 index_t(ct), C.find_triangle_vertex(index_t(ct), cv)
949 );
950
951 // If the current facet is incident to
952 // the origin vertex, then skip it (else
953 // it would generate flat tetrahedra)
954
3/4
✓ Branch 1 taken 5095721 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2364255 times.
✓ Branch 4 taken 2731466 times.
 10191442 if(facet_is_incident_to_vertex(C, c1, t0)) {
955 4728510 continue;
956 }
957
958
1/2
✓ Branch 1 taken 2731466 times.
✗ Branch 2 not taken.
 5462932 const Vertex& v1 = C.triangle_dual(c1.t);
959
960 5462932 Polyhedron::Corner c2 = c1;
961
1/2
✓ Branch 1 taken 2731466 times.
✗ Branch 2 not taken.
 5462932 C.move_to_next_around_vertex(c2);
962
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 2731466 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 5462932 geo_debug_assert(c2 != c1);
963
964 5462932 Polyhedron::Corner c3 = c2;
965
1/2
✓ Branch 1 taken 2731466 times.
✗ Branch 2 not taken.
 5462932 C.move_to_next_around_vertex(c3);
966
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 2731466 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 5462932 geo_debug_assert(c3 != c1);
967 do {
968
1/2
✓ Branch 1 taken 5833391 times.
✗ Branch 2 not taken.
 11666782 const Vertex& v2 = C.triangle_dual(c2.t);
969
1/2
✓ Branch 1 taken 5833391 times.
✗ Branch 2 not taken.
 11666782 const Vertex& v3 = C.triangle_dual(c3.t);
970
1/4
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 5 taken 5833391 times.
✗ Branch 6 not taken.
 11666782 const_cast<ACTION&> (do_it_)(
971 v, index_t(v_adj), t, index_t(t_adj),
972 *v0, v1, v2, v3
973 );
974 11666782 c2 = c3;
975
1/2
✓ Branch 1 taken 5833391 times.
✗ Branch 2 not taken.
 11666782 C.move_to_next_around_vertex(c3);
976
2/2
✓ Branch 1 taken 3101925 times.
✓ Branch 2 taken 2731466 times.
 11666782 } while(c3 != c1);
977 }
978 }
979
980 protected:
981 /**
982 * \brief Tests whether a Polyhedron facet is incident
983 * to a vertex.
984 * \param[in] C the Polyhedron
985 * \param[in] c a corner of the facet
986 * \param[in] t the index of the vertex in dual form (in other
987 * words, a triangle index).
988 * \return true if the facet incident to corner \p c
989 * is also incident to the vertex dual to \p t, false otherwise
990 */
991 10191442 bool facet_is_incident_to_vertex(
992 const Polyhedron& C, Polyhedron::Corner& c, index_t t
993 ) const {
994 10191442 Polyhedron::Corner first = c;
995 10191442 Polyhedron::Corner cur = c;
996 do {
997
2/2
✓ Branch 0 taken 2364255 times.
✓ Branch 1 taken 16220446 times.
 37169402 if(cur.t == t) {
998 4728510 return true;
999 }
1000
1/2
✓ Branch 1 taken 16220446 times.
✗ Branch 2 not taken.
 32440892 C.move_to_next_around_vertex(cur);
1001
2/2
✓ Branch 1 taken 13488980 times.
✓ Branch 2 taken 2731466 times.
 32440892 } while(cur != first);
1002 5462932 return false;
1003 }
1004
1005 protected:
1006 const ACTION& do_it_;
1007 };
1008
1009 /**
1010 * \brief Adapter class used internally to implement
1011 * for_each_primal_tetrahedron()
1012 * \details Overrides constness checks, to allow using temporaries as
1013 * argument of for_each_xxx()
1014 * \tparam ACTION the user action class
1015 */
1016 template <class ACTION>
1017 class PrimalTetrahedronAction {
1018 public:
1019 /**
1020 * \brief Constructs a new PrimalTetrahedronAction.
1021 * \param[in] do_it the user ACTION instance.
1022 */
1023 ✗ PrimalTetrahedronAction(const ACTION& do_it) :
1024 ✗ do_it_(do_it) {
1025 ✗ }
1026
1027 /**
1028 * \brief Callback called for each polyhedron
1029 * \details Routes the callback to the wrapped user action class.
1030 * \param[in] v index of current Delaunay seed
1031 * \param[in] t index of current mesh tetrahedron
1032 * \param[in] C intersection between current mesh tetrahedron
1033 * and the Voronoi cell of \p v
1034 */
1035 ✗ void operator() (
1036 index_t v,
1037 index_t t,
1038 const Polyhedron& C
1039 ) const {
1040 ✗ GEO::geo_argused(t);
1041 ✗ for(index_t it = 0; it < C.max_t(); ++it) {
1042 ✗ if(C.triangle_is_used(it)) {
1043 ✗ const SymbolicVertex& sym = C.triangle_dual(it).sym();
1044 ✗ if(sym.nb_bisectors() == 3) {
1045 ✗ index_t v1 = sym.bisector(0);
1046 ✗ index_t v2 = sym.bisector(1);
1047 ✗ index_t v3 = sym.bisector(2);
1048 // This test ensures that the tet (v,v1,v2,v3)
1049 // is generated only once.
1050 ✗ if(v < v1 && v < v2 && v < v3) {
1051 ✗ const_cast<ACTION&> (do_it_)(v, v1, v2, v3);
1052 }
1053 }
1054 }
1055 }
1056 ✗ }
1057
1058 protected:
1059 const ACTION& do_it_;
1060 };
1061
1062 public:
1063 /**
1064 * @}
1065 * \name Public interface for computation/iteration
1066 * @{
1067 */
1068
1069 /**
1070 * \brief Iterates on the facets of this RVD.
1071 * \param[in] action the user action object
1072 * \tparam ACTION needs to implement:
1073 * operator()(index_t v, index_t f, const Polygon& P) const
1074 * where v denotes the index of the current Voronoi cell
1075 * (or Delaunay vertex), f the index of the current facet
1076 * and P the computed intersection between facet f
1077 * and the Voronoi cell of v.
1078 */
1079 template <class ACTION>
1080 86 inline void for_each_polygon(const ACTION& action) {
1081
1/2
✓ Branch 1 taken 43 times.
✗ Branch 2 not taken.
 86 this->template compute_surfacic<PolygonAction<ACTION> >(
1082 86 PolygonAction<ACTION>(action)
1083 );
1084 86 }
1085
1086 /**
1087 * \brief Iterates on the facets of this RVD, triangulated on the fly.
1088 * \param[in] action the user action object
1089 * \tparam TRIACTION needs to implement:
1090 * operator()(index_t c, const TopoPolyVertex& v1, v2, v3) const
1091 * where c denotes the index of the current Voronoi cell
1092 * (or Delaunay vertex).
1093 */
1094 template <class TRIACTION>
1095 2360 inline void for_each_triangle(const TRIACTION& action) {
1096
1/2
✓ Branch 1 taken 1180 times.
✗ Branch 2 not taken.
 2360 this->template compute_surfacic<TriangleAction<TRIACTION> >(
1097 2360 TriangleAction<TRIACTION>(action)
1098 );
1099 2360 }
1100
1101 /**
1102 * \brief Iterates on the halfedges on the borders of the
1103 * restricted Voronoi cells.
1104 * \param[in] action the user action object
1105 * \tparam HEACTION needs to implement:
1106 * operator()(index_t c, const TopoPolyVertex& v1, v2) const
1107 * where c denotes the index of the current Voronoi cell
1108 * (or Delaunay vertex).
1109 */
1110 template <class HEACTION>
1111 inline void for_each_halfedge(const HEACTION& action) {
1112 this->template compute_surfacic<HalfedgeAction<HEACTION> >(
1113 HalfedgeAction<HEACTION>(action)
1114 );
1115 }
1116
1117 /**
1118 * \brief Iterates on the halfedges on the borders of the
1119 * restricted Voronoi cells that are on the boundary of the input mesh.
1120 * \param[in] action the user action object
1121 * \tparam BOACTION needs to implement:
1122 * operator()(index_t c, const TopoPolyVertexEdge& v1, v2) const
1123 * where c denotes the index of the current Voronoi cell
1124 * (or Delaunay vertex).
1125 */
1126 template <class BOACTION>
1127 inline void for_each_border_halfedge(const BOACTION& action) {
1128 this->template compute_surfacic<BorderHalfedgeAction<BOACTION> >(
1129 BorderHalfedgeAction<BOACTION>(action)
1130 );
1131 }
1132
1133 /**
1134 * \brief Iterates on the triangles of the Restricted
1135 * Delaunay Triangulation.
1136 * \param[in] action the user action object
1137 * \tparam PRIMTRIACTION needs to implement:
1138 * operator()(index_t i, unsigned j, index_t k) const
1139 * where i,j,k denote the three indices of the Delaunay vertices
1140 * that define the primal triangle.
1141 */
1142 template <class PRIMTRIACTION>
1143 ✗ inline void for_each_primal_triangle(const PRIMTRIACTION& action) {
1144 ✗ bool sym_backup = symbolic();
1145 ✗ set_symbolic(true);
1146 ✗ this->template compute_surfacic<PrimalTriangleAction<PRIMTRIACTION>>(
1147 ✗ PrimalTriangleAction<PRIMTRIACTION>(action)
1148 );
1149 ✗ set_symbolic(sym_backup);
1150 ✗ }
1151
1152 /**
1153 * \brief Iterates on the polyhedra of this RVD.
1154 * \param[in] action the user action object
1155 * \tparam ACTION needs to implement:
1156 * operator()(index_t v, index_t t, const Polyhedron& C) const
1157 * where v denotes the index of the current Voronoi cell
1158 * (or Delaunay vertex), t the index of the current tetrahedron
1159 * and C the computed intersection between tetrahedron t
1160 * and the Voronoi cell of v.
1161 */
1162 template <class ACTION>
1163 8 inline void for_each_polyhedron(const ACTION& action) {
1164
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 this->template compute_volumetric<PolyhedronAction<ACTION> >(
1165 8 PolyhedronAction<ACTION>(action)
1166 );
1167 8 }
1168
1169 /**
1170 * \brief Iterates on the polyhedra of this RVD decomposed
1171 * on the fly into tetrahedra.
1172 * \details The generated tetrahedra may be geometrically incorrect,
1173 * but they are algebraically correct. In other word, their signed
1174 * volumes sum as the volume of the restricted Voronoi cell.
1175 * \param[in] action the user action object
1176 * \param[in] visit_inner_tets if set, all the tetrahedron-cell
1177 * intersections are visited, else only tetrahedra on the border
1178 * of the restricted Voronoi cell are visited. Since all the visited
1179 * triangles are connected to the current Voronoi seed by a
1180 * tetrahedron, the computed volume is the same in both cases.
1181 * \param[in] coherent_triangles if set, this ensures that the
1182 * polygonal facets of the cells are always triangulated in a
1183 * coherent manner when seen from two different cells.
1184 * For instance, it is required if a tetrahedral mesh is
1185 * reconstructed.
1186 * \tparam ACTION needs to implement:
1187 * operator()(index_t v, signed_index_t v_adj,
1188 * index_t t, index_t t_adj,
1189 * const Vertex& v1, const Vertex& v2, const Vertex& v3
1190 * )
1191 * where the parameters are as follows:
1192 * - v is the index of the current Voronoi cell
1193 * (or Delaunay vertex)
1194 * - v_adj is the index of the Voronoi cell adjacent to t accros
1195 * facet (\p v1, \p v2, \p v3) or -1 if it does not exists
1196 * adjacent to v or -1 if current face is a tetrahedron facet
1197 * - t is the index of the current tetrahedron
1198 * - t_adj is the index of the tetrahedron adjacent to t accros
1199 * facet (\p v1, \p v2, \p v3) or -1 if it does not exists
1200 * - v1,v2 and v3 are the three vertices of the facet on the
1201 * border of the restricted Voronoi cell.
1202 */
1203 template <class ACTION>
1204 ✗ inline void for_each_volumetric_integration_simplex(
1205 const ACTION& action,
1206 bool visit_inner_tets = false, bool coherent_triangles = false
1207 ) {
1208 this->template compute_volumetric<
1209 VolumetricIntegrationSimplexAction<ACTION>
1210 ✗ >(
1211 ✗ VolumetricIntegrationSimplexAction<ACTION>(
1212 action, visit_inner_tets, coherent_triangles
1213 )
1214 );
1215 ✗ }
1216
1217 /**
1218 * \brief Iterates on the polyhedra of this RVD decomposed
1219 * on the fly into tetrahedra.
1220 * \details The tetrahedra are generated by connecting one of
1221 * the vertices of the cell to the other ones.
1222 * \param[in] action the user action object
1223 * \tparam ACTION needs to implement:
1224 * operator()(index_t v, signed_index_t v_adj,
1225 * index_t t, index_t t_adj,
1226 * const Vertex& v0, const Vertex& v1,
1227 * const Vertex& v2, const Vertex& v3
1228 * )
1229 * where the parameters are as follows:
1230 * - v is the index of the current Voronoi cell
1231 * (or Delaunay vertex)
1232 * - v_adj is the index of the Voronoi cell adjacent to t accros
1233 * facet (\p v1, \p v2, \p v3) or -1 if it does not exists
1234 * adjacent to v or -1 if current face is a tetrahedron facet
1235 * - t is the index of the current tetrahedron
1236 * - t_adj is the index of the tetrahedron adjacent to t accros
1237 * facet (\p v1, \p v2, \p v3) or -1 if it does not exists
1238 * - v0,v1,v2 and v3 are the four vertices of tetrahedron.
1239 */
1240 template <class ACTION>
1241 1280 inline void for_each_tetrahedron(
1242 const ACTION& action
1243 ) {
1244
1/2
✓ Branch 1 taken 640 times.
✗ Branch 2 not taken.
 1280 this->template compute_volumetric<TetrahedronAction<ACTION> >(
1245 1280 TetrahedronAction<ACTION>(
1246 action
1247 )
1248 );
1249 1280 }
1250
1251 /**
1252 * \brief Iterates on the primal tetrahedra of this RVD.
1253 * \details The tetrahedra are not coherently oriented,
1254 * and need a subsequent traversal operation to reorient
1255 * them. They can be also reoriented geometrically using
1256 * the orient3d() predicate.
1257 * \param[in] action the user action object
1258 * \tparam ACTION needs to implement:
1259 * operator()(index_t v0, index_t v1, index_t v2, index_t v3)
1260 * where v0,v1,v2 and v3 are the indices of the four vertices
1261 * of tetrahedron.
1262 */
1263 template <class ACTION>
1264 ✗ inline void for_each_primal_tetrahedron(const ACTION& action) {
1265 ✗ bool sym_backup = symbolic();
1266 ✗ set_symbolic(true);
1267 ✗ this->template compute_volumetric<PrimalTetrahedronAction<ACTION> >(
1268 ✗ PrimalTetrahedronAction<ACTION>(
1269 action
1270 )
1271 );
1272 ✗ set_symbolic(sym_backup);
1273 ✗ }
1274
1275 protected:
1276 /**
1277 * @}
1278 * \name Computation
1279 * @{
1280 */
1281
1282 /**
1283 * \brief Low-level API of Restricted Voronoi Diagram traversal.
1284 * \details Client code may use for_each_facet(),for_each_triangle() or
1285 * for_each_primal_triangle() instead.
1286 * \tparam ACTION needs to implement:
1287 * operator()(index_t v, index_t f, const Polygon& P) const
1288 * where v denotes the index of the current Voronoi cell
1289 * (or Delaunay vertex), f the index of the current facet
1290 * and P the computed intersection between the Voronoi cell of
1291 * v and facet f.
1292 */
1293 template <class ACTION>
1294 2446 inline void compute_surfacic(const ACTION& action) {
1295
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 1216 times.
 2446 if(connected_components_priority_) {
1296 14 this->template compute_surfacic_with_cnx_priority<ACTION>(
1297 action
1298 );
1299 } else {
1300 2432 this->template compute_surfacic_with_seeds_priority<ACTION>(
1301 action
1302 );
1303 }
1304 2446 }
1305
1306 /**
1307 * \brief Low-level API of Restricted Voronoi Diagram traversal
1308 * with seeds priority in surfacic mode.
1309 * \details Client code may use for_each_facet(),for_each_triangle() or
1310 * for_each_primal_triangle() instead.
1311 * \tparam ACTION needs to implement:
1312 * operator()(index_t v, index_t f, const Polygon& P) const
1313 * where v denotes the index of the current Voronoi cell
1314 * (or Delaunay vertex), f the index of the current facet
1315 * and P the computed intersection between the Voronoi cell of
1316 * v and facet f.
1317 */
1318 template <class ACTION>
1319 2432 inline void compute_surfacic_with_seeds_priority(const ACTION& action) {
1320 2432 if(
1321
2/2
✓ Branch 0 taken 36 times.
✓ Branch 1 taken 1180 times.
 2432 facets_begin_ == UNSPECIFIED_RANGE &&
1322
1/2
✓ Branch 0 taken 36 times.
✗ Branch 1 not taken.
 72 facets_end_ == UNSPECIFIED_RANGE
1323 ) {
1324 72 facets_begin_ = 0;
1325 72 facets_end_ = mesh_->facets.nb();
1326 }
1327 2432 current_polygon_ = nullptr;
1328
1/2
✓ Branch 1 taken 1216 times.
✗ Branch 2 not taken.
 2432 GEO::vector<index_t> seed_stamp(
1329 2432 delaunay_->nb_vertices(), index_t(-1)
1330 );
1331
1/2
✓ Branch 1 taken 1216 times.
✗ Branch 2 not taken.
 2432 GEO::vector<bool> facet_is_marked(facets_end_-facets_begin_, false);
1332
1/2
✓ Branch 1 taken 1216 times.
✗ Branch 2 not taken.
 2432 init_get_neighbors();
1333
1334
1/2
✓ Branch 1 taken 1216 times.
✗ Branch 2 not taken.
 2432 FacetSeedStack adjacent_facets;
1335
1/2
✓ Branch 1 taken 1216 times.
✗ Branch 2 not taken.
 2432 SeedStack adjacent_seeds;
1336
1/2
✓ Branch 1 taken 1216 times.
✗ Branch 2 not taken.
 2432 Polygon F;
1337
1/2
✓ Branch 1 taken 1216 times.
✗ Branch 2 not taken.
 2432 GEO::Attribute<double> vertex_weight;
1338
1/2
✓ Branch 1 taken 1216 times.
✗ Branch 2 not taken.
 2432 vertex_weight.bind_if_is_defined(
1339
1/2
✓ Branch 1 taken 1216 times.
✗ Branch 2 not taken.
 4864 mesh_->vertices.attributes(), "weight"
1340 );
1341
1342 // The algorithm propagates along both the facet-graph of
1343 // the surface and the 1-skeleton of the Delaunay triangulation,
1344 // and computes all the relevant intersections between
1345 // each Voronoi cell and facet.
1346
2/2
✓ Branch 0 taken 3808107 times.
✓ Branch 1 taken 1216 times.
 7618646 for(index_t f = facets_begin_; f < facets_end_; f++) {
1347
2/2
✓ Branch 2 taken 1567 times.
✓ Branch 3 taken 3806540 times.
 7616214 if(!facet_is_marked[f-facets_begin_]) {
1348 // Propagate along the facet-graph.
1349 3134 facet_is_marked[f-facets_begin_] = true;
1350
1/2
✓ Branch 1 taken 1567 times.
✗ Branch 2 not taken.
 3134 adjacent_facets.push(
1351
1/2
✓ Branch 1 taken 1567 times.
✗ Branch 2 not taken.
 3134 FacetSeed(f, find_seed_near_facet(f))
1352 );
1353
2/2
✓ Branch 1 taken 3808107 times.
✓ Branch 2 taken 1567 times.
 7619348 while(!adjacent_facets.empty()) {
1354 7616214 current_facet_ = adjacent_facets.top().f;
1355 7616214 current_seed_ = adjacent_facets.top().seed;
1356 7616214 adjacent_facets.pop();
1357
1358 // Copy the current facet from the Mesh into
1359 // RestrictedVoronoiDiagram's Polygon data structure
1360 // (gathers all the necessary information)
1361 7616214 F.initialize_from_mesh_facet(
1362
1/2
✓ Branch 1 taken 3808107 times.
✗ Branch 2 not taken.
 7616214 mesh_, current_facet_, symbolic_, vertex_weight
1363 );
1364
1365 // Propagate along the Delaunay 1-skeleton
1366 // This will traverse all the seeds such that their
1367 // Voronoi cell has a non-empty intersection with
1368 // the current facet.
1369
1/2
✓ Branch 1 taken 3808107 times.
✗ Branch 2 not taken.
 7616214 seed_stamp[current_seed_] = current_facet_;
1370
1/2
✓ Branch 1 taken 3808107 times.
✗ Branch 2 not taken.
 7616214 adjacent_seeds.push(current_seed_);
1371
1372
2/2
✓ Branch 1 taken 9581210 times.
✓ Branch 2 taken 3808107 times.
 26778634 while(!adjacent_seeds.empty()) {
1373 19162420 current_seed_ = adjacent_seeds.top();
1374 19162420 adjacent_seeds.pop();
1375
1376
1/2
✓ Branch 1 taken 9581210 times.
✗ Branch 2 not taken.
 19162420 current_polygon_ = intersect_cell_facet(
1377 current_seed_, F
1378 );
1379
1380
1/2
✓ Branch 2 taken 9581210 times.
✗ Branch 3 not taken.
 19162420 action(
1381 current_seed_, current_facet_, current_polygon()
1382 );
1383
1384 // Propagate to adjacent facets and adjacent seeds
1385 94465676 for(index_t v = 0;
1386
3/4
✓ Branch 2 taken 47232838 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 37651628 times.
✓ Branch 5 taken 9581210 times.
 94465676 v < current_polygon().nb_vertices(); v++
1387 ) {
1388
1/2
✓ Branch 2 taken 37651628 times.
✗ Branch 3 not taken.
 75303256 const Vertex& ve = current_polygon().vertex(v);
1389 75303256 signed_index_t neigh_f = ve.adjacent_facet();
1390 75303256 if(
1391
2/2
✓ Branch 0 taken 19499097 times.
✓ Branch 1 taken 18152531 times.
 75303256 neigh_f >= signed_index_t(facets_begin_) &&
1392
2/2
✓ Branch 0 taken 19261869 times.
✓ Branch 1 taken 237228 times.
 38998194 neigh_f < signed_index_t(facets_end_) &&
1393
1/2
✓ Branch 0 taken 19261869 times.
✗ Branch 1 not taken.
 38523738 neigh_f != signed_index_t(current_facet_)
1394 ) {
1395
2/2
✓ Branch 1 taken 3806540 times.
✓ Branch 2 taken 15455329 times.
 38523738 if(!facet_is_marked[
1396 38523738 index_t(neigh_f)-facets_begin_
1397 ]) {
1398 facet_is_marked[
1399 7613080 index_t(neigh_f)-facets_begin_
1400 7613080 ] = true;
1401
1/2
✓ Branch 1 taken 3806540 times.
✗ Branch 2 not taken.
 7613080 adjacent_facets.push(
1402 15226160 FacetSeed(
1403 index_t(neigh_f),
1404 current_seed_
1405 )
1406 );
1407 }
1408 }
1409 75303256 signed_index_t neigh_s = ve.adjacent_seed();
1410
2/2
✓ Branch 0 taken 17830376 times.
✓ Branch 1 taken 19821252 times.
 75303256 if(neigh_s != -1) {
1411 35660752 if(
1412
3/4
✓ Branch 1 taken 17830376 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5773103 times.
✓ Branch 4 taken 12057273 times.
 35660752 seed_stamp[neigh_s] != current_facet_
1413 ) {
1414
1/2
✓ Branch 1 taken 5773103 times.
✗ Branch 2 not taken.
 11546206 seed_stamp[neigh_s] = current_facet_;
1415
1/2
✓ Branch 1 taken 5773103 times.
✗ Branch 2 not taken.
 11546206 adjacent_seeds.push(index_t(neigh_s));
1416 }
1417 }
1418 }
1419 }
1420 }
1421 }
1422 }
1423 2432 current_polygon_ = nullptr;
1424 2432 }
1425
1426 /**
1427 * \brief Low-level API of Restricted Voronoi Diagram traversal .
1428 * \details Selects seed-priority or tetrahedron-priority modes
1429 * according to connected_components_priority mode.
1430 * Client code may use for_each_polyhedron() or
1431 * for_each_volumetric_integration_simplex() instead of this function.
1432 * \tparam ACTION needs to implement:
1433 * operator()(index_t v, index_t t, const Polyhedron& C) const
1434 * where v denotes the index of the current Voronoi cell
1435 * (or Delaunay vertex), t the index of the current tetrahedron
1436 * and C the computed intersection between the Voronoi cell of
1437 * v and tetrahedron t
1438 */
1439 template <class ACTION>
1440 1288 inline void compute_volumetric(const ACTION& action) {
1441
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 640 times.
 1288 if(connected_components_priority_) {
1442 8 this->template compute_volumetric_with_cnx_priority<ACTION>(
1443 action
1444 );
1445 } else {
1446 1280 this->template compute_volumetric_with_seeds_priority<ACTION>(
1447 action
1448 );
1449 }
1450 1288 }
1451
1452 /**
1453 * \brief Low-level API of Restricted Voronoi Diagram traversal
1454 * with seeds priority in volumetric mode.
1455 * \details Client code may use for_each_polyhedron() or
1456 * for_each_volumetric_integration_simplex() instead.
1457 * \tparam ACTION needs to implement:
1458 * operator()(index_t v, index_t t, const Polyhedron& C) const
1459 * where v denotes the index of the current Voronoi cell
1460 * (or Delaunay vertex), t the index of the current tetrahedron
1461 * and C the computed intersection between the Voronoi cell of
1462 * v and tetrahedron t
1463 */
1464 template <class ACTION>
1465 1280 inline void compute_volumetric_with_seeds_priority(const ACTION& action){
1466 1280 if(
1467
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 640 times.
 1280 tets_begin_ == UNSPECIFIED_RANGE &&
1468 ✗ tets_end_ == UNSPECIFIED_RANGE
1469 ) {
1470 ✗ tets_begin_ = 0;
1471 ✗ tets_end_ = mesh_->cells.nb();
1472 }
1473
1474
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 640 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 1280 geo_assert(tets_begin_ != UNSPECIFIED_RANGE);
1475
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 640 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 1280 geo_assert(tets_end_ != UNSPECIFIED_RANGE);
1476
1477
1/2
✓ Branch 1 taken 640 times.
✗ Branch 2 not taken.
 1280 GEO::vector<index_t> seed_stamp(
1478 1280 delaunay_->nb_vertices(), index_t(-1)
1479 );
1480
1/2
✓ Branch 1 taken 640 times.
✗ Branch 2 not taken.
 1280 GEO::vector<bool> tet_is_marked(tets_end_-tets_begin_, false);
1481
1/2
✓ Branch 1 taken 640 times.
✗ Branch 2 not taken.
 1280 init_get_neighbors();
1482
1483
1/2
✓ Branch 1 taken 640 times.
✗ Branch 2 not taken.
 1280 TetSeedStack adjacent_tets;
1484
1/2
✓ Branch 1 taken 640 times.
✗ Branch 2 not taken.
 1280 SeedStack adjacent_seeds;
1485
1/2
✓ Branch 2 taken 640 times.
✗ Branch 3 not taken.
 1280 Polyhedron C(dimension());
1486
1/2
✓ Branch 1 taken 640 times.
✗ Branch 2 not taken.
 1280 GEO::Attribute<double> vertex_weight;
1487
1/2
✓ Branch 1 taken 640 times.
✗ Branch 2 not taken.
 1280 vertex_weight.bind_if_is_defined(
1488
1/2
✓ Branch 1 taken 640 times.
✗ Branch 2 not taken.
 2560 mesh_->vertices.attributes(), "weight"
1489 );
1490
1491 1280 current_polyhedron_ = &C;
1492 // The algorithm propagates along both the facet-graph of
1493 // the surface and the 1-skeleton of the Delaunay triangulation,
1494 // and computes all the relevant intersections between
1495 // each Voronoi cell and facet.
1496
2/2
✓ Branch 0 taken 61440 times.
✓ Branch 1 taken 640 times.
 124160 for(index_t t = tets_begin_; t < tets_end_; ++t) {
1497
2/2
✓ Branch 2 taken 640 times.
✓ Branch 3 taken 60800 times.
 122880 if(!tet_is_marked[t-tets_begin_]) {
1498 // Propagate along the tet-graph.
1499 1280 tet_is_marked[t-tets_begin_] = true;
1500
1/2
✓ Branch 1 taken 640 times.
✗ Branch 2 not taken.
 1280 adjacent_tets.push(
1501
1/2
✓ Branch 1 taken 640 times.
✗ Branch 2 not taken.
 1280 TetSeed(t, find_seed_near_tet(t))
1502 );
1503
2/2
✓ Branch 1 taken 61440 times.
✓ Branch 2 taken 640 times.
 124160 while(!adjacent_tets.empty()) {
1504 122880 current_tet_ = adjacent_tets.top().f;
1505 122880 current_seed_ = adjacent_tets.top().seed;
1506 122880 adjacent_tets.pop();
1507
1508 // Note: current cell could be looked up here,
1509 // (from current_tet_) if we chose to keep it
1510 // and copy it right before clipping (I am
1511 // not sure that it is worth it, lookup time
1512 // will be probably fast enough)
1513
1514 // Propagate along the Delaunay 1-skeleton
1515 // This will traverse all the seeds such that their
1516 // Voronoi cell has a non-empty intersection with
1517 // the current facet.
1518
1/2
✓ Branch 1 taken 61440 times.
✗ Branch 2 not taken.
 122880 seed_stamp[current_seed_] = current_tet_;
1519
1/2
✓ Branch 1 taken 61440 times.
✗ Branch 2 not taken.
 122880 adjacent_seeds.push(current_seed_);
1520
1521
2/2
✓ Branch 1 taken 789878 times.
✓ Branch 2 taken 61440 times.
 1702636 while(!adjacent_seeds.empty()) {
1522 1579756 current_seed_ = adjacent_seeds.top();
1523 1579756 adjacent_seeds.pop();
1524
1525 1579756 C.initialize_from_mesh_tetrahedron(
1526
1/2
✓ Branch 1 taken 789878 times.
✗ Branch 2 not taken.
 1579756 mesh_, current_tet_, symbolic_, vertex_weight
1527 );
1528
1529
1/2
✓ Branch 1 taken 789878 times.
✗ Branch 2 not taken.
 1579756 intersect_cell_cell(
1530 current_seed_, C
1531 );
1532
1533
1/2
✓ Branch 2 taken 789878 times.
✗ Branch 3 not taken.
 1579756 action(
1534 current_seed_, current_tet_,
1535 current_polyhedron()
1536 );
1537
1538 // Propagate to adjacent tets and adjacent seeds
1539 // Iterate on the vertices of the cell (remember:
1540 // the cell is represented in dual form)
1541 27427544 for(index_t v = 0;
1542
3/4
✓ Branch 2 taken 13713772 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 12923894 times.
✓ Branch 5 taken 789878 times.
 27427544 v < current_polyhedron().max_v(); ++v
1543 ) {
1544
1545 // Skip clipping planes that are no longer
1546 // connected to a cell facet.
1547 41504134 if(
1548
1/2
✓ Branch 2 taken 12923894 times.
✗ Branch 3 not taken.
 25847788 current_polyhedron().vertex_triangle(v)
1549
2/2
✓ Branch 0 taken 7828173 times.
✓ Branch 1 taken 5095721 times.
 25847788 == -1
1550 ) {
1551 15656346 continue;
1552 }
1553
1554 signed_index_t id =
1555
1/2
✓ Branch 2 taken 5095721 times.
✗ Branch 3 not taken.
 10191442 current_polyhedron().vertex_id(v);
1556
2/2
✓ Branch 0 taken 2883354 times.
✓ Branch 1 taken 2212367 times.
 10191442 if(id > 0) {
1557 // Propagate to adjacent seed
1558 5766708 index_t neigh_s = index_t(id - 1);
1559
3/4
✓ Branch 1 taken 2883354 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 728438 times.
✓ Branch 4 taken 2154916 times.
 5766708 if(seed_stamp[neigh_s] != current_tet_) {
1560
1/2
✓ Branch 1 taken 728438 times.
✗ Branch 2 not taken.
 1456876 seed_stamp[neigh_s] = current_tet_;
1561
1/2
✓ Branch 1 taken 728438 times.
✗ Branch 2 not taken.
 1456876 adjacent_seeds.push(neigh_s);
1562 }
1563
2/2
✓ Branch 0 taken 2017011 times.
✓ Branch 1 taken 195356 times.
 4424734 } else if(id < 0) {
1564 // id==0 corresponds to facet on boundary
1565 // (skipped)
1566 // id<0 corresponds to adjacent tet index
1567
1568 // Propagate to adjacent tet
1569 4034022 signed_index_t neigh_t = -id - 1;
1570 4034022 if(
1571 neigh_t >=
1572
2/2
✓ Branch 0 taken 1946233 times.
✓ Branch 1 taken 70778 times.
 4034022 signed_index_t(tets_begin_) &&
1573
2/2
✓ Branch 0 taken 1875465 times.
✓ Branch 1 taken 70768 times.
 3892466 neigh_t < signed_index_t(tets_end_) &&
1574
1/2
✓ Branch 0 taken 1875465 times.
✗ Branch 1 not taken.
 3750930 neigh_t != signed_index_t(current_tet_)
1575 ) {
1576
2/2
✓ Branch 1 taken 60800 times.
✓ Branch 2 taken 1814665 times.
 3750930 if(!tet_is_marked[
1577 3750930 index_t(neigh_t)-tets_begin_
1578 ]) {
1579 tet_is_marked[
1580 121600 index_t(neigh_t)-tets_begin_
1581 121600 ] = true;
1582
1/2
✓ Branch 1 taken 60800 times.
✗ Branch 2 not taken.
 121600 adjacent_tets.push(
1583 243200 TetSeed(
1584 index_t(neigh_t),
1585 current_seed_
1586 )
1587 );
1588 }
1589 }
1590 }
1591 }
1592 }
1593 }
1594 }
1595 }
1596 1280 current_polyhedron_ = nullptr;
1597 1280 }
1598
1599
1600 /**
1601 * \brief Low-level API of Restricted Voronoi Diagram traversal
1602 * with connected components priority.
1603 * \details Client code may use for_each_cell() instead.
1604 * This version of the algorithm traverses the RVD and ensures that
1605 * the group of subfacets that belong to the same restricted Voronoi
1606 * cell will be traversed consecutively. It is used by the algorithm
1607 * that computes the final surface in CVT (i.e., the dual of the
1608 * connected components).
1609 * \note This function is less efficient than
1610 * compute_volumetric_with_seeds_priority() but
1611 * is required by some traversals that need to be done in that order.
1612 * \tparam ACTION needs to implement:
1613 * operator()(index_t v, index_t t, const Polyhedron& C) const
1614 * where v denotes the index of the current Voronoi cell
1615 * (or Delaunay vertex), c the index of the current tetrahedron
1616 * and C the computed intersection between the Voronoi cell of
1617 * v and tetrahedron t.
1618 */
1619 template <class ACTION>
1620 8 inline void compute_volumetric_with_cnx_priority(
1621 const ACTION& action
1622 ) {
1623
1624 8 if(
1625
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 8 tets_begin_ == UNSPECIFIED_RANGE &&
1626
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 8 tets_end_ == UNSPECIFIED_RANGE
1627 ) {
1628 8 tets_begin_ = 0;
1629 8 tets_end_ = mesh_->cells.nb();
1630 }
1631
1632
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 8 geo_assert(tets_begin_ != UNSPECIFIED_RANGE);
1633
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 8 geo_assert(tets_end_ != UNSPECIFIED_RANGE);
1634
1635 8 current_polyhedron_ = nullptr;
1636
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 init_get_neighbors();
1637
1638
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 std::deque<TetSeed> adjacent_seeds;
1639
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 std::stack<index_t> adjacent_tets;
1640
1641 static constexpr index_t NO_STAMP = index_t(-1);
1642 8 GEO::vector<index_t> tet_stamp(
1643
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 tets_end_ - tets_begin_, NO_STAMP
1644 );
1645
1646 16 TetSeedMarking visited(
1647
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 8 tets_end_ - tets_begin_, delaunay_->nb_vertices()
1648 );
1649
1650 // Yes, facet_seed_marking_ points to the TetSeedMarking,
1651 // (TetSeedMarking is typedef-ed as FacetSeedMarking),
1652 // ugly I know... to be revised.
1653 8 facet_seed_marking_ = &visited;
1654
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 8 Polyhedron C(dimension());
1655 8 current_polyhedron_ = &C;
1656
1657 8 current_connected_component_ = 0;
1658 // index_t C_index = tets_end_ + 1; // Unused (see comment later)
1659
1660
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 GEO::Attribute<double> vertex_weight;
1661
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 vertex_weight.bind_if_is_defined(
1662
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 16 mesh_->vertices.attributes(),"weight"
1663 );
1664
1665 // The algorithm propagates along both the facet-graph of
1666 // the surface and the 1-skeleton of the Delaunay triangulation,
1667 // and computes all the relevant intersections between
1668 // each Voronoi cell and facet.
1669
2/2
✓ Branch 0 taken 1536 times.
✓ Branch 1 taken 4 times.
 3080 for(index_t t = tets_begin_; t < tets_end_; ++t) {
1670
3/4
✓ Branch 1 taken 1536 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 1532 times.
 3072 if(tet_stamp[t - tets_begin_] == NO_STAMP) {
1671 8 current_tet_ = t;
1672
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 current_seed_ = find_seed_near_tet(t);
1673
1674
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 adjacent_seeds.push_back(
1675 8 TetSeed(current_tet_, current_seed_)
1676 );
1677
1678 // Propagate along the Delaunay-graph.
1679
2/2
✓ Branch 1 taken 35552 times.
✓ Branch 2 taken 4 times.
 71112 while(!adjacent_seeds.empty()) {
1680 // Yes, f, because TetSeed is typedef-ed as FacetSeed
1681 71104 current_tet_ = adjacent_seeds.front().f;
1682 71104 current_seed_ = adjacent_seeds.front().seed;
1683 71104 adjacent_seeds.pop_front();
1684 92086 if(
1685
1/2
✓ Branch 1 taken 35552 times.
✗ Branch 2 not taken.
 71104 tet_stamp[current_tet_ - tets_begin_] ==
1686
2/2
✓ Branch 0 taken 10491 times.
✓ Branch 1 taken 25061 times.
 71104 current_seed_
1687 ) {
1688 20982 continue;
1689 }
1690
1691
3/4
✓ Branch 1 taken 25061 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 23807 times.
✓ Branch 4 taken 1254 times.
 50122 if(visited.is_marked(current_tet_, current_seed_)) {
1692 47614 continue;
1693 }
1694
1695 2508 connected_component_changed_ = true;
1696
1/2
✓ Branch 1 taken 1254 times.
✗ Branch 2 not taken.
 2508 adjacent_tets.push(current_tet_);
1697 5016 tet_stamp[current_tet_ - tets_begin_] =
1698
1/2
✓ Branch 1 taken 1254 times.
✗ Branch 2 not taken.
 2508 current_seed_;
1699
2/2
✓ Branch 1 taken 19547 times.
✓ Branch 2 taken 1254 times.
 44110 while(!adjacent_tets.empty()) {
1700 39094 current_tet_ = adjacent_tets.top();
1701 39094 adjacent_tets.pop();
1702
1703 // Copy the current tet from the Mesh into
1704 // RestrictedVoronoiDiagram's Polyhedron
1705 // data structure (gathers all the necessary
1706 // information)
1707
1708 39094 C.initialize_from_mesh_tetrahedron(
1709
1/2
✓ Branch 1 taken 19547 times.
✗ Branch 2 not taken.
 39094 mesh_, current_tet_, symbolic_, vertex_weight
1710 );
1711
1712 // Note: difference with
1713 // compute_surfacic_with_cnx_priority():
1714 // Since intersect_cell_cell() overwrites C, we
1715 // need to initialize C from the mesh for each
1716 // visited (tet,seed) pair (and the test for
1717 // current_tet_ change with C_index is not
1718 // used here).
1719 // C_index = current_tet_;
1720
1721
1/2
✓ Branch 1 taken 19547 times.
✗ Branch 2 not taken.
 39094 intersect_cell_cell(current_seed_, C);
1722
1/2
✓ Branch 2 taken 19547 times.
✗ Branch 3 not taken.
 39094 action(
1723 current_seed_, current_tet_, current_polyhedron()
1724 );
1725 39094 connected_component_changed_ = false;
1726
1727 39094 bool touches_RVC_border = false;
1728
1729 // Propagate to adjacent tets and adjacent seeds
1730 657258 for(index_t v = 0;
1731
3/4
✓ Branch 2 taken 328629 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 309082 times.
✓ Branch 5 taken 19547 times.
 657258 v < current_polyhedron().max_v(); ++v
1732 ) {
1733
1734 // Skip clipping planes that are no longer
1735 // connected to a cell facet.
1736 984270 if(
1737
1/2
✓ Branch 2 taken 309082 times.
✗ Branch 3 not taken.
 618164 current_polyhedron().vertex_triangle(v)
1738
2/2
✓ Branch 0 taken 183053 times.
✓ Branch 1 taken 126029 times.
 618164 == -1
1739 ) {
1740 366106 continue;
1741 }
1742
1743 signed_index_t id =
1744
1/2
✓ Branch 2 taken 126029 times.
✗ Branch 3 not taken.
 252058 current_polyhedron().vertex_id(v);
1745
2/2
✓ Branch 0 taken 50032 times.
✓ Branch 1 taken 75997 times.
 252058 if(id < 0) {
1746 // id == 0 corresponds to facet on boundary
1747 // (skipped)
1748 // id < 0 corresponds to adjacent tet index
1749 100064 signed_index_t s_neigh_t = -id-1;
1750 100064 if(
1751
1/2
✓ Branch 0 taken 50032 times.
✗ Branch 1 not taken.
 100064 s_neigh_t >= signed_index_t(tets_begin_)
1752 100064 &&
1753
1/2
✓ Branch 0 taken 50032 times.
✗ Branch 1 not taken.
 100064 s_neigh_t < signed_index_t(tets_end_)
1754 ) {
1755
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 50032 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 100064 geo_debug_assert(
1756 s_neigh_t !=
1757 signed_index_t(current_tet_)
1758 );
1759 100064 index_t neigh_t = index_t(s_neigh_t);
1760 100064 if(
1761
1/2
✓ Branch 1 taken 50032 times.
✗ Branch 2 not taken.
 100064 tet_stamp[neigh_t - tets_begin_] !=
1762
2/2
✓ Branch 0 taken 18293 times.
✓ Branch 1 taken 31739 times.
 100064 current_seed_
1763 ) {
1764 73172 tet_stamp[neigh_t - tets_begin_] =
1765
1/2
✓ Branch 1 taken 18293 times.
✗ Branch 2 not taken.
 36586 current_seed_;
1766
1/2
✓ Branch 1 taken 18293 times.
✗ Branch 2 not taken.
 36586 adjacent_tets.push(neigh_t);
1767 }
1768 }
1769
2/2
✓ Branch 0 taken 71101 times.
✓ Branch 1 taken 4896 times.
 151994 } else if(id > 0) {
1770 142202 index_t neigh_s = index_t(id-1);
1771 142202 touches_RVC_border = true;
1772 142202 TetSeed ts(current_tet_, neigh_s);
1773
3/4
✓ Branch 1 taken 71101 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 35548 times.
✓ Branch 4 taken 35553 times.
 142202 if(!visited.is_marked(ts)) {
1774
1/2
✓ Branch 1 taken 35548 times.
✗ Branch 2 not taken.
 71096 adjacent_seeds.push_back(ts);
1775 }
1776 }
1777
1778 }
1779
2/2
✓ Branch 0 taken 19538 times.
✓ Branch 1 taken 9 times.
 39094 if(touches_RVC_border) {
1780
1/2
✓ Branch 1 taken 19538 times.
✗ Branch 2 not taken.
 39076 visited.mark(
1781 78152 TetSeed(current_tet_, current_seed_),
1782 current_connected_component_
1783 );
1784 }
1785 }
1786 2508 ++current_connected_component_;
1787 }
1788 }
1789 }
1790 8 facet_seed_marking_ = nullptr;
1791 8 }
1792
1793
1794 public:
1795 /**
1796 * \brief Tests whether a (facet,seed) couple was visited.
1797 * \param[in] f index of the facet
1798 * \param[in] s index of the seed
1799 */
1800 bool facet_seed_is_visited(index_t f, index_t s) const {
1801 geo_debug_assert(facet_seed_marking_ != nullptr);
1802 return facet_seed_marking_->is_marked(FacetSeed(f, s));
1803 }
1804
1805 /**
1806 * \brief Gets the index of the connected component associated
1807 * with a (facet,seed).
1808 * \param[in] f index of the facet
1809 * \param[in] s index of the seed
1810 * \return the index of the connected component or -1 if the
1811 * (\p f, \p s) couple was not visited already
1812 */
1813 836556 index_t get_facet_seed_connected_component(index_t f, index_t s) const {
1814
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 418278 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 836556 geo_debug_assert(facet_seed_marking_ != nullptr);
1815
1/2
✓ Branch 1 taken 418278 times.
✗ Branch 2 not taken.
 836556 return facet_seed_marking_->get_connected_component(
1816 1673112 FacetSeed(f, s)
1817 1673112 );
1818 }
1819
1820 /**
1821 * \brief Tests whether the current connected component changed.
1822 */
1823 358780 bool connected_component_changed() const {
1824 358780 return connected_component_changed_;
1825 }
1826
1827 /**
1828 * \brief Gets the index of the current connected component.
1829 */
1830 418278 index_t current_connected_component() const {
1831 418278 return current_connected_component_;
1832 }
1833
1834 protected:
1835 /**
1836 * \brief Low-level API of Restricted Voronoi Diagram traversal
1837 * with connected components priority.
1838 * \details Client code may use for_each_facet(),for_each_triangle() or
1839 * for_each_primal_triangle() instead.
1840 * This version of the algorithm traverses the RVD and ensures that
1841 * the group of subfacets that belong to the same restricted Voronoi
1842 * cell will be traversed consecutively. It is used by the algorithm
1843 * that computes the final surface in CVT (i.e., the dual of the
1844 * connected components).
1845 * \note This function is less efficient than
1846 * compute_surfacic_with_seeds_priority() but
1847 * is required by some traversals that need to be done in that order.
1848 * \tparam ACTION needs to implement:
1849 * operator()(index_t v, index_t f, const Polygon& P) const
1850 * where v denotes the index of the current Voronoi cell
1851 * (or Delaunay vertex), f the index of the current facet
1852 * and P the computed intersection between the Voronoi cell of
1853 * v and facet f.
1854 */
1855 template <class ACTION>
1856 14 inline void compute_surfacic_with_cnx_priority(
1857 const ACTION& action
1858 ) {
1859
1860 14 if(
1861
1/2
✓ Branch 0 taken 7 times.
✗ Branch 1 not taken.
 14 facets_begin_ == UNSPECIFIED_RANGE &&
1862
1/2
✓ Branch 0 taken 7 times.
✗ Branch 1 not taken.
 14 facets_end_ == UNSPECIFIED_RANGE
1863 ) {
1864 14 facets_begin_ = 0;
1865 14 facets_end_ = mesh_->facets.nb();
1866 }
1867
1868 14 current_polygon_ = nullptr;
1869
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 14 init_get_neighbors();
1870
1871
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 14 std::deque<FacetSeed> adjacent_seeds;
1872
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 14 std::stack<index_t> adjacent_facets;
1873
1874 static constexpr index_t NO_STAMP = index_t(-1);
1875 14 GEO::vector<index_t> facet_stamp(
1876
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 14 facets_end_ - facets_begin_, NO_STAMP
1877 );
1878
1879 28 FacetSeedMarking visited(
1880
1/2
✓ Branch 2 taken 7 times.
✗ Branch 3 not taken.
 14 facets_end_ - facets_begin_, delaunay_->nb_vertices()
1881 );
1882
1883 14 facet_seed_marking_ = &visited;
1884
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 14 Polygon F;
1885 14 current_connected_component_ = 0;
1886 14 index_t F_index = facets_end_ + 1;
1887
1888
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 14 GEO::Attribute<double> vertex_weight;
1889
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 14 vertex_weight.bind_if_is_defined(
1890
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 28 mesh_->vertices.attributes(),"weight"
1891 );
1892
1893 // The algorithm propagates along both the facet-graph of
1894 // the surface and the 1-skeleton of the Delaunay triangulation,
1895 // and computes all the relevant intersections between
1896 // each Voronoi cell and facet.
1897
2/2
✓ Branch 0 taken 70265 times.
✓ Branch 1 taken 7 times.
 140544 for(index_t f = facets_begin_; f < facets_end_; ++f) {
1898
3/4
✓ Branch 1 taken 70265 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 7 times.
✓ Branch 4 taken 70258 times.
 140530 if(facet_stamp[f - facets_begin_] == NO_STAMP) {
1899 14 current_facet_ = f;
1900
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 14 current_seed_ = find_seed_near_facet(f);
1901
1902
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 14 adjacent_seeds.push_back(
1903 14 FacetSeed(current_facet_, current_seed_)
1904 );
1905
1906 // Propagate along the Delaunay-graph.
1907 14 while(
1908
2/2
✓ Branch 1 taken 178845 times.
✓ Branch 2 taken 7 times.
 357704 !adjacent_seeds.empty()
1909 ) {
1910 357690 current_facet_ = adjacent_seeds.front().f;
1911 357690 current_seed_ = adjacent_seeds.front().seed;
1912 357690 adjacent_seeds.pop_front();
1913 506050 if(
1914
1/2
✓ Branch 1 taken 178845 times.
✗ Branch 2 not taken.
 357690 facet_stamp[current_facet_ - facets_begin_] ==
1915
2/2
✓ Branch 0 taken 74180 times.
✓ Branch 1 taken 104665 times.
 357690 current_seed_
1916 ) {
1917 148360 continue;
1918 }
1919
1920
3/4
✓ Branch 1 taken 104665 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 69667 times.
✓ Branch 4 taken 34998 times.
 209330 if(visited.is_marked(current_facet_, current_seed_)) {
1921 139334 continue;
1922 }
1923
1924 69996 connected_component_changed_ = true;
1925
1/2
✓ Branch 1 taken 34998 times.
✗ Branch 2 not taken.
 69996 adjacent_facets.push(current_facet_);
1926 139992 facet_stamp[current_facet_ - facets_begin_] =
1927
1/2
✓ Branch 1 taken 34998 times.
✗ Branch 2 not taken.
 69996 current_seed_;
1928
2/2
✓ Branch 1 taken 179390 times.
✓ Branch 2 taken 34998 times.
 428776 while(!adjacent_facets.empty()) {
1929 358780 current_facet_ = adjacent_facets.top();
1930 358780 adjacent_facets.pop();
1931
1932 // Copy the current facet from the Mesh into
1933 // RestrictedVoronoiDiagram's Polygon data structure
1934 // (gathers all the necessary information)
1935
2/2
✓ Branch 0 taken 164986 times.
✓ Branch 1 taken 14404 times.
 358780 if(F_index != current_facet_) {
1936 329972 F.initialize_from_mesh_facet(
1937
1/2
✓ Branch 1 taken 164986 times.
✗ Branch 2 not taken.
 329972 mesh_, current_facet_, symbolic_,
1938 vertex_weight
1939 );
1940 329972 F_index = current_facet_;
1941 }
1942
1943
1/2
✓ Branch 1 taken 179390 times.
✗ Branch 2 not taken.
 358780 current_polygon_ = intersect_cell_facet(
1944 current_seed_, F
1945 );
1946
1/2
✓ Branch 2 taken 179390 times.
✗ Branch 3 not taken.
 358780 action(
1947 current_seed_, current_facet_, current_polygon()
1948 );
1949 358780 connected_component_changed_ = false;
1950
1951 358780 bool touches_RVC_border = false;
1952
1953 // Propagate to adjacent facets and adjacent seeds
1954 1792796 for(index_t v = 0;
1955
3/4
✓ Branch 2 taken 896398 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 717008 times.
✓ Branch 5 taken 179390 times.
 1792796 v < current_polygon().nb_vertices(); v++
1956 ) {
1957
1/2
✓ Branch 2 taken 717008 times.
✗ Branch 3 not taken.
 1434016 const Vertex& ve = current_polygon().vertex(v);
1958 1434016 signed_index_t s_neigh_f = ve.adjacent_facet();
1959 1434016 if(
1960
2/2
✓ Branch 0 taken 357820 times.
✓ Branch 1 taken 359188 times.
 1434016 s_neigh_f >= signed_index_t(facets_begin_)
1961 715640 &&
1962
1/2
✓ Branch 0 taken 357820 times.
✗ Branch 1 not taken.
 715640 s_neigh_f < signed_index_t(facets_end_)
1963 ) {
1964
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 357820 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 715640 geo_debug_assert(
1965 s_neigh_f !=
1966 signed_index_t(current_facet_)
1967 );
1968 715640 index_t neigh_f = index_t(s_neigh_f);
1969 715640 if(
1970
1/2
✓ Branch 1 taken 357820 times.
✗ Branch 2 not taken.
 715640 facet_stamp[neigh_f - facets_begin_] !=
1971
2/2
✓ Branch 0 taken 144392 times.
✓ Branch 1 taken 213428 times.
 715640 current_seed_
1972 ) {
1973 577568 facet_stamp[neigh_f - facets_begin_] =
1974
1/2
✓ Branch 1 taken 144392 times.
✗ Branch 2 not taken.
 288784 current_seed_;
1975
1/2
✓ Branch 1 taken 144392 times.
✗ Branch 2 not taken.
 288784 adjacent_facets.push(neigh_f);
1976 }
1977 }
1978 1434016 signed_index_t neigh_s = ve.adjacent_seed();
1979
2/2
✓ Branch 0 taken 357676 times.
✓ Branch 1 taken 359332 times.
 1434016 if(neigh_s != -1) {
1980 715352 touches_RVC_border = true;
1981 715352 FacetSeed fs(
1982 current_facet_, index_t(neigh_s)
1983 );
1984
3/4
✓ Branch 1 taken 357676 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 178838 times.
✓ Branch 4 taken 178838 times.
 715352 if(!visited.is_marked(fs)) {
1985
1/2
✓ Branch 1 taken 178838 times.
✗ Branch 2 not taken.
 357676 adjacent_seeds.push_back(fs);
1986 }
1987 }
1988 }
1989
2/2
✓ Branch 0 taken 144015 times.
✓ Branch 1 taken 35375 times.
 358780 if(touches_RVC_border) {
1990
1/2
✓ Branch 1 taken 144015 times.
✗ Branch 2 not taken.
 288030 visited.mark(
1991 576060 FacetSeed(current_facet_, current_seed_),
1992 current_connected_component_
1993 );
1994 }
1995 }
1996 69996 ++current_connected_component_;
1997 }
1998 }
1999 }
2000 14 facet_seed_marking_ = nullptr;
2001 14 }
2002
2003 /**
2004 * \brief Finds a seed near a given facet.
2005 * \param[in] f index of the facet in the mesh
2006 * \return the index of a Voronoi seed such that there is a
2007 * non-empty intersection between the Voronoi cell
2008 * of the seed and facet \p f.
2009 */
2010 3148 index_t find_seed_near_facet(index_t f) {
2011 3148 const double* p = mesh_->vertices.point_ptr(
2012 3148 mesh_->facets.vertex(f,0)
2013 );
2014 3148 return find_seed_near_point(p);
2015 }
2016
2017 /**
2018 * \brief Finds a seed near a given tetrahedron.
2019 * \param[in] t index of the tetrahedron in the mesh
2020 * \return the index of a Voronoi seed such that there is a
2021 * non-empty intersection between the Voronoi cell
2022 * of the seed and tetrahedron \p t.
2023 */
2024 1288 index_t find_seed_near_tet(index_t t) {
2025 1288 index_t v = mesh_->cells.tet_vertex(t, 0);
2026 1288 const double* p = mesh_->vertices.point_ptr(v);
2027 1288 return find_seed_near_point(p);
2028 }
2029
2030 /**
2031 * \brief Finds a seed near a given point.
2032 * \param[in] p pointer to the coordinates of the point
2033 * \return the index of a Voronoi seed such that its
2034 * Voronoi cell contains the point \p p.
2035 */
2036 4436 index_t find_seed_near_point(const double* p) {
2037 // In order to be compatible with the symbolic
2038 // perturbation, if the nearest neighbor is
2039 // non-unique, we need to return the one of
2040 // lowest index (because in case of several seeds
2041 // at equal distance, the one of lowest index
2042 // is guaranteed to have the facet in its Voronoi
2043 // cell from the point of view of symbolic
2044 // perturbation).
2045
3/4
✓ Branch 0 taken 48 times.
✓ Branch 1 taken 2170 times.
✓ Branch 2 taken 48 times.
✗ Branch 3 not taken.
 4436 if(exact_ && delaunay_nn_ != nullptr) {
2046 // TODO: may need more than 10
2047 index_t neighbors[10];
2048 double neighbors_sq_dist[10];
2049 96 index_t nb = 10;
2050
2/2
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 45 times.
 96 if(delaunay_nn_->nb_vertices() < nb) {
2051 6 nb = delaunay_nn_->nb_vertices();
2052 }
2053
2/4
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 48 times.
✗ Branch 5 not taken.
 96 delaunay_nn_->nn_search()->get_nearest_neighbors(
2054 nb, p, neighbors, neighbors_sq_dist
2055 );
2056 96 index_t nearest = neighbors[0];
2057 96 double min_d = neighbors_sq_dist[0];
2058
1/2
✓ Branch 0 taken 63 times.
✗ Branch 1 not taken.
 126 for(index_t i = 1; i < nb; ++i) {
2059
2/2
✓ Branch 0 taken 48 times.
✓ Branch 1 taken 15 times.
 126 if(neighbors_sq_dist[i] != min_d) {
2060 96 break;
2061 }
2062
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 7 times.
 30 if(neighbors[i] < nearest) {
2063 16 nearest = neighbors[i];
2064 }
2065 }
2066 96 return nearest;
2067 }
2068
2069 4340 return delaunay_->nearest_vertex(p);
2070 }
2071
2072 /**
2073 * @}
2074 * \name Clipping for surfacic mode
2075 * @{
2076 */
2077
2078 /**
2079 * \brief Swaps two pointers between two polygons.
2080 * \details Used by re-entrant Sutherlang-Hogdman clipping.
2081 */
2082 93529484 void swap_polygons(Polygon*& ping, Polygon*& pong) {
2083
4/4
✓ Branch 0 taken 30203107 times.
✓ Branch 1 taken 16561635 times.
✓ Branch 2 taken 8955021 times.
✓ Branch 3 taken 21248086 times.
 93529484 if(ping != &P1 && ping != &P2) {
2084 // First clipping operation, ping points to F
2085 // (current facet copied)
2086 17910042 ping = &P2;
2087 17910042 pong = &P1;
2088 } else {
2089 75619442 std::swap(ping, pong);
2090 }
2091 93529484 }
2092
2093 /**
2094 * \brief Computes the intersection between the Voronoi cell
2095 * of a seed and a facet.
2096 * \param[in] seed the index of the seed
2097 * \param[in] F the facet represented as a Polygon
2098 * \details The result is provided in current_polygon_
2099 */
2100 19521200 Polygon* intersect_cell_facet(index_t seed, Polygon& F) {
2101 19521200 intersections_.clear();
2102
2103 // Initialize ping-pong pointers for Sutherland-Hodgman
2104 // re-entrant clipping and copy current facet into 'ping' buffer.
2105 19521200 Polygon* ping = &F;
2106 19521200 Polygon* pong = &P2;
2107
2108 // Clip current facet by current Voronoi cell (associated with seed)
2109
1/2
✓ Branch 0 taken 9760600 times.
✗ Branch 1 not taken.
 19521200 if(delaunay_nn_ != nullptr) {
2110
1/2
✓ Branch 1 taken 9760600 times.
✗ Branch 2 not taken.
 19521200 clip_by_cell_SR(seed, ping, pong); // "Security Radius" mode.
2111 } else {
2112 ✗ clip_by_cell(seed, ping, pong); // Standard mode.
2113 }
2114
2115 19521200 return ping; // Yes, 'ping', and not 'pong'
2116 // see comments in clip_by_cell()
2117 }
2118
2119 /**
2120 * \brief Computes the intersection between the Voronoi cell of a
2121 * vertex and the Mesh 'ping'.
2122 *
2123 * \details The result is returned in \p ping (Note that
2124 * \p ping and \p pong are references, and that they are swapped
2125 * after each bisector clipping, this is why the final result
2126 * is in \p ping (and not in \p pong).
2127 * This version uses the Security Radius algorithm.
2128 *
2129 * \param[in] i index of the vertex that defines the Voronoi cell
2130 * \param[in,out] ping the input polygon. On exit, contains the result.
2131 * \param[out] pong a buffer used to implement reentrant clipping.
2132 * Its content is modified by the function.
2133 */
2134 19521200 void clip_by_cell_SR(index_t i, Polygon*& ping, Polygon*& pong) {
2135 // 'Security radius' mode.
2136 // Note: the vertices of the neighborhood are returned in
2137 // increasing distance to pi. We stop the clippings as soon as the
2138 // 'security radius' is reached.
2139 19521200 const double* geo_restrict pi = delaunay_->vertex_ptr(i);
2140 geo_assume_aligned(pi, geo_dim_alignment(DIM));
2141
2142 19521200 index_t jj = 0;
2143 19521200 index_t prev_nb_neighbors = 0;
2144 19521200 neighbors_.resize(0);
2145
2/2
✓ Branch 2 taken 9786301 times.
✓ Branch 3 taken 96 times.
 19572794 while(neighbors_.size() < delaunay_nn_->nb_vertices() - 1) {
2146
2147
1/2
✓ Branch 1 taken 9786301 times.
✗ Branch 2 not taken.
 19572602 delaunay_nn_->get_neighbors(i, neighbors_);
2148
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 9786301 times.
 19572602 if(neighbors_.size() == 0) {
2149 19521008 return;
2150 }
2151
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 9786301 times.
 19572602 if(prev_nb_neighbors == neighbors_.size()) {
2152 ✗ return;
2153 }
2154
2155
2/2
✓ Branch 1 taken 56499924 times.
✓ Branch 2 taken 51119 times.
 113102086 for(; jj < neighbors_.size(); jj++) {
2156
1/2
✓ Branch 1 taken 56499924 times.
✗ Branch 2 not taken.
 112999848 index_t j = neighbors_[jj];
2157 112999848 double R2 = 0.0;
2158
3/4
✓ Branch 1 taken 259991469 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 203491545 times.
✓ Branch 4 taken 56499924 times.
 519982938 for(index_t k = 0; k < ping->nb_vertices(); k++) {
2159 geo_decl_aligned(double dik);
2160
1/2
✓ Branch 1 taken 203491545 times.
✗ Branch 2 not taken.
 406983090 const double* geo_restrict pk = ping->vertex(k).point();
2161 geo_assume_aligned(pk, geo_dim_alignment(DIM));
2162
1/2
✓ Branch 2 taken 203491545 times.
✗ Branch 3 not taken.
 406983090 dik = GEO::Geom::distance2(pi, pk, dimension());
2163 406983090 R2 = std::max(R2, dik);
2164 }
2165 geo_decl_aligned(double dij);
2166
1/2
✓ Branch 1 taken 56499924 times.
✗ Branch 2 not taken.
 112999848 const double* geo_restrict pj = delaunay_->vertex_ptr(j);
2167 geo_assume_aligned(pj, geo_dim_alignment(DIM));
2168
1/2
✓ Branch 2 taken 56499924 times.
✗ Branch 3 not taken.
 112999848 dij = GEO::Geom::distance2(pi, pj, dimension());
2169 // A little bit more than 4, because when
2170 // exact predicates are used, we need to
2171 // include tangent bisectors in the computation.
2172
2/2
✓ Branch 0 taken 9735182 times.
✓ Branch 1 taken 46764742 times.
 112999848 if(dij > 4.1 * R2) {
2173 19470364 return;
2174 }
2175
1/2
✓ Branch 1 taken 46764742 times.
✗ Branch 2 not taken.
 93529484 clip_by_plane(*ping, *pong, i, j);
2176 93529484 swap_polygons(ping, pong);
2177 }
2178
2179
2/2
✓ Branch 0 taken 25322 times.
✓ Branch 1 taken 25797 times.
 102238 if(!check_SR_) {
2180 50644 return;
2181 }
2182
2183 51594 index_t nb_neighbors = neighbors_.size();
2184 51594 prev_nb_neighbors = nb_neighbors;
2185
2186
2/2
✓ Branch 0 taken 25770 times.
✓ Branch 1 taken 27 times.
 51594 if(nb_neighbors > 8) {
2187 51540 nb_neighbors += nb_neighbors / 8;
2188 } else {
2189 54 nb_neighbors++;
2190 }
2191
2192 103188 nb_neighbors = std::min(
2193 nb_neighbors,
2194 51594 delaunay_nn_->nb_vertices() - 1
2195 );
2196
2197
1/2
✓ Branch 1 taken 25797 times.
✗ Branch 2 not taken.
 51594 delaunay_nn_->enlarge_neighborhood(i, nb_neighbors);
2198 }
2199 }
2200
2201 /**
2202 * \brief Computes the intersection between a Voronoi cell
2203 * and a polygon.
2204 *
2205 * \details The Voronoi cell is determined by vertex \p i and
2206 * the input polygon is in \p ping. The result is returned
2207 * in \p ping (Note that
2208 * \p ping and \p pong are references, and that they are swapped
2209 * after each bisector clipping, this is why the final result
2210 * is in \p ping (and not in \p pong).
2211 *
2212 * \param[in] i index of the vertex that defines the Voronoi cell
2213 * \param[in,out] ping the input polygon. On exit, contains the result.
2214 * \param[out] pong a buffer used to implement reentrant clipping.
2215 * Its content is modified by the function.
2216 */
2217 ✗ void clip_by_cell(index_t i, Polygon*& ping, Polygon*& pong) {
2218 ✗ get_neighbors(i);
2219 ✗ for(index_t jj = 0; jj < neighbors_.size(); jj++) {
2220 ✗ index_t j = neighbors_[jj];
2221 ✗ clip_by_plane(*ping, *pong, i, j);
2222 ✗ swap_polygons(ping, pong);
2223 }
2224 ✗ }
2225
2226 /**
2227 * \brief Computes the intersection between a polygon and a half-space.
2228 *
2229 * \details The input polygon is in \p ping
2230 * and the half-space is determined by the positive side
2231 * of the bisector of segment [\p i,\p j] (the side of \p i).
2232 * The result is stored into the Polygon \p pong.
2233 *
2234 * \param[in] i index of the first extremity of the bisector
2235 * \param[in] j index of the second extremity of the bisector
2236 * \param[in] ping the input polygon
2237 * \param[out] pong \p ping clipped by the bisector
2238 */
2239
2240 93529484 void clip_by_plane(
2241 Polygon& ping, Polygon& pong,
2242 index_t i, index_t j
2243 ) {
2244 93529484 ping.clip_by_plane<DIM>(
2245 93529484 pong, intersections_, mesh_, delaunay_, i, j, exact_, symbolic_
2246 );
2247 93529484 }
2248
2249 /**
2250 * @}
2251 * \name Clipping for volumetric mode
2252 * @{
2253 */
2254
2255 public:
2256 /**
2257 * \brief Computes the intersection between a Voronoi cell
2258 * and a cell with radius of security or plain mode.
2259 * \param[in] seed the index of the seed that defines the Voronoi cell
2260 * \param[in,out] C the cell to be clipped
2261 */
2262 1618850 void intersect_cell_cell(index_t seed, Polyhedron& C) {
2263 // Clip current facet by current Voronoi cell (associated with seed)
2264
1/2
✓ Branch 0 taken 809425 times.
✗ Branch 1 not taken.
 1618850 if(delaunay_nn_ != nullptr) {
2265 1618850 clip_by_cell_SR(seed, C); // "Security Radius" mode.
2266 } else {
2267 ✗ clip_by_cell(seed, C); // Standard mode.
2268 }
2269 1618850 }
2270
2271 protected:
2272 /**
2273 * \brief Computes the intersection between a Voronoi cell
2274 * and a cell in radius-of-security mode.
2275 * \param[in] seed the index of the seed that defines the Voronoi cell
2276 * \param[in,out] C the cell to be clipped
2277 */
2278 1618850 void clip_by_cell_SR(index_t seed, Polyhedron& C) {
2279
2280 // 'Security radius' mode.
2281 // Note: the vertices of the neighborhood are returned in
2282 // increasing distance to pi. We stop the clippings as soon as the
2283 // 'security radius' is reached.
2284 1618850 const double* geo_restrict pi = delaunay_->vertex_ptr(seed);
2285 geo_assume_aligned(pi, geo_dim_alignment(DIM));
2286
2287 1618850 index_t jj = 0;
2288 1618850 index_t prev_nb_neighbors = 0;
2289 1618850 neighbors_.resize(0);
2290
2291
1/2
✓ Branch 2 taken 809425 times.
✗ Branch 3 not taken.
 1618850 while(neighbors_.size() < delaunay_nn_->nb_vertices() - 1) {
2292
2293
1/2
✓ Branch 1 taken 809425 times.
✗ Branch 2 not taken.
 1618850 delaunay_nn_->get_neighbors(seed, neighbors_);
2294
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 809425 times.
 1618850 if(neighbors_.size() == 0) {
2295 1618850 return;
2296 }
2297
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 809425 times.
 1618850 if(prev_nb_neighbors == neighbors_.size()) {
2298 ✗ return;
2299 }
2300
2301
2/2
✓ Branch 1 taken 10763443 times.
✓ Branch 2 taken 65533 times.
 21657952 for(; jj < neighbors_.size(); jj++) {
2302
1/2
✓ Branch 1 taken 10763443 times.
✗ Branch 2 not taken.
 21526886 index_t j = neighbors_[jj];
2303 21526886 double R2 = 0.0;
2304
3/4
✓ Branch 1 taken 125075892 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 114312449 times.
✓ Branch 4 taken 10763443 times.
 250151784 for(index_t k = 0; k < C.max_t(); ++k) {
2305
3/4
✓ Branch 1 taken 114312449 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 32051021 times.
✓ Branch 4 taken 82261428 times.
 228624898 if(!C.triangle_is_used(k)) {
2306 64102042 continue;
2307 }
2308 geo_decl_aligned(double dik);
2309 const double* geo_restrict pk =
2310
1/2
✓ Branch 1 taken 82261428 times.
✗ Branch 2 not taken.
 164522856 C.triangle_dual(k).point();
2311 geo_assume_aligned(pk, geo_dim_alignment(DIM));
2312
1/2
✓ Branch 2 taken 82261428 times.
✗ Branch 3 not taken.
 164522856 dik = GEO::Geom::distance2(pi, pk, dimension());
2313 164522856 R2 = std::max(R2, dik);
2314 }
2315 geo_decl_aligned(double dij);
2316
1/2
✓ Branch 1 taken 10763443 times.
✗ Branch 2 not taken.
 21526886 const double* geo_restrict pj = delaunay_->vertex_ptr(j);
2317 geo_assume_aligned(pj, geo_dim_alignment(DIM));
2318
1/2
✓ Branch 2 taken 10763443 times.
✗ Branch 3 not taken.
 21526886 dij = GEO::Geom::distance2(pi, pj, dimension());
2319 // A little bit more than 4, because when
2320 // exact predicates are used, we need to
2321 // include tangent bisectors in the computation.
2322
2/2
✓ Branch 0 taken 743892 times.
✓ Branch 1 taken 10019551 times.
 21526886 if(dij > 4.1 * R2) {
2323 1487784 return;
2324 }
2325
1/2
✓ Branch 1 taken 10019551 times.
✗ Branch 2 not taken.
 20039102 clip_by_plane(C, seed, j);
2326 }
2327
2328
1/2
✓ Branch 0 taken 65533 times.
✗ Branch 1 not taken.
 131066 if(!check_SR_) {
2329 131066 return;
2330 }
2331
2332 ✗ index_t nb_neighbors = neighbors_.size();
2333 ✗ prev_nb_neighbors = nb_neighbors;
2334
2335 ✗ if(nb_neighbors > 8) {
2336 ✗ nb_neighbors += nb_neighbors / 8;
2337 } else {
2338 ✗ nb_neighbors++;
2339 }
2340
2341 ✗ nb_neighbors = std::min(
2342 nb_neighbors,
2343 ✗ delaunay_nn_->nb_vertices() - 1
2344 );
2345 ✗ delaunay_nn_->enlarge_neighborhood(seed, nb_neighbors);
2346 }
2347 }
2348
2349 /**
2350 * \brief Computes the intersection between a Voronoi cell
2351 * and a cell in plain mode.
2352 * \param[in] seed the index of the seed that defines the Voronoi cell
2353 * \param[in,out] C the cell to be clipped
2354 */
2355 ✗ void clip_by_cell(index_t seed, Polyhedron& C) {
2356 ✗ get_neighbors(seed);
2357 // Check whether cell is empty (may happen with
2358 // power diagrams)
2359 ✗ if(neighbors_.size() == 0) {
2360 ✗ C.clear();
2361 }
2362 ✗ for(index_t jj = 0; jj < neighbors_.size(); jj++) {
2363 ✗ index_t j = neighbors_[jj];
2364 ✗ clip_by_plane(C, seed, j);
2365 }
2366 ✗ }
2367
2368 /**
2369 * \brief Computes the intersection between a Voronoi cell
2370 * and a half-space determined by a bisector.
2371 * \param[in,out] C cell to be clipped
2372 * \param[in] i index of the first extremity of the bisector
2373 * \param[in] j index of the second extremity of the bisector
2374 */
2375 20039102 void clip_by_plane(Polyhedron& C, index_t i, index_t j) {
2376 20039102 C.clip_by_plane<DIM>(
2377 20039102 mesh_, delaunay_, i, j, exact_, symbolic_
2378 );
2379 20039102 }
2380
2381 /**
2382 * @}
2383 * \name Optimized get neighbors
2384 * @{
2385 */
2386
2387 /**
2388 * \brief Creates the data structure for optimized get_neighbors()
2389 * function.
2390 *
2391 * \details This function is only used when the stored delaunay
2392 * triangulation is a traditional one. When the stored delaunay
2393 * triangulation is represented by a KdTree, function is not used.
2394 */
2395 3734 void init_get_neighbors() {
2396 // In dimension 3 (and if we do not used the ANN-based algorithm),
2397 // we can use the faster 'stamp-based' algorithm for finding the
2398 // neighbors.
2399
4/6
✓ Branch 1 taken 579 times.
✓ Branch 2 taken 1288 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 579 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1867 times.
 3734 if(delaunay_->dimension() == 3 && delaunay_->nb_cells() != 0) {
2400 ✗ cur_stamp_ = 0;
2401 ✗ stamp_.assign(delaunay_->nb_vertices(), NO_INDEX);
2402 }
2403 3734 }
2404
2405 /**
2406 * \brief Caches the neighbors of a Delaunay vertex.
2407 *
2408 * \details This function is only used when the stored delaunay
2409 * triangulation is a traditional one. When the stored delaunay
2410 * triangulation is represented by a KdTree, function is not used.
2411 */
2412 ✗ void get_neighbors(index_t v) {
2413 ✗ if(stamp_.size() == 0) {
2414 // Used in ANN mode and with higher dimensions.
2415 ✗ delaunay_->get_neighbors(v, neighbors_);
2416 } else {
2417 // Used in 3D mode with standard Delaunay.
2418 // The following loop replaces
2419 // delaunay_->get_neighbors(v,neighbors_) ;
2420 // (and makes the overall algorithm 10 to 30% more efficient)
2421 ✗ neighbors_.resize(0);
2422 ✗ index_t t = index_t(delaunay_->vertex_cell(v));
2423 do {
2424 ✗ index_t lv = delaunay_->index(t, v);
2425 ✗ for(index_t lw = 0; lw < delaunay_->cell_size(); lw++) {
2426 ✗ if(lw != lv) {
2427 ✗ index_t w = index_t(delaunay_->cell_vertex(t, lw));
2428 ✗ if(stamp_[w] != cur_stamp_) {
2429 ✗ stamp_[w] = cur_stamp_;
2430 ✗ neighbors_.push_back(w);
2431 }
2432 }
2433 }
2434 ✗ t = index_t(delaunay_->next_around_vertex(t, lv));
2435 ✗ } while(t != index_t(delaunay_->vertex_cell(v)));
2436 ✗ cur_stamp_++;
2437 }
2438 ✗ }
2439
2440 /** @} */
2441
2442 protected:
2443 GEO::Mesh* mesh_;
2444 Delaunay* delaunay_;
2445 GEO::Delaunay_NearestNeighbors* delaunay_nn_;
2446
2447 PointAllocator intersections_;
2448 Polygon* current_polygon_;
2449 Polygon P1, P2;
2450 GEO::vector<index_t> neighbors_;
2451 index_t current_facet_;
2452 index_t current_seed_;
2453 Polyhedron* current_polyhedron_;
2454 index_t current_tet_;
2455
2456 // For optimized get_neighbors().
2457 index_t cur_stamp_;
2458 GEO::vector<index_t> stamp_;
2459
2460 bool symbolic_;
2461 bool check_SR_;
2462 bool exact_;
2463
2464 coord_index_t dimension_;
2465
2466 static constexpr index_t UNSPECIFIED_RANGE = index_t(-1);
2467
2468 index_t facets_begin_;
2469 index_t facets_end_;
2470
2471 index_t tets_begin_;
2472 index_t tets_end_;
2473
2474 bool connected_components_priority_;
2475 FacetSeedMarking* facet_seed_marking_;
2476 bool connected_component_changed_;
2477 index_t current_connected_component_;
2478
2479 /**
2480 * \brief Forbids construction from copy.
2481 */
2482 RestrictedVoronoiDiagram(const thisclass& rhs) = delete;
2483
2484 /**
2485 * \brief Forbids assignment.
2486 */
2487 thisclass& operator= (const thisclass& rhs) = delete;
2488 };
2489 }
2490
2491 namespace GEO {
2492
2493 /**
2494 * \brief Symbolic representation of a RestrictedVoronoiDiagram vertex.
2495 */
2496 typedef GEOGen::SymbolicVertex SymbolicVertex;
2497 }
2498
2499 #endif
2500