Merge pull request #6605 from sloriot/PMP-snap_no_deg_faces_created

Avoid creating degenerate faces in snap

BGL/include/CGAL/boost/graph/Euler_operations.h

@@ -1432,6 +1432,106 @@ add_face_to_border(typename boost::graph_traits<Graph>::halfedge_descriptor h1,
   return newh;
 }
 
+/**
+ *  \returns `true` if `e` satisfies the *link condition* \cgalCite{degn-tpec-98}, which guarantees that the surface is also 2-manifold after the edge collapse.
+ */
+template<typename Graph>
+bool
+does_satisfy_link_condition(typename boost::graph_traits<Graph>::edge_descriptor e,
+                            const Graph& g)
+{
+  typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
+  typedef typename boost::graph_traits<Graph>::halfedge_descriptor halfedge_descriptor;
+  typedef CGAL::Halfedge_around_source_iterator<Graph> out_edge_iterator;
+
+  halfedge_descriptor v0_v1 = halfedge(e,g);
+  halfedge_descriptor v1_v0 = opposite(v0_v1,g);
+
+  vertex_descriptor v0 = target(v1_v0,g), v1 = target(v0_v1,g);
+
+  vertex_descriptor vL = target(next(v0_v1,g),g);
+  vertex_descriptor vR = target(next(v1_v0,g),g);
+
+  out_edge_iterator eb1, ee1 ;
+  out_edge_iterator eb2, ee2 ;
+
+  // The following loop checks the link condition for v0_v1.
+  // Specifically, that for every vertex 'k' adjacent to both 'p and 'q', 'pkq' is a face of the mesh.
+  //
+  for ( boost::tie(eb1,ee1) = halfedges_around_source(v0,g) ;  eb1 != ee1 ; ++ eb1 )
+  {
+    halfedge_descriptor v0_k = *eb1;
+
+    if ( v0_k != v0_v1 )
+    {
+      vertex_descriptor k = target(v0_k,g);
+
+      for ( boost::tie(eb2,ee2) =  halfedges_around_source(k,g) ; eb2 != ee2 ; ++ eb2 )
+      {
+        halfedge_descriptor k_v1 = *eb2;
+
+        if ( target(k_v1,g) == v1 )
+        {
+          // At this point we know p-q-k are connected and we need to determine if this triangle is a face of the mesh.
+          //
+          // Since the mesh is known to be triangular there are at most two faces sharing the edge p-q.
+          //
+          // If p->q is NOT a border edge, the top face is p->q->t where t is target(next(p->q))
+          // If q->p is NOT a border edge, the bottom face is q->p->b where b is target(next(q->p))
+          //
+          // If k is either t or b then p-q-k *might* be a face of the mesh. It won't be if k==t but p->q is border
+          // or k==b but q->b is a border (because in that case even though there exists triangles p->q->t (or q->p->b)
+          // they are holes, not faces)
+          //
+
+          bool lIsFace =   ( vL == k && (! is_border(v0_v1,g)) )
+            || ( vR == k && (! is_border(v1_v0,g)) ) ;
+
+          if ( !lIsFace )
+          {
+            // CGAL_ECMS_TRACE(3,"  k=V" << get(Vertex_index_map,k) << " IS NOT in a face with p-q. NON-COLLAPSABLE edge." ) ;
+            return false ;
+          }
+          else
+          {
+            //CGAL_ECMS_TRACE(4,"  k=V" << get(Vertex_index_map,k) << " is in a face with p-q") ;
+          }
+        }
+      }
+    }
+  }
+
+  // detect isolated triangle (or triangle attached to a mesh with non-manifold vertices)
+  if (!is_border(v0_v1,g) && is_border(opposite(next(v0_v1,g), g), g)
+                          && is_border(opposite(prev(v0_v1,g), g), g) ) return false;
+  if (!is_border(v1_v0,g) && is_border(opposite(next(v1_v0,g), g), g)
+                          && is_border(opposite(prev(v1_v0,g), g), g) ) return false;
+
+  if ( !is_border(v0_v1,g) && !is_border(v1_v0,g) )
+  {
+    if ( is_border(v0,g) && is_border(v1,g) )
+    {
+      //CGAL_ECMS_TRACE(3,"  both p and q are boundary vertices but p-q is not. NON-COLLAPSABLE edge." ) ;
+      return false ;
+    }
+    else
+    {
+      if ( is_tetrahedron(v0_v1,g) )
+      {
+        //CGAL_ECMS_TRACE(3,"  p-q belongs to a tetrahedron. NON-COLLAPSABLE edge." ) ;
+        return false ;
+      }
+      if ( next(v0_v1, g) == opposite(prev(v1_v0, g), g) &&
+           prev(v0_v1, g) == opposite(next(v1_v0, g), g) )
+      {
+        //CGAL_ECMS_TRACE(3,"  degenerate volume." ) ;
+        return false ;
+      }
+    }
+  }
+
+  return true ;
+}
 
 /**
  * collapses an edge in a graph.
@@ -1464,6 +1564,7 @@ collapse_edge(typename boost::graph_traits<Graph>::edge_descriptor e,
   typedef typename Traits::halfedge_descriptor            halfedge_descriptor;
 
   CGAL_precondition(is_valid_edge_descriptor(e, g));
+  CGAL_precondition(does_satisfy_link_condition(e,g));
 
   halfedge_descriptor pq = halfedge(e,g);
   halfedge_descriptor qp = opposite(pq, g);
@@ -1584,6 +1685,7 @@ collapse_edge(typename boost::graph_traits<Graph>::edge_descriptor v0v1,
   typedef typename Traits::halfedge_descriptor            halfedge_descriptor;
 
   CGAL_precondition(is_valid_edge_descriptor(v0v1, g));
+  CGAL_precondition(does_satisfy_link_condition(v0v1,g));
   CGAL_precondition(!get(Edge_is_constrained_map, v0v1));
 
   halfedge_descriptor pq = halfedge(v0v1,g);
@@ -1754,109 +1856,6 @@ flip_edge(typename boost::graph_traits<Graph>::halfedge_descriptor h,
   set_halfedge(foh,oh,g);
 }
 
-/**
- *  \returns `true` if `e` satisfies the *link condition* \cgalCite{degn-tpec-98}, which guarantees that the surface is also 2-manifold after the edge collapse.
- */
-template<typename Graph>
-bool
-does_satisfy_link_condition(typename boost::graph_traits<Graph>::edge_descriptor e,
-                            const Graph& g)
-{
-  typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
-  typedef typename boost::graph_traits<Graph>::halfedge_descriptor halfedge_descriptor;
-  typedef CGAL::Halfedge_around_source_iterator<Graph> out_edge_iterator;
-
-  CGAL_precondition(is_valid_edge_descriptor(e, g));
-
-  halfedge_descriptor v0_v1 = halfedge(e,g);
-  halfedge_descriptor v1_v0 = opposite(v0_v1,g);
-
-  vertex_descriptor v0 = target(v1_v0,g), v1 = target(v0_v1,g);
-
-  vertex_descriptor vL = target(next(v0_v1,g),g);
-  vertex_descriptor vR = target(next(v1_v0,g),g);
-
-  out_edge_iterator eb1, ee1 ;
-  out_edge_iterator eb2, ee2 ;
-
-  // The following loop checks the link condition for v0_v1.
-  // Specifically, that for every vertex 'k' adjacent to both 'p and 'q', 'pkq' is a face of the mesh.
-  //
-  for ( boost::tie(eb1,ee1) = halfedges_around_source(v0,g) ;  eb1 != ee1 ; ++ eb1 )
-  {
-    halfedge_descriptor v0_k = *eb1;
-
-    if ( v0_k != v0_v1 )
-    {
-      vertex_descriptor k = target(v0_k,g);
-
-      for ( boost::tie(eb2,ee2) =  halfedges_around_source(k,g) ; eb2 != ee2 ; ++ eb2 )
-      {
-        halfedge_descriptor k_v1 = *eb2;
-
-        if ( target(k_v1,g) == v1 )
-        {
-          // At this point we know p-q-k are connected and we need to determine if this triangle is a face of the mesh.
-          //
-          // Since the mesh is known to be triangular there are at most two faces sharing the edge p-q.
-          //
-          // If p->q is NOT a border edge, the top face is p->q->t where t is target(next(p->q))
-          // If q->p is NOT a border edge, the bottom face is q->p->b where b is target(next(q->p))
-          //
-          // If k is either t or b then p-q-k *might* be a face of the mesh. It won't be if k==t but p->q is border
-          // or k==b but q->b is a border (because in that case even though there exists triangles p->q->t (or q->p->b)
-          // they are holes, not faces)
-          //
-
-          bool lIsFace =   ( vL == k && (! is_border(v0_v1,g)) )
-            || ( vR == k && (! is_border(v1_v0,g)) ) ;
-
-          if ( !lIsFace )
-          {
-            // CGAL_ECMS_TRACE(3,"  k=V" << get(Vertex_index_map,k) << " IS NOT in a face with p-q. NON-COLLAPSABLE edge." ) ;
-            return false ;
-          }
-          else
-          {
-            //CGAL_ECMS_TRACE(4,"  k=V" << get(Vertex_index_map,k) << " is in a face with p-q") ;
-          }
-        }
-      }
-    }
-  }
-
-  // detect isolated triangle (or triangle attached to a mesh with non-manifold vertices)
-  if (!is_border(v0_v1,g) && is_border(opposite(next(v0_v1,g), g), g)
-                          && is_border(opposite(prev(v0_v1,g), g), g) ) return false;
-  if (!is_border(v1_v0,g) && is_border(opposite(next(v1_v0,g), g), g)
-                          && is_border(opposite(prev(v1_v0,g), g), g) ) return false;
-
-  if ( !is_border(v0_v1,g) && !is_border(v1_v0,g) )
-  {
-    if ( is_border(v0,g) && is_border(v1,g) )
-    {
-      //CGAL_ECMS_TRACE(3,"  both p and q are boundary vertices but p-q is not. NON-COLLAPSABLE edge." ) ;
-      return false ;
-    }
-    else
-    {
-      if ( is_tetrahedron(v0_v1,g) )
-      {
-        //CGAL_ECMS_TRACE(3,"  p-q belongs to a tetrahedron. NON-COLLAPSABLE edge." ) ;
-        return false ;
-      }
-      if ( next(v0_v1, g) == opposite(prev(v1_v0, g), g) &&
-           prev(v0_v1, g) == opposite(next(v1_v0, g), g) )
-      {
-        //CGAL_ECMS_TRACE(3,"  degenerate volume." ) ;
-        return false ;
-      }
-    }
-  }
-
-  return true ;
-}
-
 #ifndef CGAL_NO_DEPRECATED_CODE
 /// \cond SKIP_IN_MANUAL
 template<typename Graph>

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/snap.h

