PMP::isotropic_remeshing - add edge_is_ protected_map

PMP_Remeshing/doc/PMP_Remeshing/examples.txt

@@ -6,6 +6,7 @@
 \example PMP_Remeshing/isotropic_remeshing_example.cpp
 \example PMP_Remeshing/isotropic_remeshing_of_patch_example.cpp
 \example PMP_Remeshing/isotropic_remeshing_with_allow_move.cpp
+\example PMP_Remeshing/isotropic_remeshing_with_constrained_edges.cpp
 \example PMP_Remeshing/isotropic_remeshing_with_custom_sizing_example.cpp
 \example PMP_Remeshing/isotropic_remeshing_with_sizing_example.cpp
 \example PMP_Remeshing/mesh_smoothing_example.cpp

PMP_Remeshing/examples/PMP_Remeshing/CMakeLists.txt

@@ -11,6 +11,7 @@ create_single_source_cgal_program("extrude.cpp" )
 create_single_source_cgal_program("hausdorff_distance_remeshing_example.cpp")
 create_single_source_cgal_program("isotropic_remeshing_example.cpp")
 create_single_source_cgal_program("isotropic_remeshing_of_patch_example.cpp")
+create_single_source_cgal_program("isotropic_remeshing_with_constrained_edges.cpp")
 create_single_source_cgal_program("isotropic_remeshing_with_custom_sizing_example.cpp")
 create_single_source_cgal_program("isotropic_remeshing_with_allow_move.cpp")
 create_single_source_cgal_program("random_perturbation_SM_example.cpp")

PMP_Remeshing/examples/PMP_Remeshing/isotropic_remeshing_with_constrained_edges.cpp

@@ -0,0 +1,72 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/Polygon_mesh_processing/detect_features.h>
+#include <CGAL/Polygon_mesh_processing/remesh.h>
+
+#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
+
+#include <iostream>
+#include <string>
+
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Mesh = CGAL::Surface_mesh<K::Point_3>;
+using face_descriptor = boost::graph_traits<Mesh>::face_descriptor;
+
+namespace PMP = CGAL::Polygon_mesh_processing;
+
+struct Border_pmap
+{
+  typedef boost::graph_traits<Mesh>::edge_descriptor key_type;
+  typedef bool value_type;
+  typedef value_type reference;
+  typedef boost::readable_property_map_tag category;
+
+  Border_pmap(const Mesh& m)
+    : m_mesh(m) {}
+
+  inline friend value_type get(const Border_pmap& bmp, const key_type& e)
+  {
+    return is_border(e, bmp.m_mesh);
+  }
+  const Mesh& m_mesh;
+};
+
+int main(int argc, char* argv[])
+{
+  const std::string filename = (argc > 1)
+                             ? argv[1]
+                             : CGAL::data_file_path("meshes/corner_tris_with_hole.off");
+
+  Mesh mesh;
+  if(!PMP::IO::read_polygon_mesh(filename, mesh)) {
+    std::cerr << "Invalid input." << std::endl;
+    return 1;
+  }
+
+  boost::property_map<Mesh, CGAL::edge_is_feature_t>::type
+    eif = get(CGAL::edge_is_feature, mesh);
+  boost::property_map<Mesh, CGAL::face_patch_id_t<int>>::type
+    pid = get(CGAL::face_patch_id_t<int>(), mesh);
+
+  // Run segmentation to fill the edge_is_feature property map and the face_patch_id property map.
+  std::size_t number_of_patches =
+      PMP::sharp_edges_segmentation(mesh, 90, eif, pid);
+  std::cout << "The input mesh has been segmented into " << number_of_patches << " patches." << std::endl;
+
+  // Border edges are protected (i.e. kept identical),
+  // Feature edges are constrained (i.e. can only be collapsed or split).
+  PMP::isotropic_remeshing(faces(mesh),
+                          0.1,
+                          mesh,
+                          PMP::parameters::edge_is_protected_map(Border_pmap(mesh))
+                          .face_patch_map(pid)
+                          .edge_is_constrained_map(eif)
+                          .number_of_iterations(3));
+  std::cout << "Remeshing done." << std::endl;
+
+  CGAL::IO::write_polygon_mesh("out.ply", mesh,
+                               CGAL::parameters::stream_precision(17));
+
+  return EXIT_SUCCESS;
+}

PMP_Remeshing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h

