Merge pull request #3069 from sloriot/PMP-isotropic_remeshing_user_projection

Fix constrained edge map update and add user projection functor as input

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

@@ -46,6 +46,7 @@ CGAL_add_named_parameter(sparse_linear_solver_t, sparse_linear_solver, sparse_li
 CGAL_add_named_parameter(number_of_relaxation_steps_t, number_of_relaxation_steps, number_of_relaxation_steps)
 CGAL_add_named_parameter(protect_constraints_t, protect_constraints, protect_constraints)
 CGAL_add_named_parameter(relax_constraints_t, relax_constraints, relax_constraints)
+CGAL_add_named_parameter(collapse_constraints_t, collapse_constraints, collapse_constraints)
 CGAL_add_named_parameter(vertex_is_constrained_t, vertex_is_constrained, vertex_is_constrained_map)
 CGAL_add_named_parameter(face_patch_t, face_patch, face_patch_map)
 CGAL_add_named_parameter(random_uniform_sampling_t, random_uniform_sampling, use_random_uniform_sampling)
@@ -64,6 +65,7 @@ CGAL_add_named_parameter(nb_points_per_distance_unit_t, nb_points_per_distance_u
 CGAL_add_named_parameter(outward_orientation_t, outward_orientation, outward_orientation)
 CGAL_add_named_parameter(overlap_test_t, overlap_test, do_overlap_test_of_bounded_sides)
 CGAL_add_named_parameter(preserve_genus_t, preserve_genus, preserve_genus)
+CGAL_add_named_parameter(projection_functor_t, projection_functor, projection_functor)
 
 // List of named parameters that we use in the package 'Surface Mesh Simplification'
 CGAL_add_named_parameter(get_cost_policy_t, get_cost_policy, get_cost)

BGL/test/BGL/test_cgal_bgl_named_params.cpp

@@ -57,6 +57,7 @@ void test(const NamedParameters& np)
   assert(get_param(np, CGAL::internal_np::number_of_relaxation_steps).v == 16);
   assert(get_param(np, CGAL::internal_np::protect_constraints).v == 17);
   assert(get_param(np, CGAL::internal_np::relax_constraints).v == 18);
+  assert(get_param(np, CGAL::internal_np::collapse_constraints).v == 43);
   assert(get_param(np, CGAL::internal_np::vertex_is_constrained).v == 19);
   assert(get_param(np, CGAL::internal_np::face_patch).v == 20);
   assert(get_param(np, CGAL::internal_np::random_uniform_sampling).v == 21);
@@ -86,6 +87,7 @@ void test(const NamedParameters& np)
   assert(get_param(np, CGAL::internal_np::weight_calculator).v == 39);
   assert(get_param(np, CGAL::internal_np::preserve_genus).v == 40);
   assert(get_param(np, CGAL::internal_np::verbosity_level).v == 41);
+  assert(get_param(np, CGAL::internal_np::projection_functor).v == 42);
 
 
   // Test types
@@ -121,6 +123,7 @@ void test(const NamedParameters& np)
   check_same_type<16>(get_param(np, CGAL::internal_np::number_of_relaxation_steps));
   check_same_type<17>(get_param(np, CGAL::internal_np::protect_constraints));
   check_same_type<18>(get_param(np, CGAL::internal_np::relax_constraints));
+  check_same_type<43>(get_param(np, CGAL::internal_np::collapse_constraints));
   check_same_type<19>(get_param(np, CGAL::internal_np::vertex_is_constrained));
   check_same_type<20>(get_param(np, CGAL::internal_np::face_patch));
   check_same_type<21>(get_param(np, CGAL::internal_np::random_uniform_sampling));
@@ -150,6 +153,7 @@ void test(const NamedParameters& np)
   check_same_type<39>(get_param(np, CGAL::internal_np::weight_calculator));
   check_same_type<40>(get_param(np, CGAL::internal_np::preserve_genus));
   check_same_type<41>(get_param(np, CGAL::internal_np::verbosity_level));
+  check_same_type<42>(get_param(np, CGAL::internal_np::projection_functor));
 }
 
 int main()
@@ -176,6 +180,7 @@ int main()
                          .number_of_relaxation_steps(A<16>(16))
                          .protect_constraints(A<17>(17))
                          .relax_constraints(A<18>(18))
+                         .collapse_constraints(A<43>(43))
                          .vertex_is_constrained_map(A<19>(19))
                          .face_patch_map(A<20>(20))
                          .use_random_uniform_sampling(A<21>(21))
@@ -199,6 +204,7 @@ int main()
                          .weight_calculator(A<39>(39))
                          .preserve_genus(A<40>(40))
                          .verbosity_level(A<41>(41))
+                         .projection_functor(A<42>(42))
        );
 
   return EXIT_SUCCESS;

Polygon_mesh_processing/doc/Polygon_mesh_processing/NamedParameters.txt

@@ -159,6 +159,15 @@ during the remeshing process.\n
 <b>Default:</b> `false`
 \cgalNPEnd
 
+\cgalNPBegin{collapse_constraints} \anchor PMP_collapse_constraints
+enables the collapse of constraints listed by \ref PMP_edge_is_constrained_map
+"edge_is_constrained_map" and boundary edges
+in `isotropic_remeshing()`. If `false`, constraint edges cannot be collapsed
+during the remeshing process.\n
+<b>Type:</b> `bool` \n
+<b>Default:</b> `true`
+\cgalNPEnd
+
 \cgalNPBegin{relax_constraints} \anchor PMP_relax_constraints
 enables the tangential relaxation step in `isotropic_remeshing()`
 to be performed on vertices that are endpoints of constraints listed
