| 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 | // Reference: S. W. Sloan, a fast algorithm for generating | ||
| 41 | // constrained Delaunay triangulation, 1992, | ||
| 42 | // Computers and Structures | ||
| 43 | // | ||
| 44 | // Specificities of this implementation: | ||
| 45 | // | ||
| 46 | // - Edges are systematically manipulated through triangles, | ||
| 47 | // and these triangles are rotated in-place in the mesh, | ||
| 48 | // in such a way that the edge we are talking about is | ||
| 49 | // systematically edge 0 (with vertices 1 and 2). | ||
| 50 | // | ||
| 51 | // - The constraint-enforcing step manipulates a queue Q | ||
| 52 | // of edges encoded this way. It examines pairs of triangles | ||
| 53 | // t1,t2=Tadj(t1,0), decides whether to swap their common | ||
| 54 | // edge (based on convexity test and intersection of | ||
| 55 | // t1's edge 0 with the constraint). In fact, this intersection | ||
| 56 | // test only depends on the combinatorics of (t1,t2) (two cases) | ||
| 57 | // and the position of t1' vertex 0 relative to the constraint | ||
| 58 | // (two cases), that makes 4 cases in total. In these cases, | ||
| 59 | // - t1 can either leave Q or be enqueued again | ||
| 60 | // - t2 was always in Q already (because it has an edge | ||
| 61 | // that has an intersection with the constraint), but | ||
| 62 | // there is one case where it leaves Q | ||
| 63 | // DList has an O(1) function to test whether an element is in the list (using | ||
| 64 | // flags associated with the elements). It is used in one case: when t2 is | ||
| 65 | // is not in Q, it means there is no intersection. | ||
| 66 | |||
| 67 | #include <geogram/delaunay/CDT_2d.h> | ||
| 68 | #include <geogram/mesh/mesh_reorder.h> | ||
| 69 | #include <geogram/basic/numeric.h> | ||
| 70 | #include <geogram/basic/boolean_expression.h> | ||
| 71 | |||
| 72 | #ifndef GEOGRAM_PSM | ||
| 73 | #include <geogram/mesh/mesh.h> | ||
| 74 | #include <geogram/mesh/mesh_io.h> | ||
| 75 | #endif | ||
| 76 | |||
| 77 | // Used by debugging functions and statistics | ||
| 78 | #include <geogram/mesh/index.h> | ||
| 79 | #include <geogram/basic/debug_stream.h> | ||
| 80 | #include <set> | ||
| 81 | #include <deque> | ||
| 82 | #include <stack> | ||
| 83 | //#define CDT_NAIVE // use naive per-edge method (kept for reference/debugging) | ||
| 84 | |||
| 85 | #ifdef GEO_DEBUG | ||
| 86 | //#define CDT_DEBUG // display *lots* of messages and activates costly checks | ||
| 87 | #endif | ||
| 88 | |||
| 89 | #ifdef CDT_DEBUG | ||
| 90 | #define CDT_LOG(X) std::cerr << X << std::endl | ||
| 91 | #else | ||
| 92 | #define CDT_LOG(X) | ||
| 93 | #endif | ||
| 94 | |||
| 95 | namespace GEO { | ||
| 96 | |||
| 97 | 652 | CDTBase2d::CDTBase2d() : | |
| 98 | 652 | nv_(0), | |
| 99 | 652 | ncnstr_(0), | |
| 100 | 652 | delaunay_(true), | |
| 101 | 652 | exact_incircle_(true), | |
| 102 | 652 | exact_intersections_(true) { | |
| 103 | 652 | } | |
| 104 | |||
| 105 |
2/2✓ Branch 0 taken 579 times.
✓ Branch 1 taken 73 times.
|
1304 | CDTBase2d::~CDTBase2d() { |
| 106 | 1304 | } | |
| 107 | |||
| 108 | 23572 | void CDTBase2d::clear() { | |
| 109 | 23572 | nv_ = 0; | |
| 110 | 23572 | ncnstr_ = 0; | |
| 111 | 23572 | T_.resize(0); | |
| 112 | 23572 | Tadj_.resize(0); | |
| 113 | 23572 | v2T_.resize(0); | |
| 114 | 23572 | Tflags_.resize(0); | |
| 115 | 23572 | Tecnstr_first_.resize(0); | |
| 116 | 23572 | ecnstr_val_.resize(0); | |
| 117 | 23572 | ecnstr_next_.resize(0); | |
| 118 | 23572 | Tnext_.resize(0); | |
| 119 | 23572 | Tprev_.resize(0); | |
| 120 | 23572 | } | |
| 121 | |||
| 122 | 15537 | void CDTBase2d::create_enclosing_triangle( | |
| 123 | index_t v0, index_t v1, index_t v2 | ||
| 124 | ) { | ||
| 125 | 15537 | nv_ = 3; | |
| 126 | 15537 | v2T_.resize(3); | |
| 127 | geo_debug_assert(v0 <= 3); | ||
| 128 | geo_debug_assert(v1 <= 3); | ||
| 129 | geo_debug_assert(v2 <= 3); | ||
| 130 | 15537 | index_t t0 = Tnew(); | |
| 131 | 15537 | Tset(t0, v0, v1, v2, NO_INDEX, NO_INDEX, NO_INDEX); | |
| 132 | 15537 | orient_012_ = orient2d(0,1,2); | |
| 133 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 15537 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
15537 | geo_assert(orient_012_ != ZERO); |
| 134 | 15537 | } | |
| 135 | |||
| 136 | 8173 | void CDTBase2d::create_enclosing_quad( | |
| 137 | index_t v0, index_t v1, index_t v2, index_t v3 | ||
| 138 | ) { | ||
| 139 | 8173 | nv_ = 4; | |
| 140 | 8173 | v2T_.resize(4); | |
| 141 | geo_debug_assert(v0 <= 4); | ||
| 142 | geo_debug_assert(v1 <= 4); | ||
| 143 | geo_debug_assert(v2 <= 4); | ||
| 144 | geo_debug_assert(v3 <= 4); | ||
| 145 | 8173 | index_t t0 = Tnew(); | |
| 146 | 8173 | index_t t1 = Tnew(); | |
| 147 | 8173 | Tset(t0, v0, v1, v3, t1, NO_INDEX, NO_INDEX); | |
| 148 | 8173 | Tset(t1, v3, v1, v2, NO_INDEX, NO_INDEX, t0); | |
| 149 | 8173 | orient_012_ = orient2d(0,1,2); | |
| 150 | geo_debug_assert(is_convex_quad(t0)); | ||
| 151 |
2/2✓ Branch 1 taken 6131 times.
✓ Branch 2 taken 2042 times.
|
8173 | if(Sign(incircle(v0,v1,v2,v3)*orient_012_) == POSITIVE) { |
| 152 | 6131 | swap_edge(t0); | |
| 153 | } | ||
| 154 | 8173 | } | |
| 155 | |||
| 156 | 222718 | void CDTBase2d::begin_insert_transaction() { | |
| 157 | 222718 | } | |
| 158 | |||
| 159 | 222712 | void CDTBase2d::commit_insert_transaction() { | |
| 160 | 222712 | } | |
| 161 | |||
| 162 | 6 | void CDTBase2d::rollback_insert_transaction() { | |
| 163 | 6 | } | |
| 164 | |||
| 165 |
2/2✓ Branch 0 taken 286740 times.
✓ Branch 1 taken 222718 times.
|
509458 | index_t CDTBase2d::insert(index_t v, index_t hint) { |
| 166 | bool keep_duplicates = false; | ||
| 167 |
2/2✓ Branch 0 taken 286740 times.
✓ Branch 1 taken 222718 times.
|
509458 | if(v == nv()) { |
| 168 |
2/2✓ Branch 0 taken 285426 times.
✓ Branch 1 taken 1314 times.
|
286740 | v2T_.push_back(NO_INDEX); |
| 169 | 286740 | ++nv_; | |
| 170 | } else { | ||
| 171 | // We are inserting a vertex in the middle of the | ||
| 172 | // list, which means we are doing batch-insertion, | ||
| 173 | // then we will not discard duplicates. | ||
| 174 | keep_duplicates = true; | ||
| 175 | geo_debug_assert(v < nv_); | ||
| 176 | } | ||
| 177 | |||
| 178 | // Phase 1: find triangle that contains vertex i | ||
| 179 | 509458 | begin_insert_transaction(); | |
| 180 | Sign o[3]; | ||
| 181 | 509458 | index_t t = locate(v,hint,o); | |
| 182 | 509458 | int nb_z = (o[0] == ZERO) + (o[1] == ZERO) + (o[2] == ZERO); | |
| 183 | geo_debug_assert(nb_z != 3); | ||
| 184 | |||
| 185 | // Duplicated vertex | ||
| 186 |
2/2✓ Branch 0 taken 151900 times.
✓ Branch 1 taken 357558 times.
|
509458 | if(nb_z == 2) { |
| 187 | CDT_LOG("duplicated vertex"); | ||
| 188 |
2/2✓ Branch 0 taken 63712 times.
✓ Branch 1 taken 88188 times.
|
151900 | v = (o[0] != ZERO) ? Tv(t,0) : |
| 189 |
2/2✓ Branch 0 taken 43616 times.
✓ Branch 1 taken 44572 times.
|
88188 | (o[1] != ZERO) ? Tv(t,1) : |
| 190 | Tv(t,2) ; | ||
| 191 |
2/2✓ Branch 0 taken 151894 times.
✓ Branch 1 taken 6 times.
|
151900 | if(!keep_duplicates) { |
| 192 | v2T_.pop_back(); | ||
| 193 | 151894 | --nv_; | |
| 194 | } | ||
| 195 | // Used by optional predicate cache management in derived classes. | ||
| 196 | // locate() computed some orient_2d() predicates, that may be | ||
| 197 | // stored in a temporary buffer, discard it. | ||
| 198 | 151900 | rollback_insert_transaction(); | |
| 199 | 151900 | return v; | |
| 200 | } | ||
| 201 | |||
| 202 | // Used by optional predicate cache management in derived classes. | ||
| 203 | // Copy the computed orient_2d() values to predicate cache. | ||
| 204 | 357558 | commit_insert_transaction(); | |
| 205 | |||
| 206 | // Stack of triangle edges to examine for flipping. Ignored in | ||
| 207 | // non-Delaunay mode (ignored if !delaunay_) | ||
| 208 | // Note: it is always edge 0 that we examine, since new | ||
| 209 | // triangles are always created with v as vertex 0. | ||
| 210 | DList S(*this); | ||
| 211 |
2/2✓ Branch 0 taken 246202 times.
✓ Branch 1 taken 111356 times.
|
357558 | if(delaunay_) { |
| 212 | S.initialize(DLIST_S_ID); | ||
| 213 | } | ||
| 214 | |||
| 215 | // Phase 2: split triangle | ||
| 216 | // Particular case: v is on edge | ||
| 217 |
2/2✓ Branch 0 taken 74642 times.
✓ Branch 1 taken 282916 times.
|
357558 | if(nb_z == 1) { |
| 218 | CDT_LOG("insert vertex on edge"); | ||
| 219 |
2/2✓ Branch 0 taken 34421 times.
✓ Branch 1 taken 40221 times.
|
74642 | index_t le = (o[0] == ZERO) ? 0 : |
| 220 |
2/2✓ Branch 0 taken 16232 times.
✓ Branch 1 taken 18189 times.
|
34421 | (o[1] == ZERO) ? 1 : |
| 221 | 2 ; | ||
| 222 |
1/2✓ Branch 1 taken 74642 times.
✗ Branch 2 not taken.
|
74642 | insert_vertex_in_edge(v,t,le,S); |
| 223 | } else { | ||
| 224 | CDT_LOG("insert vertex in triangle"); | ||
| 225 |
1/2✓ Branch 1 taken 282916 times.
✗ Branch 2 not taken.
|
282916 | insert_vertex_in_triangle(v,t,S); |
| 226 | } | ||
| 227 | |||
| 228 | // Phase 3: recursively restore Delaunay conditions for the neighbors | ||
| 229 | // of the new vertex | ||
| 230 |
2/2✓ Branch 0 taken 246202 times.
✓ Branch 1 taken 111356 times.
|
357558 | if(delaunay_) { |
| 231 |
1/2✓ Branch 1 taken 246202 times.
✗ Branch 2 not taken.
|
246202 | Delaunayize_vertex_neighbors(v,S); |
| 232 | } | ||
| 233 | |||
| 234 | #ifdef CDT_DEBUG | ||
| 235 | debug_check_consistency(); | ||
| 236 | #endif | ||
| 237 | return v; | ||
| 238 | 357558 | } | |
| 239 | |||
| 240 | |||
| 241 | 179318 | void CDTBase2d::insert_constraint(index_t i, index_t j) { | |
| 242 | geo_debug_assert(i < nv()); | ||
| 243 | geo_debug_assert(j < nv()); | ||
| 244 | CDT_LOG("insert constraint: " << i << "-" << j); | ||
| 245 | #ifdef CDT_DEBUG | ||
| 246 | debug_check_consistency(); | ||
| 247 | #endif | ||
| 248 | 179318 | ++ncnstr_; | |
| 249 | |||
| 250 | // Index of first vertex coming from constraints intersection | ||
| 251 | // (keep track of it to re-Delaunayize their neighborhoods). | ||
| 252 |
2/2✓ Branch 0 taken 178891 times.
✓ Branch 1 taken 427 times.
|
179318 | index_t first_v_isect = nv_; |
| 253 | |||
| 254 | #ifndef CDT_NAIVE | ||
| 255 | DList Q(*this, DLIST_Q_ID); // Queue of edges to constrain | ||
| 256 | DList N(*this); // New edges to re-Delaunayize (ignored if !delaunay_) | ||
| 257 |
2/2✓ Branch 0 taken 178891 times.
✓ Branch 1 taken 427 times.
|
179318 | if(delaunay_) { |
| 258 | N.initialize(DLIST_N_ID); | ||
| 259 | } | ||
| 260 |
2/2✓ Branch 0 taken 258522 times.
✓ Branch 1 taken 179318 times.
|
437840 | while(i != j) { |
| 261 | |||
| 262 | // Step 1: find all the edges that have an intersection | ||
| 263 | // with the constraint [i,j], enqueue them in Q. | ||
| 264 | // Stop at vertex on constraint or constraint intersection | ||
| 265 | // if any (returned in k) | ||
| 266 |
1/2✓ Branch 1 taken 258522 times.
✗ Branch 2 not taken.
|
258522 | index_t k = find_intersected_edges(i,j,Q); |
| 267 | |||
| 268 | // If we found a constraint intersection, | ||
| 269 | // we need to Delaunayize the neigborhood | ||
| 270 | // of the newly created vertex. Then we | ||
| 271 | // need to find the intersected edges again, | ||
| 272 | // since they may have changed. | ||
| 273 |
6/6✓ Branch 0 taken 258009 times.
✓ Branch 1 taken 513 times.
✓ Branch 2 taken 257496 times.
✓ Branch 3 taken 513 times.
✓ Branch 4 taken 5765 times.
✓ Branch 5 taken 251731 times.
|
258522 | if(delaunay_ && exact_intersections_ && k >= first_v_isect) { |
| 274 | geo_debug_assert(insert(k) == k); | ||
| 275 | Q.clear(); | ||
| 276 |
1/2✓ Branch 1 taken 5765 times.
✗ Branch 2 not taken.
|
5765 | Delaunayize_vertex_neighbors(k); |
| 277 | #ifdef CDT_DEBUG | ||
| 278 | debug_check_geometry(); | ||
| 279 | #endif | ||
| 280 |
1/2✓ Branch 1 taken 5765 times.
✗ Branch 2 not taken.
|
5765 | index_t new_k = find_intersected_edges(i,j,Q); |
| 281 |
1/8✗ Branch 0 not taken.
✓ Branch 1 taken 5765 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
5765 | geo_assert(new_k == k); |
| 282 | } | ||
| 283 | |||
| 284 | // Step 2: constrain edges | ||
| 285 |
1/2✓ Branch 1 taken 258522 times.
✗ Branch 2 not taken.
|
258522 | constrain_edges(i,k,Q,N); |
| 286 | |||
| 287 | // Step 3: restore Delaunay condition | ||
| 288 |
2/2✓ Branch 0 taken 258009 times.
✓ Branch 1 taken 513 times.
|
258522 | if(delaunay_) { |
| 289 |
1/2✓ Branch 1 taken 258009 times.
✗ Branch 2 not taken.
|
258009 | Delaunayize_new_edges(N); |
| 290 | #ifdef CDT_DEBUG | ||
| 291 | debug_check_geometry(); | ||
| 292 | #endif | ||
| 293 | } | ||
| 294 | |||
| 295 | i = k; | ||
| 296 | } | ||
| 297 | #else | ||
| 298 | DList Q(*this, DLIST_Q_ID); // Queue of edges to constrain | ||
| 299 | vector<Edge> N; // New edges to re-Delaunayize | ||
| 300 | while(i != j) { | ||
| 301 | index_t k = find_intersected_edges(i,j,Q); | ||
| 302 | |||
| 303 | // If we found a constraint intersection, | ||
| 304 | // we need to Delaunayize the neigborhood | ||
| 305 | // of the newly created vertex. Then we | ||
| 306 | // need to find the intersected edges again, | ||
| 307 | // since they may have changed. | ||
| 308 | if(delaunay_ && exact_intersections_ && k >= first_v_isect) { | ||
| 309 | geo_debug_assert(insert(k) == k); | ||
| 310 | Q.clear(); | ||
| 311 | Delaunayize_vertex_neighbors(k); | ||
| 312 | #ifdef CDT_DEBUG | ||
| 313 | debug_check_geometry(); | ||
| 314 | #endif | ||
| 315 | index_t new_k = find_intersected_edges(i,j,Q); | ||
| 316 | geo_assert(new_k == k); | ||
| 317 | } | ||
| 318 | |||
| 319 | constrain_edges_naive(i,k,Q,N); | ||
| 320 | debug_check_combinatorics(); | ||
| 321 | if(delaunay_) { | ||
| 322 | Delaunayize_new_edges_naive(N); | ||
| 323 | #ifdef CDT_DEBUG | ||
| 324 | debug_check_geometry(); | ||
| 325 | #endif | ||
| 326 | } | ||
| 327 | debug_check_combinatorics(); | ||
| 328 | i = k; | ||
| 329 | } | ||
| 330 | #endif | ||
| 331 | // Delaunayize neighborhood of vertices yielded by constraint | ||
| 332 | // intersections now if not done before (that is, if intersections | ||
| 333 | // are not exact, like in the default CDT2d class). If intersections | ||
| 334 | // are exact, it was done before (right after creating the intersection) | ||
| 335 |
4/4✓ Branch 0 taken 178891 times.
✓ Branch 1 taken 427 times.
✓ Branch 2 taken 427 times.
✓ Branch 3 taken 178464 times.
|
179318 | if(delaunay_ && !exact_intersections_) { |
| 336 |
2/2✓ Branch 0 taken 84 times.
✓ Branch 1 taken 427 times.
|
511 | for(index_t v=first_v_isect; v<nv(); ++v) { |
| 337 |
1/2✓ Branch 1 taken 84 times.
✗ Branch 2 not taken.
|
84 | Delaunayize_vertex_neighbors(v); |
| 338 | } | ||
| 339 | } | ||
| 340 | |||
| 341 | #ifdef CDT_DEBUG | ||
| 342 | debug_check_consistency(); | ||
| 343 | #endif | ||
| 344 | 179318 | } | |
| 345 | |||
| 346 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5849 times.
|
5849 | void CDTBase2d::Delaunayize_vertex_neighbors(index_t v) { |
| 347 | CDT_LOG("Delaunayize_vertex_neighbors " << v); | ||
| 348 | |||
| 349 | // Delaunayize triangles around vertices coming from | ||
| 350 | // constraint intersections | ||
| 351 | DList S(*this, DLIST_S_ID); | ||
| 352 | |||
| 353 |
1/8✗ Branch 0 not taken.
✓ Branch 1 taken 5849 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
5849 | geo_assert(vT(v) != NO_INDEX); |
| 354 | |||
| 355 | // We cannot use for_each_triangle_around_vertex() | ||
| 356 | // because we need to Trot() t during traveral, | ||
| 357 | // to have v has t's vertex 0 | ||
| 358 | // But the good news is that v is never on the border, | ||
| 359 | // (because it comes from an edge *intersection*), | ||
| 360 | // hence traversal is easier. | ||
| 361 | index_t t0 = vT(v); // Need to store it, because we Trot() | ||
| 362 | index_t t = t0; | ||
| 363 | do { | ||
| 364 | index_t lv = Tv_find(t,v); | ||
| 365 | 23858 | Trot(t,lv); | |
| 366 | geo_debug_assert(Tv(t,0) == v); | ||
| 367 | 23858 | S.push_back(t); | |
| 368 | t = Tadj(t, 1); | ||
| 369 |
1/8✗ Branch 0 not taken.
✓ Branch 1 taken 23858 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
23858 | geo_assert(t != NO_INDEX); |
| 370 |
2/2✓ Branch 0 taken 18009 times.
✓ Branch 1 taken 5849 times.
|
23858 | } while(t != t0); |
| 371 |
1/2✓ Branch 1 taken 5849 times.
✗ Branch 2 not taken.
|
5849 | Delaunayize_vertex_neighbors(v,S); |
| 372 | 5849 | } | |
| 373 | |||
| 374 | /** | ||
| 375 | * \brief Used by the implementation of find_intersected_edges() | ||
| 376 | * \details During traversal of a constrained edge [i,j], we can be on | ||
| 377 | * a vertex (then v != NO_INDEX) or | ||
| 378 | * on a triangle (then t != NO_INDEX). | ||
| 379 | * We also keep track of the previous vertex (prev_v) and previous | ||
| 380 | * triangle(prev_t) in order to make sure we do not go backwards. | ||
| 381 | */ | ||
| 382 | struct CDT2d_ConstraintWalker { | ||
| 383 | /** | ||
| 384 | * \brief ConstraintWalker constructor | ||
| 385 | * \param[in] i_in , j_in extremities of the constrained edge | ||
| 386 | */ | ||
| 387 | 264287 | CDT2d_ConstraintWalker(index_t i_in, index_t j_in) : | |
| 388 | 264287 | i(i_in), j(j_in), | |
| 389 | 264287 | t_prev(NO_INDEX), v_prev(NO_INDEX), | |
| 390 | 264287 | t(NO_INDEX), v(i_in), | |
| 391 | 264287 | v_cnstr(NO_INDEX) | |
| 392 | { | ||
| 393 | 264287 | } | |
| 394 | index_t i, j; | ||
| 395 | index_t t_prev, v_prev; | ||
| 396 | index_t t, v; | ||
| 397 | index_t v_cnstr; | ||
| 398 | }; | ||
| 399 | |||
| 400 | 264287 | index_t CDTBase2d::find_intersected_edges(index_t i, index_t j, DList& Q) { | |
| 401 | CDT_LOG("Find intersected edges: " << i << "-" << j); | ||
| 402 | CDT2d_ConstraintWalker W(i,j); | ||
| 403 | // Stop at the first encountered vertex or constraint intersection. | ||
| 404 |
4/4✓ Branch 0 taken 264287 times.
✓ Branch 1 taken 324236 times.
✓ Branch 2 taken 59949 times.
✓ Branch 3 taken 264287 times.
|
588523 | while(W.v == i || W.v == NO_INDEX) { |
| 405 | CDT_LOG( | ||
| 406 | " t=" << int(W.t) << " v=" << int(W.v) << " " | ||
| 407 | "t_prev=" << int(W.t_prev) << " v_prev=" << int(W.v_prev) | ||
| 408 | << " " | ||
| 409 | ); | ||
| 410 |
2/2✓ Branch 0 taken 264287 times.
✓ Branch 1 taken 59949 times.
|
324236 | if(W.v != NO_INDEX) { |
| 411 | 264287 | walk_constraint_v(W); | |
| 412 | } else { | ||
| 413 | 59949 | walk_constraint_t(W,Q); | |
| 414 | } | ||
| 415 | } | ||
| 416 | 264287 | return W.v; | |
| 417 | } | ||
| 418 | |||
| 419 | // The two functions below are more complicated than I wished, but are | ||
| 420 | // simpler than it looks like. There are two main different cases: | ||
| 421 | // | ||
| 422 | // - walk_constraint_v(): we are on a vertex. | ||
| 423 | // traverse all the triangles around v and find the one that | ||
| 424 | // has an intersection. For instance, when we start from vertex i, | ||
| 425 | // and also when the previous step encountered a vertex exactly on | ||
| 426 | // the constrained segment. It is the annoying case where one has | ||
| 427 | // to traverse the triangles incident to v (using the function | ||
| 428 | // for_each_T_around_v() that takes a lambda). | ||
| 429 | // | ||
| 430 | // - walk_constraint_t(): we are on an edge intersection. | ||
| 431 | // propagate to the neighbor of t accross the intersected edge. | ||
| 432 | // It is the "generic" case, simpler (the next triangle is | ||
| 433 | // determined by the edge of t that is intersected). | ||
| 434 | // | ||
| 435 | // There are three things that makes things slightly | ||
| 436 | // more complicated: | ||
| 437 | // - each case has two sub-cases, depending on whether the next | ||
| 438 | // intersection is an existing vertex. | ||
| 439 | // - if an existing edge is embedded in the constraint, one needs | ||
| 440 | // to flag that edge as a constraint. | ||
| 441 | // - we need to test whether we are arrived at vertex j | ||
| 442 | |||
| 443 | 264287 | void CDTBase2d::walk_constraint_v(CDT2d_ConstraintWalker& W) { | |
| 444 | geo_debug_assert(W.v != NO_INDEX); | ||
| 445 | geo_debug_assert(W.t == NO_INDEX); | ||
| 446 | |||
| 447 | 264287 | index_t t_next = NO_INDEX; | |
| 448 | 264287 | index_t v_next = NO_INDEX; | |
| 449 | |||
| 450 |
1/2✓ Branch 1 taken 264287 times.
✗ Branch 2 not taken.
|
264287 | for_each_T_around_v( |
| 451 | 782848 | W.v, [&](index_t t_around_v, index_t le) { | |
| 452 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 518561 times.
|
518561 | if(t_around_v == W.t_prev) { // Don't go backwards ! |
| 453 | return false; | ||
| 454 | } | ||
| 455 |
2/2✓ Branch 0 taken 429424 times.
✓ Branch 1 taken 89137 times.
|
518561 | index_t v1 = Tv(t_around_v, (le + 1)%3); |
| 456 | 518561 | index_t v2 = Tv(t_around_v, (le + 2)%3); | |
| 457 |
4/4✓ Branch 0 taken 429424 times.
✓ Branch 1 taken 89137 times.
✓ Branch 2 taken 73023 times.
✓ Branch 3 taken 356401 times.
|
518561 | if(v1 == W.j || v2 == W.j) { // Are we arrived at j ? |
| 458 | 162160 | v_next = W.j; | |
| 459 | // Edge is flagged as constraint here, because | ||
| 460 | // it will not be seen by constraint enforcement. | ||
| 461 |
2/2✓ Branch 0 taken 73023 times.
✓ Branch 1 taken 89137 times.
|
162160 | index_t le_cnstr_edge = (v1 == W.j) ? (le+2)%3 : (le+1)%3; |
| 462 | 162160 | Tadd_edge_cnstr_with_neighbor( | |
| 463 | 162160 | t_around_v, le_cnstr_edge, ncnstr_-1 | |
| 464 | ); | ||
| 465 | CDT_LOG( | ||
| 466 | " During cnstr " << W.i << "-" << W.j << ": " | ||
| 467 | << " Constrained edge " | ||
| 468 | << Tv(t_around_v, (le_cnstr_edge+1)%3) << "-" | ||
| 469 | << Tv(t_around_v, (le_cnstr_edge+2)%3) | ||
| 470 | ); | ||
| 471 | 162160 | return true; | |
| 472 | } | ||
| 473 | 356401 | Sign o1 = orient2d(W.i,W.j,v1); | |
| 474 | 356401 | Sign o2 = orient2d(W.i,W.j,v2); | |
| 475 | 356401 | Sign o3 = orient2d(v1,v2,W.j); | |
| 476 | 356401 | Sign o4 = orient_012_; // equivalent to orient2d(v1,v2,i) | |
| 477 |
4/4✓ Branch 0 taken 40775 times.
✓ Branch 1 taken 315626 times.
✓ Branch 2 taken 17024 times.
✓ Branch 3 taken 23751 times.
|
356401 | if(o1*o2 < 0 && o3*o4 < 0) { |
| 478 | 17024 | Trot(t_around_v,le); // so that le becomes edge 0 | |
| 479 | 17024 | t_next = t_around_v; // will be added to Q during next round | |
| 480 | 17024 | return true; | |
| 481 | } else { | ||
| 482 | // Special case: v1 or v2 is exactly on [i,j] | ||
| 483 | // Edge is flagged as constraint here, because | ||
| 484 | // it will not be seen by constraint enforcement. | ||
| 485 | geo_debug_assert(o1 != ZERO || o2 != ZERO); | ||
| 486 |
5/6✓ Branch 0 taken 91252 times.
✓ Branch 1 taken 248125 times.
✓ Branch 2 taken 78539 times.
✓ Branch 3 taken 12713 times.
✓ Branch 4 taken 78539 times.
✗ Branch 5 not taken.
|
339377 | if(o1 == ZERO && o3*o4 < 0 && v1 != W.v_prev) { |
| 487 | 78539 | v_next = v1; | |
| 488 | 78539 | Tadd_edge_cnstr_with_neighbor( | |
| 489 | 78539 | t_around_v, (le + 2)%3, ncnstr_-1 | |
| 490 | ); | ||
| 491 | 78539 | return true; | |
| 492 |
5/6✓ Branch 0 taken 206005 times.
✓ Branch 1 taken 54833 times.
✓ Branch 2 taken 48269 times.
✓ Branch 3 taken 6564 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 6564 times.
|
260838 | } else if(o2 == ZERO && o3*o4 < 0 && v2 != W.v_prev) { |
| 493 | 6564 | v_next = v2; | |
| 494 | 6564 | Tadd_edge_cnstr_with_neighbor( | |
| 495 | 6564 | t_around_v, (le + 1)%3, ncnstr_-1 | |
| 496 | ); | ||
| 497 | 6564 | return true; | |
| 498 | } | ||
| 499 | } | ||
| 500 | return false; | ||
| 501 | } | ||
| 502 | ); // End of for_each_T_around_v() loop | ||
| 503 | 264287 | W.t_prev = W.t; | |
| 504 | 264287 | W.v_prev = W.v; | |
| 505 | 264287 | W.t = t_next; | |
| 506 | 264287 | W.v = v_next; | |
| 507 | 264287 | } | |
| 508 | |||
| 509 | 59949 | void CDTBase2d::walk_constraint_t(CDT2d_ConstraintWalker& W, DList& Q) { | |
| 510 | geo_debug_assert(W.v == NO_INDEX); | ||
| 511 | geo_debug_assert(W.t != NO_INDEX); | ||
| 512 | |||
| 513 | index_t v_next = NO_INDEX; | ||
| 514 | index_t t_next = NO_INDEX; | ||
| 515 | |||
| 516 |
6/6✓ Branch 0 taken 57602 times.
✓ Branch 1 taken 2347 times.
✓ Branch 2 taken 53729 times.
✓ Branch 3 taken 3873 times.
✓ Branch 4 taken 50005 times.
✓ Branch 5 taken 3724 times.
|
59949 | if(Tv(W.t,0) == W.j || Tv(W.t,1) == W.j || Tv(W.t,2) == W.j) { |
| 517 | v_next = W.j; // Are we arrived at j ? | ||
| 518 | } else { | ||
| 519 | // Test the three edges of the triangle | ||
| 520 |
1/2✓ Branch 0 taken 78176 times.
✗ Branch 1 not taken.
|
78176 | for(index_t le = 0; le<3; ++le) { |
| 521 |
2/2✓ Branch 0 taken 15776 times.
✓ Branch 1 taken 62400 times.
|
78176 | if(Tadj(W.t,le) == W.t_prev) { // Do not go backwards ! |
| 522 | 15776 | continue; | |
| 523 | } | ||
| 524 | // Test whether [v1,v2] intersects the support line of (i,j). | ||
| 525 | // No need to test the *segment* [i,j]: we know the line enters | ||
| 526 | // the triangle (it is how we came here), and we know it leaves | ||
| 527 | // it, else j would have been one of the triangle's vertices. | ||
| 528 | 62400 | index_t v1 = Tv(W.t, (le + 1)%3); | |
| 529 | 62400 | index_t v2 = Tv(W.t, (le + 2)%3); | |
| 530 | 62400 | Sign o1 = orient2d(W.i,W.j,v1); | |
| 531 | 62400 | Sign o2 = orient2d(W.i,W.j,v2); | |
| 532 |
2/2✓ Branch 0 taken 48858 times.
✓ Branch 1 taken 13542 times.
|
62400 | if(o1*o2 < 0) { |
| 533 | // [v1,v2] has a frank intersection with [i,j] | ||
| 534 | 48858 | Trot(W.t,le); // So that edge 0 is intersected edge | |
| 535 |
2/2✓ Branch 0 taken 5933 times.
✓ Branch 1 taken 42925 times.
|
48858 | if(Tedge_is_constrained(W.t,0)) { |
| 536 | CDT_LOG(" Cnstr isect with:" << v1 << "-" << v2); | ||
| 537 | 5933 | v_next = create_intersection( | |
| 538 | ncnstr()-1, W.i, W.j, | ||
| 539 | edge_cnstr(Tedge_cnstr_first(W.t,0)), v1, v2 | ||
| 540 | ); | ||
| 541 | 5933 | insert_vertex_in_edge(v_next,W.t,0); | |
| 542 | // Mark new edge as constraint if walker was previously | ||
| 543 | // on a vertex. | ||
| 544 |
2/2✓ Branch 0 taken 3421 times.
✓ Branch 1 taken 2512 times.
|
5933 | if(W.v_prev != NO_INDEX) { |
| 545 | 3421 | Tadd_edge_cnstr_with_neighbor(W.t,2,ncnstr_-1); | |
| 546 | } | ||
| 547 | } else { | ||
| 548 | CDT_LOG(" Isect: t=" << W.t <<" E=" << v1 <<"-"<< v2); | ||
| 549 | 42925 | Q.push_back(W.t); | |
| 550 | 42925 | t_next = Tadj(W.t,0); | |
| 551 | } | ||
| 552 | break; | ||
| 553 | } else { // Special case: v1 or v2 is exactly on [i,j] | ||
| 554 | geo_debug_assert(o1 != ZERO || o2 != ZERO); | ||
| 555 |
2/2✓ Branch 0 taken 13527 times.
✓ Branch 1 taken 15 times.
|
13542 | if(o1 == ZERO) { |
| 556 | v_next = v1; | ||
| 557 | break; | ||
| 558 |
2/2✓ Branch 0 taken 12395 times.
✓ Branch 1 taken 1132 times.
|
13527 | } else if(o2 == ZERO) { |
| 559 | v_next = v2; | ||
| 560 | break; | ||
| 561 | } | ||
| 562 | } | ||
| 563 | } | ||
| 564 | } | ||
| 565 | 59949 | W.t_prev = W.t; | |
| 566 | 59949 | W.v_prev = W.v; | |
| 567 | 59949 | W.t = t_next; | |
| 568 | 59949 | W.v = v_next; | |
| 569 | 59949 | } | |
| 570 | |||
| 571 | 258522 | void CDTBase2d::constrain_edges(index_t i, index_t j, DList& Q, DList& N) { | |
| 572 | |||
| 573 | #ifdef CDT_DEBUG | ||
| 574 | // The function find_edge_intersections() is super complicated, | ||
| 575 | // so in debug mode I make sure it did its job correctly (by testing | ||
| 576 | // *all* edge intersections). | ||
| 577 | check_edge_intersections(i,j,Q); | ||
| 578 | #endif | ||
| 579 | // Called each time edge le of triangle t has no isect with cnstr, | ||
| 580 | // (then it is a "new edge") | ||
| 581 | 38957 | auto new_edge = [&](index_t t,index_t le) { | |
| 582 | 38957 | Trot(t,le); | |
| 583 | if( | ||
| 584 |
4/4✓ Branch 0 taken 12779 times.
✓ Branch 1 taken 26178 times.
✓ Branch 2 taken 1543 times.
✓ Branch 3 taken 11236 times.
|
38957 | (Tv(t,1) == i && Tv(t,2) == j) || |
| 585 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 27721 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
27721 | (Tv(t,1) == j && Tv(t,2) == i) |
| 586 | ) { | ||
| 587 | // Set constraint flag if the new edge is the constrained edge | ||
| 588 | 11236 | Tadd_edge_cnstr_with_neighbor(t,0,ncnstr_-1); | |
| 589 | } else { | ||
| 590 | // Memorize new edge as "to be Delaunayized" | ||
| 591 |
2/2✓ Branch 0 taken 11353 times.
✓ Branch 1 taken 16368 times.
|
27721 | if(N.initialized()) { |
| 592 | 11353 | N.push_back(t); | |
| 593 | } | ||
| 594 | } | ||
| 595 | 297479 | }; | |
| 596 | |||
| 597 | // Called each time edge le of triangle t still has an isect with cnstr | ||
| 598 | // (then it is queued again) | ||
| 599 | auto isect_edge = [&](index_t t, index_t le) { | ||
| 600 | 18292 | Trot(t,le); | |
| 601 | 18292 | Q.push_front(t); | |
| 602 | 18292 | }; | |
| 603 | |||
| 604 |
2/2✓ Branch 0 taken 78602 times.
✓ Branch 1 taken 258522 times.
|
337124 | while(!Q.empty()) { |
| 605 | 78602 | index_t t1 = Q.pop_back(); | |
| 606 |
2/2✓ Branch 1 taken 21353 times.
✓ Branch 2 taken 57249 times.
|
78602 | if(!is_convex_quad(t1)) { |
| 607 | // Sanity check: if the only remaining edge to flip does | ||
| 608 | // not form a convex quad, it means we are going to | ||
| 609 | // flip forever ! (shoud not happen) | ||
| 610 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 21353 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
21353 | geo_assert(!Q.empty()); |
| 611 | 21353 | Q.push_front(t1); | |
| 612 | } else { | ||
| 613 | index_t t2 = Tadj(t1,0); | ||
| 614 | bool no_isect = !Q.contains(t2); | ||
| 615 | index_t v0 = Tv(t1,0); | ||
| 616 |
4/4✓ Branch 0 taken 38289 times.
✓ Branch 1 taken 18960 times.
✓ Branch 2 taken 19367 times.
✓ Branch 3 taken 18922 times.
|
57249 | bool t2v0_t1v2 = (Q.contains(t2) && Tv(t2,0) == Tv(t1,2)); |
| 617 | bool t2v0_t1v1 = (Q.contains(t2) && Tv(t2,0) == Tv(t1,1)); | ||
| 618 | geo_argused(t2v0_t1v1); | ||
| 619 | |||
| 620 |
2/2✓ Branch 0 taken 18960 times.
✓ Branch 1 taken 38289 times.
|
57249 | if(no_isect) { |
| 621 | 18960 | swap_edge(t1); | |
| 622 | geo_debug_assert(!segment_edge_intersect(i,j,t1,2)); | ||
| 623 | 18960 | new_edge(t1,2); | |
| 624 | } else { | ||
| 625 | // See comment at beginning of file | ||
| 626 | // (a small variation in Sloan's | ||
| 627 | // method that makes better use of the combinatorics) | ||
| 628 | 38289 | Sign o = Sign(orient2d(i,j,v0) * orient_012_); | |
| 629 |
2/2✓ Branch 0 taken 18922 times.
✓ Branch 1 taken 19367 times.
|
38289 | if(t2v0_t1v2) { |
| 630 | 18922 | swap_edge(t1,false); // "new t1 on top" | |
| 631 |
2/2✓ Branch 0 taken 9630 times.
✓ Branch 1 taken 9292 times.
|
18922 | if(o >= 0) { |
| 632 | geo_debug_assert(!segment_edge_intersect(i,j,t1,2)); | ||
| 633 | geo_debug_assert( segment_edge_intersect(i,j,t2,0)); | ||
| 634 | 9630 | new_edge(t1,2); | |
| 635 | } else { | ||
| 636 | geo_debug_assert( segment_edge_intersect(i,j,t1,2)); | ||
| 637 | geo_debug_assert( segment_edge_intersect(i,j,t2,0)); | ||
| 638 | isect_edge(t1,2); | ||
| 639 | } | ||
| 640 | } else { | ||
| 641 | geo_debug_assert(t2v0_t1v1); | ||
| 642 | 19367 | swap_edge(t1,true); // "new t1 on bottom" | |
| 643 |
2/2✓ Branch 0 taken 9000 times.
✓ Branch 1 taken 10367 times.
|
19367 | if(o > 0) { |
| 644 | geo_debug_assert( segment_edge_intersect(i,j,t1,1)); | ||
| 645 | geo_debug_assert( segment_edge_intersect(i,j,t2,0)); | ||
| 646 | isect_edge(t1,1); | ||
| 647 | } else { | ||
| 648 | geo_debug_assert(!segment_edge_intersect(i,j,t1,1)); | ||
| 649 | geo_debug_assert( segment_edge_intersect(i,j,t2,0)); | ||
| 650 | 10367 | new_edge(t1,1); | |
| 651 | } | ||
| 652 | } | ||
| 653 | } | ||
| 654 | } | ||
| 655 | } | ||
| 656 | 258522 | } | |
| 657 | |||
| 658 | 252051 | void CDTBase2d::Delaunayize_vertex_neighbors(index_t v, DList& S) { | |
| 659 | CDT_LOG("Delaunayize_vertex_neighbors"); | ||
| 660 | index_t count = 0; | ||
| 661 |
2/2✓ Branch 0 taken 1704732 times.
✓ Branch 1 taken 252051 times.
|
1956783 | while(!S.empty()) { |
| 662 | // NASA programming style: all loops have | ||
| 663 | // a maximum number of iterations | ||
| 664 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1704732 times.
|
1704732 | ++count; |
| 665 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1704732 times.
|
1704732 | if(count > 10*nT()) { |
| 666 | ✗ | Logger::warn("CDT2d") | |
| 667 | << "Emergency exit in Delaunayize_vertex_neighbors()" | ||
| 668 | << std::endl; | ||
| 669 | S.clear(); | ||
| 670 | ✗ | geo_assert_not_reached; // For now, assert fail. | |
| 671 | } | ||
| 672 | 1704732 | index_t t1 = S.pop_back(); | |
| 673 | geo_debug_assert(Tv(t1,0) == v); | ||
| 674 |
2/2✓ Branch 0 taken 9790 times.
✓ Branch 1 taken 1694942 times.
|
1704732 | if(Tedge_is_constrained(t1,0)) { |
| 675 | 9790 | continue; | |
| 676 | } | ||
| 677 | index_t t2 = Tadj(t1,0); | ||
| 678 |
2/2✓ Branch 0 taken 204097 times.
✓ Branch 1 taken 1490845 times.
|
1694942 | if(t2 == NO_INDEX) { |
| 679 | 204097 | continue; | |
| 680 | } | ||
| 681 |
3/4✓ Branch 0 taken 1490845 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 230778 times.
✓ Branch 4 taken 1260067 times.
|
1490845 | if(!exact_incircle_ && !is_convex_quad(t1)) { |
| 682 | 230778 | continue; | |
| 683 | } | ||
| 684 | |||
| 685 | index_t v1 = Tv(t2,0); | ||
| 686 | index_t v2 = Tv(t2,1); | ||
| 687 | index_t v3 = Tv(t2,2); | ||
| 688 |
2/2✓ Branch 1 taken 506166 times.
✓ Branch 2 taken 753901 times.
|
1260067 | if(Sign(incircle(v1,v2,v3,v)*orient_012_) == POSITIVE) { |
| 689 | 506166 | swap_edge(t1); | |
| 690 | 506166 | S.push_back(t1); | |
| 691 | 506166 | S.push_back(t2); | |
| 692 | } | ||
| 693 | } | ||
| 694 | CDT_LOG("/Delaunayize_vertex_neighbors"); | ||
| 695 | 252051 | } | |
| 696 | |||
| 697 | 258009 | void CDTBase2d::Delaunayize_new_edges(DList& N) { | |
| 698 | index_t count = 0; | ||
| 699 | bool swap_occured = true; | ||
| 700 |
2/2✓ Branch 0 taken 259998 times.
✓ Branch 1 taken 258009 times.
|
518007 | while(swap_occured) { |
| 701 | swap_occured = false; | ||
| 702 | // NASA programming style: all loops have | ||
| 703 | // a maximum number of iterations | ||
| 704 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 259998 times.
|
259998 | ++count; |
| 705 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 259998 times.
|
259998 | if(count > 10*nT()) { |
| 706 | ✗ | Logger::warn("CDT2d") | |
| 707 | << "Emergency exit in Delaunayize_new_edges()" | ||
| 708 | << std::endl; | ||
| 709 | ✗ | break; | |
| 710 | } | ||
| 711 |
2/2✓ Branch 0 taken 23837 times.
✓ Branch 1 taken 259998 times.
|
283835 | for(index_t t1 = N.front(); t1 != NO_INDEX; t1 = N.next(t1)) { |
| 712 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 23837 times.
|
23837 | if(Tedge_is_constrained(t1,0)) { |
| 713 | ✗ | continue; | |
| 714 | } | ||
| 715 | index_t v1 = Tv(t1,1); | ||
| 716 | index_t v2 = Tv(t1,2); | ||
| 717 | index_t v0 = Tv(t1,0); | ||
| 718 | index_t t2 = Tadj(t1,0); | ||
| 719 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 23837 times.
|
23837 | if(t2 == NO_INDEX) { |
| 720 | ✗ | continue; | |
| 721 | } | ||
| 722 |
3/4✓ Branch 0 taken 23837 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 12865 times.
✓ Branch 4 taken 10972 times.
|
23837 | if(!exact_incircle_ && !is_convex_quad(t1)) { |
| 723 | 12865 | continue; | |
| 724 | } | ||
| 725 | index_t e2 = Tadj_find(t2,t1); | ||
| 726 | index_t v3 = Tv(t2,e2); | ||
| 727 |
2/2✓ Branch 1 taken 7988 times.
✓ Branch 2 taken 2984 times.
|
10972 | if(Sign(incircle(v0,v1,v2,v3)*orient_012_) == POSITIVE) { |
| 728 | // t2 may also encode a new edge, we need to preserve it, | ||
| 729 | // by chosing the right swap: | ||
| 730 |
2/2✓ Branch 0 taken 2463 times.
✓ Branch 1 taken 521 times.
|
2984 | if(Tv(t2,0) == Tv(t1,1)) { |
| 731 | 2463 | swap_edge(t1, true); // t2 on top | |
| 732 | 2463 | Trot(t1,1); | |
| 733 | } else { | ||
| 734 | 521 | swap_edge(t1, false); // t1 on top | |
| 735 | 521 | Trot(t1,2); | |
| 736 | } | ||
| 737 | swap_occured = true; | ||
| 738 | } | ||
| 739 | } | ||
| 740 | } | ||
| 741 | N.clear(); | ||
| 742 | 258009 | } | |
| 743 | |||
| 744 | 509458 | index_t CDTBase2d::locate(index_t v, index_t hint, Sign* o) const { | |
| 745 | Sign o_local[3]; | ||
| 746 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 509458 times.
|
509458 | if(o == nullptr) { |
| 747 | o = o_local; | ||
| 748 | } | ||
| 749 | |||
| 750 | // Efficient locate, "walking the triangulation" | ||
| 751 | 509458 | index_t t_pred = nT()+1; // Needs to be different from NO_INDEX | |
| 752 |
2/2✓ Branch 0 taken 286772 times.
✓ Branch 1 taken 222686 times.
|
509458 | index_t t = (hint == NO_INDEX) ? |
| 753 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 286772 times.
|
286772 | index_t(Numeric::random_int32()) % nT() : |
| 754 | 286772 | hint ; | |
| 755 | #ifdef GEO_DEBUG | ||
| 756 | index_t nb_traversed_t = 0; | ||
| 757 | #endif | ||
| 758 | |||
| 759 | still_walking: | ||
| 760 | { | ||
| 761 | #ifdef GEO_DEBUG | ||
| 762 | ++nb_traversed_t; | ||
| 763 | #endif | ||
| 764 | |||
| 765 | // Infinite loop are not supposed to happen, but | ||
| 766 | // let us detect them, just in case... | ||
| 767 | geo_debug_assert(nb_traversed_t <= 2*nT()); | ||
| 768 | |||
| 769 | // You will land here if we try to locate a point outside | ||
| 770 | // the boundary | ||
| 771 | bool point_outside_boundary = (t == NO_INDEX); | ||
| 772 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 6058859 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
6058859 | geo_assert(!point_outside_boundary); |
| 773 | |||
| 774 | index_t tv[3]; | ||
| 775 | 6058859 | tv[0] = Tv(t,0); | |
| 776 | 6058859 | tv[1] = Tv(t,1); | |
| 777 | 6058859 | tv[2] = Tv(t,2); | |
| 778 | |||
| 779 | // Start from a random edge | ||
| 780 | 6058859 | index_t e0 = index_t(Numeric::random_int32()) % 3; | |
| 781 |
2/2✓ Branch 0 taken 11751128 times.
✓ Branch 1 taken 509458 times.
|
12260586 | for(index_t de = 0; de < 3; ++de) { |
| 782 |
2/2✓ Branch 0 taken 2817626 times.
✓ Branch 1 taken 8933502 times.
|
11751128 | index_t le = (e0 + de) % 3; |
| 783 | |||
| 784 | index_t t_next = Tadj(t,le); | ||
| 785 | |||
| 786 | // If the candidate next triangle is the | ||
| 787 | // one we came from, then we know already that | ||
| 788 | // the orientation is positive, thus we examine | ||
| 789 | // the next candidate (or exit the loop if they | ||
| 790 | // are exhausted). | ||
| 791 | // | ||
| 792 | // (here is why intial value of t_pred needs to be | ||
| 793 | // different from NO_INDEX) | ||
| 794 |
2/2✓ Branch 0 taken 2817626 times.
✓ Branch 1 taken 8933502 times.
|
11751128 | if(t_next == t_pred) { |
| 795 | 2817626 | o[le] = POSITIVE; | |
| 796 | 2817626 | continue ; | |
| 797 | } | ||
| 798 | |||
| 799 | // To test the orientation of p w.r.t. the facet f of | ||
| 800 | // t, we replace vertex number f with p in t (same | ||
| 801 | // convention as in CGAL). | ||
| 802 | 8933502 | index_t v_bkp = tv[le]; | |
| 803 | 8933502 | tv[le] = v; | |
| 804 | 8933502 | o[le] = Sign(orient_012_ * orient2d(tv[0], tv[1], tv[2])); | |
| 805 | |||
| 806 | // If the orientation is not negative, then we cannot | ||
| 807 | // walk towards t_next, and examine the next candidate | ||
| 808 | // (or exit the loop if they are exhausted). | ||
| 809 |
2/2✓ Branch 0 taken 3384101 times.
✓ Branch 1 taken 5549401 times.
|
8933502 | if(o[le] != NEGATIVE) { |
| 810 | 3384101 | tv[le] = v_bkp; | |
| 811 | 3384101 | continue; | |
| 812 | } | ||
| 813 | |||
| 814 | // If we reach this point, then t_next is a valid | ||
| 815 | // successor, thus we are still walking. | ||
| 816 | t_pred = t; | ||
| 817 | t = t_next; | ||
| 818 | 5549401 | goto still_walking; | |
| 819 | } | ||
| 820 | } | ||
| 821 | |||
| 822 | 509458 | return t; | |
| 823 | } | ||
| 824 | |||
| 825 | 8071 | void CDTBase2d::remove_external_triangles(bool remove_internal_holes) { | |
| 826 | |||
| 827 |
2/2✓ Branch 0 taken 8067 times.
✓ Branch 1 taken 4 times.
|
8071 | if(remove_internal_holes) { |
| 828 | DList S(*this, DLIST_S_ID); | ||
| 829 | |||
| 830 | // Step 1: get triangles adjacent to the border, | ||
| 831 | // mark them as visited, classify them | ||
| 832 |
2/2✓ Branch 0 taken 103086 times.
✓ Branch 1 taken 8067 times.
|
222306 | for(index_t t=0; t<nT(); ++t) { |
| 833 |
2/2✓ Branch 0 taken 244810 times.
✓ Branch 1 taken 70818 times.
|
315628 | for(index_t le=0; le<3; ++le) { |
| 834 |
2/2✓ Branch 0 taken 32268 times.
✓ Branch 1 taken 212542 times.
|
244810 | if(Tadj(t,le) == NO_INDEX) { |
| 835 |
1/2✓ Branch 0 taken 32268 times.
✗ Branch 1 not taken.
|
32268 | bool outside = ((Tedge_cnstr_nb(t,le)%2) == 0); |
| 836 | Tset_flag(t, T_VISITED_FLAG); | ||
| 837 |
1/2✓ Branch 0 taken 32268 times.
✗ Branch 1 not taken.
|
32268 | if(outside) { |
| 838 | Tset_flag(t, T_MARKED_FLAG); | ||
| 839 | } | ||
| 840 | 32268 | S.push_back(t); | |
| 841 | break; | ||
| 842 | } | ||
| 843 | } | ||
| 844 | } | ||
| 845 | |||
| 846 | // Step 2: recursive traversal | ||
| 847 |
2/2✓ Branch 0 taken 103086 times.
✓ Branch 1 taken 8067 times.
|
111153 | while(!S.empty()) { |
| 848 | 103086 | index_t t1 = S.pop_back(); | |
| 849 | bool t1_outside = Tflag_is_set(t1, T_MARKED_FLAG); | ||
| 850 |
2/2✓ Branch 0 taken 309258 times.
✓ Branch 1 taken 103086 times.
|
412344 | for(index_t le=0; le<3; ++le) { |
| 851 | index_t t2 = Tadj(t1,le); | ||
| 852 | if( | ||
| 853 |
4/4✓ Branch 0 taken 276990 times.
✓ Branch 1 taken 32268 times.
✓ Branch 2 taken 70818 times.
✓ Branch 3 taken 206172 times.
|
309258 | t2 != NO_INDEX && |
| 854 | !Tflag_is_set(t2,T_VISITED_FLAG) | ||
| 855 | ) { | ||
| 856 | bool t2_outside = | ||
| 857 | t1_outside ^ ((Tedge_cnstr_nb(t1,le)%2) != 0); | ||
| 858 | Tset_flag(t2, T_VISITED_FLAG); | ||
| 859 |
2/2✓ Branch 0 taken 43272 times.
✓ Branch 1 taken 27546 times.
|
70818 | if(t2_outside) { |
| 860 | Tset_flag(t2, T_MARKED_FLAG); | ||
| 861 | } | ||
| 862 | 70818 | S.push_back(t2); | |
| 863 | } | ||
| 864 | } | ||
| 865 | } | ||
| 866 | |||
| 867 | // Step 3: reset visited flag | ||
| 868 |
2/2✓ Branch 0 taken 103086 times.
✓ Branch 1 taken 8067 times.
|
111153 | for(index_t t=0; t<nT(); ++t) { |
| 869 | Treset_flag(t, T_VISITED_FLAG); | ||
| 870 | } | ||
| 871 | |||
| 872 | 8067 | } else { | |
| 873 | DList S(*this, DLIST_S_ID); | ||
| 874 | |||
| 875 | // Step 1: get triangles adjacent to the border | ||
| 876 |
2/2✓ Branch 0 taken 40652 times.
✓ Branch 1 taken 4 times.
|
81312 | for(index_t t=0; t<nT(); ++t) { |
| 877 | if( | ||
| 878 | // TODO: replace with parity check ? | ||
| 879 |
3/4✓ Branch 0 taken 40396 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 40295 times.
✓ Branch 3 taken 357 times.
|
40652 | (!Tedge_is_constrained(t,0) && Tadj(t,0) == NO_INDEX) || |
| 880 |
5/6✓ Branch 0 taken 40396 times.
✓ Branch 1 taken 256 times.
✓ Branch 2 taken 40295 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 40445 times.
✓ Branch 5 taken 207 times.
|
81304 | (!Tedge_is_constrained(t,1) && Tadj(t,1) == NO_INDEX) || |
| 881 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 40445 times.
|
40445 | (!Tedge_is_constrained(t,2) && Tadj(t,2) == NO_INDEX) |
| 882 | ) { | ||
| 883 | Tset_flag(t, T_MARKED_FLAG); | ||
| 884 | ✗ | S.push_back(t); | |
| 885 | } | ||
| 886 | } | ||
| 887 | |||
| 888 | // Step 2: recursive traversal | ||
| 889 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | while(!S.empty()) { |
| 890 | ✗ | index_t t1 = S.pop_back(); | |
| 891 | ✗ | for(index_t le=0; le<3; ++le) { | |
| 892 | index_t t2 = Tadj(t1,le); | ||
| 893 | if( | ||
| 894 | ✗ | t2 != NO_INDEX && | |
| 895 | ✗ | !Tedge_is_constrained(t1,le) && | |
| 896 | !Tflag_is_set(t2,T_MARKED_FLAG) | ||
| 897 | ) { | ||
| 898 | Tset_flag(t2, T_MARKED_FLAG); | ||
| 899 | ✗ | S.push_back(t2); | |
| 900 | } | ||
| 901 | } | ||
| 902 | } | ||
| 903 | 4 | } | |
| 904 | |||
| 905 | // Step 3: remove marked triangles | ||
| 906 | 8071 | remove_marked_triangles(); | |
| 907 | 8071 | } | |
| 908 | |||
| 909 | |||
| 910 | 8177 | void CDTBase2d::remove_marked_triangles() { | |
| 911 | // Step 1: compute old2new map | ||
| 912 | // (use Tnext_'s storage, that we do not need now) | ||
| 913 | vector<index_t>& old2new = Tnext_; | ||
| 914 | index_t cur_t_new = 0; | ||
| 915 |
2/2✓ Branch 0 taken 150272 times.
✓ Branch 1 taken 8177 times.
|
316898 | for(index_t t=0; t<nT(); ++t) { |
| 916 |
2/2✓ Branch 0 taken 79335 times.
✓ Branch 1 taken 70937 times.
|
150272 | if(Tflag_is_set(t,T_MARKED_FLAG)) { |
| 917 | 79335 | old2new[t] = NO_INDEX; | |
| 918 | } else { | ||
| 919 | 70937 | old2new[t] = cur_t_new; | |
| 920 | 70937 | ++cur_t_new; | |
| 921 | } | ||
| 922 | } | ||
| 923 | index_t nT_new = cur_t_new; | ||
| 924 | |||
| 925 | // Step 2: translate adjacency and move triangles | ||
| 926 |
2/2✓ Branch 0 taken 150272 times.
✓ Branch 1 taken 8177 times.
|
308721 | for(index_t t=0; t<nT(); ++t) { |
| 927 | 150272 | index_t t_new = old2new[t]; | |
| 928 |
2/2✓ Branch 0 taken 79335 times.
✓ Branch 1 taken 70937 times.
|
150272 | if(t_new == NO_INDEX) { |
| 929 | 79335 | continue; | |
| 930 | } | ||
| 931 | index_t adj0 = Tadj(t,0); | ||
| 932 |
2/2✓ Branch 0 taken 70925 times.
✓ Branch 1 taken 12 times.
|
70937 | if(adj0 != NO_INDEX) { |
| 933 | 70925 | adj0 = old2new[adj0]; | |
| 934 | } | ||
| 935 | index_t adj1 = Tadj(t,1); | ||
| 936 |
1/2✓ Branch 0 taken 70937 times.
✗ Branch 1 not taken.
|
70937 | if(adj1 != NO_INDEX) { |
| 937 | 70937 | adj1 = old2new[adj1]; | |
| 938 | } | ||
| 939 | index_t adj2 = Tadj(t,2); | ||
| 940 |
1/2✓ Branch 0 taken 70937 times.
✗ Branch 1 not taken.
|
70937 | if(adj2 != NO_INDEX) { |
| 941 | 70937 | adj2 = old2new[adj2]; | |
| 942 | } | ||
| 943 | 70937 | Tset( | |
| 944 | t_new, | ||
| 945 | Tv(t,0), Tv(t,1), Tv(t,2), | ||
| 946 | adj0, adj1, adj2, | ||
| 947 | Tedge_cnstr_first(t,0), | ||
| 948 | Tedge_cnstr_first(t,1), | ||
| 949 | Tedge_cnstr_first(t,2) | ||
| 950 | ); | ||
| 951 | 70937 | Tflags_[t_new] = 0; | |
| 952 | } | ||
| 953 | |||
| 954 | // Step 3: resize arrays | ||
| 955 | 8177 | T_.resize(3*nT_new); | |
| 956 | 8177 | Tadj_.resize(3*nT_new); | |
| 957 | 8177 | Tflags_.resize(nT_new); | |
| 958 | 8177 | Tecnstr_first_.resize(3*nT_new); | |
| 959 | 8177 | Tnext_.resize(nT_new); | |
| 960 | 8177 | Tprev_.resize(nT_new); | |
| 961 | |||
| 962 | // Step 4: fix v2T_ | ||
| 963 |
2/2✓ Branch 0 taken 70937 times.
✓ Branch 1 taken 8177 times.
|
87291 | for(index_t t=0; t<nT(); ++t) { |
| 964 | 70937 | v2T_[Tv(t,0)] = t; | |
| 965 | 70937 | v2T_[Tv(t,1)] = t; | |
| 966 | 70937 | v2T_[Tv(t,2)] = t; | |
| 967 | } | ||
| 968 | 8177 | } | |
| 969 | |||
| 970 | |||
| 971 | /***************** Triangulation surgery (boring code ahead) *********/ | ||
| 972 | |||
| 973 |
2/2✓ Branch 0 taken 8127 times.
✓ Branch 1 taken 72448 times.
|
80575 | void CDTBase2d::insert_vertex_in_edge( |
| 974 | index_t v, index_t t, index_t le1, DList& S | ||
| 975 | ) { | ||
| 976 | index_t cnstr_first = Tedge_cnstr_first(t,le1); | ||
| 977 | index_t t1 = t; | ||
| 978 | index_t t2 = Tadj(t1,le1); | ||
| 979 | index_t v1 = Tv(t1,le1); | ||
| 980 | 80575 | index_t v2 = Tv(t1,(le1+1)%3); | |
| 981 | 80575 | index_t v3 = Tv(t1,(le1+2)%3); | |
| 982 | index_t t1_adj2 = Tadj(t1,(le1+1)%3); | ||
| 983 | index_t t1_adj3 = Tadj(t1,(le1+2)%3); | ||
| 984 |
2/2✓ Branch 0 taken 8127 times.
✓ Branch 1 taken 72448 times.
|
80575 | if(t2 != NO_INDEX) { |
| 985 | CDT_LOG(" insert vertex in internal edge"); | ||
| 986 | // New vertex is on an edge of t1 and t1 has a neighbor | ||
| 987 | // accross that edge. Discard the two triangles t1 and t2 | ||
| 988 | // adjacent to the edge, and create four new triangles | ||
| 989 | // (t1 and t2 are recycled). | ||
| 990 | index_t le2 = Tadj_find(t2,t1); | ||
| 991 | geo_debug_assert(Tv(t2, (le2+1)%3) == v3); | ||
| 992 | geo_debug_assert(Tv(t2, (le2+2)%3) == v2); | ||
| 993 | index_t v4 = Tv(t2,le2); | ||
| 994 | 8127 | index_t t2_adj2 = Tadj(t2,(le2+1)%3); | |
| 995 | 8127 | index_t t2_adj3 = Tadj(t2,(le2+2)%3); | |
| 996 | 8127 | index_t t3 = Tnew(); | |
| 997 | 8127 | index_t t4 = Tnew(); | |
| 998 | 8127 | Tset(t1,v,v1,v2,t1_adj3,t2,t4); | |
| 999 | 8127 | Tset(t2,v,v2,v4,t2_adj2,t3,t1); | |
| 1000 | 8127 | Tset(t3,v,v4,v3,t2_adj3,t4,t2); | |
| 1001 | 8127 | Tset(t4,v,v3,v1,t1_adj2,t1,t3); | |
| 1002 | 8127 | Tadj_back_connect(t1,0,t1); | |
| 1003 | 8127 | Tadj_back_connect(t2,0,t2); | |
| 1004 | 8127 | Tadj_back_connect(t3,0,t2); | |
| 1005 | 8127 | Tadj_back_connect(t4,0,t1); | |
| 1006 | Tset_edge_cnstr_first(t1,1,cnstr_first); | ||
| 1007 | Tset_edge_cnstr_first(t2,2,cnstr_first); | ||
| 1008 | Tset_edge_cnstr_first(t3,1,cnstr_first); | ||
| 1009 | Tset_edge_cnstr_first(t4,2,cnstr_first); | ||
| 1010 |
2/2✓ Branch 0 taken 2193 times.
✓ Branch 1 taken 5934 times.
|
8127 | if(S.initialized()) { |
| 1011 | 2193 | S.push_back(t1); | |
| 1012 | 2193 | S.push_back(t2); | |
| 1013 | 2193 | S.push_back(t3); | |
| 1014 | 2193 | S.push_back(t4); | |
| 1015 | } | ||
| 1016 | } else { | ||
| 1017 | CDT_LOG(" insert vertex in border edge"); | ||
| 1018 | // New vertex is on an edge of t1 and t1 has no neighbor | ||
| 1019 | // accross that edge. Discard t1 and replace it with two | ||
| 1020 | // new triangles (recycle t1). | ||
| 1021 | 72448 | t2 = Tnew(); | |
| 1022 | 72448 | Tset(t1,v,v1,v2,t1_adj3,NO_INDEX,t2); | |
| 1023 | 72448 | Tset(t2,v,v3,v1,t1_adj2,t1,NO_INDEX); | |
| 1024 | 72448 | Tadj_back_connect(t1,0,t1); | |
| 1025 | 72448 | Tadj_back_connect(t2,0,t1); | |
| 1026 | Tset_edge_cnstr_first(t1,1,cnstr_first); | ||
| 1027 | Tset_edge_cnstr_first(t2,2,cnstr_first); | ||
| 1028 |
2/2✓ Branch 0 taken 72257 times.
✓ Branch 1 taken 191 times.
|
72448 | if(S.initialized()) { |
| 1029 | 72257 | S.push_back(t1); | |
| 1030 | 72257 | S.push_back(t2); | |
| 1031 | } | ||
| 1032 | } | ||
| 1033 | 80575 | } | |
| 1034 | |||
| 1035 | 282916 | void CDTBase2d::insert_vertex_in_triangle(index_t v, index_t t, DList& S) { | |
| 1036 | // New vertex is in t1. Discard t1 and replace it with three | ||
| 1037 | // new triangles (recycle t1). | ||
| 1038 | index_t t1 = t; | ||
| 1039 | index_t v1 = Tv(t1,0); | ||
| 1040 | index_t v2 = Tv(t1,1); | ||
| 1041 | index_t v3 = Tv(t1,2); | ||
| 1042 | index_t adj1 = Tadj(t1,0); | ||
| 1043 | index_t adj2 = Tadj(t1,1); | ||
| 1044 | index_t adj3 = Tadj(t1,2); | ||
| 1045 | 282916 | index_t t2 = Tnew(); | |
| 1046 | 282916 | index_t t3 = Tnew(); | |
| 1047 | 282916 | Tset(t1,v,v2,v3,adj1,t2,t3); | |
| 1048 | 282916 | Tset(t2,v,v3,v1,adj2,t3,t1); | |
| 1049 | 282916 | Tset(t3,v,v1,v2,adj3,t1,t2); | |
| 1050 | 282916 | Tadj_back_connect(t1,0,t1); | |
| 1051 | 282916 | Tadj_back_connect(t2,0,t1); | |
| 1052 | 282916 | Tadj_back_connect(t3,0,t1); | |
| 1053 |
2/2✓ Branch 0 taken 171752 times.
✓ Branch 1 taken 111164 times.
|
282916 | if(S.initialized()) { |
| 1054 | 171752 | S.push_back(t1); | |
| 1055 | 171752 | S.push_back(t2); | |
| 1056 | 171752 | S.push_back(t3); | |
| 1057 | } | ||
| 1058 | 282916 | } | |
| 1059 | |||
| 1060 |
2/2✓ Branch 0 taken 334339 times.
✓ Branch 1 taken 238191 times.
|
572530 | void CDTBase2d::swap_edge(index_t t1, bool swap_t1_t2) { |
| 1061 | geo_debug_assert(!Tedge_is_constrained(t1,0)); | ||
| 1062 | index_t v1 = Tv(t1,0); | ||
| 1063 | index_t v2 = Tv(t1,1); | ||
| 1064 | index_t v3 = Tv(t1,2); | ||
| 1065 | index_t t1_adj2 = Tadj(t1,1); | ||
| 1066 | index_t t1_adj3 = Tadj(t1,2); | ||
| 1067 | index_t t2 = Tadj(t1,0); | ||
| 1068 | index_t le2 = Tadj_find(t2,t1); | ||
| 1069 | index_t v4 = Tv(t2,le2); | ||
| 1070 | geo_debug_assert(Tv(t2,(le2+1)%3) == v3); | ||
| 1071 | geo_debug_assert(Tv(t2,(le2+2)%3) == v2); | ||
| 1072 | |||
| 1073 | debug_Tcheck(t1); | ||
| 1074 | debug_Tcheck(t2); | ||
| 1075 | |||
| 1076 |
2/2✓ Branch 0 taken 21830 times.
✓ Branch 1 taken 550700 times.
|
572530 | index_t t2_adj2 = Tadj(t2,(le2+1)%3); |
| 1077 | 572530 | index_t t2_adj3 = Tadj(t2,(le2+2)%3); | |
| 1078 |
2/2✓ Branch 0 taken 21830 times.
✓ Branch 1 taken 550700 times.
|
572530 | if(swap_t1_t2) { |
| 1079 | 21830 | Tset(t2,v1,v4,v3,t2_adj3,t1_adj2,t1); | |
| 1080 | 21830 | Tset(t1,v1,v2,v4,t2_adj2,t2,t1_adj3); | |
| 1081 | 21830 | Tadj_back_connect(t2,0,t2); | |
| 1082 | 21830 | Tadj_back_connect(t2,1,t1); | |
| 1083 | 21830 | Tadj_back_connect(t1,0,t2); | |
| 1084 | 21830 | Tadj_back_connect(t1,2,t1); | |
| 1085 | } else { | ||
| 1086 | 550700 | Tset(t1,v1,v4,v3,t2_adj3,t1_adj2,t2); | |
| 1087 | 550700 | Tset(t2,v1,v2,v4,t2_adj2,t1,t1_adj3); | |
| 1088 | 550700 | Tadj_back_connect(t1,0,t2); | |
| 1089 | 550700 | Tadj_back_connect(t1,1,t1); | |
| 1090 | 550700 | Tadj_back_connect(t2,0,t2); | |
| 1091 | 550700 | Tadj_back_connect(t2,2,t1); | |
| 1092 | } | ||
| 1093 | |||
| 1094 | debug_Tcheck(t1); | ||
| 1095 | debug_Tcheck(t2); | ||
| 1096 | 572530 | } | |
| 1097 | |||
| 1098 | /***************** Geometry ***********************/ | ||
| 1099 | |||
| 1100 |
2/2✓ Branch 0 taken 942853 times.
✓ Branch 1 taken 650431 times.
|
1593284 | bool CDTBase2d::is_convex_quad(index_t t) const { |
| 1101 | index_t v1 = Tv(t,0); | ||
| 1102 | index_t v2 = Tv(t,1); | ||
| 1103 | index_t v3 = Tv(t,2); | ||
| 1104 | index_t t2 = Tadj(t,0); | ||
| 1105 | index_t le2 = Tadj_find(t2,t); | ||
| 1106 | index_t v4 = Tv(t2,le2); | ||
| 1107 | // t and Tadj(t,0) have the correct orientation, | ||
| 1108 | // so one just needs to check the orientation of | ||
| 1109 | // the two triangles that would be generated by | ||
| 1110 | // an edge flip. | ||
| 1111 | return | ||
| 1112 |
2/2✓ Branch 1 taken 1464800 times.
✓ Branch 2 taken 128484 times.
|
1593284 | orient2d(v1,v4,v3) == orient_012_ && |
| 1113 |
2/2✓ Branch 1 taken 136512 times.
✓ Branch 2 taken 1328288 times.
|
1464800 | orient2d(v4,v1,v2) == orient_012_ ; |
| 1114 | } | ||
| 1115 | |||
| 1116 | /******* Debugging ******************************************************/ | ||
| 1117 | |||
| 1118 | 64 | void CDTBase2d::check_geometry() const { | |
| 1119 |
3/4✓ Branch 0 taken 32 times.
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 32 times.
|
64 | if(delaunay_ && exact_incircle_) { |
| 1120 | ✗ | for(index_t t=0; t<nT(); ++t) { | |
| 1121 | for(index_t le=0; le<3; ++le) { | ||
| 1122 | // geo_assert(Tedge_is_Delaunay(t,le)); | ||
| 1123 | } | ||
| 1124 | } | ||
| 1125 | } | ||
| 1126 | 64 | } | |
| 1127 | |||
| 1128 | |||
| 1129 | ✗ | bool CDTBase2d::Tedge_is_Delaunay(index_t t1, index_t le1) const { | |
| 1130 | ✗ | if(Tedge_is_constrained(t1,le1)) { | |
| 1131 | return true; | ||
| 1132 | } | ||
| 1133 | index_t t2 = Tadj(t1,le1); | ||
| 1134 | ✗ | if(t2 == NO_INDEX) { | |
| 1135 | return true; | ||
| 1136 | } | ||
| 1137 | index_t le2 = Tadj_find(t2,t1); | ||
| 1138 | index_t v1 = Tv(t1,le1); | ||
| 1139 | ✗ | index_t v2 = Tv(t1,(le1+1)%3); | |
| 1140 | ✗ | index_t v3 = Tv(t1,(le1+2)%3); | |
| 1141 | index_t v4 = Tv(t2,le2); | ||
| 1142 | |||
| 1143 | // If we do not check that, we assert fail | ||
| 1144 | // whenever there is a vertex inserted in | ||
| 1145 | // the macroborder of the triangle | ||
| 1146 | if( | ||
| 1147 | ✗ | orient2d(v1,v4,v3) != orient_012_ || | |
| 1148 | ✗ | orient2d(v4,v1,v2) != orient_012_ | |
| 1149 | ) { | ||
| 1150 | ✗ | return true; | |
| 1151 | } | ||
| 1152 | |||
| 1153 | ✗ | return Sign(incircle(v1,v2,v3,v4)*orient_012_) <= 0; | |
| 1154 | } | ||
| 1155 | |||
| 1156 | ✗ | void CDTBase2d::check_edge_intersections( | |
| 1157 | index_t v1, index_t v2, const DList& Q | ||
| 1158 | ) { | ||
| 1159 | std::set<Edge> I; | ||
| 1160 | auto make_edge = [](index_t w1, index_t w2)->Edge { | ||
| 1161 | return std::make_pair(std::min(w1,w2), std::max(w1,w2)); | ||
| 1162 | }; | ||
| 1163 | ✗ | for(index_t t=Q.front(); t!=NO_INDEX; t = Q.next(t)) { | |
| 1164 | ✗ | geo_assert(segment_edge_intersect(v1,v2,t,0)); | |
| 1165 | ✗ | I.insert(make_edge(Tv(t,1), Tv(t,2))); | |
| 1166 | } | ||
| 1167 | ✗ | for(index_t t=0; t<nT(); ++t) { | |
| 1168 | ✗ | for(index_t le=0; le<3; ++le) { | |
| 1169 | ✗ | if(segment_edge_intersect(v1,v2,t,le)) { | |
| 1170 | ✗ | index_t w1 = Tv(t,(le+1)%3); | |
| 1171 | ✗ | index_t w2 = Tv(t,(le+2)%3); | |
| 1172 | ✗ | geo_assert(I.find(make_edge(w1,w2)) != I.end()); | |
| 1173 | } | ||
| 1174 | } | ||
| 1175 | } | ||
| 1176 | ✗ | } | |
| 1177 | |||
| 1178 | /*** Naive versions of algorithms, for reference / debugging if need be ***/ | ||
| 1179 | |||
| 1180 | ✗ | index_t CDTBase2d::locate_naive(index_t v, index_t hint, Sign* o) const { | |
| 1181 | geo_argused(hint); | ||
| 1182 | Sign o_local[3]; | ||
| 1183 | ✗ | if(o == nullptr) { | |
| 1184 | o = o_local; | ||
| 1185 | } | ||
| 1186 | |||
| 1187 | ✗ | for(index_t t=0; t<nT(); ++t) { | |
| 1188 | index_t i = Tv(t,0); | ||
| 1189 | index_t j = Tv(t,1); | ||
| 1190 | index_t k = Tv(t,2); | ||
| 1191 | ✗ | o[0] = orient2d(v,j,k); | |
| 1192 | ✗ | o[1] = orient2d(v,k,i); | |
| 1193 | ✗ | o[2] = orient2d(v,i,j); | |
| 1194 | ✗ | if(o[0]*o[1] >= 0 && o[1]*o[2] >= 0 && o[2]*o[0] >= 0) { | |
| 1195 | ✗ | return t; | |
| 1196 | } | ||
| 1197 | } | ||
| 1198 | ✗ | geo_assert_not_reached; | |
| 1199 | } | ||
| 1200 | |||
| 1201 | ✗ | void CDTBase2d::Delaunayize_new_edges_naive(vector<Edge>& N) { | |
| 1202 | CDT_LOG("Delaunayize_new_edges_naive()"); | ||
| 1203 | ✗ | for(Edge E: N) { | |
| 1204 | index_t v1 = std::min(E.first,E.second); | ||
| 1205 | index_t v2 = std::max(E.first,E.second); | ||
| 1206 | if(v2 < v1) { | ||
| 1207 | std::swap(v1,v2); | ||
| 1208 | } | ||
| 1209 | CDT_LOG("new edge: " << v1 << " " << v2); | ||
| 1210 | } | ||
| 1211 | index_t count = 0; | ||
| 1212 | bool swap_occured = true; | ||
| 1213 | ✗ | while(swap_occured) { | |
| 1214 | // NASA programming style: all loops have | ||
| 1215 | // a maximum number of iterations | ||
| 1216 | ✗ | ++count; | |
| 1217 | ✗ | if(count > 10*nT()) { | |
| 1218 | ✗ | Logger::warn("CDT2d") | |
| 1219 | << "Emergency exit in Delaunayize_new_edges_naive()" | ||
| 1220 | << std::endl; | ||
| 1221 | ✗ | return; | |
| 1222 | } | ||
| 1223 | swap_occured = false; | ||
| 1224 | ✗ | for(Edge& E: N) { | |
| 1225 | ✗ | index_t t1 = eT(E); | |
| 1226 | ✗ | if(Tedge_is_constrained(t1,0)) { | |
| 1227 | ✗ | continue; | |
| 1228 | } | ||
| 1229 | index_t v1 = Tv(t1,1); | ||
| 1230 | index_t v2 = Tv(t1,2); | ||
| 1231 | index_t v0 = Tv(t1,0); | ||
| 1232 | index_t t2 = Tadj(t1,0); | ||
| 1233 | ✗ | if(t2 == NO_INDEX) { | |
| 1234 | ✗ | continue; | |
| 1235 | } | ||
| 1236 | index_t e2 = Tadj_find(t2,t1); | ||
| 1237 | index_t v3 = Tv(t2,e2); | ||
| 1238 | ✗ | if(!exact_incircle_ && !is_convex_quad(t1)) { | |
| 1239 | ✗ | continue; | |
| 1240 | } | ||
| 1241 | ✗ | if(Sign(incircle(v0,v1,v2,v3)*orient_012_) == POSITIVE) { | |
| 1242 | CDT_LOG("swap " << v1 << " " << v2 | ||
| 1243 | << " ---> " | ||
| 1244 | << v0 << " " << v3 ); | ||
| 1245 | ✗ | swap_edge(t1); | |
| 1246 | E = std::make_pair(Tv(t1,0), Tv(t1,1)); | ||
| 1247 | swap_occured = true; | ||
| 1248 | } | ||
| 1249 | } | ||
| 1250 | } | ||
| 1251 | ✗ | N.resize(0); | |
| 1252 | CDT_LOG("/Delaunayize_new_edges_naive()"); | ||
| 1253 | } | ||
| 1254 | |||
| 1255 | ✗ | void CDTBase2d::constrain_edges_naive( | |
| 1256 | index_t i, index_t j, DList& Q_in, vector<Edge>& N | ||
| 1257 | ) { | ||
| 1258 | CDT_LOG("Q size=" << Q_in.size()); | ||
| 1259 | |||
| 1260 | std::deque<Edge> Q; | ||
| 1261 | ✗ | for(index_t t=Q_in.front(); t != NO_INDEX; t = Q_in.next(t)) { | |
| 1262 | ✗ | Q.push_back(std::make_pair(Tv(t,1), Tv(t,2))); | |
| 1263 | } | ||
| 1264 | Q_in.clear(); | ||
| 1265 | |||
| 1266 | ✗ | for(index_t t=0; t<nT(); ++t) { | |
| 1267 | geo_debug_assert(!Tis_in_list(t)); | ||
| 1268 | } | ||
| 1269 | |||
| 1270 | ✗ | while(Q.size() != 0) { | |
| 1271 | ✗ | Edge E = Q.back(); | |
| 1272 | Q.pop_back(); | ||
| 1273 | ✗ | if(!is_convex_quad(eT(E))) { | |
| 1274 | ✗ | if(Q.size() == 0) { | |
| 1275 | CDT_LOG("... infinite iteration"); | ||
| 1276 | ✗ | abort(); | |
| 1277 | } | ||
| 1278 | Q.push_front(E); | ||
| 1279 | } else { | ||
| 1280 | ✗ | index_t t = eT(E); | |
| 1281 | ✗ | swap_edge(t); | |
| 1282 | E = std::make_pair(Tv(t,0), Tv(t,1)); | ||
| 1283 | ✗ | if(segment_segment_intersect(i,j,E.first,E.second)) { | |
| 1284 | Q.push_front(E); | ||
| 1285 | } else { | ||
| 1286 | ✗ | if( | |
| 1287 | ✗ | (E.first == i && E.second == j) || | |
| 1288 | ✗ | (E.first == j && E.second == i) | |
| 1289 | ) { | ||
| 1290 | ✗ | Tadd_edge_cnstr_with_neighbor(eT(E),0,ncnstr_-1); | |
| 1291 | } else { | ||
| 1292 | ✗ | N.push_back(E); | |
| 1293 | } | ||
| 1294 | } | ||
| 1295 | } | ||
| 1296 | } | ||
| 1297 | ✗ | } | |
| 1298 | |||
| 1299 | /********************************************************************/ | ||
| 1300 | |||
| 1301 | 32 | CDT2d::CDT2d() { | |
| 1302 | 32 | exact_intersections_ = false; | |
| 1303 | 32 | exact_incircle_ = false; | |
| 1304 | 32 | } | |
| 1305 | |||
| 1306 |
1/2✓ Branch 0 taken 32 times.
✗ Branch 1 not taken.
|
64 | CDT2d::~CDT2d() { |
| 1307 | 64 | } | |
| 1308 | |||
| 1309 | ✗ | void CDT2d::clear() { | |
| 1310 | ✗ | CDTBase2d::clear(); | |
| 1311 | ✗ | point_.resize(0); | |
| 1312 | ✗ | } | |
| 1313 | |||
| 1314 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
|
32 | void CDT2d::create_enclosing_triangle( |
| 1315 | const vec2& p1, const vec2& p2, const vec2& p3 | ||
| 1316 | ) { | ||
| 1317 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
32 | geo_assert(nv() == 0); |
| 1318 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
32 | geo_assert(nT() == 0); |
| 1319 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
|
32 | point_.push_back(p1); |
| 1320 | point_.push_back(p2); | ||
| 1321 | point_.push_back(p3); | ||
| 1322 | 32 | CDTBase2d::create_enclosing_triangle(0,1,2); | |
| 1323 | 32 | } | |
| 1324 | |||
| 1325 | ✗ | void CDT2d::create_enclosing_quad( | |
| 1326 | const vec2& p1, const vec2& p2, const vec2& p3, const vec2& p4 | ||
| 1327 | ) { | ||
| 1328 | ✗ | geo_assert(nv() == 0); | |
| 1329 | ✗ | geo_assert(nT() == 0); | |
| 1330 | ✗ | point_.push_back(p1); | |
| 1331 | point_.push_back(p2); | ||
| 1332 | point_.push_back(p3); | ||
| 1333 | point_.push_back(p4); | ||
| 1334 | ✗ | CDTBase2d::create_enclosing_quad(0,1,2,3); | |
| 1335 | ✗ | } | |
| 1336 | |||
| 1337 | 7300063 | Sign CDT2d::orient2d(index_t i, index_t j, index_t k) const { | |
| 1338 | geo_debug_assert(i < nv()); | ||
| 1339 | geo_debug_assert(j < nv()); | ||
| 1340 | geo_debug_assert(k < nv()); | ||
| 1341 | 7300063 | return PCK::orient_2d( | |
| 1342 | point_[i].data(), point_[j].data(), point_[k].data() | ||
| 1343 | 7300063 | ); | |
| 1344 | } | ||
| 1345 | |||
| 1346 | 935458 | Sign CDT2d::incircle(index_t i, index_t j, index_t k, index_t l) const { | |
| 1347 | geo_debug_assert(i < nv()); | ||
| 1348 | geo_debug_assert(j < nv()); | ||
| 1349 | geo_debug_assert(k < nv()); | ||
| 1350 | geo_debug_assert(l < nv()); | ||
| 1351 | 935458 | return PCK::in_circle_2d_SOS( | |
| 1352 | point_[i].data(), point_[j].data(), point_[k].data(), | ||
| 1353 | point_[l].data() | ||
| 1354 | 935458 | ); | |
| 1355 | } | ||
| 1356 | |||
| 1357 |
2/2✓ Branch 0 taken 152 times.
✓ Branch 1 taken 16 times.
|
168 | index_t CDT2d::create_intersection( |
| 1358 | index_t E1, index_t i, index_t j, | ||
| 1359 | index_t E2, index_t k, index_t l | ||
| 1360 | ) { | ||
| 1361 | geo_argused(E1); | ||
| 1362 | geo_argused(E2); | ||
| 1363 | geo_debug_assert(i < nv()); | ||
| 1364 | geo_debug_assert(j < nv()); | ||
| 1365 | geo_debug_assert(k < nv()); | ||
| 1366 | geo_debug_assert(l < nv()); | ||
| 1367 | geo_debug_assert(E1 < ncnstr()); | ||
| 1368 | geo_debug_assert(E2 < ncnstr()); | ||
| 1369 | vec2 U = point_[j] - point_[i]; | ||
| 1370 | vec2 V = point_[l] - point_[k]; | ||
| 1371 | vec2 D = point_[k] - point_[i]; | ||
| 1372 | double delta = det(U,V); | ||
| 1373 |
2/2✓ Branch 0 taken 152 times.
✓ Branch 1 taken 16 times.
|
168 | double t = det(D,V)/delta; |
| 1374 | 168 | vec2 P = point_[i] + t*U; | |
| 1375 |
2/2✓ Branch 0 taken 152 times.
✓ Branch 1 taken 16 times.
|
168 | point_.push_back(P); |
| 1376 |
2/2✓ Branch 0 taken 152 times.
✓ Branch 1 taken 16 times.
|
168 | v2T_.push_back(NO_INDEX); |
| 1377 | 168 | index_t v = nv_; | |
| 1378 | 168 | ++nv_; | |
| 1379 | 168 | return v; | |
| 1380 | } | ||
| 1381 | |||
| 1382 |
1/2✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
|
32 | void CDT2d::insert( |
| 1383 | index_t nb_points, const double* points, | ||
| 1384 | index_t* indices, bool remove_unreferenced_vertices | ||
| 1385 | ) { | ||
| 1386 | CDT_LOG("Inserting " << nb_points << " points"); | ||
| 1387 | debug_check_consistency(); | ||
| 1388 | |||
| 1389 | // Compute spatial sort | ||
| 1390 | vector<index_t> sorted_indices; | ||
| 1391 |
1/2✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
|
32 | compute_BRIO_order(nb_points, points, sorted_indices, 2, 2); |
| 1392 | |||
| 1393 | // Insert vertices one by one, following the order given | ||
| 1394 | // by spatial sort. | ||
| 1395 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
|
32 | if(remove_unreferenced_vertices) { |
| 1396 | |||
| 1397 | // Pre-allocate memory | ||
| 1398 | ✗ | point_.reserve(point_.size()+nb_points); | |
| 1399 | ✗ | v2T_.reserve(v2T_.size()+nb_points); | |
| 1400 | |||
| 1401 | // Insert the points and vertices one by one, following | ||
| 1402 | // spatial sort order. | ||
| 1403 | index_t hint = NO_INDEX; | ||
| 1404 | ✗ | for(index_t i=0; i<nb_points; ++i) { | |
| 1405 | ✗ | point_.push_back(vec2(points+2*sorted_indices[i])); | |
| 1406 | ✗ | index_t v = CDTBase2d::insert(point_.size()-1, hint); | |
| 1407 | |||
| 1408 | // If it was a duplicated point, then remove the point | ||
| 1409 | ✗ | if(point_.size() > nv()) { | |
| 1410 | point_.pop_back(); | ||
| 1411 | } | ||
| 1412 | |||
| 1413 | ✗ | indices[sorted_indices[i]] = v; | |
| 1414 | hint = vT(v); | ||
| 1415 | } | ||
| 1416 | |||
| 1417 | } else { | ||
| 1418 | |||
| 1419 | // Insert all the points in the point_ vector | ||
| 1420 | index_t v_offset = nv(); | ||
| 1421 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
✓ Branch 3 taken 32 times.
✗ Branch 4 not taken.
|
32 | point_.reserve(point_.size()+nb_points); |
| 1422 |
2/2✓ Branch 0 taken 222718 times.
✓ Branch 1 taken 32 times.
|
222750 | for(index_t i=0; i<nb_points; ++i) { |
| 1423 |
1/4✓ Branch 1 taken 222718 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
222718 | point_.push_back(vec2(points+2*i)); |
| 1424 | } | ||
| 1425 | |||
| 1426 | // Resize vertex-to-triangle array accordingly, | ||
| 1427 | // and update number of points | ||
| 1428 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 32 times.
✓ Branch 3 taken 32 times.
✗ Branch 4 not taken.
|
32 | v2T_.resize(v2T_.size()+nb_points, NO_INDEX); |
| 1429 | 32 | nv_+=nb_points; | |
| 1430 | |||
| 1431 | // Now insert the vertices in the triangulation, | ||
| 1432 | // following the order of spatial search (but | ||
| 1433 | // this will not change the order of the points, | ||
| 1434 | // in contrast with the "remove_unreferenced_vertices" | ||
| 1435 | // alternative). In the end, each duplicated point | ||
| 1436 | // v has vT(v) == NO_INDEX (no incident triangle). | ||
| 1437 | index_t hint = NO_INDEX; | ||
| 1438 |
2/2✓ Branch 0 taken 222718 times.
✓ Branch 1 taken 32 times.
|
222750 | for(index_t i=0; i<nb_points; ++i) { |
| 1439 |
1/2✓ Branch 1 taken 222718 times.
✗ Branch 2 not taken.
|
222718 | index_t v = CDTBase2d::insert( |
| 1440 | v_offset+sorted_indices[i], hint | ||
| 1441 | ); | ||
| 1442 |
1/2✓ Branch 0 taken 222718 times.
✗ Branch 1 not taken.
|
222718 | if(indices != nullptr) { |
| 1443 | 222718 | indices[sorted_indices[i]] = v; | |
| 1444 | } | ||
| 1445 | hint = vT(v); | ||
| 1446 | } | ||
| 1447 | } | ||
| 1448 | CDT_LOG("Inserted."); | ||
| 1449 | debug_check_consistency(); | ||
| 1450 | 32 | } | |
| 1451 | |||
| 1452 | ✗ | void CDT2d::save(const std::string& filename) const { | |
| 1453 | #ifndef GEOGRAM_PSM | ||
| 1454 | ✗ | Mesh M; | |
| 1455 | ✗ | M.vertices.set_dimension(2); | |
| 1456 | ✗ | for(const vec2& P: point_) { | |
| 1457 | ✗ | M.vertices.create_vertex(P.data()); | |
| 1458 | } | ||
| 1459 | ✗ | for(index_t t=0; t<nT(); ++t) { | |
| 1460 | index_t i = Tv(t,0); | ||
| 1461 | index_t j = Tv(t,1); | ||
| 1462 | index_t k = Tv(t,2); | ||
| 1463 | ✗ | M.facets.create_triangle(i,j,k); | |
| 1464 | } | ||
| 1465 | |||
| 1466 | |||
| 1467 | Attribute<double> tex_coord; | ||
| 1468 | ✗ | tex_coord.create_vector_attribute( | |
| 1469 | ✗ | M.facet_corners.attributes(), "tex_coord", 2 | |
| 1470 | ); | ||
| 1471 | static double triangle_tex[3][2] = { | ||
| 1472 | {0.0, 0.0}, | ||
| 1473 | {1.0, 0.0}, | ||
| 1474 | {0.0, 1.0} | ||
| 1475 | }; | ||
| 1476 | ✗ | for(index_t c: M.facet_corners) { | |
| 1477 | ✗ | tex_coord[2*c] = triangle_tex[c%3][0]; | |
| 1478 | ✗ | tex_coord[2*c+1] = triangle_tex[c%3][1]; | |
| 1479 | } | ||
| 1480 | |||
| 1481 | ✗ | Attribute<bool> constraint(M.facet_corners.attributes(), "constraint"); | |
| 1482 | ✗ | for(index_t c: M.facet_corners) { | |
| 1483 | ✗ | index_t t = c/3; | |
| 1484 | ✗ | index_t lv = c%3; | |
| 1485 | constraint[c] = | ||
| 1486 | ✗ | Tedge_is_constrained(t, (lv+1)%3) || | |
| 1487 | ✗ | Tedge_is_constrained(t, (lv+2)%3) ; | |
| 1488 | } | ||
| 1489 | |||
| 1490 | ✗ | for(index_t t=0; t<nT(); ++t) { | |
| 1491 | ✗ | for(index_t le=0; le<3; ++le) { | |
| 1492 | ✗ | if(Tedge_is_constrained(t,le)) { | |
| 1493 | ✗ | index_t v1 = Tv(t, (le+1)%3); | |
| 1494 | ✗ | index_t v2 = Tv(t, (le+2)%3); | |
| 1495 | ✗ | M.edges.create_edge(v1,v2); | |
| 1496 | } | ||
| 1497 | } | ||
| 1498 | } | ||
| 1499 | |||
| 1500 | |||
| 1501 | ✗ | M.facets.connect(); | |
| 1502 | |||
| 1503 | ✗ | mesh_save(M, filename); | |
| 1504 | #else | ||
| 1505 | if(!String::string_ends_with(filename,".obj")) { | ||
| 1506 | Logger::err("CDT_2d") | ||
| 1507 | << "save() only supports .obj file format in PSM" | ||
| 1508 | << std::endl; | ||
| 1509 | return; | ||
| 1510 | } | ||
| 1511 | std::ofstream out(filename); | ||
| 1512 | for(index_t v=0; v<nv(); ++v) { | ||
| 1513 | out << "v " << point(v) << " " << 0.0 << std::endl; | ||
| 1514 | } | ||
| 1515 | for(index_t t=0; t<nT(); ++t) { | ||
| 1516 | out << "f " << Tv(t,0)+1 << " " << Tv(t,1)+1 << " " << Tv(t,2)+1 | ||
| 1517 | << std::endl; | ||
| 1518 | } | ||
| 1519 | |||
| 1520 | for(index_t t=0; t<nT(); ++t) { | ||
| 1521 | for(index_t le=0; le<3; ++le) { | ||
| 1522 | if(Tedge_is_constrained(t,le)) { | ||
| 1523 | index_t v1 = Tv(t,(le+1)%3); | ||
| 1524 | index_t v2 = Tv(t,(le+2)%3); | ||
| 1525 | out << "l " << v1+1 << " " << v2+1 << std::endl; | ||
| 1526 | } | ||
| 1527 | } | ||
| 1528 | } | ||
| 1529 | #endif | ||
| 1530 | ✗ | } | |
| 1531 | |||
| 1532 | /**************************************************************************/ | ||
| 1533 | |||
| 1534 | 391 | ExactCDT2d::ExactCDT2d(): | |
| 1535 | 391 | use_pred_cache_insert_buffer_(false) { | |
| 1536 | #ifdef GEOGRAM_USE_EXACT_NT | ||
| 1537 | CDTBase2d::exact_incircle_ = true; | ||
| 1538 | #else | ||
| 1539 | // Since incircle() with expansions computes approximated | ||
| 1540 | // lifted coordinate, we need to activate additional | ||
| 1541 | // checks for Delaunayization. | ||
| 1542 | 391 | CDTBase2d::exact_incircle_ = false; | |
| 1543 | #endif | ||
| 1544 | 391 | } | |
| 1545 | |||
| 1546 |
2/2✓ Branch 0 taken 323 times.
✓ Branch 1 taken 68 times.
|
782 | ExactCDT2d::~ExactCDT2d() { |
| 1547 | 782 | } | |
| 1548 | |||
| 1549 | 8067 | void ExactCDT2d::clear() { | |
| 1550 | 8067 | point_.resize(0); | |
| 1551 | 8067 | id_.resize(0); | |
| 1552 | pred_cache_.clear(); | ||
| 1553 | 8067 | pred_cache_insert_buffer_.resize(0); | |
| 1554 | 8067 | use_pred_cache_insert_buffer_ = false; | |
| 1555 | 8067 | cnstr_operand_bits_.resize(0); | |
| 1556 | 8067 | constraints_.resize(0); | |
| 1557 | 8067 | CDTBase2d::clear(); | |
| 1558 | #ifndef GEOGRAM_USE_EXACT_NT | ||
| 1559 | 8067 | length_.resize(0); | |
| 1560 | #endif | ||
| 1561 | 8067 | } | |
| 1562 | |||
| 1563 | ✗ | void ExactCDT2d::create_enclosing_triangle( | |
| 1564 | const ExactPoint& p1, const ExactPoint& p2, const ExactPoint& p3 | ||
| 1565 | ) { | ||
| 1566 | ✗ | geo_assert(nv() == 0); | |
| 1567 | ✗ | geo_assert(nT() == 0); | |
| 1568 | #ifndef GEOGRAM_USE_EXACT_NT | ||
| 1569 | geo_debug_assert(length_.size() == 0); | ||
| 1570 | #endif | ||
| 1571 | ✗ | add_point(p1); | |
| 1572 | ✗ | add_point(p2); | |
| 1573 | ✗ | add_point(p3); | |
| 1574 | ✗ | CDTBase2d::create_enclosing_triangle(0,1,2); | |
| 1575 | ✗ | } | |
| 1576 | |||
| 1577 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8173 times.
|
8173 | void ExactCDT2d::create_enclosing_quad( |
| 1578 | const ExactPoint& p1, const ExactPoint& p2, | ||
| 1579 | const ExactPoint& p3, const ExactPoint& p4 | ||
| 1580 | ) { | ||
| 1581 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 8173 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
8173 | geo_assert(nv() == 0); |
| 1582 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 8173 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
8173 | geo_assert(nT() == 0); |
| 1583 | #ifndef GEOGRAM_USE_EXACT_NT | ||
| 1584 | geo_debug_assert(length_.size() == 0); | ||
| 1585 | #endif | ||
| 1586 | 8173 | add_point(p1); | |
| 1587 | 8173 | add_point(p2); | |
| 1588 | 8173 | add_point(p3); | |
| 1589 | 8173 | add_point(p4); | |
| 1590 | 8173 | CDTBase2d::create_enclosing_quad(0,1,2,3); | |
| 1591 | 8173 | } | |
| 1592 | |||
| 1593 | 46603 | index_t ExactCDT2d::insert(const ExactPoint& p, index_t id, index_t hint) { | |
| 1594 | #ifndef GEOGRAM_USE_EXACT_NT | ||
| 1595 | geo_debug_assert(nv() == length_.size()); | ||
| 1596 | #endif | ||
| 1597 | debug_check_consistency(); | ||
| 1598 | 46603 | add_point(p,id); | |
| 1599 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 46603 times.
|
46603 | index_t v = CDTBase2d::insert(point_.size()-1, hint); |
| 1600 | // If inserted point already existed in | ||
| 1601 | // triangulation, then nv() did not increase | ||
| 1602 |
2/2✓ Branch 0 taken 47 times.
✓ Branch 1 taken 46556 times.
|
46603 | if(point_.size() > nv()) { |
| 1603 | point_.pop_back(); | ||
| 1604 | id_.pop_back(); | ||
| 1605 | #ifndef GEOGRAM_USE_EXACT_NT | ||
| 1606 | length_.pop_back(); | ||
| 1607 | #endif | ||
| 1608 | } | ||
| 1609 | debug_check_consistency(); | ||
| 1610 | #ifndef GEOGRAM_USE_EXACT_NT | ||
| 1611 | geo_debug_assert(nv() == length_.size()); | ||
| 1612 | #endif | ||
| 1613 | 46603 | return v; | |
| 1614 | } | ||
| 1615 | |||
| 1616 | 79295 | void ExactCDT2d::add_point(const ExactPoint& p, index_t id) { | |
| 1617 | 79295 | point_.push_back(p); | |
| 1618 |
2/2✓ Branch 0 taken 77593 times.
✓ Branch 1 taken 1702 times.
|
79295 | id_.push_back(id); |
| 1619 | #ifndef GEOGRAM_USE_EXACT_NT | ||
| 1620 | 79295 | length_.push_back( | |
| 1621 |
3/8✓ Branch 3 taken 79295 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 79295 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 79295 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
396475 | (geo_sqr(p.x) + geo_sqr(p.y)).estimate() / |
| 1622 |
1/4✓ Branch 1 taken 79295 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
158590 | geo_sqr(p.w).estimate() |
| 1623 | ); | ||
| 1624 | #endif | ||
| 1625 | 79295 | } | |
| 1626 | |||
| 1627 | 46603 | void ExactCDT2d::begin_insert_transaction() { | |
| 1628 | 46603 | use_pred_cache_insert_buffer_ = true; | |
| 1629 | 46603 | } | |
| 1630 | |||
| 1631 | 46556 | void ExactCDT2d::commit_insert_transaction() { | |
| 1632 |
2/2✓ Branch 0 taken 321416 times.
✓ Branch 1 taken 46556 times.
|
367972 | for(const auto& it: pred_cache_insert_buffer_) { |
| 1633 | 321416 | pred_cache_[it.first] = it.second; | |
| 1634 | } | ||
| 1635 | 46556 | pred_cache_insert_buffer_.resize(0); | |
| 1636 | 46556 | use_pred_cache_insert_buffer_ = false; | |
| 1637 | 46556 | } | |
| 1638 | |||
| 1639 | 47 | void ExactCDT2d::rollback_insert_transaction() { | |
| 1640 | 47 | pred_cache_insert_buffer_.resize(0); | |
| 1641 | 47 | use_pred_cache_insert_buffer_ = false; | |
| 1642 | 47 | } | |
| 1643 | |||
| 1644 | /** | ||
| 1645 | * \brief Tests the parity of the permutation of a list of | ||
| 1646 | * three distinct indices with respect to the canonical order. | ||
| 1647 | */ | ||
| 1648 | 903208 | static bool odd_order(index_t i, index_t j, index_t k) { | |
| 1649 | // Implementation: sort the elements (bubble sort is OK for | ||
| 1650 | // such a small number), and invert parity each time | ||
| 1651 | // two elements are swapped. | ||
| 1652 | 903208 | index_t tab[3] = { i, j, k}; | |
| 1653 | const int N = 3; | ||
| 1654 | bool result = false; | ||
| 1655 |
2/2✓ Branch 0 taken 1806416 times.
✓ Branch 1 taken 903208 times.
|
2709624 | for (int I = 0; I < N - 1; ++I) { |
| 1656 |
2/2✓ Branch 0 taken 2709624 times.
✓ Branch 1 taken 1806416 times.
|
4516040 | for (int J = 0; J < N - I - 1; ++J) { |
| 1657 |
2/2✓ Branch 0 taken 1600811 times.
✓ Branch 1 taken 1108813 times.
|
2709624 | if (tab[J] > tab[J + 1]) { |
| 1658 | std::swap(tab[J], tab[J + 1]); | ||
| 1659 | 1600811 | result = !result; | |
| 1660 | } | ||
| 1661 | } | ||
| 1662 | } | ||
| 1663 | 903208 | return result; | |
| 1664 | } | ||
| 1665 | |||
| 1666 |
2/2✓ Branch 0 taken 507740 times.
✓ Branch 1 taken 395468 times.
|
903208 | Sign ExactCDT2d::orient2d(index_t i, index_t j, index_t k) const { |
| 1667 | geo_debug_assert(i < nv()); | ||
| 1668 | geo_debug_assert(j < nv()); | ||
| 1669 | geo_debug_assert(k < nv()); | ||
| 1670 | |||
| 1671 | trindex K(i, j, k); | ||
| 1672 | |||
| 1673 |
2/2✓ Branch 0 taken 321731 times.
✓ Branch 1 taken 581477 times.
|
903208 | if(use_pred_cache_insert_buffer_) { |
| 1674 | 321731 | Sign result = PCK::orient_2d( | |
| 1675 | point_[K.indices[0]], | ||
| 1676 | point_[K.indices[1]], | ||
| 1677 | point_[K.indices[2]] | ||
| 1678 | ); | ||
| 1679 | 321731 | pred_cache_insert_buffer_.push_back(std::make_pair(K, result)); | |
| 1680 |
2/2✓ Branch 0 taken 137526 times.
✓ Branch 1 taken 184205 times.
|
321731 | if(odd_order(i,j,k)) { |
| 1681 | 137526 | result = Sign(-result); | |
| 1682 | } | ||
| 1683 | return result; | ||
| 1684 | } | ||
| 1685 | |||
| 1686 | bool inserted; | ||
| 1687 | std::map<trindex, Sign>::iterator it; | ||
| 1688 | 581477 | std::tie(it,inserted) = pred_cache_.insert(std::make_pair(K,ZERO)); | |
| 1689 | Sign result; | ||
| 1690 | |||
| 1691 |
2/2✓ Branch 0 taken 318333 times.
✓ Branch 1 taken 263144 times.
|
581477 | if(inserted) { |
| 1692 | 318333 | result = PCK::orient_2d( | |
| 1693 | point_[K.indices[0]], | ||
| 1694 | point_[K.indices[1]], | ||
| 1695 | point_[K.indices[2]] | ||
| 1696 | ); | ||
| 1697 | 318333 | it->second = result; | |
| 1698 | } else { | ||
| 1699 | 263144 | result = it->second; | |
| 1700 | } | ||
| 1701 | |||
| 1702 |
2/2✓ Branch 0 taken 266315 times.
✓ Branch 1 taken 315162 times.
|
581477 | if(odd_order(i,j,k)) { |
| 1703 | 266315 | result = Sign(-result); | |
| 1704 | } | ||
| 1705 | |||
| 1706 | return result; | ||
| 1707 | } | ||
| 1708 | |||
| 1709 | 156311 | Sign ExactCDT2d::incircle(index_t i,index_t j,index_t k,index_t l) const { | |
| 1710 | #ifdef GEOGRAM_USE_EXACT_NT | ||
| 1711 | return PCK::incircle_2d_SOS(point_[i], point_[j], point_[k], point_[l]); | ||
| 1712 | #else | ||
| 1713 | 156311 | return PCK::incircle_2d_SOS_with_lengths( | |
| 1714 | point_[i], point_[j], point_[k], point_[l], | ||
| 1715 | length_[i], length_[j], length_[k], length_[l] | ||
| 1716 | 156311 | ); | |
| 1717 | #endif | ||
| 1718 | } | ||
| 1719 | |||
| 1720 | 81 | index_t ExactCDT2d::create_intersection( | |
| 1721 | index_t E1, index_t i, index_t j, | ||
| 1722 | index_t E2, index_t k, index_t l | ||
| 1723 | ) { | ||
| 1724 | |||
| 1725 | geo_argused(i); | ||
| 1726 | geo_argused(j); | ||
| 1727 | geo_argused(k); | ||
| 1728 | geo_argused(l); | ||
| 1729 | |||
| 1730 | // Here we could use i,j,k,l directly, but it is *much better* to take | ||
| 1731 | // the original extremities of the constrained segments, since they have | ||
| 1732 | // simpler coordinates ! (i,j,k,l might be themselves vertices created | ||
| 1733 | // from constraints intersections, whereas constraint extremities can | ||
| 1734 | // only be initial vertices). | ||
| 1735 | 81 | i = constraints_[E1].indices[0]; | |
| 1736 | 81 | j = constraints_[E1].indices[1]; | |
| 1737 | 81 | k = constraints_[E2].indices[0]; | |
| 1738 | 81 | l = constraints_[E2].indices[1]; | |
| 1739 | |||
| 1740 | 81 | exact::vec2h U = point_[j] - point_[i]; | |
| 1741 |
1/2✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
|
81 | exact::vec2h V = point_[l] - point_[k]; |
| 1742 |
1/2✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
|
81 | exact::vec2h D = point_[k] - point_[i]; |
| 1743 | |||
| 1744 | exact::rational t( | ||
| 1745 |
2/6✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 81 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
162 | det2x2(D.x, D.y, V.x, V.y) * U.w, |
| 1746 |
1/2✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
|
81 | det2x2(U.x, U.y, V.x, V.y) * D.w |
| 1747 | ); | ||
| 1748 | |||
| 1749 |
2/4✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 81 times.
✗ Branch 6 not taken.
|
162 | point_.push_back(mix(t, point_[i], point_[j])); |
| 1750 | Numeric::optimize_number_representation(*point_.rbegin()); | ||
| 1751 | |||
| 1752 | #ifndef GEOGRAM_USE_EXACT_NT | ||
| 1753 | { | ||
| 1754 | const ExactPoint& p = *point_.rbegin(); | ||
| 1755 | 81 | length_.push_back( | |
| 1756 |
3/8✓ Branch 3 taken 81 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 81 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 81 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
405 | (geo_sqr(p.x) + geo_sqr(p.y)).estimate() / |
| 1757 |
1/4✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
162 | geo_sqr(p.w).estimate() |
| 1758 | ); | ||
| 1759 | } | ||
| 1760 | #endif | ||
| 1761 | |||
| 1762 | |||
| 1763 |
2/2✓ Branch 0 taken 77 times.
✓ Branch 1 taken 4 times.
|
81 | id_.push_back(NO_INDEX); |
| 1764 | 81 | index_t x = point_.size()-1; | |
| 1765 | |||
| 1766 |
2/2✓ Branch 0 taken 77 times.
✓ Branch 1 taken 4 times.
|
81 | CDTBase2d::v2T_.push_back(NO_INDEX); |
| 1767 | geo_debug_assert(x == CDTBase2d::nv_); | ||
| 1768 | 81 | ++CDTBase2d::nv_; | |
| 1769 | |||
| 1770 | 81 | return x; | |
| 1771 | 81 | } | |
| 1772 | |||
| 1773 | 106 | void ExactCDT2d::classify_triangles( | |
| 1774 | const std::string& expr, bool mark_only | ||
| 1775 | ) { | ||
| 1776 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 106 times.
|
106 | if(expr == "union_cnstr_operand_bits_is_operand_id") { |
| 1777 | ✗ | classify_triangles_union_cnstr_operand_bits_is_operand_id( | |
| 1778 | mark_only | ||
| 1779 | ); | ||
| 1780 | ✗ | return; | |
| 1781 | } | ||
| 1782 | |||
| 1783 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 106 times.
|
106 | facet_inclusion_bits_.assign(nT(), 0); |
| 1784 | |||
| 1785 | 106 | DList S(*this, DLIST_S_ID); | |
| 1786 | |||
| 1787 | // Step 1: get triangles adjacent to the border, | ||
| 1788 | // mark them as visited, classify them | ||
| 1789 |
2/2✓ Branch 0 taken 6534 times.
✓ Branch 1 taken 106 times.
|
13280 | for(index_t t=0; t<nT(); ++t) { |
| 1790 |
2/2✓ Branch 0 taken 18754 times.
✓ Branch 1 taken 6110 times.
|
24864 | for(index_t le=0; le<3; ++le) { |
| 1791 |
2/2✓ Branch 0 taken 424 times.
✓ Branch 1 taken 18330 times.
|
18754 | if(Tadj(t,le) == NO_INDEX) { |
| 1792 | Tset_flag(t, T_VISITED_FLAG); | ||
| 1793 | 424 | S.push_back(t); | |
| 1794 | break; | ||
| 1795 | } | ||
| 1796 | } | ||
| 1797 | } | ||
| 1798 | |||
| 1799 | // Step 2: recursive traversal | ||
| 1800 |
2/2✓ Branch 0 taken 6534 times.
✓ Branch 1 taken 106 times.
|
6640 | while(!S.empty()) { |
| 1801 | 6534 | index_t t1 = S.pop_back(); | |
| 1802 | 6534 | index_t t1_bits = facet_inclusion_bits_[t1]; | |
| 1803 |
2/2✓ Branch 0 taken 19602 times.
✓ Branch 1 taken 6534 times.
|
26136 | for(index_t le=0; le<3; ++le) { |
| 1804 | index_t t2 = Tadj(t1,le); | ||
| 1805 | if( | ||
| 1806 |
4/4✓ Branch 0 taken 19178 times.
✓ Branch 1 taken 424 times.
✓ Branch 2 taken 6110 times.
✓ Branch 3 taken 13068 times.
|
19602 | t2 != NO_INDEX && |
| 1807 | !Tflag_is_set(t2,T_VISITED_FLAG) | ||
| 1808 | ) { | ||
| 1809 | // t2 is included in the same operands as t1, | ||
| 1810 | // except for the operands that touch the boundary | ||
| 1811 | // between t1 and t2, for which inclusion changes | ||
| 1812 | index_t t2_bits = t1_bits; | ||
| 1813 | 6110 | for( | |
| 1814 | index_t ecit = Tedge_cnstr_first(t1,le); | ||
| 1815 |
2/2✓ Branch 0 taken 2390 times.
✓ Branch 1 taken 6110 times.
|
8500 | ecit != NO_INDEX; |
| 1816 | ecit = edge_cnstr_next(ecit) | ||
| 1817 | ) { | ||
| 1818 | index_t cnstr = edge_cnstr(ecit); | ||
| 1819 | 2390 | t2_bits ^= cnstr_operand_bits_[cnstr]; | |
| 1820 | } | ||
| 1821 | 6110 | facet_inclusion_bits_[t2] = t2_bits; | |
| 1822 | Tset_flag(t2, T_VISITED_FLAG); | ||
| 1823 | 6110 | S.push_back(t2); | |
| 1824 | } | ||
| 1825 | } | ||
| 1826 | } | ||
| 1827 | |||
| 1828 | // Step 3: reset visited flag | ||
| 1829 |
2/2✓ Branch 0 taken 6534 times.
✓ Branch 1 taken 106 times.
|
6640 | for(index_t t=0; t<nT(); ++t) { |
| 1830 | Treset_flag(t, T_VISITED_FLAG); | ||
| 1831 | } | ||
| 1832 | |||
| 1833 | // Step 4: mark triangles to be deleted | ||
| 1834 |
2/2✓ Branch 0 taken 9 times.
✓ Branch 1 taken 97 times.
|
106 | if(expr == "intersection") { |
| 1835 | index_t all_bits_set = 0; | ||
| 1836 |
2/2✓ Branch 0 taken 222 times.
✓ Branch 1 taken 9 times.
|
231 | for(index_t e_operand_bits: cnstr_operand_bits_) { |
| 1837 | 222 | all_bits_set |= e_operand_bits; | |
| 1838 | } | ||
| 1839 |
2/2✓ Branch 0 taken 498 times.
✓ Branch 1 taken 9 times.
|
1005 | for(index_t t=0; t<nT(); ++t) { |
| 1840 |
2/2✓ Branch 0 taken 431 times.
✓ Branch 1 taken 67 times.
|
498 | if(facet_inclusion_bits_[t] != all_bits_set) { |
| 1841 | Tset_flag(t, T_MARKED_FLAG); | ||
| 1842 | } | ||
| 1843 | } | ||
| 1844 | } else { | ||
| 1845 |
4/6✓ Branch 0 taken 88 times.
✓ Branch 1 taken 9 times.
✓ Branch 3 taken 88 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 97 times.
✗ Branch 7 not taken.
|
106 | BooleanExpression E(expr == "union" ? "*" : expr); |
| 1846 |
2/2✓ Branch 0 taken 6036 times.
✓ Branch 1 taken 97 times.
|
12266 | for(index_t t=0; t<nT(); ++t) { |
| 1847 |
3/4✓ Branch 1 taken 6036 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3364 times.
✓ Branch 4 taken 2672 times.
|
6036 | if(!E(facet_inclusion_bits_[t])) { |
| 1848 | Tset_flag(t, T_MARKED_FLAG); | ||
| 1849 | } | ||
| 1850 | } | ||
| 1851 | } | ||
| 1852 |
1/2✓ Branch 0 taken 106 times.
✗ Branch 1 not taken.
|
106 | if(!mark_only) { |
| 1853 |
1/2✓ Branch 1 taken 106 times.
✗ Branch 2 not taken.
|
106 | remove_marked_triangles(); |
| 1854 | } | ||
| 1855 | 106 | } | |
| 1856 | |||
| 1857 | typedef std::set<index_t> SparseBits; | ||
| 1858 | |||
| 1859 | ✗ | inline void sparse_bits_flip_bit(SparseBits& bits, index_t bit) { | |
| 1860 | auto it = bits.find(bit); | ||
| 1861 | ✗ | if(it == bits.end()) { | |
| 1862 | bits.insert(bit); | ||
| 1863 | } else { | ||
| 1864 | bits.erase(it); | ||
| 1865 | } | ||
| 1866 | ✗ | } | |
| 1867 | |||
| 1868 | inline bool sparse_bits_is_zero(const SparseBits& bits) { | ||
| 1869 | return (bits.size() == 0); | ||
| 1870 | } | ||
| 1871 | |||
| 1872 | ✗ | void ExactCDT2d::classify_triangles_union_cnstr_operand_bits_is_operand_id( | |
| 1873 | bool mark_only | ||
| 1874 | ) { | ||
| 1875 | |||
| 1876 | ✗ | DList S(*this, DLIST_S_ID); | |
| 1877 | std::stack<SparseBits> Sbits; | ||
| 1878 | |||
| 1879 | // Step 1: get triangles adjacent to the border, | ||
| 1880 | // mark them as visited, classify them as to-delete | ||
| 1881 | ✗ | for(index_t t=0; t<nT(); ++t) { | |
| 1882 | ✗ | for(index_t le=0; le<3; ++le) { | |
| 1883 | ✗ | if(Tadj(t,le) == NO_INDEX) { | |
| 1884 | Tset_flag(t, T_VISITED_FLAG); | ||
| 1885 | Tset_flag(t, T_MARKED_FLAG); | ||
| 1886 | ✗ | S.push_back(t); | |
| 1887 | ✗ | Sbits.push(SparseBits()); | |
| 1888 | ✗ | break; | |
| 1889 | } | ||
| 1890 | } | ||
| 1891 | } | ||
| 1892 | |||
| 1893 | // Step 2: recursive traversal | ||
| 1894 | ✗ | while(!S.empty()) { | |
| 1895 | ✗ | index_t t1 = S.pop_back(); | |
| 1896 | std::set<index_t> t1_bits = Sbits.top(); | ||
| 1897 | Sbits.pop(); | ||
| 1898 | ✗ | for(index_t le=0; le<3; ++le) { | |
| 1899 | index_t t2 = Tadj(t1,le); | ||
| 1900 | if( | ||
| 1901 | ✗ | t2 != NO_INDEX && | |
| 1902 | !Tflag_is_set(t2,T_VISITED_FLAG) | ||
| 1903 | ) { | ||
| 1904 | // t2 is included in the same operands as t1, | ||
| 1905 | // except for the operands that touch the boundary | ||
| 1906 | // between t1 and t2, for which inclusion changes | ||
| 1907 | SparseBits t2_bits = t1_bits; | ||
| 1908 | ✗ | for( | |
| 1909 | index_t ecit = Tedge_cnstr_first(t1,le); | ||
| 1910 | ✗ | ecit != NO_INDEX; | |
| 1911 | ecit = edge_cnstr_next(ecit) | ||
| 1912 | ) { | ||
| 1913 | index_t cnstr = edge_cnstr(ecit); | ||
| 1914 | ✗ | sparse_bits_flip_bit( | |
| 1915 | t2_bits, cnstr_operand_bits_[cnstr] | ||
| 1916 | ); | ||
| 1917 | } | ||
| 1918 | ✗ | if(sparse_bits_is_zero(t2_bits)) { | |
| 1919 | Tset_flag(t2, T_MARKED_FLAG); | ||
| 1920 | } | ||
| 1921 | Tset_flag(t2, T_VISITED_FLAG); | ||
| 1922 | ✗ | S.push_back(t2); | |
| 1923 | Sbits.push(t2_bits); | ||
| 1924 | } | ||
| 1925 | } | ||
| 1926 | } | ||
| 1927 | |||
| 1928 | // Step 3: reset visited flag | ||
| 1929 | ✗ | for(index_t t=0; t<nT(); ++t) { | |
| 1930 | Treset_flag(t, T_VISITED_FLAG); | ||
| 1931 | } | ||
| 1932 | |||
| 1933 | ✗ | if(!mark_only) { | |
| 1934 | ✗ | remove_marked_triangles(); | |
| 1935 | } | ||
| 1936 | ✗ | } | |
| 1937 | |||
| 1938 | ✗ | void ExactCDT2d::save(const std::string& filename) const { | |
| 1939 | #ifndef GEOGRAM_PSM | ||
| 1940 | ✗ | Mesh M; | |
| 1941 | ✗ | Attribute<index_t> nb_cnstr(M.edges.attributes(),"nb_cnstr"); | |
| 1942 | ✗ | M.vertices.set_dimension(2); | |
| 1943 | ✗ | for(const ExactPoint& P: point_) { | |
| 1944 | double w = P.w.estimate(); | ||
| 1945 | ✗ | vec2 p(P.x.estimate() / w, P.y.estimate() / w); | |
| 1946 | ✗ | M.vertices.create_vertex(p.data()); | |
| 1947 | } | ||
| 1948 | ✗ | for(index_t t=0; t<nT(); ++t) { | |
| 1949 | index_t i = Tv(t,0); | ||
| 1950 | index_t j = Tv(t,1); | ||
| 1951 | index_t k = Tv(t,2); | ||
| 1952 | ✗ | M.facets.create_triangle(i,j,k); | |
| 1953 | |||
| 1954 | ✗ | for(index_t le=0; le<3; ++le) { | |
| 1955 | ✗ | if(Tedge_is_constrained(t,le)) { | |
| 1956 | ✗ | index_t e = M.edges.create_edge( | |
| 1957 | ✗ | Tv(t,(le+1)%3), Tv(t,(le+2)%3) | |
| 1958 | ); | ||
| 1959 | ✗ | nb_cnstr[e] = Tedge_cnstr_nb(t,le); | |
| 1960 | } | ||
| 1961 | } | ||
| 1962 | } | ||
| 1963 | ✗ | M.facets.connect(); | |
| 1964 | ✗ | M.vertices.remove_isolated(); | |
| 1965 | ✗ | mesh_save(M, filename); | |
| 1966 | #else | ||
| 1967 | if(!String::string_ends_with(filename,".obj")) { | ||
| 1968 | Logger::err("CDT_2d") | ||
| 1969 | << "save() only supports .obj file format in PSM" | ||
| 1970 | << std::endl; | ||
| 1971 | return; | ||
| 1972 | } | ||
| 1973 | std::ofstream out(filename); | ||
| 1974 | for(const ExactPoint& P: point_) { | ||
| 1975 | double w = P.w.estimate(); | ||
| 1976 | vec2 p(P.x.estimate() / w, P.y.estimate() / w); | ||
| 1977 | out << "v " << p << " " << 0.0 << std::endl; | ||
| 1978 | } | ||
| 1979 | for(index_t t=0; t<nT(); ++t) { | ||
| 1980 | out << "f " << Tv(t,0)+1 << " " << Tv(t,1)+1 << " " << Tv(t,2)+1 | ||
| 1981 | << std::endl; | ||
| 1982 | } | ||
| 1983 | |||
| 1984 | for(index_t t=0; t<nT(); ++t) { | ||
| 1985 | for(index_t le=0; le<3; ++le) { | ||
| 1986 | if(Tedge_is_constrained(t,le)) { | ||
| 1987 | index_t v1 = Tv(t,(le+1)%3); | ||
| 1988 | index_t v2 = Tv(t,(le+2)%3); | ||
| 1989 | out << "l " << v1+1 << " " << v2+1 << std::endl; | ||
| 1990 | } | ||
| 1991 | } | ||
| 1992 | } | ||
| 1993 | #endif | ||
| 1994 | ✗ | } | |
| 1995 | |||
| 1996 | /***************************************************************************/ | ||
| 1997 | } | ||
| 1998 |