Merge remote-tracking branch 'cgal/master' into PMP-Local_stitch_borders-GF

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/Weights.h

@@ -141,9 +141,9 @@ public:
     Vector a = get(ppmap(), v0) - get(ppmap(), v1);
     Vector b = get(ppmap(), v2) - get(ppmap(), v1);
 
-    double dot_ab = a*b;
-    double dot_aa = a.squared_length();
-    double dot_bb = b.squared_length();
+    double dot_ab = CGAL::to_double(a * b);
+    double dot_aa = CGAL::to_double(a.squared_length());
+    double dot_bb = CGAL::to_double(b.squared_length());
     double lb = -0.999, ub = 0.999;
     double cosine = dot_ab / CGAL::sqrt(dot_aa) / CGAL::sqrt(dot_bb);
     cosine = (cosine < lb) ? lb : cosine;
@@ -853,6 +853,45 @@ public:
   }
 };
 
+// Cotangent_value_Meyer has been changed to the version:
+// Cotangent_value_Meyer_secure to avoid imprecisions from
+// the issue #4706 - https://github.com/CGAL/cgal/issues/4706.
+template<
+class PolygonMesh, class VertexPointMap = typename boost::property_map<PolygonMesh, vertex_point_t>::type>
+class Cotangent_weight_with_voronoi_area_fairing_secure {
+
+  typedef PolygonMesh PM;
+  typedef VertexPointMap VPMap;
+  Voronoi_area<PM, VPMap> voronoi_functor;
+  Cotangent_weight<PM, VPMap, Cotangent_value_Meyer_secure<PM, VPMap> > cotangent_functor;
+
+public:
+  Cotangent_weight_with_voronoi_area_fairing_secure(PM& pmesh_) :
+  voronoi_functor(pmesh_, get(CGAL::vertex_point, pmesh_)),
+  cotangent_functor(pmesh_, get(CGAL::vertex_point, pmesh_))
+  { }
+
+  Cotangent_weight_with_voronoi_area_fairing_secure(PM& pmesh_, VPMap vpmap_) :
+  voronoi_functor(pmesh_, vpmap_),
+  cotangent_functor(pmesh_, vpmap_)
+  { }
+
+  PM& pmesh() {
+    return voronoi_functor.pmesh();
+  }
+
+  typedef typename boost::graph_traits<PM>::halfedge_descriptor halfedge_descriptor;
+  typedef typename boost::graph_traits<PM>::vertex_descriptor vertex_descriptor;
+
+  double w_i(vertex_descriptor v_i) {
+    return 0.5 / voronoi_functor(v_i);
+  }
+
+  double w_ij(halfedge_descriptor he) {
+    return cotangent_functor(he) * 2.0;
+  }
+};
+
 template<class PolygonMesh>
 class Uniform_weight_fairing
 {

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/distance.h

@@ -29,7 +29,7 @@
 #include <CGAL/spatial_sort.h>
 
 #ifdef CGAL_LINKED_WITH_TBB
-#include <tbb/parallel_for.h>
+#include <tbb/parallel_reduce.h>
 #include <tbb/blocked_range.h>
 #include <atomic>
 #endif // CGAL_LINKED_WITH_TBB
@@ -81,21 +81,27 @@ struct Distance_computation{
   const AABB_tree& tree;
   const PointRange& sample_points;
   Point_3 initial_hint;
-  std::atomic<double>* distance;
-
+  double distance;
+  //constructor
   Distance_computation(
           const AABB_tree& tree,
           const Point_3& p,
-          const PointRange& sample_points,
-          std::atomic<double>* d)
+          const PointRange& sample_points)
     : tree(tree)
     , sample_points(sample_points)
     , initial_hint(p)
-    , distance(d)
+    , distance(-1)
+  {}
+  //split constructor
+  Distance_computation(Distance_computation& s, tbb::split )
+    : tree(s.tree)
+    , sample_points(s.sample_points)
+    , initial_hint(s.initial_hint)
+    , distance(-1)
   {}
 
   void
-  operator()(const tbb::blocked_range<std::size_t>& range) const
+  operator()(const tbb::blocked_range<std::size_t>& range)
   {
     Point_3 hint = initial_hint;
     double hdist = 0;
@@ -107,15 +113,11 @@ struct Distance_computation{
       if(d > hdist)
         hdist=d;
     }
-
-    // update max value stored in distance
-    double current_value = *distance;
-    while( current_value < hdist )
-    {
-      if(distance->compare_exchange_weak(current_value, hdist))
-        current_value = hdist;
-    }
+    if(hdist > distance)
+      distance = hdist;
   }
+
+  void join( Distance_computation& rhs ) {distance = (std::max)(rhs.distance, distance); }
 };
 #endif
 