@@ -152,11 +152,43 @@ void simplify_range(HalfedgeRange& halfedge_range,
             new_tolerance += CGAL::approximate_sqrt(CGAL::squared_distance(new_p, pt));
         }
 
-        if (!CGAL::Euler::does_satisfy_link_condition(edge(h, tm), tm))
+        // check that the collapse does not create a new degenerate face
+        bool do_collapse = true;
+        for(halfedge_descriptor he : halfedges_around_target(h, tm))
+        {
+          if(he != h &&
+             !is_border(he, tm) &&
+             collinear(get(vpm, source(he, tm)), new_p, get(vpm, target(next(he,tm),tm))))
+          {
+            do_collapse = false;
+            break;
+          }
+        }
+
+        if(!do_collapse)
           continue;
+
+        for(halfedge_descriptor he : halfedges_around_target(opposite(h,tm), tm))
+        {
+          if(he != opposite(h,tm) &&
+             !is_border(he, tm) &&
+             collinear(get(vpm, source(he, tm)), new_p, get(vpm, target(next(he,tm),tm))))
+          {
+            do_collapse = false;
+            break;
+          }
+        }
+
+        if(!do_collapse)
+          continue;
+
+        if(!CGAL::Euler::does_satisfy_link_condition(edge(h, tm), tm))
+          continue;
+
         const halfedge_descriptor opoh = opposite(prev(opposite(h, tm), tm), tm);
