Geogram  Version 1.9.1-rc
A programming library of geometric algorithms
periodic_delaunay_3d.h
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>
47 #include <geogram/basic/process.h>
48 #include <geogram/basic/geometry.h>
49 #include <stack>
50 
52 
53 namespace GEO {
54 
55  typedef Numeric::uint8 thread_index_t;
56  class PeriodicDelaunay3dThread;
57 
58 
68  class GEOGRAM_API PeriodicDelaunay3d : public Delaunay, public Periodic {
69  public:
70 
79  std::stack<index_t> S;
80  vector<index_t> incident_tets_set;
81 
86  incident_tets_set.resize(0);
87  }
88 
94  incident_tets_set.push_back(t);
95  }
96 
104  bool has_incident_tet(index_t t) const {
105  for(index_t i=0; i<incident_tets_set.size(); ++i) {
106  if(incident_tets_set[i] == t) {
107  return true;
108  }
109  }
110  return false;
111  }
112 
113  vector<index_t>::const_iterator begin() const {
114  return incident_tets_set.begin();
115  }
116 
117  vector<index_t>::const_iterator end() const {
118  return incident_tets_set.end();
119  }
120  };
121 
127  PeriodicDelaunay3d(bool periodic, double period=1.0);
128 
133  PeriodicDelaunay3d(const vec3& period);
134 
140  index_t nb_vertices, const double* vertices
141  ) override;
142 
151  void set_weights(const double* weights);
152 
156  void compute();
157 
167  convex_cell_exact_predicates_ = x;
168  }
169 
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 
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 
208  index_t nearest_vertex(const double* p) const override;
209 
213  void set_BRIO_levels(const vector<index_t>& levels) override;
214 
224 
236  GEO::index_t i,
237  ConvexCell& C,
239  ) const;
240 
250  GEO::index_t i,
251  ConvexCell& C
252  ) const {
254  copy_Laguerre_cell_from_Delaunay(i,C,W);
255  }
256 
265  bool has_empty_cells() const {
266  return has_empty_cells_;
267  }
268 
277  void save_cells(const std::string& basename, bool clipped);
278 
279  protected:
280 
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,
304  ) const;
305 
306 
313  index_t compress(bool shrink=true);
314 
321  void update_v_to_cell() override;
322 
326  void update_cicl() override;
327 
333 
350  index_t v,
351  ConvexCell& C,
352  bool use_instance[27],
353  bool& cell_is_on_boundary,
354  bool& cell_is_outside_cube,
356  );
357 
365 
369  void check_volume();
370 
376  PeriodicDelaunay3dThread* thread(index_t t) {
377  geo_debug_assert(t < threads_.size());
378  return reinterpret_cast<PeriodicDelaunay3dThread*>(
379  threads_[t].get()
380  );
381  }
382 
387  index_t nb_threads() const {
388  return index_t(threads_.size());
389  }
390 
391  private:
392  friend class PeriodicDelaunay3dThread;
393 
394  bool periodic_;
395  vec3 period_;
396 
397  const double* weights_;
398  vector<signed_index_t> cell_to_v_store_;
399  vector<signed_index_t> cell_to_cell_store_;
400  vector<index_t> cell_next_;
401 
402  CellStatusArray cell_status_;
403 
404  ThreadGroup threads_;
405  vector<index_t> reorder_;
406  vector<index_t> levels_;
407 
411  bool debug_mode_;
412 
416  bool verbose_debug_mode_;
417 
421  bool benchmark_mode_;
422 
423 
428  vector<Numeric::uint32> vertex_instances_;
429 
430  bool update_periodic_v_to_cell_;
431  vector<index_t> periodic_v_to_cell_rowptr_;
432  vector<index_t> periodic_v_to_cell_data_;
433 
437  bool has_empty_cells_;
438 
443  index_t nb_reallocations_;
444 
448  bool convex_cell_exact_predicates_;
449  };
450 
451 
452 }
453 
454 #endif
#define geo_debug_assert(x)
Verifies that a condition is met.
Definition: assert.h:196
An array of cell status codes associates to each tetrahedron in a Delaunay tetrahedralization.
Definition: delaunay_sync.h:62
Abstract interface for Delaunay triangulation in Nd.
Definition: delaunay.h:71
Multithreaded implementation of Delaunay in 3d with optional periodic boundary conditions.
PeriodicDelaunay3d(const vec3 &period)
Constructs a new PeriodicDelaunay3d.
PeriodicDelaunay3d(bool periodic, double period=1.0)
Constructs a new PeriodicDelaunay3d.
void check_volume()
Checks the volume of Laguerre cells.
void compute()
Computes the Delaunay triangulation.
void use_exact_predicates_for_convex_cell(bool x)
Use exact predicates in convex cell computations.
void handle_periodic_boundaries()
Duplicates the points with Voronoi cells that cross the boundary.
void save_cells(const std::string &basename, bool clipped)
Saves the cells in an Alias-Wavefront file.
void set_BRIO_levels(const vector< index_t > &levels) override
Specifies the bounds of each level to be used when hierarchic ordering is specified from outside.
void insert_vertices(index_t b, index_t e)
Insert vertices from reorder_[b] to reorder_[e-1].
index_t compress(bool shrink=true)
Removes unused tetrahedra.
double weight(index_t v) const
Gets a weight by index.
PeriodicDelaunay3dThread * thread(index_t t)
Gets a thread by index.
void update_v_to_cell() override
Stores for each vertex v a cell incident to v.
void update_cicl() override
Updates the circular incident cell lists.
void get_incident_tets(index_t v, IncidentTetrahedra &W) const
computes the set of tetrahedra that are incident to a vertex.
void copy_Laguerre_cell_from_Delaunay(GEO::index_t i, ConvexCell &C, IncidentTetrahedra &W) const
Copies a Laguerre cell from the triangulation.
index_t nb_threads() const
Gets the number of threads.
index_t nearest_vertex(const double *p) const override
Computes the nearest vertex from a query point.
void copy_Laguerre_cell_from_Delaunay(GEO::index_t i, ConvexCell &C) const
Copies a Laguerre cell from the triangulation.
index_t get_periodic_vertex_instances_to_create(index_t v, ConvexCell &C, bool use_instance[27], bool &cell_is_on_boundary, bool &cell_is_outside_cube, IncidentTetrahedra &W)
Computes the periodic vertex instances that should be generated.
vec3 vertex(index_t v) const
Gets a vertex by index.
bool has_empty_cells() const
Tests whether the Laguerre diagram has empty cells.
GEO::index_t copy_Laguerre_cell_facet_from_Delaunay(GEO::index_t i, const GEO::vec3 &Pi, double wi, double Pi_len2, GEO::index_t t, ConvexCell &C, IncidentTetrahedra &W) const
Copies a Laguerre cell facet from the triangulation.
void set_weights(const double *weights)
Sets the weights.
void set_vertices(index_t nb_vertices, const double *vertices) override
Sets the vertices of this Delaunay, and recomputes the cells.
Utilities for managing 3D periodic space.
Definition: periodic.h:57
index_t size() const
Gets the number of elements.
Definition: memory.h:674
Computes the intersection between a set of halfplanes using Bowyer-Watson algorithm.
Definition: convex_cell.h:419
Class to compute the intersection of a set of half-spaces in 3D.
Abstract interface for Delaunay.
Synchronization primitives for parallel Delaunay.
Common include file, providing basic definitions. Should be included before anything else by all head...
Geometric functions in 2d and 3d.
uint8_t uint8
Definition: numeric.h:135
Global Vorpaline namespace.
Definition: basic.h:55
std::vector< Thread_var > ThreadGroup
Collection of Threads.
Definition: process.h:155
vecng< 3, Numeric::float64 > vec3
Represents points and vectors in 3d.
Definition: geometry.h:65
geo_index_t index_t
The type for storing and manipulating indices.
Definition: numeric.h:329
Manipulation of indices for 3D periodic space.
Function and classes for process manipulation.
Gathers some structures used by some algorithms, makes multithreading more efficient by avoiding dyna...
void add_incident_tet(index_t t)
Inserts a tet into the set of incident tets.
void clear_incident_tets()
Clears the set of incident tets.
bool has_incident_tet(index_t t) const
Tests whether a tet belongs to the set of incident tets.