Fix tangency of the cut plane with the vertices of the faces

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/clip.h

@@ -1049,7 +1049,7 @@ bool clip(PolygonMesh& pm,
     clip_volume=false;
 
   if (clip_volume && !is_closed(pm)) clip_volume=false;
-  if (clip_volume && !use_compact_clipper) use_compact_clipper=true;
+  if (clip_volume && use_compact_clipper) use_compact_clipper=false;
 
   CGAL_precondition(!clip_volume || !triangulate || does_bound_a_volume(pm));
 

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/refine_with_plane.h

@@ -427,10 +427,19 @@ void refine_with_plane(PolygonMesh& pm,
     for (halfedge_descriptor h : halfedges_around_target(hv, pm))
     {
       if (is_border(h, pm)) continue;
+
+
       Oriented_side prev_ori = get(vertex_os, source(h, pm)),
                     next_ori = get(vertex_os, target(next(h, pm), pm));
-      if (prev_ori == ON_ORIENTED_BOUNDARY || next_ori == ON_ORIENTED_BOUNDARY) continue; // skip full edge
-      if (prev_ori!=next_ori) splitted_faces[face(h, pm)].push_back(h); // skip tangency point
+
+
+      // skip vertices enclosed by vertices on the plane
+      if (prev_ori==ON_ORIENTED_BOUNDARY && next_ori==ON_ORIENTED_BOUNDARY)
+        continue;
+      splitted_faces[face(h, pm)].push_back(h);
+      // we insert tangency point twice as the vertex might be use twice to split a face
+      if (prev_ori==ON_ORIENTED_BOUNDARY || next_ori==ON_ORIENTED_BOUNDARY || prev_ori==next_ori)
+        splitted_faces[face(h, pm)].push_back(h); // skip crossing points
     }
   }
 
@@ -439,14 +448,14 @@ void refine_with_plane(PolygonMesh& pm,
   {
     std::size_t nb_hedges = f_and_hs.second.size();
 
-    CGAL_assertion( nb_hedges%2 ==0 );
-
     visitor.before_subface_creations(f_and_hs.first, pm);
 
     if (nb_hedges==2)
     {
       halfedge_descriptor h1=f_and_hs.second[0], h2=f_and_hs.second[1];
-      CGAL_assertion(next(h1,pm)!=h2 && next(h2,pm)!=h1); // the edge does not already exist
+      if(next(h1,pm)==h2 || next(h2,pm)==h1 || h1==h2) // the edge does not already exist
+        continue;
+
       visitor.before_subface_created(pm);
       halfedge_descriptor res = CGAL::Euler::split_face(h1, h2, pm);
       visitor.after_subface_created(face(h2, pm), pm);
@@ -482,20 +491,55 @@ void refine_with_plane(PolygonMesh& pm,
     else
     {
       // sort hedges to make them match
-      CGAL_assertion(!triangulate);
       // TODO: need mechanism to make it robust even with EPICK
       auto less_hedge = [&pm, vpm](halfedge_descriptor h1, halfedge_descriptor h2)
       {
         return lexicographically_xyz_smaller(get(vpm,target(h1,pm)), get(vpm,target(h2,pm)));
       };
       std::sort(f_and_hs.second.begin(), f_and_hs.second.end(), less_hedge);
 
+      // remove duplicated vertex at the beginning and at the end of the sorted list
+      // in case the next/prev edge is fully included (its the entry/exit point of the line
+      // and we don't need twice the vertex in that case)
+      if (f_and_hs.second[0]==f_and_hs.second[1])
+      {
+        halfedge_descriptor h = f_and_hs.second[0];
+        if ( get(vertex_os, source(h, pm))==ON_ORIENTED_BOUNDARY || get(vertex_os, target(next(h, pm), pm))==ON_ORIENTED_BOUNDARY)
+          f_and_hs.second.erase(f_and_hs.second.begin());
+      }
+      if (f_and_hs.second.back()==*std::prev(f_and_hs.second.end(),2))
+      {
+        halfedge_descriptor h = f_and_hs.second.back();
+        if ( get(vertex_os, source(h, pm))==ON_ORIENTED_BOUNDARY || get(vertex_os, target(next(h, pm), pm))==ON_ORIENTED_BOUNDARY)
+          f_and_hs.second.pop_back();
+      }
+
+
+      nb_hedges = f_and_hs.second.size();
+      CGAL_assertion(nb_hedges%2==0);
+      CGAL_assertion(!triangulate || nb_hedges==2);
+
       for (std::size_t i=0; i<nb_hedges; i+=2)
       {
-        halfedge_descriptor h1=f_and_hs.second[i], h2=f_and_hs.second[i+1];
-        CGAL_assertion(next(h1,pm)!=h2 && next(h2,pm)!=h1); // the edge does not already exist
+        halfedge_descriptor h1=f_and_hs.second[i],
+                            h2=f_and_hs.second[i+1],
+                            h3=i+2<nb_hedges?f_and_hs.second[i+2]:boost::graph_traits<PolygonMesh>::null_halfedge();
+
+        if(next(h1,pm)==h2 || next(h2,pm)==h1 || h1==h2) // the edge does not already exist or is an outer tangency
+          continue;
+
+        bool update_h3 = h2==h3;
+
         visitor.before_subface_created(pm);
         halfedge_descriptor res = CGAL::Euler::split_face(h1, h2, pm);
+
+        if (update_h3 && face(f_and_hs.second[i+3], pm)!=face(h3,pm))
+        {
+          CGAL_assertion(target(h3, pm)==target(res,pm));
+          CGAL_assertion(face(f_and_hs.second[i+3], pm)==face(res,pm));
+          f_and_hs.second[i+2]=res;
+        }
+
         if constexpr (has_visitor)
         {
           if (face(h1,pm)!=face(h2,pm))

Polygon_mesh_processing/test/Polygon_mesh_processing/data-clip/man_torus.off

@@ -0,0 +1,36 @@
+OFF
+16 16 0
+
+-2 -2 -2
+2 -2 -2
+2 2 -2
+-2 2 -2
+-1 -1 -2
+1 -1 -2
+1 1 -2
+-1 1 -2
+-2 -2 2
+2 -2 2
+2 2 2
+-2 2 2
+-1 -1 2
+1 -1 2
+1 1 2
+-1 1 2
+4  1 0 4 5
+4  2 1 5 6
+4  3 2 6 7
+4  0 3 7 4
+4  0 1 9 8
+4  1 2 10 9
+4  2 3 11 10
+4  3 0 8 11
+4  8 9 13 12
+4  9 10 14 13
+4  10 11 15 14
+4  11 8 12 15
+4  5 4 12 13
+4  6 5 13 14
+4  7 6 14 15
+4  4 7 15 12
+

Polygon_mesh_processing/test/Polygon_mesh_processing/data-clip/tgt_case_0.off

@@ -0,0 +1,8 @@
+OFF
+5 1 0
+320 768 0
+256 768 0
+256 640 0
+320 640 0
+288 704 0
+5 0 1 2 3 4

Polygon_mesh_processing/test/Polygon_mesh_processing/data-clip/tgt_case_1.off

@@ -0,0 +1,16 @@
+OFF
+8 4 0
+
+384 768 0
+384 704 0
+240 576 0
+400 576 0
+384 640 0
+192 704 0
+240 464 0
+384 464 0
+3  2 7 3
+5  3 4 1 5 2
+3  2 6 7
+3  1 0 5
+

Polygon_mesh_processing/test/Polygon_mesh_processing/data-clip/tgt_case_2.off

@@ -0,0 +1,16 @@
+OFF
+13 1 0
+192 768 0
+320 768 0
+384 768 0
+384 704 0
+320 704 0
+320 640 0
+384 640 0
+384 576 0
+320 576 0
+320 512 0
+384 512 0
+384 448 0
+192 448 0
+13 0 1 2 3 4 5 6 7 8 9 10 11 12

Polygon_mesh_processing/test/Polygon_mesh_processing/data-clip/tgt_case_3.off

@@ -0,0 +1,20 @@
+OFF
+17 1 0
+320 768 0
+320 704 0
+288 688 0
+320 672 0
+320 624 0
+288 608 0
+320 592 0
+288 544 0
+320 528 0
+320 480 0
+288 480 0
+288 432 0
+320 432 0
+320 384 0
+288 384 0
+128 384 0
+128 768 0
+17 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

Polygon_mesh_processing/test/Polygon_mesh_processing/test_pmp_clip.cpp

@@ -427,7 +427,44 @@ void test()
     PMP::clip(tm1, K::Plane_3(0,-1,0,0), params::use_compact_clipper(false).clip_volume(true));
     assert(faces(tm1).size() == 6);
   }
-
+  {
+    TriangleMesh input;
+    std::ifstream("data-clip/man_torus.off") >> input;
+    PMP::triangulate_faces(input);
+    K::Plane_3 plane(0, -1, 0, -1);
+    PMP::clip(input, plane, CGAL::parameters::clip_volume(true));
+    assert(faces(input).size()==36);
+  }
+  {
+    TriangleMesh input;
+    std::ifstream("data-clip/man_torus.off") >> input;
+    K::Plane_3 plane(0, -1, 0, -1);
+    PMP::clip(input, plane, CGAL::parameters::clip_volume(true).do_not_triangulate_faces(true));
+    assert(faces(input).size()==14);
+  }
+  {
+    TriangleMesh input;
+    std::ifstream("data-clip/man_torus.off") >> input;
+    PMP::triangulate_faces(input);
+    K::Plane_3 plane(0, -1, 0, -2);
+    PMP::clip(input, plane, CGAL::parameters::clip_volume(true));
+    assert(faces(input).size()==32);
+  }
+  {
+    TriangleMesh input;
+    std::ifstream("data-clip/man_torus.off") >> input;
+    K::Plane_3 plane(0, -1, 0, -2);
+    PMP::clip(input, plane, CGAL::parameters::clip_volume(true).do_not_triangulate_faces(true));
+    assert(faces(input).size()==16);
+  }
+  {
+    TriangleMesh input;
+    std::ifstream("data-clip/man_torus.off") >> input;
+    PMP::triangulate_faces(input);
+    K::Plane_3 plane(0, 1, 0, 2);
+    PMP::clip(input, plane, CGAL::parameters::clip_volume(true));
+    assert(faces(input).size()==0);
+  }
   // test special case
   {
     TriangleMesh tm1, tm2;
@@ -910,15 +947,19 @@ void test_new_clip()
     assert(faces(e).size()==5);
     assert(vertices(e).size()==24);
   }
-
   {
     TriangleMesh c;
     std::ifstream("data-clip/c.off") >> c;
     PMP::refine_with_plane(c, K::Plane_3(1,0,0,-2));
     assert(faces(c).size()==2);
     assert(vertices(c).size()==8);
+    c.clear();
+    std::ifstream("data-clip/c.off") >> c;
+    PMP::reverse_face_orientations(c);
+    PMP::refine_with_plane(c, K::Plane_3(1,0,0,-2));
+    assert(faces(c).size()==2);
+    assert(vertices(c).size()==8);
   }
-
   {
     TriangleMesh e;
     std::ifstream("data-clip/ee.off") >> e;
@@ -927,7 +968,6 @@ void test_new_clip()
     assert(faces(e).size()==30);
     assert(vertices(e).size()==28);
   }
-
   {
     TriangleMesh c;
     std::ifstream("data-clip/c.off") >> c;
@@ -936,7 +976,72 @@ void test_new_clip()
     assert(faces(c).size()==8);
     assert(vertices(c).size()==9);
   }
-
+  // tangent case
+  {
+    TriangleMesh c;
+    std::ifstream("data-clip/tgt_case_0.off") >> c;
+    PMP::refine_with_plane(c, K::Plane_3(1,0,0,-288), CGAL::parameters::do_not_triangulate_faces(false));
+    assert(faces(c).size()==3);
+    assert(vertices(c).size()==7);
+    c.clear();
+    std::ifstream("data-clip/tgt_case_0.off") >> c;
+    PMP::reverse_face_orientations(c);
+    PMP::refine_with_plane(c, K::Plane_3(1,0,0,-288), CGAL::parameters::do_not_triangulate_faces(false));
+    assert(faces(c).size()==3);
+    assert(vertices(c).size()==7);
+  }
+  {
+    TriangleMesh c;
+    std::ifstream("data-clip/tgt_case_1.off") >> c;
+    PMP::refine_with_plane(c, K::Plane_3(1,0,0,-384), CGAL::parameters::do_not_triangulate_faces(false));
+    assert(faces(c).size()==6);
+    assert(vertices(c).size()==9);
+    c.clear();
+    std::ifstream("data-clip/tgt_case_1.off") >> c;
+    PMP::reverse_face_orientations(c);
+    PMP::refine_with_plane(c, K::Plane_3(1,0,0,-384), CGAL::parameters::do_not_triangulate_faces(false));
+    assert(faces(c).size()==6);
+    assert(vertices(c).size()==9);
+  }
+  {
+    TriangleMesh c;
+    std::ifstream("data-clip/tgt_case_2.off") >> c;
+    PMP::refine_with_plane(c, K::Plane_3(1,0,0,-320), CGAL::parameters::do_not_triangulate_faces(false));
+    assert(faces(c).size()==4);
+    assert(vertices(c).size()==14);
+    c.clear();
+    std::ifstream("data-clip/tgt_case_2.off") >> c;
+    PMP::reverse_face_orientations(c);
+    PMP::refine_with_plane(c, K::Plane_3(1,0,0,-320), CGAL::parameters::do_not_triangulate_faces(false));
+    assert(faces(c).size()==4);
+    assert(vertices(c).size()==14);
+  }
+  {
+    TriangleMesh c;
+    std::ifstream("data-clip/tgt_case_2.off") >> c;
+    PMP::refine_with_plane(c, K::Plane_3(1,0,0,-384), CGAL::parameters::do_not_triangulate_faces(false));
+    assert(faces(c).size()==1);
+    assert(vertices(c).size()==13);
+    c.clear();
+    std::ifstream("data-clip/tgt_case_2.off") >> c;
+    PMP::reverse_face_orientations(c);
+    PMP::refine_with_plane(c, K::Plane_3(1,0,0,-384), CGAL::parameters::do_not_triangulate_faces(false));
+    assert(faces(c).size()==1);
+    assert(vertices(c).size()==13);
+  }
+  {
+    TriangleMesh c;
+    std::ifstream("data-clip/tgt_case_3.off") >> c;
+    PMP::refine_with_plane(c, K::Plane_3(1,0,0,-288), CGAL::parameters::do_not_triangulate_faces(false));
+    assert(faces(c).size()==6);
+    assert(vertices(c).size()==18);
+    c.clear();
+    std::ifstream("data-clip/tgt_case_3.off") >> c;
+    PMP::reverse_face_orientations(c);
+    PMP::refine_with_plane(c, K::Plane_3(1,0,0,-288), CGAL::parameters::do_not_triangulate_faces(false));
+    assert(faces(c).size()==6);
+    assert(vertices(c).size()==18);
+  }
   {
     TriangleMesh ele;
     std::ifstream(CGAL::data_file_path("meshes/elephant.off")) >> ele;