| 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_2D | ||
| 41 | #define GEOGRAM_DELAUNAY_DELAUNAY_2D | ||
| 42 | |||
| 43 | #include <geogram/basic/common.h> | ||
| 44 | #include <geogram/delaunay/delaunay.h> | ||
| 45 | #include <geogram/numerics/predicates.h> | ||
| 46 | #include <geogram/basic/geometry.h> | ||
| 47 | |||
| 48 | #include <stack> | ||
| 49 | |||
| 50 | /** | ||
| 51 | * \file geogram/delaunay/delaunay_2d.h | ||
| 52 | * \brief Implementation of Delaunay in 2d. | ||
| 53 | */ | ||
| 54 | |||
| 55 | namespace GEO { | ||
| 56 | |||
| 57 | /** | ||
| 58 | * \brief Implementation of Delaunay in 2d. | ||
| 59 | * \details This package uses concepts inspired by | ||
| 60 | * two triangulation softwares, CGAL and tetgen, | ||
| 61 | * described in the following references. This package follows the | ||
| 62 | * idea used in CGAL of traversing the cavity from inside, since | ||
| 63 | * it traverses less triangles than when traversing from outside. | ||
| 64 | * - Jean-Daniel Boissonnat, Olivier Devillers, Monique Teillaud, | ||
| 65 | * and Mariette Yvinec. Triangulations in CGAL. | ||
| 66 | * In Proc. 16th Annu. ACM Sympos. Comput. Geom., pages 11-18, 2000. | ||
| 67 | * - Hang Si, Constrained Delaunay trianglesl mesh generation and | ||
| 68 | * refinement. Finite elements in Analysis and Design, | ||
| 69 | * 46 (1-2):33--46, 2010. | ||
| 70 | * | ||
| 71 | * Note that the algorithm here does not support vertex deletion nor | ||
| 72 | * degenerate input with all coplanar or all colinear points (use CGAL | ||
| 73 | * instead if you have these requirements). | ||
| 74 | * | ||
| 75 | * The core algorithm used in this code, CGAL and tetgen was | ||
| 76 | * independently and simultaneously discovered by Bowyer and Watson: | ||
| 77 | * - Adrian Bowyer, "Computing Dirichlet tessellations", | ||
| 78 | * Comput. J., vol. 24, no 2, 1981, p. 162-166 | ||
| 79 | * - David F. Watson, "Computing the n-dimensional Delaunay tessellation | ||
| 80 | * with application to Voronoi polytopes", Comput. J., vol. 24, | ||
| 81 | * no 2, 1981, p. 167-172 | ||
| 82 | * | ||
| 83 | * The spatial reordering method, that dramatically increases the | ||
| 84 | * performances, also used in this code, CGAL and tetgen was introduced | ||
| 85 | * in the following references. The second one is a smart implementation | ||
| 86 | * based on the std::nth_element() function of the STL, that inspired | ||
| 87 | * the compute_BRIO_ordering() function of this package. | ||
| 88 | * - Nina Amenta, Sunghee Choi and Gunter Rote, "Incremental constructions | ||
| 89 | * con brio", ACM Symposium on Computational Geometry 2003. | ||
| 90 | * - Christophe Delage and Olivier Devillers. Spatial Sorting. | ||
| 91 | * In CGAL User and Reference Manual. CGAL Editorial Board, | ||
| 92 | * 3.9 edition, 2011 | ||
| 93 | * | ||
| 94 | * The locate() function is based on the following two references. | ||
| 95 | * The first one randomizes the choice of the next triangle. | ||
| 96 | * The second one uses an inexact locate() function to initialize | ||
| 97 | * the exact one (it is called "structural filtering"). The first | ||
| 98 | * idea is used in both CGAL and tetgen, and the second one is used | ||
| 99 | * in CGAL. | ||
| 100 | * - Walking in a triangulation, O Devillers, S Pion, M Teillaud | ||
| 101 | * 17th Annual Symposium on Computational geometry, 106-114 | ||
| 102 | * - Stefan Funke , Kurt Mehlhorn and Stefan Naher, "Structural filtering, | ||
| 103 | * a paradigm for efficient and exact geometric programs", | ||
| 104 | * Comput. Geom., 1999 | ||
| 105 | */ | ||
| 106 | class GEOGRAM_API Delaunay2d : public Delaunay { | ||
| 107 | public: | ||
| 108 | /** | ||
| 109 | * \brief Constructs a new Delaunay2d. | ||
| 110 | * \param[in] dimension dimension of the triangulation (2 or 3). | ||
| 111 | * If dimension = 3, this creates a regular triangulation | ||
| 112 | * (dual of a power diagram). In this case: | ||
| 113 | * - the input points are 3d points, were the third coordinate | ||
| 114 | * of point \f$ i \f$ is \f$ \sqrt{W - w_i} \f$ where \f$ W \f$ is | ||
| 115 | * the maximum of the weights of all the points and \d$ w_i \$ is | ||
| 116 | * the weight associated with vertex \f$ i \f$. | ||
| 117 | * - the constructed combinatorics is a triangulated surface (2d and | ||
| 118 | * not 3d although dimension() returns 3). This triangulated surface | ||
| 119 | * corresponds to the regular triangulation of the weighted points. | ||
| 120 | */ | ||
| 121 | Delaunay2d(coord_index_t dimension = 2); | ||
| 122 | |||
| 123 | /** | ||
| 124 | * \copydoc Delaunay::set_vertices() | ||
| 125 | */ | ||
| 126 | void set_vertices(index_t nb_vertices, const double* vertices) override; | ||
| 127 | |||
| 128 | /** | ||
| 129 | * \copydoc Delaunay::nearest_vertex() | ||
| 130 | */ | ||
| 131 | index_t nearest_vertex(const double* p) const override; | ||
| 132 | |||
| 133 | /** | ||
| 134 | * \brief Tests whether the Laguerre diagram has empty cells. | ||
| 135 | * \details If the Laguerre diagram has empty cells and | ||
| 136 | * abort_if_empty_cell is set, then computation is stopped, | ||
| 137 | * and all the queries on the Laguerre diagram will not work | ||
| 138 | * (including the non-empty cells). | ||
| 139 | * \retval true if the Laguerre diagram has empty cells. | ||
| 140 | * \retval false otherwise. | ||
| 141 | */ | ||
| 142 | bool has_empty_cells() const { | ||
| 143 | return has_empty_cells_; | ||
| 144 | } | ||
| 145 | |||
| 146 | /** | ||
| 147 | * \brief Specifies behavior if an empty cell is detected. | ||
| 148 | * \param[in] x if set, then computation is aborted as soon | ||
| 149 | * as an empty cell is detected. | ||
| 150 | * \details only happens in RegularTriangulation/Laguerre diagram. | ||
| 151 | */ | ||
| 152 | void abort_if_empty_cell(bool x) { | ||
| 153 | abort_if_empty_cell_ = x; | ||
| 154 | } | ||
| 155 | |||
| 156 | protected: | ||
| 157 | |||
| 158 | /** | ||
| 159 | * \brief Symbolic constant for uninitialized hint. | ||
| 160 | * \details Locate functions can be accelerated by | ||
| 161 | * specifying a hint. This constant indicates that | ||
| 162 | * no hint is given. | ||
| 163 | */ | ||
| 164 | static constexpr index_t NO_TRIANGLE = NO_INDEX; | ||
| 165 | |||
| 166 | /** | ||
| 167 | * \brief Finds in the pointset a set of three non-colinear | ||
| 168 | * points. | ||
| 169 | * \details This function is used to initiate the incremental | ||
| 170 | * Delaunay construction. | ||
| 171 | * \param[out] iv0 index of the first vertex | ||
| 172 | * \param[out] iv1 index of the second vertex | ||
| 173 | * \param[out] iv2 index of the third vertex | ||
| 174 | * \retval true if a set of three non-colinear points was found | ||
| 175 | * \retval false if all the points are colinear | ||
| 176 | */ | ||
| 177 | bool create_first_triangle(index_t& iv0, index_t& iv1, index_t& iv2); | ||
| 178 | |||
| 179 | /** | ||
| 180 | * \brief Finds the triangle that contains a point. | ||
| 181 | * \details If the point is on an edge or vertex, | ||
| 182 | * the function returns one of the triangles incident | ||
| 183 | * to that edge or vertex. | ||
| 184 | * \param[in] p a pointer to the coordinates of the point | ||
| 185 | * \param[in] thread_safe if true, a global spinlock is | ||
| 186 | * used to protect the calls to random(), this is necessary | ||
| 187 | * if multiple threads use locate() simultaneously | ||
| 188 | * \param[out] orient a pointer to an array of three Sign%s | ||
| 189 | * or nullptr. If non-nullptr, returns the orientation with respect | ||
| 190 | * to the three edges of the triangle that contains \p p. | ||
| 191 | * \return the index of a triangle that contains \p p. | ||
| 192 | * If the point is outside the convex hull of | ||
| 193 | * the inserted so-far points, then the returned triangle | ||
| 194 | * is a virtual one (first vertex is the "vertex at infinity" | ||
| 195 | * of index -1) or NO_TRIANGLE if the virtual triangles | ||
| 196 | * were previously removed. | ||
| 197 | */ | ||
| 198 | index_t locate( | ||
| 199 | const double* p, index_t hint = NO_TRIANGLE, | ||
| 200 | bool thread_safe = false, | ||
| 201 | Sign* orient = nullptr | ||
| 202 | ) const; | ||
| 203 | |||
| 204 | /** | ||
| 205 | * \brief Finds the triangle that (approximately) | ||
| 206 | * contains a point using inexact predicates. | ||
| 207 | * \details The result of this function can be used as a hint | ||
| 208 | * for locate(). It accelerates locate as compared to calling | ||
| 209 | * it directly. This technique is referred to as "structural | ||
| 210 | * filtering". | ||
| 211 | * \param[in] p a pointer to the coordinates of the point | ||
| 212 | * \param[in] max_iter maximum number of traversed tets | ||
| 213 | * \return the index of a triangle that (approximately) | ||
| 214 | * contains \p p. | ||
| 215 | * If the point is outside the convex hull of | ||
| 216 | * the inserted so-far points, then the returned triangle | ||
| 217 | * is a virtual one (first vertex is the "vertex at infinity" | ||
| 218 | * of index -1) or NO_TRIANGLE if the virtual triangles | ||
| 219 | * were previously removed. | ||
| 220 | */ | ||
| 221 | index_t locate_inexact( | ||
| 222 | const double* p, index_t hint, index_t max_iter | ||
| 223 | ) const; | ||
| 224 | |||
| 225 | /** | ||
| 226 | * \brief Inserts a point in the triangulation. | ||
| 227 | * \param[in] v the index of the point to be inserted | ||
| 228 | * \param[in] hint the index of a triangle as near as | ||
| 229 | * possible to \p v, or -1 if unspecified | ||
| 230 | * \return the index of one of the triangles incident to | ||
| 231 | * point \p v | ||
| 232 | */ | ||
| 233 | index_t insert(index_t v, index_t hint = NO_TRIANGLE); | ||
| 234 | |||
| 235 | /** | ||
| 236 | * \brief Determines the list of triangles in conflict | ||
| 237 | * with a given point. | ||
| 238 | * \param[in] v the index of the point to be inserted | ||
| 239 | * \param[in] t the index of a triangle that contains | ||
| 240 | * \p p, as returned by locate() | ||
| 241 | * \param[in] orient an array of three signs indicating | ||
| 242 | * the orientation of \p p with respect to the three | ||
| 243 | * edges of \p t, as returned by locate() | ||
| 244 | * \param[out] t_bndry a triangle adjacent to the | ||
| 245 | * boundary of the conflict zone | ||
| 246 | * \param[out] e_bndry the edge along which t_bndry is | ||
| 247 | * adjacent to the boundary of the conflict zone | ||
| 248 | * \param[out] first the index of the first triangle in conflict | ||
| 249 | * \param[out] last the index of the last triangle in conflict | ||
| 250 | * The other triangles are linked, and can be traversed | ||
| 251 | * from \p first by using triangle_next() until \p last or END_OF_LIST | ||
| 252 | * is reached. | ||
| 253 | * The conflict zone can be empty under two circumstances: | ||
| 254 | * - the vertex \p v already exists in the triangulation | ||
| 255 | * - the triangulation is weighted and \p v is not visible | ||
| 256 | * in either cases, both \p first and \p last contain END_OF_LIST | ||
| 257 | */ | ||
| 258 | void find_conflict_zone( | ||
| 259 | index_t v, | ||
| 260 | index_t t, const Sign* orient, | ||
| 261 | index_t& t_bndry, index_t& e_bndry, | ||
| 262 | index_t& first, index_t& last | ||
| 263 | ); | ||
| 264 | |||
| 265 | /** | ||
| 266 | * \brief This function is used to implement find_conflict_zone. | ||
| 267 | * \details This function detects the neighbors of \p t that are | ||
| 268 | * in the conflict zone and calls itself recursively on them. | ||
| 269 | * \param[in] p the point to be inserted | ||
| 270 | * \param[in] t index of a triangle in the fonflict zone | ||
| 271 | * \param[out] t_bndry a triangle adjacent to the | ||
| 272 | * boundary of the conflict zone | ||
| 273 | * \param[out] e_bndry the edge along which t_bndry is | ||
| 274 | * adjacent to the boundary of the conflict zone | ||
| 275 | * \param[out] first the index of the first triangle in conflict | ||
| 276 | * \param[out] last the index of the last triangle in conflict | ||
| 277 | * \pre The triangle \p t was already marked as | ||
| 278 | * conflict (triangle_is_in_list(t)) | ||
| 279 | */ | ||
| 280 | void find_conflict_zone_iterative( | ||
| 281 | const double* p, index_t t, | ||
| 282 | index_t& t_bndry, index_t& e_bndry, | ||
| 283 | index_t& first, index_t& last | ||
| 284 | ); | ||
| 285 | |||
| 286 | /** | ||
| 287 | * \brief Creates a star of triangles filling the conflict | ||
| 288 | * zone. | ||
| 289 | * \details For each triangle edge on the border of the | ||
| 290 | * conflict zone, a new triangle is created, resting on | ||
| 291 | * the edge and incident to vertex \p v. The function is | ||
| 292 | * called recursively until the entire conflict zone is filled. | ||
| 293 | * \param[in] v the index of the point to be inserted | ||
| 294 | * \param[in] t_bndry index of a triangle on the border | ||
| 295 | * of the conflict zone. | ||
| 296 | * \param[in] e_bndry index of the facet along which \p t_bndry | ||
| 297 | * is incident to the border of the conflict zone | ||
| 298 | * \return the index of one the newly created triangles | ||
| 299 | */ | ||
| 300 | index_t stellate_conflict_zone( | ||
| 301 | index_t v, | ||
| 302 | index_t t_bndry, index_t e_bndry | ||
| 303 | ); | ||
| 304 | |||
| 305 | /*** Combinatorics - new and delete ******************************/ | ||
| 306 | |||
| 307 | /** | ||
| 308 | * \brief Maximum valid index for a triangle. | ||
| 309 | * \details This includes not only real triangles, | ||
| 310 | * but also the virtual ones on the border, the conflict | ||
| 311 | * list and the free list. | ||
| 312 | * \return the maximum valid index for a triangle. | ||
| 313 | */ | ||
| 314 | 35592 | index_t max_t() const { | |
| 315 | 35592 | return cell_to_v_store_.size() / 3; | |
| 316 | } | ||
| 317 | |||
| 318 | |||
| 319 | /** | ||
| 320 | * \brief Default symbolic value of the cell_next_ field | ||
| 321 | * that indicates that a triangle is not | ||
| 322 | * in a linked list. | ||
| 323 | * \details This is the default value. Note that it suffices | ||
| 324 | * that NOT_IN_LIST_BIT is set for a triangle | ||
| 325 | * to be not in any list. | ||
| 326 | * A triangle can be: | ||
| 327 | * - in a list (cell_next_[t] & NOT_IN_LIST_BIT == 0) | ||
| 328 | * - not in a list and not marked | ||
| 329 | * (cell_next_[t] & NOT_IN_LIST_BIT != 0) && | ||
| 330 | * (cell_next_[t] != cur_stamp_) | ||
| 331 | * - not in a list and marked | ||
| 332 | * (cell_next_[t] == cur_stamp_) | ||
| 333 | */ | ||
| 334 | static constexpr index_t NOT_IN_LIST = ~index_t(0); | ||
| 335 | |||
| 336 | /** | ||
| 337 | * \brief If cell_next_[t] & NOT_IN_LIST_BIT != 0, | ||
| 338 | * then t is not in a linked list. | ||
| 339 | * \details The other bits of cell_next_[t] are used | ||
| 340 | * to store the stamp (i.e. index of the current point | ||
| 341 | * being inserted). The stamp is used for marking triangles | ||
| 342 | * that were detected as non-conflict when inserting a point. | ||
| 343 | * A triangle can be: | ||
| 344 | * - in a list (cell_next_[t] & NOT_IN_LIST_BIT == 0) | ||
| 345 | * - not in a list and not marked | ||
| 346 | * (cell_next_[t] & NOT_IN_LIST_BIT != 0) && | ||
| 347 | * (cell_next_[t] != cur_stamp_) | ||
| 348 | * - not in a list and marked | ||
| 349 | * (cell_next_[t] == cur_stamp_) | ||
| 350 | */ | ||
| 351 | static constexpr index_t NOT_IN_LIST_BIT = | ||
| 352 | index_t(1) << (sizeof(index_t)*8-1) ; | ||
| 353 | |||
| 354 | /** | ||
| 355 | * \brief Symbolic value of the cell_next_ field | ||
| 356 | * that indicates the end of list in a linked | ||
| 357 | * list of triangles. | ||
| 358 | */ | ||
| 359 | static constexpr index_t END_OF_LIST = ~NOT_IN_LIST_BIT; | ||
| 360 | |||
| 361 | /** | ||
| 362 | * \brief Tests whether a triangle belongs to a linked | ||
| 363 | * list. | ||
| 364 | * \details Triangles can be linked, it is used to manage | ||
| 365 | * both the free list that recycles deleted triangles, | ||
| 366 | * the conflict region and the list of newly created | ||
| 367 | * triangles. In addition, a triangle that is not | ||
| 368 | * in a list can be marked. The same space is used for | ||
| 369 | * marking and chaining triangles in lists. | ||
| 370 | * A triangle can be in the following states: | ||
| 371 | * - in list | ||
| 372 | * - not in list and marked | ||
| 373 | * - not in list and not marked | ||
| 374 | * \param[in] t the index of the triangle | ||
| 375 | * \retval true if triangle \p t belongs to a linked list | ||
| 376 | * \retval false otherwise | ||
| 377 | */ | ||
| 378 | 9169 | bool triangle_is_in_list(index_t t) const { | |
| 379 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 9169 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
9169 | geo_debug_assert(t < max_t()); |
| 380 | 9169 | return (cell_next_[t] & NOT_IN_LIST_BIT) == 0; | |
| 381 | } | ||
| 382 | |||
| 383 | /** | ||
| 384 | * \brief Gets the index of a successor of a triangle. | ||
| 385 | * \details Triangles can be linked, it is used to manage | ||
| 386 | * both the free list that recycles deleted triangles. | ||
| 387 | * \param[in] t the index of the triangle | ||
| 388 | * \retval END_OF_LIST if the end of the list is reached | ||
| 389 | * \retval the index of the successor of | ||
| 390 | * triangle \t otherwise | ||
| 391 | * \pre triangle_is_in_list(t) | ||
| 392 | */ | ||
| 393 | 452 | index_t triangle_next(index_t t) const { | |
| 394 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 452 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
452 | geo_debug_assert(t < max_t()); |
| 395 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 452 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
452 | geo_debug_assert(triangle_is_in_list(t)); |
| 396 | 452 | return cell_next_[t]; | |
| 397 | } | ||
| 398 | |||
| 399 | /** | ||
| 400 | * \brief Adds a triangle to a linked list. | ||
| 401 | * \details Triangles can be linked, it is used to manage | ||
| 402 | * the free list that recycles deleted triangles. | ||
| 403 | * \param[in] t the index of the triangle | ||
| 404 | * \param[in,out] first first item of the list or END_OF_LIST if | ||
| 405 | * the list is empty | ||
| 406 | * \param[in,out] last last item of the list or END_OF_LIST if | ||
| 407 | * the list is empty | ||
| 408 | */ | ||
| 409 | 468 | void add_triangle_to_list(index_t t, index_t& first, index_t& last) { | |
| 410 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 468 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
468 | geo_debug_assert(t < max_t()); |
| 411 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 468 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
468 | geo_debug_assert(!triangle_is_in_list(t)); |
| 412 |
2/2✓ Branch 0 taken 127 times.
✓ Branch 1 taken 341 times.
|
468 | if(last == END_OF_LIST) { |
| 413 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 127 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
127 | geo_debug_assert(first == END_OF_LIST); |
| 414 | 127 | first = last = t; | |
| 415 | 127 | cell_next_[t] = END_OF_LIST; | |
| 416 | } else { | ||
| 417 | 341 | cell_next_[t] = first; | |
| 418 | 341 | first = t; | |
| 419 | } | ||
| 420 | 468 | } | |
| 421 | |||
| 422 | /** | ||
| 423 | * \brief Removes a triangle from the linked list it | ||
| 424 | * belongs to. | ||
| 425 | * \details Triangles can be linked, it is used to manage | ||
| 426 | * both the free list that recycles deleted triangles and | ||
| 427 | * the list of triangles in conflict with the inserted | ||
| 428 | * point. | ||
| 429 | * \param[in] t the index of the triangle | ||
| 430 | */ | ||
| 431 | 452 | void remove_triangle_from_list(index_t t) { | |
| 432 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 452 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
452 | geo_debug_assert(t < max_t()); |
| 433 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 452 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
452 | geo_debug_assert(triangle_is_in_list(t)); |
| 434 | 452 | cell_next_[t] = NOT_IN_LIST; | |
| 435 | 452 | } | |
| 436 | |||
| 437 | /** | ||
| 438 | * \brief Symbolic value for a vertex of a | ||
| 439 | * triangle that indicates a virtual triangle. | ||
| 440 | * \details The three other vertices then correspond to a | ||
| 441 | * facet on the convex hull of the points. | ||
| 442 | */ | ||
| 443 | static constexpr index_t VERTEX_AT_INFINITY = NO_INDEX; | ||
| 444 | |||
| 445 | /** | ||
| 446 | * \brief Tests whether a given triangle | ||
| 447 | * is a finite one. | ||
| 448 | * \details Infinite triangles are the ones | ||
| 449 | * that are incident to the infinite vertex | ||
| 450 | * (index -1) | ||
| 451 | * \param[in] t the index of the triangle | ||
| 452 | * \retval true if \p t is finite | ||
| 453 | * \retval false otherwise | ||
| 454 | */ | ||
| 455 | 262 | bool triangle_is_finite(index_t t) const { | |
| 456 | return | ||
| 457 | 262 | (cell_to_v_store_[3 * t] != NO_INDEX) && | |
| 458 |
4/4✓ Branch 0 taken 259 times.
✓ Branch 1 taken 3 times.
✓ Branch 3 taken 201 times.
✓ Branch 4 taken 58 times.
|
463 | (cell_to_v_store_[3 * t + 1] != NO_INDEX) && |
| 459 |
2/2✓ Branch 1 taken 153 times.
✓ Branch 2 taken 48 times.
|
463 | (cell_to_v_store_[3 * t + 2] != NO_INDEX) ; |
| 460 | } | ||
| 461 | |||
| 462 | /** | ||
| 463 | * \brief Tests whether a triangle is | ||
| 464 | * a real one. | ||
| 465 | * \details Real triangles are incident to | ||
| 466 | * three user-specified vertices (there are also | ||
| 467 | * virtual triangles that are incident to the | ||
| 468 | * vertex at infinity, with index -1) | ||
| 469 | * \param[in] t index of the triangle | ||
| 470 | * \retval true if triangle \p t is a real one | ||
| 471 | * \retval false otherwise | ||
| 472 | */ | ||
| 473 | 278 | bool triangle_is_real(index_t t) const { | |
| 474 |
4/4✓ Branch 1 taken 262 times.
✓ Branch 2 taken 16 times.
✓ Branch 4 taken 153 times.
✓ Branch 5 taken 109 times.
|
278 | return !triangle_is_free(t) && triangle_is_finite(t); |
| 475 | } | ||
| 476 | |||
| 477 | /** | ||
| 478 | * \brief Tests whether a triangle is | ||
| 479 | * a virtual one. | ||
| 480 | * \details Virtual triangles are triangles | ||
| 481 | * incident to the vertex at infinity. | ||
| 482 | * \param[in] t index of the triangle | ||
| 483 | * \retval true if triangle \p t is virtual | ||
| 484 | * \retval false otherwise | ||
| 485 | */ | ||
| 486 | 1675 | bool triangle_is_virtual(index_t t) const { | |
| 487 | return | ||
| 488 |
1/2✓ Branch 1 taken 1675 times.
✗ Branch 2 not taken.
|
3350 | !triangle_is_free(t) && ( |
| 489 |
2/2✓ Branch 1 taken 1664 times.
✓ Branch 2 taken 11 times.
|
1675 | cell_to_v_store_[3 * t] == VERTEX_AT_INFINITY || |
| 490 |
2/2✓ Branch 1 taken 1554 times.
✓ Branch 2 taken 110 times.
|
1664 | cell_to_v_store_[3 * t + 1] == VERTEX_AT_INFINITY || |
| 491 |
2/2✓ Branch 1 taken 54 times.
✓ Branch 2 taken 1500 times.
|
1554 | cell_to_v_store_[3 * t + 2] == VERTEX_AT_INFINITY |
| 492 | 1675 | ); | |
| 493 | } | ||
| 494 | |||
| 495 | /** | ||
| 496 | * \brief Tests whether a triangle is | ||
| 497 | * in the free list. | ||
| 498 | * \details Deleted triangles are recycled | ||
| 499 | * in a free list. | ||
| 500 | * \param[in] t index of the triangle | ||
| 501 | * \retval true if triangle \p t is in | ||
| 502 | * the free list | ||
| 503 | * \retval false otherwise | ||
| 504 | */ | ||
| 505 | 3264 | bool triangle_is_free(index_t t) const { | |
| 506 | 3264 | return triangle_is_in_list(t); | |
| 507 | } | ||
| 508 | |||
| 509 | /** | ||
| 510 | * \brief Creates a new triangle. | ||
| 511 | * \details Uses either a triangle recycled | ||
| 512 | * from the free list, or creates a new one by | ||
| 513 | * expanding the two indices arrays. | ||
| 514 | * \return the index of the newly created triangle | ||
| 515 | */ | ||
| 516 | 730 | index_t new_triangle() { | |
| 517 | index_t result; | ||
| 518 |
2/2✓ Branch 0 taken 278 times.
✓ Branch 1 taken 452 times.
|
730 | if(first_free_ == END_OF_LIST) { |
| 519 | 556 | cell_to_v_store_.resize( | |
| 520 | 278 | cell_to_v_store_.size() + 3, NO_INDEX | |
| 521 | ); | ||
| 522 | 556 | cell_to_cell_store_.resize( | |
| 523 | 278 | cell_to_cell_store_.size() + 3, NO_INDEX | |
| 524 | ); | ||
| 525 | // index_t(NOT_IN_LIST) is necessary, else with | ||
| 526 | // NOT_IN_LIST alone the compiler tries to generate a | ||
| 527 | // reference to NOT_IN_LIST resulting in a link error. | ||
| 528 |
1/2✓ Branch 1 taken 278 times.
✗ Branch 2 not taken.
|
278 | cell_next_.push_back(index_t(NOT_IN_LIST)); |
| 529 | 278 | result = max_t() - 1; | |
| 530 | } else { | ||
| 531 | 452 | result = first_free_; | |
| 532 | 452 | first_free_ = triangle_next(first_free_); | |
| 533 | 452 | remove_triangle_from_list(result); | |
| 534 | } | ||
| 535 | |||
| 536 | 730 | cell_to_cell_store_[3 * result] = NO_INDEX; | |
| 537 | 730 | cell_to_cell_store_[3 * result + 1] = NO_INDEX; | |
| 538 | 730 | cell_to_cell_store_[3 * result + 2] = NO_INDEX; | |
| 539 | |||
| 540 | 730 | return result; | |
| 541 | } | ||
| 542 | |||
| 543 | /** | ||
| 544 | * \brief Creates a new triangle. | ||
| 545 | * \details Sets the vertices. Adjacent triangles index are | ||
| 546 | * left uninitialized. Uses either a triangle recycled | ||
| 547 | * from the free list, or creates a new one by | ||
| 548 | * expanding the two indices arrays. | ||
| 549 | * \param[in] v1 index of the first vertex | ||
| 550 | * \param[in] v2 index of the second vertex | ||
| 551 | * \param[in] v3 index of the third vertex | ||
| 552 | * \return the index of the newly created triangle | ||
| 553 | */ | ||
| 554 | 730 | index_t new_triangle(index_t v1, index_t v2, index_t v3) { | |
| 555 | 730 | index_t result = new_triangle(); | |
| 556 | 730 | cell_to_v_store_[3 * result] = v1; | |
| 557 | 730 | cell_to_v_store_[3 * result + 1] = v2; | |
| 558 | 730 | cell_to_v_store_[3 * result + 2] = v3; | |
| 559 | 730 | return result; | |
| 560 | } | ||
| 561 | |||
| 562 | /** | ||
| 563 | * \brief Generates a unique stamp for marking tets. | ||
| 564 | * \details Storage is shared for list-chaining and stamp-marking | ||
| 565 | * (both are mutually exclusive), therefore the stamp has | ||
| 566 | * the NOT_IN_LIST_BIT set. | ||
| 567 | * \param[in] stamp the unique stamp for marking tets | ||
| 568 | */ | ||
| 569 | 223 | void set_triangle_mark_stamp(index_t stamp) { | |
| 570 | 223 | cur_stamp_ = (stamp | NOT_IN_LIST_BIT); | |
| 571 | 223 | } | |
| 572 | |||
| 573 | /** | ||
| 574 | * \brief Tests whether a triangle is marked. | ||
| 575 | * \details A triangle is marked whenever it is | ||
| 576 | * detected as non-conflict. The index of the | ||
| 577 | * point being inserted is used as a time-stamp | ||
| 578 | * for marking triangles. The same space is used | ||
| 579 | * for marking and for chaining the conflict list. | ||
| 580 | * A triangle can be in the following states: | ||
| 581 | * - in list | ||
| 582 | * - not in list and marked | ||
| 583 | * - not in list and not marked | ||
| 584 | * \param[in] t index of the triangle | ||
| 585 | * \retval true if triangle \p t is marked | ||
| 586 | * \retval false otherwise | ||
| 587 | */ | ||
| 588 | 1139 | bool triangle_is_marked(index_t t) const { | |
| 589 | 1139 | return cell_next_[t] == cur_stamp_; | |
| 590 | } | ||
| 591 | |||
| 592 | /** | ||
| 593 | * \brief Marks a triangle. | ||
| 594 | * \details A triangle is marked whenever it is | ||
| 595 | * detected as non-conflict. The same space is used | ||
| 596 | * for marking and for chaining the conflict list. | ||
| 597 | * The index of the point being inserted is used as a | ||
| 598 | * time-stamp for marking triangles. | ||
| 599 | * A triangle can be in the following states: | ||
| 600 | * - in list | ||
| 601 | * - not in list and marked | ||
| 602 | * - not in list and not marked | ||
| 603 | * \param[in] t index of the triangle to be marked | ||
| 604 | */ | ||
| 605 | 707 | void mark_triangle(index_t t) { | |
| 606 | 707 | cell_next_[t] = cur_stamp_; | |
| 607 | 707 | } | |
| 608 | |||
| 609 | /**** Combinatorics ****************************************/ | ||
| 610 | |||
| 611 | /** | ||
| 612 | * \brief Returns the local index of a vertex by | ||
| 613 | * edge and by local vertex index in the edge. | ||
| 614 | * \details | ||
| 615 | * tri edge vertex is such that the triangle | ||
| 616 | * formed with: | ||
| 617 | * - vertex lv | ||
| 618 | * - triangle_edge_vertex(lv,0) | ||
| 619 | * - triangle_edge_vertex(lv,1) | ||
| 620 | * has the same orientation as the original triangle for | ||
| 621 | * any vertex lv. | ||
| 622 | * \param[in] e local facet index, in (0,1,2) | ||
| 623 | * \param[in] v local vertex index, in (0,1) | ||
| 624 | * \return the local triangle vertex index of | ||
| 625 | * vertex \p v in edge \p f | ||
| 626 | */ | ||
| 627 | 24 | static index_t triangle_edge_vertex(index_t e, index_t v) { | |
| 628 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 24 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
24 | geo_debug_assert(e < 3); |
| 629 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 24 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
24 | geo_debug_assert(v < 2); |
| 630 | 24 | return index_t(triangle_edge_vertex_[e][v]); | |
| 631 | } | ||
| 632 | |||
| 633 | /** | ||
| 634 | * \brief Gets the index of a vertex of a triangle | ||
| 635 | * \param[in] t index of the triangle | ||
| 636 | * \param[in] lv local vertex (0,1,2) index in \p t | ||
| 637 | * \return the global index of the \p lv%th vertex of triangle \p t | ||
| 638 | * or -1 if the vertex is at infinity | ||
| 639 | */ | ||
| 640 | 5043 | index_t triangle_vertex(index_t t, index_t lv) const { | |
| 641 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 5043 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
5043 | geo_debug_assert(t < max_t()); |
| 642 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 5043 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
5043 | geo_debug_assert(lv < 3); |
| 643 | 5043 | return cell_to_v_store_[3 * t + lv]; | |
| 644 | } | ||
| 645 | |||
| 646 | /** | ||
| 647 | * \brief Finds the index of the vertex in a triangle. | ||
| 648 | * \param[in] t the triangle | ||
| 649 | * \param[in] v the vertex | ||
| 650 | * \return iv such that triangle_vertex(t,v)==iv | ||
| 651 | * \pre \p t is incident to \p v | ||
| 652 | */ | ||
| 653 | 682 | index_t find_triangle_vertex(index_t t, index_t v) const { | |
| 654 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 682 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
682 | geo_debug_assert(t < max_t()); |
| 655 | // Find local index of v in triangle t vertices. | ||
| 656 | 682 | const index_t* T = &(cell_to_v_store_[3 * t]); | |
| 657 | 682 | return find_3(T,v); | |
| 658 | } | ||
| 659 | |||
| 660 | |||
| 661 | /** | ||
| 662 | * \brief Gets the index of a vertex of a triangle | ||
| 663 | * \param[in] t index of the triangle | ||
| 664 | * \param[in] lv local vertex (0,1,2) index in \p t | ||
| 665 | * \return the global index of the \p lv%th vertex of triangle \p t | ||
| 666 | * \pre Vertex \p lv of triangle \p t is not at infinity | ||
| 667 | */ | ||
| 668 | 3660 | index_t finite_triangle_vertex(index_t t, index_t lv) const { | |
| 669 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 3660 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
3660 | geo_debug_assert(t < max_t()); |
| 670 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 3660 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
3660 | geo_debug_assert(lv < 3); |
| 671 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 3660 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
3660 | geo_debug_assert(cell_to_v_store_[3 * t + lv] != NO_INDEX); |
| 672 | 3660 | return cell_to_v_store_[3 * t + lv]; | |
| 673 | } | ||
| 674 | |||
| 675 | /** | ||
| 676 | * \brief Sets a triangle-to-vertex adjacency. | ||
| 677 | * \param[in] t index of the triangle | ||
| 678 | * \param[in] lv local vertex index (0,1,2) in \p t | ||
| 679 | * \param[in] v global index of the vertex | ||
| 680 | */ | ||
| 681 | void set_triangle_vertex(index_t t, index_t lv, index_t v) { | ||
| 682 | geo_debug_assert(t < max_t()); | ||
| 683 | geo_debug_assert(lv < 3); | ||
| 684 | cell_to_v_store_[3 * t + lv] = v; | ||
| 685 | } | ||
| 686 | |||
| 687 | /** | ||
| 688 | * \brief Gets the index of a triangle adjacent to another one. | ||
| 689 | * \param[in] t index of the triangle | ||
| 690 | * \param[in] le local edge (0,1,2) index in \p t | ||
| 691 | * \return the triangle adjacent to \p t accros edge \p le | ||
| 692 | */ | ||
| 693 | 7838 | index_t triangle_adjacent(index_t t, index_t le) const { | |
| 694 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 7838 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
7838 | geo_debug_assert(t < max_t()); |
| 695 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 7838 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
7838 | geo_debug_assert(le < 3); |
| 696 | 7838 | index_t result = cell_to_cell_store_[3 * t + le]; | |
| 697 | 7838 | return result; | |
| 698 | } | ||
| 699 | |||
| 700 | /** | ||
| 701 | * \brief Sets a triangle-to-triangle adjacency. | ||
| 702 | * \param[in] t1 index of the first triangle | ||
| 703 | * \param[in] le1 local facet index (0,1,2) in t1 | ||
| 704 | * \param[in] t2 index of the triangle | ||
| 705 | * adjacent to \p t1 accros \p lf1 | ||
| 706 | */ | ||
| 707 | 2912 | void set_triangle_adjacent(index_t t1, index_t le1, index_t t2) { | |
| 708 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 2912 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
2912 | geo_debug_assert(t1 < max_t()); |
| 709 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 2912 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
2912 | geo_debug_assert(t2 < max_t()); |
| 710 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 2912 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
2912 | geo_debug_assert(le1 < 3); |
| 711 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 2912 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
2912 | geo_debug_assert(t1 != t2); |
| 712 | 2912 | cell_to_cell_store_[3 * t1 + le1] = t2; | |
| 713 | 2912 | } | |
| 714 | |||
| 715 | /** | ||
| 716 | * \brief Finds the index of the edge accros which t1 is | ||
| 717 | * adjacent to t2_in. | ||
| 718 | * \param[in] t1 first triangle | ||
| 719 | * \param[in] t2 second triangle | ||
| 720 | * \return e such that triangle_adjacent(t1,e)==t2 | ||
| 721 | * \pre \p t1 and \p t2 are adjacent | ||
| 722 | */ | ||
| 723 | 722 | index_t find_triangle_adjacent(index_t t1, index_t t2) const { | |
| 724 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 722 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
722 | geo_debug_assert(t1 < max_t()); |
| 725 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 722 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
722 | geo_debug_assert(t2 < max_t()); |
| 726 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 722 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
722 | geo_debug_assert(t1 != t2); |
| 727 | |||
| 728 | // Find local index of t2 in triangle t1 adajcent tets. | ||
| 729 | 722 | const index_t* T = &(cell_to_cell_store_[3 * t1]); | |
| 730 | 722 | index_t result = find_3(T,t2); | |
| 731 | |||
| 732 | // Sanity check: make sure that t1 is adjacent to t2 | ||
| 733 | // only once! | ||
| 734 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 722 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
722 | geo_debug_assert(triangle_adjacent(t1,(result+1)%3) != t2); |
| 735 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 722 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
722 | geo_debug_assert(triangle_adjacent(t1,(result+2)%3) != t2); |
| 736 | 722 | return result; | |
| 737 | } | ||
| 738 | |||
| 739 | /** | ||
| 740 | * \brief Sets the vertices and adjacent triangles of | ||
| 741 | * a triangle. | ||
| 742 | * \param[in] t index of the triangle | ||
| 743 | * \param[in] v0 index of the first vertex | ||
| 744 | * \param[in] v1 index of the second vertex | ||
| 745 | * \param[in] v2 index of the third vertex | ||
| 746 | * \param[in] a0 index of the adjacent triangle opposite to \p v0 | ||
| 747 | * \param[in] a1 index of the adjacent triangle opposite to \p v1 | ||
| 748 | * \param[in] a2 index of the adjacent triangle opposite to \p v2 | ||
| 749 | */ | ||
| 750 | void set_triangle( | ||
| 751 | index_t t, | ||
| 752 | index_t v0, index_t v1, index_t v2, | ||
| 753 | index_t a0, index_t a1, index_t a2 | ||
| 754 | ) { | ||
| 755 | geo_debug_assert(t < max_t()); | ||
| 756 | cell_to_v_store_[3 * t] = v0; | ||
| 757 | cell_to_v_store_[3 * t + 1] = v1; | ||
| 758 | cell_to_v_store_[3 * t + 2] = v2; | ||
| 759 | cell_to_cell_store_[3 * t] = a0; | ||
| 760 | cell_to_cell_store_[3 * t + 1] = a1; | ||
| 761 | cell_to_cell_store_[3 * t + 2] = a2; | ||
| 762 | } | ||
| 763 | |||
| 764 | /******* Predicates ******************************************/ | ||
| 765 | |||
| 766 | /** | ||
| 767 | * \brief Tests whether a given triangle is in conflict with | ||
| 768 | * a given 3d point. | ||
| 769 | * \details A real triangle is in conflict with a point whenever | ||
| 770 | * the point is contained by its circumscribed sphere, and a | ||
| 771 | * virtual triangle is in conflict with a point whenever the | ||
| 772 | * triangle formed by its real face and with the point has | ||
| 773 | * positive orientation. | ||
| 774 | * \param[in] t the index of the triangle | ||
| 775 | * \param[in] p a pointer to the coordinates of the point | ||
| 776 | * \retval true if point \p p is in conflict with triangle \p t | ||
| 777 | * \retval false otherwise | ||
| 778 | */ | ||
| 779 | 1120 | bool triangle_is_conflict(index_t t, const double* p) const { | |
| 780 | |||
| 781 | // Lookup triangle vertices | ||
| 782 | const double* pv[3]; | ||
| 783 |
2/2✓ Branch 0 taken 3360 times.
✓ Branch 1 taken 1120 times.
|
4480 | for(index_t i=0; i<3; ++i) { |
| 784 |
1/2✓ Branch 1 taken 3360 times.
✗ Branch 2 not taken.
|
3360 | index_t v = triangle_vertex(t,i); |
| 785 |
3/4✓ Branch 0 taken 2967 times.
✓ Branch 1 taken 393 times.
✓ Branch 3 taken 2967 times.
✗ Branch 4 not taken.
|
3360 | pv[i] = (v == NO_INDEX) ? nullptr : vertex_ptr(v); |
| 786 | } | ||
| 787 | |||
| 788 | // Check for virtual triangles (then in_circle() | ||
| 789 | // is replaced with orient2d()) | ||
| 790 |
2/2✓ Branch 0 taken 3132 times.
✓ Branch 1 taken 727 times.
|
3859 | for(index_t le = 0; le < 3; ++le) { |
| 791 | |||
| 792 |
2/2✓ Branch 0 taken 393 times.
✓ Branch 1 taken 2739 times.
|
3132 | if(pv[le] == nullptr) { |
| 793 | |||
| 794 | // Facet of a virtual triangle opposite to | ||
| 795 | // infinite vertex corresponds to | ||
| 796 | // the triangle on the convex hull of the points. | ||
| 797 | // Orientation is obtained by replacing vertex lf | ||
| 798 | // with p. | ||
| 799 | 393 | pv[le] = p; | |
| 800 |
1/2✓ Branch 1 taken 393 times.
✗ Branch 2 not taken.
|
393 | Sign sign = PCK::orient_2d(pv[0],pv[1],pv[2]); |
| 801 | |||
| 802 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 391 times.
|
393 | if(sign > 0) { |
| 803 | 2 | return true; | |
| 804 | } | ||
| 805 | |||
| 806 |
2/2✓ Branch 0 taken 229 times.
✓ Branch 1 taken 162 times.
|
391 | if(sign < 0) { |
| 807 | 229 | return false; | |
| 808 | } | ||
| 809 | |||
| 810 | // If sign is zero, we check the real triangle | ||
| 811 | // adjacent to the facet on the convex hull. | ||
| 812 |
2/8✓ Branch 1 taken 162 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 162 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
162 | geo_debug_assert(triangle_adjacent(t, le) != NO_INDEX); |
| 813 |
1/2✓ Branch 1 taken 162 times.
✗ Branch 2 not taken.
|
162 | index_t t2 = triangle_adjacent(t, le); |
| 814 |
2/8✓ Branch 1 taken 162 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 162 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
162 | geo_debug_assert(!triangle_is_virtual(t2)); |
| 815 | |||
| 816 | // If t2 is already chained in the conflict list, | ||
| 817 | // then it is conflict | ||
| 818 |
2/4✓ Branch 1 taken 162 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 162 times.
|
162 | if(triangle_is_in_list(t2)) { |
| 819 | ✗ | return true; | |
| 820 | } | ||
| 821 | |||
| 822 | // If t2 is marked, then it is not in conflict. | ||
| 823 |
3/4✓ Branch 1 taken 162 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 158 times.
|
162 | if(triangle_is_marked(t2)) { |
| 824 | 4 | return false; | |
| 825 | } | ||
| 826 | |||
| 827 |
1/2✓ Branch 1 taken 158 times.
✗ Branch 2 not taken.
|
158 | return triangle_is_conflict(t2, p); |
| 828 | } | ||
| 829 | } | ||
| 830 | |||
| 831 | // If the triangle is a finite one, it is in conflict | ||
| 832 | // if its circumscribed sphere contains the point (this is | ||
| 833 | // the standard case). | ||
| 834 | |||
| 835 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 727 times.
|
727 | if(weighted_) { |
| 836 | ✗ | double h0 = heights_[finite_triangle_vertex(t, 0)]; | |
| 837 | ✗ | double h1 = heights_[finite_triangle_vertex(t, 1)]; | |
| 838 | ✗ | double h2 = heights_[finite_triangle_vertex(t, 2)]; | |
| 839 | index_t pindex = index_t( | ||
| 840 | ✗ | (p - vertex_ptr(0)) / int(vertex_stride_) | |
| 841 | ✗ | ); | |
| 842 | ✗ | double h = heights_[pindex]; | |
| 843 | ✗ | return (PCK::orient_2dlifted_SOS( | |
| 844 | pv[0],pv[1],pv[2],p,h0,h1,h2,h | ||
| 845 | ✗ | ) > 0) ; | |
| 846 | } | ||
| 847 | |||
| 848 |
1/2✓ Branch 1 taken 727 times.
✗ Branch 2 not taken.
|
727 | return (PCK::in_circle_2d_SOS(pv[0], pv[1], pv[2], p) > 0); |
| 849 | } | ||
| 850 | |||
| 851 | protected: | ||
| 852 | |||
| 853 | /** | ||
| 854 | * \brief Finds the index of an integer in an array of three integers. | ||
| 855 | * \param[in] T a const pointer to an array of three integers | ||
| 856 | * \param[in] v the integer to retrieve in \p T | ||
| 857 | * \return the index (0,1 or 2) of \p v in \p T | ||
| 858 | * \pre The three entries of \p T are different and one of them is | ||
| 859 | * equal to \p v. | ||
| 860 | */ | ||
| 861 | 1404 | static inline index_t find_3(const index_t* T, index_t v) { | |
| 862 | // The following expression is 10% faster than using | ||
| 863 | // if() statements. This uses the C++ norm, that | ||
| 864 | // ensures that the 'true' boolean value converted to | ||
| 865 | // an int is always 1. With most compilers, this avoids | ||
| 866 | // generating branching instructions. | ||
| 867 | // Thank to Laurent Alonso for this idea. | ||
| 868 |
2/2✓ Branch 0 taken 390 times.
✓ Branch 1 taken 1014 times.
|
1404 | index_t result = index_t( (T[1] == v) | ((T[2] == v) * 2) ); |
| 869 | // Sanity check, important if it was T[0], not explicitly | ||
| 870 | // tested (detects input that does not meet the precondition). | ||
| 871 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 1404 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1404 | geo_debug_assert(T[result] == v); |
| 872 | 1404 | return result; | |
| 873 | } | ||
| 874 | |||
| 875 | /** | ||
| 876 | * \brief Delaunay2d destructor | ||
| 877 | */ | ||
| 878 | ~Delaunay2d() override; | ||
| 879 | |||
| 880 | /** | ||
| 881 | * \brief For debugging purposes, displays a triangle. | ||
| 882 | * \param[in] t index of the triangle to display. | ||
| 883 | */ | ||
| 884 | void show_triangle(index_t t) const; | ||
| 885 | |||
| 886 | /** | ||
| 887 | * \brief For debugging purposes, displays a triangle adjacency. | ||
| 888 | * \param[in] t index of the triangle to display. | ||
| 889 | * \param[in] le local index (0,1,2) of the triangle | ||
| 890 | * facet adjacenty to display. | ||
| 891 | */ | ||
| 892 | void show_triangle_adjacent(index_t t, index_t le) const; | ||
| 893 | |||
| 894 | /** | ||
| 895 | * \brief For debugging purposes, displays a triangle. | ||
| 896 | * \param[in] first index of the first triangle in the list | ||
| 897 | * \param[in] list_name name of the list, will be displayed as well | ||
| 898 | */ | ||
| 899 | void show_list(index_t first, const std::string& list_name) const; | ||
| 900 | |||
| 901 | /** | ||
| 902 | * \brief For debugging purposes, tests some combinatorial properties. | ||
| 903 | */ | ||
| 904 | void check_combinatorics(bool verbose = false) const; | ||
| 905 | |||
| 906 | /** | ||
| 907 | * \brief For debugging purposes, test some geometrical properties. | ||
| 908 | */ | ||
| 909 | void check_geometry(bool verbose = false) const; | ||
| 910 | |||
| 911 | private: | ||
| 912 | vector<index_t> cell_to_v_store_; | ||
| 913 | vector<index_t> cell_to_cell_store_; | ||
| 914 | vector<index_t> cell_next_; | ||
| 915 | vector<index_t> reorder_; | ||
| 916 | index_t cur_stamp_; // used for marking | ||
| 917 | index_t first_free_; | ||
| 918 | bool weighted_; | ||
| 919 | vector<double> heights_; // only used in weighted mode | ||
| 920 | |||
| 921 | /** | ||
| 922 | * Performs additional checks (costly !) | ||
| 923 | */ | ||
| 924 | bool debug_mode_; | ||
| 925 | |||
| 926 | /** | ||
| 927 | * Displays the result of the additional checks. | ||
| 928 | */ | ||
| 929 | bool verbose_debug_mode_; | ||
| 930 | |||
| 931 | /** | ||
| 932 | * Displays the timing of the core algorithm. | ||
| 933 | */ | ||
| 934 | bool benchmark_mode_; | ||
| 935 | |||
| 936 | /** | ||
| 937 | * \brief Gives the indexing of triangle edge | ||
| 938 | * vertices. | ||
| 939 | * \details triangle_edge_vertex[le][lv] gives the | ||
| 940 | * local vertex index (in 0,1,2) from a | ||
| 941 | * local edge index le (in 0,1,2) and a | ||
| 942 | * local vertex index within the edge (in 0,1). | ||
| 943 | */ | ||
| 944 | static char triangle_edge_vertex_[3][2]; | ||
| 945 | |||
| 946 | /** | ||
| 947 | * \brief Used by find_conflict_zone_iterative() | ||
| 948 | */ | ||
| 949 | std::stack<index_t> S_; | ||
| 950 | |||
| 951 | /** | ||
| 952 | * \brief Regular triangulations can have empty cells. | ||
| 953 | */ | ||
| 954 | bool has_empty_cells_; | ||
| 955 | |||
| 956 | /** | ||
| 957 | * \brief Stop inserting points as soon as an empty cell | ||
| 958 | * is encountered. | ||
| 959 | */ | ||
| 960 | bool abort_if_empty_cell_; | ||
| 961 | }; | ||
| 962 | |||
| 963 | /************************************************************************/ | ||
| 964 | |||
| 965 | /** | ||
| 966 | * \brief Regular Delaunay triangulation of weighted points | ||
| 967 | * \details | ||
| 968 | * - the input points are 2d points, were the third coordinate | ||
| 969 | * of point \f$ i \f$ is \f$ \sqrt{W - w_i} \f$ where \f$ W \f$ is | ||
| 970 | * the maximum of the weights of all the points and \d$ w_i \$ is | ||
| 971 | * the weight associated with vertex \f$ i \f$. | ||
| 972 | * - the constructed combinatorics is a triangulated surface (2d and | ||
| 973 | * not 3d although dimension() returns 3). This triangulated surface | ||
| 974 | * corresponds to the regular triangulation of the weighted points. | ||
| 975 | */ | ||
| 976 | class GEOGRAM_API RegularWeightedDelaunay2d : public Delaunay2d { | ||
| 977 | public: | ||
| 978 | /** | ||
| 979 | * \brief Constructs a new Regular Delaunay2d triangulation. | ||
| 980 | * \details RegularWeightedDelaunay2d triangulations are only | ||
| 981 | * supported for dimension 3. If a different dimension is specified in | ||
| 982 | * the constructor, a InvalidDimension exception is thrown. | ||
| 983 | * \param[in] dimension dimension of the triangulation | ||
| 984 | * \throw InvalidDimension This exception is thrown if dimension is | ||
| 985 | * different from 3. | ||
| 986 | */ | ||
| 987 | RegularWeightedDelaunay2d(coord_index_t dimension = 3); | ||
| 988 | |||
| 989 | protected: | ||
| 990 | /** | ||
| 991 | * \brief RegularWeightedDelaunay2d destructor | ||
| 992 | */ | ||
| 993 | ~RegularWeightedDelaunay2d() override; | ||
| 994 | }; | ||
| 995 | } | ||
| 996 | |||
| 997 | #endif | ||
| 998 |