From b05916ee07015c30670fb6332de8e72e62318482 Mon Sep 17 00:00:00 2001 From: Bruno Levy Date: Wed, 25 Oct 2023 11:38:12 +0200 Subject: [PATCH] - Added missing "interval_nt::Rounding rounding;" in predicate filters (not needed with IntervalRN, but absolutely needed if we switch to IntervalRU) - Fix for #111: detect intersection vertices that land directly on existing vertex. --- src/lib/geogram/delaunay/CDT_2d.cpp | 23 ++++-- src/lib/geogram/delaunay/CDT_2d.h | 4 +- .../mesh/mesh_surface_intersection.cpp | 79 +++++++++++++------ .../mesh_surface_intersection_internal.cpp | 45 ++++++++++- .../mesh/mesh_surface_intersection_internal.h | 31 +++++--- src/lib/geogram/numerics/exact_geometry.cpp | 3 + src/lib/geogram/numerics/interval_nt.h | 2 +- 7 files changed, 144 insertions(+), 43 deletions(-) diff --git a/src/lib/geogram/delaunay/CDT_2d.cpp b/src/lib/geogram/delaunay/CDT_2d.cpp index fbf5db49ba8d..36cb0deff680 100644 --- a/src/lib/geogram/delaunay/CDT_2d.cpp +++ b/src/lib/geogram/delaunay/CDT_2d.cpp @@ -114,7 +114,6 @@ namespace GEO { Tecnstr_.resize(0); Tnext_.resize(0); Tprev_.resize(0); - clear_cache(); } void CDTBase2d::create_enclosing_triangle( @@ -151,7 +150,13 @@ namespace GEO { } } - void CDTBase2d::clear_cache() { + void CDTBase2d::begin_insert_transaction() { + } + + void CDTBase2d::commit_insert_transaction() { + } + + void CDTBase2d::rollback_insert_transaction() { } index_t CDTBase2d::insert(index_t v, index_t hint) { @@ -168,6 +173,7 @@ namespace GEO { } // Phase 1: find triangle that contains vertex i + begin_insert_transaction(); Sign o[3]; index_t t = locate(v,hint,o); int nb_z = (o[0] == ZERO) + (o[1] == ZERO) + (o[2] == ZERO); @@ -183,14 +189,17 @@ namespace GEO { v2T_.pop_back(); --nv_; } - // If we cached some predicates, we need to clear the cache, - // because everything cached relating to latest vertex becomes - // invalid since the same index will be used later by another - // point. - clear_cache(); + // Used by optional predicate cache management in derived classes. + // locate() computed some orient_2d() predicates, that may be + // stored in a temporary buffer, discard it. + rollback_insert_transaction(); return v; } + // Used by optional predicate cache management in derived classes. + // Copy the computed orient_2d() values to predicate cache. + commit_insert_transaction(); + // Stack of triangle edges to examine for flipping. Ignored in // non-Delaunay mode (ignored if !delaunay_) // Note: it is always edge 0 that we examine, since new diff --git a/src/lib/geogram/delaunay/CDT_2d.h b/src/lib/geogram/delaunay/CDT_2d.h index ab7fe7c775a7..dd12771d8f7c 100644 --- a/src/lib/geogram/delaunay/CDT_2d.h +++ b/src/lib/geogram/delaunay/CDT_2d.h @@ -206,7 +206,9 @@ namespace GEO { bool Tedge_is_Delaunay(index_t t, index_t le) const; protected: - virtual void clear_cache(); + virtual void begin_insert_transaction(); + virtual void commit_insert_transaction(); + virtual void rollback_insert_transaction(); /** * \brief Inserts a new point diff --git a/src/lib/geogram/mesh/mesh_surface_intersection.cpp b/src/lib/geogram/mesh/mesh_surface_intersection.cpp index afbe2239a219..ddffc568af37 100644 --- a/src/lib/geogram/mesh/mesh_surface_intersection.cpp +++ b/src/lib/geogram/mesh/mesh_surface_intersection.cpp @@ -591,15 +591,52 @@ namespace GEO { // Step 4: Epilogue // ---------------- - vector has_intersections(mesh_.facets.nb(), 0); - for(const IsectInfo& II: intersections) { - has_intersections[II.f1] = 1; - has_intersections[II.f2] = 1; + + // Vertices coming from intersections may land exactly + // on an existing vertex (see #111) + { + Stopwatch("Colocate newv"); + vector v2v(mesh_.vertices.nb()); + for(index_t v : mesh_.vertices) { + v2v[v] = v; + } + parallel_for( + 0, mesh_.vertices.nb(), + [&](index_t v) { + // If the vertex is an original vertex, not + // coming from an intersection, check whether + // it also exists as an intersection + if(vertex_to_exact_point_[v] == nullptr) { + const double* xyz = mesh_.vertices.point_ptr(v); + ExactPoint K(xyz[0], xyz[1], xyz[2], 1.0); + auto it = exact_point_to_vertex_.find(K); + if(it != exact_point_to_vertex_.end()) { + v2v[v] = it->second; + } + } + } + ); + for(index_t c : mesh_.facet_corners) { + index_t v = v2v[mesh_.facet_corners.vertex(c)]; + mesh_.facet_corners.set_vertex(c, v); + } } - - mesh_.facets.delete_elements(has_intersections); - // mesh_remove_bad_facets_no_check(mesh_); // TODO: needed here ? + // Remove original facets that have intersections. + { + vector has_intersections(mesh_.facets.nb(), 0); + for(const IsectInfo& II: intersections) { + has_intersections[II.f1] = 1; + has_intersections[II.f2] = 1; + } + mesh_.facets.delete_elements(has_intersections); + } + + // There can be duplicated facets coming from + // tesselated co-planar facets. + // Note: this updates operand_bit attribute + mesh_remove_bad_facets_no_check(mesh_); + if(use_radial_sort_) { build_Weiler_model(); } @@ -863,7 +900,7 @@ namespace GEO { if(h2 == h_ref_) { return POSITIVE; } - + for(const auto& c: refNorient_cache_) { if(c.first == h2) { return c.second; @@ -934,11 +971,6 @@ namespace GEO { void MeshSurfaceIntersection::build_Weiler_model() { - // There can be duplicated facets coming from - // tesselated co-planar facets. - // Note: this updates operand_bit attribute - mesh_remove_bad_facets_no_check(mesh_); - if(!facet_corner_alpha3_.is_bound()) { facet_corner_alpha3_.bind( mesh_.facet_corners.attributes(),"alpha3" @@ -1069,25 +1101,28 @@ namespace GEO { facet_corner_degenerate_[*it] = !OK; } - /* + // If we land here, it means we have co-planar overlapping // triangles, not supposed to happen after surface intersection, // but well, sometimes it happens ! Maybe due to underflows, // maybe due to a bug in the symbolic perturbation of the - // incircle predicate. + // incircle predicate. Maybe only due to #111 (now fixed), + // TODO: check whole Thingy10K + /* if(!OK) { std::cerr << std::endl; - for(auto it1=b; it1!=e; ++it1) { - for(auto it2=b; it2!=e; ++it2) { + for(auto it1=ib; it1!=ie; ++it1) { + for(auto it2=ib; it2!=ie; ++it2) { if(it1 != it2) { - std::cerr << (it1-b) << " " - << (it2-b) << std::endl; - RS.test(*it1, *it2); + std::cerr << (it1-ib) << " " + << (it2-ib) << std::endl; + // RS.test(*it1, *it2); } } } - save_radial("radial",b,e); - exit(-1); + if(ie-ib >= 3) { + save_radial(String::format("radial_%03d",k),ib,ie); + } } */ } diff --git a/src/lib/geogram/mesh/mesh_surface_intersection_internal.cpp b/src/lib/geogram/mesh/mesh_surface_intersection_internal.cpp index 444a67ef5c8c..274a9277b019 100644 --- a/src/lib/geogram/mesh/mesh_surface_intersection_internal.cpp +++ b/src/lib/geogram/mesh/mesh_surface_intersection_internal.cpp @@ -247,7 +247,8 @@ namespace GEO { exact_mesh_(EM), mesh_(EM.readonly_mesh()), f1_(index_t(-1)), - dry_run_(false) + dry_run_(false), + use_pred_cache_insert_buffer_(false) { #ifdef INTERSECTIONS_USE_EXACT_NT CDTBase2d::exact_incircle_ = true; @@ -259,6 +260,15 @@ namespace GEO { #endif } + void MeshInTriangle::clear() { + vertex_.resize(0); + edges_.resize(0); + f1_ = index_t(-1); + pred_cache_.clear(); + pred_cache_insert_buffer_.clear(); + use_pred_cache_insert_buffer_ = false; + CDTBase2d::clear(); + } void MeshInTriangle::begin_facet(index_t f) { f1_ = f; @@ -430,11 +440,42 @@ namespace GEO { } return result; } - + + void MeshInTriangle::begin_insert_transaction() { + use_pred_cache_insert_buffer_ = true; + } + + void MeshInTriangle::commit_insert_transaction() { + for(const auto& it: pred_cache_insert_buffer_) { + pred_cache_[it.first] = it.second; + } + pred_cache_insert_buffer_.resize(0); + use_pred_cache_insert_buffer_ = false; + } + + void MeshInTriangle::rollback_insert_transaction() { + pred_cache_insert_buffer_.resize(0); + use_pred_cache_insert_buffer_ = false; + } Sign MeshInTriangle::orient2d(index_t vx1,index_t vx2,index_t vx3) const { trindex K(vx1, vx2, vx3); + + if(use_pred_cache_insert_buffer_) { + Sign result = PCK::orient_2d_projected( + vertex_[K.indices[0]].point_exact, + vertex_[K.indices[1]].point_exact, + vertex_[K.indices[2]].point_exact, + f1_normal_axis_ + ); + pred_cache_insert_buffer_.push_back(std::make_pair(K, result)); + if(odd_order(vx1,vx2,vx3)) { + result = Sign(-result); + } + return result; + } + bool inserted; std::map::iterator it; std::tie(it,inserted) = pred_cache_.insert(std::make_pair(K,ZERO)); diff --git a/src/lib/geogram/mesh/mesh_surface_intersection_internal.h b/src/lib/geogram/mesh/mesh_surface_intersection_internal.h index c55bcf79583b..2c2e53302df7 100644 --- a/src/lib/geogram/mesh/mesh_surface_intersection_internal.h +++ b/src/lib/geogram/mesh/mesh_surface_intersection_internal.h @@ -309,17 +309,11 @@ namespace GEO { */ void commit(); - - void clear() override { - vertex_.resize(0); - edges_.resize(0); - f1_ = index_t(-1); - CDTBase2d::clear(); - } - void clear_cache() override { - pred_cache_.clear(); - } + /** + * \see CDT2d::clear() + */ + void clear() override; protected: /** @@ -425,6 +419,20 @@ namespace GEO { public: void save(const std::string& filename) const override; + protected: + void begin_insert_transaction() override; + void commit_insert_transaction() override; + void rollback_insert_transaction() override; + + void begin_pred_cache_transaction() { + use_pred_cache_insert_buffer_ = true; + } + void end_pred_cache_transaction() { + use_pred_cache_insert_buffer_ = false; + } + void commit_pred_cache_insert_buffer(); + void clear_pred_cache_insert_buffer(); + private: MeshSurfaceIntersection& exact_mesh_; const Mesh& mesh_; @@ -439,6 +447,9 @@ namespace GEO { bool has_planar_isect_; bool dry_run_; mutable std::map pred_cache_; + bool use_pred_cache_insert_buffer_; + mutable std::vector< std::pair > + pred_cache_insert_buffer_; }; /*************************************************************************/ diff --git a/src/lib/geogram/numerics/exact_geometry.cpp b/src/lib/geogram/numerics/exact_geometry.cpp index d3889b0bd623..6b08485fe2b0 100644 --- a/src/lib/geogram/numerics/exact_geometry.cpp +++ b/src/lib/geogram/numerics/exact_geometry.cpp @@ -259,6 +259,7 @@ namespace GEO { // Filter { + interval_nt::Rounding rounding; vec3HI p0I(p0); vec3HI U = vec3HI(p1)-p0I; vec3HI V = vec3HI(p2)-p0I; @@ -497,6 +498,7 @@ namespace GEO { // Filter { + interval_nt::Rounding rounding; interval_nt l3I(l3); interval_nt L1 = interval_nt(l0) - l3I; interval_nt L2 = interval_nt(l1) - l3I; @@ -597,6 +599,7 @@ namespace GEO { ) { // Filter using interval arithmetics { + interval_nt::Rounding rounding; vec3I p1I(p1); vec3I U = vec3I(p2) - p1I; vec3I V = vec3I(p3) - p1I; diff --git a/src/lib/geogram/numerics/interval_nt.h b/src/lib/geogram/numerics/interval_nt.h index fff7c2426fec..2664b8c0c384 100644 --- a/src/lib/geogram/numerics/interval_nt.h +++ b/src/lib/geogram/numerics/interval_nt.h @@ -198,7 +198,7 @@ namespace GEO { geo_argused(inf); geo_argused(sup); } -#endif +#endif };