@@ -299,15 +308,22 @@ If this parameter is not provided, the perturbation is not deterministic
 \cgalNPBegin{outward_orientation} \anchor PMP_outward_orientation
 Parameter used in orientation functions to choose between an outward or inward orientation.
 \n
-\b Type : `bool` \n
+<b>Type:</b> `bool` \n
-\b Default value is `true`
+<b>Default:</b> `true`
 
 \cgalNPBegin{do_overlap_test_of_bounded_sides} \anchor PMP_do_overlap_test_of_bounded_sides
 Parameter used in intersection test functions to indicate whether overlapping tests of bounded sides
 of close meshes are done in addition to surface intersection tests.
 \n
-\b Type : `bool` \n
+<b>Type:</b> `bool` \n
-\b Default value is `false`
+<b>Default:</b> `false`
+\cgalNPEnd
+
+\cgalNPBegin{projection_functor} \anchor PMP_projection_functor
+Parameter used in `isotropic_remeshing()` to specify an alternative vertex projection method.
+\n
+<b>Type:</b> Unary function object with vertex descriptor as argument type. \n
+<b>Default:</b> A function object projecting vertices on the input surface.
 \cgalNPEnd
 
 \cgalNPTableEnd

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

@@ -105,7 +105,7 @@ namespace internal {
     friend void put(No_constraint_pmap& , const key_type& , const bool ) {}
   };
 