-        if (is_border(opoh, tm))
+        if(is_border(opoh, tm))
-          edges_to_test.erase( opoh );
+          edges_to_test.erase(opoh);
+
         vertex_descriptor v = Euler::collapse_edge(edge(h, tm), tm);
 
         put(vpm, v, new_p);
@@ -164,7 +196,7 @@ void simplify_range(HalfedgeRange& halfedge_range,
 
         if(get(range_halfedges, prev_h))
           edges_to_test.insert(prev_h);
-        if(next_h!=opoh && get(range_halfedges, next_h))
+        if(next_h != opoh && get(range_halfedges, next_h))
           edges_to_test.insert(next_h);
         ++collapsed_n;
       }
@@ -637,7 +669,6 @@ std::size_t split_edges(EdgesToSplitContainer& edges_to_split,
 
   typedef typename boost::property_traits<VPMS>::value_type                       Point;
   typedef typename boost::property_traits<VPMT>::reference                        Point_ref;
-  typedef typename GeomTraits::Vector_3                                           Vector;
 
   typedef std::pair<halfedge_descriptor, Point>                                   Vertex_with_new_position;
   typedef std::vector<Vertex_with_new_position>                                   Vertices_with_new_position;
@@ -687,7 +718,12 @@ std::size_t split_edges(EdgesToSplitContainer& edges_to_split,
 
       bool do_split = true;
 
-      // Some splits can create degenerate faces, avoid that
+      // In case of self-snapping, avoid degenerate caps
+      const bool is_same_mesh = (&tm_T == &tm_S);
+      if(is_same_mesh && target(next(opposite(h_to_split, tm_T), tm_T), tm_T) == splitter_v)
+        do_split = false;
+
+      // Do not split if it would create a degenerate needle
       if((new_position == get(vpm_T, target(h_to_split, tm_T))) ||
          (new_position == get(vpm_T, source(h_to_split, tm_T))))
         do_split = false;
@@ -695,11 +731,58 @@ std::size_t split_edges(EdgesToSplitContainer& edges_to_split,
       if(!first_split && new_position == previous_split_position)
         do_split = false;
 
+      // check if the new faces after split will not be degenerate
+      const Point& p0 = new_position;
+      Point_ref p1 = get(vpm_T, source(h_to_split, tm_T));
+      Point_ref p2 = get(vpm_T, target(next(opposite(h_to_split, tm_T), tm_T), tm_T));
+      Point_ref p3 = get(vpm_T, target(h_to_split, tm_T));
+
+      /* Chooses the diagonal that will split the quad in two triangles that maximizes
+       * the scalar product of the un-normalized normals of the two triangles.
+       *
+       * The lengths of the un-normalized normals (computed using cross-products of two vectors)
+       * are proportional to the area of the triangles.
+       * Maximizing the scalar product of the two normals will avoid skinny triangles,
+       * and will also take into account the cosine of the angle between the two normals.
+       *
+       * In particular, if the two triangles are oriented in different directions,
+       * the scalar product will be negative.
+       */
+      auto p1p3 = CGAL::cross_product(p2-p1, p3-p2) * CGAL::cross_product(p0-p3, p1-p0);
+      auto p0p2 = CGAL::cross_product(p1-p0, p1-p2) * CGAL::cross_product(p3-p2, p3-p0);
+      bool first_split_face = (p0p2 > p1p3);
+
+      if(first_split_face)
+      {
+        if(p0p2 <= 0 || collinear(p0,p1,p2) || collinear(p0,p2,p3))
+          do_split = false;
+      }
+      else
+      {
+        if(p1p3 <= 0 || collinear(p0,p1,p3) || collinear(p1,p2,p3))
+          do_split = false;
+      }
+
+      if(do_split && !is_source_mesh_fixed)
+      {
+        for(halfedge_descriptor h : halfedges_around_target(splitter_v, tm_S))
+        {
+          if(!is_border(h,tm_S) && collinear(get(vpm_S, source(h,tm_S)), new_position, get(vpm_S, target(next(h,tm_S),tm_S))))
+          {
+            do_split = false;
+            break;
+          }
+        }
+
+        if(do_split)
+          put(vpm_S, splitter_v, new_position);
+      }
+
 #ifdef CGAL_PMP_SNAP_DEBUG_PP
       std::cout << " -.-.-. Splitting " << edge(h_to_split, tm_T) << " |||| "
                 << " Vs " << source(h_to_split, tm_T) << " (" << tm_T.point(source(h_to_split, tm_T)) << ")"
                 << " --- Vt " << target(h_to_split, tm_T) << " (" << tm_T.point(target(h_to_split, tm_T)) << ")" << std::endl;
-      std::cout << "With point: " << vnp.second << std::endl;
+      std::cout << "With point: " << new_position << " (init: " << vnp.second << ")" << std::endl;
       std::cout << "Actually split? " << do_split << std::endl;
 #endif
 
@@ -715,75 +798,27 @@ std::size_t split_edges(EdgesToSplitContainer& edges_to_split,
 
         visitor.after_vertex_edge_snap(new_v, tm_T);
       }
-
+      else
-      if(!is_source_mesh_fixed)
+      {
-        put(vpm_S, splitter_v, new_position);
+        continue;
+      }
 
       first_split = false;
       previous_split_position = new_position;
       ++snapped_n;
 
-      // Everything below is choosing the diagonal to triangulate the quad formed by the edge split
+      halfedge_descriptor v0, v1, v2, v3;
-      // So, it's only relevant if splitting has been performed
+      v0 = opposite(h_to_split, tm_T);
-      if(!do_split)
+      v1 = next(v0, tm_T);
-        continue;
+      v2 = next(v1, tm_T);
-
+      v3 = next(v2, tm_T);
-      /*          new_p
+      // halfedge_descriptor new_hd =
-       *         /   \
+      first_split_face ? CGAL::Euler::split_face(v0, v2, tm_T)
-       *    res /     \ h_to_split
+                       : CGAL::Euler::split_face(v1, v3, tm_T);
-       *       /       \
+      if(first_split_face)
-       *      /         \
+        visitor.after_split_face(v0, v2, tm_T);
-       *    left       right
-       *     |         /
-       *     |        /
-       *     |       /
-       *     |      /
-       *     |     /
-       *     |    /
-       *      opp
-       */
-
-      const halfedge_descriptor res = prev(h_to_split, tm_T);
-      const Point_ref left_pt = get(vpm_T, source(res, tm_T));
-      const Point_ref right_pt = get(vpm_T, target(h_to_split, tm_T));
-      const Point_ref opp = get(vpm_T, target(next(opposite(res, tm_T), tm_T), tm_T));
-
-      // Check if 'p' is "visible" from 'opp' (i.e. its projection on the plane 'Pl(left, opp, right)'
-      // falls in the cone with apex 'opp' and sides given by 'left' and 'right')
-      const Vector n = gt.construct_orthogonal_vector_3_object()(right_pt, left_pt, opp);
-      const Point trans_left_pt = gt.construct_translated_point_3_object()(left_pt, n);
-      const Point trans_right_pt = gt.construct_translated_point_3_object()(right_pt, n);
-
-      const Point_ref new_p = get(vpm_T, new_v);
-      const bool left_of_left = (gt.orientation_3_object()(trans_left_pt, left_pt, opp, new_p) == CGAL::POSITIVE);
-      const bool right_of_right = (gt.orientation_3_object()(right_pt, trans_right_pt, opp, new_p) == CGAL::POSITIVE);
-
-      const bool is_visible = (!left_of_left && !right_of_right);
-
-#ifdef CGAL_PMP_SNAP_DEBUG_PP
-      std::cout << "Left/Right: " << left_of_left << " " << right_of_right << std::endl;
-      std::cout << "visible from " << opp << " ? " << is_visible << std::endl;
-#endif
-
-      // h_to_split is equal to 'next(res)' after splitting
-      const halfedge_descriptor h_to_split_opp = opposite(h_to_split, tm_T);
-
-      if(is_visible)
-      {
-        halfedge_descriptor h2 = prev(prev(h_to_split_opp, tm_T), tm_T);
-        halfedge_descriptor new_hd = CGAL::Euler::split_face(h_to_split_opp,
-                                                             h2, tm_T);
-        h_to_split = opposite(prev(new_hd, tm_T), tm_T);
-        visitor.after_split_face(h_to_split_opp, h2, tm_T);
-      }
       else
-      {
+        visitor.after_split_face(v1, v3, tm_T);
-        halfedge_descriptor h2 = prev(h_to_split_opp, tm_T);
-        halfedge_descriptor new_hd = CGAL::Euler::split_face(opposite(res, tm_T),
-                                                             h2, tm_T);
-        h_to_split = opposite(next(new_hd, tm_T), tm_T);
-        visitor.after_split_face(opposite(res, tm_T), h2, tm_T);
-      }
     }
   }
 
@@ -1406,6 +1441,8 @@ std::size_t snap_borders(TriangleMesh& tm,
                                                       true /*self snapping*/, np, np);
 }
 
