| 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 | #include <geogram/delaunay/delaunay_2d.h> | ||
| 41 | #include <geogram/basic/logger.h> | ||
| 42 | #include <geogram/basic/geometry_nd.h> | ||
| 43 | #include <geogram/basic/process.h> | ||
| 44 | #include <geogram/basic/command_line.h> | ||
| 45 | #include <geogram/basic/stopwatch.h> | ||
| 46 | #include <geogram/basic/matrix.h> | ||
| 47 | #include <geogram/basic/permutation.h> | ||
| 48 | #include <geogram/mesh/mesh_reorder.h> | ||
| 49 | #include <geogram/bibliography/bibliography.h> | ||
| 50 | |||
| 51 | #include <stack> | ||
| 52 | #include <algorithm> | ||
| 53 | |||
| 54 | // TODO: optimizations: | ||
| 55 | // - convex hull traversal for nearest_vertex() | ||
| 56 | |||
| 57 | namespace { | ||
| 58 | using namespace GEO; | ||
| 59 | |||
| 60 | /** | ||
| 61 | * \brief Computes the (approximate) orientation predicate in 2d. | ||
| 62 | * \details Computes the sign of the (approximate) signed volume of | ||
| 63 | * the triangle p0, p1, p2 | ||
| 64 | * \param[in] p0 first vertex of the triangle | ||
| 65 | * \param[in] p1 second vertex of the triangle | ||
| 66 | * \param[in] p2 third vertex of the triangle | ||
| 67 | * \retval POSITIVE if the triangle is oriented positively | ||
| 68 | * \retval ZERO if the triangle is flat | ||
| 69 | * \retval NEGATIVE if the triangle is oriented negatively | ||
| 70 | * \todo check whether orientation is inverted as compared to | ||
| 71 | * Shewchuk's version. | ||
| 72 | */ | ||
| 73 | inline Sign orient_2d_inexact( | ||
| 74 | const double* p0, const double* p1, const double* p2 | ||
| 75 | ) { | ||
| 76 | 1799 | double a11 = p1[0] - p0[0] ; | |
| 77 | 1799 | double a12 = p1[1] - p0[1] ; | |
| 78 | |||
| 79 | 1799 | double a21 = p2[0] - p0[0] ; | |
| 80 | 1799 | double a22 = p2[1] - p0[1] ; | |
| 81 | |||
| 82 | double Delta = det2x2( | ||
| 83 | a11,a12, | ||
| 84 | a21,a22 | ||
| 85 | ); | ||
| 86 | |||
| 87 | return geo_sgn(Delta); | ||
| 88 | } | ||
| 89 | } | ||
| 90 | |||
| 91 | namespace GEO { | ||
| 92 | |||
| 93 | // triangle edge vertex is such that the triangle | ||
| 94 | // formed with: | ||
| 95 | // vertex lv | ||
| 96 | // triangle_edge_vertex[lv][0] | ||
| 97 | // triangle_edge_vertex[lv][1] | ||
| 98 | // has the same orientation as the original triangle for | ||
| 99 | // any vertex lv. | ||
| 100 | |||
| 101 | char Delaunay2d::triangle_edge_vertex_[3][2] = { | ||
| 102 | {1,2}, | ||
| 103 | {2,0}, | ||
| 104 | {0,1} | ||
| 105 | }; | ||
| 106 | |||
| 107 | 2 | Delaunay2d::Delaunay2d(coord_index_t dimension) : | |
| 108 |
1/2✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
|
2 | Delaunay(dimension) |
| 109 | { | ||
| 110 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | geo_cite_with_info( |
| 111 | "DBLP:journals/cj/Bowyer81", | ||
| 112 | "One of the two initial references to the algorithm, " | ||
| 113 | "discovered independently and simultaneously by Bowyer and Watson." | ||
| 114 | ); | ||
| 115 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | geo_cite_with_info( |
| 116 | "journals/cj/Watson81", | ||
| 117 | "One of the two initial references to the algorithm, " | ||
| 118 | "discovered independently and simultaneously by Bowyer and Watson." | ||
| 119 | ); | ||
| 120 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | geo_cite_with_info( |
| 121 | "DBLP:conf/compgeom/AmentaCR03", | ||
| 122 | "Using spatial sorting has a dramatic impact on the performances." | ||
| 123 | ); | ||
| 124 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | geo_cite_with_info( |
| 125 | "DBLP:journals/comgeo/FunkeMN05", | ||
| 126 | "Initializing \\verb|locate()| with a non-exact version " | ||
| 127 | " (structural filtering) gains (a bit of) performance." | ||
| 128 | ); | ||
| 129 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | geo_cite_with_info( |
| 130 | "DBLP:journals/comgeo/BoissonnatDPTY02", | ||
| 131 | "The idea of traversing the cavity from inside " | ||
| 132 | " used in GEOGRAM is inspired by the implementation of " | ||
| 133 | " \\verb|Delaunay_triangulation_3| in CGAL." | ||
| 134 | ); | ||
| 135 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | geo_cite_with_info( |
| 136 | "DBLP:conf/imr/Si06", | ||
| 137 | "The triangulation data structure used in GEOGRAM is inspired " | ||
| 138 | "by Tetgen." | ||
| 139 | ); | ||
| 140 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | geo_cite_with_info( |
| 141 | "DBLP:journals/ijfcs/DevillersPT02", | ||
| 142 | "Analysis of the different versions of the line walk algorithm " | ||
| 143 | " used by \\verb|locate()|." | ||
| 144 | ); | ||
| 145 | |||
| 146 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if(dimension != 2 && dimension != 3) { |
| 147 | ✗ | throw InvalidDimension(dimension, "Delaunay2d", "2 or 3"); | |
| 148 | } | ||
| 149 | 2 | first_free_ = END_OF_LIST; | |
| 150 | 2 | weighted_ = (dimension == 3); | |
| 151 | // In weighted mode, vertices are 3d but combinatorics is 2d. | ||
| 152 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if(weighted_) { |
| 153 | ✗ | cell_size_ = 3; | |
| 154 | ✗ | cell_v_stride_ = 3; | |
| 155 | ✗ | cell_neigh_stride_ = 3; | |
| 156 | } | ||
| 157 | 2 | cur_stamp_ = 0; | |
| 158 |
3/6✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
|
2 | debug_mode_ = CmdLine::get_arg_bool("dbg:delaunay"); |
| 159 |
3/6✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
2 | verbose_debug_mode_ = CmdLine::get_arg_bool("dbg:delaunay_verbose"); |
| 160 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | debug_mode_ = (debug_mode_ || verbose_debug_mode_); |
| 161 |
3/6✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
2 | benchmark_mode_ = CmdLine::get_arg_bool("dbg:delaunay_benchmark"); |
| 162 | 2 | has_empty_cells_ = false; | |
| 163 | 2 | abort_if_empty_cell_ = false; | |
| 164 | 2 | } | |
| 165 | |||
| 166 | 8 | Delaunay2d::~Delaunay2d() { | |
| 167 | 8 | } | |
| 168 | |||
| 169 | 2 | void Delaunay2d::set_vertices(index_t nb_vertices, const double* vertices) { | |
| 170 | 2 | has_empty_cells_ = false; | |
| 171 |
1/2✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
|
2 | Stopwatch W("DelInternal", benchmark_mode_); |
| 172 | |||
| 173 | 2 | cur_stamp_ = 0; | |
| 174 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if(weighted_) { |
| 175 | ✗ | heights_.resize(nb_vertices); | |
| 176 | ✗ | for(index_t i = 0; i < nb_vertices; ++i) { | |
| 177 | // Client code uses 3d embedding with ti = sqrt(W - wi) | ||
| 178 | // where W = max(wi) | ||
| 179 | // We recompute the standard "shifted" lifting on | ||
| 180 | // the paraboloid from it. | ||
| 181 | // (we use wi - W, everything is shifted by W, but | ||
| 182 | // we do not care since the power diagram is invariant | ||
| 183 | // by a translation of all weights). | ||
| 184 | ✗ | double w = -geo_sqr(vertices[3 * i + 2]); | |
| 185 | ✗ | heights_[i] = -w + | |
| 186 | ✗ | geo_sqr(vertices[3 * i]) + | |
| 187 | ✗ | geo_sqr(vertices[3 * i + 1]); | |
| 188 | } | ||
| 189 | } | ||
| 190 | |||
| 191 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | Delaunay::set_vertices(nb_vertices, vertices); |
| 192 | |||
| 193 | 2 | index_t expected_triangles = nb_vertices * 2; | |
| 194 | |||
| 195 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | cell_to_v_store_.reserve(expected_triangles * 3); |
| 196 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | cell_to_cell_store_.reserve(expected_triangles * 3); |
| 197 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | cell_next_.reserve(expected_triangles); |
| 198 | |||
| 199 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | cell_to_v_store_.resize(0); |
| 200 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | cell_to_cell_store_.resize(0); |
| 201 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | cell_next_.resize(0); |
| 202 | 2 | first_free_ = END_OF_LIST; | |
| 203 | |||
| 204 | // Sort the vertices spatially. This makes localisation | ||
| 205 | // faster. | ||
| 206 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | if(do_reorder_) { |
| 207 | 2 | compute_BRIO_order( | |
| 208 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | nb_vertices, vertex_ptr(0), reorder_, 2, dimension_ |
| 209 | ); | ||
| 210 | } else { | ||
| 211 | ✗ | reorder_.resize(nb_vertices); | |
| 212 | ✗ | for(index_t i = 0; i < nb_vertices; ++i) { | |
| 213 | ✗ | reorder_[i] = i; | |
| 214 | } | ||
| 215 | } | ||
| 216 | |||
| 217 | double sorting_time = 0; | ||
| 218 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | if(benchmark_mode_) { |
| 219 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | sorting_time = W.elapsed_time(); |
| 220 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
4 | Logger::out("DelInternal1") << "BRIO sorting:" |
| 221 | << sorting_time | ||
| 222 | << std::endl; | ||
| 223 | } | ||
| 224 | |||
| 225 | // The indices of the vertices of the first triangle. | ||
| 226 | index_t v0, v1, v2; | ||
| 227 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
|
2 | if(!create_first_triangle(v0, v1, v2)) { |
| 228 | ✗ | Logger::warn("Delaunay2d") << "All the points are colinear" | |
| 229 | << std::endl; | ||
| 230 | ✗ | return; | |
| 231 | } | ||
| 232 | |||
| 233 | index_t hint = NO_TRIANGLE; | ||
| 234 | // Insert all the vertices incrementally. | ||
| 235 |
2/2✓ Branch 0 taken 229 times.
✓ Branch 1 taken 2 times.
|
231 | for(index_t i = 0; i < nb_vertices; ++i) { |
| 236 | 229 | index_t v = reorder_[i]; | |
| 237 | // Do not re-insert the first four vertices. | ||
| 238 |
6/6✓ Branch 0 taken 227 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 225 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 223 times.
✓ Branch 5 taken 2 times.
|
229 | if(v != v0 && v != v1 && v != v2) { |
| 239 |
1/2✓ Branch 1 taken 223 times.
✗ Branch 2 not taken.
|
223 | index_t new_hint = insert(v, hint); |
| 240 |
2/2✓ Branch 0 taken 96 times.
✓ Branch 1 taken 127 times.
|
223 | if(new_hint == NO_TRIANGLE) { |
| 241 | 96 | has_empty_cells_ = true; | |
| 242 |
1/2✓ Branch 0 taken 96 times.
✗ Branch 1 not taken.
|
96 | if(abort_if_empty_cell_) { |
| 243 | return; | ||
| 244 | } | ||
| 245 | } else { | ||
| 246 | hint = new_hint; | ||
| 247 | } | ||
| 248 | } | ||
| 249 | } | ||
| 250 | |||
| 251 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | if(benchmark_mode_) { |
| 252 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
4 | Logger::out("DelInternal2") << "Core insertion algo:" |
| 253 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
2 | << W.elapsed_time() - sorting_time |
| 254 | << std::endl; | ||
| 255 | } | ||
| 256 | |||
| 257 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if(debug_mode_) { |
| 258 | ✗ | check_combinatorics(verbose_debug_mode_); | |
| 259 | ✗ | check_geometry(verbose_debug_mode_); | |
| 260 | } | ||
| 261 | |||
| 262 | // Compress cell_to_v_store_ and cell_to_cell_store_ | ||
| 263 | // (remove free and virtual tetrahedra). | ||
| 264 | // Since cell_next_ is not used at this point, | ||
| 265 | // we reuse it for storing the conversion array that | ||
| 266 | // maps old trgl indices to new trgl indices | ||
| 267 | // Note: trgl_is_real() uses the previous value of | ||
| 268 | // cell_next(), but we are processing indices | ||
| 269 | // in increasing order and since old2new[t] is always | ||
| 270 | // smaller or equal to t, we never overwrite a value | ||
| 271 | // before needing it. | ||
| 272 | |||
| 273 | vector<index_t>& old2new = cell_next_; | ||
| 274 | index_t nb_triangles = 0; | ||
| 275 | index_t nb_triangles_to_delete = 0; | ||
| 276 | |||
| 277 | { | ||
| 278 |
2/2✓ Branch 0 taken 278 times.
✓ Branch 1 taken 2 times.
|
558 | for(index_t t = 0; t < max_t(); ++t) { |
| 279 | if( | ||
| 280 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 278 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 153 times.
✓ Branch 5 taken 109 times.
|
540 | (keep_infinite_ && !triangle_is_free(t)) || |
| 281 | triangle_is_real(t) | ||
| 282 | ) { | ||
| 283 |
2/2✓ Branch 0 taken 152 times.
✓ Branch 1 taken 1 times.
|
153 | if(t != nb_triangles) { |
| 284 | Memory::copy( | ||
| 285 | 152 | &cell_to_v_store_[nb_triangles * 3], | |
| 286 | 152 | &cell_to_v_store_[t * 3], | |
| 287 | 3 * sizeof(index_t) | ||
| 288 | ); | ||
| 289 | Memory::copy( | ||
| 290 | &cell_to_cell_store_[nb_triangles * 3], | ||
| 291 | &cell_to_cell_store_[t * 3], | ||
| 292 | 3 * sizeof(index_t) | ||
| 293 | ); | ||
| 294 | } | ||
| 295 | 153 | old2new[t] = nb_triangles; | |
| 296 | 153 | ++nb_triangles; | |
| 297 | } else { | ||
| 298 | 125 | old2new[t] = NO_INDEX; | |
| 299 | 125 | ++nb_triangles_to_delete; | |
| 300 | } | ||
| 301 | } | ||
| 302 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | cell_to_v_store_.resize(3 * nb_triangles); |
| 303 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | cell_to_cell_store_.resize(3 * nb_triangles); |
| 304 |
2/2✓ Branch 0 taken 459 times.
✓ Branch 1 taken 2 times.
|
461 | for(index_t i = 0; i < 3 * nb_triangles; ++i) { |
| 305 | 459 | index_t t = cell_to_cell_store_[i]; | |
| 306 | geo_debug_assert(t != NO_INDEX); | ||
| 307 | 459 | t = old2new[t]; | |
| 308 | // Note: t can be equal to -1 when a real trgl is | ||
| 309 | // adjacent to a virtual one (and this is how the | ||
| 310 | // rest of Vorpaline expects to see trgls on the | ||
| 311 | // border). | ||
| 312 | 459 | cell_to_cell_store_[i] = t; | |
| 313 | } | ||
| 314 | } | ||
| 315 | |||
| 316 | // In "keep_infinite" mode, we reorder the cells in such | ||
| 317 | // a way that finite cells have indices [0..nb_finite_cells_-1] | ||
| 318 | // and infinite cells have indices [nb_finite_cells_ .. nb_cells_-1] | ||
| 319 | |||
| 320 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if(keep_infinite_) { |
| 321 | ✗ | nb_finite_cells_ = 0; | |
| 322 | index_t finite_ptr = 0; | ||
| 323 | ✗ | index_t infinite_ptr = nb_triangles - 1; | |
| 324 | for(;;) { | ||
| 325 | ✗ | while(triangle_is_finite(finite_ptr)) { | |
| 326 | ✗ | old2new[finite_ptr] = finite_ptr; | |
| 327 | ✗ | ++finite_ptr; | |
| 328 | ✗ | ++nb_finite_cells_; | |
| 329 | } | ||
| 330 | ✗ | while(!triangle_is_finite(infinite_ptr)) { | |
| 331 | ✗ | old2new[infinite_ptr] = infinite_ptr; | |
| 332 | ✗ | --infinite_ptr; | |
| 333 | } | ||
| 334 | ✗ | if(finite_ptr > infinite_ptr) { | |
| 335 | break; | ||
| 336 | } | ||
| 337 | ✗ | old2new[finite_ptr] = infinite_ptr; | |
| 338 | ✗ | old2new[infinite_ptr] = finite_ptr; | |
| 339 | ✗ | ++nb_finite_cells_; | |
| 340 | ✗ | for(index_t lf=0; lf<3; ++lf) { | |
| 341 | ✗ | std::swap( | |
| 342 | ✗ | cell_to_cell_store_[3*finite_ptr + lf], | |
| 343 | ✗ | cell_to_cell_store_[3*infinite_ptr + lf] | |
| 344 | ); | ||
| 345 | } | ||
| 346 | ✗ | for(index_t lv=0; lv<3; ++lv) { | |
| 347 | ✗ | std::swap( | |
| 348 | ✗ | cell_to_v_store_[3*finite_ptr + lv], | |
| 349 | ✗ | cell_to_v_store_[3*infinite_ptr + lv] | |
| 350 | ); | ||
| 351 | } | ||
| 352 | ✗ | ++finite_ptr; | |
| 353 | ✗ | --infinite_ptr; | |
| 354 | ✗ | } | |
| 355 | ✗ | for(index_t i = 0; i < 3 * nb_triangles; ++i) { | |
| 356 | ✗ | index_t t = cell_to_cell_store_[i]; | |
| 357 | geo_debug_assert(t != NO_INDEX); | ||
| 358 | ✗ | t = old2new[t]; | |
| 359 | geo_debug_assert(t != NO_INDEX); | ||
| 360 | ✗ | cell_to_cell_store_[i] = t; | |
| 361 | } | ||
| 362 | } | ||
| 363 | |||
| 364 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | if(benchmark_mode_) { |
| 365 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if(keep_infinite_) { |
| 366 | ✗ | Logger::out("DelCompress") | |
| 367 | << "Removed " << nb_triangles_to_delete | ||
| 368 | << " triangles (free list)" << std::endl; | ||
| 369 | } else { | ||
| 370 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
4 | Logger::out("DelCompress") |
| 371 | << "Removed " << nb_triangles_to_delete | ||
| 372 | << " triangles (free list and infinite)" << std::endl; | ||
| 373 | } | ||
| 374 | } | ||
| 375 | |||
| 376 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | set_arrays( |
| 377 | nb_triangles, | ||
| 378 | cell_to_v_store_.data(), cell_to_cell_store_.data() | ||
| 379 | ); | ||
| 380 | |||
| 381 | // Not mandatory, but doing so makes it possible to | ||
| 382 | // use locate() in derived classes outside of | ||
| 383 | // set_vertices(). | ||
| 384 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | cell_next_.assign(cell_next_.size(),NO_INDEX); |
| 385 | 2 | } | |
| 386 | |||
| 387 | ✗ | index_t Delaunay2d::nearest_vertex(const double* p) const { | |
| 388 | |||
| 389 | // TODO: For the moment, we fallback to the (unefficient) | ||
| 390 | // baseclass implementation when in weighted mode. | ||
| 391 | ✗ | if(weighted_) { | |
| 392 | ✗ | return Delaunay::nearest_vertex(p); | |
| 393 | } | ||
| 394 | |||
| 395 | // Find a triangle (real or virtual) that contains p | ||
| 396 | ✗ | index_t t = locate(p, NO_TRIANGLE, thread_safe()); | |
| 397 | |||
| 398 | // If p is outside the convex hull of the inserted points, | ||
| 399 | // a special traversal is required (not implemented yet). | ||
| 400 | // TODO: implement convex hull boundary traversal | ||
| 401 | // (for now we fallback to linear search implemented | ||
| 402 | // in baseclass) | ||
| 403 | ✗ | if(t == NO_TRIANGLE || triangle_is_virtual(t)) { | |
| 404 | ✗ | return Delaunay::nearest_vertex(p); | |
| 405 | } | ||
| 406 | |||
| 407 | double sq_dist = 1e30; | ||
| 408 | index_t result = NO_TRIANGLE; | ||
| 409 | |||
| 410 | // Find the nearest vertex among t's vertices | ||
| 411 | ✗ | for(index_t lv = 0; lv < 3; ++lv) { | |
| 412 | index_t v = triangle_vertex(t, lv); | ||
| 413 | // If the triangle is virtual, then the first vertex | ||
| 414 | // is the vertex at infinity and is skipped. | ||
| 415 | ✗ | if(v == NO_INDEX) { | |
| 416 | ✗ | continue; | |
| 417 | } | ||
| 418 | double cur_sq_dist = Geom::distance2(p, vertex_ptr(v), 2); | ||
| 419 | ✗ | if(cur_sq_dist < sq_dist) { | |
| 420 | sq_dist = cur_sq_dist; | ||
| 421 | result = v; | ||
| 422 | } | ||
| 423 | } | ||
| 424 | return result; | ||
| 425 | } | ||
| 426 | |||
| 427 | 223 | index_t Delaunay2d::locate_inexact( | |
| 428 | const double* p, index_t hint, index_t max_iter | ||
| 429 | ) const { | ||
| 430 | |||
| 431 | // If no hint specified, find a triangle randomly | ||
| 432 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 221 times.
|
223 | while(hint == NO_TRIANGLE) { |
| 433 |
2/4✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
|
2 | hint = index_t(Numeric::random_int32()) % max_t(); |
| 434 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if(triangle_is_free(hint)) { |
| 435 | hint = NO_TRIANGLE; | ||
| 436 | } | ||
| 437 | } | ||
| 438 | |||
| 439 | geo_debug_assert(!triangle_is_free(hint)); | ||
| 440 | geo_debug_assert(!triangle_is_in_list(hint)); | ||
| 441 | |||
| 442 | // Always start from a real trgl. If the trgl is virtual, | ||
| 443 | // find its real neighbor (always opposite to the | ||
| 444 | // infinite vertex) | ||
| 445 |
2/2✓ Branch 0 taken 153 times.
✓ Branch 1 taken 70 times.
|
223 | if(triangle_is_virtual(hint)) { |
| 446 |
1/2✓ Branch 0 taken 144 times.
✗ Branch 1 not taken.
|
144 | for(index_t lf = 0; lf < 3; ++lf) { |
| 447 |
2/2✓ Branch 0 taken 70 times.
✓ Branch 1 taken 74 times.
|
144 | if(triangle_vertex(hint, lf) == VERTEX_AT_INFINITY) { |
| 448 | hint = triangle_adjacent(hint, lf); | ||
| 449 | geo_debug_assert(hint != NO_TRIANGLE); | ||
| 450 | 70 | break; | |
| 451 | } | ||
| 452 | } | ||
| 453 | } | ||
| 454 | |||
| 455 | index_t t = hint; | ||
| 456 | index_t t_pred = NO_TRIANGLE; | ||
| 457 | |||
| 458 | 223 | still_walking: | |
| 459 | { | ||
| 460 | const double* pv[3]; | ||
| 461 | 997 | pv[0] = vertex_ptr(finite_triangle_vertex(t,0)); | |
| 462 | 997 | pv[1] = vertex_ptr(finite_triangle_vertex(t,1)); | |
| 463 | 997 | pv[2] = vertex_ptr(finite_triangle_vertex(t,2)); | |
| 464 | |||
| 465 |
2/2✓ Branch 0 taken 2218 times.
✓ Branch 1 taken 188 times.
|
2406 | for(index_t le = 0; le < 3; ++le) { |
| 466 | |||
| 467 | index_t t_next = triangle_adjacent(t,le); | ||
| 468 | |||
| 469 | // If the opposite trgl is -1, then it means that | ||
| 470 | // we are trying to locate() (e.g. called from | ||
| 471 | // nearest_vertex) within a triangulation | ||
| 472 | // from which the infinite trgls were removed. | ||
| 473 |
1/2✓ Branch 0 taken 2218 times.
✗ Branch 1 not taken.
|
2218 | if(t_next == NO_INDEX) { |
| 474 | 35 | return NO_TRIANGLE; | |
| 475 | } | ||
| 476 | |||
| 477 | // If the candidate next triangle is the | ||
| 478 | // one we came from, then we know already that | ||
| 479 | // the orientation is positive, thus we examine | ||
| 480 | // the next candidate (or exit the loop if they | ||
| 481 | // are exhausted). | ||
| 482 |
2/2✓ Branch 0 taken 419 times.
✓ Branch 1 taken 1799 times.
|
2218 | if(t_next == t_pred) { |
| 483 | 419 | continue ; | |
| 484 | } | ||
| 485 | |||
| 486 | // To test the orientation of p w.r.t. the facet f of | ||
| 487 | // t, we replace vertex number f with p in t (same | ||
| 488 | // convention as in CGAL). | ||
| 489 | 1799 | const double* pv_bkp = pv[le]; | |
| 490 | 1799 | pv[le] = p; | |
| 491 |
2/2✓ Branch 0 taken 990 times.
✓ Branch 1 taken 809 times.
|
1799 | Sign ori = orient_2d_inexact(pv[0], pv[1], pv[2]); |
| 492 | |||
| 493 | // If the orientation is not negative, then we cannot | ||
| 494 | // walk towards t_next, and examine the next candidate | ||
| 495 | // (or exit the loop if they are exhausted). | ||
| 496 |
2/2✓ Branch 0 taken 990 times.
✓ Branch 1 taken 809 times.
|
1799 | if(ori != NEGATIVE) { |
| 497 | 990 | pv[le] = pv_bkp; | |
| 498 | 990 | continue; | |
| 499 | } | ||
| 500 | |||
| 501 | // If the opposite trgl is a virtual trgl, then | ||
| 502 | // the point has a positive orientation relative | ||
| 503 | // to the facet on the border of the convex hull, | ||
| 504 | // thus t_next is a trgl in conflict and we are | ||
| 505 | // done. | ||
| 506 |
2/2✓ Branch 0 taken 774 times.
✓ Branch 1 taken 35 times.
|
809 | if(triangle_is_virtual(t_next)) { |
| 507 | return t_next; | ||
| 508 | } | ||
| 509 | |||
| 510 | // If we reach this point, then t_next is a valid | ||
| 511 | // successor, thus we are still walking. | ||
| 512 | t_pred = t; | ||
| 513 | t = t_next; | ||
| 514 |
1/2✓ Branch 0 taken 774 times.
✗ Branch 1 not taken.
|
774 | if(--max_iter != 0) { |
| 515 | 774 | goto still_walking; | |
| 516 | } | ||
| 517 | } | ||
| 518 | } | ||
| 519 | |||
| 520 | // If we reach this point, we did not find a valid successor | ||
| 521 | // for walking (a face for which p has negative orientation), | ||
| 522 | // thus we reached the trgl for which p has all positive | ||
| 523 | // face orientations (i.e. the trgl that contains p). | ||
| 524 | |||
| 525 | 188 | return t; | |
| 526 | } | ||
| 527 | |||
| 528 | |||
| 529 | 223 | index_t Delaunay2d::locate( | |
| 530 | const double* p, index_t hint, bool thread_safe, Sign* orient | ||
| 531 | ) const { | ||
| 532 | |||
| 533 | // Try improving the hint by using the | ||
| 534 | // inexact locate function. This gains | ||
| 535 | // (a little bit) performance (a few | ||
| 536 | // percent in total Delaunay computation | ||
| 537 | // time), but it is better than nothing... | ||
| 538 | // Note: there is a maximum number of trgls | ||
| 539 | // traversed by locate_inexact() (2500) | ||
| 540 | // since there exists configurations in which | ||
| 541 | // locate_inexact() loops forever ! | ||
| 542 | |||
| 543 | 223 | hint = locate_inexact(p, hint, 2500); | |
| 544 | |||
| 545 | static Process::spinlock locate_lock = GEOGRAM_SPINLOCK_INIT; | ||
| 546 | |||
| 547 | // We need to have this spinlock because | ||
| 548 | // of random() that is not thread-safe | ||
| 549 | // (TODO: implement a random() function with | ||
| 550 | // thread local storage) | ||
| 551 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 223 times.
|
223 | if(thread_safe) { |
| 552 | Process::acquire_spinlock(locate_lock); | ||
| 553 | } | ||
| 554 | |||
| 555 | // If no hint specified, find a triangle randomly | ||
| 556 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 223 times.
|
223 | while(hint == NO_TRIANGLE) { |
| 557 | ✗ | hint = index_t(Numeric::random_int32()) % max_t(); | |
| 558 | ✗ | if(triangle_is_free(hint)) { | |
| 559 | hint = NO_TRIANGLE; | ||
| 560 | } | ||
| 561 | } | ||
| 562 | |||
| 563 | geo_debug_assert(!triangle_is_free(hint)); | ||
| 564 | geo_debug_assert(!triangle_is_in_list(hint)); | ||
| 565 | |||
| 566 | // Always start from a real trgl. If the trgl is virtual, | ||
| 567 | // find its real neighbor (always opposite to the | ||
| 568 | // infinite vertex) | ||
| 569 |
2/2✓ Branch 0 taken 188 times.
✓ Branch 1 taken 35 times.
|
223 | if(triangle_is_virtual(hint)) { |
| 570 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 83 times.
|
83 | for(index_t le = 0; le < 3; ++le) { |
| 571 |
2/2✓ Branch 0 taken 35 times.
✓ Branch 1 taken 48 times.
|
83 | if(triangle_vertex(hint, le) == VERTEX_AT_INFINITY) { |
| 572 | hint = triangle_adjacent(hint, le); | ||
| 573 | geo_debug_assert(hint != NO_TRIANGLE); | ||
| 574 | 35 | break; | |
| 575 | } | ||
| 576 | } | ||
| 577 | } | ||
| 578 | |||
| 579 | geo_debug_assert(!triangle_is_free(hint)); | ||
| 580 | geo_debug_assert(!triangle_is_in_list(hint)); | ||
| 581 | geo_debug_assert(!triangle_is_virtual(hint)); | ||
| 582 | |||
| 583 | index_t t = hint; | ||
| 584 | index_t t_pred = NO_TRIANGLE; | ||
| 585 | Sign orient_local[3]; | ||
| 586 |
1/2✓ Branch 0 taken 223 times.
✗ Branch 1 not taken.
|
223 | if(orient == nullptr) { |
| 587 | orient = orient_local; | ||
| 588 | } | ||
| 589 | |||
| 590 | |||
| 591 | 223 | still_walking: | |
| 592 | { | ||
| 593 | const double* pv[3]; | ||
| 594 | 223 | pv[0] = vertex_ptr(finite_triangle_vertex(t,0)); | |
| 595 | 223 | pv[1] = vertex_ptr(finite_triangle_vertex(t,1)); | |
| 596 | 223 | pv[2] = vertex_ptr(finite_triangle_vertex(t,2)); | |
| 597 | |||
| 598 | // Start from a random facet | ||
| 599 | 223 | index_t e0 = index_t(Numeric::random_int32()) % 3; | |
| 600 |
2/2✓ Branch 0 taken 640 times.
✓ Branch 1 taken 188 times.
|
828 | for(index_t de = 0; de < 3; ++de) { |
| 601 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 640 times.
|
640 | index_t le = (e0 + de) % 3; |
| 602 | |||
| 603 | index_t t_next = triangle_adjacent(t,le); | ||
| 604 | |||
| 605 | // If the opposite triangle is -1, then it means that | ||
| 606 | // we are trying to locate() (e.g. called from | ||
| 607 | // nearest_vertex) within a triangulation | ||
| 608 | // from which the infinite trgls were removed. | ||
| 609 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 640 times.
|
640 | if(t_next == NO_INDEX) { |
| 610 | ✗ | if(thread_safe) { | |
| 611 | Process::release_spinlock(locate_lock); | ||
| 612 | } | ||
| 613 | 35 | return NO_TRIANGLE; | |
| 614 | } | ||
| 615 | |||
| 616 | geo_debug_assert(!triangle_is_free(t_next)); | ||
| 617 | geo_debug_assert(!triangle_is_in_list(t_next)); | ||
| 618 | |||
| 619 | // If the candidate next triangle is the | ||
| 620 | // one we came from, then we know already that | ||
| 621 | // the orientation is positive, thus we examine | ||
| 622 | // the next candidate (or exit the loop if they | ||
| 623 | // are exhausted). | ||
| 624 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 640 times.
|
640 | if(t_next == t_pred) { |
| 625 | ✗ | orient[le] = POSITIVE ; | |
| 626 | ✗ | continue ; | |
| 627 | } | ||
| 628 | |||
| 629 | // To test the orientation of p w.r.t. the facet f of | ||
| 630 | // t, we replace vertex number f with p in t (same | ||
| 631 | // convention as in CGAL). | ||
| 632 | // This is equivalent to trgl_facet_point_orient3d(t,f,p) | ||
| 633 | // (but less costly, saves a couple of lookups) | ||
| 634 | 640 | const double* pv_bkp = pv[le]; | |
| 635 | 640 | pv[le] = p; | |
| 636 | 640 | orient[le] = PCK::orient_2d(pv[0], pv[1], pv[2]); | |
| 637 | |||
| 638 | // If the orientation is not negative, then we cannot | ||
| 639 | // walk towards t_next, and examine the next candidate | ||
| 640 | // (or exit the loop if they are exhausted). | ||
| 641 |
2/2✓ Branch 0 taken 605 times.
✓ Branch 1 taken 35 times.
|
640 | if(orient[le] != NEGATIVE) { |
| 642 | 605 | pv[le] = pv_bkp; | |
| 643 | 605 | continue; | |
| 644 | } | ||
| 645 | |||
| 646 | // If the opposite trgl is a virtual trgl, then | ||
| 647 | // the point has a positive orientation relative | ||
| 648 | // to the facet on the border of the convex hull, | ||
| 649 | // thus t_next is a trgl in conflict and we are | ||
| 650 | // done. | ||
| 651 |
1/2✓ Branch 0 taken 35 times.
✗ Branch 1 not taken.
|
35 | if(triangle_is_virtual(t_next)) { |
| 652 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 35 times.
|
35 | if(thread_safe) { |
| 653 | Process::release_spinlock(locate_lock); | ||
| 654 | } | ||
| 655 |
2/2✓ Branch 0 taken 105 times.
✓ Branch 1 taken 35 times.
|
140 | for(index_t tle = 0; tle < 3; ++tle) { |
| 656 | 105 | orient[tle] = POSITIVE; | |
| 657 | } | ||
| 658 | return t_next; | ||
| 659 | } | ||
| 660 | |||
| 661 | // If we reach this point, then t_next is a valid | ||
| 662 | // successor, thus we are still walking. | ||
| 663 | t_pred = t; | ||
| 664 | t = t_next; | ||
| 665 | ✗ | goto still_walking; | |
| 666 | } | ||
| 667 | } | ||
| 668 | |||
| 669 | // If we reach this point, we did not find a valid successor | ||
| 670 | // for walking (a face for which p has negative orientation), | ||
| 671 | // thus we reached the trgl for which p has all positive | ||
| 672 | // face orientations (i.e. the trgl that contains p). | ||
| 673 | |||
| 674 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 188 times.
|
188 | if(thread_safe) { |
| 675 | Process::release_spinlock(locate_lock); | ||
| 676 | } | ||
| 677 | return t; | ||
| 678 | } | ||
| 679 | |||
| 680 | 223 | void Delaunay2d::find_conflict_zone( | |
| 681 | index_t v, | ||
| 682 | index_t t, const Sign* orient, | ||
| 683 | index_t& t_bndry, index_t& e_bndry, | ||
| 684 | index_t& first, index_t& last | ||
| 685 | ) { | ||
| 686 |
2/2✓ Branch 0 taken 127 times.
✓ Branch 1 taken 96 times.
|
223 | first = last = END_OF_LIST; |
| 687 | |||
| 688 | // Generate a unique stamp from current vertex index, | ||
| 689 | // used for marking triangles | ||
| 690 | set_triangle_mark_stamp(v); | ||
| 691 | |||
| 692 | // Pointer to the coordinates of the point to be inserted | ||
| 693 | const double* p = vertex_ptr(v); | ||
| 694 | |||
| 695 | geo_debug_assert(t != NO_TRIANGLE); | ||
| 696 | |||
| 697 | // Test whether the point already exists in | ||
| 698 | // the triangulation. The point already exists | ||
| 699 | // if it's located on three faces of the | ||
| 700 | // triangle returned by locate(). | ||
| 701 | 223 | int nb_zero = | |
| 702 | 223 | (orient[0] == ZERO) + | |
| 703 | 223 | (orient[1] == ZERO) + | |
| 704 | 223 | (orient[2] == ZERO) ; | |
| 705 | |||
| 706 |
2/2✓ Branch 0 taken 127 times.
✓ Branch 1 taken 96 times.
|
223 | if(nb_zero >= 2) { |
| 707 | return; | ||
| 708 | } | ||
| 709 | |||
| 710 | // Weighted triangulations can have dangling | ||
| 711 | // vertices. Such vertices p are characterized by | ||
| 712 | // the fact that p is not in conflict with the | ||
| 713 | // triangle returned by locate(). | ||
| 714 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 127 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
127 | if(weighted_ && !triangle_is_conflict(t, p)) { |
| 715 | return; | ||
| 716 | } | ||
| 717 | |||
| 718 | // Note: points on edges and on facets are | ||
| 719 | // handled by the way triangle_is_in_conflict() | ||
| 720 | // is implemented, that naturally inserts | ||
| 721 | // the correct triangles in the conflict list. | ||
| 722 | |||
| 723 | |||
| 724 | // Mark t as conflict | ||
| 725 | add_triangle_to_list(t, first, last); | ||
| 726 | |||
| 727 | // A small optimization: if the point to be inserted | ||
| 728 | // is on some faces of the located triangle, insert | ||
| 729 | // the neighbors accros those edges in the conflict list. | ||
| 730 | // It saves a couple of calls to the predicates in this | ||
| 731 | // specific case (combinatorics are in general less | ||
| 732 | // expensive than the predicates). | ||
| 733 |
3/4✓ Branch 0 taken 127 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 86 times.
✓ Branch 3 taken 41 times.
|
127 | if(!weighted_ && nb_zero != 0) { |
| 734 |
2/2✓ Branch 0 taken 258 times.
✓ Branch 1 taken 86 times.
|
344 | for(index_t le = 0; le < 3; ++le) { |
| 735 |
2/2✓ Branch 0 taken 86 times.
✓ Branch 1 taken 172 times.
|
258 | if(orient[le] == ZERO) { |
| 736 | index_t t2 = triangle_adjacent(t, le); | ||
| 737 | add_triangle_to_list(t2, first, last); | ||
| 738 | } | ||
| 739 | } | ||
| 740 |
2/2✓ Branch 0 taken 258 times.
✓ Branch 1 taken 86 times.
|
344 | for(index_t le = 0; le < 3; ++le) { |
| 741 |
2/2✓ Branch 0 taken 86 times.
✓ Branch 1 taken 172 times.
|
258 | if(orient[le] == ZERO) { |
| 742 | index_t t2 = triangle_adjacent(t, le); | ||
| 743 | 86 | find_conflict_zone_iterative( | |
| 744 | p,t2,t_bndry,e_bndry,first,last | ||
| 745 | ); | ||
| 746 | } | ||
| 747 | } | ||
| 748 | } | ||
| 749 | |||
| 750 | // Determine the conflict list by greedy propagation from t. | ||
| 751 | 127 | find_conflict_zone_iterative(p,t,t_bndry,e_bndry,first,last); | |
| 752 | } | ||
| 753 | |||
| 754 |
1/2✓ Branch 0 taken 213 times.
✗ Branch 1 not taken.
|
213 | void Delaunay2d::find_conflict_zone_iterative( |
| 755 | const double* p, index_t t_in, | ||
| 756 | index_t& t_bndry, index_t& e_bndry, | ||
| 757 | index_t& first, index_t& last | ||
| 758 | ) { | ||
| 759 | |||
| 760 | S_.push(t_in); | ||
| 761 | |||
| 762 |
2/2✓ Branch 0 taken 468 times.
✓ Branch 1 taken 213 times.
|
681 | while(!S_.empty()) { |
| 763 | |||
| 764 |
1/2✓ Branch 0 taken 468 times.
✗ Branch 1 not taken.
|
468 | index_t t = S_.top(); |
| 765 | S_.pop(); | ||
| 766 | |||
| 767 |
2/2✓ Branch 0 taken 1404 times.
✓ Branch 1 taken 468 times.
|
1872 | for(index_t le = 0; le < 3; ++le) { |
| 768 |
2/2✓ Branch 0 taken 977 times.
✓ Branch 1 taken 427 times.
|
1404 | index_t t2 = triangle_adjacent(t, le); |
| 769 | |||
| 770 | 442 | if( | |
| 771 |
4/4✓ Branch 0 taken 977 times.
✓ Branch 1 taken 427 times.
✓ Branch 2 taken 15 times.
✓ Branch 3 taken 962 times.
|
1404 | triangle_is_in_list(t2) || // known as conflict |
| 772 | triangle_is_marked(t2) // known as non-conflict | ||
| 773 | ) { | ||
| 774 | 697 | continue; | |
| 775 | } | ||
| 776 | |||
| 777 |
2/2✓ Branch 1 taken 255 times.
✓ Branch 2 taken 707 times.
|
962 | if(triangle_is_conflict(t2, p)) { |
| 778 | // Chain t2 in conflict list | ||
| 779 | add_triangle_to_list(t2, first, last); | ||
| 780 | S_.push(t2); | ||
| 781 | 255 | continue; | |
| 782 | } | ||
| 783 | |||
| 784 | // At this point, t is in conflict | ||
| 785 | // and t2 is not in conflict. | ||
| 786 | // We keep a reference to a trgl on the boundary | ||
| 787 | 707 | t_bndry = t; | |
| 788 | 707 | e_bndry = le; | |
| 789 | // Mark t2 as visited (but not conflict) | ||
| 790 | mark_triangle(t2); | ||
| 791 | } | ||
| 792 | } | ||
| 793 | 213 | } | |
| 794 | |||
| 795 | 127 | index_t Delaunay2d::stellate_conflict_zone( | |
| 796 | index_t v_in, index_t t1, index_t t1ebord | ||
| 797 | ) { | ||
| 798 | |||
| 799 | index_t t = t1; | ||
| 800 | index_t e = t1ebord; | ||
| 801 | index_t t_adj = triangle_adjacent(t,e); | ||
| 802 | |||
| 803 | geo_debug_assert(t_adj != NO_INDEX); | ||
| 804 | |||
| 805 | geo_debug_assert(triangle_is_in_list(t)); | ||
| 806 | geo_debug_assert(!triangle_is_in_list(t_adj)); | ||
| 807 | |||
| 808 | |||
| 809 | index_t new_t_first = NO_INDEX; | ||
| 810 | index_t new_t_prev = NO_INDEX; | ||
| 811 | |||
| 812 | do { | ||
| 813 | |||
| 814 | 722 | index_t v1 = triangle_vertex(t, (e+1)%3); | |
| 815 | 722 | index_t v2 = triangle_vertex(t, (e+2)%3); | |
| 816 | |||
| 817 | // Create new triangle | ||
| 818 | 722 | index_t new_t = new_triangle(v_in, v1, v2); | |
| 819 | |||
| 820 | // Connect new triangle to triangle on the other | ||
| 821 | // side of the conflict zone. | ||
| 822 | set_triangle_adjacent(new_t, 0, t_adj); | ||
| 823 | index_t adj_e = find_triangle_adjacent(t_adj, t); | ||
| 824 | set_triangle_adjacent(t_adj, adj_e, new_t); | ||
| 825 | |||
| 826 | |||
| 827 | // Move to next triangle | ||
| 828 | e = (e + 1)%3; | ||
| 829 | t_adj = triangle_adjacent(t,e); | ||
| 830 |
2/2✓ Branch 0 taken 682 times.
✓ Branch 1 taken 722 times.
|
1404 | while(triangle_is_in_list(t_adj)) { |
| 831 | t = t_adj; | ||
| 832 | 682 | e = (find_triangle_vertex(t,v2) + 2)%3; | |
| 833 | t_adj = triangle_adjacent(t,e); | ||
| 834 | geo_debug_assert(t_adj != NO_INDEX); | ||
| 835 | } | ||
| 836 | |||
| 837 |
2/2✓ Branch 0 taken 595 times.
✓ Branch 1 taken 127 times.
|
722 | if(new_t_prev == NO_INDEX) { |
| 838 | new_t_first = new_t; | ||
| 839 | } else { | ||
| 840 | set_triangle_adjacent(new_t_prev, 1, new_t); | ||
| 841 | set_triangle_adjacent(new_t, 2, new_t_prev); | ||
| 842 | } | ||
| 843 | |||
| 844 | new_t_prev = new_t; | ||
| 845 | |||
| 846 |
2/2✓ Branch 0 taken 595 times.
✓ Branch 1 taken 127 times.
|
722 | } while((t != t1) || (e != t1ebord)); |
| 847 | |||
| 848 | // Connect last triangle to first triangle | ||
| 849 | set_triangle_adjacent(new_t_prev, 1, new_t_first); | ||
| 850 | set_triangle_adjacent(new_t_first, 2, new_t_prev); | ||
| 851 | |||
| 852 | 127 | return new_t_prev; | |
| 853 | } | ||
| 854 | |||
| 855 | 223 | index_t Delaunay2d::insert(index_t v, index_t hint) { | |
| 856 | 223 | index_t t_bndry = NO_TRIANGLE; | |
| 857 | 223 | index_t e_bndry = NO_INDEX; | |
| 858 | 223 | index_t first_conflict = NO_TRIANGLE; | |
| 859 | 223 | index_t last_conflict = NO_TRIANGLE; | |
| 860 | |||
| 861 | const double* p = vertex_ptr(v); | ||
| 862 | |||
| 863 | Sign orient[3]; | ||
| 864 | 223 | index_t t = locate(p, hint, false, orient); | |
| 865 | 223 | find_conflict_zone( | |
| 866 | v,t,orient,t_bndry,e_bndry,first_conflict,last_conflict | ||
| 867 | ); | ||
| 868 | |||
| 869 | // The conflict list can be empty if: | ||
| 870 | // - Vertex v already exists in the triangulation | ||
| 871 | // - The triangulation is weighted and v is not visible | ||
| 872 |
2/2✓ Branch 0 taken 127 times.
✓ Branch 1 taken 96 times.
|
223 | if(first_conflict == END_OF_LIST) { |
| 873 | return NO_TRIANGLE; | ||
| 874 | } | ||
| 875 | |||
| 876 | 127 | index_t new_triangle = stellate_conflict_zone(v,t_bndry,e_bndry); | |
| 877 | |||
| 878 | // Recycle the triangles of the conflict zone. | ||
| 879 | 127 | cell_next_[last_conflict] = first_free_; | |
| 880 | 127 | first_free_ = first_conflict; | |
| 881 | |||
| 882 | // Return one of the newly created triangles | ||
| 883 | 127 | return new_triangle; | |
| 884 | } | ||
| 885 | |||
| 886 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | bool Delaunay2d::create_first_triangle( |
| 887 | index_t& iv0, index_t& iv1, index_t& iv2 | ||
| 888 | ) { | ||
| 889 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if(nb_vertices() < 3) { |
| 890 | return false; | ||
| 891 | } | ||
| 892 | |||
| 893 | 2 | iv0 = 0; | |
| 894 | |||
| 895 | 2 | iv1 = 1; | |
| 896 | 2 | while( | |
| 897 |
3/4✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 2 times.
|
6 | iv1 < nb_vertices() && |
| 898 | 3 | PCK::points_are_identical_2d( | |
| 899 | vertex_ptr(iv0), vertex_ptr(iv1) | ||
| 900 | ) | ||
| 901 | ) { | ||
| 902 | 1 | ++iv1; | |
| 903 | } | ||
| 904 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if(iv1 == nb_vertices()) { |
| 905 | return false; | ||
| 906 | } | ||
| 907 | |||
| 908 | 2 | iv2 = iv1 + 1; | |
| 909 | Sign s = ZERO; | ||
| 910 | 2 | while( | |
| 911 |
3/4✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 1 times.
|
6 | iv2 < nb_vertices() && |
| 912 | 3 | (s = PCK::orient_2d( | |
| 913 | vertex_ptr(iv0), vertex_ptr(iv1), vertex_ptr(iv2) | ||
| 914 | )) == ZERO | ||
| 915 | ) { | ||
| 916 | 1 | ++iv2; | |
| 917 | } | ||
| 918 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if(iv2 == nb_vertices()) { |
| 919 | return false; | ||
| 920 | } | ||
| 921 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
2 | if(s == NEGATIVE) { |
| 922 | std::swap(iv1,iv2); | ||
| 923 | } | ||
| 924 | |||
| 925 | // Create the first triangle | ||
| 926 | 2 | index_t t0 = new_triangle(iv0, iv1, iv2); | |
| 927 | |||
| 928 | // Create the first three virtual triangles surrounding it | ||
| 929 | index_t t[3]; | ||
| 930 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
|
8 | for(index_t e = 0; e < 3; ++e) { |
| 931 | // In reverse order since it is an adjacent triangle | ||
| 932 | index_t v1 = triangle_vertex(t0, triangle_edge_vertex(e,1)); | ||
| 933 | index_t v2 = triangle_vertex(t0, triangle_edge_vertex(e,0)); | ||
| 934 | 6 | t[e] = new_triangle(VERTEX_AT_INFINITY, v1, v2); | |
| 935 | } | ||
| 936 | |||
| 937 | // Connect the virtual triangles to the real one | ||
| 938 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
|
8 | for(index_t e=0; e<3; ++e) { |
| 939 | 6 | set_triangle_adjacent(t[e], 0, t0); | |
| 940 | set_triangle_adjacent(t0, e, t[e]); | ||
| 941 | } | ||
| 942 | |||
| 943 | // Interconnect the three virtual triangles along their common | ||
| 944 | // edges | ||
| 945 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
|
8 | for(index_t e = 0; e < 3; ++e) { |
| 946 | // In reverse order since it is an adjacent triangle | ||
| 947 | index_t lv1 = triangle_edge_vertex(e,1); | ||
| 948 | index_t lv2 = triangle_edge_vertex(e,0); | ||
| 949 | 6 | set_triangle_adjacent(t[e], 1, t[lv1]); | |
| 950 | 6 | set_triangle_adjacent(t[e], 2, t[lv2]); | |
| 951 | } | ||
| 952 | |||
| 953 | return true; | ||
| 954 | } | ||
| 955 | |||
| 956 | /************************************************************************/ | ||
| 957 | |||
| 958 | ✗ | void Delaunay2d::show_triangle(index_t t) const { | |
| 959 | std::cerr << "tri" | ||
| 960 | ✗ | << (triangle_is_in_list(t) ? '*' : ' ') | |
| 961 | << t | ||
| 962 | << ", v=[" | ||
| 963 | << triangle_vertex(t, 0) | ||
| 964 | ✗ | << ' ' | |
| 965 | << triangle_vertex(t, 1) | ||
| 966 | ✗ | << ' ' | |
| 967 | << triangle_vertex(t, 2) | ||
| 968 | ✗ | << "] adj=["; | |
| 969 | ✗ | show_triangle_adjacent(t, 0); | |
| 970 | ✗ | show_triangle_adjacent(t, 1); | |
| 971 | ✗ | show_triangle_adjacent(t, 2); | |
| 972 | ✗ | std::cerr << "] "; | |
| 973 | |||
| 974 | ✗ | for(index_t e = 0; e < 3; ++e) { | |
| 975 | ✗ | std::cerr << 'e' << e << ':'; | |
| 976 | ✗ | for(index_t v = 0; v < 2; ++v) { | |
| 977 | std::cerr << triangle_vertex(t, triangle_edge_vertex(e,v)) | ||
| 978 | ✗ | << ','; | |
| 979 | } | ||
| 980 | ✗ | std::cerr << ' '; | |
| 981 | } | ||
| 982 | std::cerr << std::endl; | ||
| 983 | ✗ | } | |
| 984 | |||
| 985 | ✗ | void Delaunay2d::show_triangle_adjacent(index_t t, index_t le) const { | |
| 986 | index_t adj = triangle_adjacent(t, le); | ||
| 987 | ✗ | if(adj != NO_INDEX) { | |
| 988 | ✗ | std::cerr << (triangle_is_in_list(adj) ? '*' : ' '); | |
| 989 | } | ||
| 990 | std::cerr << adj; | ||
| 991 | ✗ | std::cerr << ' '; | |
| 992 | ✗ | } | |
| 993 | |||
| 994 | ✗ | void Delaunay2d::show_list( | |
| 995 | index_t first, const std::string& list_name | ||
| 996 | ) const { | ||
| 997 | index_t t = first; | ||
| 998 | std::cerr << "tri list: " << list_name << std::endl; | ||
| 999 | ✗ | while(t != END_OF_LIST) { | |
| 1000 | ✗ | show_triangle(t); | |
| 1001 | t = triangle_next(t); | ||
| 1002 | } | ||
| 1003 | std::cerr << "-------------" << std::endl; | ||
| 1004 | ✗ | } | |
| 1005 | |||
| 1006 | ✗ | void Delaunay2d::check_combinatorics(bool verbose) const { | |
| 1007 | ✗ | if(verbose) { | |
| 1008 | std::cerr << std::endl; | ||
| 1009 | } | ||
| 1010 | bool ok = true; | ||
| 1011 | ✗ | std::vector<bool> v_has_triangle(nb_vertices(), false); | |
| 1012 | ✗ | for(index_t t = 0; t < max_t(); ++t) { | |
| 1013 | ✗ | if(triangle_is_free(t)) { | |
| 1014 | /* | ||
| 1015 | if(verbose) { | ||
| 1016 | std::cerr << "-Deleted tri: "; | ||
| 1017 | show_tri(t); | ||
| 1018 | } | ||
| 1019 | */ | ||
| 1020 | } else { | ||
| 1021 | /* | ||
| 1022 | if(verbose) { | ||
| 1023 | std::cerr << "Checking tri: "; | ||
| 1024 | show_tri(t); | ||
| 1025 | } | ||
| 1026 | */ | ||
| 1027 | ✗ | for(index_t le = 0; le < 3; ++le) { | |
| 1028 | ✗ | if(triangle_adjacent(t, le) == NO_INDEX) { | |
| 1029 | std::cerr << le << ":Missing adjacent tri" | ||
| 1030 | << std::endl; | ||
| 1031 | ok = false; | ||
| 1032 | ✗ | } else if(triangle_adjacent(t, le) == t) { | |
| 1033 | std::cerr << le << ":Tri is adjacent to itself" | ||
| 1034 | << std::endl; | ||
| 1035 | ok = false; | ||
| 1036 | } else { | ||
| 1037 | index_t t2 = triangle_adjacent(t, le); | ||
| 1038 | bool found = false; | ||
| 1039 | ✗ | for(index_t le2 = 0; le2 < 3; ++le2) { | |
| 1040 | ✗ | if(triangle_adjacent(t2, le2) == t) { | |
| 1041 | found = true; | ||
| 1042 | } | ||
| 1043 | } | ||
| 1044 | ✗ | if(!found) { | |
| 1045 | std::cerr | ||
| 1046 | << le << ":Adjacent link is not bidirectional" | ||
| 1047 | << std::endl; | ||
| 1048 | ok = false; | ||
| 1049 | } | ||
| 1050 | } | ||
| 1051 | } | ||
| 1052 | index_t nb_infinite = 0; | ||
| 1053 | ✗ | for(index_t lv = 0; lv < 3; ++lv) { | |
| 1054 | ✗ | if(triangle_vertex(t, lv) == NO_INDEX) { | |
| 1055 | ✗ | ++nb_infinite; | |
| 1056 | } | ||
| 1057 | } | ||
| 1058 | ✗ | if(nb_infinite > 1) { | |
| 1059 | ok = false; | ||
| 1060 | std::cerr << "More than one infinite vertex" | ||
| 1061 | << std::endl; | ||
| 1062 | } | ||
| 1063 | } | ||
| 1064 | ✗ | for(index_t lv = 0; lv < 3; ++lv) { | |
| 1065 | index_t v = triangle_vertex(t, lv); | ||
| 1066 | ✗ | if(v != NO_INDEX) { | |
| 1067 | ✗ | v_has_triangle[v] = true; | |
| 1068 | } | ||
| 1069 | } | ||
| 1070 | } | ||
| 1071 | ✗ | for(index_t v = 0; v < nb_vertices(); ++v) { | |
| 1072 | ✗ | if(!v_has_triangle[v]) { | |
| 1073 | ✗ | if(verbose) { | |
| 1074 | std::cerr << "Vertex " << v | ||
| 1075 | << " is isolated (duplicated ?)" << std::endl; | ||
| 1076 | } | ||
| 1077 | } | ||
| 1078 | } | ||
| 1079 | ✗ | geo_assert(ok); | |
| 1080 | ✗ | if(verbose) { | |
| 1081 | std::cerr << std::endl; | ||
| 1082 | } | ||
| 1083 | std::cerr << std::endl << "Delaunay Combi OK" << std::endl; | ||
| 1084 | ✗ | } | |
| 1085 | |||
| 1086 | ✗ | void Delaunay2d::check_geometry(bool verbose) const { | |
| 1087 | bool ok = true; | ||
| 1088 | ✗ | for(index_t t = 0; t < max_t(); ++t) { | |
| 1089 | ✗ | if(!triangle_is_free(t)) { | |
| 1090 | index_t v0 = triangle_vertex(t, 0); | ||
| 1091 | index_t v1 = triangle_vertex(t, 1); | ||
| 1092 | index_t v2 = triangle_vertex(t, 2); | ||
| 1093 | ✗ | for(index_t v = 0; v < nb_vertices(); ++v) { | |
| 1094 | ✗ | if(v == v0 || v == v1 || v == v2) { | |
| 1095 | ✗ | continue; | |
| 1096 | } | ||
| 1097 | ✗ | if(triangle_is_conflict(t, vertex_ptr(v))) { | |
| 1098 | ok = false; | ||
| 1099 | ✗ | if(verbose) { | |
| 1100 | std::cerr << "Tri " << t << | ||
| 1101 | " is in conflict with vertex " << v | ||
| 1102 | << std::endl; | ||
| 1103 | |||
| 1104 | ✗ | std::cerr << " offending tri: "; | |
| 1105 | ✗ | show_triangle(t); | |
| 1106 | } | ||
| 1107 | } | ||
| 1108 | } | ||
| 1109 | } | ||
| 1110 | } | ||
| 1111 | ✗ | geo_assert(ok); | |
| 1112 | std::cerr << std::endl << "Delaunay Geo OK" << std::endl; | ||
| 1113 | ✗ | } | |
| 1114 | |||
| 1115 | /************************************************************************/ | ||
| 1116 | |||
| 1117 | ✗ | RegularWeightedDelaunay2d::RegularWeightedDelaunay2d( | |
| 1118 | coord_index_t dimension | ||
| 1119 | ✗ | ) : | |
| 1120 | ✗ | Delaunay2d(3) | |
| 1121 | { | ||
| 1122 | ✗ | if(dimension != 3) { | |
| 1123 | ✗ | throw InvalidDimension(dimension, "RegularWeightedDelaunay2d", "3"); | |
| 1124 | } | ||
| 1125 | ✗ | } | |
| 1126 | |||
| 1127 | ✗ | RegularWeightedDelaunay2d::~RegularWeightedDelaunay2d() { | |
| 1128 | ✗ | } | |
| 1129 | } | ||
| 1130 |