-  template <typename PM, typename FaceRange, typename FaceIndexMap>
+  template <typename PM, typename FaceIndexMap>
   struct Border_constraint_pmap
   {
     typedef typename boost::graph_traits<PM>::halfedge_descriptor halfedge_descriptor;
@@ -125,6 +125,8 @@ namespace internal {
       : border_edges_ptr(new std::set<edge_descriptor>() )
       , pmesh_ptr_(NULL)
     {}
+
+    template <class FaceRange>
     Border_constraint_pmap(const PM& pmesh
                          , const FaceRange& faces
                          , const FIMap& fimap)
@@ -139,14 +141,14 @@ namespace internal {
         border_edges_ptr->insert(edge(h, *pmesh_ptr_));
     }
 
-    friend bool get(const Border_constraint_pmap<PM, FaceRange, FIMap>& map,
+    friend bool get(const Border_constraint_pmap<PM, FIMap>& map,
                     const edge_descriptor& e)
     {
       CGAL_assertion(map.pmesh_ptr_!=NULL);
       return map.border_edges_ptr->count(e)!=0;
     }
 
-    friend void put(Border_constraint_pmap<PM, FaceRange, FIMap>& map,
+    friend void put(Border_constraint_pmap<PM, FIMap>& map,
                     const edge_descriptor& e,
                     const bool is)
     {
@@ -160,21 +162,33 @@ namespace internal {
 
 
   template <typename PM,
-            typename EdgeIsConstrainedMap,
             typename FaceIndexMap>
   struct Connected_components_pmap
   {
     typedef typename boost::graph_traits<PM>::face_descriptor   face_descriptor;
     typedef std::size_t                                         Patch_id;
     typedef FaceIndexMap                                        FIMap;
-    typedef EdgeIsConstrainedMap                                ECMap;
+    typedef Connected_components_pmap<PM, FIMap>                CCMap;
-    typedef Connected_components_pmap<PM, ECMap, FIMap>         CCMap;
 
     typedef CGAL::dynamic_face_property_t<Patch_id> Face_property_tag;
     typedef typename boost::property_map<PM, Face_property_tag >::type Patch_ids_map;
     Patch_ids_map patch_ids_map;
     std::size_t nb_cc;
 
+    template <class Range>
+    bool same_range(const Range& r1, const Range& r2)
+    {
+      return boost::begin(r1)==boost::begin(r2) &&
+             boost::end(r1)==boost::end(r2);
+    }
+
+    template <class Range1, class Range2>
+    bool same_range(const Range1& r1, const Range2& r2)
+    {
+      return std::distance(boost::begin(r1), boost::end(r1)) ==
+             std::distance(boost::begin(r2), boost::end(r2));
+    }
+
   public:
     typedef face_descriptor                     key_type;
     typedef Patch_id                            value_type;
@@ -183,7 +197,9 @@ namespace internal {
 
     //note pmesh is a non-const ref because properties are added and removed
     //modify the mesh data structure, but not the mesh itself
-    Connected_components_pmap(PM& pmesh
+    template <class FaceRange, class EdgeIsConstrainedMap>
+    Connected_components_pmap(const FaceRange& face_range
+                            , PM& pmesh
                             , EdgeIsConstrainedMap ecmap
                             , FIMap fimap
                             , const bool do_init = true)
@@ -194,15 +210,29 @@ namespace internal {
 #ifdef CGAL_PMP_REMESHING_VERBOSE
         std::cout << "Compute connected components property map." << std::endl;
 #endif
-        nb_cc
+        if ( same_range(face_range, (faces(pmesh))) )
-          = PMP::connected_components(pmesh,
+        {
-                                      patch_ids_map,
+          // applied on the whole mesh
-                                      PMP::parameters::edge_is_constrained_map(ecmap)
+          nb_cc
-                                     .face_index_map(fimap));
+            = PMP::connected_components(pmesh,
-        if(nb_cc == 1){
+                                        patch_ids_map,
-          patch_ids_map = Patch_ids_map();
+                                        PMP::parameters::edge_is_constrained_map(ecmap)
+                                       .face_index_map(fimap));
+        }
+        else
+        {
+          // applied on a subset of the mesh
+          nb_cc
+            = PMP::connected_components(pmesh,
+                                        patch_ids_map,
+                                        PMP::parameters::edge_is_constrained_map(
+                                          make_OR_property_map(ecmap
+                                          , internal::Border_constraint_pmap<PM, FIMap>(pmesh, face_range, fimap) ) )
+                                       .face_index_map(fimap));
         }
       }
+      else
+        nb_cc=0; // default value
     }
 
 
@@ -221,10 +251,12 @@ namespace internal {
 
   template<typename PM,
            typename EdgeConstraintMap,
-           typename VertexPointMap>
+           typename VertexPointMap,
+           typename FacePatchMap>
   bool constraints_are_short_enough(const PM& pmesh,
                                     EdgeConstraintMap ecmap,
                                     VertexPointMap vpmap,
+                                    const FacePatchMap& fpm,
                                     const double& high)
   {
     double sqh = high*high;
@@ -232,9 +264,11 @@ namespace internal {
     typedef typename boost::graph_traits<PM>::edge_descriptor     edge_descriptor;
     BOOST_FOREACH(edge_descriptor e, edges(pmesh))
     {
-      if (get(ecmap, e))
+      halfedge_descriptor h = halfedge(e, pmesh);
+      if (  is_border(e, pmesh) ||
+            get(ecmap, e) ||
+            get(fpm, face(h,pmesh))!=get(fpm, face(opposite(h,pmesh),pmesh)) )
       {
-        halfedge_descriptor h = halfedge(e, pmesh);
         if (sqh < CGAL::squared_distance(get(vpmap, source(h, pmesh)),
                                          get(vpmap, target(h, pmesh))))
         {
@@ -311,9 +345,6 @@ namespace internal {
     {
       halfedge_status_pmap_ = get(CGAL::dynamic_halfedge_property_t<Halfedge_status>(),
                                   pmesh);
-      BOOST_FOREACH(halfedge_descriptor h, halfedges(mesh_))
-        put(halfedge_status_pmap_, h, MESH);
-
       CGAL_assertion(CGAL::is_triangle_mesh(mesh_));
     }
 
@@ -368,6 +399,7 @@ namespace internal {
 
 
     // split edges of edge_range that have their length > high
+    // Note: only used to split a range of edges provided as input
     template<typename EdgeRange>
     void split_long_edges(const EdgeRange& edge_range,
                           const double& high)
@@ -389,12 +421,7 @@ namespace internal {
       {
         double sqlen = sqlength(e);
         if (sqlen > sq_high)
-        {
           long_edges.insert(long_edge(halfedge(e, mesh_), sqlen));
-          put(ecmap_, e, false);
-        }
-        else
-          put(ecmap_, e, true);
       }
 
       //split long edges
@@ -410,6 +437,8 @@ namespace internal {
         //split edge
         Point refinement_point = this->midpoint(he);
         halfedge_descriptor hnew = CGAL::Euler::split_edge(he, mesh_);
+        // propagate the constrained status
+        put(ecmap_, edge(hnew, mesh_), get(ecmap_, edge(he, mesh_)));
         CGAL_assertion(he == next(hnew, mesh_));
         ++nb_splits;
 
@@ -428,23 +457,22 @@ namespace internal {
           long_edges.insert(long_edge(hnew, sqlen_new));
           long_edges.insert(long_edge(next(hnew, mesh_), sqlen_new));
         }
-        else
-        {
-          put(ecmap_, edge(hnew, mesh_), true);
-          put(ecmap_, edge(next(hnew, mesh_), mesh_), true);
-        }
 
         //insert new edges to keep triangular faces, and update long_edges
         if (!is_border(hnew, mesh_))
         {
-          CGAL::Euler::split_face(hnew, next(next(hnew, mesh_), mesh_), mesh_);
+          halfedge_descriptor hnew2 =
+            CGAL::Euler::split_face(hnew, next(next(hnew, mesh_), mesh_), mesh_);
+          put(ecmap_, edge(hnew2, mesh_), false);
         }
 
         //do it again on the other side if we're not on boundary
         halfedge_descriptor hnew_opp = opposite(hnew, mesh_);
         if (!is_border(hnew_opp, mesh_))
         {
-          CGAL::Euler::split_face(prev(hnew_opp, mesh_), next(hnew_opp, mesh_), mesh_);
+          halfedge_descriptor hnew2 =
+            CGAL::Euler::split_face(prev(hnew_opp, mesh_), next(hnew_opp, mesh_), mesh_);
+          put(ecmap_, edge(hnew2, mesh_), false);
         }
       }
 #ifdef CGAL_PMP_REMESHING_VERBOSE
@@ -509,6 +537,7 @@ namespace internal {
         Point refinement_point = this->midpoint(he);
         halfedge_descriptor hnew = CGAL::Euler::split_edge(he, mesh_);
         CGAL_assertion(he == next(hnew, mesh_));
+        put(ecmap_, edge(hnew, mesh_), get(ecmap_, edge(he, mesh_)) );
         ++nb_splits;
 
         //move refinement point
@@ -538,6 +567,7 @@ namespace internal {
           halfedge_descriptor hnew2 = CGAL::Euler::split_face(hnew,
                                                               next(next(hnew, mesh_), mesh_),
                                                               mesh_);
+          put(ecmap_, edge(hnew2, mesh_), false);
           Halfedge_status snew = (is_on_patch(hnew) || is_on_patch_border(hnew))
             ? PATCH
             : MESH;
@@ -560,6 +590,7 @@ namespace internal {
           halfedge_descriptor hnew2 = CGAL::Euler::split_face(prev(hnew_opp, mesh_),
                                                               next(hnew_opp, mesh_),
                                                               mesh_);
+          put(ecmap_, edge(hnew2, mesh_), false);
           Halfedge_status snew = (is_on_patch(hnew_opp) || is_on_patch_border(hnew_opp))
              ? PATCH
             : MESH;
@@ -595,7 +626,9 @@ namespace internal {
     // "collapses and thus removes all edges that are shorter than a
     // threshold `low`. [...] testing before each collapse whether the collapse
     // would produce an edge that is longer than `high`"
-    void collapse_short_edges(const double& low, const double& high)
+    void collapse_short_edges(const double& low,
+                              const double& high,
+                              const bool collapse_constraints)
     {
       typedef boost::bimap<
         boost::bimaps::set_of<halfedge_descriptor>,
@@ -617,7 +650,7 @@ namespace internal {
       BOOST_FOREACH(edge_descriptor e, edges(mesh_))
       {
         double sqlen = sqlength(e);
-        if( (sqlen < sq_low) && is_collapse_allowed(e) )
+        if ((sqlen < sq_low) && is_collapse_allowed(e, collapse_constraints))
           short_edges.insert(short_edge(halfedge(e, mesh_), sqlen));
       }
 #ifdef CGAL_PMP_REMESHING_VERBOSE_PROGRESS
@@ -639,7 +672,7 @@ namespace internal {
 #endif
 
         edge_descriptor e = edge(he, mesh_);
-        if (!is_collapse_allowed(e))
+        if (!is_collapse_allowed(e, collapse_constraints))
           continue; //situation could have changed since it was added to the bimap
 
         //handle the boundary case :
@@ -690,7 +723,7 @@ namespace internal {
             continue;//both directions invert a face
         }
         CGAL_assertion(collapse_does_not_invert_face(he));
-        CGAL_assertion(is_collapse_allowed(e));
+        CGAL_assertion(is_collapse_allowed(e, collapse_constraints));
 
         if (!CGAL::Euler::does_satisfy_link_condition(e, mesh_))//necessary to collapse
           continue;
@@ -729,46 +762,32 @@ namespace internal {
           bool mesh_border_case     = is_on_border(he);
           bool mesh_border_case_opp = is_on_border(he_opp);
           halfedge_descriptor ep_p  = prev(he_opp, mesh_);
-          halfedge_descriptor epo_p = opposite(ep_p, mesh_);
           halfedge_descriptor en    = next(he, mesh_);
+          halfedge_descriptor ep    = prev(he, mesh_);
           halfedge_descriptor en_p  = next(he_opp, mesh_);
-          Halfedge_status s_ep_p    = status(ep_p);
-          Halfedge_status s_epo_p   = status(epo_p);
-          Halfedge_status s_ep      = status(prev(he, mesh_));
-          Halfedge_status s_epo     = status(opposite(prev(he, mesh_), mesh_));
 
           // merge halfedge_status to keep the more important on both sides
           //do it before collapse is performed to be sure everything is valid
           if (!mesh_border_case)
-            merge_status(en, s_epo, s_ep);
+            merge_and_update_status(en, ep);
-          if (!mesh_border_case_opp && s_ep_p!=PATCH_BORDER)
-            merge_status(en_p, s_epo_p, s_ep_p);
-
-          halfedge_and_opp_removed(he);
-          if (!mesh_border_case)
-            halfedge_and_opp_removed(prev(he, mesh_));
           if (!mesh_border_case_opp)
-          {
+            merge_and_update_status(en_p, ep_p);
-            if (s_ep_p!=PATCH_BORDER)
+
-              halfedge_and_opp_removed(prev(he_opp, mesh_));
+          if (!protect_constraints_)
-            else{
+            put(ecmap_, e, false);
-              // swap edges so as to keep constained edges:
+          else
-              // replace en_p by epo_p and ep_p by eno_p
+            CGAL_assertion( !get(ecmap_, e) );
-              ::CGAL::internal::swap_edges(en_p, epo_p, mesh_);
-              CGAL_assertion(is_valid(mesh_));
-            }
-          }
 
           //perform collapse
           CGAL_assertion(target(halfedge(e, mesh_), mesh_) == vb);
-          vertex_descriptor vkept = CGAL::Euler::collapse_edge(e, mesh_);
+          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
           ++nb_collapses;
 
           //fix constrained case
           CGAL_assertion((is_constrained(vkept) || is_on_patch_border(vb)) == (is_va_constrained || is_vb_constrained));
-          fix_degenerate_faces(vkept, short_edges, sq_low);
+          fix_degenerate_faces(vkept, short_edges, sq_low, collapse_constraints);
 
 #ifdef CGAL_PMP_REMESHING_DEBUG
           debug_status_map();
@@ -779,7 +798,7 @@ namespace internal {
           BOOST_FOREACH(halfedge_descriptor ht, halfedges_around_target(vkept, mesh_))
           {
             double sqlen = sqlength(ht);
-            if( (sqlen < sq_low) && is_collapse_allowed(edge(ht, mesh_)) )
+            if ((sqlen < sq_low) && is_collapse_allowed(edge(ht, mesh_), collapse_constraints))
               short_edges.insert(short_edge(ht, sqlen));
           }
         }//end if(collapse_ok)
@@ -839,6 +858,8 @@ namespace internal {
 
         Patch_id pid = get_patch_id(face(he, mesh_));
 
+        CGAL_assertion( is_flip_topologically_allowed(edge(he, mesh_)) );
+        CGAL_assertion( !get(ecmap_, edge(he, mesh_)) );
         CGAL::Euler::flip_edge(he, mesh_);
         vva -= 1;
         vvb -= 1;
@@ -876,6 +897,8 @@ namespace internal {
           || !check_normals(target(he, mesh_))
           || !check_normals(source(he, mesh_)))
         {
+          CGAL_assertion( is_flip_topologically_allowed(edge(he, mesh_)) );
+          CGAL_assertion( !get(ecmap_, edge(he, mesh_)) );
           CGAL::Euler::flip_edge(he, mesh_);
           --nb_flips;
 
@@ -1035,7 +1058,7 @@ namespace internal {
 
     // PMP book :
     // "maps the vertices back to the surface"
-    void project_to_surface()
+    void project_to_surface(boost::param_not_found)
     {
       //todo : handle the case of boundary vertices
 #ifdef CGAL_PMP_REMESHING_VERBOSE
@@ -1061,6 +1084,35 @@ namespace internal {
       std::cout << "done." << std::endl;
 #endif
 
+#ifdef CGAL_DUMP_REMESHING_STEPS
+      dump("5-project.off");
+#endif
+    }
+
+    template <class ProjectionFunctor>
+    void project_to_surface(const ProjectionFunctor& proj)
+    {
+      //todo : handle the case of boundary vertices
+#ifdef CGAL_PMP_REMESHING_VERBOSE
+      std::cout << "Project to surface...";
+      std::cout.flush();
+#endif
+      BOOST_FOREACH(vertex_descriptor v, vertices(mesh_))
+      {
+        if (is_constrained(v) || is_isolated(v) || !is_on_patch(v))
+          continue;
+        //note if v is constrained, it has not moved
+        put(vpmap_, v,  proj(v));
+      }
+      CGAL_assertion(is_valid(mesh_));
+      CGAL_assertion(is_triangle_mesh(mesh_));
+#ifdef CGAL_PMP_REMESHING_DEBUG
+      debug_self_intersections();
+#endif
+#ifdef CGAL_PMP_REMESHING_VERBOSE
+      std::cout << "done." << std::endl;
+#endif
+
 #ifdef CGAL_DUMP_REMESHING_STEPS
       dump("5-project.off");
 #endif
@@ -1203,7 +1255,8 @@ private:
       }
     }
 
-    bool is_collapse_allowed(const edge_descriptor& e) const
+    bool is_collapse_allowed(const edge_descriptor& e
+                           , const bool collapse_constraints) const
     {
       halfedge_descriptor he = halfedge(e, mesh_);
       halfedge_descriptor hopp = opposite(he, mesh_);
@@ -1211,7 +1264,7 @@ private:
       if (is_on_mesh(he) && is_on_mesh(hopp))
         return false;
 
-      if (protect_constraints_ && is_constrained(e))
+      if ( (protect_constraints_ || !collapse_constraints) && is_constrained(e))
         return false;
       if (is_on_patch(he)) //hopp is also on patch
       {
@@ -1275,10 +1328,23 @@ private:
       return true;//we already checked we're not pinching a hole in the patch
     }
 
+    bool is_flip_topologically_allowed(const edge_descriptor& e) const
+    {
+      halfedge_descriptor h=halfedge(e, mesh_);
+      return !halfedge(target(next(h, mesh_), mesh_),
+               target(next(opposite(h, mesh_), mesh_), mesh_),
+               mesh_).second;
+    }
+
     bool is_flip_allowed(const edge_descriptor& e) const
     {
-      return is_flip_allowed(halfedge(e, mesh_))
+      bool flip_possible = is_flip_allowed(halfedge(e, mesh_))
-          && is_flip_allowed(opposite(halfedge(e, mesh_), mesh_));
+                        && is_flip_allowed(opposite(halfedge(e, mesh_), mesh_));
+
+      if (!flip_possible) return false;
+
+      // the flip is not possible if the edge already exists
+      return is_flip_topologically_allowed(e);
     }
 
     bool is_flip_allowed(const halfedge_descriptor& h) const
@@ -1411,8 +1477,9 @@ private:
     template<typename FaceRange>
     void tag_halfedges_status(const FaceRange& face_range)
     {
-      //tag MESH,        //h and hopp belong to the mesh, not the patch
+      //init halfedges as:
-      //tag MESH_BORDER  //h belongs to the mesh, face(hopp, pmesh) == null_face()
+      //  - MESH,        //h and hopp belong to the mesh, not the patch
+      //  - MESH_BORDER  //h belongs to the mesh, face(hopp, pmesh) == null_face()
       BOOST_FOREACH(halfedge_descriptor h, halfedges(mesh_))
       {
         //being part of the border of the mesh is predominant
@@ -1434,35 +1501,41 @@ private:
         }
       }
 
-      internal::Border_constraint_pmap<PM, FaceRange, FaceIndexMap>
+      // tag patch border halfedges
-        border_map(mesh_, face_range, fimap_);
+      BOOST_FOREACH(halfedge_descriptor h, halfedges(mesh_))
-      //override the border of PATCH
-      //tag PATCH_BORDER,//h belongs to the patch, hopp doesn't
-      BOOST_FOREACH(edge_descriptor e, edges(mesh_))
       {
-        if (get(ecmap_, e)
+        if (status(h)==PATCH && status(opposite(h, mesh_))!=PATCH)
-          || get(border_map, e)
-          || get_patch_id(face(halfedge(e, mesh_), mesh_))
-              != get_patch_id(face(opposite(halfedge(e, mesh_), mesh_), mesh_)))
         {
-          //deal with h and hopp for borders that are sharp edges to be preserved
+          set_status(h, PATCH_BORDER);
-          halfedge_descriptor h = halfedge(e, mesh_);
+          has_border_ = true;
-          if (status(h) == PATCH){
+        }
-            set_status(h, PATCH_BORDER);
+      }
-            has_border_ = true;
-          }
 
-          halfedge_descriptor hopp = opposite(h, mesh_);
+      // update status using constrained edge map
-          if (status(hopp) == PATCH){
+      if (!boost::is_same<EdgeIsConstrainedMap,
-            set_status(hopp, PATCH_BORDER);
+                          No_constraint_pmap<edge_descriptor> >::value)
-            has_border_ = true;
+      {
+        BOOST_FOREACH(edge_descriptor e, edges(mesh_))
+        {
+          if (get(ecmap_, e))
+          {
+            //deal with h and hopp for borders that are sharp edges to be preserved
+            halfedge_descriptor h = halfedge(e, mesh_);
+            if (status(h) == PATCH){
+              set_status(h, PATCH_BORDER);
+              has_border_ = true;
+            }
+
+            halfedge_descriptor hopp = opposite(h, mesh_);
+            if (status(hopp) == PATCH){
+              set_status(hopp, PATCH_BORDER);
+              has_border_ = true;
+            }
           }
         }
       }
     }
 
-
-    // for a Surface_mesh::Property_map we make MESH the default value
     Halfedge_status status(const halfedge_descriptor& h) const
     {
       return get(halfedge_status_pmap_,h);
@@ -1474,17 +1547,16 @@ private:
       put(halfedge_status_pmap_,h,s);
     }
 
-    void merge_status(const halfedge_descriptor& en,
+    void merge_and_update_status(halfedge_descriptor en,
-                      const Halfedge_status& s_epo,
+                                 halfedge_descriptor ep)
-                      const Halfedge_status& s_ep)
     {
-      //get missing data
+
       halfedge_descriptor eno = opposite(en, mesh_);
+      halfedge_descriptor epo = opposite(ep, mesh_);
       Halfedge_status s_eno = status(eno);
+      Halfedge_status s_epo = status(epo);
 
-      if (s_epo == MESH_BORDER && s_eno == MESH_BORDER)
+      Halfedge_status s_ep = status(ep);
-        return;
-
       if(s_epo == MESH_BORDER
         || s_ep == MESH_BORDER
         || s_epo == PATCH_BORDER
@@ -1493,13 +1565,26 @@ private:
         set_status(en, s_epo);
         set_status(eno, s_ep);
       }
+      else
+      {
+        Halfedge_status s_en = status(en);
+        if(s_eno == MESH_BORDER
+          || s_en == MESH_BORDER
+          || s_eno == PATCH_BORDER
+          || s_en == PATCH_BORDER)
+        {
+          set_status(ep, s_epo);
+          set_status(epo, s_ep);
+        }
+      }
       // else keep current status for en and eno
     }
 
     template<typename Bimap>
     void fix_degenerate_faces(const vertex_descriptor& v,
                               Bimap& short_edges,
-                              const double& sq_low)
+                              const double& sq_low,
+                              const bool collapse_constraints)
     {
       CGAL_assertion_code(std::size_t nb_done = 0);
       boost::unordered_set<halfedge_descriptor> degenerate_faces;
@@ -1543,7 +1628,8 @@ private:
 
             short_edges.left.erase(hf);
             short_edges.left.erase(hfo);
-
+            CGAL_assertion( is_flip_topologically_allowed(edge(hf, mesh_)) );
+            CGAL_assertion( !get(ecmap_, edge(hf, mesh_)) );
             CGAL::Euler::flip_edge(hf, mesh_);
             CGAL_assertion_code(++nb_done);
 
@@ -1560,7 +1646,7 @@ private:
 #endif
 
             //insert new edges in 'short_edges'
-            if (is_collapse_allowed(edge(hf, mesh_)))
+            if (is_collapse_allowed(edge(hf, mesh_), collapse_constraints))
             {
               double sqlen = sqlength(hf);
               if (sqlen < sq_low)
@@ -1703,12 +1789,6 @@ private:
         set_status(h, s);
     }
 
-    void halfedge_and_opp_removed(const halfedge_descriptor& h)
-    {
-      set_status(h,MESH);
-      set_status(opposite(h, mesh_),MESH);
-    }
-
     std::size_t nb_valid_halfedges() const
     {
       return static_cast<std::size_t>(
@@ -1833,18 +1913,6 @@ private:
       return input_patch_ids_;
     }
 
-    void update_constraints_property_map()
-    {
-      BOOST_FOREACH(edge_descriptor e, edges(mesh_))
-      {
-        if (is_on_patch_border(halfedge(e, mesh_))
-          || is_on_patch_border(opposite(halfedge(e, mesh_), mesh_)))
-          put(ecmap_, e, true);
-        else
-          put(ecmap_, e, false);
-      }
-    }
-
   private:
     PolygonMesh& mesh_;
     VertexPointMap& vpmap_;

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h

@@ -81,6 +81,7 @@ namespace Polygon_mesh_processing {
 *  \cgalParamBegin{edge_is_constrained_map} a property map containing the
 *    constrained-or-not status of each edge of `pmesh`. A constrained edge can be split
 *    or collapsed, but not flipped, nor its endpoints moved by smoothing.
+*    Sub-edges generated by splitting are set to be constrained.
 *    Note that patch boundary edges (i.e. incident to only one face in the range)
 *    are always considered as constrained edges.
 *  \cgalParamEnd
@@ -96,6 +97,10 @@ namespace Polygon_mesh_processing {
 *     good quality results. It can even fail to terminate because of cascading vertex
 *     insertions.
 *  \cgalParamEnd
+*  \cgalParamBegin{collapse_constraints} If `true`, the edges set as constrained
+*     in `edge_is_constrained_map` (or by default the boundary edges)
+*     are collapsed during remeshing. This value is ignored if `protect_constraints` is true;
+*  \cgalParamEnd
 *  \cgalParamBegin{face_patch_map} a property map with the patch id's associated to the
      faces of `faces`. Instance of a class model of `ReadWritePropertyMap`. It gets
      updated during the remeshing process while new faces are created.
@@ -107,11 +112,19 @@ namespace Polygon_mesh_processing {
 *    constrained in `edge_is_constrained_map` and boundary edges move along the
 *    constrained polylines they belong to.
 *  \cgalParamEnd
+*  \cgalParamBegin{do_project} a boolean that sets whether vertices should be reprojected
+*    on the input surface after creation or displacement.
+*  \cgalParamEnd
+*  \cgalParamBegin{projection_functor}
+*  A function object used to project input vertices (moved by the smoothing) and created vertices.
+*  It must have `%Point_3 operator()(vertex_descriptor)`, `%Point_3` being the value type
+*  of the vertex point map.
+*  If not provided, vertices are projected on the input surface mesh.
+*  \cgalParamEnd
 * \cgalNamedParamsEnd
 *
 * @sa `split_long_edges()`
 *
-*@todo add possibility to provide a functor that projects to a prescribed surface
 *@todo Deal with exact constructions Kernel. The only thing that makes sense is to
 *      guarantee that the output vertices are exactly on the input surface.
 *      To do so, we can do every construction in `double`, and use an exact process for
@@ -134,6 +147,7 @@ void isotropic_remeshing(const FaceRange& faces
 
   typedef PolygonMesh PM;
   typedef typename boost::graph_traits<PM>::vertex_descriptor vertex_descriptor;
+  typedef typename boost::graph_traits<PM>::edge_descriptor edge_descriptor;
   using boost::get_param;
   using boost::choose_param;
 
@@ -145,6 +159,9 @@ void isotropic_remeshing(const FaceRange& faces
   t.start();
 #endif
 
+  static const bool need_aabb_tree =
+    boost::is_default_param(get_param(np, internal_np::projection_functor));
+
   typedef typename GetGeomTraits<PM, NamedParameters>::type GT;
 
   typedef typename GetVertexPointMap<PM, NamedParameters>::type VPMap;
@@ -158,14 +175,10 @@ void isotropic_remeshing(const FaceRange& faces
   typedef typename boost::lookup_named_param_def <
       internal_np::edge_is_constrained_t,
       NamedParameters,
-      internal::Border_constraint_pmap<PM, FaceRange, FIMap>//default
+      internal::No_constraint_pmap<edge_descriptor>//default
     > ::type ECMap;
-  ECMap ecmap = (boost::is_same<ECMap, internal::Border_constraint_pmap<PM, FaceRange, FIMap> >::value)
+  ECMap ecmap = choose_param(get_param(np, internal_np::edge_is_constrained)
-     //avoid constructing the Border_constraint_pmap if it's not used
+                            , internal::No_constraint_pmap<edge_descriptor>());
-    ? choose_param(get_param(np, internal_np::edge_is_constrained)
-                 , internal::Border_constraint_pmap<PM, FaceRange, FIMap>(pmesh, faces, fimap))
-    : choose_param(get_param(np, internal_np::edge_is_constrained)
-                 , internal::Border_constraint_pmap<PM, FaceRange, FIMap>());
 
   typedef typename boost::lookup_named_param_def <
       internal_np::vertex_is_constrained_t,
@@ -175,29 +188,35 @@ void isotropic_remeshing(const FaceRange& faces
   VCMap vcmap = choose_param(get_param(np, internal_np::vertex_is_constrained),
                              internal::No_constraint_pmap<vertex_descriptor>());
 
+  bool protect = choose_param(get_param(np, internal_np::protect_constraints), false);
   typedef typename boost::lookup_named_param_def <
       internal_np::face_patch_t,
       NamedParameters,
-      internal::Connected_components_pmap<PM, ECMap, FIMap>//default
+      internal::Connected_components_pmap<PM, FIMap>//default
     > ::type FPMap;
   FPMap fpmap = choose_param(
     get_param(np, internal_np::face_patch),
-    internal::Connected_components_pmap<PM, ECMap, FIMap>(pmesh, ecmap, fimap,
+    internal::Connected_components_pmap<PM, FIMap>(faces, pmesh, ecmap, fimap,
-        boost::is_default_param(get_param(np, internal_np::face_patch))));
+      boost::is_default_param(get_param(np, internal_np::face_patch)) && (need_aabb_tree
+#if !defined(CGAL_NO_PRECONDITIONS)
+      || protect // face patch map is used to identify patch border edges to check protected edges are short enough
+#endif
+    ) ) );
 
   double low = 4. / 5. * target_edge_length;
   double high = 4. / 3. * target_edge_length;
 
-  bool protect = choose_param(get_param(np, internal_np::protect_constraints), false);
+#if !defined(CGAL_NO_PRECONDITIONS)
   if(protect)
   {
     std::string msg("Isotropic remeshing : protect_constraints cannot be set to");
     msg.append(" true with constraints larger than 4/3 * target_edge_length.");
     msg.append(" Remeshing aborted.");
     CGAL_precondition_msg(
-      internal::constraints_are_short_enough(pmesh, ecmap, vpmap, high),
+      internal::constraints_are_short_enough(pmesh, ecmap, vpmap, fpmap, high),
       msg.c_str());
   }
+#endif
 
 #ifdef CGAL_PMP_REMESHING_VERBOSE
   t.stop();
@@ -208,7 +227,7 @@ void isotropic_remeshing(const FaceRange& faces
 #endif
 
   typename internal::Incremental_remesher<PM, VPMap, GT, ECMap, VCMap, FPMap, FIMap>
-    remesher(pmesh, vpmap, protect, ecmap, vcmap, fpmap, fimap);
+    remesher(pmesh, vpmap, protect, ecmap, vcmap, fpmap, fimap, need_aabb_tree);
   remesher.init_remeshing(faces);
 
 #ifdef CGAL_PMP_REMESHING_VERBOSE
@@ -216,6 +235,7 @@ void isotropic_remeshing(const FaceRange& faces
   std::cout << " done ("<< t.time() <<" sec)." << std::endl;
 #endif
 
+  bool collapse_constraints = choose_param(get_param(np, internal_np::collapse_constraints), true);
   unsigned int nb_iterations = choose_param(get_param(np, internal_np::number_of_iterations), 1);
   bool smoothing_1d = choose_param(get_param(np, internal_np::relax_constraints), false);
   unsigned int nb_laplacian = choose_param(get_param(np, internal_np::number_of_relaxation_steps), 1);
@@ -235,19 +255,17 @@ void isotropic_remeshing(const FaceRange& faces
     if (target_edge_length>0)
     {
       remesher.split_long_edges(high);
-      remesher.collapse_short_edges(low, high);
+      remesher.collapse_short_edges(low, high, collapse_constraints);
     }
     remesher.equalize_valences();
     remesher.tangential_relaxation(smoothing_1d, nb_laplacian);
-    remesher.project_to_surface();
+    if ( choose_param(get_param(np, internal_np::do_project), true) )
-
+      remesher.project_to_surface(get_param(np, internal_np::projection_functor));
 #ifdef CGAL_PMP_REMESHING_VERBOSE
     std::cout << std::endl;
 #endif
   }
 
-  remesher.update_constraints_property_map();
-
 #ifdef CGAL_PMP_REMESHING_VERBOSE
   t.stop();
   std::cout << "Remeshing done (size = " << target_edge_length;
@@ -296,6 +314,12 @@ void isotropic_remeshing(
 *  \cgalParamBegin{vertex_point_map} the property map with the points associated
 *    to the vertices of `pmesh`. Instance of a class model of `ReadWritePropertyMap`.
 *  \cgalParamEnd
+*  \cgalParamBegin{face_index_map} a property map containing the index of each face of `pmesh`
+*  \cgalParamEnd
+*  \cgalParamBegin{edge_is_constrained_map} a property map containing the
+*    constrained-or-not status of each edge of `pmesh`. A constrained edge can be split,
+*    and the sub-edges are set to be constrained.
+*  \cgalParamEnd
 * \cgalNamedParamsEnd
 *
 * @sa `isotropic_remeshing()`
@@ -334,13 +358,13 @@ void split_long_edges(const EdgeRange& edges
   
   typename internal::Incremental_remesher<PM, VPMap, GT, ECMap,
     internal::No_constraint_pmap<vertex_descriptor>,
-    internal::Connected_components_pmap<PM, ECMap, FIMap>,
+    internal::Connected_components_pmap<PM, FIMap>,
     FIMap
   >
     remesher(pmesh, vpmap, false/*protect constraints*/
            , ecmap
            , internal::No_constraint_pmap<vertex_descriptor>()
-           , internal::Connected_components_pmap<PM, ECMap, FIMap>(pmesh, ecmap, fimap, false)
+           , internal::Connected_components_pmap<PM, FIMap>(faces(pmesh), pmesh, ecmap, fimap, false)
            , fimap
            , false/*need aabb_tree*/);
 

Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_plugin.cpp

@@ -400,6 +400,7 @@ public Q_SLOTS:
                  *selection_item->polyhedron(),
                  selection_item->constrained_edges_pmap(),
                  get(CGAL::vertex_point, *selection_item->polyhedron()),
+                 CGAL::Static_property_map<face_descriptor, std::size_t>(1),
                  4. / 3. * target_length))
             {
               QApplication::restoreOverrideCursor();

Property_map/include/CGAL/property_map.h

@@ -85,6 +85,8 @@ class OR_property_map {
   PM2 pm2;
 
  public:
+  OR_property_map() {} // required by boost::connected_components
+
   OR_property_map(PM1 pm1, PM2 pm2)
     : pm1(pm1),pm2(pm2)
   {}
@@ -103,9 +105,15 @@ class OR_property_map {
     put(pm.pm1,k, v);
     put(pm.pm2,k, v);
   }
-
 };
 
+template <class PM1, class PM2>
+OR_property_map<PM1, PM2>
+make_OR_property_map(const PM1& pm1, const PM2& pm2)
+{
+  return OR_property_map<PM1, PM2>(pm1, pm2);
+}
+
 // A property map that uses the result of a property map as key.
 template <class KeyMap, class ValueMap>
 struct Property_map_binder{