diff --git a/BGL/include/CGAL/boost/graph/parameters_interface.h b/BGL/include/CGAL/boost/graph/parameters_interface.h index fd75acf47d3..d152be492c8 100644 --- a/BGL/include/CGAL/boost/graph/parameters_interface.h +++ b/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) diff --git a/BGL/test/BGL/test_cgal_bgl_named_params.cpp b/BGL/test/BGL/test_cgal_bgl_named_params.cpp index 91cd3c5ff4a..fd1f631ea6b 100644 --- a/BGL/test/BGL/test_cgal_bgl_named_params.cpp +++ b/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; diff --git a/Polygon_mesh_processing/doc/Polygon_mesh_processing/NamedParameters.txt b/Polygon_mesh_processing/doc/Polygon_mesh_processing/NamedParameters.txt index 4202b10e136..22101de1caf 100644 --- a/Polygon_mesh_processing/doc/Polygon_mesh_processing/NamedParameters.txt +++ b/Polygon_mesh_processing/doc/Polygon_mesh_processing/NamedParameters.txt @@ -159,6 +159,15 @@ during the remeshing process.\n Default: `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 +Type: `bool` \n +Default: `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 Default value is `true` +Type: `bool` \n +Default: `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 Default value is `false` +Type: `bool` \n +Default: `false` +\cgalNPEnd + +\cgalNPBegin{projection_functor} \anchor PMP_projection_functor +Parameter used in `isotropic_remeshing()` to specify an alternative vertex projection method. +\n +Type: Unary function object with vertex descriptor as argument type. \n +Default: A function object projecting vertices on the input surface. \cgalNPEnd \cgalNPTableEnd diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h index acf556b71e3..bfbe181b40e 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h +++ b/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 + template struct Border_constraint_pmap { typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; @@ -125,6 +125,8 @@ namespace internal { : border_edges_ptr(new std::set() ) , pmesh_ptr_(NULL) {} + + template 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& map, + friend bool get(const Border_constraint_pmap& 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& map, + friend void put(Border_constraint_pmap& map, const edge_descriptor& e, const bool is) { @@ -160,21 +162,33 @@ namespace internal { template struct Connected_components_pmap { typedef typename boost::graph_traits::face_descriptor face_descriptor; typedef std::size_t Patch_id; typedef FaceIndexMap FIMap; - typedef EdgeIsConstrainedMap ECMap; - typedef Connected_components_pmap CCMap; + typedef Connected_components_pmap CCMap; typedef CGAL::dynamic_face_property_t Face_property_tag; typedef typename boost::property_map::type Patch_ids_map; Patch_ids_map patch_ids_map; std::size_t nb_cc; + template + bool same_range(const Range& r1, const Range& r2) + { + return boost::begin(r1)==boost::begin(r2) && + boost::end(r1)==boost::end(r2); + } + + template + 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 + 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 - = PMP::connected_components(pmesh, - patch_ids_map, - PMP::parameters::edge_is_constrained_map(ecmap) - .face_index_map(fimap)); - if(nb_cc == 1){ - patch_ids_map = Patch_ids_map(); + if ( same_range(face_range, (faces(pmesh))) ) + { + // applied on the whole mesh + nb_cc + = PMP::connected_components(pmesh, + 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(pmesh, face_range, fimap) ) ) + .face_index_map(fimap)); } } + else + nb_cc=0; // default value } @@ -221,10 +251,12 @@ namespace internal { template + 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::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(), 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 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, @@ -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); - 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_)); + merge_and_update_status(en, ep); if (!mesh_border_case_opp) - { - if (s_ep_p!=PATCH_BORDER) - halfedge_and_opp_removed(prev(he_opp, mesh_)); - else{ - // swap edges so as to keep constained edges: - // replace en_p by epo_p and ep_p by eno_p - ::CGAL::internal::swap_edges(en_p, epo_p, mesh_); - CGAL_assertion(is_valid(mesh_)); - } - } + merge_and_update_status(en_p, ep_p); + + if (!protect_constraints_) + put(ecmap_, e, false); + else + CGAL_assertion( !get(ecmap_, e) ); //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 + 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_)) - && is_flip_allowed(opposite(halfedge(e, mesh_), mesh_)); + bool flip_possible = is_flip_allowed(halfedge(e, 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 void tag_halfedges_status(const FaceRange& face_range) { - //tag MESH, //h and hopp belong to the mesh, not the patch - //tag MESH_BORDER //h belongs to the mesh, face(hopp, pmesh) == null_face() + //init halfedges as: + // - 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 - border_map(mesh_, face_range, fimap_); - //override the border of PATCH - //tag PATCH_BORDER,//h belongs to the patch, hopp doesn't - BOOST_FOREACH(edge_descriptor e, edges(mesh_)) + // tag patch border halfedges + BOOST_FOREACH(halfedge_descriptor h, halfedges(mesh_)) { - if (get(ecmap_, e) - || get(border_map, e) - || get_patch_id(face(halfedge(e, mesh_), mesh_)) - != get_patch_id(face(opposite(halfedge(e, mesh_), mesh_), mesh_))) + if (status(h)==PATCH && status(opposite(h, mesh_))!=PATCH) { - //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; - } + 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; + // update status using constrained edge map + if (!boost::is_same >::value) + { + 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, - const Halfedge_status& s_epo, - const Halfedge_status& s_ep) + void merge_and_update_status(halfedge_descriptor en, + halfedge_descriptor 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) - return; - + Halfedge_status s_ep = status(ep); 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 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 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( @@ -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_; diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h index d689e3db724..d87cd651b6a 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h +++ b/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::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::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::type GT; typedef typename GetVertexPointMap::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//default + internal::No_constraint_pmap//default > ::type ECMap; - ECMap ecmap = (boost::is_same >::value) - //avoid constructing the Border_constraint_pmap if it's not used - ? choose_param(get_param(np, internal_np::edge_is_constrained) - , internal::Border_constraint_pmap(pmesh, faces, fimap)) - : choose_param(get_param(np, internal_np::edge_is_constrained) - , internal::Border_constraint_pmap()); + ECMap ecmap = choose_param(get_param(np, internal_np::edge_is_constrained) + , internal::No_constraint_pmap()); 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()); + 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//default + internal::Connected_components_pmap//default > ::type FPMap; FPMap fpmap = choose_param( get_param(np, internal_np::face_patch), - internal::Connected_components_pmap(pmesh, ecmap, fimap, - boost::is_default_param(get_param(np, internal_np::face_patch)))); + internal::Connected_components_pmap(faces, pmesh, ecmap, fimap, + 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 - 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, - internal::Connected_components_pmap, + internal::Connected_components_pmap, FIMap > remesher(pmesh, vpmap, false/*protect constraints*/ , ecmap , internal::No_constraint_pmap() - , internal::Connected_components_pmap(pmesh, ecmap, fimap, false) + , internal::Connected_components_pmap(faces(pmesh), pmesh, ecmap, fimap, false) , fimap , false/*need aabb_tree*/); diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_plugin.cpp index 12b3e251c11..45ce18f78fc 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_plugin.cpp +++ b/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(1), 4. / 3. * target_length)) { QApplication::restoreOverrideCursor(); diff --git a/Property_map/include/CGAL/property_map.h b/Property_map/include/CGAL/property_map.h index d025f3962e2..f589496dcba 100644 --- a/Property_map/include/CGAL/property_map.h +++ b/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 +OR_property_map +make_OR_property_map(const PM1& pm1, const PM2& pm2) +{ + return OR_property_map(pm1, pm2); +} + // A property map that uses the result of a property map as key. template struct Property_map_binder{