| 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_DELAUNAY_DELAUNAY_3D | ||
| 41 | #define GEOGRAM_DELAUNAY_DELAUNAY_3D | ||
| 42 | |||
| 43 | #include <geogram/basic/common.h> | ||
| 44 | #include <geogram/delaunay/delaunay.h> | ||
| 45 | #include <geogram/delaunay/cavity.h> | ||
| 46 | #include <geogram/numerics/predicates.h> | ||
| 47 | #include <geogram/basic/geometry.h> | ||
| 48 | |||
| 49 | #include <stack> | ||
| 50 | |||
| 51 | /** | ||
| 52 | * \file geogram/delaunay/delaunay_3d.h | ||
| 53 | * \brief Implementation of Delaunay in 3d. | ||
| 54 | */ | ||
| 55 | |||
| 56 | namespace GEO { | ||
| 57 | |||
| 58 | /** | ||
| 59 | * \brief Implementation of Delaunay in 3d. | ||
| 60 | * \details This package uses concepts inspired by | ||
| 61 | * two triangulation softwares, CGAL and tetgen, | ||
| 62 | * described in the following references. This package follows the | ||
| 63 | * idea used in CGAL of traversing the cavity from inside, since | ||
| 64 | * it traverses less tetrahedra than when traversing from outside. | ||
| 65 | * - Jean-Daniel Boissonnat, Olivier Devillers, Monique Teillaud, | ||
| 66 | * and Mariette Yvinec. Triangulations in CGAL. | ||
| 67 | * In Proc. 16th Annu. ACM Sympos. Comput. Geom., pages 11-18, 2000. | ||
| 68 | * - Hang Si, Constrained Delaunay tetrahedral mesh generation and | ||
| 69 | * refinement. Finite elements in Analysis and Design, | ||
| 70 | * 46 (1-2):33--46, 2010. | ||
| 71 | * | ||
| 72 | * Note that the algorithm here does not support vertex deletion nor | ||
| 73 | * degenerate input with all coplanar or all colinear points (use CGAL | ||
| 74 | * instead if you have these requirements). | ||
| 75 | * | ||
| 76 | * The core algorithm used in this code, CGAL and tetgen was | ||
| 77 | * independently and simultaneously discovered by Bowyer and Watson: | ||
| 78 | * - Adrian Bowyer, "Computing Dirichlet tessellations", | ||
| 79 | * Comput. J., vol. 24, no 2, 1981, p. 162-166 | ||
| 80 | * - David F. Watson, "Computing the n-dimensional Delaunay tessellation | ||
| 81 | * with application to Voronoi polytopes", Comput. J., vol. 24, | ||
| 82 | * no 2, 1981, p. 167-172 | ||
| 83 | * | ||
| 84 | * The spatial reordering method, that dramatically increases the | ||
| 85 | * performances, also used in this code, CGAL and tetgen was introduced | ||
| 86 | * in the following references. The second one is a smart implementation | ||
| 87 | * based on the std::nth_element() function of the STL, that inspired | ||
| 88 | * the compute_BRIO_ordering() function of this package. | ||
| 89 | * - Nina Amenta, Sunghee Choi and Gunter Rote, "Incremental constructions | ||
| 90 | * con brio", ACM Symposium on Computational Geometry 2003. | ||
| 91 | * - Christophe Delage and Olivier Devillers. Spatial Sorting. | ||
| 92 | * In CGAL User and Reference Manual. CGAL Editorial Board, | ||
| 93 | * 3.9 edition, 2011 | ||
| 94 | * | ||
| 95 | * The locate() function is based on the following two references. | ||
| 96 | * The first one randomizes the choice of the next tetrahedron. | ||
| 97 | * The second one uses an inexact locate() function to initialize | ||
| 98 | * the exact one (it is called "structural filtering"). The first | ||
| 99 | * idea is used in both CGAL and tetgen, and the second one is used | ||
| 100 | * in CGAL. | ||
| 101 | * - Walking in a triangulation, O Devillers, S Pion, M Teillaud | ||
| 102 | * 17th Annual Symposium on Computational geometry, 106-114 | ||
| 103 | * - Stefan Funke , Kurt Mehlhorn and Stefan Naher, "Structural filtering, | ||
| 104 | * a paradigm for efficient and exact geometric programs", | ||
| 105 | * Comput. Geom., 1999 | ||
| 106 | */ | ||
| 107 | class GEOGRAM_API Delaunay3d : public Delaunay { | ||
| 108 | public: | ||
| 109 | /** | ||
| 110 | * \brief Constructs a new Delaunay3d. | ||
| 111 | * \param[in] dimension dimension of the triangulation (3 or 4). | ||
| 112 | * If dimension = 4, this creates a regular triangulation | ||
| 113 | * (dual of a power diagram). In this case: | ||
| 114 | * - the input points are 4d points, were the fourth coordinate | ||
| 115 | * of point \f$ i \f$ is \f$ \sqrt{W - w_i} \f$ where \f$ W \f$ is | ||
| 116 | * the maximum of the weights of all the points and \d$ w_i \$ is | ||
| 117 | * the weight associated with vertex \f$ i \f$. | ||
| 118 | * - the constructed combinatorics is a tetrahedralized volume (3d and | ||
| 119 | * not 4d although dimension() returns 4). This tetrahedralized volume | ||
| 120 | * corresponds to the regular triangulation of the weighted points. | ||
| 121 | */ | ||
| 122 | Delaunay3d(coord_index_t dimension = 3); | ||
| 123 | |||
| 124 | /** | ||
| 125 | * \copydoc Delaunay::set_vertices() | ||
| 126 | */ | ||
| 127 | void set_vertices(index_t nb_vertices, const double* vertices) override; | ||
| 128 | |||
| 129 | /** | ||
| 130 | * \copydoc Delaunay::nearest_vertex() | ||
| 131 | */ | ||
| 132 | index_t nearest_vertex(const double* p) const override; | ||
| 133 | |||
| 134 | protected: | ||
| 135 | |||
| 136 | /** | ||
| 137 | * \brief Symbolic constant for uninitialized hint. | ||
| 138 | * \details Locate functions can be accelerated by | ||
| 139 | * specifying a hint. This constant indicates that | ||
| 140 | * no hint is given. | ||
| 141 | */ | ||
| 142 | static constexpr index_t NO_TETRAHEDRON = NO_INDEX; | ||
| 143 | |||
| 144 | /** | ||
| 145 | * \brief Finds in the pointset a set of four non-coplanar | ||
| 146 | * points. | ||
| 147 | * \details This function is used to initiate the incremental | ||
| 148 | * Delaunay construction. | ||
| 149 | * \param[out] iv0 index of the first vertex | ||
| 150 | * \param[out] iv1 index of the second vertex | ||
| 151 | * \param[out] iv2 index of the third vertex | ||
| 152 | * \param[out] iv3 index of the fourth vertex | ||
| 153 | * \retval true if a set of four non-coplanar points was found | ||
| 154 | * \retval false if all the points are coplanar | ||
| 155 | */ | ||
| 156 | bool create_first_tetrahedron( | ||
| 157 | index_t& iv0, index_t& iv1, index_t& iv2, index_t& iv3 | ||
| 158 | ); | ||
| 159 | |||
| 160 | /** | ||
| 161 | * \brief Finds the tetrahedron that contains a point. | ||
| 162 | * \details If the point is on a face, edge or vertex, | ||
| 163 | * the function returns one of the tetrahedra incident | ||
| 164 | * to that face, edge or vertex. | ||
| 165 | * \param[in] p a pointer to the coordinates of the point | ||
| 166 | * \param[in] thread_safe if true, a global spinlock is | ||
| 167 | * used to protect the calls to random(), this is necessary | ||
| 168 | * if multiple threads use locate() simultaneously | ||
| 169 | * \param[out] orient a pointer to an array of four Sign%s | ||
| 170 | * or nullptr. If non-nullptr, returns the orientation with respect | ||
| 171 | * to the four facets of the tetrahedron that contains \p p. | ||
| 172 | * \return the index of a tetrahedron that contains \p p. | ||
| 173 | * If the point is outside the convex hull of | ||
| 174 | * the inserted so-far points, then the returned tetrahedron | ||
| 175 | * is a virtual one (first vertex is the "vertex at infinity" | ||
| 176 | * of index -1) or NO_TETRAHEDRON if the virtual tetrahedra | ||
| 177 | * were previously removed. | ||
| 178 | */ | ||
| 179 | index_t locate( | ||
| 180 | const double* p, index_t hint = NO_TETRAHEDRON, | ||
| 181 | bool thread_safe = false, | ||
| 182 | Sign* orient = nullptr | ||
| 183 | ) const; | ||
| 184 | |||
| 185 | /** | ||
| 186 | * \brief Finds the tetrahedron that (approximately) | ||
| 187 | * contains a point using inexact predicates. | ||
| 188 | * \details The result of this function can be used as a hint | ||
| 189 | * for locate(). It accelerates locate as compared to calling | ||
| 190 | * it directly. This technique is referred to as "structural | ||
| 191 | * filtering". | ||
| 192 | * \param[in] p a pointer to the coordinates of the point | ||
| 193 | * \param[in] max_iter maximum number of traversed tets | ||
| 194 | * \return the index of a tetrahedron that (approximately) | ||
| 195 | * contains \p p. | ||
| 196 | * If the point is outside the convex hull of | ||
| 197 | * the inserted so-far points, then the returned tetrahedron | ||
| 198 | * is a virtual one (first vertex is the "vertex at infinity" | ||
| 199 | * of index -1) or NO_TETRAHEDRON if the virtual tetrahedra | ||
| 200 | * were previously removed. | ||
| 201 | */ | ||
| 202 | index_t locate_inexact( | ||
| 203 | const double* p, index_t hint, index_t max_iter | ||
| 204 | ) const; | ||
| 205 | |||
| 206 | /** | ||
| 207 | * \brief Inserts a point in the triangulation. | ||
| 208 | * \param[in] v the index of the point to be inserted | ||
| 209 | * \param[in] hint the index of a tetrahedron as near as | ||
| 210 | * possible to \p v, or -1 if unspecified | ||
| 211 | * \return the index of one of the tetrahedra incident to | ||
| 212 | * point \p v | ||
| 213 | */ | ||
| 214 | index_t insert(index_t v, index_t hint = NO_TETRAHEDRON); | ||
| 215 | |||
| 216 | /** | ||
| 217 | * \brief Determines the list of tetrahedra in conflict | ||
| 218 | * with a given point. | ||
| 219 | * \param[in] v the index of the point to be inserted | ||
| 220 | * \param[in] t the index of a tetrahedron that contains | ||
| 221 | * \p p, as returned by locate() | ||
| 222 | * \param[in] orient an array of four signs indicating | ||
| 223 | * the orientation of \p p with respect to the four | ||
| 224 | * faces of \p t, as returned by locate() | ||
| 225 | * \param[out] t_bndry a tetrahedron adjacent to the | ||
| 226 | * boundary of the conflict zone | ||
| 227 | * \param[out] f_bndry the facet along which t_bndry is | ||
| 228 | * adjacent to the boundary of the conflict zone | ||
| 229 | * \param[out] first the index of the first tetrahedron in conflict | ||
| 230 | * \param[out] last the index of the last tetrahedron in conflict | ||
| 231 | * The other tetrahedra are linked, and can be traversed | ||
| 232 | * from \p first by using tet_next() until \p last or END_OF_LIST | ||
| 233 | * is reached. | ||
| 234 | * The conflict zone can be empty under two circumstances: | ||
| 235 | * - the vertex \p v already exists in the triangulation | ||
| 236 | * - the triangulation is weighted and \p v is not visible | ||
| 237 | * in either cases, both \p first and \p last contain END_OF_LIST | ||
| 238 | */ | ||
| 239 | void find_conflict_zone( | ||
| 240 | index_t v, | ||
| 241 | index_t t, const Sign* orient, | ||
| 242 | index_t& t_bndry, index_t& f_bndry, | ||
| 243 | index_t& first, index_t& last | ||
| 244 | ); | ||
| 245 | |||
| 246 | /** | ||
| 247 | * \brief This function is used to implement find_conflict_zone. | ||
| 248 | * \details This function detects the neighbors of \p t that are | ||
| 249 | * in the conflict zone and calls itself recursively on them. | ||
| 250 | * \param[in] p the point to be inserted | ||
| 251 | * \param[in] t index of a tetrahedron in the fonflict zone | ||
| 252 | * \param[out] t_bndry a tetrahedron adjacent to the | ||
| 253 | * boundary of the conflict zone | ||
| 254 | * \param[out] f_bndry the facet along which t_bndry is | ||
| 255 | * adjacent to the boundary of the conflict zone | ||
| 256 | * \param[out] first the index of the first tetrahedron in conflict | ||
| 257 | * \param[out] last the index of the last tetrahedron in conflict | ||
| 258 | * \pre The tetrahedron \p t was alredy marked as | ||
| 259 | * conflict (tet_is_in_list(t)) | ||
| 260 | */ | ||
| 261 | void find_conflict_zone_iterative( | ||
| 262 | const double* p, index_t t, | ||
| 263 | index_t& t_bndry, index_t& f_bndry, | ||
| 264 | index_t& first, index_t& last | ||
| 265 | ); | ||
| 266 | |||
| 267 | /** | ||
| 268 | * \brief Creates a star of tetrahedra filling the conflict | ||
| 269 | * zone. | ||
| 270 | * \param[in] v the index of the point to be inserted | ||
| 271 | * \details This function is used when the Cavity computed | ||
| 272 | * when traversing the conflict zone is OK, that is to say | ||
| 273 | * when its array sizes were not exceeded. | ||
| 274 | * \return the index of one the newly created tetrahedron | ||
| 275 | */ | ||
| 276 | index_t stellate_cavity(index_t v); | ||
| 277 | |||
| 278 | |||
| 279 | /** | ||
| 280 | * \brief Creates a star of tetrahedra filling the conflict | ||
| 281 | * zone. | ||
| 282 | * \details For each tetrahedron facet on the border of the | ||
| 283 | * conflict zone, a new tetrahedron is created, resting on | ||
| 284 | * the facet and incident to vertex \p v. The function is | ||
| 285 | * called recursively until the entire conflict zone is filled. | ||
| 286 | * \param[in] v the index of the point to be inserted | ||
| 287 | * \param[in] t_bndry index of a tetrahedron on the border | ||
| 288 | * of the conflict zone. | ||
| 289 | * \param[in] f_bndry index of the facet along which \p t_bndry | ||
| 290 | * is incident to the border of the conflict zone | ||
| 291 | * \param[in] prev_f the facet of \p t_bndry connected to the | ||
| 292 | * tetrahedron that \p t_bndry was reached from, or NO_INDEX | ||
| 293 | * if it is the first tetrahedron. | ||
| 294 | * \return the index of one the newly created tetrahedron | ||
| 295 | */ | ||
| 296 | index_t stellate_conflict_zone_iterative( | ||
| 297 | index_t v, | ||
| 298 | index_t t_bndry, index_t f_bndry, | ||
| 299 | index_t prev_f=NO_INDEX | ||
| 300 | ); | ||
| 301 | |||
| 302 | /** | ||
| 303 | * \brief Finds the neighbor of a tetrahedron on the border of the | ||
| 304 | * conflict zone. | ||
| 305 | * \details This function is used by stellate_conflict_zone_iterative() | ||
| 306 | * \param[in] t1 a tetrahedron on the border of the conflict zone | ||
| 307 | * \param[in] t1fborder the local facet index of \p t1 along which it | ||
| 308 | * is on the border of the conflict zone | ||
| 309 | * \param[in] t1ft2 the local facet index of \p t1 that will be | ||
| 310 | * traversed | ||
| 311 | * \param[out] t2 a tetrahedron on the border of the conflict zone, | ||
| 312 | * with an edge common to facets \p t1fborder and \p t1ft2 of | ||
| 313 | * tetrahedron \p t1 | ||
| 314 | * \param[out] t2fborder the local facet index of \p t2 along which it | ||
| 315 | * is on the border of the conflict zone | ||
| 316 | * \param[out] t2ft1 the local index of the facet of \p t2 that has a | ||
| 317 | * common edge with facets \p t1fborder and \p t1ft2 of tetrahedron | ||
| 318 | * \p t1 | ||
| 319 | * \retval true if \p t2 is a newly created tetrahedron | ||
| 320 | * \retval false if \p t2 is an old tetrahedron in conflict | ||
| 321 | */ | ||
| 322 | 750 | bool get_neighbor_along_conflict_zone_border( | |
| 323 | index_t t1, | ||
| 324 | index_t t1fborder, | ||
| 325 | index_t t1ft2, | ||
| 326 | index_t& t2, | ||
| 327 | index_t& t2fborder, | ||
| 328 | index_t& t2ft1 | ||
| 329 | ) const { | ||
| 330 | |||
| 331 | // Note: this function is a bit long for an inline function, | ||
| 332 | // but I observed a (modest) performance gain doing so. | ||
| 333 | |||
| 334 | // Find two vertices that are both on facets new_f and f1 | ||
| 335 | // (the edge around which we are turning) | ||
| 336 | // This uses duality as follows: | ||
| 337 | // Primal form (not used here): | ||
| 338 | // halfedge_facet_[v1][v2] returns a facet that is incident | ||
| 339 | // to both v1 and v2. | ||
| 340 | // Dual form (used here): | ||
| 341 | // halfedge_facet_[f1][f2] returns a vertex that both | ||
| 342 | // f1 and f2 are incident to. | ||
| 343 | index_t ev1 = | ||
| 344 |
1/2✓ Branch 1 taken 750 times.
✗ Branch 2 not taken.
|
750 | tet_vertex(t1, index_t(halfedge_facet_[t1ft2][t1fborder])); |
| 345 | index_t ev2 = | ||
| 346 |
1/2✓ Branch 1 taken 750 times.
✗ Branch 2 not taken.
|
750 | tet_vertex(t1, index_t(halfedge_facet_[t1fborder][t1ft2])); |
| 347 | |||
| 348 | // Turn around edge [ev1,ev2] inside the conflict zone | ||
| 349 | // until we reach again the boundary of the conflict zone. | ||
| 350 | // Traversing inside the conflict zone is faster (as compared | ||
| 351 | // to outside) since it traverses a smaller number of tets. | ||
| 352 | 750 | index_t cur_t = t1; | |
| 353 | 750 | index_t cur_f = t1ft2; | |
| 354 |
1/2✓ Branch 1 taken 750 times.
✗ Branch 2 not taken.
|
750 | index_t next_t = tet_adjacent(cur_t,cur_f); |
| 355 |
3/4✓ Branch 1 taken 1893 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1143 times.
✓ Branch 4 taken 750 times.
|
1893 | while(tet_is_in_list(next_t)) { |
| 356 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 1143 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1143 | geo_debug_assert(next_t != t1); |
| 357 | 1143 | cur_t = next_t; | |
| 358 |
1/2✓ Branch 1 taken 1143 times.
✗ Branch 2 not taken.
|
1143 | cur_f = get_facet_by_halfedge(cur_t,ev1,ev2); |
| 359 |
1/2✓ Branch 1 taken 1143 times.
✗ Branch 2 not taken.
|
1143 | next_t = tet_adjacent(cur_t, cur_f); |
| 360 | } | ||
| 361 | |||
| 362 | // At this point, cur_t is in conflict zone and | ||
| 363 | // next_t is outside the conflict zone. | ||
| 364 | index_t f12,f21; | ||
| 365 |
1/2✓ Branch 1 taken 750 times.
✗ Branch 2 not taken.
|
750 | get_facets_by_halfedge(next_t, ev1, ev2, f12, f21); |
| 366 |
1/2✓ Branch 1 taken 750 times.
✗ Branch 2 not taken.
|
750 | t2 = tet_adjacent(next_t,f21); |
| 367 |
1/2✓ Branch 1 taken 750 times.
✗ Branch 2 not taken.
|
750 | index_t v_neigh_opposite = tet_vertex(next_t,f12); |
| 368 |
1/2✓ Branch 1 taken 750 times.
✗ Branch 2 not taken.
|
750 | t2ft1 = find_tet_vertex(t2, v_neigh_opposite); |
| 369 | 750 | t2fborder = cur_f; | |
| 370 | |||
| 371 | // Test whether the found neighboring tet was created | ||
| 372 | // (then return true) or is an old tet in conflict | ||
| 373 | // (then return false). | ||
| 374 | 750 | return(t2 != cur_t); | |
| 375 | } | ||
| 376 | |||
| 377 | /****** Combinatorics - new and delete ***************************/ | ||
| 378 | |||
| 379 | /** | ||
| 380 | * \brief Maximum valid index for a tetrahedron. | ||
| 381 | * \details This includes not only real tetrahedra, | ||
| 382 | * but also the virtual ones on the border, the conflict | ||
| 383 | * list and the free list. | ||
| 384 | * \return the maximum valid index for a tetrahedron | ||
| 385 | */ | ||
| 386 | 198779 | index_t max_t() const { | |
| 387 | 198779 | return cell_to_v_store_.size() / 4; | |
| 388 | } | ||
| 389 | |||
| 390 | /** | ||
| 391 | * \brief Default symbolic value of the cell_next_ field | ||
| 392 | * that indicates that a tetrahedron is not | ||
| 393 | * in a linked list. | ||
| 394 | * \details This is the default value. Note that it suffices | ||
| 395 | * that NOT_IN_LIST_BIT is set for a tetrahedron | ||
| 396 | * to be not in any list. | ||
| 397 | * A tetrahedron can be: | ||
| 398 | * - in a list (cell_next_[t] & NOT_IN_LIST_BIT == 0) | ||
| 399 | * - not in a list and not marked | ||
| 400 | * (cell_next_[t] & NOT_IN_LIST_BIT != 0) && | ||
| 401 | * (cell_next_[t] != cur_stamp_) | ||
| 402 | * - not in a list and marked | ||
| 403 | * (cell_next_[t] == cur_stamp_) | ||
| 404 | */ | ||
| 405 | static constexpr index_t NOT_IN_LIST = ~index_t(0); | ||
| 406 | |||
| 407 | /** | ||
| 408 | * \brief If cell_next_[t] & NOT_IN_LIST_BIT != 0, | ||
| 409 | * then t is not in a linked list. | ||
| 410 | * \details The other bits of cell_next_[t] are used | ||
| 411 | * to store the stamp (i.e. index of the current point | ||
| 412 | * being inserted). The stamp is used for marking tetrahedra | ||
| 413 | * that were detected as non-conflict when inserting a point. | ||
| 414 | * A tetrahedron can be: | ||
| 415 | * - in a list (cell_next_[t] & NOT_IN_LIST_BIT == 0) | ||
| 416 | * - not in a list and not marked | ||
| 417 | * (cell_next_[t] & NOT_IN_LIST_BIT != 0) && | ||
| 418 | * (cell_next_[t] != cur_stamp_) | ||
| 419 | * - not in a list and marked | ||
| 420 | * (cell_next_[t] == cur_stamp_) | ||
| 421 | */ | ||
| 422 | static constexpr index_t NOT_IN_LIST_BIT = | ||
| 423 | index_t(1) << (sizeof(index_t)*8-1) ; | ||
| 424 | |||
| 425 | /** | ||
| 426 | * \brief Symbolic value of the cell_next_ field | ||
| 427 | * that indicates the end of list in a linked | ||
| 428 | * list of tetrahedra. | ||
| 429 | */ | ||
| 430 | static constexpr index_t END_OF_LIST = ~NOT_IN_LIST_BIT; | ||
| 431 | |||
| 432 | |||
| 433 | /** | ||
| 434 | * \brief Tests whether a tetrahedron belongs to a linked | ||
| 435 | * list. | ||
| 436 | * \details Tetrahedra can be linked, it is used to manage | ||
| 437 | * both the free list that recycles deleted tetrahedra, | ||
| 438 | * the conflict region and the list of newly created | ||
| 439 | * tetrahedra. In addition, a tetrahedron that is not | ||
| 440 | * in a list can be marked. The same space is used for | ||
| 441 | * marking and chaining tetrahedra in lists. | ||
| 442 | * A tetrahedron can be in the following states: | ||
| 443 | * - in list | ||
| 444 | * - not in list and marked | ||
| 445 | * - not in list and not marked | ||
| 446 | * \param[in] t the index of the tetrahedron | ||
| 447 | * \retval true if tetrahedron \p t belongs to a linked list | ||
| 448 | * \retval false otherwise | ||
| 449 | */ | ||
| 450 | 31104 | bool tet_is_in_list(index_t t) const { | |
| 451 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 31104 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
31104 | geo_debug_assert(t < max_t()); |
| 452 | 31104 | return (cell_next_[t] & NOT_IN_LIST_BIT) == 0; | |
| 453 | } | ||
| 454 | |||
| 455 | /** | ||
| 456 | * \brief Gets the index of a successor of a tetrahedron. | ||
| 457 | * \details Tetrahedra can be linked, it is used to manage | ||
| 458 | * both the free list that recycles deleted tetrahedra. | ||
| 459 | * \param[in] t the index of the tetrahedron | ||
| 460 | * \retval END_OF_LIST if the end of the list is reached | ||
| 461 | * \retval the index of the successor of | ||
| 462 | * tetrahedron \t otherwise | ||
| 463 | * \pre tet_is_in_list(t) | ||
| 464 | */ | ||
| 465 | 3272 | index_t tet_next(index_t t) const { | |
| 466 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 3272 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
3272 | geo_debug_assert(t < max_t()); |
| 467 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 3272 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
3272 | geo_debug_assert(tet_is_in_list(t)); |
| 468 | 3272 | return cell_next_[t]; | |
| 469 | } | ||
| 470 | |||
| 471 | /** | ||
| 472 | * \brief Adds a tetrahedron to a linked list. | ||
| 473 | * \details Tetrahedra can be linked, it is used to manage | ||
| 474 | * the free list that recycles deleted tetrahedra. | ||
| 475 | * \param[in] t the index of the tetrahedron | ||
| 476 | * \param[in,out] first first item of the list or END_OF_LIST if | ||
| 477 | * the list is empty | ||
| 478 | * \param[in,out] last last item of the list or END_OF_LIST if | ||
| 479 | * the list is empty | ||
| 480 | */ | ||
| 481 | 3313 | void add_tet_to_list(index_t t, index_t& first, index_t& last) { | |
| 482 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 3313 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
3313 | geo_debug_assert(t < max_t()); |
| 483 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 3313 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
3313 | geo_debug_assert(!tet_is_in_list(t)); |
| 484 |
2/2✓ Branch 0 taken 221 times.
✓ Branch 1 taken 3092 times.
|
3313 | if(last == END_OF_LIST) { |
| 485 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 221 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
221 | geo_debug_assert(first == END_OF_LIST); |
| 486 | 221 | first = last = t; | |
| 487 | 221 | cell_next_[t] = END_OF_LIST; | |
| 488 | } else { | ||
| 489 | 3092 | cell_next_[t] = first; | |
| 490 | 3092 | first = t; | |
| 491 | } | ||
| 492 | 3313 | } | |
| 493 | |||
| 494 | /** | ||
| 495 | * \brief Removes a tetrahedron from the linked list it | ||
| 496 | * belongs to. | ||
| 497 | * \details Tetrahedra can be linked, it is used to manage | ||
| 498 | * both the free list that recycles deleted tetrahedra and | ||
| 499 | * the list of tetrahedra in conflict with the inserted | ||
| 500 | * point. | ||
| 501 | * \param[in] t the index of the tetrahedron | ||
| 502 | */ | ||
| 503 | 3272 | void remove_tet_from_list(index_t t) { | |
| 504 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 3272 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
3272 | geo_debug_assert(t < max_t()); |
| 505 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 3272 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
3272 | geo_debug_assert(tet_is_in_list(t)); |
| 506 | 3272 | cell_next_[t] = NOT_IN_LIST; | |
| 507 | 3272 | } | |
| 508 | |||
| 509 | /** | ||
| 510 | * \brief Symbolic value for a vertex of a | ||
| 511 | * tetrahedron that indicates a virtual tetrahedron. | ||
| 512 | * \details The three other vertices then correspond to a | ||
| 513 | * facet on the convex hull of the points. | ||
| 514 | */ | ||
| 515 | static constexpr index_t VERTEX_AT_INFINITY = NO_INDEX; | ||
| 516 | |||
| 517 | /** | ||
| 518 | * \brief Tests whether a given tetrahedron | ||
| 519 | * is a finite one. | ||
| 520 | * \details Infinite tetrahedra are the ones | ||
| 521 | * that are incident to the infinite vertex | ||
| 522 | * (index -1) | ||
| 523 | * \param[in] t the index of the tetrahedron | ||
| 524 | * \retval true if \p t is finite | ||
| 525 | * \retval false otherwise | ||
| 526 | */ | ||
| 527 | 1145 | bool tet_is_finite(index_t t) const { | |
| 528 | return | ||
| 529 | 1145 | cell_to_v_store_[4 * t] != NO_INDEX && | |
| 530 |
2/2✓ Branch 1 taken 996 times.
✓ Branch 2 taken 148 times.
|
1144 | cell_to_v_store_[4 * t + 1] != NO_INDEX && |
| 531 |
4/4✓ Branch 0 taken 1144 times.
✓ Branch 1 taken 1 times.
✓ Branch 3 taken 914 times.
✓ Branch 4 taken 82 times.
|
3203 | cell_to_v_store_[4 * t + 2] != NO_INDEX && |
| 532 |
2/2✓ Branch 1 taken 779 times.
✓ Branch 2 taken 135 times.
|
2059 | cell_to_v_store_[4 * t + 3] != NO_INDEX; |
| 533 | } | ||
| 534 | |||
| 535 | /** | ||
| 536 | * \brief Tests whether a tetrahedron is | ||
| 537 | * a real one. | ||
| 538 | * \details Real tetrahedra are incident to | ||
| 539 | * four user-specified vertices (there are also | ||
| 540 | * virtual tetrahedra that are incident to the | ||
| 541 | * vertex at infinity, with index -1) | ||
| 542 | * \param[in] t index of the tetrahedron | ||
| 543 | * \retval true if tetrahedron \p t is a real one | ||
| 544 | * \retval false otherwise | ||
| 545 | */ | ||
| 546 | 1186 | bool tet_is_real(index_t t) const { | |
| 547 |
4/4✓ Branch 1 taken 1145 times.
✓ Branch 2 taken 41 times.
✓ Branch 4 taken 779 times.
✓ Branch 5 taken 366 times.
|
1186 | return !tet_is_free(t) && tet_is_finite(t); |
| 548 | } | ||
| 549 | |||
| 550 | /** | ||
| 551 | * \brief Tests whether a tetrahedron is | ||
| 552 | * a virtual one. | ||
| 553 | * \details Virtual tetrahedra are tetrahedra | ||
| 554 | * incident to the vertex at infinity. | ||
| 555 | * \param[in] t index of the tetrahedron | ||
| 556 | * \retval true if tetrahedron \p t is virtual | ||
| 557 | * \retval false otherwise | ||
| 558 | */ | ||
| 559 | 2946 | bool tet_is_virtual(index_t t) const { | |
| 560 | return | ||
| 561 |
2/2✓ Branch 1 taken 2824 times.
✓ Branch 2 taken 122 times.
|
5770 | !tet_is_free(t) && ( |
| 562 |
2/2✓ Branch 1 taken 2818 times.
✓ Branch 2 taken 6 times.
|
2824 | cell_to_v_store_[4 * t] == VERTEX_AT_INFINITY || |
| 563 |
2/2✓ Branch 1 taken 2655 times.
✓ Branch 2 taken 163 times.
|
2818 | cell_to_v_store_[4 * t + 1] == VERTEX_AT_INFINITY || |
| 564 |
2/2✓ Branch 1 taken 2561 times.
✓ Branch 2 taken 94 times.
|
2655 | cell_to_v_store_[4 * t + 2] == VERTEX_AT_INFINITY || |
| 565 |
2/2✓ Branch 1 taken 121 times.
✓ Branch 2 taken 2440 times.
|
5507 | cell_to_v_store_[4 * t + 3] == VERTEX_AT_INFINITY) ; |
| 566 | } | ||
| 567 | |||
| 568 | /** | ||
| 569 | * \brief Tests whether a tetrahedron is | ||
| 570 | * in the free list. | ||
| 571 | * \details Deleted tetrahedra are recycled | ||
| 572 | * in a free list. | ||
| 573 | * \param[in] t index of the tetrahedron | ||
| 574 | * \retval true if tetrahedron \p t is in | ||
| 575 | * the free list | ||
| 576 | * \retval false otherwise | ||
| 577 | */ | ||
| 578 | 4134 | bool tet_is_free(index_t t) const { | |
| 579 | 4134 | return tet_is_in_list(t); | |
| 580 | } | ||
| 581 | |||
| 582 | /** | ||
| 583 | * \brief Creates a new tetrahedron. | ||
| 584 | * \details Uses either a tetrahedron recycled | ||
| 585 | * from the free list, or creates a new one by | ||
| 586 | * expanding the two indices arrays. | ||
| 587 | * \return the index of the newly created tetrahedron | ||
| 588 | */ | ||
| 589 | 4458 | index_t new_tetrahedron() { | |
| 590 | index_t result; | ||
| 591 |
2/2✓ Branch 0 taken 1186 times.
✓ Branch 1 taken 3272 times.
|
4458 | if(first_free_ == END_OF_LIST) { |
| 592 | 1186 | cell_to_v_store_.resize( | |
| 593 | 1186 | cell_to_v_store_.size() + 4, NO_INDEX | |
| 594 | ); | ||
| 595 | 1186 | cell_to_cell_store_.resize( | |
| 596 | 1186 | cell_to_cell_store_.size() + 4, NO_INDEX | |
| 597 | ); | ||
| 598 | // index_t(NOT_IN_LIST) is necessary, else with | ||
| 599 | // NOT_IN_LIST alone the compiler tries to generate a | ||
| 600 | // reference to NOT_IN_LIST resulting in a link error. | ||
| 601 |
1/2✓ Branch 1 taken 1186 times.
✗ Branch 2 not taken.
|
1186 | cell_next_.push_back(index_t(NOT_IN_LIST)); |
| 602 | 1186 | result = max_t() - 1; | |
| 603 | } else { | ||
| 604 | 3272 | result = first_free_; | |
| 605 | 3272 | first_free_ = tet_next(first_free_); | |
| 606 | 3272 | remove_tet_from_list(result); | |
| 607 | } | ||
| 608 | |||
| 609 | 4458 | cell_to_cell_store_[4 * result] = NO_INDEX; | |
| 610 | 4458 | cell_to_cell_store_[4 * result + 1] = NO_INDEX; | |
| 611 | 4458 | cell_to_cell_store_[4 * result + 2] = NO_INDEX; | |
| 612 | 4458 | cell_to_cell_store_[4 * result + 3] = NO_INDEX; | |
| 613 | |||
| 614 | 4458 | return result; | |
| 615 | } | ||
| 616 | |||
| 617 | /** | ||
| 618 | * \brief Creates a new tetrahedron. | ||
| 619 | * \details Sets the vertices. Adjacent tetrahedra index are | ||
| 620 | * left uninitialized. Uses either a tetrahedron recycled | ||
| 621 | * from the free list, or creates a new one by | ||
| 622 | * expanding the two indices arrays. | ||
| 623 | * \param[in] v1 index of the first vertex | ||
| 624 | * \param[in] v2 index of the second vertex | ||
| 625 | * \param[in] v3 index of the third vertex | ||
| 626 | * \param[in] v4 index of the fourth vertex | ||
| 627 | * \return the index of the newly created tetrahedron | ||
| 628 | */ | ||
| 629 | 4458 | index_t new_tetrahedron( | |
| 630 | index_t v1, index_t v2, | ||
| 631 | index_t v3, index_t v4 | ||
| 632 | ) { | ||
| 633 | 4458 | index_t result = new_tetrahedron(); | |
| 634 | 4458 | cell_to_v_store_[4 * result] = v1; | |
| 635 | 4458 | cell_to_v_store_[4 * result + 1] = v2; | |
| 636 | 4458 | cell_to_v_store_[4 * result + 2] = v3; | |
| 637 | 4458 | cell_to_v_store_[4 * result + 3] = v4; | |
| 638 | 4458 | return result; | |
| 639 | } | ||
| 640 | |||
| 641 | /** | ||
| 642 | * \brief Generates a unique stamp for marking tets. | ||
| 643 | * \details Storage is shared for list-chaining and stamp-marking | ||
| 644 | * (both are mutually exclusive), therefore the stamp has | ||
| 645 | * the NOT_IN_LIST_BIT set. | ||
| 646 | * \param[in] stamp the unique stamp for marking tets | ||
| 647 | */ | ||
| 648 | 221 | void set_tet_mark_stamp(index_t stamp) { | |
| 649 | 221 | cur_stamp_ = (stamp | NOT_IN_LIST_BIT); | |
| 650 | 221 | } | |
| 651 | |||
| 652 | /** | ||
| 653 | * \brief Tests whether a tetrahedron is marked. | ||
| 654 | * \details A tetrahedron is marked whenever it is | ||
| 655 | * detected as non-conflict. The index of the | ||
| 656 | * point being inserted is used as a time-stamp | ||
| 657 | * for marking tetrahedra. The same space is used | ||
| 658 | * for marking and for chaining the conflict list. | ||
| 659 | * A tetrahedron can be in the following states: | ||
| 660 | * - in list | ||
| 661 | * - not in list and marked | ||
| 662 | * - not in list and not marked | ||
| 663 | * \param[in] t index of the tetrahedron | ||
| 664 | * \retval true if tetrahedron \p t is marked | ||
| 665 | * \retval false otherwise | ||
| 666 | */ | ||
| 667 | 8290 | bool tet_is_marked(index_t t) const { | |
| 668 | 8290 | return cell_next_[t] == cur_stamp_; | |
| 669 | } | ||
| 670 | |||
| 671 | /** | ||
| 672 | * \brief Marks a tetrahedron. | ||
| 673 | * \details A tetrahedron is marked whenever it is | ||
| 674 | * detected as non-conflict. The same space is used | ||
| 675 | * for marking and for chaining the conflict list. | ||
| 676 | * The index of the point being inserted is used as a | ||
| 677 | * time-stamp for marking tetrahedra. | ||
| 678 | * A tetrahedron can be in the following states: | ||
| 679 | * - in list | ||
| 680 | * - not in list and marked | ||
| 681 | * - not in list and not marked | ||
| 682 | * \param[in] t index of the tetrahedron to be marked | ||
| 683 | */ | ||
| 684 | 3952 | void mark_tet(index_t t) { | |
| 685 | 3952 | cell_next_[t] = cur_stamp_; | |
| 686 | 3952 | } | |
| 687 | |||
| 688 | /********* Combinatorics ******************************************/ | ||
| 689 | |||
| 690 | /** | ||
| 691 | * \brief Returns the local index of a vertex by | ||
| 692 | * facet and by local vertex index in the facet. | ||
| 693 | * \details | ||
| 694 | * tet facet vertex is such that the tetrahedron | ||
| 695 | * formed with: | ||
| 696 | * - vertex lv | ||
| 697 | * - tet_facet_vertex(lv,0) | ||
| 698 | * - tet_facet_vertex(lv,1) | ||
| 699 | * - tet_facet_vertex(lv,2) | ||
| 700 | * has the same orientation as the original tetrahedron for | ||
| 701 | * any vertex lv. | ||
| 702 | * \param[in] f local facet index, in (0,1,2,3) | ||
| 703 | * \param[in] v local vertex index, in (0,1,2) | ||
| 704 | * \return the local tetrahedron vertex index of | ||
| 705 | * vertex \p v in facet \p f | ||
| 706 | */ | ||
| 707 | 13392 | static index_t tet_facet_vertex(index_t f, index_t v) { | |
| 708 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 13392 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
13392 | geo_debug_assert(f < 4); |
| 709 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 13392 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
13392 | geo_debug_assert(v < 3); |
| 710 | 13392 | return index_t(tet_facet_vertex_[f][v]); | |
| 711 | } | ||
| 712 | |||
| 713 | /** | ||
| 714 | * \brief Gets the index of a vertex of a tetrahedron | ||
| 715 | * \param[in] t index of the tetrahedron | ||
| 716 | * \param[in] lv local vertex (0,1,2 or 3) index in \p t | ||
| 717 | * \return the global index of the \p lv%th vertex of tetrahedron \p t | ||
| 718 | * or -1 if the vertex is at infinity | ||
| 719 | */ | ||
| 720 | 48791 | index_t tet_vertex(index_t t, index_t lv) const { | |
| 721 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 48791 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
48791 | geo_debug_assert(t < max_t()); |
| 722 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 48791 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
48791 | geo_debug_assert(lv < 4); |
| 723 | 48791 | return cell_to_v_store_[4 * t + lv]; | |
| 724 | } | ||
| 725 | |||
| 726 | /** | ||
| 727 | * \brief Finds the index of the vertex in a tetrahedron. | ||
| 728 | * \param[in] t the tetrahedron | ||
| 729 | * \param[in] v the vertex | ||
| 730 | * \return iv such that tet_vertex(t,v)==iv | ||
| 731 | * \pre \p t is incident to \p v | ||
| 732 | */ | ||
| 733 | 750 | index_t find_tet_vertex(index_t t, index_t v) const { | |
| 734 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 750 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
750 | geo_debug_assert(t < max_t()); |
| 735 | // Find local index of v in tetrahedron t vertices. | ||
| 736 | 750 | const index_t* T = &(cell_to_v_store_[4 * t]); | |
| 737 | 750 | return find_4(T,v); | |
| 738 | } | ||
| 739 | |||
| 740 | /** | ||
| 741 | * \brief Gets the index of a vertex of a tetrahedron | ||
| 742 | * \param[in] t index of the tetrahedron | ||
| 743 | * \param[in] lv local vertex (0,1,2 or 3) index in \p t | ||
| 744 | * \return the global index of the \p lv%th vertex of tetrahedron \p t | ||
| 745 | * \pre Vertex \p lv of tetrahedron \p t is not at infinity | ||
| 746 | */ | ||
| 747 | 7072 | index_t finite_tet_vertex(index_t t, index_t lv) const { | |
| 748 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 7072 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
7072 | geo_debug_assert(t < max_t()); |
| 749 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 7072 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
7072 | geo_debug_assert(lv < 4); |
| 750 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 7072 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
7072 | geo_debug_assert(cell_to_v_store_[4 * t + lv] != NO_INDEX); |
| 751 | 7072 | return cell_to_v_store_[4 * t + lv]; | |
| 752 | } | ||
| 753 | |||
| 754 | /** | ||
| 755 | * \brief Sets a tetrahedron-to-vertex adjacency. | ||
| 756 | * \param[in] t index of the tetrahedron | ||
| 757 | * \param[in] lv local vertex index (0,1,2 or 3) in \p t | ||
| 758 | * \param[in] v global index of the vertex | ||
| 759 | */ | ||
| 760 | 500 | void set_tet_vertex(index_t t, index_t lv, index_t v) { | |
| 761 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 500 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
500 | geo_debug_assert(t < max_t()); |
| 762 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 500 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
500 | geo_debug_assert(lv < 4); |
| 763 | 500 | cell_to_v_store_[4 * t + lv] = v; | |
| 764 | 500 | } | |
| 765 | |||
| 766 | /** | ||
| 767 | * \brief Gets the index of a tetrahedron adjacent to another one. | ||
| 768 | * \param[in] t index of the tetrahedron | ||
| 769 | * \param[in] lf local facet (0,1,2 or 3) index in \p t | ||
| 770 | * \return the tetrahedron adjacent to \p t accorss facet \p lf | ||
| 771 | */ | ||
| 772 | 42980 | index_t tet_adjacent(index_t t, index_t lf) const { | |
| 773 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 42980 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
42980 | geo_debug_assert(t < max_t()); |
| 774 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 42980 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
42980 | geo_debug_assert(lf < 4); |
| 775 | 42980 | index_t result = cell_to_cell_store_[4 * t + lf]; | |
| 776 | 42980 | return result; | |
| 777 | } | ||
| 778 | |||
| 779 | /** | ||
| 780 | * \brief Sets a tetrahedron-to-tetrahedron adjacency. | ||
| 781 | * \param[in] t1 index of the first tetrahedron | ||
| 782 | * \param[in] lf1 local facet index (0,1,2 or 3) in t1 | ||
| 783 | * \param[in] t2 index of the tetrahedron | ||
| 784 | * adjacent to \p t1 accros \p lf1 | ||
| 785 | */ | ||
| 786 | 22280 | void set_tet_adjacent(index_t t1, index_t lf1, index_t t2) { | |
| 787 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 22280 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
22280 | geo_debug_assert(t1 < max_t()); |
| 788 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 22280 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
22280 | geo_debug_assert(t2 < max_t()); |
| 789 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 22280 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
22280 | geo_debug_assert(lf1 < 4); |
| 790 | 22280 | cell_to_cell_store_[4 * t1 + lf1] = t2; | |
| 791 | 22280 | } | |
| 792 | |||
| 793 | /** | ||
| 794 | * \brief Finds the index of the facet accros which t1 is | ||
| 795 | * adjacent to t2_in. | ||
| 796 | * \param[in] t1 first tetrahedron | ||
| 797 | * \param[in] t2 second tetrahedron | ||
| 798 | * \return f such that tet_adjacent(t1,f)==t2 | ||
| 799 | * \pre \p t1 and \p t2 are adjacent | ||
| 800 | */ | ||
| 801 | 4448 | index_t find_tet_adjacent(index_t t1, index_t t2) const { | |
| 802 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 4448 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
4448 | geo_debug_assert(t1 < max_t()); |
| 803 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 4448 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
4448 | geo_debug_assert(t2 < max_t()); |
| 804 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 4448 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
4448 | geo_debug_assert(t1 != t2); |
| 805 | |||
| 806 | // Find local index of t2 in tetrahedron t1 adajcent tets. | ||
| 807 | 4448 | const index_t* T = &(cell_to_cell_store_[4 * t1]); | |
| 808 | 4448 | index_t result = find_4(T,t2); | |
| 809 | |||
| 810 | // Sanity check: make sure that t1 is adjacent to t2 | ||
| 811 | // only once! | ||
| 812 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 4448 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
4448 | geo_debug_assert(tet_adjacent(t1,(result+1)%4) != t2); |
| 813 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 4448 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
4448 | geo_debug_assert(tet_adjacent(t1,(result+2)%4) != t2); |
| 814 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 4448 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
4448 | geo_debug_assert(tet_adjacent(t1,(result+3)%4) != t2); |
| 815 | 4448 | return result; | |
| 816 | } | ||
| 817 | |||
| 818 | /** | ||
| 819 | * \brief Sets the vertices and adjacent tetrahedra of | ||
| 820 | * a tetrahedron. | ||
| 821 | * \param[in] t index of the tetrahedron | ||
| 822 | * \param[in] v0 index of the first vertex | ||
| 823 | * \param[in] v1 index of the second vertex | ||
| 824 | * \param[in] v2 index of the third vertex | ||
| 825 | * \param[in] v3 index of the fourth vertex | ||
| 826 | * \param[in] a0 index of the adjacent tetrahedron opposite to \p v0 | ||
| 827 | * \param[in] a1 index of the adjacent tetrahedron opposite to \p v1 | ||
| 828 | * \param[in] a2 index of the adjacent tetrahedron opposite to \p v2 | ||
| 829 | * \param[in] a3 index of the adjacent tetrahedron opposite to \p v3 | ||
| 830 | */ | ||
| 831 | void set_tet( | ||
| 832 | index_t t, | ||
| 833 | index_t v0, index_t v1, index_t v2, index_t v3, | ||
| 834 | index_t a0, index_t a1, index_t a2, index_t a3 | ||
| 835 | ) { | ||
| 836 | geo_debug_assert(t < max_t()); | ||
| 837 | cell_to_v_store_[4 * t] = v0; | ||
| 838 | cell_to_v_store_[4 * t + 1] = v1; | ||
| 839 | cell_to_v_store_[4 * t + 2] = v2; | ||
| 840 | cell_to_v_store_[4 * t + 3] = v3; | ||
| 841 | cell_to_cell_store_[4 * t] = a0; | ||
| 842 | cell_to_cell_store_[4 * t + 1] = a1; | ||
| 843 | cell_to_cell_store_[4 * t + 2] = a2; | ||
| 844 | cell_to_cell_store_[4 * t + 3] = a3; | ||
| 845 | } | ||
| 846 | |||
| 847 | /****** Combinatorics - traversals ************************/ | ||
| 848 | |||
| 849 | /** | ||
| 850 | * Gets the local facet index incident to an | ||
| 851 | * oriented halfedge. | ||
| 852 | * \param[in] t index of the tetrahedron | ||
| 853 | * \param[in] v1 global index of the first extremity | ||
| 854 | * \param[in] v2 global index of the second extremity | ||
| 855 | * \return the local index of the facet incident to | ||
| 856 | * the oriented edge \p v1, \p v2. | ||
| 857 | */ | ||
| 858 | 1143 | index_t get_facet_by_halfedge(index_t t, index_t v1, index_t v2) const { | |
| 859 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 1143 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
1143 | geo_debug_assert(t < max_t()); |
| 860 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 1143 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1143 | geo_debug_assert(v1 != v2); |
| 861 | // Find local index of v1 and v2 in tetrahedron t | ||
| 862 | 1143 | const index_t* T = &(cell_to_v_store_[4 * t]); | |
| 863 | 1143 | index_t lv1 = find_4(T,v1); | |
| 864 | 1143 | index_t lv2 = find_4(T,v2); | |
| 865 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 1143 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1143 | geo_debug_assert(lv1 != lv2); |
| 866 | 1143 | return index_t(halfedge_facet_[lv1][lv2]); | |
| 867 | } | ||
| 868 | |||
| 869 | /** | ||
| 870 | * Gets the local facet indices incident to an | ||
| 871 | * oriented halfedge. | ||
| 872 | * \param[in] t index of the tetrahedron | ||
| 873 | * \param[in] v1 global index of the first extremity | ||
| 874 | * \param[in] v2 global index of the second extremity | ||
| 875 | * \param[out] f12 the local index of the facet | ||
| 876 | * indicent to the halfedge [v1,v2] | ||
| 877 | * \param[out] f21 the local index of the facet | ||
| 878 | * indicent to the halfedge [v2,v1] | ||
| 879 | */ | ||
| 880 | 750 | void get_facets_by_halfedge( | |
| 881 | index_t t, index_t v1, index_t v2, | ||
| 882 | index_t& f12, index_t& f21 | ||
| 883 | ) const { | ||
| 884 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 750 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
750 | geo_debug_assert(t < max_t()); |
| 885 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 750 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
750 | geo_debug_assert(v1 != v2); |
| 886 | |||
| 887 | // Find local index of v1 and v2 in tetrahedron t | ||
| 888 | // The following expression is 10% faster than using | ||
| 889 | // if() statements (multiply by boolean result of test). | ||
| 890 | // Thank to Laurent Alonso for this idea. | ||
| 891 | 750 | const index_t* T = &(cell_to_v_store_[4 * t]); | |
| 892 | |||
| 893 | 750 | index_t lv1 = index_t( | |
| 894 |
4/4✓ Branch 0 taken 188 times.
✓ Branch 1 taken 562 times.
✓ Branch 2 taken 194 times.
✓ Branch 3 taken 556 times.
|
750 | (T[1] == v1) | ((T[2] == v1) * 2) | ((T[3] == v1) * 3) |
| 895 | ); | ||
| 896 | |||
| 897 | 750 | index_t lv2 = index_t( | |
| 898 |
4/4✓ Branch 0 taken 189 times.
✓ Branch 1 taken 561 times.
✓ Branch 2 taken 147 times.
✓ Branch 3 taken 603 times.
|
750 | (T[1] == v2) | ((T[2] == v2) * 2) | ((T[3] == v2) * 3) |
| 899 | ); | ||
| 900 | |||
| 901 |
3/8✓ Branch 0 taken 163 times.
✓ Branch 1 taken 587 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 163 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
750 | geo_debug_assert(lv1 != 0 || T[0] == v1); |
| 902 |
3/8✓ Branch 0 taken 277 times.
✓ Branch 1 taken 473 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 277 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
750 | geo_debug_assert(lv2 != 0 || T[0] == v2); |
| 903 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 750 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
750 | geo_debug_assert(lv1 != lv2); |
| 904 | |||
| 905 | 750 | f12 = index_t(halfedge_facet_[lv1][lv2]); | |
| 906 | 750 | f21 = index_t(halfedge_facet_[lv2][lv1]); | |
| 907 | 750 | } | |
| 908 | |||
| 909 | |||
| 910 | /** | ||
| 911 | * \brief Gets the next tetrahedron around an oriented edge of | ||
| 912 | * a tetrahedron. | ||
| 913 | * \param[in,out] t the tetrahedron | ||
| 914 | * \param[in] v1 global index of the first extremity of the edge | ||
| 915 | * \param[in] v2 global index of the second extremity of the edge | ||
| 916 | * \return the next tetrahedron from \p t around the oriented edge | ||
| 917 | * (\p v1 \p v2). | ||
| 918 | */ | ||
| 919 | index_t next_around_halfedge(index_t& t, index_t v1, index_t v2) const { | ||
| 920 | return (index_t)tet_adjacent( | ||
| 921 | t, get_facet_by_halfedge(t, v1, v2) | ||
| 922 | ); | ||
| 923 | } | ||
| 924 | |||
| 925 | /****** Predicates **********************************************/ | ||
| 926 | |||
| 927 | /** | ||
| 928 | * \brief Tests whether a given tetrahedron is in conflict with | ||
| 929 | * a given 3d point. | ||
| 930 | * \details A real tetrahedron is in conflict with a point whenever | ||
| 931 | * the point is contained by its circumscribed sphere, and a | ||
| 932 | * virtual tetrahedron is in conflict with a point whenever the | ||
| 933 | * tetrahedron formed by its real face and with the point has | ||
| 934 | * positive orientation. | ||
| 935 | * \param[in] t the index of the tetrahedron | ||
| 936 | * \param[in] p a pointer to the coordinates of the point | ||
| 937 | * \retval true if point \p p is in conflict with tetrahedron \p t | ||
| 938 | * \retval false otherwise | ||
| 939 | */ | ||
| 940 | 7668 | bool tet_is_conflict(index_t t, const double* p) const { | |
| 941 | |||
| 942 | // Lookup tetrahedron vertices | ||
| 943 | const double* pv[4]; | ||
| 944 |
2/2✓ Branch 0 taken 30672 times.
✓ Branch 1 taken 7668 times.
|
38340 | for(index_t i=0; i<4; ++i) { |
| 945 |
1/2✓ Branch 1 taken 30672 times.
✗ Branch 2 not taken.
|
30672 | index_t v = tet_vertex(t,i); |
| 946 |
3/4✓ Branch 0 taken 28710 times.
✓ Branch 1 taken 1962 times.
✓ Branch 3 taken 28710 times.
✗ Branch 4 not taken.
|
30672 | pv[i] = (v == NO_INDEX) ? nullptr : vertex_ptr(v); |
| 947 | } | ||
| 948 | |||
| 949 | // Check for virtual tetrahedra (then in_sphere() | ||
| 950 | // is replaced with orient3d()) | ||
| 951 |
2/2✓ Branch 0 taken 28588 times.
✓ Branch 1 taken 5706 times.
|
34294 | for(index_t lf = 0; lf < 4; ++lf) { |
| 952 | |||
| 953 |
2/2✓ Branch 0 taken 1962 times.
✓ Branch 1 taken 26626 times.
|
28588 | if(pv[lf] == nullptr) { |
| 954 | |||
| 955 | // Facet of a virtual tetrahedron opposite to | ||
| 956 | // infinite vertex corresponds to | ||
| 957 | // the triangle on the convex hull of the points. | ||
| 958 | // Orientation is obtained by replacing vertex lf | ||
| 959 | // with p. | ||
| 960 | 1962 | pv[lf] = p; | |
| 961 |
1/2✓ Branch 1 taken 1962 times.
✗ Branch 2 not taken.
|
1962 | Sign sign = PCK::orient_3d(pv[0],pv[1],pv[2],pv[3]); |
| 962 | |||
| 963 |
2/2✓ Branch 0 taken 203 times.
✓ Branch 1 taken 1759 times.
|
1962 | if(sign > 0) { |
| 964 | 203 | return true; | |
| 965 | } | ||
| 966 | |||
| 967 |
2/2✓ Branch 0 taken 791 times.
✓ Branch 1 taken 968 times.
|
1759 | if(sign < 0) { |
| 968 | 791 | return false; | |
| 969 | } | ||
| 970 | |||
| 971 | // If sign is zero, we check the real tetrahedron | ||
| 972 | // adjacent to the facet on the convex hull. | ||
| 973 |
2/8✓ Branch 1 taken 968 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 968 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
968 | geo_debug_assert(tet_adjacent(t, lf) != NO_INDEX); |
| 974 |
1/2✓ Branch 1 taken 968 times.
✗ Branch 2 not taken.
|
968 | index_t t2 = tet_adjacent(t, lf); |
| 975 |
2/8✓ Branch 1 taken 968 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 968 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
968 | geo_debug_assert(!tet_is_virtual(t2)); |
| 976 | |||
| 977 | // If t2 is already chained in the conflict list, | ||
| 978 | // then it is conflict | ||
| 979 |
3/4✓ Branch 1 taken 968 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 122 times.
✓ Branch 4 taken 846 times.
|
968 | if(tet_is_in_list(t2)) { |
| 980 | 122 | return true; | |
| 981 | } | ||
| 982 | |||
| 983 | // If t2 is marked, then it is not in conflict. | ||
| 984 |
3/4✓ Branch 1 taken 846 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 126 times.
✓ Branch 4 taken 720 times.
|
846 | if(tet_is_marked(t2)) { |
| 985 | 126 | return false; | |
| 986 | } | ||
| 987 | |||
| 988 |
1/2✓ Branch 1 taken 720 times.
✗ Branch 2 not taken.
|
720 | return tet_is_conflict(t2, p); |
| 989 | } | ||
| 990 | } | ||
| 991 | |||
| 992 | // If the tetrahedron is a finite one, it is in conflict | ||
| 993 | // if its circumscribed sphere contains the point (this is | ||
| 994 | // the standard case). | ||
| 995 | |||
| 996 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5706 times.
|
5706 | if(weighted_) { |
| 997 | ✗ | double h0 = heights_[finite_tet_vertex(t, 0)]; | |
| 998 | ✗ | double h1 = heights_[finite_tet_vertex(t, 1)]; | |
| 999 | ✗ | double h2 = heights_[finite_tet_vertex(t, 2)]; | |
| 1000 | ✗ | double h3 = heights_[finite_tet_vertex(t, 3)]; | |
| 1001 | index_t pindex = index_t( | ||
| 1002 | ✗ | (p - vertex_ptr(0)) / int(vertex_stride_) | |
| 1003 | ✗ | ); | |
| 1004 | ✗ | double h = heights_[pindex]; | |
| 1005 | ✗ | return (PCK::orient_3dlifted_SOS( | |
| 1006 | pv[0],pv[1],pv[2],pv[3],p,h0,h1,h2,h3,h | ||
| 1007 | ✗ | ) > 0) ; | |
| 1008 | } | ||
| 1009 | |||
| 1010 |
1/2✓ Branch 1 taken 5706 times.
✗ Branch 2 not taken.
|
5706 | return (PCK::in_sphere_3d_SOS(pv[0], pv[1], pv[2], pv[3], p) > 0); |
| 1011 | } | ||
| 1012 | |||
| 1013 | protected: | ||
| 1014 | |||
| 1015 | /** | ||
| 1016 | * \brief Finds the index of an integer in an array of four integers. | ||
| 1017 | * \param[in] T a const pointer to an array of four integers | ||
| 1018 | * \param[in] v the integer to retrieve in \p T | ||
| 1019 | * \return the index (0,1,2 or 3) of \p v in \p T | ||
| 1020 | * \pre The four entries of \p T are different and one of them is | ||
| 1021 | * equal to \p v. | ||
| 1022 | */ | ||
| 1023 | 7484 | static index_t find_4(const index_t* T, index_t v) { | |
| 1024 | // The following expression is 10% faster than using | ||
| 1025 | // if() statements. This uses the C++ norm, that | ||
| 1026 | // ensures that the 'true' boolean value converted to | ||
| 1027 | // an int is always 1. With most compilers, this avoids | ||
| 1028 | // generating branching instructions. | ||
| 1029 | // Thank to Laurent Alonso for this idea. | ||
| 1030 | // Note: Laurent also has this version: | ||
| 1031 | // (T[0] != v)+(T[2]==v)+2*(T[3]==v) | ||
| 1032 | // that avoids a *3 multiply, but it is not faster in | ||
| 1033 | // practice. | ||
| 1034 | 7484 | index_t result = index_t( | |
| 1035 |
4/4✓ Branch 0 taken 2121 times.
✓ Branch 1 taken 5363 times.
✓ Branch 2 taken 1896 times.
✓ Branch 3 taken 5588 times.
|
7484 | (T[1] == v) | ((T[2] == v) * 2) | ((T[3] == v) * 3) |
| 1036 | ); | ||
| 1037 | // Sanity check, important if it was T[0], not explicitly | ||
| 1038 | // tested (detects input that does not meet the precondition). | ||
| 1039 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 7484 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
7484 | geo_debug_assert(T[result] == v); |
| 1040 | 7484 | return result; | |
| 1041 | } | ||
| 1042 | |||
| 1043 | /** | ||
| 1044 | * \brief Delaunay3d destructor | ||
| 1045 | */ | ||
| 1046 | ~Delaunay3d() override; | ||
| 1047 | |||
| 1048 | /** | ||
| 1049 | * \brief For debugging purposes, displays a tetrahedron. | ||
| 1050 | * \param[in] t index of the tetrahedron to display. | ||
| 1051 | */ | ||
| 1052 | void show_tet(index_t t) const; | ||
| 1053 | |||
| 1054 | /** | ||
| 1055 | * \brief For debugging purposes, displays a tetrahedron adjacency. | ||
| 1056 | * \param[in] t index of the tetrahedron to display. | ||
| 1057 | * \param[in] lf local index (0,1,2 or 3) of the tetrahedron | ||
| 1058 | * facet adjacenty to display. | ||
| 1059 | */ | ||
| 1060 | void show_tet_adjacent(index_t t, index_t lf) const; | ||
| 1061 | |||
| 1062 | /** | ||
| 1063 | * \brief For debugging purposes, displays a tetrahedron. | ||
| 1064 | * \param[in] first index of the first tetrahedron in the list | ||
| 1065 | * \param[in] list_name name of the list, will be displayed as well | ||
| 1066 | */ | ||
| 1067 | void show_list(index_t first, const std::string& list_name) const; | ||
| 1068 | |||
| 1069 | /** | ||
| 1070 | * \brief For debugging purposes, tests some combinatorial properties. | ||
| 1071 | */ | ||
| 1072 | void check_combinatorics(bool verbose = false) const; | ||
| 1073 | |||
| 1074 | /** | ||
| 1075 | * \brief For debugging purposes, test some geometrical properties. | ||
| 1076 | */ | ||
| 1077 | void check_geometry(bool verbose = false) const; | ||
| 1078 | |||
| 1079 | private: | ||
| 1080 | vector<index_t> cell_to_v_store_; | ||
| 1081 | vector<index_t> cell_to_cell_store_; | ||
| 1082 | vector<index_t> cell_next_; | ||
| 1083 | vector<index_t> reorder_; | ||
| 1084 | index_t cur_stamp_; // used for marking | ||
| 1085 | index_t first_free_; | ||
| 1086 | bool weighted_; | ||
| 1087 | vector<double> heights_; // only used in weighted mode | ||
| 1088 | |||
| 1089 | /** | ||
| 1090 | * Performs additional checks (costly !) | ||
| 1091 | */ | ||
| 1092 | bool debug_mode_; | ||
| 1093 | |||
| 1094 | /** | ||
| 1095 | * Displays the result of the additional checks. | ||
| 1096 | */ | ||
| 1097 | bool verbose_debug_mode_; | ||
| 1098 | |||
| 1099 | /** | ||
| 1100 | * Displays the timing of the core algorithm. | ||
| 1101 | */ | ||
| 1102 | bool benchmark_mode_; | ||
| 1103 | |||
| 1104 | /** | ||
| 1105 | * \brief Gives the indexing of tetrahedron facet | ||
| 1106 | * vertices. | ||
| 1107 | * \details tet_facet_vertex[lf][lv] gives the | ||
| 1108 | * local vertex index (in 0,1,2,3) from a | ||
| 1109 | * local facet index lf (in 0,1,2,3) and a | ||
| 1110 | * local vertex index within the facet (in 0,1,2). | ||
| 1111 | */ | ||
| 1112 | static char tet_facet_vertex_[4][3]; | ||
| 1113 | |||
| 1114 | /** | ||
| 1115 | * \brief Gives a local facet index by | ||
| 1116 | * halfedge extremities local indices. | ||
| 1117 | */ | ||
| 1118 | static char halfedge_facet_[4][4]; | ||
| 1119 | |||
| 1120 | /** | ||
| 1121 | * \brief Used by the (de-recursified) | ||
| 1122 | * find_conflict_zone_iterative() function. | ||
| 1123 | */ | ||
| 1124 | std::stack<index_t> S_; | ||
| 1125 | |||
| 1126 | /** | ||
| 1127 | * \brief Used to represent the stack in the | ||
| 1128 | * (de-recursified) stellate_conflict_zone_iterative() | ||
| 1129 | * function. | ||
| 1130 | */ | ||
| 1131 | class StellateConflictStack { | ||
| 1132 | public: | ||
| 1133 | |||
| 1134 | /** | ||
| 1135 | * \brief Pushes a new frame onto the stack. | ||
| 1136 | * \param[in] t1 index of a tetrahedron on the border of | ||
| 1137 | * the conflict zone | ||
| 1138 | * \param[in] t1fbord index of the facet of \p t1 that is | ||
| 1139 | * on the border of the conflict zone | ||
| 1140 | * \param[in] t1fprev index of the facet of \p t1 that we | ||
| 1141 | * come from, or NO_INDEX if \p t1 is the first tetrahedron | ||
| 1142 | */ | ||
| 1143 | 500 | void push(index_t t1, index_t t1fbord, index_t t1fprev) { | |
| 1144 | 500 | store_.resize(store_.size()+1); | |
| 1145 | 500 | top().t1 = t1; | |
| 1146 | 500 | top().t1fbord = Numeric::uint8(t1fbord); | |
| 1147 | 500 | top().t1fprev = Numeric::uint8(t1fprev); | |
| 1148 | 500 | } | |
| 1149 | |||
| 1150 | /** | ||
| 1151 | * \brief Saves local variables into the current stack frame. | ||
| 1152 | * \param[in] new_t the index of the newly created tetrahedron | ||
| 1153 | * \param[in] t1ft2 the facet of t1 that is adjacent to t2 | ||
| 1154 | * \param[in] t2ft1 the facet of t2 that is adjacent to t1 | ||
| 1155 | */ | ||
| 1156 | 497 | void save_locals(index_t new_t, index_t t1ft2, index_t t2ft1) { | |
| 1157 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 497 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
497 | geo_debug_assert(!empty()); |
| 1158 | 497 | top().new_t = new_t; | |
| 1159 | 497 | top().t1ft2 = Numeric::uint8(t1ft2); | |
| 1160 | 497 | top().t2ft1 = Numeric::uint8(t2ft1); | |
| 1161 | 497 | } | |
| 1162 | |||
| 1163 | /** | ||
| 1164 | * \brief Gets the parameters from the current stack frame. | ||
| 1165 | * \param[out] t1 index of a tetrahedron on the border of | ||
| 1166 | * the conflict zone | ||
| 1167 | * \param[out] t1fbord index of the facet of \p t1 that is | ||
| 1168 | * on the border of the conflict zone | ||
| 1169 | * \param[out] t1fprev index of the facet of \p t1 that we | ||
| 1170 | * come from, or NO_INDEX if \p t1 is the first tetrahedron | ||
| 1171 | */ | ||
| 1172 | 997 | void get_parameters( | |
| 1173 | index_t& t1, index_t& t1fbord, index_t& t1fprev | ||
| 1174 | ) const { | ||
| 1175 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 997 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
997 | geo_debug_assert(!empty()); |
| 1176 | 997 | t1 = top().t1; | |
| 1177 | 997 | t1fbord = index_t(top().t1fbord); | |
| 1178 | 997 | t1fprev = index_t(top().t1fprev); | |
| 1179 | 997 | } | |
| 1180 | |||
| 1181 | |||
| 1182 | /** | ||
| 1183 | * \brief Gets the local variables from the current stack frame. | ||
| 1184 | * \param[out] new_t the index of the newly created tetrahedron | ||
| 1185 | * \param[out] t1ft2 the facet of t1 that is adjacent to t2 | ||
| 1186 | * \param[out] t2ft1 the facet of t2 that is adjacent to t1 | ||
| 1187 | */ | ||
| 1188 | 497 | void get_locals( | |
| 1189 | index_t& new_t, index_t& t1ft2, index_t& t2ft1 | ||
| 1190 | ) const { | ||
| 1191 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 497 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
497 | geo_debug_assert(!empty()); |
| 1192 | 497 | new_t = top().new_t; | |
| 1193 | 497 | t1ft2 = index_t(top().t1ft2); | |
| 1194 | 497 | t2ft1 = index_t(top().t2ft1); | |
| 1195 | 497 | } | |
| 1196 | |||
| 1197 | /** | ||
| 1198 | * \brief Pops a stack frame. | ||
| 1199 | */ | ||
| 1200 | 500 | void pop() { | |
| 1201 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 500 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
500 | geo_debug_assert(!empty()); |
| 1202 | 500 | store_.pop_back(); | |
| 1203 | 500 | } | |
| 1204 | |||
| 1205 | /** | ||
| 1206 | * \brief Tests whether the stack is empty. | ||
| 1207 | * \retval true if the stack is empty | ||
| 1208 | * \retval false otherwise | ||
| 1209 | */ | ||
| 1210 | 10464 | bool empty() const { | |
| 1211 | 10464 | return store_.empty(); | |
| 1212 | } | ||
| 1213 | |||
| 1214 | private: | ||
| 1215 | |||
| 1216 | /** | ||
| 1217 | * \brief The parameters and local | ||
| 1218 | * variables stored in a stack frame. | ||
| 1219 | */ | ||
| 1220 | struct Frame { | ||
| 1221 | // Parameters | ||
| 1222 | index_t t1; | ||
| 1223 | index_t new_t; | ||
| 1224 | Numeric::uint8 t1fbord ; | ||
| 1225 | |||
| 1226 | // Local variables | ||
| 1227 | Numeric::uint8 t1fprev ; | ||
| 1228 | Numeric::uint8 t1ft2 ; | ||
| 1229 | Numeric::uint8 t2ft1 ; | ||
| 1230 | }; | ||
| 1231 | |||
| 1232 | /** | ||
| 1233 | * \brief Gets the top of the stack. | ||
| 1234 | * \return a modifiable reference to the Frame on | ||
| 1235 | * the top of the stack | ||
| 1236 | * \pre !empty() | ||
| 1237 | */ | ||
| 1238 | 2991 | Frame& top() { | |
| 1239 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 2991 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
2991 | geo_debug_assert(!empty()); |
| 1240 | 2991 | return *store_.rbegin(); | |
| 1241 | } | ||
| 1242 | |||
| 1243 | /** | ||
| 1244 | * \brief Gets the top of the stack. | ||
| 1245 | * \return a const reference to the Frame on | ||
| 1246 | * the top of the stack | ||
| 1247 | * \pre !empty() | ||
| 1248 | */ | ||
| 1249 | 4482 | const Frame& top() const { | |
| 1250 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 4482 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
4482 | geo_debug_assert(!empty()); |
| 1251 | 4482 | return *store_.rbegin(); | |
| 1252 | } | ||
| 1253 | |||
| 1254 | std::vector<Frame> store_; | ||
| 1255 | }; | ||
| 1256 | |||
| 1257 | /** | ||
| 1258 | * \brief Used by the (de-recursified) | ||
| 1259 | * stellate_conflict_zone_iterative() function. | ||
| 1260 | */ | ||
| 1261 | StellateConflictStack S2_; | ||
| 1262 | |||
| 1263 | Cavity cavity_; | ||
| 1264 | }; | ||
| 1265 | |||
| 1266 | /************************************************************************/ | ||
| 1267 | |||
| 1268 | /** | ||
| 1269 | * \brief Regular Delaunay triangulation of weighted points | ||
| 1270 | * \details | ||
| 1271 | * - the input points are 4d points, were the fourth coordinate | ||
| 1272 | * of point \f$ i \f$ is \f$ \sqrt{W - w_i} \f$ where \f$ W \f$ is | ||
| 1273 | * the maximum of the weights of all the points and \d$ w_i \$ is | ||
| 1274 | * the weight associated with vertex \f$ i \f$. | ||
| 1275 | * - the constructed combinatorics is a tetrahedralized volume (3d and | ||
| 1276 | * not 4d although dimension() returns 4). This tetrahedralized volume | ||
| 1277 | * corresponds to the regular triangulation of the weighted points. | ||
| 1278 | */ | ||
| 1279 | class GEOGRAM_API RegularWeightedDelaunay3d : public Delaunay3d { | ||
| 1280 | public: | ||
| 1281 | /** | ||
| 1282 | * \brief Constructs a new Regular Delaunay3d triangulation. | ||
| 1283 | * \details RegularWeightedDelaunay3d triangulations are only | ||
| 1284 | * supported for dimension 3. If a different dimension is specified in | ||
| 1285 | * the constructor, a InvalidDimension exception is thrown. | ||
| 1286 | * \param[in] dimension dimension of the triangulation | ||
| 1287 | * \throw InvalidDimension This exception is thrown if dimension is | ||
| 1288 | * different than 3. | ||
| 1289 | */ | ||
| 1290 | RegularWeightedDelaunay3d(coord_index_t dimension = 4); | ||
| 1291 | |||
| 1292 | protected: | ||
| 1293 | /** | ||
| 1294 | * \brief RegularWeightedDelaunay3d destructor | ||
| 1295 | */ | ||
| 1296 | ~RegularWeightedDelaunay3d() override; | ||
| 1297 | }; | ||
| 1298 | } | ||
| 1299 | |||
| 1300 | #endif | ||
| 1301 |