GCC Code Coverage Report

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			void clear_incident_tets() {
85		✗	incident_tets_set.resize(0);
86			}
87
88			/**
89			* \brief Inserts a tet into the set of incident tets.
90			* \param[in] t the tet to be inserted.
91			*/
92			void add_incident_tet(index_t t) {
93			incident_tets_set.push_back(t);
94			}
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			bool has_incident_tet(index_t t) const {
104		✗	for(index_t i=0; i<incident_tets_set.size(); ++i) {
105		✗	if(incident_tets_set[i] == t) {
106			return true;
107			}
108			}
109			return false;
110			}
111
112			vector<index_t>::const_iterator begin() const {
113			return incident_tets_set.begin();
114			}
115
116			vector<index_t>::const_iterator end() const {
117			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 Use exact predicates in convex cell computations.
160			* \details Convex cell computations are used in periodic
161			* mode for determining the cells that straddle the domain
162			* boundary.
163			* \param[in] x true if exact predicates should be used
164			* (default), false otherwise.
165			*/
166			void use_exact_predicates_for_convex_cell(bool x) {
167			convex_cell_exact_predicates_ = x;
168			}
169
170			/**
171			* \brief Gets a vertex by index.
172			* \param[in] v a vertex index. Can be a virtual
173			* vertex index when in periodic mode.
174			* \return the 3d point associated with the vertex.
175			* In periodic mode, if \p v is a virtual vertex,
176			* then the translation is applied to the real vertex.
177			*/
178		✗	vec3 vertex(index_t v) const {
179		✗	if(!periodic_) {
180			geo_debug_assert(v < nb_vertices());
181		✗	return vec3(vertices_ + 3*v);
182			}
183		✗	index_t instance = v/nb_vertices_non_periodic_;
184		✗	v = v%nb_vertices_non_periodic_;
185		✗	vec3 result(vertices_ + 3*v);
186		✗	result.x += double(translation[instance][0]) * period_.x;
187		✗	result.y += double(translation[instance][1]) * period_.y;
188		✗	result.z += double(translation[instance][2]) * period_.z;
189		✗	return result;
190			}
191
192			/**
193			* \brief Gets a weight by index.
194			* \param[in] v a vertex index. Can be a virtual
195			* vertex index in periodic mode.
196			* \return the weight associated with the vertex.
197			*/
198			double weight(index_t v) const {
199		✗	if(weights_ == nullptr) {
200			return 0.0;
201			}
202		✗	return periodic_ ? weights_[periodic_vertex_real(v)] : weights_[v] ;
203			}
204
205			/**
206			* \copydoc Delaunay::nearest_vertex()
207			*/
208			index_t nearest_vertex(const double* p) const override;
209
210			/**
211			* \copydoc Delaunay::set_BRIO_levels()
212			*/
213			void set_BRIO_levels(const vector<index_t>& levels) override;
214
215			/**
216			* \brief computes the set of tetrahedra that are incident to
217			* a vertex.
218			* \param[in] v the index of the vertex.
219			* \param[in,out] W a reference to a
220			* PeriodicDelaunay3d::IncidentTetrahedra.
221			* On exit it contains the list of incident tets.
222			*/
223			void get_incident_tets(index_t v, IncidentTetrahedra& W) const;
224
225			/**
226			* \brief Copies a Laguerre cell from the triangulation.
227			* \details Delaunay neigbhors are stored in ConvexCell vertex global
228			* indices.
229			* \param[in] i the index of the vertex of which the Laguerre cell
230			* should be computed.
231			* \param[out] C the Laguerre cell.
232			* \param[in,out] W a reference to a
233			* PeriodicDelaunay3d::IncidentTetrahedra
234			*/
235			void copy_Laguerre_cell_from_Delaunay(
236			GEO::index_t i,
237			ConvexCell& C,
238			IncidentTetrahedra& W
239			) const;
240
241			/**
242			* \brief Copies a Laguerre cell from the triangulation.
243			* \details Delaunay neigbhors are stored in ConvexCell vertex global
244			* indices.
245			* \param[in] i the index of the vertex of which the Laguerre cell
246			* should be computed.
247			* \param[out] C the Laguerre cell.
248			*/
249			void copy_Laguerre_cell_from_Delaunay(
250			GEO::index_t i,
251			ConvexCell& C
252			) const {
253			IncidentTetrahedra W;
254			copy_Laguerre_cell_from_Delaunay(i,C,W);
255			}
256
257			/**
258			* \brief Tests whether the Laguerre diagram has empty cells.
259			* \details If the Laguerre diagram has empty cells, then
260			* computation is stopped, and all the queries on the Laguerre
261			* diagram will not work (including the non-empty cells).
262			* \retval true if the Laguerre diagram has empty cells.
263			* \retval false otherwise.
264			*/
265			bool has_empty_cells() const {
266			return has_empty_cells_;
267			}
268
269			/**
270			* \brief Saves the cells in an Alias-Wavefront file.
271			* \param[in] basename the basename of the file. Filename
272			* is basename_callid.obj, where callid is the number of
273			* times the function was invoked.
274			* \param[in] clipped if true, clip the cells by the domain
275			* without saving, else keep the cells as is.
276			*/
277			void save_cells(const std::string& basename, bool clipped);
278
279			protected:
280
281			/**
282			* \brief Copies a Laguerre cell facet from the triangulation.
283			* \param[in] i the index of the vertex of which the Laguerre cell
284			* should be computed.
285			* \param[in] Pi the coordinates of vertex \p i
286			* \param[in] wi the weight associated to vertex \p i
287			* \param[in] Pi_len2 the squared length of vertex \p i (considered as
288			* a vector).
289			* \param[in] t a tetrahedron of the Delaunay triangulation,
290			* incident to vertex i
291			* \param[out] C the Laguerre cell.
292			* \param[in,out] W a reference to a
293			* PeriodicDelaunay3d::IncidentTetrahedra
294			* \return the local index of vertex \p i within tetrahedron \p t
295			*/
296			GEO::index_t copy_Laguerre_cell_facet_from_Delaunay(
297			GEO::index_t i,
298			const GEO::vec3& Pi,
299			double wi,
300			double Pi_len2,
301			GEO::index_t t,
302			ConvexCell& C,
303			IncidentTetrahedra& W
304			) const;
305
306
307			/**
308			* \brief Removes unused tetrahedra.
309			* \return the final number of tetrahedra.
310			* \param[in] shrink if true, then array space is shrunk
311			* to fit the new number of tetrahedra.
312			*/
313			index_t compress(bool shrink=true);
314
315			/**
316			* \copydoc Delaunay::update_v_to_cell()
317			* \details if update_periodic_v_to_cell_ is set to true,
318			* also updates the map periodic_v_to_cell_ that maps
319			* each virtual vertex to a tet incident to it.
320			*/
321			void update_v_to_cell() override;
322
323			/**
324			* \copydoc Delaunay::update_cicl()
325			*/
326			void update_cicl() override;
327
328			/**
329			* \brief Duplicates the points with Voronoi cells
330			* that cross the boundary.
331			*/
332			void handle_periodic_boundaries();
333
334
335			/**
336			* \brief Phase I of periodic boundaries handling
337			* \details For each cell that traverses the boundary, generates
338			* the periodic vertices instances corresponding to the crossed
339			* faces of the domain, and inserts them in the list of vertices
340			* to be inserted (the reorder_ vector).
341			*/
342			void handle_periodic_boundaries_phase_I();
343
344			/**
345			* \brief Tests the position of a Laguerre vertex w.r.t. a plane
346			* \details The positive side of the plane equation corresponds to
347			* what is kept. In other words, the normal vector P.x, P.y, P.z
348			* points towards the interior of this ConvexCell.
349			* \param[in] t a tetrahedron index. The considered Laguerre vertex
350			* is the dual of this tetrahedron.
351			* \param[in] P the plane equation.
352			* \retval true if the Laguerre vertex would be clipped by plane P
353			* \retval false otherwise
354			*/
355			bool Laguerre_vertex_is_in_conflict_with_plane(index_t t, vec4 P) const;
356
357			/**
358			* \brief Phase II of periodic boundaries handling
359			* \details Adds the newly discovered neighbors of
360			* the vertices inserted during phase I, back-translated
361			* to their original domain instances in the list of
362			* vertices to be inserted (in the reorder_ member).
363			*/
364			void handle_periodic_boundaries_phase_II();
365
366			/**
367			* \brief Inserts vertices from reorder_[b] to reorder_[e-1] using
368			* multithreaded Delaunay. Called by insert_vertices() if there
369			* are many vertices to insert.
370			* \details If an empty cells is detected, has_empty_cells_ is
371			* set and the function exits.
372			*/
373			void insert_vertices(const char* phase, index_t b, index_t e);
374
375			/**
376			* \brief Inserts vertices as indicated by a reordering vector
377			* and a vector of BRIO levels, as obtained using
378			* compute_BRIO_order_periodic() (internal function, in the .cpp)
379			* \details used by insert_vertices() (and maybe also compute(), we shall
380			* see if we can). Returns immediatly if an empty cell is encountered,
381			* then it sets has_empty_cells_.
382			* \param[in] levels a const reference to the vector of BRIO levels, as
383			* offsets in the member reordering vector reorder_.
384			*/
385			void insert_vertices_with_BRIO(
386			const char* phase, const vector<index_t>& levels
387			);
388
389			/**
390			* \brief Checks the volume of Laguerre cells.
391			*/
392			void check_volume();
393
394			/**
395			* \brief Gets a thread by index.
396			* \param[in] t the index of the thread, in 0..nb_threads()-1
397			* \return a pointer to the t-th thread.
398			*/
399			PeriodicDelaunay3dThread* thread(index_t t) {
400			geo_debug_assert(t < threads_.size());
401			return reinterpret_cast<PeriodicDelaunay3dThread*>(
402		✗	threads_[t].get()
403			);
404			}
405
406			/**
407			* \brief Gets the number of threads.
408			* \return the number of threads.
409			*/
410			index_t nb_threads() const {
411		✗	return index_t(threads_.size());
412			}
413
414			void check_max_t();
415
416			private:
417			friend class PeriodicDelaunay3dThread;
418
419			bool periodic_;
420			vec3 period_;
421
422			const double* weights_;
423			vector<index_t> cell_to_v_store_;
424			vector<index_t> cell_to_cell_store_;
425			vector<index_t> cell_next_;
426
427			CellStatusArray cell_status_;
428
429			ThreadGroup threads_;
430			vector<index_t> reorder_;
431			vector<index_t> levels_;
432
433			/**
434			* Performs additional checks (costly !)
435			*/
436			bool debug_mode_;
437
438			/**
439			* Displays the result of the additional checks.
440			*/
441			bool verbose_debug_mode_;
442
443			/**
444			* Displays a synthetic summary of timings at the end of the algorithm
445			*/
446			bool benchmark_mode_;
447
448			/**
449			* Displays the detailed timing of all the phases of the algorithm
450			*/
451			bool detailed_benchmark_mode_;
452
453			/**
454			* \brief Bitmask that indicates for each real vertex
455			* the virtual vertices that were created.
456			*/
457			vector<Numeric::uint32> vertex_instances_;
458
459			bool update_periodic_v_to_cell_;
460			vector<index_t> periodic_v_to_cell_rowptr_;
461			vector<index_t> periodic_v_to_cell_data_;
462
463			/**
464			* \brief Early detection of empty cells.
465			*/
466			bool has_empty_cells_;
467
468			/**
469			* \brief Number of reallocations when inserting vertices
470			* in sequential mode.
471			*/
472			index_t nb_reallocations_;
473
474			/**
475			* \brief Use exact predicates in convex cell.
476			*/
477			bool convex_cell_exact_predicates_;
478
479			struct Stats {
480
481			Stats();
482
483			void reset();
484
485		✗	std::string to_string() {
486		✗	return raw_ ? to_string_raw() : to_string_pretty();
487			}
488
489			std::string to_string_raw() const;
490			std::string to_string_pretty() const;
491
492			/**
493			* If set, displays numbers without any formatting.
494			* Set by constructor if command line argument dbg:raw_logs is set.
495			*/
496			bool raw_;
497
498			double total_t_;
499
500			double phase_0_t_;
501
502			double phase_I_t_;
503			double phase_I_classify_t_;
504			index_t phase_I_nb_inside_;
505			index_t phase_I_nb_cross_;
506			index_t phase_I_nb_outside_;
507			double phase_I_insert_t_;
508			index_t phase_I_insert_nb_;
509
510			double phase_II_t_;
511			double phase_II_classify_t_;
512			double phase_II_insert_t_;
513			index_t phase_II_insert_nb_;
514			} stats_;
515
516			friend class LaguerreDiagramOmegaSimple3d;
517			};
518
519
520			}
521
522			#endif
523