Fix Coref Boolean for identical nested patches

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

@@ -743,6 +743,10 @@ class Face_graph_output_builder
     boost::dynamic_bitset<> coplanar_patches_of_tm2_for_union_and_intersection(nb_patches_tm2,false);
     patch_status_not_set_tm1.set();
     patch_status_not_set_tm2.set();
+    // extra containers used when entire connected components are identical (filled only if needed)
+    std::vector<std::size_t> coplanar_tm1_to_coplanar_tm2;
+    std::vector<vertex_descriptor> extreme_vertex_per_cc_1;
+    std::vector<vertex_descriptor> extreme_vertex_per_cc_2;
 
     // first set coplanar status of patches using the coplanar faces collected during the
     // extra intersection edges collected. This is important in the case of full connected components
@@ -1499,8 +1503,61 @@ class Face_graph_output_builder
           {
             if (coplanar_patches_of_tm1.test(patch_id))
             {
-              if (is_tm1_inside_out == is_tm2_inside_out)
+              // Two "identical" coplanar patches that are entire connected components
+              //we have the correspondance between cpln patches thanks to faces in tm1_coplanar_faces and tm2_coplanar_faces
+              CGAL_assertion(tm1_coplanar_faces.size()==tm2_coplanar_faces.size());
+              if (coplanar_tm1_to_coplanar_tm2.empty()) // fill container only once
+              {
+                coplanar_tm1_to_coplanar_tm2.resize(nb_patches_tm1, NID);
+                for (std::size_t i=0; i<tm1_coplanar_faces.size(); ++i)
+                {
+                  std::size_t pid1 = tm1_patch_ids[get(fids1, tm1_coplanar_faces[i])];
+                  std::size_t pid2 = tm2_patch_ids[get(fids2, tm2_coplanar_faces[i])];
+                  coplanar_tm1_to_coplanar_tm2[pid1]=pid2;
+                }
+
+                const vertex_descriptor null_v = boost::graph_traits<TriangleMesh>::null_vertex();
+                extreme_vertex_per_cc_1.assign(nb_patches_tm1, null_v);
+                for (face_descriptor fd : faces(tm1))
+                {
+                  std::size_t patch_id = tm1_patch_ids[get(fids1, fd)];
+                  if (!coplanar_patches_of_tm1.test(patch_id)) continue;
+                  halfedge_descriptor hd=halfedge(fd, tm1);
+                  for (halfedge_descriptor h : CGAL::halfedges_around_face(hd, tm1))
+                  {
+                    vertex_descriptor vd = target(h, tm1);
+                    if (extreme_vertex_per_cc_1[patch_id]==null_v || get(vpm1,extreme_vertex_per_cc_1[patch_id])<get(vpm1,vd))
+                      extreme_vertex_per_cc_1[patch_id]=vd;
+                  }
+                }
+
+                extreme_vertex_per_cc_2.assign(nb_patches_tm2, null_v);
+                for (face_descriptor fd : faces(tm2))
+                {
+                  std::size_t patch_id = tm2_patch_ids[get(fids2, fd)];
+                  if (!coplanar_patches_of_tm2.test(patch_id)) continue;
+                  halfedge_descriptor hd=halfedge(fd, tm2);
+                  for (halfedge_descriptor h : CGAL::halfedges_around_face(hd, tm2))
+                  {
+                    vertex_descriptor vd = target(h, tm2);
+                    if (extreme_vertex_per_cc_2[patch_id]==null_v || get(vpm2,extreme_vertex_per_cc_2[patch_id])<get(vpm2,vd))
+                      extreme_vertex_per_cc_2[patch_id]=vd;
+                  }
+                }
+              }
+
+              const std::size_t patch_id2=coplanar_tm1_to_coplanar_tm2[patch_id];
+              CGAL_assertion(patch_id2!=NID);
+
+              bool is_oo_tm1 = ::CGAL::Polygon_mesh_processing::internal::is_outward_oriented(extreme_vertex_per_cc_1[patch_id], tm1, parameters::vertex_point_map(vpm1)),
+                   is_oo_tm2 = ::CGAL::Polygon_mesh_processing::internal::is_outward_oriented(extreme_vertex_per_cc_2[patch_id2], tm2, parameters::vertex_point_map(vpm2));
+
+              if (is_oo_tm1==is_oo_tm2)
+              {
                 coplanar_patches_of_tm1_for_union_and_intersection.set(patch_id);
+                coplanar_patches_of_tm2_for_union_and_intersection.set(patch_id2);
+                patch_status_not_set_tm2.reset( patch_id2 );
+              }
             }
             else
             {

Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt

@@ -65,6 +65,7 @@ create_single_source_cgal_program("triangulate_hole_with_cdt_2_test.cpp")
 create_single_source_cgal_program("test_pmp_polyhedral_envelope.cpp")
 create_single_source_cgal_program("test_pmp_np_function.cpp")
 create_single_source_cgal_program("test_degenerate_pmp_clip_split_corefine.cpp")
+create_single_source_cgal_program("test_corefinement_cavities.cpp")
 #  create_single_source_cgal_program("test_pmp_repair_self_intersections.cpp")
 
 find_package(Eigen3 3.2.0 QUIET) #(requires 3.2.0 or greater)

Polygon_mesh_processing/test/Polygon_mesh_processing/test_corefinement_cavities.cpp

@@ -0,0 +1,123 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/Polygon_mesh_processing/corefinement.h>
+
+template <class Point, class Mesh>
+void make_hexa(double x, double y, double z,
+               double X, double Y, double Z,
+               Mesh& mesh, int t)
+{
+  CGAL::make_hexahedron(
+    Point(x,y,Z),
+    Point(X,y,Z),
+    Point(X,y,z),
+    Point(x,y,z),
+    Point(x,Y,z),
+    Point(x,Y,Z),
+    Point(X,Y,Z),
+    Point(X,Y,z),
+  mesh);
+
+  using face_descriptor = typename boost::graph_traits<Mesh>::face_descriptor;
+  using halfedge_descriptor = typename boost::graph_traits<Mesh>::halfedge_descriptor;
+
+  std::vector<face_descriptor> fcs(faces(mesh).begin(), faces(mesh).end());
+  for (face_descriptor f : fcs)
+  {
+    halfedge_descriptor h = halfedge(f, mesh);
+    if (t==1) h=next(h,mesh);
+    halfedge_descriptor h2=next(next(h, mesh), mesh);
+    CGAL::Euler::split_face(h, h2, mesh);
+  }
+}
+
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Point_3 = K::Point_3;
+using Mesh = CGAL::Surface_mesh<Point_3>;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+
+void test_operations(Mesh A, Mesh B,
+                     bool reverse_A, bool reverse_B,
+                     std::string round,
+                     std::size_t union_v, std::size_t inter_v, std::size_t diff1_v, std::size_t diff2_v)
+{
+#ifndef VERBOSE
+  CGAL_USE(round);
+#endif
+  if (reverse_A) PMP::reverse_face_orientations(A);
+  if (reverse_B) PMP::reverse_face_orientations(B);
+
+  Mesh out_union, out_inter, out_diff1, out_diff2;
+  std::array<boost::optional<Mesh*>, 4> output;
+  output[PMP::Corefinement::UNION] = &out_union;
+  output[PMP::Corefinement::INTERSECTION] = &out_inter;
+  output[PMP::Corefinement::TM1_MINUS_TM2] = &out_diff1;
+  output[PMP::Corefinement::TM2_MINUS_TM1] = &out_diff2;
+
+  Mesh lA=A, lB=B;
+  PMP::corefine_and_compute_boolean_operations(lA,lB,output);
+#ifdef VERBOSE
+  std::ofstream("out_union_"+round+".off") << out_union;
+  std::ofstream("out_inter_"+round+".off") << out_inter;
+  std::ofstream("out_diff1_"+round+".off") << out_diff1;
+  std::ofstream("out_diff2_"+round+".off") << out_diff2;
+#endif
+  assert(vertices(out_union).size()==union_v);
+  assert(vertices(out_inter).size()==inter_v);
+  assert(vertices(out_diff1).size()==diff1_v);
+  assert(vertices(out_diff2).size()==diff2_v);
+}
+
+int main()
+{
+
+  Mesh A, mh, B;
+  make_hexa<Point_3>(0, 0, 0,
+                     4, 4, 4,
+                     A, 0);
+  make_hexa<Point_3>(1, 1, 1,
+                     2, 2, 2,
+                     mh, 0);
+  make_hexa<Point_3>(1, 1, 1,
+                     2, 2, 2,
+                     B, 1);
+
+  Mesh A2, mh2, B2;
+  make_hexa<Point_3>(5, 0, 0,
+                     9, 4, 4,
+                     A2, 0);
+  make_hexa<Point_3>(6, 1, 1,
+                     7, 2, 2,
+                     mh2, 0);
+  make_hexa<Point_3>(6, 1, 1,
+                     7, 2, 2,
+                     B2, 1);
+
+  A.join(A2);
+  mh.join(mh2);
+  PMP::reverse_face_orientations(mh);
+  A.join(mh);
+  B.join(B2);
+
+#ifdef VERBOSE
+  std::ofstream("A.off") << A;
+  std::ofstream("B.off") << B;
+#endif
+
+  test_operations(A, B, false, false, "r00", 16, 0, 44, 28);
+  test_operations(A, B, false, true,  "r01", 28, 44, 0, 16);
+  test_operations(A, B, true,  false, "r10", 44, 28, 16, 0);
+  test_operations(A, B, true,  true,  "r11", 0, 16, 28, 44);
+
+  test_operations(A, A, false, false, "a00", 32, 32, 0, 0);
+  test_operations(A, A, false, true,  "a01", 0, 0, 32, 32);
+  test_operations(A, A, true,  false, "a10", 0, 0, 32, 32);
+  test_operations(A, A, true,  true,  "a11", 32, 32, 0, 0);
+
+  test_operations(B, B, false, false, "b00", 16, 16, 0, 0);
+  test_operations(B, B, false, true,  "b01", 0, 0, 16, 16);
+  test_operations(B, B, true,  false, "b10", 0, 0, 16, 16);
+  test_operations(B, B, true,  true,  "b11", 16, 16, 0, 0);
+}