GCC Code Coverage Report


Directory: ./
File: delaunay/periodic_delaunay_3d.h
Date: 2026-09-27 03:24:14
Exec Total Coverage
Lines: 40 43 93.0%
Functions: 11 12 91.7%
Branches: 9 24 37.5%

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 PERIODIC_DELAUNAY_TRIANGULATION_3D
41 #define PERIODIC_DELAUNAY_TRIANGULATION_3D
42
43 #include <geogram/basic/common.h>
44 #include <geogram/delaunay/delaunay.h>
45 #include <geogram/delaunay/periodic.h>
46 #include <geogram/voronoi/convex_cell.h>
47 #include <geogram/basic/process.h>
48 #include <geogram/basic/geometry.h>
49 #include <stack>
50
51 #include <geogram/delaunay/delaunay_sync.h>
52
53 namespace GEO {
54
55 class PeriodicDelaunay3dThread;
56
57
58 /**
59 * \brief Multithreaded implementation of Delaunay in 3d with
60 * optional periodic boundary conditions.
61 * \details Periodicity is taken into account by detecting the
62 * vertices with Voronoi cells that straddle the boundary and
63 * duplicating them as need be (with an additional propagation to
64 * have the correct neighborhoods).
65 * \see Delaunay3d, ParallelDelaunay3d
66 */
67 class GEOGRAM_API PeriodicDelaunay3d : public Delaunay, public Periodic {
68 public:
69
70 /**
71 * \brief Gathers some structures used by some algorithms, makes
72 * multithreading more efficient by avoiding dynamic reallocations.
73 * \details It is used to compute the set of tetrahedra incident to
74 * a vertex. It gathers a stack and the vector of incident tets
75 * obtained so far.
76 */
77 struct IncidentTetrahedra {
78 std::stack<index_t> S;
79 vector<index_t> incident_tets_set;
80
81 /**
82 * \brief Clears the set of incident tets.
83 */
84 442 void clear_incident_tets() {
85 442 incident_tets_set.resize(0);
86 442 }
87
88 /**
89 * \brief Inserts a tet into the set of incident tets.
90 * \param[in] t the tet to be inserted.
91 */
92 9762 void add_incident_tet(index_t t) {
93 9762 incident_tets_set.push_back(t);
94 9762 }
95
96 /**
97 * \brief Tests whether a tet belongs to the set of incident
98 * tets.
99 * \param[in] t the tet to be tested
100 * \retval true if the tet belongs to the set of incident tets
101 * \retval false otherwise
102 */
103 29286 bool has_incident_tet(index_t t) const {
104
2/2
✓ Branch 1 taken 653533 times.
✓ Branch 2 taken 9320 times.
662853 for(index_t i=0; i<incident_tets_set.size(); ++i) {
105
2/2
✓ Branch 1 taken 19966 times.
✓ Branch 2 taken 633567 times.
653533 if(incident_tets_set[i] == t) {
106 19966 return true;
107 }
108 }
109 9320 return false;
110 }
111
112 442 vector<index_t>::const_iterator begin() const {
113 442 return incident_tets_set.begin();
114 }
115
116 442 vector<index_t>::const_iterator end() const {
117 442 return incident_tets_set.end();
118 }
119 };
120
121 /**
122 * \brief Constructs a new PeriodicDelaunay3d.
123 * \param[in] periodic if true, constructs a periodic triangulation.
124 * \param[in] period the edge length of the periodic domain.
125 */
126 PeriodicDelaunay3d(bool periodic, double period=1.0);
127
128 /**
129 * \brief Constructs a new PeriodicDelaunay3d.
130 * \param[in] period the edge lengths along x,y,z in the periodic domain
131 * \param[in] periodic if true, constructs a periodic triangulation.
132 */
133 PeriodicDelaunay3d(const vec3& period, bool periodic = true);
134
135 /**
136 * \copydoc Delaunay::set_vertices()
137 * \note compute() needs to be called after.
138 */
139 void set_vertices(
140 index_t nb_vertices, const double* vertices
141 ) override;
142
143 /**
144 * \brief Sets the weights.
145 * \param[in] weights pointer to the array of
146 * weights. Size is the number of real vertices,
147 * i.e., the parameter nb_vertices passed to
148 * set_vertices().
149 * \note compute() needs to be called after.
150 */
151 void set_weights(const double* weights);
152
153 /**
154 * \brief Computes the Delaunay triangulation.
155 */
156 void compute();
157
158 /**
159 * \brief Tests whether this PeriodicDelaunay3d is in periodic mode
160 * \retval true if in periodic mode, that is, uses periodic boundary
161 * conditions
162 * \retval false otherwise
163 */
164 2 bool periodic() const {
165 2 return periodic_;
166 }
167
168 /**
169 * \brief Use exact predicates in convex cell computations.
170 * \details Convex cell computations are used when checking volume
171 * and when saving the cells
172 * \param[in] x true if exact predicates should be used
173 * (default), false otherwise.
174 */
175 2 void use_exact_predicates_for_convex_cell(bool x) {
176 2 convex_cell_exact_predicates_ = x;
177 2 }
178
179 /**
180 * \brief Gets a vertex by index.
181 * \param[in] v a vertex index. Can be a virtual
182 * vertex index when in periodic mode.
183 * \return the 3d point associated with the vertex.
184 * In periodic mode, if \p v is a virtual vertex,
185 * then the translation is applied to the real vertex.
186 */
187 29944 vec3 vertex(index_t v) const {
188
2/2
✓ Branch 0 taken 2477 times.
✓ Branch 1 taken 27467 times.
29944 if(!periodic_) {
189
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 2477 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
2477 geo_debug_assert(v < nb_vertices());
190 2477 return vec3(vertices_ + 3*v);
191 }
192 27467 index_t instance = v/nb_vertices_non_periodic_;
193 27467 v = v%nb_vertices_non_periodic_;
194 27467 vec3 result(vertices_ + 3*v);
195 27467 result.x += double(translation[instance][0]) * period_.x;
196 27467 result.y += double(translation[instance][1]) * period_.y;
197 27467 result.z += double(translation[instance][2]) * period_.z;
198 27467 return result;
199 }
200
201 /**
202 * \brief Gets a weight by index.
203 * \param[in] v a vertex index. Can be a virtual
204 * vertex index in periodic mode.
205 * \return the weight associated with the vertex.
206 */
207 23572 double weight(index_t v) const {
208
1/2
✓ Branch 0 taken 23572 times.
✗ Branch 1 not taken.
23572 if(weights_ == nullptr) {
209 23572 return 0.0;
210 }
211 ✗ return periodic_ ? weights_[periodic_vertex_real(v)] : weights_[v] ;
212 }
213
214 /**
215 * \copydoc Delaunay::nearest_vertex()
216 */
217 index_t nearest_vertex(const double* p) const override;
218
219 /**
220 * \copydoc Delaunay::set_BRIO_levels()
221 */
222 void set_BRIO_levels(const vector<index_t>& levels) override;
223
224 /**
225 * \brief computes the set of tetrahedra that are incident to
226 * a vertex.
227 * \param[in] v the index of the vertex.
228 * \param[in,out] W a reference to a
229 * PeriodicDelaunay3d::IncidentTetrahedra.
230 * On exit it contains the list of incident tets.
231 */
232 void get_incident_tets(index_t v, IncidentTetrahedra& W) const;
233
234 /**
235 * \brief Copies a Laguerre cell from the triangulation.
236 * \details Delaunay neigbhors are stored in ConvexCell vertex global
237 * indices.
238 * \param[in] i the index of the vertex of which the Laguerre cell
239 * should be computed.
240 * \param[out] C the Laguerre cell.
241 * \param[in,out] W a reference to a
242 * PeriodicDelaunay3d::IncidentTetrahedra
243 */
244 void copy_Laguerre_cell_from_Delaunay(
245 GEO::index_t i,
246 ConvexCell& C,
247 IncidentTetrahedra& W
248 ) const;
249
250 /**
251 * \brief Copies a Laguerre cell from the triangulation.
252 * \details Delaunay neigbhors are stored in ConvexCell vertex global
253 * indices.
254 * \param[in] i the index of the vertex of which the Laguerre cell
255 * should be computed.
256 * \param[out] C the Laguerre cell.
257 */
258 void copy_Laguerre_cell_from_Delaunay(
259 GEO::index_t i,
260 ConvexCell& C
261 ) const {
262 IncidentTetrahedra W;
263 copy_Laguerre_cell_from_Delaunay(i,C,W);
264 }
265
266 /**
267 * \brief Sets whether computation should be stopped when an
268 * empty cell is encountered
269 * \param[in] x if set, computation is aborted as soon as an
270 * empty cell is encountered. Default behavior is off.
271 */
272 void abort_if_empty_cell(bool x) {
273 abort_on_empty_cell_ = x;
274 }
275
276 /**
277 * \brief Tests whether the Laguerre diagram has empty cells.
278 * \details If the Laguerre diagram has empty cells, then
279 * computation is stopped, and all the queries on the Laguerre
280 * diagram will not work (including the non-empty cells).
281 * \retval true if the Laguerre diagram has empty cells.
282 * \retval false otherwise.
283 */
284 bool has_empty_cells() const {
285 return has_empty_cells_;
286 }
287
288 /**
289 * \brief Saves the cells in an Alias-Wavefront file.
290 * \param[in] basename the basename of the file. Filename
291 * is basename_callid.obj, where callid is the number of
292 * times the function was invoked.
293 * \param[in] clipped if true, clip the cells by the domain
294 * without saving, else keep the cells as is.
295 */
296 void save_cells(const std::string& basename, bool clipped);
297
298 protected:
299
300 /**
301 * \brief Copies a Laguerre cell facet from the triangulation.
302 * \param[in] i the index of the vertex of which the Laguerre cell
303 * should be computed.
304 * \param[in] Pi the coordinates of vertex \p i
305 * \param[in] wi the weight associated to vertex \p i
306 * \param[in] Pi_len2 the squared length of vertex \p i (considered as
307 * a vector).
308 * \param[in] t a tetrahedron of the Delaunay triangulation,
309 * incident to vertex i
310 * \param[out] C the Laguerre cell.
311 * \param[in,out] W a reference to a
312 * PeriodicDelaunay3d::IncidentTetrahedra
313 * \return the local index of vertex \p i within tetrahedron \p t
314 */
315 GEO::index_t copy_Laguerre_cell_facet_from_Delaunay(
316 GEO::index_t i,
317 const GEO::vec3& Pi,
318 double wi,
319 double Pi_len2,
320 GEO::index_t t,
321 ConvexCell& C,
322 IncidentTetrahedra& W
323 ) const;
324
325
326 /**
327 * \brief Removes unused tetrahedra.
328 * \return the final number of tetrahedra.
329 * \param[in] shrink if true, then array space is shrunk
330 * to fit the new number of tetrahedra.
331 */
332 index_t compress(bool shrink=true);
333
334 /**
335 * \copydoc Delaunay::update_v_to_cell()
336 * \details if update_periodic_v_to_cell_ is set to true,
337 * also updates the map periodic_v_to_cell_ that maps
338 * each virtual vertex to a tet incident to it.
339 */
340 void update_v_to_cell() override;
341
342 /**
343 * \copydoc Delaunay::update_cicl()
344 */
345 void update_cicl() override;
346
347 /**
348 * \brief Duplicates the points with Voronoi cells
349 * that cross the boundary.
350 */
351 void handle_periodic_boundaries();
352
353
354 /**
355 * \brief Phase I of periodic boundaries handling
356 * \details For each cell that traverses the boundary, generates
357 * the periodic vertices instances corresponding to the crossed
358 * faces of the domain, and inserts them in the list of vertices
359 * to be inserted (the reorder_ vector).
360 */
361 void handle_periodic_boundaries_phase_I();
362
363 /**
364 * \brief Tests the position of a Laguerre vertex w.r.t. a plane
365 * \details The positive side of the plane equation corresponds to
366 * what is kept. In other words, the normal vector P.x, P.y, P.z
367 * points towards the interior of this ConvexCell.
368 * \param[in] t a tetrahedron index. The considered Laguerre vertex
369 * is the dual of this tetrahedron.
370 * \param[in] P the plane equation.
371 * \retval true if the Laguerre vertex would be clipped by plane P
372 * \retval false otherwise
373 */
374 bool Laguerre_vertex_is_in_conflict_with_plane(index_t t, vec4 P) const;
375
376 /**
377 * \brief Phase II of periodic boundaries handling
378 * \details Adds the newly discovered neighbors of
379 * the vertices inserted during phase I, back-translated
380 * to their original domain instances in the list of
381 * vertices to be inserted (in the reorder_ member).
382 */
383 void handle_periodic_boundaries_phase_II();
384
385 /**
386 * \brief Inserts vertices from reorder_[b] to reorder_[e-1] using
387 * multithreaded Delaunay. Called by insert_vertices() if there
388 * are many vertices to insert.
389 * \details If an empty cells is detected, has_empty_cells_ is
390 * set and the function exits.
391 */
392 void insert_vertices(const char* phase, index_t b, index_t e);
393
394 /**
395 * \brief Inserts vertices as indicated by a reordering vector
396 * and a vector of BRIO levels, as obtained using
397 * compute_BRIO_order_periodic() (internal function, in the .cpp)
398 * \details used by insert_vertices() (and maybe also compute(), we shall
399 * see if we can). Returns immediatly if an empty cell is encountered,
400 * then it sets has_empty_cells_.
401 * \param[in] levels a const reference to the vector of BRIO levels, as
402 * offsets in the member reordering vector reorder_.
403 */
404 void insert_vertices_with_BRIO(
405 const char* phase, const vector<index_t>& levels
406 );
407
408 /**
409 * \brief Checks the volume of Laguerre cells.
410 */
411 void check_volume();
412
413 /**
414 * \brief Gets a thread by index.
415 * \param[in] t the index of the thread, in 0..nb_threads()-1
416 * \return a pointer to the t-th thread.
417 */
418 114 PeriodicDelaunay3dThread* thread(index_t t) {
419
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 114 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
114 geo_debug_assert(t < threads_.size());
420 return reinterpret_cast<PeriodicDelaunay3dThread*>(
421 114 threads_[t].get()
422 114 );
423 }
424
425 /**
426 * \brief Gets the number of threads.
427 * \return the number of threads.
428 */
429 20 index_t nb_threads() const {
430 20 return index_t(threads_.size());
431 }
432
433 void check_max_t();
434
435 private:
436 friend class PeriodicDelaunay3dThread;
437
438 bool periodic_;
439 vec3 period_;
440
441 const double* weights_;
442 vector<index_t> cell_to_v_store_;
443 vector<index_t> cell_to_cell_store_;
444 vector<index_t> cell_next_;
445
446 CellStatusArray cell_status_;
447
448 ThreadGroup threads_;
449 vector<index_t> reorder_;
450 vector<index_t> levels_;
451
452 /**
453 * Performs additional checks (costly !)
454 */
455 bool debug_mode_;
456
457 /**
458 * Displays the result of the additional checks.
459 */
460 bool verbose_debug_mode_;
461
462 /**
463 * Displays a synthetic summary of timings at the end of the algorithm
464 */
465 bool benchmark_mode_;
466
467 /**
468 * Displays the detailed timing of all the phases of the algorithm
469 */
470 bool detailed_benchmark_mode_;
471
472 /**
473 * \brief Bitmask that indicates for each real vertex
474 * the virtual vertices that were created.
475 */
476 vector<Numeric::uint32> vertex_instances_;
477
478 bool update_periodic_v_to_cell_;
479 vector<index_t> periodic_v_to_cell_rowptr_;
480 vector<index_t> periodic_v_to_cell_data_;
481
482 /**
483 * \brief If set, abort as soon as an empty cell
484 * is encountered.
485 */
486 bool abort_on_empty_cell_;
487
488 /**
489 * \brief Early detection of empty cells.
490 */
491 bool has_empty_cells_;
492
493 /**
494 * \brief Number of reallocations when inserting vertices
495 * in sequential mode.
496 */
497 index_t nb_reallocations_;
498
499 /**
500 * \brief Use exact predicates in convex cell.
501 */
502 bool convex_cell_exact_predicates_;
503
504 struct Stats {
505
506 Stats();
507
508 void reset();
509
510 ✗ std::string to_string() {
511 ✗ return raw_ ? to_string_raw() : to_string_pretty();
512 }
513
514 std::string to_string_raw() const;
515 std::string to_string_pretty() const;
516
517 /**
518 * If set, displays numbers without any formatting.
519 * Set by constructor if command line argument dbg:raw_logs is set.
520 */
521 bool raw_;
522
523 double total_t_;
524
525 double phase_0_t_;
526
527 double phase_I_t_;
528 double phase_I_classify_t_;
529 index_t phase_I_nb_inside_;
530 index_t phase_I_nb_cross_;
531 index_t phase_I_nb_outside_;
532 double phase_I_insert_t_;
533 index_t phase_I_insert_nb_;
534
535 double phase_II_t_;
536 double phase_II_classify_t_;
537 double phase_II_insert_t_;
538 index_t phase_II_insert_nb_;
539 } stats_;
540
541 friend class LaguerreDiagramOmegaSimple3d;
542 };
543
544
545 }
546
547 #endif
548