GCC Code Coverage Report

Directory: ./
File: lib/geogram/voronoi/generic_RVD.h
Date: 2026-09-07 02:28:19
Exec Total Coverage
Lines: 448 526 85.2%
Functions: 84 749 11.2%
Branches: 348 1526 22.8%
Line Branch Exec Source
1 /*
2 * Copyright (c) 2000-2022 Inria
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 * * Redistributions in binary form must reproduce the above copyright notice,
11 * this list of conditions and the following disclaimer in the documentation
12 * and/or other materials provided with the distribution.
13 * * Neither the name of the ALICE Project-Team nor the names of its
14 * contributors may be used to endorse or promote products derived from this
15 * software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 * POSSIBILITY OF SUCH DAMAGE.
28 *
29 * Contact: Bruno Levy
30 *
31 * https://www.inria.fr/fr/bruno-levy
32 *
33 * Inria,
34 * Domaine de Voluceau,
35 * 78150 Le Chesnay - Rocquencourt
36 * FRANCE
37 *
38 */
39
40 #ifndef GEOGRAM_VORONOI_GENERIC_RVD
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 static coord_index_t dimension() {
88 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 172 RestrictedVoronoiDiagram(
119 Delaunay* delaunay,
120 GEO::Mesh* mesh
121 ) :
122 172 mesh_(mesh),
123
2/2
✓ Branch 0 taken 46 times.
✓ Branch 1 taken 40 times.
 172 delaunay_(delaunay),
124 intersections_(DIM),
125 172 symbolic_(false),
126 172 check_SR_(true),
127 172 exact_(false)
128 {
129
2/2
✓ Branch 0 taken 46 times.
✓ Branch 1 taken 40 times.
 172 delaunay_nn_ = dynamic_cast<GEO::Delaunay_NearestNeighbors*>(
130 delaunay_
131 );
132 172 dimension_ = DIM;
133 172 facets_begin_ = UNSPECIFIED_RANGE;
134 172 facets_end_ = UNSPECIFIED_RANGE;
135 172 tets_begin_ = UNSPECIFIED_RANGE;
136 172 tets_end_ = UNSPECIFIED_RANGE;
137 172 connected_components_priority_ = false;
138 172 facet_seed_marking_ = nullptr;
139 172 connected_component_changed_ = false;
140 172 current_connected_component_ = 0;
141 172 cur_stamp_ = NO_INDEX;
142 172 current_facet_ = GEO::max_index_t();
143 172 current_seed_ = GEO::max_index_t();
144 172 current_polygon_ = nullptr;
145 172 current_tet_ = GEO::max_index_t();
146 172 current_polyhedron_ = nullptr;
147 172 }
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 void set_connected_components_priority(bool x) {
157 ✗ connected_components_priority_ = x;
158 6 }
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 GEO::Mesh* mesh() {
186
4/14
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
 4 return mesh_;
187 }
188
189 /**
190 * \brief Gets the Delaunay triangulation.
191 */
192 const Delaunay* delaunay() const {
193
1/14
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✓ Branch 54 taken 4 times.
✗ Branch 55 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 62 not taken.
✗ Branch 63 not taken.
✗ Branch 66 not taken.
✗ Branch 67 not taken.
 11893313 return delaunay_;
194 }
195
196 /**
197 * \brief Gets the Delaunay triangulation.
198 */
199 Delaunay* delaunay() {
200 80 return delaunay_;
201 }
202
203 /**
204 * \brief Sets the Delaunay triangulation.
205 */
206 void set_delaunay(Delaunay* delaunay) {
207 40 delaunay_ = delaunay;
208
4/14
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 20 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 4 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 40 delaunay_nn_ = dynamic_cast<GEO::Delaunay_NearestNeighbors*>(
209 delaunay_
210 );
211 }
212
213 /**
214 * \brief Sets the input mesh.
215 */
216 void set_mesh(GEO::Mesh* mesh) {
217
0/14
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
 40 mesh_ = mesh;
218 }
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 void set_facets_range(index_t facets_begin, index_t facets_end) {
227 geo_debug_assert(facets_end >= facets_begin);
228 40 facets_begin_ = facets_begin;
229 40 facets_end_ = facets_end;
230 }
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 void set_tetrahedra_range(index_t tets_begin, index_t tets_end) {
241 geo_debug_assert(tets_end >= tets_begin);
242 16 tets_begin_ = tets_begin;
243 16 tets_end_ = tets_end;
244 }
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 index_t current_facet() const {
268
2/42
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 60502 times.
✓ Branch 3 taken 149907 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
 210409 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 const Polygon& current_polygon() const {
288
7/364
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✓ Branch 57 taken 1051420 times.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✗ Branch 63 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 66 not taken.
✓ Branch 67 taken 112972 times.
✗ Branch 68 not taken.
✗ Branch 69 not taken.
✗ Branch 70 not taken.
✗ Branch 71 not taken.
✗ Branch 72 not taken.
✗ Branch 73 not taken.
✗ Branch 74 not taken.
✗ Branch 75 not taken.
✗ Branch 76 not taken.
✗ Branch 77 not taken.
✗ Branch 78 not taken.
✗ Branch 79 not taken.
✗ Branch 80 not taken.
✓ Branch 81 taken 3740244 times.
✗ Branch 82 not taken.
✗ Branch 83 not taken.
✗ Branch 84 not taken.
✗ Branch 85 not taken.
✗ Branch 86 not taken.
✗ Branch 87 not taken.
✗ Branch 88 not taken.
✗ Branch 89 not taken.
✗ Branch 90 not taken.
✗ Branch 91 not taken.
✗ Branch 92 not taken.
✗ Branch 93 not taken.
✗ Branch 94 not taken.
✗ Branch 95 not taken.
✗ Branch 96 not taken.
✓ Branch 97 taken 573384 times.
✗ Branch 98 not taken.
✗ Branch 99 not taken.
✗ Branch 100 not taken.
✗ Branch 101 not taken.
✗ Branch 102 not taken.
✗ Branch 103 not taken.
✗ Branch 104 not taken.
✗ Branch 105 not taken.
✗ Branch 106 not taken.
✗ Branch 107 not taken.
✗ Branch 108 not taken.
✗ Branch 109 not taken.
✗ Branch 110 not taken.
✗ Branch 111 not taken.
✗ Branch 112 not taken.
✗ Branch 113 not taken.
✗ Branch 114 not taken.
✗ Branch 115 not taken.
✗ Branch 116 not taken.
✗ Branch 117 not taken.
✗ Branch 118 not taken.
✗ Branch 119 not taken.
✗ Branch 120 not taken.
✗ Branch 121 not taken.
✗ Branch 122 not taken.
✗ Branch 123 not taken.
✗ Branch 124 not taken.
✗ Branch 125 not taken.
✗ Branch 126 not taken.
✗ Branch 127 not taken.
✗ Branch 128 not taken.
✗ Branch 129 not taken.
✗ Branch 130 not taken.
✗ Branch 131 not taken.
✗ Branch 132 not taken.
✗ Branch 133 not taken.
✗ Branch 134 not taken.
✗ Branch 135 not taken.
✗ Branch 136 not taken.
✗ Branch 137 not taken.
✗ Branch 138 not taken.
✗ Branch 139 not taken.
✗ Branch 140 not taken.
✗ Branch 141 not taken.
✗ Branch 142 not taken.
✗ Branch 143 not taken.
✗ Branch 144 not taken.
✗ Branch 145 not taken.
✗ Branch 146 not taken.
✗ Branch 147 not taken.
✗ Branch 148 not taken.
✗ Branch 149 not taken.
✗ Branch 150 not taken.
✗ Branch 151 not taken.
✗ Branch 152 not taken.
✗ Branch 153 not taken.
✗ Branch 154 not taken.
✗ Branch 155 not taken.
✗ Branch 156 not taken.
✗ Branch 157 not taken.
✗ Branch 158 not taken.
✗ Branch 159 not taken.
✗ Branch 160 not taken.
✗ Branch 161 not taken.
✗ Branch 162 not taken.
✗ Branch 163 not taken.
✗ Branch 164 not taken.
✗ Branch 165 not taken.
✗ Branch 166 not taken.
✗ Branch 167 not taken.
✗ Branch 168 not taken.
✗ Branch 169 not taken.
✗ Branch 170 not taken.
✓ Branch 171 taken 312339 times.
✗ Branch 172 not taken.
✗ Branch 173 not taken.
✗ Branch 174 not taken.
✗ Branch 175 not taken.
✗ Branch 176 not taken.
✗ Branch 177 not taken.
✗ Branch 178 not taken.
✗ Branch 179 not taken.
✗ Branch 180 not taken.
✗ Branch 181 not taken.
✗ Branch 182 not taken.
✗ Branch 183 not taken.
✗ Branch 184 not taken.
✓ Branch 185 taken 10382596 times.
✗ Branch 186 not taken.
✗ Branch 187 not taken.
✗ Branch 188 not taken.
✗ Branch 189 not taken.
✗ Branch 190 not taken.
✗ Branch 191 not taken.
✗ Branch 192 not taken.
✗ Branch 193 not taken.
✗ Branch 194 not taken.
✗ Branch 195 not taken.
✗ Branch 196 not taken.
✗ Branch 197 not taken.
✗ Branch 198 not taken.
✗ Branch 199 not taken.
✗ Branch 200 not taken.
✓ Branch 201 taken 1591980 times.
✗ Branch 202 not taken.
✗ Branch 203 not taken.
✗ Branch 204 not taken.
✗ Branch 205 not taken.
✗ Branch 206 not taken.
✗ Branch 207 not taken.
✗ Branch 208 not taken.
✗ Branch 209 not taken.
✗ Branch 210 not taken.
✗ Branch 211 not taken.
✗ Branch 212 not taken.
✗ Branch 213 not taken.
✗ Branch 214 not taken.
✗ Branch 215 not taken.
✗ Branch 216 not taken.
✗ Branch 217 not taken.
✗ Branch 218 not taken.
✗ Branch 219 not taken.
✗ Branch 220 not taken.
✗ Branch 221 not taken.
✗ Branch 222 not taken.
✗ Branch 223 not taken.
✗ Branch 224 not taken.
✗ Branch 225 not taken.
✗ Branch 226 not taken.
✗ Branch 227 not taken.
✗ Branch 228 not taken.
✗ Branch 229 not taken.
✗ Branch 230 not taken.
✗ Branch 231 not taken.
✗ Branch 232 not taken.
✗ Branch 233 not taken.
✗ Branch 234 not taken.
✗ Branch 235 not taken.
✗ Branch 236 not taken.
✗ Branch 237 not taken.
✗ Branch 238 not taken.
✗ Branch 239 not taken.
✗ Branch 240 not taken.
✗ Branch 241 not taken.
✗ Branch 242 not taken.
✗ Branch 243 not taken.
✗ Branch 244 not taken.
✗ Branch 245 not taken.
✗ Branch 246 not taken.
✗ Branch 247 not taken.
✗ Branch 248 not taken.
✗ Branch 249 not taken.
✗ Branch 250 not taken.
✗ Branch 251 not taken.
✗ Branch 252 not taken.
✗ Branch 253 not taken.
✗ Branch 254 not taken.
✗ Branch 255 not taken.
✗ Branch 256 not taken.
✗ Branch 257 not taken.
✗ Branch 258 not taken.
✗ Branch 259 not taken.
✗ Branch 260 not taken.
✗ Branch 261 not taken.
✗ Branch 262 not taken.
✗ Branch 263 not taken.
✗ Branch 264 not taken.
✗ Branch 265 not taken.
✗ Branch 266 not taken.
✗ Branch 267 not taken.
✗ Branch 268 not taken.
✗ Branch 269 not taken.
✗ Branch 270 not taken.
✗ Branch 271 not taken.
✗ Branch 272 not taken.
✗ Branch 273 not taken.
✗ Branch 274 not taken.
✗ Branch 275 not taken.
✗ Branch 276 not taken.
✗ Branch 277 not taken.
✗ Branch 278 not taken.
✗ Branch 279 not taken.
✗ Branch 280 not taken.
✗ Branch 281 not taken.
✗ Branch 282 not taken.
✗ Branch 283 not taken.
✗ Branch 284 not taken.
✗ Branch 285 not taken.
✗ Branch 286 not taken.
✗ Branch 287 not taken.
✗ Branch 288 not taken.
✗ Branch 289 not taken.
✗ Branch 290 not taken.
✗ Branch 291 not taken.
✗ Branch 292 not taken.
✗ Branch 293 not taken.
✗ Branch 294 not taken.
✗ Branch 295 not taken.
✗ Branch 296 not taken.
✗ Branch 297 not taken.
✗ Branch 298 not taken.
✗ Branch 299 not taken.
✗ Branch 300 not taken.
✗ Branch 301 not taken.
✗ Branch 302 not taken.
✗ Branch 303 not taken.
✗ Branch 304 not taken.
✗ Branch 305 not taken.
✗ Branch 306 not taken.
✗ Branch 307 not taken.
✗ Branch 308 not taken.
✗ Branch 309 not taken.
✗ Branch 310 not taken.
✗ Branch 311 not taken.
✗ Branch 312 not taken.
✗ Branch 313 not taken.
✗ Branch 314 not taken.
✗ Branch 315 not taken.
✗ Branch 316 not taken.
✗ Branch 317 not taken.
✗ Branch 318 not taken.
✗ Branch 319 not taken.
✗ Branch 320 not taken.
✗ Branch 321 not taken.
✗ Branch 322 not taken.
✗ Branch 323 not taken.
✗ Branch 324 not taken.
✗ Branch 325 not taken.
✗ Branch 326 not taken.
✗ Branch 327 not taken.
✗ Branch 328 not taken.
✗ Branch 329 not taken.
✗ Branch 330 not taken.
✗ Branch 331 not taken.
✗ Branch 332 not taken.
✗ Branch 333 not taken.
✗ Branch 334 not taken.
✗ Branch 335 not taken.
✗ Branch 336 not taken.
✗ Branch 337 not taken.
✗ Branch 338 not taken.
✗ Branch 339 not taken.
✗ Branch 340 not taken.
✗ Branch 341 not taken.
✗ Branch 342 not taken.
✗ Branch 343 not taken.
✗ Branch 344 not taken.
✗ Branch 345 not taken.
✗ Branch 346 not taken.
✗ Branch 347 not taken.
✗ Branch 348 not taken.
✗ Branch 349 not taken.
✗ Branch 350 not taken.
✗ Branch 351 not taken.
✗ Branch 352 not taken.
✗ Branch 353 not taken.
✗ Branch 354 not taken.
✗ Branch 355 not taken.
✗ Branch 356 not taken.
✗ Branch 357 not taken.
✗ Branch 358 not taken.
✗ Branch 359 not taken.
✗ Branch 360 not taken.
✗ Branch 361 not taken.
✗ Branch 362 not taken.
✗ Branch 363 not taken.
 17764935 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 const Polyhedron& current_polyhedron() const {
309
12/252
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✓ Branch 39 taken 57806 times.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 57 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✓ Branch 61 taken 2048838 times.
✗ Branch 62 not taken.
✗ Branch 63 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 66 not taken.
✗ Branch 67 not taken.
✗ Branch 68 not taken.
✓ Branch 69 taken 356901 times.
✗ Branch 70 not taken.
✗ Branch 71 not taken.
✗ Branch 72 not taken.
✗ Branch 73 not taken.
✗ Branch 74 not taken.
✓ Branch 75 taken 75602 times.
✗ Branch 76 not taken.
✗ Branch 77 not taken.
✗ Branch 78 not taken.
✗ Branch 79 not taken.
✗ Branch 80 not taken.
✗ Branch 81 not taken.
✗ Branch 82 not taken.
✗ Branch 83 not taken.
✗ Branch 84 not taken.
✗ Branch 85 not taken.
✗ Branch 86 not taken.
✗ Branch 87 not taken.
✗ Branch 88 not taken.
✗ Branch 89 not taken.
✗ Branch 90 not taken.
✗ Branch 91 not taken.
✗ Branch 92 not taken.
✗ Branch 93 not taken.
✗ Branch 94 not taken.
✗ Branch 95 not taken.
✗ Branch 96 not taken.
✓ Branch 97 taken 3023729 times.
✗ Branch 98 not taken.
✗ Branch 99 not taken.
✗ Branch 100 not taken.
✗ Branch 101 not taken.
✗ Branch 102 not taken.
✗ Branch 103 not taken.
✗ Branch 104 not taken.
✓ Branch 105 taken 460817 times.
✗ Branch 106 not taken.
✗ Branch 107 not taken.
✗ Branch 108 not taken.
✗ Branch 109 not taken.
✗ Branch 110 not taken.
✓ Branch 111 taken 99143 times.
✗ Branch 112 not taken.
✗ Branch 113 not taken.
✗ Branch 114 not taken.
✗ Branch 115 not taken.
✗ Branch 116 not taken.
✗ Branch 117 not taken.
✗ Branch 118 not taken.
✗ Branch 119 not taken.
✗ Branch 120 not taken.
✗ Branch 121 not taken.
✗ Branch 122 not taken.
✗ Branch 123 not taken.
✗ Branch 124 not taken.
✗ Branch 125 not taken.
✗ Branch 126 not taken.
✗ Branch 127 not taken.
✗ Branch 128 not taken.
✗ Branch 129 not taken.
✗ Branch 130 not taken.
✗ Branch 131 not taken.
✗ Branch 132 not taken.
✓ Branch 133 taken 3350910 times.
✗ Branch 134 not taken.
✗ Branch 135 not taken.
✗ Branch 136 not taken.
✗ Branch 137 not taken.
✗ Branch 138 not taken.
✗ Branch 139 not taken.
✗ Branch 140 not taken.
✓ Branch 141 taken 544856 times.
✗ Branch 142 not taken.
✗ Branch 143 not taken.
✗ Branch 144 not taken.
✗ Branch 145 not taken.
✗ Branch 146 not taken.
✓ Branch 147 taken 99838 times.
✗ Branch 148 not taken.
✗ Branch 149 not taken.
✗ Branch 150 not taken.
✗ Branch 151 not taken.
✗ Branch 152 not taken.
✗ Branch 153 not taken.
✗ Branch 154 not taken.
✗ Branch 155 not taken.
✗ Branch 156 not taken.
✗ Branch 157 not taken.
✗ Branch 158 not taken.
✗ Branch 159 not taken.
✗ Branch 160 not taken.
✗ Branch 161 not taken.
✗ Branch 162 not taken.
✗ Branch 163 not taken.
✗ Branch 164 not taken.
✗ Branch 165 not taken.
✗ Branch 166 not taken.
✗ Branch 167 not taken.
✗ Branch 168 not taken.
✓ Branch 169 taken 3418339 times.
✗ Branch 170 not taken.
✗ Branch 171 not taken.
✗ Branch 172 not taken.
✗ Branch 173 not taken.
✗ Branch 174 not taken.
✗ Branch 175 not taken.
✗ Branch 176 not taken.
✓ Branch 177 taken 582979 times.
✗ Branch 178 not taken.
✗ Branch 179 not taken.
✗ Branch 180 not taken.
✗ Branch 181 not taken.
✗ Branch 182 not taken.
✗ Branch 183 not taken.
✗ Branch 184 not taken.
✗ Branch 185 not taken.
✗ Branch 186 not taken.
✗ Branch 187 not taken.
✗ Branch 188 not taken.
✗ Branch 189 not taken.
✗ Branch 190 not taken.
✗ Branch 191 not taken.
✗ Branch 192 not taken.
✗ Branch 193 not taken.
✗ Branch 194 not taken.
✗ Branch 195 not taken.
✗ Branch 196 not taken.
✗ Branch 197 not taken.
✗ Branch 198 not taken.
✗ Branch 199 not taken.
✗ Branch 200 not taken.
✗ Branch 201 not taken.
✗ Branch 202 not taken.
✗ Branch 203 not taken.
✗ Branch 204 not taken.
✗ Branch 205 not taken.
✗ Branch 206 not taken.
✗ Branch 207 not taken.
✗ Branch 208 not taken.
✗ Branch 209 not taken.
✗ Branch 210 not taken.
✗ Branch 211 not taken.
✗ Branch 212 not taken.
✗ Branch 213 not taken.
✗ Branch 214 not taken.
✗ Branch 215 not taken.
✗ Branch 216 not taken.
✗ Branch 217 not taken.
✗ Branch 218 not taken.
✗ Branch 219 not taken.
✗ Branch 220 not taken.
✗ Branch 221 not taken.
✗ Branch 222 not taken.
✗ Branch 223 not taken.
✗ Branch 224 not taken.
✗ Branch 225 not taken.
✗ Branch 226 not taken.
✗ Branch 227 not taken.
✗ Branch 228 not taken.
✗ Branch 229 not taken.
✗ Branch 230 not taken.
✗ Branch 231 not taken.
✗ Branch 232 not taken.
✗ Branch 233 not taken.
✗ Branch 234 not taken.
✗ Branch 235 not taken.
✗ Branch 236 not taken.
✗ Branch 237 not taken.
✗ Branch 238 not taken.
✗ Branch 239 not taken.
✗ Branch 240 not taken.
✗ Branch 241 not taken.
✗ Branch 242 not taken.
✗ Branch 243 not taken.
✗ Branch 244 not taken.
✗ Branch 245 not taken.
✗ Branch 246 not taken.
✗ Branch 247 not taken.
✗ Branch 248 not taken.
✗ Branch 249 not taken.
✗ Branch 250 not taken.
✗ Branch 251 not taken.
 14119758 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 void set_symbolic(bool x) {
319 92 symbolic_ = x;
320 // exact mode requires symbolic mode.
321
8/70
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 35 times.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✓ Branch 32 taken 1 times.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✓ Branch 37 taken 4 times.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✓ Branch 42 taken 1 times.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 57 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✗ Branch 63 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 66 not taken.
✗ Branch 67 not taken.
✗ Branch 68 not taken.
✗ Branch 69 not taken.
 46 if(exact_) {
322 43 symbolic_ = true;
323 }
324 }
325
326 /**
327 * \brief Tests whether symbolic mode is active.
328 */
329 bool symbolic() const {
330
5/98
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✓ Branch 28 taken 36 times.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✓ Branch 38 taken 1 times.
✗ Branch 39 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✓ Branch 54 taken 1 times.
✗ Branch 55 not taken.
✗ Branch 57 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✗ Branch 62 not taken.
✗ Branch 63 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 66 not taken.
✗ Branch 67 not taken.
✗ Branch 68 not taken.
✗ Branch 69 not taken.
✓ Branch 70 taken 1 times.
✗ Branch 71 not taken.
✗ Branch 73 not taken.
✗ Branch 74 not taken.
✗ Branch 75 not taken.
✗ Branch 76 not taken.
✗ Branch 78 not taken.
✗ Branch 79 not taken.
✗ Branch 80 not taken.
✗ Branch 81 not taken.
✗ Branch 82 not taken.
✗ Branch 83 not taken.
✗ Branch 84 not taken.
✗ Branch 85 not taken.
✗ Branch 86 not taken.
✗ Branch 87 not taken.
✗ Branch 89 not taken.
✗ Branch 90 not taken.
✗ Branch 91 not taken.
✗ Branch 92 not taken.
✗ Branch 94 not taken.
✗ Branch 95 not taken.
✗ Branch 96 not taken.
✗ Branch 97 not taken.
✗ Branch 98 not taken.
✗ Branch 99 not taken.
✗ Branch 100 not taken.
✗ Branch 101 not taken.
✗ Branch 102 not taken.
✗ Branch 103 not taken.
✗ Branch 105 not taken.
✗ Branch 106 not taken.
✗ Branch 107 not taken.
✗ Branch 108 not taken.
✗ Branch 110 not taken.
✗ Branch 111 not taken.
 46 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 void set_exact_predicates(bool x) {
339 95 exact_ = x;
340 // exact mode requires symbolic mode.
341
8/14
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 40 times.
✓ Branch 3 taken 12 times.
✓ Branch 4 taken 5 times.
✓ Branch 5 taken 4 times.
✓ Branch 6 taken 5 times.
✓ Branch 7 taken 20 times.
✓ Branch 8 taken 5 times.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 95 if(exact_) {
342 55 symbolic_ = true;
343 }
344 }
345
346 /**
347 * \brief Tests whether exact predicates are used.
348 */
349 bool exact_predicates() const {
350 40 return exact_;
351 }
352
353 /**
354 * \brief Specifies whether radius of security should be enforced.
355 */
356 void set_check_SR(bool x) {
357 106 check_SR_ = x;
358 }
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 bool check_SR() const {
366 40 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 42 PolygonAction(const ACTION& do_it) :
399 42 do_it_(do_it) {
400 }
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 void operator() (
411 index_t v,
412 index_t f,
413 const Polygon& P
414 ) const {
415 GEO::geo_argused(f);
416 285070 const_cast<ACTION&> ( do_it_)(v, P);
417 285070 }
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 912 TriangleAction(const ACTION& do_it) :
439 912 do_it_(do_it) {
440 }
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 5725870 void operator() (
453 index_t v,
454 index_t f,
455 const Polygon& P
456 ) const {
457 GEO::geo_argused(f);
458
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10562750 times.
✓ Branch 2 taken 7699815 times.
✓ Branch 3 taken 2862935 times.
 21125500 for(index_t i = 1; i + 1 < P.nb_vertices(); i++) {
459 15399630 const_cast<ACTION&> (do_it_)(
460 v, P.vertex(0), P.vertex(i), P.vertex(i + 1)
461 );
462 }
463 5725870 }
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 4 PolyhedronAction(const ACTION& do_it) :
643 4 do_it_(do_it) {
644 }
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 void operator() (
655 index_t v,
656 index_t t,
657 const Polyhedron& C
658 ) const {
659 19667 const_cast<ACTION&> ( do_it_)(v, t, C);
660 19667 }
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 548 VolumetricIntegrationSimplexAction(
708 const ACTION& do_it,
709 bool visit_inner_tets = false,
710 bool coherent_triangles = false
711 ) :
712 548 do_it_(do_it),
713 548 visit_inner_tets_(visit_inner_tets),
714 548 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 1338616 void operator() (index_t v, index_t t, const Polyhedron& C) const {
727
2/2
✓ Branch 0 taken 11172508 times.
✓ Branch 1 taken 669308 times.
 47367264 for(index_t cv = 0; cv < C.max_v(); ++cv) {
728 signed_index_t ct = C.vertex_triangle(cv);
729
2/2
✓ Branch 0 taken 6848074 times.
✓ Branch 1 taken 4324434 times.
 22345016 if(ct == -1) {
730 17121288 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
2/2
✓ Branch 0 taken 1712570 times.
✓ Branch 1 taken 2611864 times.
 8648868 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
1/2
✓ Branch 0 taken 1712570 times.
✗ Branch 1 not taken.
 3425140 if(!visit_inner_tets_) {
744 3425140 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 5223728 v_adj = adjacent - 1;
751 }
752 // and adjacent indicex equal to zero corresponds
753 // to tet on border.
754
755
2/2
✓ Branch 0 taken 2177675 times.
✓ Branch 1 taken 434189 times.
 5223728 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
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2611864 times.
 5223728 if(coherent_triangles_) {
764 ✗ move_to_first_corner_of_facet(C, c1, v);
765 }
766
767 5223728 const Vertex& v1 = C.triangle_dual(c1.t);
768
769 5223728 Polyhedron::Corner c2 = c1;
770 5223728 C.move_to_next_around_vertex(c2);
771 geo_debug_assert(c2 != c1);
772
773 5223728 Polyhedron::Corner c3 = c2;
774 5223728 C.move_to_next_around_vertex(c3);
775 geo_debug_assert(c3 != c1);
776 do {
777 10369872 const Vertex& v2 = C.triangle_dual(c2.t);
778 10369872 const Vertex& v3 = C.triangle_dual(c3.t);
779 10369872 const_cast<ACTION&> (do_it_)(
780 v, index_t(v_adj), t, index_t(t_adj), v1, v2, v3
781 );
782 10369872 c2 = c3;
783 10369872 C.move_to_next_around_vertex(c3);
784
2/2
✓ Branch 0 taken 2573072 times.
✓ Branch 1 taken 2611864 times.
 10369872 } while(c3 != c1);
785 }
786 1338616 }
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 80 TetrahedronAction(
887 const ACTION& do_it
888 ) :
889 80 do_it_(do_it)
890 {
891 }
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 204510 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 const Vertex* v0 = nullptr;
911 index_t t0;
912
2/2
✓ Branch 0 taken 151921 times.
✓ Branch 1 taken 1365 times.
 613144 for(t0 = 0; t0 < C.max_t(); ++t0) {
913
2/2
✓ Branch 0 taken 100890 times.
✓ Branch 1 taken 51031 times.
 303842 if(C.triangle_is_used(t0)) {
914 v0 = &C.triangle_dual(t0);
915 201780 break;
916 }
917 }
918
919 // If current cell is empty, return
920
2/2
✓ Branch 0 taken 100890 times.
✓ Branch 1 taken 1365 times.
 204510 if(v0 == nullptr) {
921 return;
922 }
923
924
2/2
✓ Branch 0 taken 1843298 times.
✓ Branch 1 taken 100890 times.
 7574972 for(index_t cv = 0; cv < C.max_v(); ++cv) {
925 signed_index_t ct = C.vertex_triangle(cv);
926
2/2
✓ Branch 0 taken 1188589 times.
✓ Branch 1 taken 654709 times.
 3686596 if(ct == -1) {
927 2982518 continue;
928 }
929 geo_debug_assert(C.triangle_is_used(index_t(ct)));
930
931 signed_index_t adjacent = C.vertex_id(cv);
932 signed_index_t v_adj = -1;
933 signed_index_t t_adj = -1;
934
935
2/2
✓ Branch 0 taken 257959 times.
✓ Branch 1 taken 396750 times.
 1309418 if(adjacent < 0) {
936 // Negative adjacent indices correspond to
937 // tet-tet links
938 515918 t_adj = -adjacent - 1;
939 } else if(adjacent > 0) {
940 // Positive adjacent indices correspond to
941 // Voronoi seed - Voroni seed link
942 793500 v_adj = adjacent - 1;
943 }
944 // and adjacent indicex equal to zero corresponds
945 // to tet on border.
946
947
2/2
✓ Branch 0 taken 472813 times.
✓ Branch 1 taken 181896 times.
 1309418 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
2/2
✓ Branch 0 taken 302670 times.
✓ Branch 1 taken 352039 times.
 1309418 if(facet_is_incident_to_vertex(C, c1, t0)) {
955 605340 continue;
956 }
957
958 const Vertex& v1 = C.triangle_dual(c1.t);
959
960 704078 Polyhedron::Corner c2 = c1;
961 704078 C.move_to_next_around_vertex(c2);
962 geo_debug_assert(c2 != c1);
963
964 704078 Polyhedron::Corner c3 = c2;
965 704078 C.move_to_next_around_vertex(c3);
966 geo_debug_assert(c3 != c1);
967 do {
968 1506142 const Vertex& v2 = C.triangle_dual(c2.t);
969 1506142 const Vertex& v3 = C.triangle_dual(c3.t);
970 1506142 const_cast<ACTION&> (do_it_)(
971 v, index_t(v_adj), t, index_t(t_adj),
972 *v0, v1, v2, v3
973 );
974 1506142 c2 = c3;
975 1506142 C.move_to_next_around_vertex(c3);
976
2/2
✓ Branch 0 taken 401032 times.
✓ Branch 1 taken 352039 times.
 1506142 } 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 1309418 bool facet_is_incident_to_vertex(
992 const Polyhedron& C, Polyhedron::Corner& c, index_t t
993 ) const {
994 1309418 Polyhedron::Corner first = c;
995 1309418 Polyhedron::Corner cur = c;
996 do {
997
2/2
✓ Branch 0 taken 2088497 times.
✓ Branch 1 taken 302670 times.
 4782334 if(cur.t == t) {
998 return true;
999 }
1000 4176994 C.move_to_next_around_vertex(cur);
1001
2/2
✓ Branch 0 taken 1736458 times.
✓ Branch 1 taken 352039 times.
 4176994 } while(cur != first);
1002 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 inline void for_each_polygon(const ACTION& action) {
1081
3/28
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 11 taken 36 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✓ Branch 30 taken 4 times.
✗ Branch 31 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
 42 this->template compute_surfacic<PolygonAction<ACTION> >(
1082 42 PolygonAction<ACTION>(action)
1083 );
1084 }
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 inline void for_each_triangle(const TRIACTION& action) {
1096 912 this->template compute_surfacic<TriangleAction<TRIACTION> >(
1097 912 TriangleAction<TRIACTION>(action)
1098 );
1099 ✗ }
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 inline void for_each_polyhedron(const ACTION& action) {
1164 4 this->template compute_volumetric<PolyhedronAction<ACTION> >(
1165 4 PolyhedronAction<ACTION>(action)
1166 );
1167 4 }
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
0/28
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
 548 >(
1211 548 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 inline void for_each_tetrahedron(
1242 const ACTION& action
1243 ) {
1244
0/14
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
 80 this->template compute_volumetric<TetrahedronAction<ACTION> >(
1245 80 TetrahedronAction<ACTION>(
1246 action
1247 )
1248 );
1249 }
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 1908 inline void compute_surfacic(const ACTION& action) {
1295
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 948 times.
 1908 if(connected_components_priority_) {
1296 12 this->template compute_surfacic_with_cnx_priority<ACTION>(
1297 action
1298 );
1299 } else {
1300 1896 this->template compute_surfacic_with_seeds_priority<ACTION>(
1301 action
1302 );
1303 }
1304 1908 }
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 1896 inline void compute_surfacic_with_seeds_priority(const ACTION& action) {
1320 1896 if(
1321
2/2
✓ Branch 0 taken 36 times.
✓ Branch 1 taken 912 times.
 1896 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 1896 current_polygon_ = nullptr;
1328 GEO::vector<index_t> seed_stamp(
1329 1896 delaunay_->nb_vertices(), index_t(-1)
1330 );
1331
1/2
✓ Branch 1 taken 948 times.
✗ Branch 2 not taken.
 1896 GEO::vector<bool> facet_is_marked(facets_end_-facets_begin_, false);
1332
1/2
✓ Branch 1 taken 948 times.
✗ Branch 2 not taken.
 1896 init_get_neighbors();
1333
1334 FacetSeedStack adjacent_facets;
1335 SeedStack adjacent_seeds;
1336 Polygon F;
1337 GEO::Attribute<double> vertex_weight;
1338
1/2
✓ Branch 1 taken 948 times.
✗ Branch 2 not taken.
 1896 vertex_weight.bind_if_is_defined(
1339
2/4
✓ Branch 1 taken 948 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 948 times.
✗ Branch 5 not taken.
 1896 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 357826 times.
✓ Branch 1 taken 948 times.
 717548 for(index_t f = facets_begin_; f < facets_end_; f++) {
1347
2/2
✓ Branch 0 taken 1160 times.
✓ Branch 1 taken 356666 times.
 715652 if(!facet_is_marked[f-facets_begin_]) {
1348 // Propagate along the facet-graph.
1349 facet_is_marked[f-facets_begin_] = true;
1350 adjacent_facets.push(
1351
1/2
✓ Branch 1 taken 1160 times.
✗ Branch 2 not taken.
 2320 FacetSeed(f, find_seed_near_facet(f))
1352 );
1353
2/2
✓ Branch 0 taken 357826 times.
✓ Branch 1 taken 1160 times.
 717972 while(!adjacent_facets.empty()) {
1354 715652 current_facet_ = adjacent_facets.top().f;
1355
2/2
✓ Branch 0 taken 352940 times.
✓ Branch 1 taken 4886 times.
 715652 current_seed_ = adjacent_facets.top().seed;
1356 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 715652 F.initialize_from_mesh_facet(
1362
1/2
✓ Branch 1 taken 357826 times.
✗ Branch 2 not taken.
 715652 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 0 taken 357826 times.
✗ Branch 1 not taken.
 715652 seed_stamp[current_seed_] = current_facet_;
1370
1/2
✓ Branch 0 taken 357826 times.
✗ Branch 1 not taken.
 715652 adjacent_seeds.push(current_seed_);
1371
1372
2/2
✓ Branch 0 taken 3073344 times.
✓ Branch 1 taken 357826 times.
 6862340 while(!adjacent_seeds.empty()) {
1373
2/2
✓ Branch 0 taken 3071432 times.
✓ Branch 1 taken 1912 times.
 6146688 current_seed_ = adjacent_seeds.top();
1374 adjacent_seeds.pop();
1375
1376
1/2
✓ Branch 1 taken 3073344 times.
✗ Branch 2 not taken.
 6146688 current_polygon_ = intersect_cell_facet(
1377 current_seed_, F
1378 );
1379
1380
1/2
✓ Branch 1 taken 210409 times.
✗ Branch 2 not taken.
 6146688 action(
1381 current_seed_, current_facet_, current_polygon()
1382 );
1383
1384 // Propagate to adjacent facets and adjacent seeds
1385 34258430 for(index_t v = 0;
1386
2/2
✓ Branch 0 taken 14266280 times.
✓ Branch 1 taken 3073344 times.
 34679248 v < current_polygon().nb_vertices(); v++
1387 ) {
1388 const Vertex& ve = current_polygon().vertex(v);
1389 signed_index_t neigh_f = ve.adjacent_facet();
1390 28532560 if(
1391
2/2
✓ Branch 0 taken 3981660 times.
✓ Branch 1 taken 10284620 times.
 28532560 neigh_f >= signed_index_t(facets_begin_) &&
1392
2/2
✓ Branch 0 taken 3811923 times.
✓ Branch 1 taken 169737 times.
 7963320 neigh_f < signed_index_t(facets_end_) &&
1393
1/2
✓ Branch 0 taken 3811923 times.
✗ Branch 1 not taken.
 7623846 neigh_f != signed_index_t(current_facet_)
1394 ) {
1395 7623846 if(!facet_is_marked[
1396
2/2
✓ Branch 0 taken 356666 times.
✓ Branch 1 taken 3455257 times.
 7623846 index_t(neigh_f)-facets_begin_
1397 ]) {
1398 facet_is_marked[
1399 index_t(neigh_f)-facets_begin_
1400 ] = true;
1401 adjacent_facets.push(
1402
1/2
✓ Branch 1 taken 356666 times.
✗ Branch 2 not taken.
 1426664 FacetSeed(
1403 index_t(neigh_f),
1404 current_seed_
1405 )
1406 );
1407 }
1408 }
1409 signed_index_t neigh_s = ve.adjacent_seed();
1410
2/2
✓ Branch 0 taken 10096707 times.
✓ Branch 1 taken 4169573 times.
 28532560 if(neigh_s != -1) {
1411 20193414 if(
1412
2/2
✓ Branch 0 taken 2715518 times.
✓ Branch 1 taken 7381189 times.
 20193414 seed_stamp[neigh_s] != current_facet_
1413 ) {
1414 5431036 seed_stamp[neigh_s] = current_facet_;
1415
1/4
✓ Branch 1 taken 2715518 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 10862072 adjacent_seeds.push(index_t(neigh_s));
1416 }
1417 }
1418 }
1419 }
1420 }
1421 }
1422 }
1423
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 948 times.
 1896 current_polygon_ = nullptr;
1424 1896 }
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 1264 inline void compute_volumetric(const ACTION& action) {
1441
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 628 times.
 1264 if(connected_components_priority_) {
1442 8 this->template compute_volumetric_with_cnx_priority<ACTION>(
1443 action
1444 );
1445 } else {
1446 1256 this->template compute_volumetric_with_seeds_priority<ACTION>(
1447 action
1448 );
1449 }
1450 1264 }
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 1256 inline void compute_volumetric_with_seeds_priority(const ACTION& action){
1466 1256 if(
1467
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 628 times.
 1256 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 628 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 1256 geo_assert(tets_begin_ != UNSPECIFIED_RANGE);
1475
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 628 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 1256 geo_assert(tets_end_ != UNSPECIFIED_RANGE);
1476
1477 GEO::vector<index_t> seed_stamp(
1478 1256 delaunay_->nb_vertices(), index_t(-1)
1479 );
1480
1/2
✓ Branch 1 taken 628 times.
✗ Branch 2 not taken.
 1256 GEO::vector<bool> tet_is_marked(tets_end_-tets_begin_, false);
1481
1/2
✓ Branch 1 taken 628 times.
✗ Branch 2 not taken.
 1256 init_get_neighbors();
1482
1483 TetSeedStack adjacent_tets;
1484 SeedStack adjacent_seeds;
1485 Polyhedron C(dimension());
1486 GEO::Attribute<double> vertex_weight;
1487
1/2
✓ Branch 1 taken 628 times.
✗ Branch 2 not taken.
 1256 vertex_weight.bind_if_is_defined(
1488
2/4
✓ Branch 1 taken 628 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 628 times.
✗ Branch 5 not taken.
 1256 mesh_->vertices.attributes(), "weight"
1489 );
1490
1491 1256 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 60288 times.
✓ Branch 1 taken 628 times.
 121832 for(index_t t = tets_begin_; t < tets_end_; ++t) {
1497
2/2
✓ Branch 0 taken 628 times.
✓ Branch 1 taken 59660 times.
 120576 if(!tet_is_marked[t-tets_begin_]) {
1498 // Propagate along the tet-graph.
1499 tet_is_marked[t-tets_begin_] = true;
1500 adjacent_tets.push(
1501
1/2
✓ Branch 1 taken 628 times.
✗ Branch 2 not taken.
 1256 TetSeed(t, find_seed_near_tet(t))
1502 );
1503
2/2
✓ Branch 0 taken 60288 times.
✓ Branch 1 taken 628 times.
 121832 while(!adjacent_tets.empty()) {
1504 120576 current_tet_ = adjacent_tets.top().f;
1505
1/2
✓ Branch 0 taken 60288 times.
✗ Branch 1 not taken.
 120576 current_seed_ = adjacent_tets.top().seed;
1506 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 0 taken 60288 times.
✗ Branch 1 not taken.
 120576 seed_stamp[current_seed_] = current_tet_;
1519
1/2
✓ Branch 0 taken 60288 times.
✗ Branch 1 not taken.
 120576 adjacent_seeds.push(current_seed_);
1520
1521
2/2
✓ Branch 0 taken 771563 times.
✓ Branch 1 taken 60288 times.
 1663702 while(!adjacent_seeds.empty()) {
1522
1/2
✓ Branch 0 taken 771563 times.
✗ Branch 1 not taken.
 1543126 current_seed_ = adjacent_seeds.top();
1523 adjacent_seeds.pop();
1524
1525 1543126 C.initialize_from_mesh_tetrahedron(
1526
1/2
✓ Branch 1 taken 771563 times.
✗ Branch 2 not taken.
 1543126 mesh_, current_tet_, symbolic_, vertex_weight
1527 );
1528
1529
1/2
✓ Branch 1 taken 771563 times.
✗ Branch 2 not taken.
 1543126 intersect_cell_cell(
1530 current_seed_, C
1531 );
1532
1533
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 1543126 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 27574738 for(index_t v = 0;
1542
2/2
✓ Branch 0 taken 13015806 times.
✓ Branch 1 taken 771563 times.
 27574738 v < current_polyhedron().max_v(); ++v
1543 ) {
1544
1545 // Skip clipping planes that are no longer
1546 // connected to a cell facet.
1547
2/2
✓ Branch 0 taken 8036663 times.
✓ Branch 1 taken 4979143 times.
 26031612 if(
1548 current_polyhedron().vertex_triangle(v)
1549 == -1
1550 ) {
1551 16073326 continue;
1552 }
1553
1554 signed_index_t id =
1555 current_polyhedron().vertex_id(v);
1556
2/2
✓ Branch 0 taken 2816983 times.
✓ Branch 1 taken 2162160 times.
 9958286 if(id > 0) {
1557 // Propagate to adjacent seed
1558
2/2
✓ Branch 0 taken 711275 times.
✓ Branch 1 taken 2105708 times.
 5633966 index_t neigh_s = index_t(id - 1);
1559
2/2
✓ Branch 0 taken 711275 times.
✓ Branch 1 taken 2105708 times.
 5633966 if(seed_stamp[neigh_s] != current_tet_) {
1560
1/2
✓ Branch 0 taken 711275 times.
✗ Branch 1 not taken.
 1422550 seed_stamp[neigh_s] = current_tet_;
1561 adjacent_seeds.push(neigh_s);
1562 }
1563
2/2
✓ Branch 0 taken 1970529 times.
✓ Branch 1 taken 191631 times.
 4324320 } 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 3941058 signed_index_t neigh_t = -id - 1;
1570 3941058 if(
1571 neigh_t >=
1572
2/2
✓ Branch 0 taken 1901345 times.
✓ Branch 1 taken 69184 times.
 3941058 signed_index_t(tets_begin_) &&
1573
2/2
✓ Branch 0 taken 1832170 times.
✓ Branch 1 taken 69175 times.
 3802690 neigh_t < signed_index_t(tets_end_) &&
1574
1/2
✓ Branch 0 taken 1832170 times.
✗ Branch 1 not taken.
 3664340 neigh_t != signed_index_t(current_tet_)
1575 ) {
1576 3664340 if(!tet_is_marked[
1577
2/2
✓ Branch 0 taken 59660 times.
✓ Branch 1 taken 1772510 times.
 3664340 index_t(neigh_t)-tets_begin_
1578 ]) {
1579 tet_is_marked[
1580 index_t(neigh_t)-tets_begin_
1581 ] = true;
1582 adjacent_tets.push(
1583
1/4
✓ Branch 1 taken 59660 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 238640 TetSeed(
1584 index_t(neigh_t),
1585 current_seed_
1586 )
1587 );
1588 }
1589 }
1590 }
1591 }
1592 }
1593 }
1594 }
1595 }
1596
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 628 times.
 1256 current_polyhedron_ = nullptr;
1597 2512 }
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 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 8 geo_assert(tets_begin_ != UNSPECIFIED_RANGE);
1633
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 8 geo_assert(tets_end_ != UNSPECIFIED_RANGE);
1634
1635 8 current_polyhedron_ = nullptr;
1636 8 init_get_neighbors();
1637
1638 std::deque<TetSeed> adjacent_seeds;
1639 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
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 TetSeedMarking visited(
1647
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 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
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 facet_seed_marking_ = &visited;
1654 Polyhedron C(dimension());
1655 8 current_polyhedron_ = &C;
1656
1657
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 current_connected_component_ = 0;
1658 // index_t C_index = tets_end_ + 1; // Unused (see comment later)
1659
1660 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
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 8 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
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 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 adjacent_seeds.push_back(
1675
1/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 8 TetSeed(current_tet_, current_seed_)
1676 );
1677
1678 // Propagate along the Delaunay-graph.
1679
2/2
✓ Branch 0 taken 35662 times.
✓ Branch 1 taken 4 times.
 71332 while(!adjacent_seeds.empty()) {
1680 // Yes, f, because TetSeed is typedef-ed as FacetSeed
1681 71324 current_tet_ = adjacent_seeds.front().f;
1682 71324 current_seed_ = adjacent_seeds.front().seed;
1683 71324 adjacent_seeds.pop_front();
1684 91998 if(
1685
2/2
✓ Branch 0 taken 10337 times.
✓ Branch 1 taken 25325 times.
 71324 tet_stamp[current_tet_ - tets_begin_] ==
1686
2/2
✓ Branch 0 taken 10337 times.
✓ Branch 1 taken 25325 times.
 71324 current_seed_
1687 ) {
1688 20674 continue;
1689 }
1690
1691
2/2
✓ Branch 0 taken 24064 times.
✓ Branch 1 taken 1261 times.
 50650 if(visited.is_marked(current_tet_, current_seed_)) {
1692 48128 continue;
1693 }
1694
1695 2522 connected_component_changed_ = true;
1696
1/2
✓ Branch 0 taken 1261 times.
✗ Branch 1 not taken.
 2522 adjacent_tets.push(current_tet_);
1697 2522 tet_stamp[current_tet_ - tets_begin_] =
1698 2522 current_seed_;
1699
2/2
✓ Branch 0 taken 19667 times.
✓ Branch 1 taken 1261 times.
 44378 while(!adjacent_tets.empty()) {
1700
1/2
✓ Branch 0 taken 19667 times.
✗ Branch 1 not taken.
 39334 current_tet_ = adjacent_tets.top();
1701 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 39334 C.initialize_from_mesh_tetrahedron(
1709
1/2
✓ Branch 1 taken 19667 times.
✗ Branch 2 not taken.
 39334 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 19667 times.
✗ Branch 2 not taken.
 39334 intersect_cell_cell(current_seed_, C);
1722
1/2
✓ Branch 1 taken 19667 times.
✗ Branch 2 not taken.
 39334 action(
1723 current_seed_, current_tet_, current_polyhedron()
1724 );
1725 39334 connected_component_changed_ = false;
1726
1727 bool touches_RVC_border = false;
1728
1729 // Propagate to adjacent tets and adjacent seeds
1730 664778 for(index_t v = 0;
1731
2/2
✓ Branch 0 taken 312722 times.
✓ Branch 1 taken 19667 times.
 664778 v < current_polyhedron().max_v(); ++v
1732 ) {
1733
1734 // Skip clipping planes that are no longer
1735 // connected to a cell facet.
1736
2/2
✓ Branch 0 taken 186074 times.
✓ Branch 1 taken 126648 times.
 625444 if(
1737 current_polyhedron().vertex_triangle(v)
1738 == -1
1739 ) {
1740 372148 continue;
1741 }
1742
1743 signed_index_t id =
1744 current_polyhedron().vertex_id(v);
1745
2/2
✓ Branch 0 taken 50404 times.
✓ Branch 1 taken 76244 times.
 253296 if(id < 0) {
1746 // id == 0 corresponds to facet on boundary
1747 // (skipped)
1748 // id < 0 corresponds to adjacent tet index
1749 100808 signed_index_t s_neigh_t = -id-1;
1750 100808 if(
1751
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 50404 times.
 100808 s_neigh_t >= signed_index_t(tets_begin_)
1752 100808 &&
1753
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 50404 times.
 100808 s_neigh_t < signed_index_t(tets_end_)
1754 ) {
1755 geo_debug_assert(
1756 s_neigh_t !=
1757 signed_index_t(current_tet_)
1758 );
1759 100808 index_t neigh_t = index_t(s_neigh_t);
1760 100808 if(
1761
2/2
✓ Branch 0 taken 18406 times.
✓ Branch 1 taken 31998 times.
 100808 tet_stamp[neigh_t - tets_begin_] !=
1762
2/2
✓ Branch 0 taken 18406 times.
✓ Branch 1 taken 31998 times.
 100808 current_seed_
1763 ) {
1764
1/2
✓ Branch 0 taken 18406 times.
✗ Branch 1 not taken.
 36812 tet_stamp[neigh_t - tets_begin_] =
1765 current_seed_;
1766 adjacent_tets.push(neigh_t);
1767 }
1768 }
1769
2/2
✓ Branch 0 taken 4928 times.
✓ Branch 1 taken 71316 times.
 152488 } else if(id > 0) {
1770 142632 index_t neigh_s = index_t(id-1);
1771 touches_RVC_border = true;
1772 142632 TetSeed ts(current_tet_, neigh_s);
1773
2/2
✓ Branch 0 taken 35658 times.
✓ Branch 1 taken 35658 times.
 142632 if(!visited.is_marked(ts)) {
1774
1/2
✓ Branch 1 taken 35658 times.
✗ Branch 2 not taken.
 71316 adjacent_seeds.push_back(ts);
1775 }
1776 }
1777
1778 }
1779
1/2
✓ Branch 0 taken 19667 times.
✗ Branch 1 not taken.
 39334 if(touches_RVC_border) {
1780 39334 visited.mark(
1781 39334 TetSeed(current_tet_, current_seed_),
1782 current_connected_component_
1783 );
1784 }
1785 }
1786 2522 ++current_connected_component_;
1787 }
1788 }
1789 }
1790
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 8 facet_seed_marking_ = nullptr;
1791 16 }
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 index_t get_facet_seed_connected_component(index_t f, index_t s) const {
1814 geo_debug_assert(facet_seed_marking_ != nullptr);
1815 179412 return facet_seed_marking_->get_connected_component(
1816 FacetSeed(f, s)
1817 );
1818 }
1819
1820 /**
1821 * \brief Tests whether the current connected component changed.
1822 */
1823 bool connected_component_changed() const {
1824
4/14
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 10000 times.
✓ Branch 3 taken 8686 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 20000 times.
✓ Branch 7 taken 35975 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 74661 return connected_component_changed_;
1825 }
1826
1827 /**
1828 * \brief Gets the index of the current connected component.
1829 */
1830 index_t current_connected_component() const {
1831
4/14
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 29648 times.
✓ Branch 3 taken 30352 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 59354 times.
✓ Branch 7 taken 60058 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 179412 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 12 inline void compute_surfacic_with_cnx_priority(
1857 const ACTION& action
1858 ) {
1859
1860 12 if(
1861
1/2
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
 12 facets_begin_ == UNSPECIFIED_RANGE &&
1862
1/2
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
 12 facets_end_ == UNSPECIFIED_RANGE
1863 ) {
1864 12 facets_begin_ = 0;
1865 12 facets_end_ = mesh_->facets.nb();
1866 }
1867
1868 12 current_polygon_ = nullptr;
1869 12 init_get_neighbors();
1870
1871 std::deque<FacetSeed> adjacent_seeds;
1872 std::stack<index_t> adjacent_facets;
1873
1874 static constexpr index_t NO_STAMP = index_t(-1);
1875 12 GEO::vector<index_t> facet_stamp(
1876
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 12 facets_end_ - facets_begin_, NO_STAMP
1877 );
1878
1879
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 12 FacetSeedMarking visited(
1880
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 12 facets_end_ - facets_begin_, delaunay_->nb_vertices()
1881 );
1882
1883
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 12 facet_seed_marking_ = &visited;
1884 Polygon F;
1885 12 current_connected_component_ = 0;
1886
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 12 index_t F_index = facets_end_ + 1;
1887
1888 GEO::Attribute<double> vertex_weight;
1889
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 12 vertex_weight.bind_if_is_defined(
1890
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 12 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 7798 times.
✓ Branch 1 taken 6 times.
 15608 for(index_t f = facets_begin_; f < facets_end_; ++f) {
1898
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 7792 times.
 15596 if(facet_stamp[f - facets_begin_] == NO_STAMP) {
1899 12 current_facet_ = f;
1900
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 12 current_seed_ = find_seed_near_facet(f);
1901
1902 adjacent_seeds.push_back(
1903
1/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 12 FacetSeed(current_facet_, current_seed_)
1904 );
1905
1906 // Propagate along the Delaunay-graph.
1907 12 while(
1908
2/2
✓ Branch 0 taken 126673 times.
✓ Branch 1 taken 6 times.
 253358 !adjacent_seeds.empty()
1909 ) {
1910 253346 current_facet_ = adjacent_seeds.front().f;
1911 253346 current_seed_ = adjacent_seeds.front().seed;
1912 253346 adjacent_seeds.pop_front();
1913 330332 if(
1914
2/2
✓ Branch 0 taken 38493 times.
✓ Branch 1 taken 88180 times.
 253346 facet_stamp[current_facet_ - facets_begin_] ==
1915
2/2
✓ Branch 0 taken 38493 times.
✓ Branch 1 taken 88180 times.
 253346 current_seed_
1916 ) {
1917 76986 continue;
1918 }
1919
1920
2/2
✓ Branch 0 taken 58180 times.
✓ Branch 1 taken 30000 times.
 176360 if(visited.is_marked(current_facet_, current_seed_)) {
1921 116360 continue;
1922 }
1923
1924 60000 connected_component_changed_ = true;
1925
1/2
✓ Branch 0 taken 30000 times.
✗ Branch 1 not taken.
 60000 adjacent_facets.push(current_facet_);
1926 60000 facet_stamp[current_facet_ - facets_begin_] =
1927 60000 current_seed_;
1928
2/2
✓ Branch 0 taken 74661 times.
✓ Branch 1 taken 30000 times.
 209322 while(!adjacent_facets.empty()) {
1929
1/2
✓ Branch 0 taken 74661 times.
✗ Branch 1 not taken.
 149322 current_facet_ = adjacent_facets.top();
1930 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 61451 times.
✓ Branch 1 taken 13210 times.
 149322 if(F_index != current_facet_) {
1936 122902 F.initialize_from_mesh_facet(
1937
1/2
✓ Branch 1 taken 61451 times.
✗ Branch 2 not taken.
 122902 mesh_, current_facet_, symbolic_,
1938 vertex_weight
1939 );
1940 122902 F_index = current_facet_;
1941 }
1942
1943
1/2
✓ Branch 1 taken 74661 times.
✗ Branch 2 not taken.
 149322 current_polygon_ = intersect_cell_facet(
1944 current_seed_, F
1945 );
1946
1/2
✓ Branch 1 taken 74661 times.
✗ Branch 2 not taken.
 149322 action(
1947 current_seed_, current_facet_, current_polygon()
1948 );
1949 149322 connected_component_changed_ = false;
1950
1951 bool touches_RVC_border = false;
1952
1953 // Propagate to adjacent facets and adjacent seeds
1954 850622 for(index_t v = 0;
1955
2/2
✓ Branch 0 taken 350650 times.
✓ Branch 1 taken 74661 times.
 850622 v < current_polygon().nb_vertices(); v++
1956 ) {
1957 const Vertex& ve = current_polygon().vertex(v);
1958 signed_index_t s_neigh_f = ve.adjacent_facet();
1959 701300 if(
1960
2/2
✓ Branch 0 taken 96944 times.
✓ Branch 1 taken 253706 times.
 701300 s_neigh_f >= signed_index_t(facets_begin_)
1961 193888 &&
1962
1/2
✓ Branch 0 taken 96944 times.
✗ Branch 1 not taken.
 193888 s_neigh_f < signed_index_t(facets_end_)
1963 ) {
1964 geo_debug_assert(
1965 s_neigh_f !=
1966 signed_index_t(current_facet_)
1967 );
1968 193888 index_t neigh_f = index_t(s_neigh_f);
1969 193888 if(
1970
2/2
✓ Branch 0 taken 44661 times.
✓ Branch 1 taken 52283 times.
 193888 facet_stamp[neigh_f - facets_begin_] !=
1971
2/2
✓ Branch 0 taken 44661 times.
✓ Branch 1 taken 52283 times.
 193888 current_seed_
1972 ) {
1973
1/2
✓ Branch 0 taken 44661 times.
✗ Branch 1 not taken.
 89322 facet_stamp[neigh_f - facets_begin_] =
1974 current_seed_;
1975 adjacent_facets.push(neigh_f);
1976 }
1977 }
1978 signed_index_t neigh_s = ve.adjacent_seed();
1979
2/2
✓ Branch 0 taken 253334 times.
✓ Branch 1 taken 97316 times.
 701300 if(neigh_s != -1) {
1980 touches_RVC_border = true;
1981 506668 FacetSeed fs(
1982 current_facet_, index_t(neigh_s)
1983 );
1984
2/2
✓ Branch 0 taken 126667 times.
✓ Branch 1 taken 126667 times.
 506668 if(!visited.is_marked(fs)) {
1985
1/2
✓ Branch 1 taken 126667 times.
✗ Branch 2 not taken.
 253334 adjacent_seeds.push_back(fs);
1986 }
1987 }
1988 }
1989
2/2
✓ Branch 0 taken 74647 times.
✓ Branch 1 taken 14 times.
 149322 if(touches_RVC_border) {
1990 149294 visited.mark(
1991 149294 FacetSeed(current_facet_, current_seed_),
1992 current_connected_component_
1993 );
1994 }
1995 }
1996 60000 ++current_connected_component_;
1997 }
1998 }
1999 }
2000
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
 12 facet_seed_marking_ = nullptr;
2001 24 }
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 2332 index_t find_seed_near_facet(index_t f) {
2011
1/2
✓ Branch 0 taken 1166 times.
✗ Branch 1 not taken.
 2332 const double* p = mesh_->vertices.point_ptr(
2012 mesh_->facets.vertex(f,0)
2013 );
2014 2332 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 1264 index_t find_seed_near_tet(index_t t) {
2025 1264 index_t v = mesh_->cells.tet_vertex(t, 0);
2026 const double* p = mesh_->vertices.point_ptr(v);
2027 1264 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 3596 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 1750 times.
✓ Branch 2 taken 48 times.
✗ Branch 3 not taken.
 3596 if(exact_ && delaunay_nn_ != nullptr) {
2046 // TODO: may need more than 10
2047 index_t neighbors[10];
2048 double neighbors_sq_dist[10];
2049 index_t nb = 10;
2050 96 if(delaunay_nn_->nb_vertices() < nb) {
2051 nb = delaunay_nn_->nb_vertices();
2052 }
2053 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 15 times.
✓ Branch 1 taken 48 times.
 126 if(neighbors_sq_dist[i] != min_d) {
2060 break;
2061 }
2062 30 if(neighbors[i] < nearest) {
2063 nearest = neighbors[i];
2064 }
2065 }
2066 return nearest;
2067 }
2068
2069 3500 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 void swap_polygons(Polygon*& ping, Polygon*& pong) {
2083
8/56
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 4031330 times.
✓ Branch 13 taken 2471209 times.
✓ Branch 14 taken 943413 times.
✓ Branch 15 taken 3087917 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✓ Branch 28 taken 8304740 times.
✓ Branch 29 taken 4664299 times.
✓ Branch 30 taken 2204489 times.
✓ Branch 31 taken 6100251 times.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
 19471578 if(ping != &P1 && ping != &P2) {
2084 // First clipping operation, ping points to F
2085 // (current facet copied)
2086 3147902 ping = &P2;
2087 3147902 pong = &P1;
2088 } else {
2089 std::swap(ping, pong);
2090 }
2091 }
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
1/2
✓ Branch 0 taken 3148005 times.
✗ Branch 1 not taken.
 6296010 Polygon* intersect_cell_facet(index_t seed, Polygon& F) {
2101 intersections_.clear();
2102
2103 // Initialize ping-pong pointers for Sutherland-Hodgman
2104 // re-entrant clipping and copy current facet into 'ping' buffer.
2105 6296010 Polygon* ping = &F;
2106 6296010 Polygon* pong = &P2;
2107
2108 // Clip current facet by current Voronoi cell (associated with seed)
2109
1/2
✓ Branch 0 taken 3148005 times.
✗ Branch 1 not taken.
 6296010 if(delaunay_nn_ != nullptr) {
2110 6296010 clip_by_cell_SR(seed, ping, pong); // "Security Radius" mode.
2111 } else {
2112 ✗ clip_by_cell(seed, ping, pong); // Standard mode.
2113 }
2114
2115 6296010 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 6296010 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 6296010 const double* geo_restrict pi = delaunay_->vertex_ptr(i);
2140 geo_assume_aligned(pi, geo_dim_alignment(DIM));
2141
2142 index_t jj = 0;
2143 index_t prev_nb_neighbors = 0;
2144 6296010 neighbors_.resize(0);
2145
2/2
✓ Branch 0 taken 3173729 times.
✓ Branch 1 taken 96 times.
 12643660 while(neighbors_.size() < delaunay_nn_->nb_vertices() - 1) {
2146
2147 6347458 delaunay_nn_->get_neighbors(i, neighbors_);
2148
1/2
✓ Branch 0 taken 3173729 times.
✗ Branch 1 not taken.
 6347458 if(neighbors_.size() == 0) {
2149 6295818 return;
2150 }
2151
1/2
✓ Branch 0 taken 3173729 times.
✗ Branch 1 not taken.
 6347458 if(prev_nb_neighbors == neighbors_.size()) {
2152 return;
2153 }
2154
2155
2/2
✓ Branch 0 taken 22614645 times.
✓ Branch 1 taken 30662 times.
 84233770 for(; jj < neighbors_.size(); jj++) {
2156 45229290 index_t j = neighbors_[jj];
2157 45229290 double R2 = 0.0;
2158
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 110666494 times.
✓ Branch 2 taken 88051849 times.
✓ Branch 3 taken 22614645 times.
 221332988 for(index_t k = 0; k < ping->nb_vertices(); k++) {
2159 geo_decl_aligned(double dik);
2160 const double* geo_restrict pk = ping->vertex(k).point();
2161 geo_assume_aligned(pk, geo_dim_alignment(DIM));
2162 dik = GEO::Geom::distance2(pi, pk, dimension());
2163 176103698 R2 = std::max(R2, dik);
2164 }
2165 geo_decl_aligned(double dij);
2166 45229290 const double* geo_restrict pj = delaunay_->vertex_ptr(j);
2167 geo_assume_aligned(pj, geo_dim_alignment(DIM));
2168 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 3143067 times.
✓ Branch 1 taken 19471578 times.
 45229290 if(dij > 4.1 * R2) {
2173 6286134 return;
2174 }
2175 38943156 clip_by_plane(*ping, *pong, i, j);
2176 swap_polygons(ping, pong);
2177 }
2178
2179
2/2
✓ Branch 0 taken 25820 times.
✓ Branch 1 taken 4842 times.
 61324 if(!check_SR_) {
2180 return;
2181 }
2182
2183 index_t nb_neighbors = neighbors_.size();
2184 prev_nb_neighbors = nb_neighbors;
2185
2186
2/2
✓ Branch 0 taken 25793 times.
✓ Branch 1 taken 27 times.
 51640 if(nb_neighbors > 8) {
2187 51586 nb_neighbors += nb_neighbors / 8;
2188 } else {
2189 54 nb_neighbors++;
2190 }
2191
2192 nb_neighbors = std::min(
2193 nb_neighbors,
2194 51640 delaunay_nn_->nb_vertices() - 1
2195 );
2196
2197 51640 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 void clip_by_plane(
2241 Polygon& ping, Polygon& pong,
2242 index_t i, index_t j
2243 ) {
2244 19471578 ping.clip_by_plane<DIM>(
2245 19471578 pong, intersections_, mesh_, delaunay_, i, j, exact_, symbolic_
2246 );
2247 }
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 1582460 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 791230 times.
✗ Branch 1 not taken.
 1582460 if(delaunay_nn_ != nullptr) {
2265 1582460 clip_by_cell_SR(seed, C); // "Security Radius" mode.
2266 } else {
2267 ✗ clip_by_cell(seed, C); // Standard mode.
2268 }
2269 1582460 }
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 1582460 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 1582460 const double* geo_restrict pi = delaunay_->vertex_ptr(seed);
2285 geo_assume_aligned(pi, geo_dim_alignment(DIM));
2286
2287 index_t jj = 0;
2288 index_t prev_nb_neighbors = 0;
2289 1582460 neighbors_.resize(0);
2290
2291
1/2
✓ Branch 0 taken 798390 times.
✗ Branch 1 not taken.
 3179240 while(neighbors_.size() < delaunay_nn_->nb_vertices() - 1) {
2292
2293 1596780 delaunay_nn_->get_neighbors(seed, neighbors_);
2294
1/2
✓ Branch 0 taken 798390 times.
✗ Branch 1 not taken.
 1596780 if(neighbors_.size() == 0) {
2295 1582460 return;
2296 }
2297
1/2
✓ Branch 0 taken 798390 times.
✗ Branch 1 not taken.
 1596780 if(prev_nb_neighbors == neighbors_.size()) {
2298 return;
2299 }
2300
2301
2/2
✓ Branch 0 taken 10947472 times.
✓ Branch 1 taken 33192 times.
 42325876 for(; jj < neighbors_.size(); jj++) {
2302 21894944 index_t j = neighbors_[jj];
2303 21894944 double R2 = 0.0;
2304
2/2
✓ Branch 0 taken 117449460 times.
✓ Branch 1 taken 10947472 times.
 513587728 for(index_t k = 0; k < C.max_t(); ++k) {
2305
2/2
✓ Branch 0 taken 33088866 times.
✓ Branch 1 taken 84360594 times.
 234898920 if(!C.triangle_is_used(k)) {
2306 continue;
2307 }
2308 geo_decl_aligned(double dik);
2309 const double* geo_restrict pk =
2310 C.triangle_dual(k).point();
2311 geo_assume_aligned(pk, geo_dim_alignment(DIM));
2312 dik = GEO::Geom::distance2(pi, pk, dimension());
2313 168721188 R2 = std::max(R2, dik);
2314 }
2315 geo_decl_aligned(double dij);
2316 21894944 const double* geo_restrict pj = delaunay_->vertex_ptr(j);
2317 geo_assume_aligned(pj, geo_dim_alignment(DIM));
2318 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 765198 times.
✓ Branch 1 taken 10182274 times.
 21894944 if(dij > 4.1 * R2) {
2323 1530396 return;
2324 }
2325 clip_by_plane(C, seed, j);
2326 }
2327
2328
2/2
✓ Branch 0 taken 7160 times.
✓ Branch 1 taken 26032 times.
 66384 if(!check_SR_) {
2329 return;
2330 }
2331
2332 index_t nb_neighbors = neighbors_.size();
2333 prev_nb_neighbors = nb_neighbors;
2334
2335
1/2
✓ Branch 0 taken 7160 times.
✗ Branch 1 not taken.
 14320 if(nb_neighbors > 8) {
2336 14320 nb_neighbors += nb_neighbors / 8;
2337 } else {
2338 ✗ nb_neighbors++;
2339 }
2340
2341 nb_neighbors = std::min(
2342 nb_neighbors,
2343 14320 delaunay_nn_->nb_vertices() - 1
2344 );
2345 14320 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 void clip_by_plane(Polyhedron& C, index_t i, index_t j) {
2376 10182274 C.clip_by_plane<DIM>(
2377 10182274 mesh_, delaunay_, i, j, exact_, symbolic_
2378 );
2379 }
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 3172 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
3/4
✓ Branch 0 taken 499 times.
✓ Branch 1 taken 1087 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 499 times.
 3172 if(delaunay_->dimension() == 3 && delaunay_->nb_cells() != 0) {
2400 ✗ cur_stamp_ = 0;
2401 ✗ stamp_.assign(delaunay_->nb_vertices(), NO_INDEX);
2402 }
2403 3172 }
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