+//TODO:add an option to preserve orientation?
+
 } // end namespace experimental
 } // end namespace Polygon_mesh_processing
 } // end namespace CGAL

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/snap_vertices.h

@@ -706,6 +706,7 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A,
   typedef typename GetGeomTraits<PolygonMesh, NamedParameters_A>::type         GT;
   typedef typename GT::FT                                                      FT;
   typedef typename boost::property_traits<VPM_B>::value_type                   Point;
+  typedef typename boost::property_traits<VPM_B>::reference                    Point_ref;
 
   typedef std::vector<halfedge_descriptor>                                     Vertex_container;
   typedef std::pair<Vertex_container, FT>                                      Unique_vertex;
@@ -955,8 +956,23 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A,
     {
       if(is_second_mesh_fixed)
       {
+        const Point_ref new_p = get(vpm_B, vb);
         for(const halfedge_descriptor ha : vs_a)
-          put(vpm_A, target(ha, tm_A), get(vpm_B, vb));
+        {
+          bool skip = false;
+          for(halfedge_descriptor haa : halfedges_around_target(ha, tm_A))
+          {
+            if(!is_border(haa,tm_A) &&
+               collinear(get(vpm_A, source(haa,tm_A)), new_p, get(vpm_A, target(next(haa,tm_A),tm_A))))
+            {
+              skip = true;
+              break;
+            }
+          }
+
+          if(!skip)
+            put(vpm_A, target(ha, tm_A), new_p);
+        }
       }
       else
       {
@@ -972,12 +988,40 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A,
 #endif
 
         for(const halfedge_descriptor ha : vs_a)
-          put(vpm_A, target(ha, tm_A), new_p);
+        {
+          bool skip = false;
+          for(halfedge_descriptor haa : halfedges_around_target(ha, tm_A))
+          {
+            if(!is_border(haa,tm_A) &&
+               collinear(get(vpm_A, source(haa,tm_A)), new_p, get(vpm_A, target(next(haa,tm_A),tm_A))))
+            {
+              skip = true;
+              break;
+            }
+          }
+
+          if(!skip)
+            put(vpm_A, target(ha, tm_A), new_p);
+        }
 
         for(const halfedge_descriptor hb : vs_b)
-          put(vpm_B, target(hb, tm_B), new_p);
+        {
-      }
+          bool skip = false;
+          for(halfedge_descriptor hbb : halfedges_around_target(hb, tm_B))
+          {
+            if(!is_border(hbb,tm_B) &&
+               collinear(get(vpm_B, source(hbb,tm_B)), new_p, get(vpm_B, target(next(hbb,tm_B),tm_B))))
+            {
+              skip = true;
+              break;
+            }
+          }
 
+          if(!skip)
+            put(vpm_B, target(hb, tm_B), new_p);
+        }
+      }
+      //TODO: the counter shall depend on skip?
       ++counter;
     }