1 /* 2 * Copyright (c) 2012-2014, Bruno Levy 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 * If you modify this software, you should include a notice giving the 30 * name of the person performing the modification, the date of modification, 31 * and the reason for such modification. 32 * 33 * Contact: Bruno Levy 34 * 35 * Bruno.Levy@inria.fr 36 * http://www.loria.fr/~levy 37 * 38 * ALICE Project 39 * LORIA, INRIA Lorraine, 40 * Campus Scientifique, BP 239 41 * 54506 VANDOEUVRE LES NANCY CEDEX 42 * FRANCE 43 * 44 */ 45 46 #ifndef PERIODIC_DELAUNAY_TRIANGULATION_3D 47 #define PERIODIC_DELAUNAY_TRIANGULATION_3D 48 49 #include <geogram/basic/common.h> 50 #include <geogram/delaunay/delaunay.h> 51 #include <geogram/delaunay/periodic.h> 52 #include <geogram/voronoi/convex_cell.h> 53 #include <geogram/basic/process.h> 54 #include <geogram/basic/geometry.h> 55 #include <stack> 56 57 namespace GEO { 58 59 typedef Numeric::uint8 thread_index_t; 60 class PeriodicDelaunay3dThread; 61 62 63 /** 64 * \brief Multithreaded implementation of Delaunay in 3d with 65 * optional periodic boundary conditions. 66 * \details Periodicity is taken into account by detecting the 67 * vertices with Voronoi cells that straddle the boundary and 68 * duplicating them as need be (with an additional propagation to 69 * have the correct neighborhoods). 70 * \see Delaunay3d, ParallelDelaunay3d 71 */ 72 class GEOGRAM_API PeriodicDelaunay3d : public Delaunay, public Periodic { 73 public: 74 75 /** 76 * \brief Gathers some structures used by some algorithms, makes 77 * multithreading more efficient by avoiding dynamic reallocations. 78 * \details It is used to compute the set of tetrahedra incident to 79 * a vertex. It gathers a stack and the vector of incident tets 80 * obtained so far. 81 */ 82 struct IncidentTetrahedra { 83 std::stack<index_t> S; 84 vector<index_t> incident_tets_set; 85 86 /** 87 * \brief Clears the set of incident tets. 88 */ clear_incident_tetsIncidentTetrahedra89 void clear_incident_tets() { 90 incident_tets_set.resize(0); 91 } 92 93 /** 94 * \brief Inserts a tet into the set of incident tets. 95 * \param[in] t the tet to be inserted. 96 */ add_incident_tetIncidentTetrahedra97 void add_incident_tet(index_t t) { 98 incident_tets_set.push_back(t); 99 } 100 101 /** 102 * \brief Tests whether a tet belongs to the set of incident 103 * tets. 104 * \param[in] t the tet to be tested 105 * \retval true if the tet belongs to the set of incident tets 106 * \retval false otherwise 107 */ has_incident_tetIncidentTetrahedra108 bool has_incident_tet(index_t t) const { 109 for(index_t i=0; i<incident_tets_set.size(); ++i) { 110 if(incident_tets_set[i] == t) { 111 return true; 112 } 113 } 114 return false; 115 } 116 beginIncidentTetrahedra117 vector<index_t>::const_iterator begin() const { 118 return incident_tets_set.begin(); 119 } 120 endIncidentTetrahedra121 vector<index_t>::const_iterator end() const { 122 return incident_tets_set.end(); 123 } 124 }; 125 126 127 /** 128 * \brief Constructs a new PeriodicDelaunay3d. 129 * \param[in] periodic if true, constructs a periodic triangulation. 130 * \param[in] period the edge length of the periodic domain. 131 */ 132 PeriodicDelaunay3d(bool periodic, double period=1.0); 133 134 /** 135 * \copydoc Delaunay::set_vertices() 136 * \note compute() needs to be called after. 137 */ 138 virtual void set_vertices( 139 index_t nb_vertices, const double* vertices 140 ); 141 142 /** 143 * \brief Sets the weights. 144 * \param[in] weights pointer to the array of 145 * weights. Size is the number of real vertices, 146 * i.e., the parameter nb_vertices passed to 147 * set_vertices(). 148 * \note compute() needs to be called after. 149 */ 150 void set_weights(const double* weights); 151 152 /** 153 * \brief Computes the Delaunay triangulation. 154 */ 155 void compute(); 156 157 /** 158 * \brief Use exact predicates in convex cell computations. 159 * \details Convex cell computations are used in periodic 160 * mode for determining the cells that straddle the domain 161 * boundary. 162 * \param[in] x true if exact predicates should be used 163 * (default), false otherwise. 164 */ use_exact_predicates_for_convex_cell(bool x)165 void use_exact_predicates_for_convex_cell(bool x) { 166 convex_cell_exact_predicates_ = x; 167 } 168 169 /** 170 * \brief Gets a vertex by index. 171 * \param[in] v a vertex index. Can be a virtual 172 * vertex index when in periodic mode. 173 * \return the 3d point associated with the vertex. 174 * In periodic mode, if \p v is a virtual vertex, 175 * then the translation is applied to the real vertex. 176 */ vertex(index_t v)177 vec3 vertex(index_t v) const { 178 if(!periodic_) { 179 geo_debug_assert(v < nb_vertices()); 180 return vec3(vertices_ + 3*v); 181 } 182 index_t instance = v/nb_vertices_non_periodic_; 183 v = v%nb_vertices_non_periodic_; 184 vec3 result(vertices_ + 3*v); 185 result.x += double(translation[instance][0]) * period_; 186 result.y += double(translation[instance][1]) * period_; 187 result.z += double(translation[instance][2]) * period_; 188 return result; 189 } 190 191 /** 192 * \brief Gets a weight by index. 193 * \param[in] v a vertex index. Can be a virtual 194 * vertex index in periodic mode. 195 * \return the weight associated with the vertex. 196 */ weight(index_t v)197 double weight(index_t v) const { 198 if(weights_ == nullptr) { 199 return 0.0; 200 } 201 return periodic_ ? weights_[periodic_vertex_real(v)] : weights_[v] ; 202 } 203 204 /** 205 * \copydoc Delaunay::nearest_vertex() 206 */ 207 virtual index_t nearest_vertex(const double* p) const; 208 209 /** 210 * \copydoc Delaunay::set_BRIO_levels() 211 */ 212 virtual void set_BRIO_levels(const vector<index_t>& levels); 213 214 /** 215 * \brief computes the set of tetrahedra that are incident to 216 * a vertex. 217 * \param[in] v the index of the vertex. 218 * \param[in,out] W a reference to a 219 * PeriodicDelaunay3d::IncidentTetrahedra. 220 * On exit it contains the list of incident tets. 221 */ 222 void get_incident_tets(index_t v, IncidentTetrahedra& W) const; 223 224 /** 225 * \brief Copies a Laguerre cell from the triangulation. 226 * \details Delaunay neigbhors are stored in ConvexCell vertex global 227 * indices. 228 * \param[in] i the index of the vertex of which the Laguerre cell 229 * should be computed. 230 * \param[out] C the Laguerre cell. 231 * \param[in,out] W a reference to a 232 * PeriodicDelaunay3d::IncidentTetrahedra 233 */ 234 void copy_Laguerre_cell_from_Delaunay( 235 GEO::index_t i, 236 ConvexCell& C, 237 IncidentTetrahedra& W 238 ) const; 239 240 /** 241 * \brief Copies a Laguerre cell from the triangulation. 242 * \details Delaunay neigbhors are stored in ConvexCell vertex global 243 * indices. 244 * \param[in] i the index of the vertex of which the Laguerre cell 245 * should be computed. 246 * \param[out] C the Laguerre cell. 247 */ copy_Laguerre_cell_from_Delaunay(GEO::index_t i,ConvexCell & C)248 void copy_Laguerre_cell_from_Delaunay( 249 GEO::index_t i, 250 ConvexCell& C 251 ) const { 252 IncidentTetrahedra W; 253 copy_Laguerre_cell_from_Delaunay(i,C,W); 254 } 255 256 /** 257 * \brief Tests whether the Laguerre diagram has empty cells. 258 * \details If the Laguerre diagram has empty cells, then 259 * computation is stopped, and all the queries on the Laguerre 260 * diagram will not work (including the non-empty cells). 261 * \retval true if the Laguerre diagram has empty cells. 262 * \retval false otherwise. 263 */ has_empty_cells()264 bool has_empty_cells() const { 265 return has_empty_cells_; 266 } 267 268 /** 269 * \brief Saves the cells in an Alias-Wavefront file. 270 * \param[in] basename the basename of the file. Filename 271 * is basename_callid.obj, where callid is the number of 272 * times the function was invoked. 273 * \param[in] clipped if true, clip the cells by the domain 274 * without saving, else keep the cells as is. 275 */ 276 void save_cells(const std::string& basename, bool clipped); 277 278 protected: 279 280 /** 281 * \brief Copies a Laguerre cell facet from the triangulation. 282 * \param[in] i the index of the vertex of which the Laguerre cell 283 * should be computed. 284 * \param[in] Pi the coordinates of vertex \p i 285 * \param[in] wi the weight associated to vertex \p i 286 * \param[in] Pi_len2 the squared length of vertex \p i (considered as 287 * a vector). 288 * \param[in] t a tetrahedron of the Delaunay triangulation, 289 * incident to vertex i 290 * \param[out] C the Laguerre cell. 291 * \param[in,out] W a reference to a 292 * PeriodicDelaunay3d::IncidentTetrahedra 293 * \return the local index of vertex \p i within tetrahedron \p t 294 */ 295 GEO::index_t copy_Laguerre_cell_facet_from_Delaunay( 296 GEO::index_t i, 297 const GEO::vec3& Pi, 298 double wi, 299 double Pi_len2, 300 GEO::index_t t, 301 ConvexCell& C, 302 IncidentTetrahedra& W 303 ) const; 304 305 306 /** 307 * \brief Removes unused tetrahedra. 308 * \return the final number of tetrahedra. 309 * \param[in] shrink if true, then array space is shrunk 310 * to fit the new number of tetrahedra. 311 */ 312 index_t compress(bool shrink=true); 313 314 /** 315 * \copydoc Delaunay::update_v_to_cell() 316 * \details if update_periodic_v_to_cell_ is set to true, 317 * also updates the map periodic_v_to_cell_ that maps 318 * each virtual vertex to a tet incident to it. 319 */ 320 virtual void update_v_to_cell(); 321 322 /** 323 * \copydoc Delaunay::update_cicl() 324 */ 325 virtual void update_cicl(); 326 327 /** 328 * \brief Duplicates the points with Voronoi cells 329 * that cross the boundary. 330 */ 331 void handle_periodic_boundaries(); 332 333 /** 334 * \brief Computes the periodic vertex instances 335 * that should be generated. 336 * \param[in] v vertex index, in 0..nb_vertices_non_periodic_-1 337 * \param[out] C the clipped Laguerre cell 338 * \param[out] use_instance the array of booleans that indicates which 339 * instance should be generated. 340 * \param[out] cell_is_on_boundary true if the cell 341 * has an intersection with the cube, false otherwise. 342 * \param[out] cell_is_outside_cube true if the cell 343 * is completely outside the cube, false otherwise. 344 * \param[in,out] W a reference to a 345 * PeriodicDelaunay3d::IncidentTetrahedra 346 * \return the number of instances to generate. 347 */ 348 index_t get_periodic_vertex_instances_to_create( 349 index_t v, 350 ConvexCell& C, 351 bool use_instance[27], 352 bool& cell_is_on_boundary, 353 bool& cell_is_outside_cube, 354 IncidentTetrahedra& W 355 ); 356 357 /** 358 * \brief Insert vertices from 359 * reorder_[b] to reorder_[e-1] 360 * \details If an empty cells is detected, has_empty_cells_ is 361 * set and the function exits. 362 */ 363 void insert_vertices(index_t b, index_t e); 364 365 /** 366 * \brief Checks the volume of Laguerre cells. 367 */ 368 void check_volume(); 369 370 /** 371 * \brief Gets a thread by index. 372 * \param[in] t the index of the thread, in 0..nb_threads()-1 373 * \return a pointer to the t-th thread. 374 */ thread(index_t t)375 PeriodicDelaunay3dThread* thread(index_t t) { 376 geo_debug_assert(t < threads_.size()); 377 return reinterpret_cast<PeriodicDelaunay3dThread*>( 378 threads_[t].get() 379 ); 380 } 381 382 /** 383 * \brief Gets the number of threads. 384 * \return the number of threads. 385 */ nb_threads()386 index_t nb_threads() const { 387 return index_t(threads_.size()); 388 } 389 390 private: 391 friend class PeriodicDelaunay3dThread; 392 393 bool periodic_; 394 double period_; 395 396 const double* weights_; 397 vector<signed_index_t> cell_to_v_store_; 398 vector<signed_index_t> cell_to_cell_store_; 399 vector<index_t> cell_next_; 400 vector<thread_index_t> cell_thread_; 401 ThreadGroup threads_; 402 vector<index_t> reorder_; 403 vector<index_t> levels_; 404 405 /** 406 * Performs additional checks (costly !) 407 */ 408 bool debug_mode_; 409 410 /** 411 * Displays the result of the additional checks. 412 */ 413 bool verbose_debug_mode_; 414 415 /** 416 * Displays the timing of the core algorithm. 417 */ 418 bool benchmark_mode_; 419 420 421 /** 422 * \brief Bitmask that indicates for each real vertex 423 * the virtual vertices that were created. 424 */ 425 vector<Numeric::uint32> vertex_instances_; 426 427 bool update_periodic_v_to_cell_; 428 vector<index_t> periodic_v_to_cell_rowptr_; 429 vector<index_t> periodic_v_to_cell_data_; 430 431 /** 432 * \brief Early detection of empty cells. 433 */ 434 bool has_empty_cells_; 435 436 /** 437 * \brief Number of reallocations when inserting vertices 438 * in sequential mode. 439 */ 440 index_t nb_reallocations_; 441 442 /** 443 * \brief Use exact predicates in convex cell. 444 */ 445 bool convex_cell_exact_predicates_; 446 }; 447 448 449 } 450 451 #endif 452