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 83f580237b8..df404a2f9ba 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 @@ -6,7 +6,10 @@ #include #include -#include + +#include +#include +#include namespace CGAL { namespace Polygon_mesh_processing { @@ -34,30 +37,36 @@ namespace internal { void split_long_edges(const double& high) { - //collect long edges + 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; - std::map long_edges; + + //collect long edges + Boost_bimap long_edges; BOOST_FOREACH(edge_descriptor e, edges(mesh_)) { double sqlen = sqlength(e); if(sqlen > sq_high) - long_edges[halfedge(e, mesh_)] = sqlen; + long_edges.insert(long_edge(halfedge(e, mesh_), sqlen)); } //split long edges while (!long_edges.empty()) { - typename std::map::iterator eit = long_edges.begin(); - edge_descriptor e = eit->first; - double sqlen = eit->second; - long_edges.erase(eit); - Point refinement_point = this->midpoint(e); + 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); + Point refinement_point = this->midpoint(he); //split edge - bool is_border = (face(halfedge(e, mesh_), mesh_) == boost::graph_traits::null_face()); - halfedge_descriptor hnew = (!is_border) - ? CGAL::Euler::split_edge(halfedge(e, mesh_), mesh_) - : CGAL::Euler::split_edge(opposite(halfedge(e, mesh_), mesh_), mesh_); + halfedge_descriptor hnew = (!is_border(he, mesh_)) + ? CGAL::Euler::split_edge(he, mesh_) + : CGAL::Euler::split_edge(opposite(he, mesh_), mesh_); vertex_descriptor vnew = target(hnew, mesh_); vpmap_[vnew] = refinement_point; @@ -67,37 +76,146 @@ namespace internal { if (sqlen_new > sq_high) { //if it was more than twice the "long" threshold, insert them - long_edges[hnew] = sqlen_new; - long_edges[next(hnew, mesh_)] = sqlen_new; + 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 (face(hnew, mesh_) != boost::graph_traits::null_face()) + if (!is_border(hnew, mesh_)) { halfedge_descriptor hnew2 = CGAL::Euler::split_face(hnew, next(next(hnew, mesh_), mesh_), mesh_); double sql = sqlength(hnew2); if (sql > sq_high) - long_edges[hnew2] = sql; + long_edges.insert(long_edge(hnew2, sql)); } //do it again on the other side if we're not on boundary halfedge_descriptor hnew_opp = opposite(hnew, mesh_); - if (face(hnew_opp, mesh_) != boost::graph_traits::null_face()) + if (!is_border(hnew_opp, mesh_)) { halfedge_descriptor hnew2 = CGAL::Euler::split_face(prev(hnew_opp, mesh_), next(hnew_opp, mesh_), mesh_); double sql = sqlength(hnew2); if (sql > sq_high) - long_edges[hnew2] = sql; + long_edges.insert(long_edge(hnew2, sql)); } } + std::cout << " done." << std::endl; +#ifdef CGAL_DUMP_REMESHING_STEPS + dump("1-edge_split.off"); +#endif } + void collapse_short_edges(const double& low, const double& high) { - ; + typedef boost::bimap< + boost::bimaps::set_of, + boost::bimaps::multiset_of > > Boost_bimap; + typedef typename Boost_bimap::value_type short_edge; + + std::cout << "Collapse short edges (" << low << ", " << high << ")..."; + double sq_low = low*low; + double sq_high = high*high; + + Boost_bimap short_edges; + BOOST_FOREACH(edge_descriptor e, edges(mesh_)) + { + double sqlen = sqlength(e); + if (sqlen < sq_low) + short_edges.insert(short_edge(halfedge(e, mesh_), sqlen)); + } + + unsigned int nb_collapses = 0; + while (!short_edges.empty()) + { + typename Boost_bimap::right_map::iterator eit = short_edges.right.begin(); + halfedge_descriptor he = eit->second; + double sqlen = eit->first; + short_edges.right.erase(eit); + + vertex_descriptor va = target(he, mesh_); + vertex_descriptor vb = source(he, mesh_); + + //avoid collapsing away from boundary a halfedge incident to boundary + if (is_border(vb, mesh_) || is_border(va, mesh_)) + continue; //temporary + + if (is_border(vb, mesh_)) + { + if (is_border(va, mesh_)) + continue; //happens when he crosses the surface + he = opposite(he, mesh_); + std::swap(va, vb); + } + + if (degree(va, mesh_) < 3 + || degree(vb, mesh_) < 3 + || is_border_edge(he, mesh_) //other collapses could have changed that + || !CGAL::Euler::satisfies_link_condition(he, mesh_))//necessary to collapse + continue; + + //check that collapse would not create an edge with length > high + //iterate on vertices va_i of the one-ring of va + bool collapse_ok = true; + BOOST_FOREACH(halfedge_descriptor ha, halfedges_around_target(va, mesh_)) + { + vertex_descriptor va_i = source(ha, mesh_); + if (sqlength(vb, va_i) > sq_high) + { + collapse_ok = false; + break; + } + } + //if it is allowed, perform the collapse + if (collapse_ok) + { + //"collapse va into vb along e" + // remove edges incident to va and vb, because their lengths will change + BOOST_FOREACH(halfedge_descriptor ha, halfedges_around_target(va, mesh_)) + { + short_edges.left.erase(ha); + short_edges.left.erase(opposite(ha, mesh_)); + } + BOOST_FOREACH(halfedge_descriptor hb, halfedges_around_target(vb, mesh_)) + { + short_edges.left.erase(hb); + short_edges.left.erase(opposite(hb, mesh_)); + } + + CGAL_assertion_code( + halfedge_descriptor en = next(he, mesh_); + halfedge_descriptor enp = next(opposite(he, mesh_), mesh_); + ); + + //perform collapse + Point target_point = vpmap_[vb]; + vertex_descriptor vkept = CGAL::Euler::collapse_edge(edge(he, mesh_), mesh_); + vpmap_[vkept] = target_point; + ++nb_collapses; + + CGAL_assertion(source(en, mesh_) == source(enp, mesh_)); + CGAL_expensive_assertion(is_triangle_mesh(mesh_)); + + //insert new/remaining short edges + BOOST_FOREACH(halfedge_descriptor ht, halfedges_around_target(vkept, mesh_)) + { + double sqlen = sqlength(ht); + if (sqlen < sq_low) + short_edges.insert(short_edge(ht, sqlen)); + } + + std::cout << "."; + std::cout.flush(); + } + } + std::cout << " done (" << nb_collapses << " collapses)." << std::endl; +#ifdef CGAL_DUMP_REMESHING_STEPS + dump("2-edge_collapse.off"); +#endif } + void equalize_valences() { ; @@ -112,12 +230,19 @@ namespace internal { } private: + double sqlength(const vertex_descriptor& v1, + const vertex_descriptor& v2) const + { + return CGAL::squared_distance(vpmap_[v1], vpmap_[v2]); + } + double sqlength(const halfedge_descriptor& h) const { vertex_descriptor v1 = target(h, mesh_); vertex_descriptor v2 = source(h, mesh_); - return CGAL::squared_distance(vpmap_[v1], vpmap_[v2]); + return sqlength(v1, v2); } + double sqlength(const edge_descriptor& e) const { return sqlength(halfedge(e, mesh_)); @@ -130,6 +255,13 @@ namespace internal { return CGAL::midpoint(p1, p2); } + void dump(const char* filename) const + { + std::ofstream out(filename); + out << mesh_; + out.close(); + } + 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 578d714d123..c1414d28c92 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/remesh.h @@ -21,6 +21,8 @@ #ifndef CGAL_POLYGON_MESH_PROCESSING_REMESH_H #define CGAL_POLYGON_MESH_PROCESSING_REMESH_H +#define CGAL_DUMP_REMESHING_STEPS + #include #include @@ -57,7 +59,7 @@ void incremental_triangle_based_remeshing(PolygonMesh& pmesh typename internal::Incremental_remesher remesher(pmesh, vpmap); - unsigned int nb_iterations = choose_param(get_param(np, number_of_iterations), 10); + unsigned int nb_iterations = choose_param(get_param(np, number_of_iterations), 1); double low = 4. / 5. * target_edge_length; double high = 4. / 3. * target_edge_length; 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 d2c0139340c..4fa01cf653f 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/remeshing_test.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/remeshing_test.cpp @@ -25,6 +25,8 @@ void test_edge_lengths(const Mesh& pmesh, vertex_descriptor v1 = target(halfedge(e, pmesh), pmesh); vertex_descriptor v2 = source(halfedge(e, pmesh), pmesh); double sql = CGAL::squared_distance(vpmap[v1], vpmap[v2]); + if (sqhigh < sql) + std::cout << "sqhigh = " << sqhigh << "\t sql = " << sql << std::endl; CGAL_assertion(sqhigh >= sql); } } @@ -39,18 +41,18 @@ int main() return 1; } - double target_edge_length = 0.03; + double target_edge_length = 0.01; double low = 4. / 5. * target_edge_length; double high = 4. / 3. * target_edge_length; CGAL::Polygon_mesh_processing::incremental_triangle_based_remeshing(m, target_edge_length, - CGAL::Polygon_mesh_processing::parameters::number_of_iterations(10)); + CGAL::Polygon_mesh_processing::parameters::number_of_iterations(1)); boost::property_map::const_type vpmap = boost::get(CGAL::vertex_point, m); - test_edge_lengths(m, high, vpmap); +// test_edge_lengths(m, high, vpmap); std::ofstream out("U_remeshed.off"); out << m;