| 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 | 651 | CDTBase2d::CDTBase2d() : | |
| 98 | 651 | nv_(0), | |
| 99 | 651 | ncnstr_(0), | |
| 100 | 651 | delaunay_(true), | |
| 101 | 651 | exact_incircle_(true), | |
| 102 | 651 | exact_intersections_(true) { | |
| 103 | 651 | } | |
| 104 | |||
| 105 | 1302 | CDTBase2d::~CDTBase2d() { | |
| 106 | 1302 | } | |
| 107 | |||
| 108 | 23573 | void CDTBase2d::clear() { | |
| 109 | 23573 | nv_ = 0; | |
| 110 | 23573 | ncnstr_ = 0; | |
| 111 | 23573 | T_.resize(0); | |
| 112 | 23573 | Tadj_.resize(0); | |
| 113 | 23573 | v2T_.resize(0); | |
| 114 | 23573 | Tflags_.resize(0); | |
| 115 | 23573 | Tecnstr_first_.resize(0); | |
| 116 | 23573 | ecnstr_val_.resize(0); | |
| 117 | 23573 | ecnstr_next_.resize(0); | |
| 118 | 23573 | Tnext_.resize(0); | |
| 119 | 23573 | Tprev_.resize(0); | |
| 120 | 23573 | } | |
| 121 | |||
| 122 | 15538 | void CDTBase2d::create_enclosing_triangle( | |
| 123 | index_t v0, index_t v1, index_t v2 | ||
| 124 | ) { | ||
| 125 | 15538 | nv_ = 3; | |
| 126 | 15538 | v2T_.resize(3); | |
| 127 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 15538 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
15538 | geo_debug_assert(v0 <= 3); |
| 128 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 15538 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
15538 | geo_debug_assert(v1 <= 3); |
| 129 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 15538 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
15538 | geo_debug_assert(v2 <= 3); |
| 130 | 15538 | index_t t0 = Tnew(); | |
| 131 | 15538 | Tset(t0, v0, v1, v2, NO_INDEX, NO_INDEX, NO_INDEX); | |
| 132 | 15538 | orient_012_ = orient2d(0,1,2); | |
| 133 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 15538 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
15538 | geo_assert(orient_012_ != ZERO); |
| 134 | 15538 | } | |
| 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 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 8173 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
8173 | geo_debug_assert(v0 <= 4); |
| 142 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 8173 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
8173 | geo_debug_assert(v1 <= 4); |
| 143 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 8173 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
8173 | geo_debug_assert(v2 <= 4); |
| 144 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 8173 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
8173 | 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 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 8173 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
8173 | geo_debug_assert(is_convex_quad(t0)); |
| 151 |
2/2✓ Branch 1 taken 6130 times.
✓ Branch 2 taken 2043 times.
|
8173 | if(Sign(incircle(v0,v1,v2,v3)*orient_012_) == POSITIVE) { |
| 152 | 6130 | 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 | 515224 | index_t CDTBase2d::insert(index_t v, index_t hint) { | |
| 166 | 515224 | bool keep_duplicates = false; | |
| 167 |
2/2✓ Branch 1 taken 286741 times.
✓ Branch 2 taken 228483 times.
|
515224 | if(v == nv()) { |
| 168 |
1/2✓ Branch 1 taken 286741 times.
✗ Branch 2 not taken.
|
286741 | v2T_.push_back(NO_INDEX); |
| 169 | 286741 | ++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 | 228483 | keep_duplicates = true; | |
| 175 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 228483 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
228483 | geo_debug_assert(v < nv_); |
| 176 | } | ||
| 177 | |||
| 178 | // Phase 1: find triangle that contains vertex i | ||
| 179 |
1/2✓ Branch 1 taken 515224 times.
✗ Branch 2 not taken.
|
515224 | begin_insert_transaction(); |
| 180 | Sign o[3]; | ||
| 181 |
1/2✓ Branch 1 taken 515224 times.
✗ Branch 2 not taken.
|
515224 | index_t t = locate(v,hint,o); |
| 182 | 515224 | int nb_z = (o[0] == ZERO) + (o[1] == ZERO) + (o[2] == ZERO); | |
| 183 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 515224 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
515224 | geo_debug_assert(nb_z != 3); |
| 184 | |||
| 185 | // Duplicated vertex | ||
| 186 |
2/2✓ Branch 0 taken 157666 times.
✓ Branch 1 taken 357558 times.
|
515224 | if(nb_z == 2) { |
| 187 | CDT_LOG("duplicated vertex"); | ||
| 188 |
3/4✓ Branch 0 taken 69330 times.
✓ Branch 1 taken 88336 times.
✓ Branch 3 taken 69330 times.
✗ Branch 4 not taken.
|
157666 | v = (o[0] != ZERO) ? Tv(t,0) : |
| 189 |
3/4✓ Branch 0 taken 44171 times.
✓ Branch 1 taken 44165 times.
✓ Branch 3 taken 44171 times.
✗ Branch 4 not taken.
|
88336 | (o[1] != ZERO) ? Tv(t,1) : |
| 190 |
1/2✓ Branch 1 taken 44165 times.
✗ Branch 2 not taken.
|
44165 | Tv(t,2) ; |
| 191 |
2/2✓ Branch 0 taken 151895 times.
✓ Branch 1 taken 5771 times.
|
157666 | if(!keep_duplicates) { |
| 192 | 151895 | v2T_.pop_back(); | |
| 193 | 151895 | --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 |
1/2✓ Branch 1 taken 157666 times.
✗ Branch 2 not taken.
|
157666 | rollback_insert_transaction(); |
| 199 | 157666 | 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 |
1/2✓ Branch 1 taken 357558 times.
✗ Branch 2 not taken.
|
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 | 357558 | DList S(*this); | |
| 211 |
2/2✓ Branch 0 taken 246202 times.
✓ Branch 1 taken 111356 times.
|
357558 | if(delaunay_) { |
| 212 |
1/2✓ Branch 1 taken 246202 times.
✗ Branch 2 not taken.
|
246202 | 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 74588 times.
✓ Branch 1 taken 282970 times.
|
357558 | if(nb_z == 1) { |
| 218 | CDT_LOG("insert vertex on edge"); | ||
| 219 |
2/2✓ Branch 0 taken 34409 times.
✓ Branch 1 taken 40179 times.
|
108997 | index_t le = (o[0] == ZERO) ? 0 : |
| 220 |
2/2✓ Branch 0 taken 18211 times.
✓ Branch 1 taken 16198 times.
|
34409 | (o[1] == ZERO) ? 1 : |
| 221 | 2 ; | ||
| 222 |
1/2✓ Branch 1 taken 74588 times.
✗ Branch 2 not taken.
|
74588 | insert_vertex_in_edge(v,t,le,S); |
| 223 | } else { | ||
| 224 | CDT_LOG("insert vertex in triangle"); | ||
| 225 |
1/2✓ Branch 1 taken 282970 times.
✗ Branch 2 not taken.
|
282970 | 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 | 357558 | return v; | |
| 238 | 357558 | } | |
| 239 | |||
| 240 | |||
| 241 | 179299 | void CDTBase2d::insert_constraint(index_t i, index_t j) { | |
| 242 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 179299 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
179299 | geo_debug_assert(i < nv()); |
| 243 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 179299 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
179299 | geo_debug_assert(j < nv()); |
| 244 | CDT_LOG("insert constraint: " << i << "-" << j); | ||
| 245 | #ifdef CDT_DEBUG | ||
| 246 | debug_check_consistency(); | ||
| 247 | #endif | ||
| 248 | 179299 | ++ncnstr_; | |
| 249 | |||
| 250 | // Index of first vertex coming from constraints intersection | ||
| 251 | // (keep track of it to re-Delaunayize their neighborhoods). | ||
| 252 | 179299 | index_t first_v_isect = nv_; | |
| 253 | |||
| 254 | #ifndef CDT_NAIVE | ||
| 255 |
1/2✓ Branch 1 taken 179299 times.
✗ Branch 2 not taken.
|
179299 | DList Q(*this, DLIST_Q_ID); // Queue of edges to constrain |
| 256 | 179299 | DList N(*this); // New edges to re-Delaunayize (ignored if !delaunay_) | |
| 257 |
2/2✓ Branch 0 taken 178872 times.
✓ Branch 1 taken 427 times.
|
179299 | if(delaunay_) { |
| 258 |
1/2✓ Branch 1 taken 178872 times.
✗ Branch 2 not taken.
|
178872 | N.initialize(DLIST_N_ID); |
| 259 | } | ||
| 260 |
2/2✓ Branch 0 taken 258641 times.
✓ Branch 1 taken 179299 times.
|
437940 | 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 258641 times.
✗ Branch 2 not taken.
|
258641 | 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 258128 times.
✓ Branch 1 taken 513 times.
✓ Branch 2 taken 257615 times.
✓ Branch 3 taken 513 times.
✓ Branch 4 taken 5765 times.
✓ Branch 5 taken 251850 times.
|
258641 | if(delaunay_ && exact_intersections_ && k >= first_v_isect) { |
| 274 |
2/8✓ Branch 1 taken 5765 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 5765 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
5765 | geo_debug_assert(insert(k) == k); |
| 275 |
1/2✓ Branch 1 taken 5765 times.
✗ Branch 2 not taken.
|
5765 | 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/6✗ 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.
|
5765 | geo_assert(new_k == k); |
| 282 | } | ||
| 283 | |||
| 284 | // Step 2: constrain edges | ||
| 285 |
1/2✓ Branch 1 taken 258641 times.
✗ Branch 2 not taken.
|
258641 | constrain_edges(i,k,Q,N); |
| 286 | |||
| 287 | // Step 3: restore Delaunay condition | ||
| 288 |
2/2✓ Branch 0 taken 258128 times.
✓ Branch 1 taken 513 times.
|
258641 | if(delaunay_) { |
| 289 |
1/2✓ Branch 1 taken 258128 times.
✗ Branch 2 not taken.
|
258128 | Delaunayize_new_edges(N); |
| 290 | #ifdef CDT_DEBUG | ||
| 291 | debug_check_geometry(); | ||
| 292 | #endif | ||
| 293 | } | ||
| 294 | |||
| 295 | 258641 | 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 178872 times.
✓ Branch 1 taken 427 times.
✓ Branch 2 taken 427 times.
✓ Branch 3 taken 178445 times.
|
179299 | if(delaunay_ && !exact_intersections_) { |
| 336 |
2/2✓ Branch 1 taken 84 times.
✓ Branch 2 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 | 179299 | } | |
| 345 | |||
| 346 | 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 |
1/2✓ Branch 1 taken 5849 times.
✗ Branch 2 not taken.
|
5849 | DList S(*this, DLIST_S_ID); |
| 352 | |||
| 353 |
2/8✓ Branch 1 taken 5849 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 5849 times.
✗ 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 |
1/2✓ Branch 1 taken 5849 times.
✗ Branch 2 not taken.
|
5849 | index_t t0 = vT(v); // Need to store it, because we Trot() |
| 362 | 5849 | index_t t = t0; | |
| 363 | do { | ||
| 364 |
1/2✓ Branch 1 taken 23858 times.
✗ Branch 2 not taken.
|
23858 | index_t lv = Tv_find(t,v); |
| 365 |
1/2✓ Branch 1 taken 23858 times.
✗ Branch 2 not taken.
|
23858 | Trot(t,lv); |
| 366 |
2/8✓ Branch 1 taken 23858 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 23858 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
23858 | geo_debug_assert(Tv(t,0) == v); |
| 367 |
1/2✓ Branch 1 taken 23858 times.
✗ Branch 2 not taken.
|
23858 | S.push_back(t); |
| 368 |
1/2✓ Branch 1 taken 23858 times.
✗ Branch 2 not taken.
|
23858 | t = Tadj(t, 1); |
| 369 |
1/6✗ 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.
|
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 | 264406 | CDT2d_ConstraintWalker(index_t i_in, index_t j_in) : | |
| 388 | 264406 | i(i_in), j(j_in), | |
| 389 | 264406 | t_prev(NO_INDEX), v_prev(NO_INDEX), | |
| 390 | 264406 | t(NO_INDEX), v(i_in), | |
| 391 | 264406 | v_cnstr(NO_INDEX) | |
| 392 | { | ||
| 393 | 264406 | } | |
| 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 | 264406 | index_t CDTBase2d::find_intersected_edges(index_t i, index_t j, DList& Q) { | |
| 401 | CDT_LOG("Find intersected edges: " << i << "-" << j); | ||
| 402 | 264406 | CDT2d_ConstraintWalker W(i,j); | |
| 403 | // Stop at the first encountered vertex or constraint intersection. | ||
| 404 |
4/4✓ Branch 0 taken 264406 times.
✓ Branch 1 taken 324212 times.
✓ Branch 2 taken 59806 times.
✓ Branch 3 taken 264406 times.
|
588618 | 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 264406 times.
✓ Branch 1 taken 59806 times.
|
324212 | if(W.v != NO_INDEX) { |
| 411 |
1/2✓ Branch 1 taken 264406 times.
✗ Branch 2 not taken.
|
264406 | walk_constraint_v(W); |
| 412 | } else { | ||
| 413 |
1/2✓ Branch 1 taken 59806 times.
✗ Branch 2 not taken.
|
59806 | walk_constraint_t(W,Q); |
| 414 | } | ||
| 415 | } | ||
| 416 | 264406 | 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 | 264406 | void CDTBase2d::walk_constraint_v(CDT2d_ConstraintWalker& W) { | |
| 444 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 264406 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
264406 | geo_debug_assert(W.v != NO_INDEX); |
| 445 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 264406 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
264406 | geo_debug_assert(W.t == NO_INDEX); |
| 446 | |||
| 447 | 264406 | index_t t_next = NO_INDEX; | |
| 448 | 264406 | index_t v_next = NO_INDEX; | |
| 449 | |||
| 450 |
1/2✓ Branch 1 taken 264406 times.
✗ Branch 2 not taken.
|
264406 | for_each_T_around_v( |
| 451 |
1/2✓ Branch 1 taken 264406 times.
✗ Branch 2 not taken.
|
528812 | W.v, [&](index_t t_around_v, index_t le) { |
| 452 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 518601 times.
|
518601 | if(t_around_v == W.t_prev) { // Don't go backwards ! |
| 453 | ✗ | return false; | |
| 454 | } | ||
| 455 | 518601 | index_t v1 = Tv(t_around_v, (le + 1)%3); | |
| 456 | 518601 | index_t v2 = Tv(t_around_v, (le + 2)%3); | |
| 457 |
4/4✓ Branch 0 taken 429396 times.
✓ Branch 1 taken 89205 times.
✓ Branch 2 taken 72977 times.
✓ Branch 3 taken 356419 times.
|
518601 | if(v1 == W.j || v2 == W.j) { // Are we arrived at j ? |
| 458 | 162182 | 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 89205 times.
✓ Branch 1 taken 72977 times.
|
162182 | index_t le_cnstr_edge = (v1 == W.j) ? (le+2)%3 : (le+1)%3; |
| 462 | 162182 | Tadd_edge_cnstr_with_neighbor( | |
| 463 | 162182 | 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 | 162182 | return true; | |
| 472 | } | ||
| 473 | 356419 | Sign o1 = orient2d(W.i,W.j,v1); | |
| 474 | 356419 | Sign o2 = orient2d(W.i,W.j,v2); | |
| 475 | 356419 | Sign o3 = orient2d(v1,v2,W.j); | |
| 476 | 356419 | Sign o4 = orient_012_; // equivalent to orient2d(v1,v2,i) | |
| 477 |
4/4✓ Branch 0 taken 40644 times.
✓ Branch 1 taken 315775 times.
✓ Branch 2 taken 17018 times.
✓ Branch 3 taken 23626 times.
|
356419 | if(o1*o2 < 0 && o3*o4 < 0) { |
| 478 | 17018 | Trot(t_around_v,le); // so that le becomes edge 0 | |
| 479 | 17018 | t_next = t_around_v; // will be added to Q during next round | |
| 480 | 17018 | 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 |
3/8✓ Branch 0 taken 91560 times.
✓ Branch 1 taken 247841 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 91560 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
339401 | geo_debug_assert(o1 != ZERO || o2 != ZERO); |
| 486 |
5/6✓ Branch 0 taken 91560 times.
✓ Branch 1 taken 247841 times.
✓ Branch 2 taken 78540 times.
✓ Branch 3 taken 13020 times.
✓ Branch 4 taken 78540 times.
✗ Branch 5 not taken.
|
339401 | if(o1 == ZERO && o3*o4 < 0 && v1 != W.v_prev) { |
| 487 | 78540 | v_next = v1; | |
| 488 | 78540 | Tadd_edge_cnstr_with_neighbor( | |
| 489 | 78540 | t_around_v, (le + 2)%3, ncnstr_-1 | |
| 490 | ); | ||
| 491 | 78540 | return true; | |
| 492 |
5/6✓ Branch 0 taken 54947 times.
✓ Branch 1 taken 205914 times.
✓ Branch 2 taken 6666 times.
✓ Branch 3 taken 48281 times.
✓ Branch 4 taken 6666 times.
✗ Branch 5 not taken.
|
260861 | } else if(o2 == ZERO && o3*o4 < 0 && v2 != W.v_prev) { |
| 493 | 6666 | v_next = v2; | |
| 494 | 6666 | Tadd_edge_cnstr_with_neighbor( | |
| 495 | 6666 | t_around_v, (le + 1)%3, ncnstr_-1 | |
| 496 | ); | ||
| 497 | 6666 | return true; | |
| 498 | } | ||
| 499 | } | ||
| 500 | 254195 | return false; | |
| 501 | } | ||
| 502 | ); // End of for_each_T_around_v() loop | ||
| 503 | 264406 | W.t_prev = W.t; | |
| 504 | 264406 | W.v_prev = W.v; | |
| 505 | 264406 | W.t = t_next; | |
| 506 | 264406 | W.v = v_next; | |
| 507 | 264406 | } | |
| 508 | |||
| 509 | 59806 | void CDTBase2d::walk_constraint_t(CDT2d_ConstraintWalker& W, DList& Q) { | |
| 510 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 59806 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
59806 | geo_debug_assert(W.v == NO_INDEX); |
| 511 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 59806 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
59806 | geo_debug_assert(W.t != NO_INDEX); |
| 512 | |||
| 513 | 59806 | index_t v_next = NO_INDEX; | |
| 514 | 59806 | index_t t_next = NO_INDEX; | |
| 515 | |||
| 516 |
8/8✓ Branch 1 taken 57478 times.
✓ Branch 2 taken 2328 times.
✓ Branch 4 taken 53619 times.
✓ Branch 5 taken 3859 times.
✓ Branch 7 taken 3737 times.
✓ Branch 8 taken 49882 times.
✓ Branch 9 taken 9924 times.
✓ Branch 10 taken 49882 times.
|
59806 | if(Tv(W.t,0) == W.j || Tv(W.t,1) == W.j || Tv(W.t,2) == W.j) { |
| 517 | 9924 | 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 78028 times.
✗ Branch 1 not taken.
|
78028 | for(index_t le = 0; le<3; ++le) { |
| 521 |
2/2✓ Branch 1 taken 15717 times.
✓ Branch 2 taken 62311 times.
|
78028 | if(Tadj(W.t,le) == W.t_prev) { // Do not go backwards ! |
| 522 | 15717 | 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 | 62311 | index_t v1 = Tv(W.t, (le + 1)%3); | |
| 529 | 62311 | index_t v2 = Tv(W.t, (le + 2)%3); | |
| 530 | 62311 | Sign o1 = orient2d(W.i,W.j,v1); | |
| 531 | 62311 | Sign o2 = orient2d(W.i,W.j,v2); | |
| 532 |
2/2✓ Branch 0 taken 48721 times.
✓ Branch 1 taken 13590 times.
|
62311 | if(o1*o2 < 0) { |
| 533 | // [v1,v2] has a frank intersection with [i,j] | ||
| 534 | 48721 | Trot(W.t,le); // So that edge 0 is intersected edge | |
| 535 |
2/2✓ Branch 1 taken 5933 times.
✓ Branch 2 taken 42788 times.
|
48721 | if(Tedge_is_constrained(W.t,0)) { |
| 536 | CDT_LOG(" Cnstr isect with:" << v1 << "-" << v2); | ||
| 537 | 11866 | v_next = create_intersection( | |
| 538 | 5933 | 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 3392 times.
✓ Branch 1 taken 2541 times.
|
5933 | if(W.v_prev != NO_INDEX) { |
| 545 | 3392 | Tadd_edge_cnstr_with_neighbor(W.t,2,ncnstr_-1); | |
| 546 | } | ||
| 547 | } else { | ||
| 548 | CDT_LOG(" Isect: t=" << W.t <<" E=" << v1 <<"-"<< v2); | ||
| 549 | 42788 | Q.push_back(W.t); | |
| 550 | 42788 | t_next = Tadj(W.t,0); | |
| 551 | } | ||
| 552 | 48721 | break; | |
| 553 | } else { // Special case: v1 or v2 is exactly on [i,j] | ||
| 554 |
3/8✓ Branch 0 taken 8 times.
✓ Branch 1 taken 13582 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 8 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
13590 | geo_debug_assert(o1 != ZERO || o2 != ZERO); |
| 555 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 13582 times.
|
13590 | if(o1 == ZERO) { |
| 556 | 8 | v_next = v1; | |
| 557 | 8 | break; | |
| 558 |
2/2✓ Branch 0 taken 1153 times.
✓ Branch 1 taken 12429 times.
|
13582 | } else if(o2 == ZERO) { |
| 559 | 1153 | v_next = v2; | |
| 560 | 1153 | break; | |
| 561 | } | ||
| 562 | } | ||
| 563 | } | ||
| 564 | } | ||
| 565 | 59806 | W.t_prev = W.t; | |
| 566 | 59806 | W.v_prev = W.v; | |
| 567 | 59806 | W.t = t_next; | |
| 568 | 59806 | W.v = v_next; | |
| 569 | 59806 | } | |
| 570 | |||
| 571 | 258641 | 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 | 38938 | auto new_edge = [&](index_t t,index_t le) { | |
| 582 | 38938 | Trot(t,le); | |
| 583 | 38938 | if( | |
| 584 |
6/6✓ Branch 1 taken 12736 times.
✓ Branch 2 taken 26202 times.
✓ Branch 4 taken 1506 times.
✓ Branch 5 taken 11230 times.
✓ Branch 6 taken 11230 times.
✓ Branch 7 taken 27708 times.
|
66646 | (Tv(t,1) == i && Tv(t,2) == j) || |
| 585 |
1/4✗ Branch 1 not taken.
✓ Branch 2 taken 27708 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
27708 | (Tv(t,1) == j && Tv(t,2) == i) |
| 586 | ) { | ||
| 587 | // Set constraint flag if the new edge is the constrained edge | ||
| 588 | 11230 | Tadd_edge_cnstr_with_neighbor(t,0,ncnstr_-1); | |
| 589 | } else { | ||
| 590 | // Memorize new edge as "to be Delaunayized" | ||
| 591 |
2/2✓ Branch 1 taken 11340 times.
✓ Branch 2 taken 16368 times.
|
27708 | if(N.initialized()) { |
| 592 | 11340 | N.push_back(t); | |
| 593 | } | ||
| 594 | } | ||
| 595 | 297579 | }; | |
| 596 | |||
| 597 | // Called each time edge le of triangle t still has an isect with cnstr | ||
| 598 | // (then it is queued again) | ||
| 599 | 18311 | auto isect_edge = [&](index_t t, index_t le) { | |
| 600 | 18311 | Trot(t,le); | |
| 601 | 18311 | Q.push_front(t); | |
| 602 | 276952 | }; | |
| 603 | |||
| 604 |
3/4✓ Branch 1 taken 337816 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 79175 times.
✓ Branch 4 taken 258641 times.
|
337816 | while(!Q.empty()) { |
| 605 |
1/2✓ Branch 1 taken 79175 times.
✗ Branch 2 not taken.
|
79175 | index_t t1 = Q.pop_back(); |
| 606 |
3/4✓ Branch 1 taken 79175 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 21926 times.
✓ Branch 4 taken 57249 times.
|
79175 | 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 |
2/8✓ Branch 1 taken 21926 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 21926 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
21926 | geo_assert(!Q.empty()); |
| 611 |
1/2✓ Branch 1 taken 21926 times.
✗ Branch 2 not taken.
|
21926 | Q.push_front(t1); |
| 612 | } else { | ||
| 613 |
1/2✓ Branch 1 taken 57249 times.
✗ Branch 2 not taken.
|
57249 | index_t t2 = Tadj(t1,0); |
| 614 |
1/2✓ Branch 1 taken 57249 times.
✗ Branch 2 not taken.
|
57249 | bool no_isect = !Q.contains(t2); |
| 615 |
1/2✓ Branch 1 taken 57249 times.
✗ Branch 2 not taken.
|
57249 | index_t v0 = Tv(t1,0); |
| 616 |
7/10✓ Branch 1 taken 57249 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 38501 times.
✓ Branch 4 taken 18748 times.
✓ Branch 6 taken 38501 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 38501 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 18900 times.
✓ Branch 12 taken 19601 times.
|
57249 | bool t2v0_t1v2 = (Q.contains(t2) && Tv(t2,0) == Tv(t1,2)); |
| 617 |
7/10✓ Branch 1 taken 57249 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 38501 times.
✓ Branch 4 taken 18748 times.
✓ Branch 6 taken 38501 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 38501 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 19601 times.
✓ Branch 12 taken 18900 times.
|
57249 | bool t2v0_t1v1 = (Q.contains(t2) && Tv(t2,0) == Tv(t1,1)); |
| 618 | 57249 | geo_argused(t2v0_t1v1); | |
| 619 | |||
| 620 |
2/2✓ Branch 0 taken 18748 times.
✓ Branch 1 taken 38501 times.
|
57249 | if(no_isect) { |
| 621 |
1/2✓ Branch 1 taken 18748 times.
✗ Branch 2 not taken.
|
18748 | swap_edge(t1); |
| 622 |
2/8✓ Branch 1 taken 18748 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 18748 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
18748 | geo_debug_assert(!segment_edge_intersect(i,j,t1,2)); |
| 623 |
1/2✓ Branch 1 taken 18748 times.
✗ Branch 2 not taken.
|
18748 | 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 |
1/2✓ Branch 1 taken 38501 times.
✗ Branch 2 not taken.
|
38501 | Sign o = Sign(orient2d(i,j,v0) * orient_012_); |
| 629 |
2/2✓ Branch 0 taken 18900 times.
✓ Branch 1 taken 19601 times.
|
38501 | if(t2v0_t1v2) { |
| 630 |
1/2✓ Branch 1 taken 18900 times.
✗ Branch 2 not taken.
|
18900 | swap_edge(t1,false); // "new t1 on top" |
| 631 |
2/2✓ Branch 0 taken 9614 times.
✓ Branch 1 taken 9286 times.
|
18900 | if(o >= 0) { |
| 632 |
2/8✓ Branch 1 taken 9614 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9614 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
9614 | geo_debug_assert(!segment_edge_intersect(i,j,t1,2)); |
| 633 |
2/8✓ Branch 1 taken 9614 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9614 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
9614 | geo_debug_assert( segment_edge_intersect(i,j,t2,0)); |
| 634 |
1/2✓ Branch 1 taken 9614 times.
✗ Branch 2 not taken.
|
9614 | new_edge(t1,2); |
| 635 | } else { | ||
| 636 |
2/8✓ Branch 1 taken 9286 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9286 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
9286 | geo_debug_assert( segment_edge_intersect(i,j,t1,2)); |
| 637 |
2/8✓ Branch 1 taken 9286 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9286 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
9286 | geo_debug_assert( segment_edge_intersect(i,j,t2,0)); |
| 638 |
1/2✓ Branch 1 taken 9286 times.
✗ Branch 2 not taken.
|
9286 | isect_edge(t1,2); |
| 639 | } | ||
| 640 | } else { | ||
| 641 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 19601 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
19601 | geo_debug_assert(t2v0_t1v1); |
| 642 |
1/2✓ Branch 1 taken 19601 times.
✗ Branch 2 not taken.
|
19601 | swap_edge(t1,true); // "new t1 on bottom" |
| 643 |
2/2✓ Branch 0 taken 9025 times.
✓ Branch 1 taken 10576 times.
|
19601 | if(o > 0) { |
| 644 |
2/8✓ Branch 1 taken 9025 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9025 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
9025 | geo_debug_assert( segment_edge_intersect(i,j,t1,1)); |
| 645 |
2/8✓ Branch 1 taken 9025 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9025 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
9025 | geo_debug_assert( segment_edge_intersect(i,j,t2,0)); |
| 646 |
1/2✓ Branch 1 taken 9025 times.
✗ Branch 2 not taken.
|
9025 | isect_edge(t1,1); |
| 647 | } else { | ||
| 648 |
2/8✓ Branch 1 taken 10576 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 10576 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
10576 | geo_debug_assert(!segment_edge_intersect(i,j,t1,1)); |
| 649 |
2/8✓ Branch 1 taken 10576 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 10576 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
10576 | geo_debug_assert( segment_edge_intersect(i,j,t2,0)); |
| 650 |
1/2✓ Branch 1 taken 10576 times.
✗ Branch 2 not taken.
|
10576 | new_edge(t1,1); |
| 651 | } | ||
| 652 | } | ||
| 653 | } | ||
| 654 | } | ||
| 655 | } | ||
| 656 | 258641 | } | |
| 657 | |||
| 658 | 252051 | void CDTBase2d::Delaunayize_vertex_neighbors(index_t v, DList& S) { | |
| 659 | CDT_LOG("Delaunayize_vertex_neighbors"); | ||
| 660 | 252051 | index_t count = 0; | |
| 661 |
2/2✓ Branch 1 taken 1704836 times.
✓ Branch 2 taken 252051 times.
|
1956887 | while(!S.empty()) { |
| 662 | // NASA programming style: all loops have | ||
| 663 | // a maximum number of iterations | ||
| 664 | 1704836 | ++count; | |
| 665 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 1704836 times.
|
1704836 | 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 | 1704836 | index_t t1 = S.pop_back(); | |
| 673 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 1704836 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
1704836 | geo_debug_assert(Tv(t1,0) == v); |
| 674 |
2/2✓ Branch 1 taken 9667 times.
✓ Branch 2 taken 1695169 times.
|
1704836 | if(Tedge_is_constrained(t1,0)) { |
| 675 | 9667 | continue; | |
| 676 | } | ||
| 677 | 1695169 | index_t t2 = Tadj(t1,0); | |
| 678 |
2/2✓ Branch 0 taken 204183 times.
✓ Branch 1 taken 1490986 times.
|
1695169 | if(t2 == NO_INDEX) { |
| 679 | 204183 | continue; | |
| 680 | } | ||
| 681 |
5/6✓ Branch 0 taken 1490986 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 230536 times.
✓ Branch 4 taken 1260450 times.
✓ Branch 5 taken 230536 times.
✓ Branch 6 taken 1260450 times.
|
1490986 | if(!exact_incircle_ && !is_convex_quad(t1)) { |
| 682 | 230536 | continue; | |
| 683 | } | ||
| 684 | |||
| 685 | 1260450 | index_t v1 = Tv(t2,0); | |
| 686 | 1260450 | index_t v2 = Tv(t2,1); | |
| 687 | 1260450 | index_t v3 = Tv(t2,2); | |
| 688 |
2/2✓ Branch 1 taken 506247 times.
✓ Branch 2 taken 754203 times.
|
1260450 | if(Sign(incircle(v1,v2,v3,v)*orient_012_) == POSITIVE) { |
| 689 | 506247 | swap_edge(t1); | |
| 690 | 506247 | S.push_back(t1); | |
| 691 | 506247 | S.push_back(t2); | |
| 692 | } | ||
| 693 | } | ||
| 694 | CDT_LOG("/Delaunayize_vertex_neighbors"); | ||
| 695 | 252051 | } | |
| 696 | |||
| 697 | 258128 | void CDTBase2d::Delaunayize_new_edges(DList& N) { | |
| 698 | 258128 | index_t count = 0; | |
| 699 | 258128 | bool swap_occured = true; | |
| 700 |
2/2✓ Branch 0 taken 260093 times.
✓ Branch 1 taken 258128 times.
|
518221 | while(swap_occured) { |
| 701 | 260093 | swap_occured = false; | |
| 702 | // NASA programming style: all loops have | ||
| 703 | // a maximum number of iterations | ||
| 704 | 260093 | ++count; | |
| 705 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 260093 times.
|
260093 | 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 2 taken 22921 times.
✓ Branch 3 taken 260093 times.
|
283014 | for(index_t t1 = N.front(); t1 != NO_INDEX; t1 = N.next(t1)) { |
| 712 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 22921 times.
|
22921 | if(Tedge_is_constrained(t1,0)) { |
| 713 | ✗ | continue; | |
| 714 | } | ||
| 715 | 22921 | index_t v1 = Tv(t1,1); | |
| 716 | 22921 | index_t v2 = Tv(t1,2); | |
| 717 | 22921 | index_t v0 = Tv(t1,0); | |
| 718 | 22921 | index_t t2 = Tadj(t1,0); | |
| 719 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 22921 times.
|
22921 | if(t2 == NO_INDEX) { |
| 720 | ✗ | continue; | |
| 721 | } | ||
| 722 |
5/6✓ Branch 0 taken 22921 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 12281 times.
✓ Branch 4 taken 10640 times.
✓ Branch 5 taken 12281 times.
✓ Branch 6 taken 10640 times.
|
22921 | if(!exact_incircle_ && !is_convex_quad(t1)) { |
| 723 | 12281 | continue; | |
| 724 | } | ||
| 725 | 10640 | index_t e2 = Tadj_find(t2,t1); | |
| 726 | 10640 | index_t v3 = Tv(t2,e2); | |
| 727 |
2/2✓ Branch 1 taken 2944 times.
✓ Branch 2 taken 7696 times.
|
10640 | 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 2 taken 2393 times.
✓ Branch 3 taken 551 times.
|
2944 | if(Tv(t2,0) == Tv(t1,1)) { |
| 731 | 2393 | swap_edge(t1, true); // t2 on top | |
| 732 | 2393 | Trot(t1,1); | |
| 733 | } else { | ||
| 734 | 551 | swap_edge(t1, false); // t1 on top | |
| 735 | 551 | Trot(t1,2); | |
| 736 | } | ||
| 737 | 2944 | swap_occured = true; | |
| 738 | } | ||
| 739 | } | ||
| 740 | } | ||
| 741 | 258128 | N.clear(); | |
| 742 | 258128 | } | |
| 743 | |||
| 744 | 515224 | 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 515224 times.
|
515224 | if(o == nullptr) { |
| 747 | ✗ | o = o_local; | |
| 748 | } | ||
| 749 | |||
| 750 | // Efficient locate, "walking the triangulation" | ||
| 751 |
1/2✓ Branch 1 taken 515224 times.
✗ Branch 2 not taken.
|
515224 | index_t t_pred = nT()+1; // Needs to be different from NO_INDEX |
| 752 |
2/2✓ Branch 0 taken 292538 times.
✓ Branch 1 taken 222686 times.
|
515224 | index_t t = (hint == NO_INDEX) ? |
| 753 |
2/4✓ Branch 1 taken 292538 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 292538 times.
✗ Branch 5 not taken.
|
292538 | index_t(Numeric::random_int32()) % nT() : |
| 754 | 515224 | hint ; | |
| 755 | #ifdef GEO_DEBUG | ||
| 756 | 515224 | index_t nb_traversed_t = 0; | |
| 757 | #endif | ||
| 758 | |||
| 759 | 6121594 | still_walking: | |
| 760 | { | ||
| 761 | #ifdef GEO_DEBUG | ||
| 762 | 6121594 | ++nb_traversed_t; | |
| 763 | #endif | ||
| 764 | |||
| 765 | // Infinite loop are not supposed to happen, but | ||
| 766 | // let us detect them, just in case... | ||
| 767 |
2/8✓ Branch 1 taken 6121594 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 6121594 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
6121594 | 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 | 6121594 | bool point_outside_boundary = (t == NO_INDEX); | |
| 772 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 6121594 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
6121594 | geo_assert(!point_outside_boundary); |
| 773 | |||
| 774 | index_t tv[3]; | ||
| 775 |
1/2✓ Branch 1 taken 6121594 times.
✗ Branch 2 not taken.
|
6121594 | tv[0] = Tv(t,0); |
| 776 |
1/2✓ Branch 1 taken 6121594 times.
✗ Branch 2 not taken.
|
6121594 | tv[1] = Tv(t,1); |
| 777 |
1/2✓ Branch 1 taken 6121594 times.
✗ Branch 2 not taken.
|
6121594 | tv[2] = Tv(t,2); |
| 778 | |||
| 779 | // Start from a random edge | ||
| 780 |
1/2✓ Branch 1 taken 6121594 times.
✗ Branch 2 not taken.
|
6121594 | index_t e0 = index_t(Numeric::random_int32()) % 3; |
| 781 |
2/2✓ Branch 0 taken 11877250 times.
✓ Branch 1 taken 515224 times.
|
12392474 | for(index_t de = 0; de < 3; ++de) { |
| 782 | 11877250 | index_t le = (e0 + de) % 3; | |
| 783 | |||
| 784 |
1/2✓ Branch 1 taken 11877250 times.
✗ Branch 2 not taken.
|
11877250 | 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 2846392 times.
✓ Branch 1 taken 9030858 times.
|
11877250 | if(t_next == t_pred) { |
| 795 | 2846392 | o[le] = POSITIVE; | |
| 796 | 2846392 | 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 | 9030858 | index_t v_bkp = tv[le]; | |
| 803 | 9030858 | tv[le] = v; | |
| 804 |
1/2✓ Branch 1 taken 9030858 times.
✗ Branch 2 not taken.
|
9030858 | 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 3424488 times.
✓ Branch 1 taken 5606370 times.
|
9030858 | if(o[le] != NEGATIVE) { |
| 810 | 3424488 | tv[le] = v_bkp; | |
| 811 | 3424488 | continue; | |
| 812 | } | ||
| 813 | |||
| 814 | // If we reach this point, then t_next is a valid | ||
| 815 | // successor, thus we are still walking. | ||
| 816 | 5606370 | t_pred = t; | |
| 817 | 5606370 | t = t_next; | |
| 818 | 5606370 | goto still_walking; | |
| 819 | } | ||
| 820 | } | ||
| 821 | |||
| 822 | 515224 | 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 |
1/2✓ Branch 1 taken 8067 times.
✗ Branch 2 not taken.
|
8067 | DList S(*this, DLIST_S_ID); |
| 829 | |||
| 830 | // Step 1: get triangles adjacent to the border, | ||
| 831 | // mark them as visited, classify them | ||
| 832 |
3/4✓ Branch 1 taken 111147 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 103080 times.
✓ Branch 4 taken 8067 times.
|
111147 | for(index_t t=0; t<nT(); ++t) { |
| 833 |
2/2✓ Branch 0 taken 244778 times.
✓ Branch 1 taken 70812 times.
|
315590 | for(index_t le=0; le<3; ++le) { |
| 834 |
3/4✓ Branch 1 taken 244778 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 32268 times.
✓ Branch 4 taken 212510 times.
|
244778 | if(Tadj(t,le) == NO_INDEX) { |
| 835 |
1/2✓ Branch 1 taken 32268 times.
✗ Branch 2 not taken.
|
32268 | bool outside = ((Tedge_cnstr_nb(t,le)%2) == 0); |
| 836 |
1/2✓ Branch 1 taken 32268 times.
✗ Branch 2 not taken.
|
32268 | Tset_flag(t, T_VISITED_FLAG); |
| 837 |
1/2✓ Branch 0 taken 32268 times.
✗ Branch 1 not taken.
|
32268 | if(outside) { |
| 838 |
1/2✓ Branch 1 taken 32268 times.
✗ Branch 2 not taken.
|
32268 | Tset_flag(t, T_MARKED_FLAG); |
| 839 | } | ||
| 840 |
1/2✓ Branch 1 taken 32268 times.
✗ Branch 2 not taken.
|
32268 | S.push_back(t); |
| 841 | 32268 | break; | |
| 842 | } | ||
| 843 | } | ||
| 844 | } | ||
| 845 | |||
| 846 | // Step 2: recursive traversal | ||
| 847 |
3/4✓ Branch 1 taken 111147 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 103080 times.
✓ Branch 4 taken 8067 times.
|
111147 | while(!S.empty()) { |
| 848 |
1/2✓ Branch 1 taken 103080 times.
✗ Branch 2 not taken.
|
103080 | index_t t1 = S.pop_back(); |
| 849 |
1/2✓ Branch 1 taken 103080 times.
✗ Branch 2 not taken.
|
103080 | bool t1_outside = Tflag_is_set(t1, T_MARKED_FLAG); |
| 850 |
2/2✓ Branch 0 taken 309240 times.
✓ Branch 1 taken 103080 times.
|
412320 | for(index_t le=0; le<3; ++le) { |
| 851 |
1/2✓ Branch 1 taken 309240 times.
✗ Branch 2 not taken.
|
309240 | index_t t2 = Tadj(t1,le); |
| 852 | 309240 | if( | |
| 853 |
4/4✓ Branch 0 taken 276972 times.
✓ Branch 1 taken 32268 times.
✓ Branch 2 taken 70812 times.
✓ Branch 3 taken 238428 times.
|
586212 | t2 != NO_INDEX && |
| 854 |
3/4✓ Branch 1 taken 276972 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 70812 times.
✓ Branch 4 taken 206160 times.
|
276972 | !Tflag_is_set(t2,T_VISITED_FLAG) |
| 855 | ) { | ||
| 856 | bool t2_outside = | ||
| 857 |
1/2✓ Branch 1 taken 70812 times.
✗ Branch 2 not taken.
|
70812 | t1_outside ^ ((Tedge_cnstr_nb(t1,le)%2) != 0); |
| 858 |
1/2✓ Branch 1 taken 70812 times.
✗ Branch 2 not taken.
|
70812 | Tset_flag(t2, T_VISITED_FLAG); |
| 859 |
2/2✓ Branch 0 taken 43269 times.
✓ Branch 1 taken 27543 times.
|
70812 | if(t2_outside) { |
| 860 |
1/2✓ Branch 1 taken 43269 times.
✗ Branch 2 not taken.
|
43269 | Tset_flag(t2, T_MARKED_FLAG); |
| 861 | } | ||
| 862 |
1/2✓ Branch 1 taken 70812 times.
✗ Branch 2 not taken.
|
70812 | S.push_back(t2); |
| 863 | } | ||
| 864 | } | ||
| 865 | } | ||
| 866 | |||
| 867 | // Step 3: reset visited flag | ||
| 868 |
3/4✓ Branch 1 taken 111147 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 103080 times.
✓ Branch 4 taken 8067 times.
|
111147 | for(index_t t=0; t<nT(); ++t) { |
| 869 |
1/2✓ Branch 1 taken 103080 times.
✗ Branch 2 not taken.
|
103080 | Treset_flag(t, T_VISITED_FLAG); |
| 870 | } | ||
| 871 | |||
| 872 | 8067 | } else { | |
| 873 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
4 | DList S(*this, DLIST_S_ID); |
| 874 | |||
| 875 | // Step 1: get triangles adjacent to the border | ||
| 876 |
3/4✓ Branch 1 taken 40656 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 40652 times.
✓ Branch 4 taken 4 times.
|
40656 | for(index_t t=0; t<nT(); ++t) { |
| 877 | 40652 | if( | |
| 878 | // TODO: replace with parity check ? | ||
| 879 |
3/6✓ Branch 1 taken 40652 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 40396 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 40396 times.
✗ Branch 7 not taken.
|
81048 | (!Tedge_is_constrained(t,0) && Tadj(t,0) == NO_INDEX) || |
| 880 |
8/12✓ Branch 0 taken 40396 times.
✓ Branch 1 taken 256 times.
✓ Branch 3 taken 40652 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 40295 times.
✓ Branch 6 taken 357 times.
✓ Branch 8 taken 40295 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 40295 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 40652 times.
|
121700 | (!Tedge_is_constrained(t,1) && Tadj(t,1) == NO_INDEX) || |
| 881 |
5/8✓ Branch 1 taken 40652 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 40445 times.
✓ Branch 4 taken 207 times.
✓ Branch 6 taken 40445 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 40445 times.
|
40652 | (!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 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 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 | 8177 | vector<index_t>& old2new = Tnext_; | |
| 914 | 8177 | index_t cur_t_new = 0; | |
| 915 |
2/2✓ Branch 1 taken 150266 times.
✓ Branch 2 taken 8177 times.
|
158443 | for(index_t t=0; t<nT(); ++t) { |
| 916 |
2/2✓ Branch 1 taken 79332 times.
✓ Branch 2 taken 70934 times.
|
150266 | if(Tflag_is_set(t,T_MARKED_FLAG)) { |
| 917 | 79332 | old2new[t] = NO_INDEX; | |
| 918 | } else { | ||
| 919 | 70934 | old2new[t] = cur_t_new; | |
| 920 | 70934 | ++cur_t_new; | |
| 921 | } | ||
| 922 | } | ||
| 923 | 8177 | index_t nT_new = cur_t_new; | |
| 924 | |||
| 925 | // Step 2: translate adjacency and move triangles | ||
| 926 |
2/2✓ Branch 1 taken 150266 times.
✓ Branch 2 taken 8177 times.
|
158443 | for(index_t t=0; t<nT(); ++t) { |
| 927 | 150266 | index_t t_new = old2new[t]; | |
| 928 |
2/2✓ Branch 0 taken 79332 times.
✓ Branch 1 taken 70934 times.
|
150266 | if(t_new == NO_INDEX) { |
| 929 | 79332 | continue; | |
| 930 | } | ||
| 931 | 70934 | index_t adj0 = Tadj(t,0); | |
| 932 |
2/2✓ Branch 0 taken 70922 times.
✓ Branch 1 taken 12 times.
|
70934 | if(adj0 != NO_INDEX) { |
| 933 | 70922 | adj0 = old2new[adj0]; | |
| 934 | } | ||
| 935 | 70934 | index_t adj1 = Tadj(t,1); | |
| 936 |
1/2✓ Branch 0 taken 70934 times.
✗ Branch 1 not taken.
|
70934 | if(adj1 != NO_INDEX) { |
| 937 | 70934 | adj1 = old2new[adj1]; | |
| 938 | } | ||
| 939 | 70934 | index_t adj2 = Tadj(t,2); | |
| 940 |
1/2✓ Branch 0 taken 70934 times.
✗ Branch 1 not taken.
|
70934 | if(adj2 != NO_INDEX) { |
| 941 | 70934 | adj2 = old2new[adj2]; | |
| 942 | } | ||
| 943 | 70934 | 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 | 70934 | 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 1 taken 70934 times.
✓ Branch 2 taken 8177 times.
|
79111 | for(index_t t=0; t<nT(); ++t) { |
| 964 | 70934 | v2T_[Tv(t,0)] = t; | |
| 965 | 70934 | v2T_[Tv(t,1)] = t; | |
| 966 | 70934 | v2T_[Tv(t,2)] = t; | |
| 967 | } | ||
| 968 | 8177 | } | |
| 969 | |||
| 970 | |||
| 971 | /***************** Triangulation surgery (boring code ahead) *********/ | ||
| 972 | |||
| 973 | 80521 | void CDTBase2d::insert_vertex_in_edge( | |
| 974 | index_t v, index_t t, index_t le1, DList& S | ||
| 975 | ) { | ||
| 976 | 80521 | index_t cnstr_first = Tedge_cnstr_first(t,le1); | |
| 977 | 80521 | index_t t1 = t; | |
| 978 | 80521 | index_t t2 = Tadj(t1,le1); | |
| 979 | 80521 | index_t v1 = Tv(t1,le1); | |
| 980 | 80521 | index_t v2 = Tv(t1,(le1+1)%3); | |
| 981 | 80521 | index_t v3 = Tv(t1,(le1+2)%3); | |
| 982 | 80521 | index_t t1_adj2 = Tadj(t1,(le1+1)%3); | |
| 983 | 80521 | index_t t1_adj3 = Tadj(t1,(le1+2)%3); | |
| 984 |
2/2✓ Branch 0 taken 8071 times.
✓ Branch 1 taken 72450 times.
|
80521 | 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 | 8071 | index_t le2 = Tadj_find(t2,t1); | |
| 991 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 8071 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
8071 | geo_debug_assert(Tv(t2, (le2+1)%3) == v3); |
| 992 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 8071 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
8071 | geo_debug_assert(Tv(t2, (le2+2)%3) == v2); |
| 993 | 8071 | index_t v4 = Tv(t2,le2); | |
| 994 | 8071 | index_t t2_adj2 = Tadj(t2,(le2+1)%3); | |
| 995 | 8071 | index_t t2_adj3 = Tadj(t2,(le2+2)%3); | |
| 996 | 8071 | index_t t3 = Tnew(); | |
| 997 | 8071 | index_t t4 = Tnew(); | |
| 998 | 8071 | Tset(t1,v,v1,v2,t1_adj3,t2,t4); | |
| 999 | 8071 | Tset(t2,v,v2,v4,t2_adj2,t3,t1); | |
| 1000 | 8071 | Tset(t3,v,v4,v3,t2_adj3,t4,t2); | |
| 1001 | 8071 | Tset(t4,v,v3,v1,t1_adj2,t1,t3); | |
| 1002 | 8071 | Tadj_back_connect(t1,0,t1); | |
| 1003 | 8071 | Tadj_back_connect(t2,0,t2); | |
| 1004 | 8071 | Tadj_back_connect(t3,0,t2); | |
| 1005 | 8071 | Tadj_back_connect(t4,0,t1); | |
| 1006 | 8071 | Tset_edge_cnstr_first(t1,1,cnstr_first); | |
| 1007 | 8071 | Tset_edge_cnstr_first(t2,2,cnstr_first); | |
| 1008 | 8071 | Tset_edge_cnstr_first(t3,1,cnstr_first); | |
| 1009 | 8071 | Tset_edge_cnstr_first(t4,2,cnstr_first); | |
| 1010 |
2/2✓ Branch 1 taken 2137 times.
✓ Branch 2 taken 5934 times.
|
8071 | if(S.initialized()) { |
| 1011 | 2137 | S.push_back(t1); | |
| 1012 | 2137 | S.push_back(t2); | |
| 1013 | 2137 | S.push_back(t3); | |
| 1014 | 2137 | 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 | 72450 | t2 = Tnew(); | |
| 1022 | 72450 | Tset(t1,v,v1,v2,t1_adj3,NO_INDEX,t2); | |
| 1023 | 72450 | Tset(t2,v,v3,v1,t1_adj2,t1,NO_INDEX); | |
| 1024 | 72450 | Tadj_back_connect(t1,0,t1); | |
| 1025 | 72450 | Tadj_back_connect(t2,0,t1); | |
| 1026 | 72450 | Tset_edge_cnstr_first(t1,1,cnstr_first); | |
| 1027 | 72450 | Tset_edge_cnstr_first(t2,2,cnstr_first); | |
| 1028 |
2/2✓ Branch 1 taken 72259 times.
✓ Branch 2 taken 191 times.
|
72450 | if(S.initialized()) { |
| 1029 | 72259 | S.push_back(t1); | |
| 1030 | 72259 | S.push_back(t2); | |
| 1031 | } | ||
| 1032 | } | ||
| 1033 | 80521 | } | |
| 1034 | |||
| 1035 | 282970 | 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 | 282970 | index_t t1 = t; | |
| 1039 | 282970 | index_t v1 = Tv(t1,0); | |
| 1040 | 282970 | index_t v2 = Tv(t1,1); | |
| 1041 | 282970 | index_t v3 = Tv(t1,2); | |
| 1042 | 282970 | index_t adj1 = Tadj(t1,0); | |
| 1043 | 282970 | index_t adj2 = Tadj(t1,1); | |
| 1044 | 282970 | index_t adj3 = Tadj(t1,2); | |
| 1045 | 282970 | index_t t2 = Tnew(); | |
| 1046 | 282970 | index_t t3 = Tnew(); | |
| 1047 | 282970 | Tset(t1,v,v2,v3,adj1,t2,t3); | |
| 1048 | 282970 | Tset(t2,v,v3,v1,adj2,t3,t1); | |
| 1049 | 282970 | Tset(t3,v,v1,v2,adj3,t1,t2); | |
| 1050 | 282970 | Tadj_back_connect(t1,0,t1); | |
| 1051 | 282970 | Tadj_back_connect(t2,0,t1); | |
| 1052 | 282970 | Tadj_back_connect(t3,0,t1); | |
| 1053 |
2/2✓ Branch 1 taken 171806 times.
✓ Branch 2 taken 111164 times.
|
282970 | if(S.initialized()) { |
| 1054 | 171806 | S.push_back(t1); | |
| 1055 | 171806 | S.push_back(t2); | |
| 1056 | 171806 | S.push_back(t3); | |
| 1057 | } | ||
| 1058 | 282970 | } | |
| 1059 | |||
| 1060 | 572570 | void CDTBase2d::swap_edge(index_t t1, bool swap_t1_t2) { | |
| 1061 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 572570 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
572570 | geo_debug_assert(!Tedge_is_constrained(t1,0)); |
| 1062 | 572570 | index_t v1 = Tv(t1,0); | |
| 1063 | 572570 | index_t v2 = Tv(t1,1); | |
| 1064 | 572570 | index_t v3 = Tv(t1,2); | |
| 1065 | 572570 | index_t t1_adj2 = Tadj(t1,1); | |
| 1066 | 572570 | index_t t1_adj3 = Tadj(t1,2); | |
| 1067 | 572570 | index_t t2 = Tadj(t1,0); | |
| 1068 | 572570 | index_t le2 = Tadj_find(t2,t1); | |
| 1069 | 572570 | index_t v4 = Tv(t2,le2); | |
| 1070 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 572570 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
572570 | geo_debug_assert(Tv(t2,(le2+1)%3) == v3); |
| 1071 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 572570 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
572570 | geo_debug_assert(Tv(t2,(le2+2)%3) == v2); |
| 1072 | |||
| 1073 | 572570 | debug_Tcheck(t1); | |
| 1074 | 572570 | debug_Tcheck(t2); | |
| 1075 | |||
| 1076 | 572570 | index_t t2_adj2 = Tadj(t2,(le2+1)%3); | |
| 1077 | 572570 | index_t t2_adj3 = Tadj(t2,(le2+2)%3); | |
| 1078 |
2/2✓ Branch 0 taken 21994 times.
✓ Branch 1 taken 550576 times.
|
572570 | if(swap_t1_t2) { |
| 1079 | 21994 | Tset(t2,v1,v4,v3,t2_adj3,t1_adj2,t1); | |
| 1080 | 21994 | Tset(t1,v1,v2,v4,t2_adj2,t2,t1_adj3); | |
| 1081 | 21994 | Tadj_back_connect(t2,0,t2); | |
| 1082 | 21994 | Tadj_back_connect(t2,1,t1); | |
| 1083 | 21994 | Tadj_back_connect(t1,0,t2); | |
| 1084 | 21994 | Tadj_back_connect(t1,2,t1); | |
| 1085 | } else { | ||
| 1086 | 550576 | Tset(t1,v1,v4,v3,t2_adj3,t1_adj2,t2); | |
| 1087 | 550576 | Tset(t2,v1,v2,v4,t2_adj2,t1,t1_adj3); | |
| 1088 | 550576 | Tadj_back_connect(t1,0,t2); | |
| 1089 | 550576 | Tadj_back_connect(t1,1,t1); | |
| 1090 | 550576 | Tadj_back_connect(t2,0,t2); | |
| 1091 | 550576 | Tadj_back_connect(t2,2,t1); | |
| 1092 | } | ||
| 1093 | |||
| 1094 | 572570 | debug_Tcheck(t1); | |
| 1095 | 572570 | debug_Tcheck(t2); | |
| 1096 | 572570 | } | |
| 1097 | |||
| 1098 | /***************** Geometry ***********************/ | ||
| 1099 | |||
| 1100 | 1601255 | bool CDTBase2d::is_convex_quad(index_t t) const { | |
| 1101 | 1601255 | index_t v1 = Tv(t,0); | |
| 1102 | 1601255 | index_t v2 = Tv(t,1); | |
| 1103 | 1601255 | index_t v3 = Tv(t,2); | |
| 1104 | 1601255 | index_t t2 = Tadj(t,0); | |
| 1105 | 1601255 | index_t le2 = Tadj_find(t2,t); | |
| 1106 | 1601255 | 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 1473330 times.
✓ Branch 2 taken 127925 times.
|
3074585 | orient2d(v1,v4,v3) == orient_012_ && |
| 1113 |
2/2✓ Branch 1 taken 1336512 times.
✓ Branch 2 taken 136818 times.
|
3074585 | orient2d(v4,v1,v2) == orient_012_ ; |
| 1114 | } | ||
| 1115 | |||
| 1116 | /******* Debugging ******************************************************/ | ||
| 1117 | |||
| 1118 | 93328 | void CDTBase2d::check_geometry() const { | |
| 1119 |
3/4✓ Branch 0 taken 93264 times.
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 93264 times.
|
93328 | 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 | 93328 | } | |
| 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 | 64 | CDT2d::~CDT2d() { | |
| 1307 | 64 | } | |
| 1308 | |||
| 1309 | ✗ | void CDT2d::clear() { | |
| 1310 | ✗ | CDTBase2d::clear(); | |
| 1311 | ✗ | point_.resize(0); | |
| 1312 | ✗ | } | |
| 1313 | |||
| 1314 | 32 | void CDT2d::create_enclosing_triangle( | |
| 1315 | const vec2& p1, const vec2& p2, const vec2& p3 | ||
| 1316 | ) { | ||
| 1317 |
1/6✗ Branch 1 not taken.
✓ Branch 2 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 1 not taken.
✓ Branch 2 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 | 32 | point_.push_back(p1); | |
| 1320 | 32 | point_.push_back(p2); | |
| 1321 | 32 | 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 | 9081215 | Sign CDT2d::orient2d(index_t i, index_t j, index_t k) const { | |
| 1338 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 9081215 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
9081215 | geo_debug_assert(i < nv()); |
| 1339 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 9081215 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
9081215 | geo_debug_assert(j < nv()); |
| 1340 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 9081215 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
9081215 | geo_debug_assert(k < nv()); |
| 1341 | 9081215 | return PCK::orient_2d( | |
| 1342 | 9081215 | point_[i].data(), point_[j].data(), point_[k].data() | |
| 1343 | 9081215 | ); | |
| 1344 | } | ||
| 1345 | |||
| 1346 | 935458 | Sign CDT2d::incircle(index_t i, index_t j, index_t k, index_t l) const { | |
| 1347 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 935458 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
935458 | geo_debug_assert(i < nv()); |
| 1348 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 935458 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
935458 | geo_debug_assert(j < nv()); |
| 1349 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 935458 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
935458 | geo_debug_assert(k < nv()); |
| 1350 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 935458 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
935458 | geo_debug_assert(l < nv()); |
| 1351 | 935458 | return PCK::in_circle_2d_SOS( | |
| 1352 | 935458 | point_[i].data(), point_[j].data(), point_[k].data(), | |
| 1353 | 935458 | point_[l].data() | |
| 1354 | 935458 | ); | |
| 1355 | } | ||
| 1356 | |||
| 1357 | 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 | 168 | geo_argused(E1); | |
| 1362 | 168 | geo_argused(E2); | |
| 1363 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 168 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
168 | geo_debug_assert(i < nv()); |
| 1364 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 168 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
168 | geo_debug_assert(j < nv()); |
| 1365 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 168 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
168 | geo_debug_assert(k < nv()); |
| 1366 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 168 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
168 | geo_debug_assert(l < nv()); |
| 1367 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 168 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
168 | geo_debug_assert(E1 < ncnstr()); |
| 1368 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 168 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
168 | geo_debug_assert(E2 < ncnstr()); |
| 1369 |
2/4✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 168 times.
✗ Branch 5 not taken.
|
168 | vec2 U = point_[j] - point_[i]; |
| 1370 |
2/4✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 168 times.
✗ Branch 5 not taken.
|
168 | vec2 V = point_[l] - point_[k]; |
| 1371 |
2/4✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 168 times.
✗ Branch 5 not taken.
|
168 | vec2 D = point_[k] - point_[i]; |
| 1372 | 168 | double delta = det(U,V); | |
| 1373 | 168 | double t = det(D,V)/delta; | |
| 1374 |
1/2✓ Branch 2 taken 168 times.
✗ Branch 3 not taken.
|
168 | vec2 P = point_[i] + t*U; |
| 1375 |
1/2✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
|
168 | point_.push_back(P); |
| 1376 |
1/2✓ Branch 1 taken 168 times.
✗ Branch 2 not taken.
|
168 | v2T_.push_back(NO_INDEX); |
| 1377 | 168 | index_t v = nv_; | |
| 1378 | 168 | ++nv_; | |
| 1379 | 168 | return v; | |
| 1380 | } | ||
| 1381 | |||
| 1382 | 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 |
1/2✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
|
32 | debug_check_consistency(); |
| 1388 | |||
| 1389 | // Compute spatial sort | ||
| 1390 | 32 | 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 | 32 | index_t v_offset = nv(); | |
| 1421 |
1/2✓ Branch 2 taken 32 times.
✗ Branch 3 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/2✓ Branch 2 taken 222718 times.
✗ Branch 3 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 |
1/2✓ Branch 2 taken 32 times.
✗ Branch 3 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 | 32 | 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 | 668154 | index_t v = CDTBase2d::insert( | |
| 1440 |
2/4✓ Branch 1 taken 222718 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 222718 times.
✗ Branch 5 not taken.
|
222718 | v_offset+sorted_indices[i], hint |
| 1441 | ); | ||
| 1442 |
1/2✓ Branch 0 taken 222718 times.
✗ Branch 1 not taken.
|
222718 | if(indices != nullptr) { |
| 1443 |
1/2✓ Branch 1 taken 222718 times.
✗ Branch 2 not taken.
|
222718 | indices[sorted_indices[i]] = v; |
| 1444 | } | ||
| 1445 |
1/2✓ Branch 1 taken 222718 times.
✗ Branch 2 not taken.
|
222718 | hint = vT(v); |
| 1446 | } | ||
| 1447 | } | ||
| 1448 | CDT_LOG("Inserted."); | ||
| 1449 |
1/2✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
|
32 | 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 | 782 | ExactCDT2d::~ExactCDT2d() { | |
| 1547 | 782 | } | |
| 1548 | |||
| 1549 | 8067 | void ExactCDT2d::clear() { | |
| 1550 | 8067 | point_.resize(0); | |
| 1551 | 8067 | id_.resize(0); | |
| 1552 | 8067 | 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 | 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 1 not taken.
✓ Branch 2 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 1 not taken.
✓ Branch 2 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 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 8173 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
8173 | 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 | 46600 | index_t ExactCDT2d::insert(const ExactPoint& p, index_t id, index_t hint) { | |
| 1594 | #ifndef GEOGRAM_USE_EXACT_NT | ||
| 1595 |
1/6✗ Branch 2 not taken.
✓ Branch 3 taken 46600 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
46600 | geo_debug_assert(nv() == length_.size()); |
| 1596 | #endif | ||
| 1597 | 46600 | debug_check_consistency(); | |
| 1598 | 46600 | add_point(p,id); | |
| 1599 | 46600 | 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 2 taken 47 times.
✓ Branch 3 taken 46553 times.
|
46600 | if(point_.size() > nv()) { |
| 1603 | 47 | point_.pop_back(); | |
| 1604 | 47 | id_.pop_back(); | |
| 1605 | #ifndef GEOGRAM_USE_EXACT_NT | ||
| 1606 | 47 | length_.pop_back(); | |
| 1607 | #endif | ||
| 1608 | } | ||
| 1609 | 46600 | debug_check_consistency(); | |
| 1610 | #ifndef GEOGRAM_USE_EXACT_NT | ||
| 1611 |
1/6✗ Branch 2 not taken.
✓ Branch 3 taken 46600 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
46600 | geo_debug_assert(nv() == length_.size()); |
| 1612 | #endif | ||
| 1613 | 46600 | return v; | |
| 1614 | } | ||
| 1615 | |||
| 1616 | 79292 | void ExactCDT2d::add_point(const ExactPoint& p, index_t id) { | |
| 1617 | 79292 | point_.push_back(p); | |
| 1618 | 79292 | id_.push_back(id); | |
| 1619 | #ifndef GEOGRAM_USE_EXACT_NT | ||
| 1620 |
2/4✓ Branch 1 taken 79292 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 79292 times.
✗ Branch 5 not taken.
|
79292 | length_.push_back(squared_length(p.x, p.y, p.w)); |
| 1621 | #endif | ||
| 1622 | 79292 | } | |
| 1623 | |||
| 1624 | 46681 | void ExactCDT2d::begin_insert_transaction() { | |
| 1625 | 46681 | use_pred_cache_insert_buffer_ = true; | |
| 1626 | 46681 | } | |
| 1627 | |||
| 1628 | 46553 | void ExactCDT2d::commit_insert_transaction() { | |
| 1629 |
2/2✓ Branch 2 taken 323266 times.
✓ Branch 3 taken 46553 times.
|
416372 | for(const auto& it: pred_cache_insert_buffer_) { |
| 1630 |
1/2✓ Branch 1 taken 323266 times.
✗ Branch 2 not taken.
|
323266 | pred_cache_[it.first] = it.second; |
| 1631 | } | ||
| 1632 | 46553 | pred_cache_insert_buffer_.resize(0); | |
| 1633 | 46553 | use_pred_cache_insert_buffer_ = false; | |
| 1634 | 46553 | } | |
| 1635 | |||
| 1636 | 128 | void ExactCDT2d::rollback_insert_transaction() { | |
| 1637 | 128 | pred_cache_insert_buffer_.resize(0); | |
| 1638 | 128 | use_pred_cache_insert_buffer_ = false; | |
| 1639 | 128 | } | |
| 1640 | |||
| 1641 | /** | ||
| 1642 | * \brief Tests the parity of the permutation of a list of | ||
| 1643 | * three distinct indices with respect to the canonical order. | ||
| 1644 | */ | ||
| 1645 | 1288897 | static bool odd_order(index_t i, index_t j, index_t k) { | |
| 1646 | // Implementation: sort the elements (bubble sort is OK for | ||
| 1647 | // such a small number), and invert parity each time | ||
| 1648 | // two elements are swapped. | ||
| 1649 | 1288897 | index_t tab[3] = { i, j, k}; | |
| 1650 | 1288897 | const int N = 3; | |
| 1651 | 1288897 | bool result = false; | |
| 1652 |
2/2✓ Branch 0 taken 2577794 times.
✓ Branch 1 taken 1288897 times.
|
3866691 | for (int I = 0; I < N - 1; ++I) { |
| 1653 |
2/2✓ Branch 0 taken 3866691 times.
✓ Branch 1 taken 2577794 times.
|
6444485 | for (int J = 0; J < N - I - 1; ++J) { |
| 1654 |
2/2✓ Branch 0 taken 2448819 times.
✓ Branch 1 taken 1417872 times.
|
3866691 | if (tab[J] > tab[J + 1]) { |
| 1655 | 2448819 | std::swap(tab[J], tab[J + 1]); | |
| 1656 | 2448819 | result = !result; | |
| 1657 | } | ||
| 1658 | } | ||
| 1659 | } | ||
| 1660 | 1288897 | return result; | |
| 1661 | } | ||
| 1662 | |||
| 1663 | 1288897 | Sign ExactCDT2d::orient2d(index_t i, index_t j, index_t k) const { | |
| 1664 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 1288897 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
1288897 | geo_debug_assert(i < nv()); |
| 1665 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 1288897 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
1288897 | geo_debug_assert(j < nv()); |
| 1666 |
1/6✗ Branch 1 not taken.
✓ Branch 2 taken 1288897 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
1288897 | geo_debug_assert(k < nv()); |
| 1667 | |||
| 1668 |
1/2✓ Branch 1 taken 1288897 times.
✗ Branch 2 not taken.
|
1288897 | trindex K(i, j, k); |
| 1669 | |||
| 1670 |
2/2✓ Branch 0 taken 324811 times.
✓ Branch 1 taken 964086 times.
|
1288897 | if(use_pred_cache_insert_buffer_) { |
| 1671 | 649622 | Sign result = PCK::orient_2d( | |
| 1672 |
1/2✓ Branch 1 taken 324811 times.
✗ Branch 2 not taken.
|
324811 | point_[K.indices[0]], |
| 1673 |
1/2✓ Branch 1 taken 324811 times.
✗ Branch 2 not taken.
|
324811 | point_[K.indices[1]], |
| 1674 |
2/4✓ Branch 1 taken 324811 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 324811 times.
✗ Branch 5 not taken.
|
324811 | point_[K.indices[2]] |
| 1675 | 324811 | ); | |
| 1676 |
2/4✓ Branch 1 taken 324811 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 324811 times.
✗ Branch 5 not taken.
|
324811 | pred_cache_insert_buffer_.push_back(std::make_pair(K, result)); |
| 1677 |
2/2✓ Branch 1 taken 139041 times.
✓ Branch 2 taken 185770 times.
|
324811 | if(odd_order(i,j,k)) { |
| 1678 | 139041 | result = Sign(-result); | |
| 1679 | } | ||
| 1680 | 324811 | return result; | |
| 1681 | } | ||
| 1682 | |||
| 1683 | bool inserted; | ||
| 1684 | 964086 | std::map<trindex, Sign>::iterator it; | |
| 1685 |
2/4✓ Branch 1 taken 964086 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 964086 times.
✗ Branch 5 not taken.
|
964086 | std::tie(it,inserted) = pred_cache_.insert(std::make_pair(K,ZERO)); |
| 1686 | Sign result; | ||
| 1687 | |||
| 1688 |
2/2✓ Branch 0 taken 346204 times.
✓ Branch 1 taken 617882 times.
|
964086 | if(inserted) { |
| 1689 | 692408 | result = PCK::orient_2d( | |
| 1690 |
1/2✓ Branch 1 taken 346204 times.
✗ Branch 2 not taken.
|
346204 | point_[K.indices[0]], |
| 1691 |
1/2✓ Branch 1 taken 346204 times.
✗ Branch 2 not taken.
|
346204 | point_[K.indices[1]], |
| 1692 |
2/4✓ Branch 1 taken 346204 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 346204 times.
✗ Branch 5 not taken.
|
346204 | point_[K.indices[2]] |
| 1693 | ); | ||
| 1694 | 346204 | it->second = result; | |
| 1695 | } else { | ||
| 1696 | 617882 | result = it->second; | |
| 1697 | } | ||
| 1698 | |||
| 1699 |
2/2✓ Branch 1 taken 415574 times.
✓ Branch 2 taken 548512 times.
|
964086 | if(odd_order(i,j,k)) { |
| 1700 | 415574 | result = Sign(-result); | |
| 1701 | } | ||
| 1702 | |||
| 1703 | 964086 | return result; | |
| 1704 | } | ||
| 1705 | |||
| 1706 | 156493 | Sign ExactCDT2d::incircle(index_t i,index_t j,index_t k,index_t l) const { | |
| 1707 | #ifdef GEOGRAM_USE_EXACT_NT | ||
| 1708 | return PCK::incircle_2d_SOS(point_[i], point_[j], point_[k], point_[l]); | ||
| 1709 | #else | ||
| 1710 | 938958 | return PCK::incircle_2d_SOS_with_lengths( | |
| 1711 | 156493 | point_[i], point_[j], point_[k], point_[l], | |
| 1712 | 156493 | length_[i], length_[j], length_[k], length_[l] | |
| 1713 | 156493 | ); | |
| 1714 | #endif | ||
| 1715 | } | ||
| 1716 | |||
| 1717 | 81 | index_t ExactCDT2d::create_intersection( | |
| 1718 | index_t E1, index_t i, index_t j, | ||
| 1719 | index_t E2, index_t k, index_t l | ||
| 1720 | ) { | ||
| 1721 | |||
| 1722 | 81 | geo_argused(i); | |
| 1723 | 81 | geo_argused(j); | |
| 1724 | 81 | geo_argused(k); | |
| 1725 | 81 | geo_argused(l); | |
| 1726 | |||
| 1727 | // Here we could use i,j,k,l directly, but it is *much better* to take | ||
| 1728 | // the original extremities of the constrained segments, since they have | ||
| 1729 | // simpler coordinates ! (i,j,k,l might be themselves vertices created | ||
| 1730 | // from constraints intersections, whereas constraint extremities can | ||
| 1731 | // only be initial vertices). | ||
| 1732 |
1/2✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
|
81 | i = constraints_[E1].indices[0]; |
| 1733 |
1/2✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
|
81 | j = constraints_[E1].indices[1]; |
| 1734 |
1/2✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
|
81 | k = constraints_[E2].indices[0]; |
| 1735 |
1/2✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
|
81 | l = constraints_[E2].indices[1]; |
| 1736 | |||
| 1737 |
3/6✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 81 times.
✗ Branch 8 not taken.
|
81 | exact::vec2h U = point_[j] - point_[i]; |
| 1738 |
3/6✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 81 times.
✗ Branch 8 not taken.
|
81 | exact::vec2h V = point_[l] - point_[k]; |
| 1739 |
3/6✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 81 times.
✗ Branch 8 not taken.
|
81 | exact::vec2h D = point_[k] - point_[i]; |
| 1740 | |||
| 1741 | exact::rational t( | ||
| 1742 |
2/4✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
|
162 | det2x2(D.x, D.y, V.x, V.y) * U.w, |
| 1743 |
2/4✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
|
162 | det2x2(U.x, U.y, V.x, V.y) * D.w |
| 1744 |
1/2✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
|
81 | ); |
| 1745 | |||
| 1746 |
4/8✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 81 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 81 times.
✗ Branch 11 not taken.
|
81 | point_.push_back(mix(t, point_[i], point_[j])); |
| 1747 |
1/2✓ Branch 3 taken 81 times.
✗ Branch 4 not taken.
|
81 | Numeric::optimize_number_representation(*point_.rbegin()); |
| 1748 | |||
| 1749 | #ifndef GEOGRAM_USE_EXACT_NT | ||
| 1750 | { | ||
| 1751 | 81 | const ExactPoint& p = *point_.rbegin(); | |
| 1752 | 162 | length_.push_back( | |
| 1753 |
6/12✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 81 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 81 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 81 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 81 times.
✗ Branch 17 not taken.
|
162 | (geo_sqr(p.x) + geo_sqr(p.y)).estimate() / |
| 1754 |
4/8✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 81 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 81 times.
✗ Branch 11 not taken.
|
81 | geo_sqr(p.w).estimate() |
| 1755 | ); | ||
| 1756 | } | ||
| 1757 | #endif | ||
| 1758 | |||
| 1759 | |||
| 1760 |
1/2✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
|
81 | id_.push_back(NO_INDEX); |
| 1761 | 81 | index_t x = point_.size()-1; | |
| 1762 | |||
| 1763 |
1/2✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
|
81 | CDTBase2d::v2T_.push_back(NO_INDEX); |
| 1764 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 81 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
81 | geo_debug_assert(x == CDTBase2d::nv_); |
| 1765 | 81 | ++CDTBase2d::nv_; | |
| 1766 | |||
| 1767 | 81 | return x; | |
| 1768 | 81 | } | |
| 1769 | |||
| 1770 | 106 | void ExactCDT2d::classify_triangles( | |
| 1771 | const std::string& expr, bool mark_only | ||
| 1772 | ) { | ||
| 1773 |
2/4✓ Branch 1 taken 106 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 106 times.
|
106 | if(expr == "union_cnstr_operand_bits_is_operand_id") { |
| 1774 | ✗ | classify_triangles_union_cnstr_operand_bits_is_operand_id( | |
| 1775 | mark_only | ||
| 1776 | ); | ||
| 1777 | ✗ | return; | |
| 1778 | } | ||
| 1779 | |||
| 1780 |
2/4✓ Branch 1 taken 106 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 106 times.
✗ Branch 5 not taken.
|
106 | facet_inclusion_bits_.assign(nT(), 0); |
| 1781 | |||
| 1782 |
1/2✓ Branch 1 taken 106 times.
✗ Branch 2 not taken.
|
106 | DList S(*this, DLIST_S_ID); |
| 1783 | |||
| 1784 | // Step 1: get triangles adjacent to the border, | ||
| 1785 | // mark them as visited, classify them | ||
| 1786 |
3/4✓ Branch 1 taken 6640 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6534 times.
✓ Branch 4 taken 106 times.
|
6640 | for(index_t t=0; t<nT(); ++t) { |
| 1787 |
2/2✓ Branch 0 taken 18754 times.
✓ Branch 1 taken 6110 times.
|
24864 | for(index_t le=0; le<3; ++le) { |
| 1788 |
3/4✓ Branch 1 taken 18754 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 424 times.
✓ Branch 4 taken 18330 times.
|
18754 | if(Tadj(t,le) == NO_INDEX) { |
| 1789 |
1/2✓ Branch 1 taken 424 times.
✗ Branch 2 not taken.
|
424 | Tset_flag(t, T_VISITED_FLAG); |
| 1790 |
1/2✓ Branch 1 taken 424 times.
✗ Branch 2 not taken.
|
424 | S.push_back(t); |
| 1791 | 424 | break; | |
| 1792 | } | ||
| 1793 | } | ||
| 1794 | } | ||
| 1795 | |||
| 1796 | // Step 2: recursive traversal | ||
| 1797 |
3/4✓ Branch 1 taken 6640 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6534 times.
✓ Branch 4 taken 106 times.
|
6640 | while(!S.empty()) { |
| 1798 |
1/2✓ Branch 1 taken 6534 times.
✗ Branch 2 not taken.
|
6534 | index_t t1 = S.pop_back(); |
| 1799 |
1/2✓ Branch 1 taken 6534 times.
✗ Branch 2 not taken.
|
6534 | index_t t1_bits = facet_inclusion_bits_[t1]; |
| 1800 |
2/2✓ Branch 0 taken 19602 times.
✓ Branch 1 taken 6534 times.
|
26136 | for(index_t le=0; le<3; ++le) { |
| 1801 |
1/2✓ Branch 1 taken 19602 times.
✗ Branch 2 not taken.
|
19602 | index_t t2 = Tadj(t1,le); |
| 1802 | 19602 | if( | |
| 1803 |
4/4✓ Branch 0 taken 19178 times.
✓ Branch 1 taken 424 times.
✓ Branch 2 taken 6110 times.
✓ Branch 3 taken 13492 times.
|
38780 | t2 != NO_INDEX && |
| 1804 |
3/4✓ Branch 1 taken 19178 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6110 times.
✓ Branch 4 taken 13068 times.
|
19178 | !Tflag_is_set(t2,T_VISITED_FLAG) |
| 1805 | ) { | ||
| 1806 | // t2 is included in the same operands as t1, | ||
| 1807 | // except for the operands that touch the boundary | ||
| 1808 | // between t1 and t2, for which inclusion changes | ||
| 1809 | 6110 | index_t t2_bits = t1_bits; | |
| 1810 | 6110 | for( | |
| 1811 |
1/2✓ Branch 1 taken 6110 times.
✗ Branch 2 not taken.
|
6110 | index_t ecit = Tedge_cnstr_first(t1,le); |
| 1812 |
2/2✓ Branch 0 taken 2391 times.
✓ Branch 1 taken 6110 times.
|
8501 | ecit != NO_INDEX; |
| 1813 | 2391 | ecit = edge_cnstr_next(ecit) | |
| 1814 | ) { | ||
| 1815 |
1/2✓ Branch 1 taken 2391 times.
✗ Branch 2 not taken.
|
2391 | index_t cnstr = edge_cnstr(ecit); |
| 1816 |
2/4✓ Branch 1 taken 2391 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2391 times.
✗ Branch 5 not taken.
|
2391 | t2_bits ^= cnstr_operand_bits_[cnstr]; |
| 1817 | } | ||
| 1818 |
1/2✓ Branch 1 taken 6110 times.
✗ Branch 2 not taken.
|
6110 | facet_inclusion_bits_[t2] = t2_bits; |
| 1819 |
1/2✓ Branch 1 taken 6110 times.
✗ Branch 2 not taken.
|
6110 | Tset_flag(t2, T_VISITED_FLAG); |
| 1820 |
1/2✓ Branch 1 taken 6110 times.
✗ Branch 2 not taken.
|
6110 | S.push_back(t2); |
| 1821 | } | ||
| 1822 | } | ||
| 1823 | } | ||
| 1824 | |||
| 1825 | // Step 3: reset visited flag | ||
| 1826 |
3/4✓ Branch 1 taken 6640 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6534 times.
✓ Branch 4 taken 106 times.
|
6640 | for(index_t t=0; t<nT(); ++t) { |
| 1827 |
1/2✓ Branch 1 taken 6534 times.
✗ Branch 2 not taken.
|
6534 | Treset_flag(t, T_VISITED_FLAG); |
| 1828 | } | ||
| 1829 | |||
| 1830 | // Step 4: mark triangles to be deleted | ||
| 1831 |
3/4✓ Branch 1 taken 106 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 9 times.
✓ Branch 4 taken 97 times.
|
106 | if(expr == "intersection") { |
| 1832 | 9 | index_t all_bits_set = 0; | |
| 1833 |
2/2✓ Branch 2 taken 222 times.
✓ Branch 3 taken 9 times.
|
240 | for(index_t e_operand_bits: cnstr_operand_bits_) { |
| 1834 | 222 | all_bits_set |= e_operand_bits; | |
| 1835 | } | ||
| 1836 |
3/4✓ Branch 1 taken 507 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 498 times.
✓ Branch 4 taken 9 times.
|
507 | for(index_t t=0; t<nT(); ++t) { |
| 1837 |
3/4✓ Branch 1 taken 498 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 431 times.
✓ Branch 4 taken 67 times.
|
498 | if(facet_inclusion_bits_[t] != all_bits_set) { |
| 1838 |
1/2✓ Branch 1 taken 431 times.
✗ Branch 2 not taken.
|
431 | Tset_flag(t, T_MARKED_FLAG); |
| 1839 | } | ||
| 1840 | } | ||
| 1841 | } else { | ||
| 1842 |
8/14✓ Branch 1 taken 97 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 88 times.
✓ Branch 4 taken 9 times.
✓ Branch 6 taken 88 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 9 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 97 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 88 times.
✓ Branch 16 taken 9 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
|
185 | BooleanExpression E(expr == "union" ? "*" : expr); |
| 1843 |
3/4✓ Branch 1 taken 6133 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6036 times.
✓ Branch 4 taken 97 times.
|
6133 | for(index_t t=0; t<nT(); ++t) { |
| 1844 |
4/6✓ Branch 1 taken 6036 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6036 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 3364 times.
✓ Branch 7 taken 2672 times.
|
6036 | if(!E(facet_inclusion_bits_[t])) { |
| 1845 |
1/2✓ Branch 1 taken 3364 times.
✗ Branch 2 not taken.
|
3364 | Tset_flag(t, T_MARKED_FLAG); |
| 1846 | } | ||
| 1847 | } | ||
| 1848 | 97 | } | |
| 1849 |
1/2✓ Branch 0 taken 106 times.
✗ Branch 1 not taken.
|
106 | if(!mark_only) { |
| 1850 |
1/2✓ Branch 1 taken 106 times.
✗ Branch 2 not taken.
|
106 | remove_marked_triangles(); |
| 1851 | } | ||
| 1852 | 106 | } | |
| 1853 | |||
| 1854 | typedef std::set<index_t> SparseBits; | ||
| 1855 | |||
| 1856 | ✗ | inline void sparse_bits_flip_bit(SparseBits& bits, index_t bit) { | |
| 1857 | ✗ | auto it = bits.find(bit); | |
| 1858 | ✗ | if(it == bits.end()) { | |
| 1859 | ✗ | bits.insert(bit); | |
| 1860 | } else { | ||
| 1861 | ✗ | bits.erase(it); | |
| 1862 | } | ||
| 1863 | ✗ | } | |
| 1864 | |||
| 1865 | ✗ | inline bool sparse_bits_is_zero(const SparseBits& bits) { | |
| 1866 | ✗ | return (bits.size() == 0); | |
| 1867 | } | ||
| 1868 | |||
| 1869 | ✗ | void ExactCDT2d::classify_triangles_union_cnstr_operand_bits_is_operand_id( | |
| 1870 | bool mark_only | ||
| 1871 | ) { | ||
| 1872 | |||
| 1873 | ✗ | DList S(*this, DLIST_S_ID); | |
| 1874 | ✗ | std::stack<SparseBits> Sbits; | |
| 1875 | |||
| 1876 | // Step 1: get triangles adjacent to the border, | ||
| 1877 | // mark them as visited, classify them as to-delete | ||
| 1878 | ✗ | for(index_t t=0; t<nT(); ++t) { | |
| 1879 | ✗ | for(index_t le=0; le<3; ++le) { | |
| 1880 | ✗ | if(Tadj(t,le) == NO_INDEX) { | |
| 1881 | ✗ | Tset_flag(t, T_VISITED_FLAG); | |
| 1882 | ✗ | Tset_flag(t, T_MARKED_FLAG); | |
| 1883 | ✗ | S.push_back(t); | |
| 1884 | ✗ | Sbits.push(SparseBits()); | |
| 1885 | ✗ | break; | |
| 1886 | } | ||
| 1887 | } | ||
| 1888 | } | ||
| 1889 | |||
| 1890 | // Step 2: recursive traversal | ||
| 1891 | ✗ | while(!S.empty()) { | |
| 1892 | ✗ | index_t t1 = S.pop_back(); | |
| 1893 | ✗ | std::set<index_t> t1_bits = Sbits.top(); | |
| 1894 | ✗ | Sbits.pop(); | |
| 1895 | ✗ | for(index_t le=0; le<3; ++le) { | |
| 1896 | ✗ | index_t t2 = Tadj(t1,le); | |
| 1897 | ✗ | if( | |
| 1898 | ✗ | t2 != NO_INDEX && | |
| 1899 | ✗ | !Tflag_is_set(t2,T_VISITED_FLAG) | |
| 1900 | ) { | ||
| 1901 | // t2 is included in the same operands as t1, | ||
| 1902 | // except for the operands that touch the boundary | ||
| 1903 | // between t1 and t2, for which inclusion changes | ||
| 1904 | ✗ | SparseBits t2_bits = t1_bits; | |
| 1905 | ✗ | for( | |
| 1906 | ✗ | index_t ecit = Tedge_cnstr_first(t1,le); | |
| 1907 | ✗ | ecit != NO_INDEX; | |
| 1908 | ✗ | ecit = edge_cnstr_next(ecit) | |
| 1909 | ) { | ||
| 1910 | ✗ | index_t cnstr = edge_cnstr(ecit); | |
| 1911 | ✗ | sparse_bits_flip_bit( | |
| 1912 | ✗ | t2_bits, cnstr_operand_bits_[cnstr] | |
| 1913 | ); | ||
| 1914 | } | ||
| 1915 | ✗ | if(sparse_bits_is_zero(t2_bits)) { | |
| 1916 | ✗ | Tset_flag(t2, T_MARKED_FLAG); | |
| 1917 | } | ||
| 1918 | ✗ | Tset_flag(t2, T_VISITED_FLAG); | |
| 1919 | ✗ | S.push_back(t2); | |
| 1920 | ✗ | Sbits.push(t2_bits); | |
| 1921 | ✗ | } | |
| 1922 | } | ||
| 1923 | ✗ | } | |
| 1924 | |||
| 1925 | // Step 3: reset visited flag | ||
| 1926 | ✗ | for(index_t t=0; t<nT(); ++t) { | |
| 1927 | ✗ | Treset_flag(t, T_VISITED_FLAG); | |
| 1928 | } | ||
| 1929 | |||
| 1930 | ✗ | if(!mark_only) { | |
| 1931 | ✗ | remove_marked_triangles(); | |
| 1932 | } | ||
| 1933 | ✗ | } | |
| 1934 | |||
| 1935 | ✗ | void ExactCDT2d::save(const std::string& filename) const { | |
| 1936 | #ifndef GEOGRAM_PSM | ||
| 1937 | ✗ | Mesh M; | |
| 1938 | ✗ | Attribute<index_t> nb_cnstr(M.edges.attributes(),"nb_cnstr"); | |
| 1939 | ✗ | M.vertices.set_dimension(2); | |
| 1940 | ✗ | for(const ExactPoint& P: point_) { | |
| 1941 | ✗ | double w = P.w.estimate(); | |
| 1942 | ✗ | vec2 p(P.x.estimate() / w, P.y.estimate() / w); | |
| 1943 | ✗ | M.vertices.create_vertex(p.data()); | |
| 1944 | } | ||
| 1945 | ✗ | for(index_t t=0; t<nT(); ++t) { | |
| 1946 | ✗ | index_t i = Tv(t,0); | |
| 1947 | ✗ | index_t j = Tv(t,1); | |
| 1948 | ✗ | index_t k = Tv(t,2); | |
| 1949 | ✗ | M.facets.create_triangle(i,j,k); | |
| 1950 | |||
| 1951 | ✗ | for(index_t le=0; le<3; ++le) { | |
| 1952 | ✗ | if(Tedge_is_constrained(t,le)) { | |
| 1953 | ✗ | index_t e = M.edges.create_edge( | |
| 1954 | ✗ | Tv(t,(le+1)%3), Tv(t,(le+2)%3) | |
| 1955 | ); | ||
| 1956 | ✗ | nb_cnstr[e] = Tedge_cnstr_nb(t,le); | |
| 1957 | } | ||
| 1958 | } | ||
| 1959 | } | ||
| 1960 | ✗ | M.facets.connect(); | |
| 1961 | ✗ | M.vertices.remove_isolated(); | |
| 1962 | ✗ | mesh_save(M, filename); | |
| 1963 | #else | ||
| 1964 | if(!String::string_ends_with(filename,".obj")) { | ||
| 1965 | Logger::err("CDT_2d") | ||
| 1966 | << "save() only supports .obj file format in PSM" | ||
| 1967 | << std::endl; | ||
| 1968 | return; | ||
| 1969 | } | ||
| 1970 | std::ofstream out(filename); | ||
| 1971 | for(const ExactPoint& P: point_) { | ||
| 1972 | double w = P.w.estimate(); | ||
| 1973 | vec2 p(P.x.estimate() / w, P.y.estimate() / w); | ||
| 1974 | out << "v " << p << " " << 0.0 << std::endl; | ||
| 1975 | } | ||
| 1976 | for(index_t t=0; t<nT(); ++t) { | ||
| 1977 | out << "f " << Tv(t,0)+1 << " " << Tv(t,1)+1 << " " << Tv(t,2)+1 | ||
| 1978 | << std::endl; | ||
| 1979 | } | ||
| 1980 | |||
| 1981 | for(index_t t=0; t<nT(); ++t) { | ||
| 1982 | for(index_t le=0; le<3; ++le) { | ||
| 1983 | if(Tedge_is_constrained(t,le)) { | ||
| 1984 | index_t v1 = Tv(t,(le+1)%3); | ||
| 1985 | index_t v2 = Tv(t,(le+2)%3); | ||
| 1986 | out << "l " << v1+1 << " " << v2+1 << std::endl; | ||
| 1987 | } | ||
| 1988 | } | ||
| 1989 | } | ||
| 1990 | #endif | ||
| 1991 | ✗ | } | |
| 1992 | |||
| 1993 | 371635 | double ExactCDT2d::squared_length( | |
| 1994 | const expansion_nt& x, const expansion_nt& y, const expansion_nt& w | ||
| 1995 | ) { | ||
| 1996 | 371635 | double result = 0.0; | |
| 1997 | #if defined(__SIZEOF_FLOAT128__) || defined(__FLOAT128__) | ||
| 1998 | { | ||
| 1999 |
5/10✓ Branch 1 taken 371635 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 371635 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 371635 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 371635 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 371635 times.
✗ Branch 14 not taken.
|
371635 | expansion_nt Num = geo_sqr(x) + geo_sqr(y); |
| 2000 |
2/4✓ Branch 1 taken 371635 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 371635 times.
✗ Branch 5 not taken.
|
371635 | expansion_nt Denom = geo_sqr(w); |
| 2001 | |||
| 2002 | 371635 | __float128 num = 0.0; | |
| 2003 |
2/2✓ Branch 2 taken 3351738 times.
✓ Branch 3 taken 371635 times.
|
3723373 | for(index_t i=0; i<Num.rep().length(); ++i) { |
| 2004 |
1/2✓ Branch 2 taken 3351738 times.
✗ Branch 3 not taken.
|
3351738 | num += Num.rep()[i]; |
| 2005 | } | ||
| 2006 | 371635 | __float128 denom = 0.0; | |
| 2007 |
2/2✓ Branch 2 taken 2016179 times.
✓ Branch 3 taken 371635 times.
|
2387814 | for(index_t i=0; i<Denom.rep().length(); ++i) { |
| 2008 |
1/2✓ Branch 2 taken 2016179 times.
✗ Branch 3 not taken.
|
2016179 | denom += Denom.rep()[i]; |
| 2009 | } | ||
| 2010 | 371635 | result = double(num/denom); | |
| 2011 | 371635 | } | |
| 2012 | #else | ||
| 2013 | result = (geo_sqr(x) + geo_sqr(y)).estimate() / geo_sqr(w).estimate() ; | ||
| 2014 | #endif | ||
| 2015 | 371635 | return result; | |
| 2016 | } | ||
| 2017 | |||
| 2018 | /***************************************************************************/ | ||
| 2019 | } | ||
| 2020 |