@@ -258,6 +258,7 @@ namespace internal {
          , typename VertexPointMap
          , typename GeomTraits
          , typename EdgeIsConstrainedMap
+         , typename EdgeIsProtectedMap
          , typename VertexIsConstrainedMap
          , typename FacePatchMap
          , typename FaceIndexMap
@@ -279,6 +280,7 @@ namespace internal {
     typedef Incremental_remesher<PM, VertexPointMap
                                , GeomTraits
                                , EdgeIsConstrainedMap
+                               , EdgeIsProtectedMap
                                , VertexIsConstrainedMap
                                , FacePatchMap
                                , FaceIndexMap
@@ -304,6 +306,7 @@ namespace internal {
                        , const GeomTraits& gt
                        , const bool protect_constraints
                        , EdgeIsConstrainedMap ecmap
+                       , EdgeIsProtectedMap epmap
                        , VertexIsConstrainedMap vcmap
                        , FacePatchMap fpmap
                        , FaceIndexMap fimap
@@ -318,8 +321,10 @@ namespace internal {
       , protect_constraints_(protect_constraints)
       , patch_ids_map_(fpmap)
       , ecmap_(ecmap)
+      , epmap_(epmap)
       , vcmap_(vcmap)
       , fimap_(fimap)
+      , feature_edges_pmap_(ecmap, epmap)
     {
       halfedge_status_pmap_ = get(CGAL::dynamic_halfedge_property_t<Halfedge_status>(),
                                   pmesh);
@@ -338,6 +343,24 @@ namespace internal {
       }
     }
 
+    static constexpr bool has_constrained_edges()
+    {
+      using False_pmap = CGAL::Static_boolean_property_map<edge_descriptor, false>;
+      return !std::is_same_v<EdgeIsProtectedMap,   False_pmap>
+          || !std::is_same_v<EdgeIsConstrainedMap, False_pmap>;
+    }
+    bool has_protected_edges()
+    {
+      using False_pmap = CGAL::Static_boolean_property_map<edge_descriptor, false>;
+      return protect_constraints_
+          || !std::is_same_v<EdgeIsProtectedMap, False_pmap>;
+    }
+
+    auto all_feature_edges_pmap() const
+    {
+      return feature_edges_pmap_;
+    }
+
     template<typename FaceRange>
     void init_remeshing(const FaceRange& face_range)
     {
@@ -399,6 +422,8 @@ namespace internal {
         );
       for(edge_descriptor e : edge_range)
       {
+        if(!is_split_allowed(e))
+          continue;
         const halfedge_descriptor he = halfedge(e, mesh_);
         std::optional<double> sqlen = sizing.is_too_long(source(he, mesh_), target(he, mesh_), mesh_);
         if(sqlen != std::nullopt)
@@ -522,7 +547,7 @@ namespace internal {
         std::cout.flush();
 #endif
 
-        if (protect_constraints_ && !is_longest_on_faces(edge(he, mesh_)))
+        if(has_protected_edges() && !is_longest_on_faces(edge(he, mesh_)))
           continue;
 
         //collect patch_ids
@@ -794,6 +819,7 @@ namespace internal {
 
           //perform collapse
           CGAL_assertion(target(halfedge(e, mesh_), mesh_) == vb);
+          CGAL_assertion(!get(epmap_, e));
           vertex_descriptor vkept = CGAL::Euler::collapse_edge(e, mesh_, ecmap_);
           CGAL_assertion(is_valid(mesh_));
           CGAL_assertion(vkept == vb);//is the constrained point still here
@@ -878,15 +904,15 @@ namespace internal {
 
         std::array<halfedge_descriptor, 2> r1 = internal::is_badly_shaped(
             face(he, mesh_),
-            mesh_, vpmap_, vcmap_, ecmap_, gt_,
+            mesh_, vpmap_, vcmap_, feature_edges_pmap_, gt_,
             4, // bound on shortest/longest edge above 4 => needle
             cap_threshold, // bound on the angle: above 160 deg => cap
             0,// collapse length threshold : not needed here
             0); // flip triangle height threshold
 
         std::array<halfedge_descriptor, 2> r2 = internal::is_badly_shaped(
             face(opposite(he, mesh_), mesh_),
-            mesh_, vpmap_, vcmap_, ecmap_, gt_, 4, cap_threshold, 0, 0);
+            mesh_, vpmap_, vcmap_, feature_edges_pmap_, gt_, 4, cap_threshold, 0, 0);
 
         const bool badly_shaped = (r1[0] != boost::graph_traits<PolygonMesh>::null_halfedge()//needle
                                 || r1[1] != boost::graph_traits<PolygonMesh>::null_halfedge()//cap
@@ -912,7 +938,7 @@ namespace internal {
         CGAL_assertion_code(Halfedge_status s1o = status(opposite(he, mesh_)));
 
         CGAL_assertion( is_flip_topologically_allowed(edge(he, mesh_)) );
-        CGAL_assertion( !get(ecmap_, edge(he, mesh_)) );
+        CGAL_assertion(!get(feature_edges_pmap_, edge(he, mesh_)));
         CGAL::Euler::flip_edge(he, mesh_);
 
         if (!badly_shaped)
@@ -1041,7 +1067,7 @@ namespace internal {
       auto constrained_vertices_pmap
         = boost::make_function_property_map<vertex_descriptor>(vertex_constraint);
 
-      if constexpr (std::is_same_v<SizingFunction, Uniform_sizing_field<PM, VertexPointMap>>)
+      if constexpr(std::is_same_v<SizingFunction, Uniform_sizing_field<PM, VertexPointMap>>)
       {
 #ifdef CGAL_PMP_REMESHING_VERBOSE
         std::cout << " using tangential relaxation with weights equal to 1";
@@ -1265,14 +1291,24 @@ namespace internal {
       return is_on_border(e) || is_on_patch_border(e);
     }
 
+    bool is_protected(const edge_descriptor& e) const
+    {
+      return get(epmap_, e) || (protect_constraints_ && is_constrained(e));
+    }
+    bool is_protected(const halfedge_descriptor& h) const
+    {
+      return is_protected(edge(h, mesh_));
+    }
+
     bool is_split_allowed(const edge_descriptor& e) const
     {
+      if(is_protected(e))
+        return false;
+
       halfedge_descriptor h = halfedge(e, mesh_);
       halfedge_descriptor hopp = opposite(h, mesh_);
 
-      if (protect_constraints_ && is_constrained(e))
-        return false;
-      else //allow splitting constraints
+      //allow splitting constraints
       {
         if (is_on_mesh(h) && is_on_mesh(hopp))
           return false;
@@ -1299,7 +1335,10 @@ namespace internal {
       if (is_an_isolated_constraint(he) || is_an_isolated_constraint(hopp))
         return false;
 
-      if ( (protect_constraints_ || !collapse_constraints) && is_constrained(e))
+      if (!collapse_constraints && is_constrained(e))
+        return false;
+
+      if(is_protected(e))
         return false;
       if (is_on_patch(he)) //hopp is also on patch
       {
@@ -1373,6 +1412,9 @@ namespace internal {
 
     bool is_flip_allowed(const edge_descriptor& e) const
     {
+      if(is_protected(e))
+        return false;
+
       bool flip_possible = is_flip_allowed(halfedge(e, mesh_))
                         && is_flip_allowed(opposite(halfedge(e, mesh_), mesh_));
 
@@ -1403,6 +1445,8 @@ namespace internal {
       {
         if (is_on_patch(h))
           continue;
+        else if(is_protected(h))
+          return false;
         else if (is_on_patch_border(h) && relax_constraints)
           continue;
         else
@@ -1595,12 +1639,11 @@ namespace internal {
       }
 
       // update status using constrained edge map
-      if (!std::is_same<EdgeIsConstrainedMap,
-                          Static_boolean_property_map<edge_descriptor, false> >::value)
+      if constexpr (has_constrained_edges())
       {
         for(edge_descriptor e : edges(mesh_))
         {
-          if (get(ecmap_, e))
+          if(get(feature_edges_pmap_, e))
           {
             //deal with h and hopp for borders that are sharp edges to be preserved
             halfedge_descriptor h = halfedge(e, mesh_);
@@ -1763,7 +1806,8 @@ namespace internal {
             short_edges.left.erase(hf);
             short_edges.left.erase(hfo);
 
-            CGAL_assertion( !get(ecmap_, edge(hf, mesh_)) );
+            CGAL_assertion_code(edge_descriptor ef = edge(hf, mesh_));
+            CGAL_assertion( !get(ecmap_, ef) && !get(epmap_, ef));
 
             if (!is_flip_topologically_allowed(edge(hf, mesh_)))
               continue;
@@ -2063,11 +2107,13 @@ namespace internal {
     Triangle_list input_triangles_;
     Patch_id_list input_patch_ids_;
     Halfedge_status_pmap halfedge_status_pmap_;
-    bool protect_constraints_;
+    const bool protect_constraints_;
     FacePatchMap patch_ids_map_;
     EdgeIsConstrainedMap ecmap_;
+    EdgeIsProtectedMap epmap_;
     VertexIsConstrainedMap vcmap_;
     FaceIndexMap fimap_;
+    CGAL::OR_property_map<EdgeIsConstrainedMap, EdgeIsProtectedMap> feature_edges_pmap_;
     CGAL_assertion_code(bool input_mesh_is_valid_;)
 
   };//end class Incremental_remesher