| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /* | ||
| 2 | * Copyright (c) 2000-2022 Inria | ||
| 3 | * All rights reserved. | ||
| 4 | * | ||
| 5 | * Redistribution and use in source and binary forms, with or without | ||
| 6 | * modification, are permitted provided that the following conditions are met: | ||
| 7 | * | ||
| 8 | * * Redistributions of source code must retain the above copyright notice, | ||
| 9 | * this list of conditions and the following disclaimer. | ||
| 10 | * * Redistributions in binary form must reproduce the above copyright notice, | ||
| 11 | * this list of conditions and the following disclaimer in the documentation | ||
| 12 | * and/or other materials provided with the distribution. | ||
| 13 | * * Neither the name of the ALICE Project-Team nor the names of its | ||
| 14 | * contributors may be used to endorse or promote products derived from this | ||
| 15 | * software without specific prior written permission. | ||
| 16 | * | ||
| 17 | * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" | ||
| 18 | * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE | ||
| 19 | * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE | ||
| 20 | * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE | ||
| 21 | * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR | ||
| 22 | * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF | ||
| 23 | * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS | ||
| 24 | * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN | ||
| 25 | * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) | ||
| 26 | * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE | ||
| 27 | * POSSIBILITY OF SUCH DAMAGE. | ||
| 28 | * | ||
| 29 | * Contact: Bruno Levy | ||
| 30 | * | ||
| 31 | * https://www.inria.fr/fr/bruno-levy | ||
| 32 | * | ||
| 33 | * Inria, | ||
| 34 | * Domaine de Voluceau, | ||
| 35 | * 78150 Le Chesnay - Rocquencourt | ||
| 36 | * FRANCE | ||
| 37 | * | ||
| 38 | */ | ||
| 39 | |||
| 40 | #ifndef GEOGRAM_MESH_MESH_SURFACE_INTERSECTION_INTERNAL | ||
| 41 | #define GEOGRAM_MESH_MESH_SURFACE_INTERSECTION_INTERNAL | ||
| 42 | |||
| 43 | /** | ||
| 44 | * \file mesh_surface_intersection_internal.h | ||
| 45 | * \brief Classes used by MeshSurfaceIntersection | ||
| 46 | */ | ||
| 47 | |||
| 48 | #include <geogram/basic/common.h> | ||
| 49 | #include <geogram/mesh/mesh_surface_intersection.h> | ||
| 50 | #include <geogram/mesh/mesh.h> | ||
| 51 | #include <geogram/mesh/mesh_io.h> | ||
| 52 | #include <geogram/mesh/index.h> | ||
| 53 | #include <geogram/mesh/triangle_intersection.h> | ||
| 54 | #include <geogram/delaunay/CDT_2d.h> | ||
| 55 | #include <geogram/numerics/exact_geometry.h> | ||
| 56 | |||
| 57 | #include <map> | ||
| 58 | |||
| 59 | namespace GEO { | ||
| 60 | |||
| 61 | /** | ||
| 62 | * \brief Meshes a single triangle with the constraints that come from | ||
| 63 | * the intersections with the other triangles. | ||
| 64 | * \details Inherits CDTBase2d (constrained Delaunay triangulation), and | ||
| 65 | * redefines orient2d(), incircle2d() and create_intersection() using | ||
| 66 | * vectors with homogeneous coordinates stored as arithmetic expansions | ||
| 67 | * (vec2HE) or arbitrary-precision floating point numbers (vec2HEx) if | ||
| 68 | * compiled with Tessael's geogramplus extension package. | ||
| 69 | */ | ||
| 70 | class MeshInTriangle : public CDTBase2d { | ||
| 71 | public: | ||
| 72 | |||
| 73 | typedef exact::vec3h ExactPoint; | ||
| 74 | |||
| 75 | /***************************************************************/ | ||
| 76 | |||
| 77 | /** | ||
| 78 | * \brief An edge of the mesh. | ||
| 79 | * \details It represents the constraints to be used by the | ||
| 80 | * constrained triangulation to remesh the facet. Makes a maximum | ||
| 81 | * use of the combinatorial information to reduce the complexity | ||
| 82 | * (degree) of the constructed coordinates as much as possible. | ||
| 83 | */ | ||
| 84 | class Edge { | ||
| 85 | public: | ||
| 86 | Edge( | ||
| 87 | index_t v1_in = NO_INDEX, | ||
| 88 | index_t v2_in = NO_INDEX, | ||
| 89 | index_t f2 = NO_INDEX, | ||
| 90 | TriangleRegion R2 = T2_RGN_T | ||
| 91 | 124650 | ) : v1(v1_in), | |
| 92 | 124650 | v2(v2_in) { | |
| 93 | 124650 | sym.f2 = f2; | |
| 94 | 124650 | sym.R2 = R2; | |
| 95 | } | ||
| 96 | index_t v1,v2; // The two extremities of the edge | ||
| 97 | struct { // Symbolic information: this edge = f1 /\ f2.R2 | ||
| 98 | index_t f2; | ||
| 99 | TriangleRegion R2; | ||
| 100 | } sym; | ||
| 101 | }; | ||
| 102 | |||
| 103 | /***************************************************************/ | ||
| 104 | |||
| 105 | /** | ||
| 106 | * \brief A vertex of the triangulation | ||
| 107 | * \details Stores geometric information in exact precision, both | ||
| 108 | * in 3D and in local 2D coordinates. It also stores symbolic | ||
| 109 | * information, that is, facet indices and regions that generated | ||
| 110 | * the vertex. | ||
| 111 | */ | ||
| 112 | 305311 | class Vertex { | |
| 113 | public: | ||
| 114 | |||
| 115 | enum Type { | ||
| 116 | UNINITIALIZED, MESH_VERTEX, PRIMARY_ISECT, SECONDARY_ISECT | ||
| 117 | }; | ||
| 118 | |||
| 119 | /** | ||
| 120 | * \brief Constructor for macro-triangle vertices. | ||
| 121 | * \param[in] f facet index, supposed to correspond to | ||
| 122 | * MeshInTriangle's current facet | ||
| 123 | * \param[in] lv local vertex index in \p f | ||
| 124 | */ | ||
| 125 |
1/2✓ Branch 1 taken 46518 times.
✗ Branch 2 not taken.
|
46518 | Vertex(MeshInTriangle* M, index_t f, index_t lv) { |
| 126 | geo_debug_assert(f == M->f1_); | ||
| 127 | 46518 | type = MESH_VERTEX; | |
| 128 | 46518 | mit = M; | |
| 129 | 46518 | init_sym(f, NO_INDEX, TriangleRegion(lv), T2_RGN_T); | |
| 130 |
1/2✓ Branch 1 taken 46518 times.
✗ Branch 2 not taken.
|
46518 | init_geometry(compute_geometry()); |
| 131 | 46518 | } | |
| 132 | |||
| 133 | /** | ||
| 134 | * \brief Constructor for intersections with other facets. | ||
| 135 | * \param[in] f1 , f2 the two facets. \p f1 is suposed to | ||
| 136 | * correspond to MeshInTriangle's current facet | ||
| 137 | * \param[in] R1 , R2 the two facet regions. | ||
| 138 | */ | ||
| 139 | 240141 | Vertex( | |
| 140 | MeshInTriangle* M, | ||
| 141 | index_t f1, index_t f2, | ||
| 142 | TriangleRegion R1, TriangleRegion R2 | ||
| 143 |
1/2✓ Branch 1 taken 240141 times.
✗ Branch 2 not taken.
|
240141 | ) { |
| 144 | geo_debug_assert(f1 == M->f1_); | ||
| 145 | 240141 | type = PRIMARY_ISECT; | |
| 146 | 240141 | mit = M; | |
| 147 | init_sym(f1,f2,R1,R2); | ||
| 148 |
1/2✓ Branch 1 taken 240141 times.
✗ Branch 2 not taken.
|
240141 | init_geometry(compute_geometry()); |
| 149 | 240141 | } | |
| 150 | |||
| 151 | /** | ||
| 152 | * \brief Constructor for intersections between constraints. | ||
| 153 | * \param[in] point_exact_in exact 3D coordinates | ||
| 154 | * of the intersection | ||
| 155 | */ | ||
| 156 | 5684 | Vertex( | |
| 157 | MeshInTriangle* M, const ExactPoint& point_exact_in | ||
| 158 |
1/2✓ Branch 1 taken 5684 times.
✗ Branch 2 not taken.
|
5684 | ) { |
| 159 | 5684 | type = SECONDARY_ISECT; | |
| 160 | 5684 | mit = M; | |
| 161 | init_sym(NO_INDEX, NO_INDEX, T1_RGN_T, T2_RGN_T); | ||
| 162 | 5684 | init_geometry(point_exact_in); | |
| 163 | 5684 | } | |
| 164 | |||
| 165 | /** | ||
| 166 | * \brief Default constructor | ||
| 167 | */ | ||
| 168 | ✗ | Vertex() { | |
| 169 | ✗ | type = UNINITIALIZED; | |
| 170 | ✗ | mit = nullptr; | |
| 171 | init_sym(NO_INDEX, NO_INDEX, T1_RGN_T, T2_RGN_T); | ||
| 172 | mesh_vertex_index = NO_INDEX; | ||
| 173 | } | ||
| 174 | |||
| 175 | /** | ||
| 176 | * \brief Gets the mesh | ||
| 177 | * \return a reference to the mesh | ||
| 178 | */ | ||
| 179 | const Mesh& mesh() const { | ||
| 180 | 68173 | return mit->mesh(); | |
| 181 | } | ||
| 182 | |||
| 183 | /** | ||
| 184 | * \brief Prints this vertex | ||
| 185 | * \details Displays the combinatorial information | ||
| 186 | * \param[out] out an optional stream where to print | ||
| 187 | */ | ||
| 188 | void print(std::ostream& out=std::cerr) const; | ||
| 189 | |||
| 190 | /** | ||
| 191 | * \brief Gets a string representation of this Vertex | ||
| 192 | * \return a string with the combinatorial information | ||
| 193 | * of this Vertex | ||
| 194 | */ | ||
| 195 | std::string to_string() const { | ||
| 196 | std::ostringstream out; | ||
| 197 | print(out); | ||
| 198 | return out.str(); | ||
| 199 | } | ||
| 200 | |||
| 201 | ✗ | vec2 get_UV_approx() const { | |
| 202 | ✗ | double u = point_exact[mit->u_].estimate(); | |
| 203 | ✗ | double v = point_exact[mit->v_].estimate(); | |
| 204 | double w = point_exact.w.estimate(); | ||
| 205 | ✗ | return vec2(u/w,v/w); | |
| 206 | } | ||
| 207 | |||
| 208 | protected: | ||
| 209 | |||
| 210 | /** | ||
| 211 | * \brief Initializes the symbolic information of this Vertex | ||
| 212 | * \param[in] f1 , f2 the two facets. \p f1 is suposed to | ||
| 213 | * correspond to MeshInTriangle's current facet | ||
| 214 | * \param[in] R1 , R2 the two facet regions. | ||
| 215 | */ | ||
| 216 | void init_sym( | ||
| 217 | index_t f1, index_t f2, TriangleRegion R1, TriangleRegion R2 | ||
| 218 | ) { | ||
| 219 | 292343 | sym.f1 = f1; | |
| 220 | 292343 | sym.f2 = f2; | |
| 221 | 292343 | sym.R1 = R1; | |
| 222 | 292343 | sym.R2 = R2; | |
| 223 |
3/6✓ Branch 1 taken 5684 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 240141 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 46518 times.
✗ Branch 8 not taken.
|
292343 | mesh_vertex_index = NO_INDEX; |
| 224 | } | ||
| 225 | |||
| 226 | /** | ||
| 227 | * \brief Gets the geometry of this vertex | ||
| 228 | * \details Computes the exact 3D position of this vertex | ||
| 229 | * based on the mesh and the combinatorial information | ||
| 230 | */ | ||
| 231 | ExactPoint compute_geometry(); | ||
| 232 | |||
| 233 | /** | ||
| 234 | * \brief Optimizes exact numbers in generated | ||
| 235 | * points and computes approximate coordinates. | ||
| 236 | */ | ||
| 237 | void init_geometry(const ExactPoint& P); | ||
| 238 | |||
| 239 | public: | ||
| 240 | MeshInTriangle* mit; | ||
| 241 | ExactPoint point_exact; // Exact homogeneous coords using expansions | ||
| 242 | Type type; // MESH_VERTEX, PRIMARY_ISECT or SECONDARY_ISECT | ||
| 243 | index_t mesh_vertex_index; // Global mesh vertex index once created | ||
| 244 | struct { // Symbolic information - tri-tri isect | ||
| 245 | index_t f1,f2; // global facet indices in mesh | ||
| 246 | TriangleRegion R1,R2; // triangle regions | ||
| 247 | } sym; | ||
| 248 | #ifndef GEOGRAM_USE_EXACT_NT | ||
| 249 | double l; // precomputed approximated (p[u]^2 + p[v]^2) / p.w^2 | ||
| 250 | #endif | ||
| 251 | }; | ||
| 252 | |||
| 253 | /***************************************************************/ | ||
| 254 | |||
| 255 | MeshInTriangle(MeshSurfaceIntersection& EM); | ||
| 256 | |||
| 257 | /** | ||
| 258 | * \brief Gets the readonly initial mesh | ||
| 259 | * \return a const reference to a copy of the initial mesh | ||
| 260 | */ | ||
| 261 | const Mesh& mesh() const { | ||
| 262 |
2/4✓ Branch 0 taken 46518 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 21655 times.
✗ Branch 3 not taken.
|
68173 | return mesh_; |
| 263 | } | ||
| 264 | |||
| 265 | /** | ||
| 266 | * \brief Gets the target mesh | ||
| 267 | * \return a reference to the target mesh | ||
| 268 | */ | ||
| 269 | Mesh& target_mesh() { | ||
| 270 | 131392 | return exact_mesh_.target_mesh(); | |
| 271 | } | ||
| 272 | |||
| 273 | /** | ||
| 274 | * \brief In dry run mode, the computed local triangulations | ||
| 275 | * are not inserted in the global mesh. This is for benchmarking. | ||
| 276 | * Default is off. | ||
| 277 | */ | ||
| 278 | void set_dry_run(bool x) { | ||
| 279 |
1/2✓ Branch 1 taken 228 times.
✗ Branch 2 not taken.
|
228 | dry_run_ = x; |
| 280 | } | ||
| 281 | |||
| 282 | /** | ||
| 283 | * \brief For debugging, save constraints to a file | ||
| 284 | * \param[in] filename a mesh filename where to solve the constraints | ||
| 285 | * (.obj or .geogram) | ||
| 286 | */ | ||
| 287 | ✗ | void save_constraints(const std::string& filename) { | |
| 288 | ✗ | Mesh M; | |
| 289 | ✗ | get_constraints(M); | |
| 290 | ✗ | mesh_save(M,filename); | |
| 291 | ✗ | } | |
| 292 | |||
| 293 | void begin_facet(index_t f); | ||
| 294 | |||
| 295 | index_t add_vertex(index_t f2, TriangleRegion R1, TriangleRegion R2); | ||
| 296 | |||
| 297 | void add_edge( | ||
| 298 | index_t f2, | ||
| 299 | TriangleRegion AR1, TriangleRegion AR2, | ||
| 300 | TriangleRegion BR1, TriangleRegion BR2 | ||
| 301 | ); | ||
| 302 | |||
| 303 | /** | ||
| 304 | * \brief Creates new vertices and new triangles in target mesh | ||
| 305 | */ | ||
| 306 | void commit(); | ||
| 307 | |||
| 308 | |||
| 309 | /** | ||
| 310 | * \see CDT2d::clear() | ||
| 311 | */ | ||
| 312 | void clear() override; | ||
| 313 | |||
| 314 | protected: | ||
| 315 | /** | ||
| 316 | * \brief For debugging, copies the constraints to a mesh | ||
| 317 | */ | ||
| 318 | void get_constraints(Mesh& M, bool with_edges=true) const; | ||
| 319 | |||
| 320 | vec3 mesh_vertex(index_t v) const { | ||
| 321 | return vec3(mesh().vertices.point_ptr(v)); | ||
| 322 | } | ||
| 323 | |||
| 324 |
1/2✓ Branch 0 taken 1500959 times.
✗ Branch 1 not taken.
|
1500959 | vec3 mesh_facet_vertex(index_t f, index_t lv) const { |
| 325 | index_t v = mesh().facets.vertex(f,lv); | ||
| 326 | 1500959 | return mesh_vertex(v); | |
| 327 | } | ||
| 328 | |||
| 329 | vec2 mesh_vertex_UV(index_t v) const { | ||
| 330 | const double* p = mesh().vertices.point_ptr(v); | ||
| 331 | 621120 | return vec2(p[u_], p[v_]); | |
| 332 | } | ||
| 333 | |||
| 334 |
1/2✓ Branch 0 taken 621120 times.
✗ Branch 1 not taken.
|
621120 | vec2 mesh_facet_vertex_UV(index_t f, index_t lv) const { |
| 335 | index_t v = mesh().facets.vertex(f,lv); | ||
| 336 | 621120 | return mesh_vertex_UV(v); | |
| 337 | } | ||
| 338 | |||
| 339 | |||
| 340 | void log_err() const { | ||
| 341 | std::cerr << "Houston, we got a problem (while remeshing facet " | ||
| 342 | << f1_ << "):" << std::endl; | ||
| 343 | } | ||
| 344 | |||
| 345 | protected: | ||
| 346 | |||
| 347 | /********************** CDTBase2d overrides ***********************/ | ||
| 348 | |||
| 349 | /** | ||
| 350 | * \brief Tests the orientation of three vertices | ||
| 351 | * \param[in] v1 , v2 , v3 the three vertices | ||
| 352 | * \retval POSITIVE if they are in the trigonometric order | ||
| 353 | * \retval ZERO if they are aligned | ||
| 354 | * \retval NEGATIVE if they are in the anti-trigonometric order | ||
| 355 | */ | ||
| 356 | Sign orient2d(index_t v1,index_t v2,index_t v3) const override; | ||
| 357 | |||
| 358 | /** | ||
| 359 | * \brief Tests the relative position of a point with respect | ||
| 360 | * to the circumscribed circle of a triangle | ||
| 361 | * \param[in] v1 , v2 , v3 the three vertices of the triangle | ||
| 362 | * oriented anticlockwise | ||
| 363 | * \param[in] v4 the point to be tested | ||
| 364 | * \retval POSITIVE if the point is inside the circle | ||
| 365 | * \retval ZERO if the point is on the circle | ||
| 366 | * \retval NEGATIVE if the point is outside the circle | ||
| 367 | */ | ||
| 368 | Sign incircle( | ||
| 369 | index_t v1,index_t v2,index_t v3,index_t v4 | ||
| 370 | ) const override; | ||
| 371 | |||
| 372 | /** | ||
| 373 | * \brief Given two segments that have an intersection, create the | ||
| 374 | * intersection | ||
| 375 | * \details The intersection is given both as the indices of segment | ||
| 376 | * extremities (i,j) and (k,l), that one can use to retreive the | ||
| 377 | * points in derived classes, and constraint indices E1 and E2, that | ||
| 378 | * derived classes may use to retreive symbolic information attached | ||
| 379 | * to the constraint | ||
| 380 | * \param[in] e1 the index of the first edge, corresponding to the | ||
| 381 | * value of ncnstr() when insert_constraint() was called for | ||
| 382 | * that edge | ||
| 383 | * \param[in] i , j the vertices of the first segment | ||
| 384 | * \param[in] e2 the index of the second edge, corresponding to the | ||
| 385 | * value of ncnstr() when insert_constraint() was called for | ||
| 386 | * that edge | ||
| 387 | * \param[in] k , l the vertices of the second segment | ||
| 388 | * \return the index of a newly created vertex that corresponds to | ||
| 389 | * the intersection between [\p i , \p j] and [\p k , \p l] | ||
| 390 | */ | ||
| 391 | index_t create_intersection( | ||
| 392 | index_t e1, index_t i, index_t j, | ||
| 393 | index_t e2, index_t k, index_t l | ||
| 394 | ) override; | ||
| 395 | |||
| 396 | /** | ||
| 397 | * \brief Computes the intersection between two edges | ||
| 398 | * \param[in] e1 , e2 the two edges | ||
| 399 | * \param[out] I the intersection | ||
| 400 | */ | ||
| 401 | void get_edge_edge_intersection( | ||
| 402 | index_t e1, index_t e2, ExactPoint& I | ||
| 403 | ) const; | ||
| 404 | |||
| 405 | /** | ||
| 406 | * \brief Auxilliary function used by get_edge_edge_intersection() | ||
| 407 | * for the special case when the two edges are coplanar | ||
| 408 | * \param[in] e1 , e2 the two edges | ||
| 409 | * \param[out] I the intersection | ||
| 410 | */ | ||
| 411 | void get_edge_edge_intersection_2D( | ||
| 412 | index_t e1, index_t e2, ExactPoint& I | ||
| 413 | ) const; | ||
| 414 | |||
| 415 | public: | ||
| 416 | void save(const std::string& filename) const override; | ||
| 417 | |||
| 418 | protected: | ||
| 419 | void begin_insert_transaction() override; | ||
| 420 | void commit_insert_transaction() override; | ||
| 421 | void rollback_insert_transaction() override; | ||
| 422 | |||
| 423 | private: | ||
| 424 | MeshSurfaceIntersection& exact_mesh_; | ||
| 425 | const Mesh& mesh_; | ||
| 426 | index_t f1_; | ||
| 427 | coord_index_t f1_normal_axis_; | ||
| 428 | coord_index_t u_; // = (f1_normal_axis_ + 1)%3 | ||
| 429 | coord_index_t v_; // = (f1_normal_axis_ + 2)%3 | ||
| 430 | vector<Vertex> vertex_; | ||
| 431 | vector<Edge> edges_; | ||
| 432 | bool dry_run_; | ||
| 433 | mutable std::map<trindex, Sign> pred_cache_; | ||
| 434 | bool use_pred_cache_insert_buffer_; | ||
| 435 | mutable std::vector< std::pair<trindex, Sign> > | ||
| 436 | pred_cache_insert_buffer_; | ||
| 437 | }; | ||
| 438 | |||
| 439 | /*************************************************************************/ | ||
| 440 | |||
| 441 | /** | ||
| 442 | * \brief Stores information about a triangle-triangle intersection. | ||
| 443 | * \details The intersection is a segment A-B. Its extremities A and B | ||
| 444 | * are indicated by the regions in f1 and f2 that created the | ||
| 445 | * intersection. If the intersection is just a point, | ||
| 446 | * then A and B regions are the same. | ||
| 447 | */ | ||
| 448 | struct IsectInfo { | ||
| 449 | public: | ||
| 450 | |||
| 451 | /** | ||
| 452 | * Swaps the two facets and updates the combinatorial | ||
| 453 | * information accordingly. | ||
| 454 | */ | ||
| 455 | 67883 | void flip() { | |
| 456 | std::swap(f1,f2); | ||
| 457 | 67883 | A_rgn_f1 = swap_T1_T2(A_rgn_f1); | |
| 458 | 67883 | A_rgn_f2 = swap_T1_T2(A_rgn_f2); | |
| 459 | std::swap(A_rgn_f1, A_rgn_f2); | ||
| 460 | 67883 | B_rgn_f1 = swap_T1_T2(B_rgn_f1); | |
| 461 | 67883 | B_rgn_f2 = swap_T1_T2(B_rgn_f2); | |
| 462 | std::swap(B_rgn_f1, B_rgn_f2); | ||
| 463 | 67883 | } | |
| 464 | |||
| 465 | /** | ||
| 466 | * \brief Tests whether intersection is just a point. | ||
| 467 | * \details Points are encoded as segments with the | ||
| 468 | * same symbolic information for both vertices. | ||
| 469 | */ | ||
| 470 | bool is_point() const { | ||
| 471 | return | ||
| 472 |
2/2✓ Branch 0 taken 54379 times.
✓ Branch 1 taken 81387 times.
|
135766 | A_rgn_f1 == B_rgn_f1 && |
| 473 |
2/2✓ Branch 0 taken 8606 times.
✓ Branch 1 taken 45773 times.
|
54379 | A_rgn_f2 == B_rgn_f2 ; |
| 474 | } | ||
| 475 | |||
| 476 | index_t f1; | ||
| 477 | index_t f2; | ||
| 478 | TriangleRegion A_rgn_f1; | ||
| 479 | TriangleRegion A_rgn_f2; | ||
| 480 | TriangleRegion B_rgn_f1; | ||
| 481 | TriangleRegion B_rgn_f2; | ||
| 482 | }; | ||
| 483 | |||
| 484 | /**********************************************************************/ | ||
| 485 | |||
| 486 | /** | ||
| 487 | * \brief Detects and retriangulates a set of coplanar facets for | ||
| 488 | * MeshSurfaceIntersection. | ||
| 489 | */ | ||
| 490 | class CoplanarFacets { | ||
| 491 | public: | ||
| 492 | static constexpr index_t NON_MANIFOLD = index_t(-2); | ||
| 493 | typedef MeshSurfaceIntersection::ExactPoint ExactPoint; | ||
| 494 | |||
| 495 | /** | ||
| 496 | * \brief Constructs a CoplanarFacets object associated with a | ||
| 497 | * MeshSurfaceIntersection | ||
| 498 | * \details No set of facets is identified. One needs to call get(). | ||
| 499 | * \param[in] I a reference to the MeshSurfaceIntersection | ||
| 500 | * \param[in] clear_attributes if set, resets facet_chart and | ||
| 501 | * keep_vertex | ||
| 502 | * \param[in] angle_tolerance angle tolerance for detecting coplanar | ||
| 503 | * facets and colinear edges (in degrees) | ||
| 504 | */ | ||
| 505 | CoplanarFacets( | ||
| 506 | MeshSurfaceIntersection& I, bool clear_attributes, | ||
| 507 | double angle_tolerance = 0.0 | ||
| 508 | ); | ||
| 509 | |||
| 510 | /** | ||
| 511 | * \brief Gets the set of coplanar facets from a given facet and | ||
| 512 | * group id. | ||
| 513 | * \details Uses the "group" facets attribute. If \p f's group is | ||
| 514 | * uninitialized (NO_INDEX), determines the facets of the group | ||
| 515 | * geometrically and initializes the attribute, else gets the | ||
| 516 | * facets based on the attribute. | ||
| 517 | * \param[in] f the facet | ||
| 518 | * \param[in] group_id the facet group id | ||
| 519 | */ | ||
| 520 | void get(index_t f, index_t group_id); | ||
| 521 | |||
| 522 | /** | ||
| 523 | * \brief Marks the vertices that need to be kept in the | ||
| 524 | * simplified facets. | ||
| 525 | * \details A vertex is kept if it is incident to at least | ||
| 526 | * two non-colinear | ||
| 527 | * edges on the border. The status of the vertices is stored in the | ||
| 528 | * "keep" vertex attribute. | ||
| 529 | */ | ||
| 530 | void mark_vertices_to_keep(); | ||
| 531 | |||
| 532 | /** | ||
| 533 | * \brief For debugging purposes, saves border edges to a file. | ||
| 534 | * \param[in] filename the file where to store the borders. | ||
| 535 | */ | ||
| 536 | void save_borders(const std::string& filename); | ||
| 537 | |||
| 538 | /** | ||
| 539 | * \brief For debugging purposes, saves all the facets of the group | ||
| 540 | * to a file. | ||
| 541 | * \param[in] filename the file where to store the facets of the group. | ||
| 542 | */ | ||
| 543 | void save_facet_group(const std::string& filename); | ||
| 544 | |||
| 545 | /** | ||
| 546 | * \brief Triangulates the kept vertices. | ||
| 547 | * \details One can get the triangle through the (public) CDT member | ||
| 548 | * (ExactCDT2d). | ||
| 549 | */ | ||
| 550 | void triangulate(); | ||
| 551 | |||
| 552 | protected: | ||
| 553 | |||
| 554 | /** | ||
| 555 | * \brief Finds all the pairs of coplanar facets | ||
| 556 | * \details Initializes c_is_coplanar_[], a vector of booleans indexed | ||
| 557 | * by facet corners. | ||
| 558 | */ | ||
| 559 | void find_coplanar_facets(); | ||
| 560 | |||
| 561 | /** | ||
| 562 | * \brief Tests whether two triangles are coplanar | ||
| 563 | * \details This is used to determine the facets that can be | ||
| 564 | * merged | ||
| 565 | * \param[in] p1 , p2 , p3 the vertices of the first triangle | ||
| 566 | * \param[in] q1 , q2 , q3 the vertices of the second triangle | ||
| 567 | * \retval true if the two triangles are coplanar | ||
| 568 | * \retval false otherwise | ||
| 569 | * \details uses angle_tolerance specified to the constructor (if set | ||
| 570 | * to zero, uses exact computation) | ||
| 571 | */ | ||
| 572 | bool triangles_are_coplanar( | ||
| 573 | const vec3& p1, const vec3& p2, const vec3& p3, | ||
| 574 | const vec3& q1, const vec3& q2, const vec3& q3 | ||
| 575 | ) const; | ||
| 576 | |||
| 577 | |||
| 578 | /** | ||
| 579 | * \brief Tests whether two edges are co-linear | ||
| 580 | * \param[in] P1 , P2 , P3 the vertices of the two edges | ||
| 581 | * \retval true if [P1,P2] and [P2,P3] are co-linear, and P2 is between | ||
| 582 | * P1 and p3 | ||
| 583 | * \retval false otherwise | ||
| 584 | * \details uses angle_tolerance specified to the constructor (if set | ||
| 585 | * to zero, uses exact computation) | ||
| 586 | */ | ||
| 587 | bool edges_are_colinear( | ||
| 588 | const ExactPoint& P1, const ExactPoint& P2, const ExactPoint& P3 | ||
| 589 | ) const; | ||
| 590 | |||
| 591 | |||
| 592 | |||
| 593 | public: | ||
| 594 | ExactCDT2d CDT; | ||
| 595 | |||
| 596 | /** | ||
| 597 | * \brief Gets the number of coplanar facets | ||
| 598 | * \return the number of coplanar facets present in the mesh | ||
| 599 | */ | ||
| 600 | index_t nb_facets() { | ||
| 601 | return facets_.size(); | ||
| 602 | } | ||
| 603 | |||
| 604 | /** | ||
| 605 | * \brief Marks the facets | ||
| 606 | * \param[out] facet_is_marked on exit, set to 1 for facets present | ||
| 607 | * in the list of coplanar facets. Needs to be of size | ||
| 608 | * mesh_.facets.nb(). | ||
| 609 | */ | ||
| 610 | void mark_facets(vector<index_t>& facet_is_marked) { | ||
| 611 |
2/2✓ Branch 0 taken 58475 times.
✓ Branch 1 taken 8069 times.
|
66544 | for(index_t f: facets_) { |
| 612 | 58475 | facet_is_marked[f] = 1; | |
| 613 | } | ||
| 614 | } | ||
| 615 | |||
| 616 | private: | ||
| 617 | MeshSurfaceIntersection& I_; | ||
| 618 | Mesh& mesh_; | ||
| 619 | const Mesh& mesh_copy_; | ||
| 620 | double angle_tolerance_; | ||
| 621 | index_t group_id_; | ||
| 622 | Attribute<index_t> facet_group_; | ||
| 623 | Attribute<bool> keep_vertex_; | ||
| 624 | Attribute<bool> c_is_coplanar_; | ||
| 625 | Attribute<bool> f_is_flipped_; | ||
| 626 | vector<bool> f_visited_; | ||
| 627 | vector<bool> h_visited_; | ||
| 628 | vector<bool> v_visited_; | ||
| 629 | vector<index_t> v_idx_; | ||
| 630 | coord_index_t u_; | ||
| 631 | coord_index_t v_; | ||
| 632 | |||
| 633 | /***********************************************************/ | ||
| 634 | |||
| 635 | vector<index_t> vertices_; | ||
| 636 | vector<index_t> facets_; | ||
| 637 | |||
| 638 | /** | ||
| 639 | * \brief A 2d Incident Edge Lists data structure | ||
| 640 | */ | ||
| 641 | class Halfedges { | ||
| 642 | public: | ||
| 643 | |||
| 644 | /** | ||
| 645 | * \brief Halfedges constructor | ||
| 646 | * \param[in] coplanar_facets a reference to the CoplanarFacets | ||
| 647 | */ | ||
| 648 | Halfedges( | ||
| 649 | CoplanarFacets& coplanar_facets | ||
| 650 | 285 | ) : mesh_(coplanar_facets.mesh_) { | |
| 651 | } | ||
| 652 | |||
| 653 | /** | ||
| 654 | * \brief Initializes this Halfedges | ||
| 655 | * \details Clears the list of halfedges and incident edge lists | ||
| 656 | */ | ||
| 657 | 146176 | void initialize() { | |
| 658 | // Use resize rather than assign so that we do not traverse | ||
| 659 | // all the halfedges of the mesh_ | ||
| 660 | 146176 | v_first_halfedge_.resize(mesh_.vertices.nb(), NO_INDEX); | |
| 661 | 146176 | h_next_around_v_.resize(mesh_.facet_corners.nb(), NO_INDEX); | |
| 662 | // We only need to reset the halfedges of this set of coplanar | ||
| 663 | // facets. | ||
| 664 |
2/2✓ Branch 0 taken 491821 times.
✓ Branch 1 taken 146176 times.
|
637997 | for(index_t h: halfedges_) { |
| 665 | 491821 | v_first_halfedge_[vertex(h,0)] = NO_INDEX; | |
| 666 | 491821 | h_next_around_v_[h] = NO_INDEX; | |
| 667 | } | ||
| 668 | 146176 | halfedges_.resize(0); | |
| 669 | 146176 | } | |
| 670 | |||
| 671 | /** | ||
| 672 | * \brief Gets a vertex of a halfedge | ||
| 673 | * \param[in] h the halfedge | ||
| 674 | * \param[in] dlv 0 for origin, 1 for destination, 2 for opposite | ||
| 675 | * \return the vertex | ||
| 676 | */ | ||
| 677 | index_t vertex(index_t h, index_t dlv) const { | ||
| 678 | 2168759 | index_t f = h/3; | |
| 679 |
1/4✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 4 taken 246922 times.
✗ Branch 5 not taken.
|
246922 | index_t lv = (h+dlv)%3; |
| 680 |
3/8✓ Branch 0 taken 51865 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 246922 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 491821 times.
✗ Branch 7 not taken.
|
790608 | return mesh_.facets.vertex(f,lv); |
| 681 | } | ||
| 682 | |||
| 683 | /** | ||
| 684 | * \brief Gets the incident facet | ||
| 685 | * \param[in] h a halfedge | ||
| 686 | * \return the facet incident to the halfedge | ||
| 687 | */ | ||
| 688 | index_t facet(index_t h) const { | ||
| 689 | return h/3; | ||
| 690 | } | ||
| 691 | |||
| 692 | /** | ||
| 693 | * \brief Gets the opposite halfedge | ||
| 694 | * \param[in] h a halfedge | ||
| 695 | * \return the halfedge opposite to \p h, or NO_INDEX if there | ||
| 696 | * is no such halfedge | ||
| 697 | */ | ||
| 698 | index_t alpha2(index_t h) const { | ||
| 699 | index_t t1 = facet(h); | ||
| 700 | index_t t2 = mesh_.facet_corners.adjacent_facet(h); | ||
| 701 | if(t2 == NO_INDEX) { | ||
| 702 | return NO_INDEX; | ||
| 703 | } | ||
| 704 | for(index_t h2: mesh_.facets.corners(t2)) { | ||
| 705 | if(mesh_.facet_corners.adjacent_facet(h2) == t1) { | ||
| 706 | return h2; | ||
| 707 | } | ||
| 708 | } | ||
| 709 | geo_assert_not_reached; | ||
| 710 | } | ||
| 711 | |||
| 712 | /** | ||
| 713 | * \brief Adds a halfedge | ||
| 714 | * \param[in] h the halfedge | ||
| 715 | */ | ||
| 716 | 493844 | void add(index_t h) { | |
| 717 |
2/2✓ Branch 0 taken 492401 times.
✓ Branch 1 taken 1443 times.
|
493844 | halfedges_.push_back(h); |
| 718 |
1/2✓ Branch 0 taken 493844 times.
✗ Branch 1 not taken.
|
493844 | index_t v1 = vertex(h,0); |
| 719 | 493844 | h_next_around_v_[h] = v_first_halfedge_[v1]; | |
| 720 | 493844 | v_first_halfedge_[v1] = h; | |
| 721 | 493844 | } | |
| 722 | |||
| 723 | /** | ||
| 724 | * \brief used by range-based for | ||
| 725 | * \return an iterator to the first halfedge | ||
| 726 | */ | ||
| 727 | vector<index_t>::const_iterator begin() const { | ||
| 728 | return halfedges_.begin(); | ||
| 729 | } | ||
| 730 | |||
| 731 | /** | ||
| 732 | * \brief used by range-based for | ||
| 733 | * \return an iterator to one position past the last halfedge | ||
| 734 | */ | ||
| 735 | vector<index_t>::const_iterator end() const { | ||
| 736 | return halfedges_.end(); | ||
| 737 | } | ||
| 738 | |||
| 739 | /** | ||
| 740 | * \brief Gets the first halfedge starting from a vertex | ||
| 741 | * \param[in] v the vertex | ||
| 742 | * \return the first halfedge starting from \p v | ||
| 743 | */ | ||
| 744 | index_t vertex_first_halfedge(index_t v) const { | ||
| 745 | 987686 | return v_first_halfedge_[v]; | |
| 746 | } | ||
| 747 | |||
| 748 | /** | ||
| 749 | * \brief Gets the next halfedge in the incident edge list | ||
| 750 | * \param[in] h a halfedge | ||
| 751 | * \return the next halfedge around the origin of \p h or | ||
| 752 | * NO_INDEX if there is no such halfedge | ||
| 753 | */ | ||
| 754 | index_t next_around_vertex(index_t h) const { | ||
| 755 | 987696 | return h_next_around_v_[h]; | |
| 756 | } | ||
| 757 | |||
| 758 | /** | ||
| 759 | * \brief Gets the number of halfedges around a vertex | ||
| 760 | * \partam[in] v a vertex | ||
| 761 | * \return the number of halfedges starting from \p v | ||
| 762 | */ | ||
| 763 | index_t nb_halfedges_around_vertex(index_t v) const { | ||
| 764 | index_t result = 0; | ||
| 765 | 987686 | for( | |
| 766 | index_t h = vertex_first_halfedge(v); | ||
| 767 |
4/4✓ Branch 0 taken 493844 times.
✓ Branch 1 taken 493842 times.
✓ Branch 2 taken 493848 times.
✓ Branch 3 taken 493844 times.
|
1975378 | h != NO_INDEX; |
| 768 | h = next_around_vertex(h) | ||
| 769 | ) { | ||
| 770 | 987692 | ++result; | |
| 771 | } | ||
| 772 | return result; | ||
| 773 | } | ||
| 774 | |||
| 775 | /** | ||
| 776 | * \brief Gets the next halfedge along a polyline | ||
| 777 | * \param[in] h a halfedge | ||
| 778 | * \return the halfedge on the same polyline as \p h starting | ||
| 779 | * from \p h destination or NO_INDEX if there is no such | ||
| 780 | * halfedge. Polyline stops where it encounters a vertex that does | ||
| 781 | * not have exactly 1 incident halfedge, that is, | ||
| 782 | * where the halfedges graph is non-manifold. | ||
| 783 | */ | ||
| 784 |
1/2✓ Branch 0 taken 493844 times.
✗ Branch 1 not taken.
|
493844 | index_t next_along_polyline(index_t h) const { |
| 785 | index_t v2 = vertex(h,1); | ||
| 786 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 493840 times.
|
493844 | if(nb_halfedges_around_vertex(v2) != 1) { |
| 787 | 4 | return NO_INDEX; | |
| 788 | } | ||
| 789 | return vertex_first_halfedge(v2); | ||
| 790 | } | ||
| 791 | |||
| 792 | private: | ||
| 793 | Mesh& mesh_; | ||
| 794 | vector<index_t> halfedges_; | ||
| 795 | vector<index_t> v_first_halfedge_; | ||
| 796 | vector<index_t> h_next_around_v_; | ||
| 797 | } halfedges_; | ||
| 798 | |||
| 799 | |||
| 800 | /** | ||
| 801 | * \brief Organizes halfedges as a set of chains starting and ending | ||
| 802 | * at non-manifold vertices | ||
| 803 | */ | ||
| 804 | class Polylines { | ||
| 805 | public: | ||
| 806 | |||
| 807 | /** | ||
| 808 | * \brief Polylines constructor | ||
| 809 | * \param[in] CF a reference to the CoplanarFacets | ||
| 810 | */ | ||
| 811 |
2/2✓ Branch 0 taken 57 times.
✓ Branch 1 taken 228 times.
|
285 | Polylines(CoplanarFacets& CF) : CF_(CF) { |
| 812 | } | ||
| 813 | |||
| 814 | /** | ||
| 815 | * \brief Initializes this Polylines | ||
| 816 | * \details Resets all the stored polylines | ||
| 817 | */ | ||
| 818 | 146176 | void initialize() { | |
| 819 | 146176 | H_.resize(0); | |
| 820 | 146176 | polyline_start_.resize(0); | |
| 821 | 146176 | polyline_start_.push_back(0); | |
| 822 | 146176 | } | |
| 823 | |||
| 824 | /** | ||
| 825 | * \brief Gets the number of polylines | ||
| 826 | * \return the number of polylines | ||
| 827 | */ | ||
| 828 | index_t nb() const { | ||
| 829 | 81157 | return polyline_start_.size() - 1; | |
| 830 | } | ||
| 831 | |||
| 832 | /** | ||
| 833 | * \brief used by range-based for | ||
| 834 | * \return a non-iterator corresponding to the first polyline index. | ||
| 835 | */ | ||
| 836 | index_as_iterator begin() const { | ||
| 837 | return index_as_iterator(0); | ||
| 838 | } | ||
| 839 | |||
| 840 | /** | ||
| 841 | * \brief used by range-based for | ||
| 842 | * \return a non-iterator to one position past the | ||
| 843 | * last polyline index. | ||
| 844 | */ | ||
| 845 | index_as_iterator end() const { | ||
| 846 | 81157 | return index_as_iterator(nb()); | |
| 847 | } | ||
| 848 | |||
| 849 | /** | ||
| 850 | * \brief Gets the halfedges in a polyline | ||
| 851 | * \param[in] polyline the polyline index | ||
| 852 | * \return an iteratable sequence of halfedges | ||
| 853 | */ | ||
| 854 | const_index_ptr_range halfedges(index_t polyline) const { | ||
| 855 | geo_debug_assert(polyline < nb()); | ||
| 856 | return const_index_ptr_range( | ||
| 857 | 81361 | H_, polyline_start_[polyline], polyline_start_[polyline+1] | |
| 858 | 81361 | ); | |
| 859 | } | ||
| 860 | |||
| 861 | /** | ||
| 862 | * \brief Creates a new polyline | ||
| 863 | */ | ||
| 864 | void begin_polyline() { | ||
| 865 | } | ||
| 866 | |||
| 867 | /** | ||
| 868 | * \brief Finishes a polyline creation | ||
| 869 | */ | ||
| 870 | void end_polyline() { | ||
| 871 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 146376 times.
|
146380 | polyline_start_.push_back(H_.size()); |
| 872 | 146380 | } | |
| 873 | |||
| 874 | /** | ||
| 875 | * \brief Adds a halfedge to the current polyline | ||
| 876 | * \details Needs to be called between begin_polyline() and | ||
| 877 | * end_polyline() | ||
| 878 | * \param[in] h the halfedge to be added to the current polyline | ||
| 879 | */ | ||
| 880 | void add_halfedge(index_t h) { | ||
| 881 |
3/4✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 492377 times.
✓ Branch 3 taken 1443 times.
|
493844 | H_.push_back(h); |
| 882 | } | ||
| 883 | |||
| 884 | /** | ||
| 885 | * \brief Gets the first vertex of a polyline | ||
| 886 | * \param[in] polyline a polyline index | ||
| 887 | * \return the index of the first vertex of \p polyline | ||
| 888 | */ | ||
| 889 |
1/2✓ Branch 0 taken 162722 times.
✗ Branch 1 not taken.
|
162722 | index_t first_vertex(index_t polyline) const { |
| 890 | 162722 | index_t h = H_[polyline_start_[polyline]]; | |
| 891 |
1/2✓ Branch 0 taken 162722 times.
✗ Branch 1 not taken.
|
162722 | return CF_.halfedges_.vertex(h,0); |
| 892 | } | ||
| 893 | |||
| 894 | /** | ||
| 895 | * \brief Gets the last vertex of a polyline | ||
| 896 | * \details if the polyline is closed, last_vertex() is the same as | ||
| 897 | * first_vertex() | ||
| 898 | * \param[in] polyline a polyline index | ||
| 899 | * \return the index of the first vertex of \p polyline | ||
| 900 | */ | ||
| 901 | 154551 | index_t last_vertex(index_t polyline) const { | |
| 902 |
1/2✓ Branch 0 taken 154551 times.
✗ Branch 1 not taken.
|
154551 | index_t h = H_[polyline_start_[polyline+1]-1]; |
| 903 |
1/2✓ Branch 0 taken 154551 times.
✗ Branch 1 not taken.
|
154551 | return CF_.halfedges_.vertex(h,1); |
| 904 | } | ||
| 905 | |||
| 906 | /** | ||
| 907 | * \brief Gets the predecessor of the first vertex | ||
| 908 | * \param[in] polyline a polyline | ||
| 909 | * \return if the polyline is closed, the predecessor | ||
| 910 | * of the first vertex, otherwise NO_INDEX | ||
| 911 | */ | ||
| 912 | 73190 | index_t prev_first_vertex(index_t polyline) const { | |
| 913 |
1/2✓ Branch 0 taken 73190 times.
✗ Branch 1 not taken.
|
73190 | if(first_vertex(polyline) != last_vertex(polyline)) { |
| 914 | return NO_INDEX; | ||
| 915 | } | ||
| 916 |
1/2✓ Branch 0 taken 73190 times.
✗ Branch 1 not taken.
|
73190 | index_t h = H_[polyline_start_[polyline+1]-1]; |
| 917 |
1/2✓ Branch 0 taken 73190 times.
✗ Branch 1 not taken.
|
146380 | return CF_.halfedges_.vertex(h,0); |
| 918 | } | ||
| 919 | |||
| 920 | private: | ||
| 921 | CoplanarFacets& CF_; | ||
| 922 | vector<index_t> H_; | ||
| 923 | vector<index_t> polyline_start_; | ||
| 924 | } polylines_; | ||
| 925 | |||
| 926 | }; | ||
| 927 | |||
| 928 | /**********************************************************************/ | ||
| 929 | |||
| 930 | } | ||
| 931 | |||
| 932 | #endif | ||
| 933 |