- 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

@@ -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

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 

src/lib/geogram/mesh/mesh_surface_intersection.cpp

@@ -591,15 +591,52 @@ namespace GEO {
         // Step 4: Epilogue
         // ----------------
 
-        vector<index_t> 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<index_t> 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<index_t> 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);
+                            }
                         }
                         */
                     }

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<trindex, Sign>::iterator it;
         std::tie(it,inserted) = pred_cache_.insert(std::make_pair(K,ZERO));

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<trindex, Sign> pred_cache_;
+        bool use_pred_cache_insert_buffer_;
+        mutable std::vector< std::pair<trindex, Sign> >
+             pred_cache_insert_buffer_;
     };
 
     /*************************************************************************/

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;

src/lib/geogram/numerics/interval_nt.h

@@ -198,7 +198,7 @@ namespace GEO {
             geo_argused(inf);
             geo_argused(sup);
         }
-#endif        
+#endif
     };