diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/remesh_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/remesh_impl.h index b1653c1f8cd..53837fb2dc0 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/remesh_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/remesh_impl.h @@ -68,10 +68,8 @@ namespace internal { }; template class Incremental_remesher { @@ -94,9 +92,7 @@ namespace internal { public: Incremental_remesher(PolygonMesh& pmesh - , const FaceRange& face_range , VertexPointMap& vpmap - , const EdgeIsConstrainedMap& ecmap , const bool protect_constraints) : mesh_(pmesh) , vpmap_(vpmap) @@ -120,8 +116,6 @@ namespace internal { Triangle_3(get(vpmap, v1), get(vpmap, v2), get(vpmap, v3))); } tree_ptr_ = new AABB_tree(input_triangles_.begin(), input_triangles_.end()); - - tag_halfedges_status(face_range, ecmap); } ~Incremental_remesher() @@ -130,6 +124,85 @@ namespace internal { delete tree_ptr_; } + template + void init_faces_remeshing(const FaceRange& face_range + , const EdgeIsConstrainedMap& ecmap) + { + tag_halfedges_status(face_range, ecmap); + } + + // split edges of edge_range that have their length > high + template + void split_long_edges(const EdgeRange& edge_range, + const double& high) + { + typedef boost::bimap< + boost::bimaps::set_of, + boost::bimaps::multiset_of > > Boost_bimap; + typedef typename Boost_bimap::value_type long_edge; + + std::cout << "Split long edges (" << high << ")..."; + double sq_high = high*high; + + //collect long edges + Boost_bimap long_edges; + BOOST_FOREACH(edge_descriptor e, edge_range) + { + double sqlen = sqlength(e); + if (sqlen > sq_high) + long_edges.insert(long_edge(halfedge(e, mesh_), sqlen)); + } + + //split long edges + unsigned int nb_splits = 0; + while (!long_edges.empty()) + { + //the edge with longest length + typename Boost_bimap::right_map::iterator eit = long_edges.right.begin(); + halfedge_descriptor he = eit->second; + double sqlen = eit->first; + long_edges.right.erase(eit); + + //split edge + Point refinement_point = this->midpoint(he); + halfedge_descriptor hnew = CGAL::Euler::split_edge(he, mesh_); + CGAL_assertion(he == next(hnew, mesh_)); + ++nb_splits; + + //move refinement point + vertex_descriptor vnew = target(hnew, mesh_); + put(vpmap_, vnew, refinement_point); + + //after splitting + halfedge_descriptor hnew_opp = opposite(hnew, mesh_); + + //check sub-edges + double sqlen_new = 0.25 * sqlen; + if (sqlen_new > sq_high) + { + //if it was more than twice the "long" threshold, insert them + long_edges.insert(long_edge(hnew, sqlen_new)); + long_edges.insert(long_edge(next(hnew, mesh_), sqlen_new)); + } + + //insert new edges to keep triangular faces, and update long_edges + if (!is_border(hnew, mesh_)) + { + halfedge_descriptor hnew2 + = CGAL::Euler::split_face(hnew, next(next(hnew, mesh_), mesh_), mesh_); + } + + //do it again on the other side if we're not on boundary + if (!is_border(hnew_opp, mesh_)) + { + halfedge_descriptor hnew2 + = CGAL::Euler::split_face(prev(hnew_opp, mesh_), next(hnew_opp, mesh_), mesh_); + } + } + std::cout << " done (" << nb_splits << " splits)." << std::endl; + } + // PMP book : // "visits all edges of the mesh //if an edge is longer than the given threshold `high`, the edge @@ -693,6 +766,7 @@ namespace internal { || (is_on_border(hopp) && is_on_patch_border(he)); } + template void tag_halfedges_status(const FaceRange& face_range , const EdgeIsConstrainedMap& ecmap) { 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 a1459f129ef..fb054794d7b 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h @@ -67,6 +67,8 @@ namespace Polygon_mesh_processing { * \cgalParamEnd * \cgalNamedParamsEnd * +* @sa `split_long_edges` +* *@todo we suppose `faces` describe only one patch. Handle several patches. *@todo document `number_of_iterations` */ @@ -78,8 +80,6 @@ void incremental_triangle_based_remeshing(PolygonMesh& pmesh , const double& target_edge_length , const NamedParameters& np) { - std::cout.precision(18); - typedef PolygonMesh PM; using boost::choose_pmap; using boost::get_param; @@ -103,8 +103,9 @@ void incremental_triangle_based_remeshing(PolygonMesh& pmesh bool protect = choose_param(get_param(np, protect_constraints), false); - typename internal::Incremental_remesher - remesher(pmesh, faces, vpmap, ecmap, protect); + typename internal::Incremental_remesher + remesher(pmesh, vpmap, protect); + remesher.init_faces_remeshing(faces, ecmap); unsigned int nb_iterations = choose_param(get_param(np, number_of_iterations), 1); @@ -134,6 +135,48 @@ void incremental_triangle_based_remeshing(PolygonMesh& pmesh parameters::all_default()); } +/*! +* \ingroup PkgPolygonMeshProcessing +* @brief splits the edges listed in `edges` +* +* @tparam PolygonMesh model of `MutableFaceGraph` that +* has an internal property map for `CGAL::vertex_point_t` +*/ +template +void split_long_edges(PolygonMesh& pmesh + , EdgeRange& edge_range + , const double& max_length + , const NamedParameters& np) +{ + typedef PolygonMesh PM; + using boost::choose_pmap; + using boost::get_param; + + typedef typename GetGeomTraits::type GT; + typedef typename GetVertexPointMap::type VPMap; + VPMap vpmap = choose_pmap(get_param(np, boost::vertex_point), + pmesh, + boost::vertex_point); + + typename internal::Incremental_remesher + remesher(pmesh, vpmap, false/*protect constraints*/); + + remesher.split_long_edges(edge_range, max_length); +} + +template +void split_long_edges(PolygonMesh& pmesh + , EdgeRange& edges + , const double& max_length) +{ + split_long_edges(pmesh, + edges, + max_length, + parameters::all_default()); +} + } //end namespace Polygon_mesh_processing } //end namespace CGAL diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/remeshing_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/remeshing_test.cpp index 7104a6f6188..99109836bc3 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/remeshing_test.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/remeshing_test.cpp @@ -6,6 +6,7 @@ #include #include +#include #include #include @@ -103,13 +104,21 @@ int main(int argc, char* argv[]) } const std::set& patch = k_ring(patch_center, 3, m); + std::vector border; + PMP::get_border(m, patch, std::back_inserter(border)); + + CGAL::Polygon_mesh_processing::split_long_edges( + m, border, 1.5 * target_edge_length); + CGAL::Timer t; t.start(); CGAL::Polygon_mesh_processing::incremental_triangle_based_remeshing(m, - faces(m),//patch, + patch, target_edge_length, - CGAL::Polygon_mesh_processing::parameters::number_of_iterations(nb_iter)); + CGAL::Polygon_mesh_processing::parameters::number_of_iterations(nb_iter) + .protect_constraints(true) + ); t.stop(); std::cout << "Remeshing took " << t.time() << std::endl;