@@ -136,9 +138,9 @@ double approximate_Hausdorff_distance_impl(
   {
     std::atomic<double> distance;
     distance=0;
-    Distance_computation<AABBTree, PointRange> f(tree, hint, sample_points, &distance);
-    tbb::parallel_for(tbb::blocked_range<std::size_t>(0, sample_points.size()), f);
-    return distance;
+    Distance_computation<AABBTree, PointRange> f(tree, hint, sample_points);
+    tbb::parallel_reduce(tbb::blocked_range<std::size_t>(0, sample_points.size()), f);
+    return f.distance;
   }
   else
 #endif

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/fair.h

@@ -137,7 +137,11 @@ namespace internal {
 #endif
 
     typedef typename GetVertexPointMap < TriangleMesh, NamedParameters>::type VPMap;
-    typedef CGAL::internal::Cotangent_weight_with_voronoi_area_fairing<TriangleMesh, VPMap>
+
+    // Cotangent_weight_with_voronoi_area_fairing has been changed to the version:
+    // Cotangent_weight_with_voronoi_area_fairing_secure to avoid imprecisions from
+    // the issue #4706 - https://github.com/CGAL/cgal/issues/4706.
+    typedef CGAL::internal::Cotangent_weight_with_voronoi_area_fairing_secure<TriangleMesh, VPMap>
       Default_Weight_calculator;
 
     VPMap vpmap_ = choose_parameter(get_parameter(np, internal_np::vertex_point),

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Face_graph_output_builder.h

@@ -682,6 +682,10 @@ public:
       halfedge_descriptor h1 = it->second.first[&tm1];
       halfedge_descriptor h2 = it->second.first[&tm2];
 
+#ifdef CGAL_COREFINEMENT_DEBUG
+      std::cout << "Looking at triangles around edge " << tm1.point(source(h1, tm1)) << " " << tm1.point(target(h1, tm1)) << "\n";
+#endif
+
       CGAL_assertion(ids.first==vertex_to_node_id1[source(h1,tm1)]);
       CGAL_assertion(ids.second==vertex_to_node_id1[target(h1,tm1)]);
       CGAL_assertion(ids.first==vertex_to_node_id2[source(h2,tm2)]);
@@ -713,6 +717,9 @@ public:
             //Nothing allowed
             if (!used_to_clip_a_surface)
             {
+#ifdef CGAL_COREFINEMENT_DEBUG
+              std::cout << "  Non-manifold edge case 1\n";
+#endif
               impossible_operation.set();
               return;
             }
@@ -723,6 +730,9 @@ public:
           //Ambiguous, we can do nothing
           if (!used_to_clip_a_surface)
           {
+#ifdef CGAL_COREFINEMENT_DEBUG
+              std::cout << "  Non-manifold edge case 2\n";
+#endif
             impossible_operation.set();
             return;
           }
@@ -776,6 +786,9 @@ public:
         {
           CGAL_assertion(!used_to_clip_a_surface);
           //Ambiguous, we do nothing
+#ifdef CGAL_COREFINEMENT_DEBUG
+              std::cout << "  Non-manifold edge case 3\n";
+#endif
           impossible_operation.set();
           return;
         }
@@ -971,6 +984,9 @@ public:
                 // poly_second - poly_first             = {0}
                 // poly_first \cap poly_second          = q1q2
                 // opposite( poly_first U poly_second ) = p2p1
+#ifdef CGAL_COREFINEMENT_DEBUG
+              std::cout << "  Non-manifold edge case 4\n";
+#endif
                 impossible_operation.set(TM1_MINUS_TM2); // tm1-tm2 is non-manifold
               }
               else{
@@ -982,7 +998,12 @@ public:
                 is_patch_inside_tm2.set(patch_id_p1);
                 is_patch_inside_tm2.set(patch_id_p2);
                 if (!used_to_clip_a_surface)
+                {
+#ifdef CGAL_COREFINEMENT_DEBUG
+              std::cout << "  Non-manifold edge case 5\n";
+#endif
                   impossible_operation.set(INTERSECTION); // tm1 n tm2 is non-manifold
+                }
               }
             }
             else
@@ -997,6 +1018,9 @@ public:
               {
                 if (!used_to_clip_a_surface)
                 {
+#ifdef CGAL_COREFINEMENT_DEBUG
+              std::cout << "  Non-manifold edge case 6\n";
+#endif
                   impossible_operation.set();
                   return;
                 }
@@ -1018,6 +1042,9 @@ public:
               {
                 if (!used_to_clip_a_surface)
                 {
+#ifdef CGAL_COREFINEMENT_DEBUG
+              std::cout << "  Non-manifold edge case 7\n";
+#endif
                   impossible_operation.set();
                   return;
                 }
@@ -1039,6 +1066,9 @@ public:
                 // poly_second - poly_first             = q1q2
                 // poly_first \cap poly_second          = {0}
                 // opposite( poly_first U poly_second ) = p2q1 U q2p1
+#ifdef CGAL_COREFINEMENT_DEBUG
+              std::cout << "  Non-manifold edge case 8\n";
+#endif
                 impossible_operation.set(UNION); // tm1 U tm2 is non-manifold
               }
               else{
@@ -1049,6 +1079,9 @@ public:
                 // poly_second - poly_first             = q1p1 U p2q2
                 // poly_first \cap poly_second          = p1p2
                 // opposite( poly_first U poly_second ) = q2q1
+#ifdef CGAL_COREFINEMENT_DEBUG
+              std::cout << "  Non-manifold edge case 9\n";
+#endif
                 impossible_operation.set(TM2_MINUS_TM1); // tm2 - tm1 is non-manifold
               }
             }

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_impl.h

@@ -35,6 +35,7 @@
 #include <iostream>
 #include <utility>
 #include <vector>
+#include <unordered_map>
 
 namespace CGAL {
 namespace Polygon_mesh_processing {

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_self_intersections.h

@@ -30,6 +30,7 @@
 #include <CGAL/boost/graph/Named_function_parameters.h>
 #include <CGAL/boost/graph/named_params_helper.h>
 #include <CGAL/boost/graph/selection.h>
+#include <CGAL/box_intersection_d.h>
 #include <CGAL/utility.h>
 
 #include <array>
@@ -62,6 +63,7 @@ namespace Polygon_mesh_processing {
 namespace internal {
 
 #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
+static int unsolved_self_intersections = 0;
 static int self_intersections_solved_by_constrained_smoothing = 0;
 static int self_intersections_solved_by_unconstrained_smoothing = 0;
 static int self_intersections_solved_by_constrained_hole_filling = 0;
@@ -252,8 +254,7 @@ FaceOutputIterator replace_faces_with_patch(const std::vector<typename boost::gr
   CGAL_postcondition(is_valid_polygon_mesh(pmesh));
 
 #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
-  if(Polygon_mesh_processing::does_self_intersect(new_faces, pmesh))
-    std::cout << "!! NEW FACES SELF INTERSECT !!" << std::endl;
+  CGAL_postcondition(!does_self_intersect(new_faces, pmesh));
 #endif
 
   return out;
@@ -948,101 +949,6 @@ bool construct_tentative_sub_hole_patch(std::vector<std::vector<typename boost::
 
 // -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
 
-// This function is only called when the hole is NOT subdivided into smaller holes
-template <typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
-bool fill_hole(std::vector<typename boost::graph_traits<TriangleMesh>::halfedge_descriptor>& cc_border_hedges,
-               std::set<typename boost::graph_traits<TriangleMesh>::face_descriptor>& cc_faces,
-               std::set<typename boost::graph_traits<TriangleMesh>::face_descriptor>& working_face_range,
-               TriangleMesh& tmesh,
-               VertexPointMap vpm,
-               const GeomTraits& gt)
-{
-  typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor     vertex_descriptor;
-  typedef typename boost::graph_traits<TriangleMesh>::edge_descriptor       edge_descriptor;
-  typedef typename boost::graph_traits<TriangleMesh>::face_descriptor       face_descriptor;
-
-  typedef typename boost::property_traits<VertexPointMap>::value_type       Point;
-
-#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
-  std::cout << "  DEBUG: Attempting hole-filling (no constraints), " << cc_faces.size() << " faces\n";
-#endif
-
-  if(!order_border_halfedge_range(cc_border_hedges, tmesh))
-  {
-    CGAL_assertion(false); // we shouldn't fail to orient the boundary cycle of the complete hole
-    return false;
-  }
-
-  std::set<vertex_descriptor> cc_interior_vertices;
-  std::set<edge_descriptor> cc_interior_edges;
-
-  std::vector<vertex_descriptor> cc_border_vertices;
-  cc_border_vertices.reserve(cc_border_hedges.size());
-
-  std::vector<std::vector<Point> > patch;
-  if(!construct_tentative_hole_patch(cc_border_vertices, cc_interior_vertices, cc_interior_edges,
-                                     cc_border_hedges, cc_faces, patch, tmesh, vpm, gt))
-  {
-#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
-    std::cout << "  DEBUG: Failed to find acceptable hole patch\n";
-#endif
-
-    return false;
-  }
-
-  // Could renew the range directly within the patch replacement function
-  // to avoid erasing and re-adding the same face
-  for(const face_descriptor f : cc_faces)
-    working_face_range.erase(f);
-
-  // Plug the new triangles in the mesh, reusing previous edges and faces
-  replace_faces_with_patch(cc_border_vertices, cc_interior_vertices,
-                           cc_border_hedges, cc_interior_edges,
-                           cc_faces, patch, tmesh, vpm,
-                           std::inserter(working_face_range, working_face_range.end()));
-
-#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT
-  static int filed_hole_id = 0;
-  std::stringstream oss;
-  oss << "results/filled_basic_" << filed_hole_id++ << ".off" << std::ends;
-  std::ofstream(oss.str().c_str()) << std::setprecision(17) << tmesh;
-#endif
-
-  CGAL_postcondition(is_valid_polygon_mesh(tmesh));
-
-  return true;
-}
-
-// Same function as above but border of the hole is not known
-template <typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
-bool fill_hole(std::set<typename boost::graph_traits<TriangleMesh>::face_descriptor>& cc_faces,
-               std::set<typename boost::graph_traits<TriangleMesh>::face_descriptor>& working_face_range,
-               TriangleMesh& tmesh,
-               VertexPointMap vpm,
-               const GeomTraits& gt)
-{
-  typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor   halfedge_descriptor;
-  typedef typename boost::graph_traits<TriangleMesh>::face_descriptor       face_descriptor;
-
-  std::vector<halfedge_descriptor> cc_border_hedges;
-  for(face_descriptor fd : cc_faces)
-  {
-    halfedge_descriptor h = halfedge(fd, tmesh);
-    for(int i=0; i<3; ++i)
-    {
-      if(is_border(opposite(h, tmesh), tmesh) || cc_faces.count(face(opposite(h, tmesh), tmesh)) == 0)
-        cc_border_hedges.push_back(h);
-
-      h = next(h, tmesh);
-    }
-  }
-
-  if(order_border_halfedge_range(cc_border_hedges, tmesh))
-    return fill_hole(cc_border_hedges, cc_faces, working_face_range, tmesh, vpm, gt);
-  else
-    return false;
-}
-
 // Patch is not valid if:
 // - we insert the same face more than once
 // - insert (geometric) non-manifold edges
@@ -1120,6 +1026,102 @@ bool check_patch_sanity(const std::vector<std::vector<Point> >& patch)
   return true;
 }
 
+// This function is only called when the hole is NOT subdivided into smaller holes
+template <typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
+bool fill_hole(std::vector<typename boost::graph_traits<TriangleMesh>::halfedge_descriptor>& cc_border_hedges,
+               std::set<typename boost::graph_traits<TriangleMesh>::face_descriptor>& cc_faces,
+               std::set<typename boost::graph_traits<TriangleMesh>::face_descriptor>& working_face_range,
+               TriangleMesh& tmesh,
+               VertexPointMap vpm,
+               const GeomTraits& gt)
+{
+  typedef typename boost::graph_traits<TriangleMesh>::vertex_descriptor     vertex_descriptor;
+  typedef typename boost::graph_traits<TriangleMesh>::edge_descriptor       edge_descriptor;
+  typedef typename boost::graph_traits<TriangleMesh>::face_descriptor       face_descriptor;
+
+  typedef typename boost::property_traits<VertexPointMap>::value_type       Point;
+
+#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
+  std::cout << "  DEBUG: Attempting hole-filling (no constraints), " << cc_faces.size() << " faces\n";
+#endif
+
+  if(!order_border_halfedge_range(cc_border_hedges, tmesh))
+  {
+    CGAL_assertion(false); // we shouldn't fail to orient the boundary cycle of the complete hole
+    return false;
+  }
+
+  std::set<vertex_descriptor> cc_interior_vertices;
+  std::set<edge_descriptor> cc_interior_edges;
+
+  std::vector<vertex_descriptor> cc_border_vertices;
+  cc_border_vertices.reserve(cc_border_hedges.size());
+
+  std::vector<std::vector<Point> > patch;
+  if(!construct_tentative_hole_patch(cc_border_vertices, cc_interior_vertices, cc_interior_edges,
+                                     cc_border_hedges, cc_faces, patch, tmesh, vpm, gt) ||
+     !check_patch_sanity<TriangleMesh>(patch))
+  {
+#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
+    std::cout << "  DEBUG: Failed to find acceptable hole patch\n";
+#endif
+
+    return false;
+  }
+
+  // Could renew the range directly within the patch replacement function
+  // to avoid erasing and re-adding the same face
+  for(const face_descriptor f : cc_faces)
+    working_face_range.erase(f);
+
+  // Plug the new triangles in the mesh, reusing previous edges and faces
+  replace_faces_with_patch(cc_border_vertices, cc_interior_vertices,
+                           cc_border_hedges, cc_interior_edges,
+                           cc_faces, patch, tmesh, vpm,
+                           std::inserter(working_face_range, working_face_range.end()));
+
+#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT
+  static int filed_hole_id = 0;
+  std::stringstream oss;
+  oss << "results/filled_basic_" << filed_hole_id++ << ".off" << std::ends;
+  std::ofstream(oss.str().c_str()) << std::setprecision(17) << tmesh;
+#endif
+
+  CGAL_postcondition(is_valid_polygon_mesh(tmesh));
+
+  return true;
+}
+
+// Same function as above but border of the hole is not known
+template <typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
+bool fill_hole(std::set<typename boost::graph_traits<TriangleMesh>::face_descriptor>& cc_faces,
+               std::set<typename boost::graph_traits<TriangleMesh>::face_descriptor>& working_face_range,
+               TriangleMesh& tmesh,
+               VertexPointMap vpm,
+               const GeomTraits& gt)
+{
+  typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor   halfedge_descriptor;
+  typedef typename boost::graph_traits<TriangleMesh>::face_descriptor       face_descriptor;
+
+  std::vector<halfedge_descriptor> cc_border_hedges;
+  for(face_descriptor fd : cc_faces)
+  {
+    halfedge_descriptor h = halfedge(fd, tmesh);
+    for(int i=0; i<3; ++i)
+    {
+      if(is_border(opposite(h, tmesh), tmesh) || cc_faces.count(face(opposite(h, tmesh), tmesh)) == 0)
+        cc_border_hedges.push_back(h);
+
+      h = next(h, tmesh);
+    }
+  }
+
+  if(order_border_halfedge_range(cc_border_hedges, tmesh))
+    return fill_hole(cc_border_hedges, cc_faces, working_face_range, tmesh, vpm, gt);
+  else
+    return false;
+}
+
 template <typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
 bool fill_hole_with_constraints(std::vector<typename boost::graph_traits<TriangleMesh>::halfedge_descriptor>& cc_border_hedges,
                                 std::set<typename boost::graph_traits<TriangleMesh>::face_descriptor>& cc_faces,
@@ -1161,8 +1163,9 @@ bool fill_hole_with_constraints(std::vector<typename boost::graph_traits<Triangl
                                                  CGAL::parameters::edge_is_constrained_map(eif));
 
     visited_faces.insert(sub_cc.begin(), sub_cc.end());
-
+#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
     std::cout << "CC of size " << sub_cc.size() << " (total: " << cc_faces.size() << ")" << std::endl;
+#endif
     ++cc_counter;
 
 #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT
@@ -1176,9 +1179,9 @@ bool fill_hole_with_constraints(std::vector<typename boost::graph_traits<Triangl
       return fill_hole(cc_border_hedges, cc_faces, working_face_range, tmesh, vpm, gt);
     }
   }
-
+#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
   std::cout << cc_counter << " independent sub holes" << std::endl;
-
+#endif
 #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT
   std::ofstream out("results/hole_fillers.off");
   out.precision(17);
@@ -1204,7 +1207,9 @@ bool fill_hole_with_constraints(std::vector<typename boost::graph_traits<Triangl
   // We're assembling multiple patches so we could have the same face appearing multiple times...
   if(!check_patch_sanity<TriangleMesh>(patch))
   {
+#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
     std::cout << "Unhealthy patch, use base fill_hole" << std::endl;
+#endif
     return fill_hole(cc_border_hedges, cc_faces, working_face_range, tmesh, vpm, gt);
   }
 
@@ -1224,6 +1229,94 @@ bool fill_hole_with_constraints(std::vector<typename boost::graph_traits<Triangl
   return true;
 }
 
+template <class Box, class TM, class VPM, class GT, class OutputIterator>
+struct Strict_intersect_edges // "strict" as in "not sharing a vertex"
+{
+  typedef typename boost::graph_traits<TM>::halfedge_descriptor               halfedge_descriptor;
+  typedef typename GT::Segment_3                                              Segment;
+
+  mutable OutputIterator m_iterator;
+  const TM& m_tmesh;
+  const VPM m_vpmap;
+
+  typename GT::Construct_segment_3 m_construct_segment;
+  typename GT::Do_intersect_3 m_do_intersect;
+
+  Strict_intersect_edges(const TM& tmesh, VPM vpmap, const GT& gt, OutputIterator it)
+    :
+      m_iterator(it),
+      m_tmesh(tmesh),
+      m_vpmap(vpmap),
+      m_construct_segment(gt.construct_segment_3_object()),
+      m_do_intersect(gt.do_intersect_3_object())
+  {}
+
+  void operator()(const Box* b, const Box* c) const
+  {
+    const halfedge_descriptor h = b->info();
+    const halfedge_descriptor g = c->info();
+
+    if(source(h, m_tmesh) == target(g, m_tmesh) || target(h, m_tmesh) == source(g, m_tmesh))
+      return;
+
+    const Segment s1 = m_construct_segment(get(m_vpmap, source(h, m_tmesh)), get(m_vpmap, target(h, m_tmesh)));
+    const Segment s2 = m_construct_segment(get(m_vpmap, source(g, m_tmesh)), get(m_vpmap, target(g, m_tmesh)));
+
+    if(m_do_intersect(s1, s2))
+      *m_iterator++ = std::make_pair(b->info(), c->info());
+  }
+};
+
+template <typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
+bool is_simple_3(const std::vector<typename boost::graph_traits<TriangleMesh>::halfedge_descriptor>& cc_border_hedges,
+                 const TriangleMesh& tmesh,
+                 VertexPointMap vpm,
+                 const GeomTraits& gt)
+{
+  typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor                       halfedge_descriptor;
+
+  typedef typename boost::property_traits<VertexPointMap>::reference                            Point_ref;
+
+  typedef CGAL::Box_intersection_d::ID_FROM_BOX_ADDRESS                                         Box_policy;
+  typedef CGAL::Box_intersection_d::Box_with_info_d<double, 3, halfedge_descriptor, Box_policy> Box;
+
+  std::vector<Box> boxes;
+  boxes.reserve(cc_border_hedges.size());
+
+  for(halfedge_descriptor h : cc_border_hedges)
+  {
+    const Point_ref p = get(vpm, source(h, tmesh));
+    const Point_ref q = get(vpm, target(h, tmesh));
+    CGAL_assertion(!gt.equal_3_object()(p, q));
+
+    boxes.emplace_back(p.bbox() + q.bbox(), h);
+  }
+
+  // generate box pointers
+  std::vector<const Box*> box_ptr;
+  box_ptr.reserve(boxes.size());
+
+  for(Box& b : boxes)
+    box_ptr.push_back(&b);
+
+  typedef boost::function_output_iterator<CGAL::internal::Throw_at_output>          Throwing_output_iterator;
+  typedef internal::Strict_intersect_edges<Box, TriangleMesh, VertexPointMap,
+                                           GeomTraits, Throwing_output_iterator>    Throwing_filter;
+  Throwing_filter throwing_filter(tmesh, vpm, gt, Throwing_output_iterator());
+
+  try
+  {
+    const std::ptrdiff_t cutoff = 2000;
+    CGAL::box_self_intersection_d<Parallel_if_available_tag>(box_ptr.begin(), box_ptr.end(), throwing_filter, cutoff);
+  }
+  catch(CGAL::internal::Throw_at_output_exception&)
+  {
+    return false;
+  }
+
+  return true;
+}
+
 template <typename TriangleMesh, typename VertexPointMap, typename GeomTraits>
 bool remove_self_intersections_with_hole_filling(std::vector<typename boost::graph_traits<TriangleMesh>::halfedge_descriptor>& cc_border_hedges,
                                                  std::set<typename boost::graph_traits<TriangleMesh>::face_descriptor>& cc_faces,
@@ -1243,6 +1336,14 @@ bool remove_self_intersections_with_hole_filling(std::vector<typename boost::gra
   out.close();
 #endif
 
+  if(!is_simple_3(cc_border_hedges, tmesh, vpm, gt))
+  {
+#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
+    std::cout << "Hole filling cannot handle non-simple border" << std::endl;
+#endif
+    return false;
+  }
+
 #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTIONS_NO_CONSTRAINTS_IN_HOLE_FILLING
   // Do not try to impose sharp edge constraints if we are not doing local-only self intersections removal
   local_self_intersection_removal = false;
@@ -1294,7 +1395,7 @@ remove_self_intersections_one_step(std::set<typename boost::graph_traits<Triangl
   std::set<face_descriptor> faces_to_remove_copy = faces_to_remove;
 
 #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
-  std::cout << "DEBUG: running remove_self_intersections_one_step, step " << step
+  std::cout << "##### running remove_self_intersections_one_step, step " << step
             << " with " << faces_to_remove.size() << " intersecting faces\n";
 #endif
 
@@ -1307,9 +1408,12 @@ remove_self_intersections_one_step(std::set<typename boost::graph_traits<Triangl
 #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
   std::cout << "  DEBUG: is_valid in one_step(tmesh)? ";
   std::cout.flush();
-  std::cout << is_valid_polygon_mesh(tmesh) << "\n";
+
+  unsolved_self_intersections = 0;
 #endif
 
+  CGAL_precondition(is_valid_polygon_mesh(tmesh));
+
   while(!faces_to_remove.empty())
   {
 #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
@@ -1505,6 +1609,8 @@ remove_self_intersections_one_step(std::set<typename boost::graph_traits<Triangl
       topology_issue = true;
 #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
       std::cout << "  DEBUG: CC not handled due to the presence at least one non-manifold vertex\n";
+
+      ++unsolved_self_intersections;
 #endif
 
       continue; // cannot replace a patch containing a nm vertex by a disk
@@ -1530,7 +1636,12 @@ remove_self_intersections_one_step(std::set<typename boost::graph_traits<Triangl
     }
 
     if(cc_faces.size() == 1) // it is a triangle nothing better can be done
+    {
+#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
+     ++unsolved_self_intersections;
+#endif
       continue;
+    }
 
     working_face_range.insert(cc_faces.begin(), cc_faces.end());
 
@@ -1633,8 +1744,9 @@ remove_self_intersections_one_step(std::set<typename boost::graph_traits<Triangl
       if(nb_cycles > (only_border_edges ? 1 : 0))
       {
 #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
-        std::cout << "  DEBUG: CC not handled due to the presence of  "
+        std::cout << "  DEBUG: CC not handled due to the presence of "
                   << nb_cycles << " of boundary edges\n";
+     ++unsolved_self_intersections;
 #endif
 
         topology_issue = true;
@@ -1646,6 +1758,7 @@ remove_self_intersections_one_step(std::set<typename boost::graph_traits<Triangl
         {
 #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
           std::cout << "  DEBUG: CC not handled because it is not a topological disk (preserve_genus=true)\n";
+          ++unsolved_self_intersections;
 #endif
 
           all_fixed = false;
@@ -1681,6 +1794,7 @@ remove_self_intersections_one_step(std::set<typename boost::graph_traits<Triangl
 #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
           std::cout << "  DEBUG: CC not handled because it is not a topological disk("
                     << nbc << " boundary cycles)\n";
+          ++unsolved_self_intersections;
 #endif
 
           all_fixed = false;
@@ -1702,6 +1816,7 @@ remove_self_intersections_one_step(std::set<typename boost::graph_traits<Triangl
     {
 #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
       std::cout << "  DEBUG: Failed to fill hole\n";
+      ++unsolved_self_intersections;
 #endif
 
       all_fixed = false;
@@ -1796,9 +1911,9 @@ bool remove_self_intersections(const FaceRange& face_range,
       // TODO : possible optimization to reduce the range to check with the bbox
       // of the previous patches or something.
       self_intersections(working_face_range, tmesh, std::back_inserter(self_inter));
-
+#ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
       std::cout << self_inter.size() << " intersecting pairs" << std::endl;
-
+#endif
       for(const Face_pair& fp : self_inter)
       {
         faces_to_remove.insert(fp.first);
@@ -1831,6 +1946,7 @@ bool remove_self_intersections(const FaceRange& face_range,
   std::cout << "solved by unconstrained smoothing: " << internal::self_intersections_solved_by_unconstrained_smoothing << std::endl;
   std::cout << "solved by constrained hole-filling: " << internal::self_intersections_solved_by_constrained_hole_filling << std::endl;
   std::cout << "solved by unconstrained hole-filling: " << internal::self_intersections_solved_by_unconstrained_hole_filling << std::endl;
+  std::cout << "unsolved: " << internal::unsolved_self_intersections << std::endl;
 #endif
 
 #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT