diff --git a/BGL/doc/BGL/Concepts/HalfedgeGraph.h b/BGL/doc/BGL/Concepts/HalfedgeGraph.h index 5a5fdfa2dfe..19844256407 100644 --- a/BGL/doc/BGL/Concepts/HalfedgeGraph.h +++ b/BGL/doc/BGL/Concepts/HalfedgeGraph.h @@ -49,7 +49,7 @@ class HalfedgeGraph { }; /*! \relates HalfedgeGraph -returns the edge corresponding to halfedges `h` and `opposite(h,g)`. +returns the edge corresponding to halfedges `h` and `opposite(h,g)`, with the following invariant `halfedge(edge(h,g),g)==h`. */ template boost::graph_traits::edge_descriptor diff --git a/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h b/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h index 8f008b0ffb6..e82007d8f85 100644 --- a/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h +++ b/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_3.h @@ -67,8 +67,10 @@ namespace CGAL typedef typename Polyhedron_dual::Vertex_const_iterator Vertex_const_iterator; - // typedefs for primal + // typedef and type for primal typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typename boost::property_map::type vpm = + get(CGAL::vertex_point, primal); // Typedefs for intersection typedef typename Kernel::Plane_3 Plane_3; @@ -115,8 +117,9 @@ namespace CGAL origin.y() + pp->y(), origin.z() + pp->z()); - - primal_vertices[it] = add_vertex(ppp, primal); + vertex_descriptor vd = add_vertex(primal); + primal_vertices[it] = vd; + put(vpm, vd, ppp); } // Then, add facets to the primal polyhedron diff --git a/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_with_constructions_3.h b/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_with_constructions_3.h index 5ec0cbd16cd..45cb470effb 100644 --- a/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_with_constructions_3.h +++ b/Convex_hull_3/include/CGAL/Convex_hull_3/dual/halfspace_intersection_with_constructions_3.h @@ -56,22 +56,25 @@ namespace CGAL typedef typename boost::graph_traits::face_descriptor face_descriptor; typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typename boost::property_map::const_type vpmap = get(CGAL::vertex_point, primal); + typename boost::property_map::const_type vpm_primal = get(CGAL::vertex_point, primal); + typename boost::property_map::type vpm_dual = get(CGAL::vertex_point, dual); // compute coordinates of extreme vertices in the dual polyhedron // from primal faces boost::unordered_map extreme_points; BOOST_FOREACH (face_descriptor fd , faces( primal)){ halfedge_descriptor h = halfedge(fd,primal); - Plane_3 p (get(vpmap, target(h, primal)), - get(vpmap, target(next(h, primal), primal)), - get(vpmap, target(next(next(h, primal), primal), primal))); + Plane_3 p (get(vpm_primal, target(h, primal)), + get(vpm_primal, target(next(h, primal), primal)), + get(vpm_primal, target(next(next(h, primal), primal), primal))); // translate extreme vertex Point_3 extreme_p = CGAL::ORIGIN + p.orthogonal_vector () / (-p.d()); Point_3 translated_extreme_p(extreme_p.x() + origin.x(), extreme_p.y() + origin.y(), extreme_p.z() + origin.z()); - extreme_points[fd] = add_vertex(translated_extreme_p,dual); + vertex_descriptor vd = add_vertex(dual); + extreme_points[fd] = vd; + put(vpm_dual, vd, translated_extreme_p); } // build faces diff --git a/Convex_hull_3/include/CGAL/convex_hull_3.h b/Convex_hull_3/include/CGAL/convex_hull_3.h index 14b0d0fde44..ea8771f3ef4 100644 --- a/Convex_hull_3/include/CGAL/convex_hull_3.h +++ b/Convex_hull_3/include/CGAL/convex_hull_3.h @@ -316,13 +316,26 @@ void coplanar_3_hull(InputIterator first, InputIterator beyond, typename boost::property_map::type vpm = get(CGAL::vertex_point, P); - std::vector::vertex_descriptor> vertices; + typedef boost::graph_traits Graph_traits; + typedef typename Graph_traits::vertex_descriptor vertex_descriptor; + typedef typename Graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename Graph_traits::face_descriptor face_descriptor; + std::vector vertices; vertices.reserve(CH_2.size()); BOOST_FOREACH(const Point_3& p, CH_2){ vertices.push_back(add_vertex(P)); put(vpm, vertices.back(),p); } - Euler::add_face(vertices, P); + face_descriptor f = Euler::add_face(vertices, P); + + // Then triangulate that face + const halfedge_descriptor he = halfedge(f, P); + halfedge_descriptor other_he = next(next(he, P), P); + for(std::size_t i = 3, end = vertices.size(); i < end; ++i) { + const halfedge_descriptor next_he = next(other_he, P); + Euler::split_face(other_he, he, P); + other_he = next_he; + } } diff --git a/Convex_hull_3/test/Convex_hull_3/quickhull_test_3.cpp b/Convex_hull_3/test/Convex_hull_3/quickhull_test_3.cpp index c5c4f22a255..61756f28d79 100644 --- a/Convex_hull_3/test/Convex_hull_3/quickhull_test_3.cpp +++ b/Convex_hull_3/test/Convex_hull_3/quickhull_test_3.cpp @@ -106,6 +106,22 @@ void test_small_hull() polyhedron2.size_of_facets() == 6 ); } +void test_coplanar_hull() +{ + std::vector points; + points.push_back(Point_3(0,0,0)); + points.push_back(Point_3(1,0,0)); + points.push_back(Point_3(0,1,0)); + points.push_back(Point_3(1,1,0)); + points.push_back(Point_3(1.5,0.5,0)); + points.push_back(Point_3(0.5,1.1,0)); + Polyhedron_3 polyhedron; + CGAL::convex_hull_3(points.begin(), points.end(), polyhedron, Traits()); + assert( polyhedron.is_valid() ); + assert( polyhedron.size_of_facets() == 4 ); + assert( polyhedron.is_pure_triangle() ); +} + int main() { @@ -115,6 +131,8 @@ int main() test_tetrahedron_convexity(); std::cerr << "Testing small hull" << std::endl; test_small_hull(); + std::cerr << "Testing coplanar hull" << std::endl; + test_coplanar_hull(); std::cerr << "Testing 500 random points" << std::endl; std::vector points; diff --git a/Mesh_2/doc/Mesh_2/CGAL/lloyd_optimize_mesh_2.h b/Mesh_2/doc/Mesh_2/CGAL/lloyd_optimize_mesh_2.h index 329d0b12b5b..e9cbce4e5d5 100644 --- a/Mesh_2/doc/Mesh_2/CGAL/lloyd_optimize_mesh_2.h +++ b/Mesh_2/doc/Mesh_2/CGAL/lloyd_optimize_mesh_2.h @@ -101,13 +101,13 @@ lloyd_optimize_mesh_2(cdt, \endcode -\sa `Mesh_optimization_return_code` +\sa `CGAL::Mesh_optimization_return_code` \sa `CGAL::refine_Delaunay_mesh_2()` */ template -Mesh_optimization_return_code +CGAL::Mesh_optimization_return_code lloyd_optimize_mesh_2(CDT& cdt, double parameters::time_limit=0, std::size_t parameters::max_iteration_number=0, diff --git a/Mesh_2/include/CGAL/Mesh_2/Mesh_global_optimizer_2.h b/Mesh_2/include/CGAL/Mesh_2/Mesh_global_optimizer_2.h index 704bccffd72..346b0f85baf 100644 --- a/Mesh_2/include/CGAL/Mesh_2/Mesh_global_optimizer_2.h +++ b/Mesh_2/include/CGAL/Mesh_2/Mesh_global_optimizer_2.h @@ -308,7 +308,7 @@ private: // Find the minimum value do { - if(face->is_in_domain()) + if (!cdt_.is_infinite(face)) min_sqr = (std::min)(min_sqr, sq_circumradius(face)); face++; } diff --git a/Point_set_3/include/CGAL/Point_set_3.h b/Point_set_3/include/CGAL/Point_set_3.h index db6726c8ae4..b9d6cb3d473 100644 --- a/Point_set_3/include/CGAL/Point_set_3.h +++ b/Point_set_3/include/CGAL/Point_set_3.h @@ -413,7 +413,9 @@ public: else { -- m_nb_removed; - return m_indices.end() - m_nb_removed - 1; + iterator out = m_indices.end() - m_nb_removed - 1; + m_base.reset(*out); + return out; } } diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h index e94fa36c967..ecad70671ce 100644 --- a/Point_set_processing_3/include/CGAL/estimate_scale.h +++ b/Point_set_processing_3/include/CGAL/estimate_scale.h @@ -71,7 +71,8 @@ class Quick_multiscale_approximate_knn_distance::reference + operator() (const ValueType& v) const { return get(point_pmap, v); } }; std::size_t m_cluster_size; @@ -227,7 +228,8 @@ class Quick_multiscale_approximate_knn_distance::reference + operator() (const ValueType& v) const { return get(point_pmap, v); } }; template @@ -243,7 +245,8 @@ class Quick_multiscale_approximate_knn_distance::reference + p2 = get(ppmap.point_pmap, i); return value_type (p2.x(), p2.y(), 0.); } @@ -368,7 +371,8 @@ public: FT nb = 0.; std::size_t index = 0; - const typename Kernel::Point_2& pquery = get(point_pmap, *query); + typename boost::property_traits::reference + pquery = get(point_pmap, *query); for (std::size_t t = 0; t < m_point_sets.size(); ++ t) { std::size_t size = ((t == m_point_sets.size() - 1) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/parameters_interface.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/parameters_interface.h index 9247d934a7f..04e9e59b591 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/parameters_interface.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/parameters_interface.h @@ -32,9 +32,9 @@ CGAL_add_named_parameter(protect_constraints_t, protect_constraints, protect_con CGAL_add_named_parameter(relax_constraints_t, relax_constraints, relax_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, random_uniform_sampling) -CGAL_add_named_parameter(grid_sampling_t, grid_sampling, grid_sampling) -CGAL_add_named_parameter(monte_carlo_sampling_t, monte_carlo_sampling, monte_carlo_sampling) +CGAL_add_named_parameter(random_uniform_sampling_t, random_uniform_sampling, use_random_uniform_sampling) +CGAL_add_named_parameter(grid_sampling_t, grid_sampling, use_grid_sampling) +CGAL_add_named_parameter(monte_carlo_sampling_t, monte_carlo_sampling, use_monte_carlo_sampling) CGAL_add_named_parameter(do_sample_edges_t, do_sample_edges, do_sample_edges) CGAL_add_named_parameter(do_sample_vertices_t, do_sample_vertices, do_sample_vertices) CGAL_add_named_parameter(do_sample_faces_t, do_sample_faces, do_sample_faces) diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h index fad2c1986b2..5c3e63307d3 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/measure.h @@ -36,6 +36,8 @@ #include #include +#include // needed for CGAL::exact(FT)/CGAL::exact(Lazy_exact_nt) + #include #ifdef DOXYGEN_RUNNING diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/measures_test.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/measures_test.cpp index 3565c074914..a63da169ce1 100644 --- a/Polygon_mesh_processing/test/Polygon_mesh_processing/measures_test.cpp +++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/measures_test.cpp @@ -1,10 +1,10 @@ +#include #include #include #include #include -#include #include #include diff --git a/Ridges_3/doc/Ridges_3/CGAL/Ridges.h b/Ridges_3/doc/Ridges_3/CGAL/Ridges.h index b08aa4753a7..1c6cec82ad6 100644 --- a/Ridges_3/doc/Ridges_3/CGAL/Ridges.h +++ b/Ridges_3/doc/Ridges_3/CGAL/Ridges.h @@ -170,21 +170,21 @@ Ridge_approximation(const TriangleMesh &tm, Outputs ridges of types `MAX_ELLIPTIC_RIDGE` and `MAX_HYPERBOLIC_RIDGE`. \tparam OutputIterator an output iterator wìth value type `Ridge_line*`. */ -template OutputIterator compute_max_ridges(OutputIterator it, Rdge_order ord = Ridge_order_3); +template OutputIterator compute_max_ridges(OutputIterator it, Ridge_order ord = Ridge_order_3); /*! Outputs ridges of types `MIN_ELLIPTIC_RIDGE` and `MIN_HYPERBOLIC_RIDGE`. \tparam OutputIterator an output iterator with value type `Ridge_line*`. */ -template OutputIterator compute_min_ridges(OutputIterator it, Rdge_order ord = Ridge_order_3); +template OutputIterator compute_min_ridges(OutputIterator it, Ridge_order ord = Ridge_order_3); /*! Outputs ridges of types `MAX_CREST_RIDGE` and `MIN_CREST_RIDGE`. \tparam OutputIterator is an output iterator with value type `Ridge_line*`. */ -template OutputIterator compute_crest_ridges(OutputIterator it, Rdge_order ord = Ridge_order_3); +template OutputIterator compute_crest_ridges(OutputIterator it, Ridge_order ord = Ridge_order_3); /// @} @@ -226,7 +226,7 @@ typedef typename TriangleMesh::Traits::FT FT; A halfedge crossed by a ridge is paired with the barycentric coordinate of the crossing point. */ -typedef std::pair< halfedge_descriptor, FT> Ridge_halfhedge; +typedef std::pair< halfedge_descriptor, FT> Ridge_halfedge; /// @} @@ -261,7 +261,7 @@ FT sharpness() const; /*! */ -const std::list* line() const; +const std::list* line() const; /// @} diff --git a/Ridges_3/include/CGAL/Ridges.h b/Ridges_3/include/CGAL/Ridges.h index 0da44d6852f..1aaaaec915a 100644 --- a/Ridges_3/include/CGAL/Ridges.h +++ b/Ridges_3/include/CGAL/Ridges.h @@ -69,7 +69,8 @@ public: typedef typename Kernel::FT FT; typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef std::pair< halfedge_descriptor, FT> ridge_halfhedge; + typedef std::pair< halfedge_descriptor, FT> Ridge_halfedge; + typedef Ridge_halfedge ridge_halfhedge; //kept for backward compatibility Ridge_type line_type() const {return m_line_type;} Ridge_type& line_type() {return m_line_type;} @@ -80,8 +81,8 @@ public: const FT sharpness() const {return m_sharpness;} FT& sharpness() {return m_sharpness;} - const std::list* line() const { return &m_line;} - std::list* line() { return &m_line;} + const std::list* line() const { return &m_line;} + std::list* line() { return &m_line;} //constructor Ridge_line(const TriangleMesh& P); @@ -100,7 +101,7 @@ protected: //one of MAX_ELLIPTIC_RIDGE, MAX_HYPERBOLIC_RIDGE, MAX_CREST_RIDGE, //MIN_ELLIPTIC_RIDGE, MIN_HYPERBOLIC_RIDGE or MIN_CREST_RIDGE Ridge_type m_line_type; - std::list m_line; + std::list m_line; FT m_strength;// = integral of ppal curvature along the line FT m_sharpness;// = (integral of second derivative of curvature // along the line) multiplied by the squared of @@ -131,7 +132,7 @@ dump_4ogl(std::ostream& out_stream, << strength() << " " << sharpness() << " "; typedef typename boost::property_traits::value_type Point_3; - typename std::list::const_iterator + typename std::list::const_iterator iter = line()->begin(), ite = line()->end(); for (;iter!=ite;iter++){ @@ -156,7 +157,7 @@ dump_verbose(std::ostream& out_stream, VertexPointMap vpm) const << "Sharpness is : " << sharpness() << std::endl << "Polyline point coordinates are : " << std::endl; - typename std::list::const_iterator + typename std::list::const_iterator iter = line()->begin(), ite = line()->end(); for (;iter!=ite;iter++){ @@ -204,7 +205,8 @@ class Ridge_approximation CGAL_static_assertion((boost::is_same::value)); CGAL_static_assertion((boost::is_same::value)); - typedef std::pair< halfedge_descriptor, FT> Ridge_halfhedge; + typedef std::pair< halfedge_descriptor, FT> Ridge_halfedge; + typedef Ridge_halfedge Ridge_halfhedge; // kept for backward compatibility typedef CGAL::Ridge_line Ridge_line; Ridge_approximation(const TriangleMesh &P, @@ -502,7 +504,7 @@ facet_ridge_type(const face_descriptor f, halfedge_descriptor& he1, halfedge_des //check for regular facet //i.e. if there is a coherent orientation of ppal dir at the facet vertices - if ( d1[v1]*d1[v2] * d1[v1]*d1[v3] * d1[v2]*d1[v3] < 0 ) + if ( get(d1,v1)*get(d1,v2) * get(d1,v1)*get(d1,v3) * get(d1,v2)*get(d1,v3) < 0 ) return NO_RIDGE; //determine potential crest color @@ -511,11 +513,11 @@ facet_ridge_type(const face_descriptor f, halfedge_descriptor& he1, halfedge_des Ridge_type crest_color = NO_RIDGE; if (r_type == CREST_RIDGE) { - if ( CGAL::abs(k1[v1]+k1[v2]+k1[v3]) > CGAL::abs(k2[v1]+k2[v2]+k2[v3]) ) + if ( CGAL::abs(get(k1,v1)+get(k1,v2)+get(k1,v3)) > CGAL::abs(get(k2, v1)+get(k2,v2)+get(k2,v3)) ) crest_color = MAX_CREST_RIDGE; - if ( CGAL::abs(k1[v1]+k1[v2]+k1[v3]) < CGAL::abs(k2[v1]+k2[v2]+k2[v3]) ) + if ( CGAL::abs(get(k1,v1)+get(k1,v2)+get(k1,v3)) < CGAL::abs(get(k2,v1)+get(k2,v2)+get(k2,v3)) ) crest_color = MIN_CREST_RIDGE; - if ( CGAL::abs(k1[v1]+k1[v2]+k1[v3]) == CGAL::abs(k2[v1]+k2[v2]+k2[v3]) ) + if ( CGAL::abs(get(k1,v1)+get(k1,v2)+get(k1,v3)) == CGAL::abs(get(k2,v1)+get(k2,v2)+get(k2,v3)) ) return NO_RIDGE; } @@ -587,15 +589,15 @@ xing_on_edge(const halfedge_descriptor he, bool& is_crossed, Ridge_interrogation is_crossed = false; FT sign = 0; FT b_p, b_q; // extremalities at p and q for he: p->q - Vector_3 d_p = d1[target(opposite(he,P),P)], - d_q = d1[target(he,P)]; //ppal dir + Vector_3 d_p = get(d1,target(opposite(he,P),P)), + d_q = get(d1,target(he,P)); //ppal dir if ( color == MAX_RIDGE ) { - b_p = b0[target(opposite(he,P),P)]; - b_q = b0[target(he,P)]; + b_p = get(b0,target(opposite(he,P),P)); + b_q = get(b0,target(he,P)); } else { - b_p = b3[target(opposite(he,P),P)]; - b_q = b3[target(he,P)]; + b_p = get(b3,target(opposite(he,P),P)); + b_q = get(b3,target(he,P)); } if ( b_p == 0 && b_q == 0 ) return; if ( b_p == 0 && b_q !=0 ) sign = d_p*d_q * b_q; @@ -620,13 +622,13 @@ tag_as_elliptic_hyperbolic(const Ridge_interrogation_type color, FT coord1, coord2; if (color == MAX_RIDGE) { - coord1 = CGAL::abs(b0[v_q1]) / ( CGAL::abs(b0[v_p1]) + CGAL::abs(b0[v_q1]) ); - coord2 = CGAL::abs(b0[v_q2]) / ( CGAL::abs(b0[v_p2]) + CGAL::abs(b0[v_q2]) ); + coord1 = CGAL::abs(get(b0,v_q1)) / ( CGAL::abs(get(b0,v_p1)) + CGAL::abs(get(b0,v_q1)) ); + coord2 = CGAL::abs(get(b0,v_q2)) / ( CGAL::abs(get(b0,v_p2)) + CGAL::abs(get(b0,v_q2)) ); } else { - coord1 = CGAL::abs(b3[v_q1]) / ( CGAL::abs(b3[v_p1]) + CGAL::abs(b3[v_q1]) ); - coord2 = CGAL::abs(b3[v_q2]) / ( CGAL::abs(b3[v_p2]) + CGAL::abs(b3[v_q2]) ); + coord1 = CGAL::abs(get(b3,v_q1)) / ( CGAL::abs(get(b3,v_p1)) + CGAL::abs(get(b3,v_q1)) ); + coord2 = CGAL::abs(get(b3,v_q2)) / ( CGAL::abs(get(b3,v_p2)) + CGAL::abs(get(b3,v_q2)) ); } if ( tag_order == Ridge_order_3 ) { @@ -649,10 +651,10 @@ tag_as_elliptic_hyperbolic(const Ridge_interrogation_type color, // of Pi at the two crossing points FT sign_P; if (color == MAX_RIDGE) - sign_P = P1[v_p1]*coord1 + P1[v_q1]*(1-coord1) - + P1[v_p2]*coord2 + P1[v_q2]*(1-coord2); - else sign_P = P2[v_p1]*coord1 + P2[v_q1]*(1-coord1) - + P2[v_p2]*coord2 + P2[v_q2]*(1-coord2); + sign_P = get(P1,v_p1)*coord1 + get(P1,v_q1)*(1-coord1) + + get(P1,v_p2)*coord2 + get(P1,v_q2)*(1-coord2); + else sign_P = get(P2,v_p1)*coord1 + get(P2,v_q1)*(1-coord1) + + get(P2,v_p2)*coord2 + get(P2,v_q2)*(1-coord2); if ( sign_P < 0 ) return true; else return false; } @@ -672,20 +674,20 @@ b_sign_pointing_to_ridge(const vertex_descriptor v1, Vector_3 r = r2 - r1, dv1, dv2, dv3; FT bv1, bv2, bv3; if ( color == MAX_RIDGE ) { - bv1 = b0[v1]; - bv2 = b0[v2]; - bv3 = b0[v3]; - dv1 = d1[v1]; - dv2 = d1[v2]; - dv3 = d1[v3]; + bv1 = get(b0,v1); + bv2 = get(b0,v2); + bv3 = get(b0,v3); + dv1 = get(d1,v1); + dv2 = get(d1,v2); + dv3 = get(d1,v3); } else { - bv1 = b3[v1]; - bv2 = b3[v2]; - bv3 = b3[v3]; - dv1 = d2[v1]; - dv2 = d2[v2]; - dv3 = d2[v3]; + bv1 = get(b3,v1); + bv2 = get(b3,v2); + bv3 = get(b3,v3); + dv1 = get(d2,v1); + dv2 = get(d2,v2); + dv3 = get(d2,v3); } if ( r != CGAL::NULL_VECTOR ) r = r/CGAL::sqrt(r*r); FT sign1, sign2, sign3; @@ -712,7 +714,7 @@ init_ridge_line(Ridge_line* ridge_line, const Ridge_type r_type) { ridge_line->line_type() = r_type; - ridge_line->line()->push_back(Ridge_halfhedge(h1, bary_coord(h1,r_type))); + ridge_line->line()->push_back(Ridge_halfedge(h1, bary_coord(h1,r_type))); addback(ridge_line, h2, r_type); } @@ -735,8 +737,8 @@ addback(Ridge_line* ridge_line, const halfedge_descriptor he, FT k1x, k2x; //abs value of the ppal curvatures at the Xing point on he. FT k_second = 0; // abs value of the second derivative of the curvature // along the line of curvature - k1x = CGAL::abs(k1[v_p]) * coord + CGAL::abs(k1[v_q]) * (1-coord) ; - k2x = CGAL::abs(k2[v_p]) * coord + CGAL::abs(k2[v_q]) * (1-coord) ; + k1x = CGAL::abs(get(k1,v_p)) * coord + CGAL::abs(get(k1,v_q)) * (1-coord) ; + k2x = CGAL::abs(get(k2,v_p)) * coord + CGAL::abs(get(k2,v_q)) * (1-coord) ; if ( (ridge_line->line_type() == MAX_ELLIPTIC_RIDGE) || (ridge_line->line_type() == MAX_HYPERBOLIC_RIDGE) @@ -744,7 +746,7 @@ addback(Ridge_line* ridge_line, const halfedge_descriptor he, ridge_line->strength() += k1x * CGAL::sqrt(segment * segment); if (tag_order == Ridge_order_4) { if (k1x != k2x) - k_second =CGAL::abs(( CGAL::abs(P1[v_p]) * coord + CGAL::abs(P1[v_q]) * (1-coord) )/(k1x-k2x)); + k_second =CGAL::abs(( CGAL::abs(get(P1,v_p)) * coord + CGAL::abs(get(P1,v_q)) * (1-coord) )/(k1x-k2x)); ridge_line->sharpness() += k_second * CGAL::sqrt(segment * segment) * squared_model_size; } } if ( (ridge_line->line_type() == MIN_ELLIPTIC_RIDGE) @@ -753,10 +755,10 @@ addback(Ridge_line* ridge_line, const halfedge_descriptor he, ridge_line->strength() += k2x * CGAL::sqrt(segment * segment); if (tag_order == Ridge_order_4) { if (k1x != k2x) - k_second =CGAL::abs(( CGAL::abs(P2[v_p]) * coord + CGAL::abs(P2[v_q]) * (1-coord) )/(k1x-k2x)); + k_second =CGAL::abs(( CGAL::abs(get(P2,v_p)) * coord + CGAL::abs(get(P2,v_q)) * (1-coord) )/(k1x-k2x)); ridge_line->sharpness() += k_second * CGAL::sqrt(segment * segment) * squared_model_size; } } - ridge_line->line()->push_back( Ridge_halfhedge(he, coord)); + ridge_line->line()->push_back( Ridge_halfedge(he, coord)); } @@ -779,8 +781,8 @@ addfront(Ridge_line* ridge_line, FT k1x, k2x; //abs value of the ppal curvatures at the Xing point on he. FT k_second = 0.; // abs value of the second derivative of the curvature // along the line of curvature - k1x = CGAL::abs(k1[v_p]) * coord + CGAL::abs(k1[v_q]) * (1-coord) ; - k2x = CGAL::abs(k2[v_p]) * coord + CGAL::abs(k2[v_q]) * (1-coord) ; + k1x = CGAL::abs(get(k1,v_p)) * coord + CGAL::abs(get(k1,v_q)) * (1-coord) ; + k2x = CGAL::abs(get(k2,v_p)) * coord + CGAL::abs(get(k2,v_q)) * (1-coord) ; if ( (ridge_line->line_type() == MAX_ELLIPTIC_RIDGE) || (ridge_line->line_type() == MAX_HYPERBOLIC_RIDGE) @@ -788,7 +790,7 @@ addfront(Ridge_line* ridge_line, ridge_line->strength() += k1x * CGAL::sqrt(segment * segment); if (tag_order == Ridge_order_4) { if (k1x != k2x) - k_second =CGAL::abs(( CGAL::abs(P1[v_p]) * coord + CGAL::abs(P1[v_q]) * (1-coord) )/(k1x-k2x)); + k_second =CGAL::abs(( CGAL::abs(get(P1,v_p)) * coord + CGAL::abs(get(P1,v_q)) * (1-coord) )/(k1x-k2x)); ridge_line->sharpness() += k_second * CGAL::sqrt(segment * segment) * squared_model_size; } } if ( (ridge_line->line_type() == MIN_ELLIPTIC_RIDGE) @@ -797,10 +799,10 @@ addfront(Ridge_line* ridge_line, ridge_line->strength() += k2x * CGAL::sqrt(segment * segment); if (tag_order == Ridge_order_4) { if (k1x != k2x) - k_second =CGAL::abs(( CGAL::abs(P2[v_p]) * coord + CGAL::abs(P2[v_q]) * (1-coord) )/(k1x-k2x)); + k_second =CGAL::abs(( CGAL::abs(get(P2,v_p)) * coord + CGAL::abs(get(P2,v_q)) * (1-coord) )/(k1x-k2x)); ridge_line->sharpness() += k_second * CGAL::sqrt(segment * segment) * squared_model_size; } } - ridge_line->line()->push_front( Ridge_halfhedge(he, coord)); + ridge_line->line()->push_front( Ridge_halfedge(he, coord)); } @@ -815,14 +817,14 @@ bary_coord(const halfedge_descriptor he, const Ridge_type r_type) if ( (r_type == MAX_ELLIPTIC_RIDGE) || (r_type == MAX_HYPERBOLIC_RIDGE) || (r_type == MAX_CREST_RIDGE) ) { - b_p = b0[target(opposite(he,P),P)]; - b_q = b0[target(he,P)]; + b_p = get(b0,target(opposite(he,P),P)); + b_q = get(b0,target(he,P)); } if ( (r_type == MIN_ELLIPTIC_RIDGE) || (r_type == MIN_HYPERBOLIC_RIDGE) || (r_type == MIN_CREST_RIDGE) ) { - b_p = b3[target(opposite(he,P),P)]; - b_q = b3[target(he,P)]; + b_p = get(b3,target(opposite(he,P),P)); + b_q = get(b3,target(he,P)); } //denominator cannot be 0 since there is no crossing when both extremalities are 0 return CGAL::abs(b_q) / ( CGAL::abs(b_q) + CGAL::abs(b_p) ); diff --git a/Ridges_3/include/CGAL/Umbilics.h b/Ridges_3/include/CGAL/Umbilics.h index 0b7af0874c8..96431b1c170 100644 --- a/Ridges_3/include/CGAL/Umbilics.h +++ b/Ridges_3/include/CGAL/Umbilics.h @@ -194,7 +194,7 @@ compute(OutputIterator umbilics_it, FT size) boost::tie(itb,ite) = vertices(P); for (;itb != ite; itb++) { vertex_descriptor vh = *itb; - umbilicEstimatorVertex = cgal_abs(k1[vh]-k2[vh]); + umbilicEstimatorVertex = cgal_abs(get(k1,vh)-get(k2,vh)); //reset vector, list and bool vces.clear(); contour.clear(); @@ -218,7 +218,7 @@ compute(OutputIterator umbilics_it, FT size) itev = vces.end(); itbv++; for (; itbv != itev; itbv++) - { umbilicEstimatorNeigh = cgal_abs( k1[*itbv] - k2[*itbv] ); + { umbilicEstimatorNeigh = cgal_abs( get(k1,*itbv) - get(k2,*itbv) ); if ( umbilicEstimatorNeigh < umbilicEstimatorVertex ) {is_umbilic = false; break;} } @@ -245,14 +245,14 @@ compute_type(Umbilic& umb) itlast = --umb.contour_list().end(); v = target(*itb, P); - dir = d1[v]; - normal = CGAL::cross_product(d1[v], d2[v]); + dir = get(d1,v); + normal = CGAL::cross_product(get(d1,v), get(d2,v)); //sum angles along the contour do{ itb++; v = target(*itb, P); - dirnext = d1[v]; + dirnext = get(d1,v); cosinus = To_double(dir*dirnext); if (cosinus < 0) {dirnext = dirnext*(-1); cosinus *= -1;} if (cosinus>1) cosinus = 1; @@ -261,13 +261,13 @@ compute_type(Umbilic& umb) else angle = -acos(cosinus); angleSum += angle; dir = dirnext; - normal = CGAL::cross_product(d1[v], d2[v]); + normal = CGAL::cross_product(get(d1,v), get(d2,v)); } while (itb != (itlast)); //angle (v_last, v_0) v = target(*umb.contour_list().begin(), P); - dirnext = d1[v]; + dirnext = get(d1,v); cosinus = To_double(dir*dirnext); if (cosinus < 0) {dirnext = dirnext*(-1); cosinus *= -1;} if (cosinus>1) cosinus = 1; diff --git a/Subdivision_method_3/doc/Subdivision_method_3/Subdivision_method_3.txt b/Subdivision_method_3/doc/Subdivision_method_3/Subdivision_method_3.txt index 39d503841f4..48a87621908 100644 --- a/Subdivision_method_3/doc/Subdivision_method_3/Subdivision_method_3.txt +++ b/Subdivision_method_3/doc/Subdivision_method_3/Subdivision_method_3.txt @@ -127,12 +127,12 @@ There is only one line deserving a detailed explanation: \code{.cpp} -Subdivision_method_3::CatmullClark_subdivision(pmesh, d); +Subdivision_method_3::CatmullClark_subdivision(pmesh, params::number_of_iterations(d)); \endcode `Subdivision_method_3` specifies the namespace of the -subdivision functions. `CatmullClark_subdivision(P, d)` computes the +subdivision functions. `CatmullClark_subdivision(P, params::number_of_iterations(d))` computes the Catmull-Clark subdivision surface of the polygon mesh `pmesh` after `d` iterations of the refinements. The polygon mesh `pmesh` is passed by reference, and is modified (i.e.\ subdivided) by the @@ -150,12 +150,6 @@ in the internal container are sequentially ordered (e.g. This implies that the iterators traverse the primitives in the order of their creations/insertions. -\attention The four subdivision techniques have been made to work with -`Surface_mesh` despite `Surface_mesh` not necessarily inserting new elements at -the end of its containers (see Section \ref sectionSurfaceMesh_properties of the -`Surface_mesh`). However, it is required that the `Surface_mesh` passed as input -is garbage-free, which can be achieved by first calling `Surface_mesh::collect_garbage()`. - Section \ref secRefHost gives detailed explanations on this two restrictions. @@ -187,8 +181,8 @@ three components of the Catmull-Clark subdivision method. \code{.cpp} -template -void PQQ(PolygonMesh& p, Mask mask, int depth) +template +void PQQ(PolygonMesh& p, Mask mask, NamedParameters np) \endcode @@ -203,8 +197,7 @@ new points by cooperating with the `mask`. To implement Catmull-Clark subdivision, `Mask`, the geometry policy, has to realize the geometry masks of Catmull-Clark subdivision. -The parameter `depth` specifies the iterations of the refinement -on the control mesh. +The number of iterations as well as the vertex point map can be specified using the named parameter `np`. To implement the geometry masks, we need to know how a refinement host communicates with its geometry masks. @@ -258,8 +251,8 @@ class CatmullClark_mask_3 { }; void edge_node(halfedge_descriptor edge, Point_3& pt) { - Point_3& p1 = get(vpm,target(edge, pmesh)); - Point_3& p2 = get(vpm,source(edge, pmesh)); + Point_3 p1 = get(vpm,target(edge, pmesh)); + Point_3 p2 = get(vpm,source(edge, pmesh)); Point_3 f1, f2; face_node(face(edge,pmesh), f1); face_node(face(opposite(edge,pmesh),pmesh), f2); @@ -273,11 +266,11 @@ class CatmullClark_mask_3 { typename boost::graph_traits::degree_size_type n = degree(vertex,pmesh); FT Q[] = {0.0, 0.0, 0.0}, R[] = {0.0, 0.0, 0.0}; - Point_3& S = get(vpm,vertex); + Point_3 S = get(vpm,vertex); Point_3 q; for (int i = 0; i < n; i++, ++vcir) { - Point_3& p2 = get(vpm, target(opposite(*vcir,pmesh),pmesh)); + Point_3 p2 = get(vpm, target(opposite(*vcir,pmesh),pmesh)); R[0] += (S[0]+p2[0])/2; R[1] += (S[1]+p2[1])/2; R[2] += (S[2]+p2[2])/2; @@ -305,7 +298,7 @@ call `PQQ()` with the Catmull-Clark masks that we have just defined. \code{.cpp} -PQQ(pmesh, CatmullClark_mask_3(pmesh), depth); +PQQ(pmesh, CatmullClark_mask_3(pmesh), params::number_of_iterations(depth)); \endcode @@ -334,17 +327,17 @@ and \f$ \sqrt{3}\f$ subdivision. \code{.cpp} namespace Subdivision_method_3 { -template -void PQQ(PolygonMesh& pmesh, Mask mask, int step); +template +void PQQ(PolygonMesh& pmesh, Mask mask, NamedParameters np); -template -void PTQ(PolygonMesh& pmesh, Mask mask, int step); +template +void PTQ(PolygonMesh& pmesh, Mask mask, NamedParameters np); -template -void DQQ(PolygonMesh& pmesh, Mask mask, int step) +template +void DQQ(PolygonMesh& pmesh, Mask mask, NamedParameters np) -template -void Sqrt3(PolygonMesh& pmesh, Mask mask, int step) +template +void Sqrt3(PolygonMesh& pmesh, Mask mask, NamedParameters np) } \endcode @@ -461,15 +454,15 @@ is listed below, where `ept` returns the new point splitting \code{.cpp} void border_node(halfedge_descriptor edge, Point_3& ept, Point_3& vpt) { - Point_3& ep1 = get(vpm, target(edge, pmesh)); - Point_3& ep2 = get(vpm, target(opposite(edge, pmesh), pmesh)); + Point_3 ep1 = get(vpm, target(edge, pmesh)); + Point_3 ep2 = get(vpm, target(opposite(edge, pmesh), pmesh)); ept = Point_3((ep1[0]+ep2[0])/2, (ep1[1]+ep2[1])/2, (ep1[2]+ep2[2])/2); Halfedge_around_target_circulator vcir(edge, pmesh); - Point_3& vp1 = get(vpm,target(opposite(*vcir, pmesh ), pmesh)); - Point_3& vp0 = get(vpm, target(*vcir, pmesh)); + Point_3 vp1 = get(vpm,target(opposite(*vcir, pmesh ), pmesh)); + Point_3 vp0 = get(vpm, target(*vcir, pmesh)); --vcir; - Point_3& vp_1 = get(vpm, target(opposite(*vcir, pmesh), pmesh)); + Point_3 vp_1 = get(vpm, target(opposite(*vcir, pmesh), pmesh)); vpt = Point_3((vp_1[0] + 6*vp0[0] + vp1[0])/8, (vp_1[1] + 6*vp0[1] + vp1[1])/8, (vp_1[2] + 6*vp0[2] + vp1[2])/8 ); @@ -533,24 +526,24 @@ based on that kernel. \code{.cpp} namespace Subdivision_method_3 { - template - void CatmullClark_subdivision(PolygonMesh& pmesh, int step = 1) { - PQQ(pmesh, CatmullClark_mask_3(pmesh), step); + template + void CatmullClark_subdivision(PolygonMesh& pmesh, NamedParameters np) { + PQQ(pmesh, CatmullClark_mask_3(pmesh), np); } - template - void Loop_subdivision(PolygonMesh& pmesh, int step = 1) { - PTQ(pmesh, Loop_mask_3(pmesh) , step); + template + void Loop_subdivision(PolygonMesh& pmesh, NamedParameters np) { + PTQ(pmesh, Loop_mask_3(pmesh) , np); } - template - void DooSabin_subdivision(PolygonMesh& pmesh, int step = 1) { - DQQ(pmesh, DooSabin_mask_3(pmesh), step); + template + void DooSabin_subdivision(PolygonMesh& pmesh, NamedParameters np) { + DQQ(pmesh, DooSabin_mask_3(pmesh), np); } - template - void Sqrt3_subdivision(PolygonMesh& pmesh, int step = 1) { - Sqrt3(pmesh, Sqrt3_mask_3(pmesh), step); + template + void Sqrt3_subdivision(PolygonMesh& pmesh, NamedParameters np) { + Sqrt3(pmesh, Sqrt3_mask_3(pmesh), np); } } diff --git a/Subdivision_method_3/examples/Subdivision_method_3/CatmullClark_subdivision.cpp b/Subdivision_method_3/examples/Subdivision_method_3/CatmullClark_subdivision.cpp index b0cc8a9138c..cfaf8938158 100644 --- a/Subdivision_method_3/examples/Subdivision_method_3/CatmullClark_subdivision.cpp +++ b/Subdivision_method_3/examples/Subdivision_method_3/CatmullClark_subdivision.cpp @@ -16,6 +16,7 @@ typedef CGAL::Surface_mesh PolygonMesh; using namespace std; using namespace CGAL; +namespace params = CGAL::Polygon_mesh_processing::parameters; int main(int argc, char** argv) { if (argc > 4) { @@ -40,7 +41,7 @@ int main(int argc, char** argv) { Timer t; t.start(); - Subdivision_method_3::CatmullClark_subdivision(pmesh, d); + Subdivision_method_3::CatmullClark_subdivision(pmesh, params::number_of_iterations(d)); std::cerr << "Done (" << t.time() << " s)" << std::endl; std::ofstream out(out_file); diff --git a/Subdivision_method_3/examples/Subdivision_method_3/Customized_subdivision.cpp b/Subdivision_method_3/examples/Subdivision_method_3/Customized_subdivision.cpp index 5dbdaf8e198..dc93f81bc89 100644 --- a/Subdivision_method_3/examples/Subdivision_method_3/Customized_subdivision.cpp +++ b/Subdivision_method_3/examples/Subdivision_method_3/Customized_subdivision.cpp @@ -16,6 +16,7 @@ typedef CGAL::Simple_cartesian Kernel; typedef CGAL::Surface_mesh PolygonMesh; using namespace std; using namespace CGAL; +namespace params = CGAL::Polygon_mesh_processing::parameters; // ====================================================================== template @@ -27,6 +28,7 @@ class WLoop_mask_3 { typedef typename boost::property_map::type Vertex_pmap; typedef typename boost::property_traits::value_type Point; + typedef typename boost::property_traits::reference Point_ref; PolygonMesh& pmesh; Vertex_pmap vpm; @@ -37,10 +39,10 @@ public: {} void edge_node(halfedge_descriptor hd, Point& pt) { - Point& p1 = get(vpm, target(hd,pmesh)); - Point& p2 = get(vpm, target(opposite(hd,pmesh),pmesh)); - Point& f1 = get(vpm, target(next(hd,pmesh),pmesh)); - Point& f2 = get(vpm, target(next(opposite(hd,pmesh),pmesh),pmesh)); + Point_ref p1 = get(vpm, target(hd,pmesh)); + Point_ref p2 = get(vpm, target(opposite(hd,pmesh),pmesh)); + Point_ref f1 = get(vpm, target(next(hd,pmesh),pmesh)); + Point_ref f2 = get(vpm, target(next(opposite(hd,pmesh),pmesh),pmesh)); pt = Point((3*(p1[0]+p2[0])+f1[0]+f2[0])/8, (3*(p1[1]+p2[1])+f1[1]+f2[1])/8, @@ -48,12 +50,12 @@ public: } void vertex_node(vertex_descriptor vd, Point& pt) { double R[] = {0.0, 0.0, 0.0}; - Point& S = get(vpm,vd); + Point_ref S = get(vpm,vd); std::size_t n = 0; BOOST_FOREACH(halfedge_descriptor hd, halfedges_around_target(vd, pmesh)){ ++n; - Point& p = get(vpm, target(opposite(hd,pmesh),pmesh)); + Point_ref p = get(vpm, target(opposite(hd,pmesh),pmesh)); R[0] += p[0]; R[1] += p[1]; R[2] += p[2]; } @@ -71,15 +73,15 @@ public: } void border_node(halfedge_descriptor hd, Point& ept, Point& vpt) { - Point& ep1 = get(vpm, target(hd,pmesh)); - Point& ep2 = get(vpm, target(opposite(hd,pmesh),pmesh)); + Point_ref ep1 = get(vpm, target(hd,pmesh)); + Point_ref ep2 = get(vpm, target(opposite(hd,pmesh),pmesh)); ept = Point((ep1[0]+ep2[0])/2, (ep1[1]+ep2[1])/2, (ep1[2]+ep2[2])/2); Halfedge_around_target_circulator vcir(hd,pmesh); - Point& vp1 = get(vpm, target(opposite(*vcir,pmesh),pmesh)); - Point& vp0 = get(vpm, target(*vcir,pmesh)); + Point_ref vp1 = get(vpm, target(opposite(*vcir,pmesh),pmesh)); + Point_ref vp0 = get(vpm, target(*vcir,pmesh)); --vcir; - Point& vp_1 = get(vpm,target(opposite(*vcir,pmesh),pmesh)); + Point_ref vp_1 = get(vpm,target(opposite(*vcir,pmesh),pmesh)); vpt = Point((vp_1[0] + 6*vp0[0] + vp1[0])/8, (vp_1[1] + 6*vp0[1] + vp1[1])/8, (vp_1[2] + 6*vp0[2] + vp1[2])/8 ); @@ -109,7 +111,7 @@ int main(int argc, char **argv) { Timer t; t.start(); - Subdivision_method_3::PTQ(pmesh, WLoop_mask_3(pmesh), d); + Subdivision_method_3::PTQ(pmesh, WLoop_mask_3(pmesh), params::number_of_iterations(d)); std::cerr << "Done (" << t.time() << " s)" << std::endl; std::ofstream out(out_file); diff --git a/Subdivision_method_3/examples/Subdivision_method_3/DooSabin_subdivision.cpp b/Subdivision_method_3/examples/Subdivision_method_3/DooSabin_subdivision.cpp index 72cbe53e6f1..016bf1f8d33 100644 --- a/Subdivision_method_3/examples/Subdivision_method_3/DooSabin_subdivision.cpp +++ b/Subdivision_method_3/examples/Subdivision_method_3/DooSabin_subdivision.cpp @@ -18,6 +18,7 @@ typedef CGAL::Polyhedron_3 PolygonMesh; using namespace std; using namespace CGAL; +namespace params = CGAL::Polygon_mesh_processing::parameters; int main(int argc, char **argv) { if (argc > 4) { @@ -42,7 +43,7 @@ int main(int argc, char **argv) { Timer t; t.start(); - Subdivision_method_3::DooSabin_subdivision(pmesh, d); + Subdivision_method_3::DooSabin_subdivision(pmesh, params::number_of_iterations(d)); std::cerr << "Done (" << t.time() << " s)" << std::endl; std::ofstream out(out_file); diff --git a/Subdivision_method_3/examples/Subdivision_method_3/Loop_subdivision.cpp b/Subdivision_method_3/examples/Subdivision_method_3/Loop_subdivision.cpp index 88526ed999a..8298450fdda 100644 --- a/Subdivision_method_3/examples/Subdivision_method_3/Loop_subdivision.cpp +++ b/Subdivision_method_3/examples/Subdivision_method_3/Loop_subdivision.cpp @@ -16,6 +16,7 @@ typedef CGAL::Surface_mesh PolygonMesh; using namespace std; using namespace CGAL; +namespace params = CGAL::Polygon_mesh_processing::parameters; int main(int argc, char **argv) { if (argc > 4) { @@ -40,7 +41,7 @@ int main(int argc, char **argv) { Timer t; t.start(); - Subdivision_method_3::Loop_subdivision(pmesh,d); + Subdivision_method_3::Loop_subdivision(pmesh, params::number_of_iterations(d)); std::cerr << "Done (" << t.time() << " s)" << std::endl; std::ofstream out(out_file); diff --git a/Subdivision_method_3/examples/Subdivision_method_3/Sqrt3_subdivision.cpp b/Subdivision_method_3/examples/Subdivision_method_3/Sqrt3_subdivision.cpp index cec6f63b4e5..447ae6e1671 100644 --- a/Subdivision_method_3/examples/Subdivision_method_3/Sqrt3_subdivision.cpp +++ b/Subdivision_method_3/examples/Subdivision_method_3/Sqrt3_subdivision.cpp @@ -16,6 +16,7 @@ typedef CGAL::Surface_mesh PolygonMesh; using namespace std; using namespace CGAL; +namespace params = CGAL::Polygon_mesh_processing::parameters; int main(int argc, char **argv) { if (argc > 4) { @@ -40,7 +41,7 @@ int main(int argc, char **argv) { Timer t; t.start(); - Subdivision_method_3::Sqrt3_subdivision(pmesh,d); + Subdivision_method_3::Sqrt3_subdivision(pmesh, params::number_of_iterations(d)); std::cerr << "Done (" << t.time() << " s)" << std::endl; std::ofstream out(out_file); diff --git a/Subdivision_method_3/include/CGAL/Subdivision_method_3/internal/subdivision_hosts_impl_3.h b/Subdivision_method_3/include/CGAL/Subdivision_method_3/internal/subdivision_hosts_impl_3.h index df791dd8a48..cf467a3112d 100644 --- a/Subdivision_method_3/include/CGAL/Subdivision_method_3/internal/subdivision_hosts_impl_3.h +++ b/Subdivision_method_3/include/CGAL/Subdivision_method_3/internal/subdivision_hosts_impl_3.h @@ -33,11 +33,8 @@ #include #include -#include #include -#include #include -#include #include #include @@ -49,18 +46,34 @@ namespace Subdivision_method_3 { namespace internal { +template +void call_reserve(PolygonMesh& pm, std::size_t v, std::size_t e, std::size_t f) +{ + typedef boost::graph_traits GT; + reserve(pm, + static_cast(v), + static_cast(e), + static_cast(f)); +} + template void PQQ_1step(Poly& p, VertexPointMap vpm, Mask mask) { typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; typedef typename boost::graph_traits::edge_descriptor edge_descriptor; - - typedef typename boost::graph_traits::vertex_iterator vertex_iterator; - typedef typename boost::graph_traits::edge_iterator edge_iterator; - typedef typename boost::graph_traits::face_iterator face_iterator; + typedef typename boost::graph_traits::face_descriptor face_descriptor; typedef Halfedge_around_face_circulator Halfedge_around_facet_circulator; + //First back up initial vertices/faces/edges + std::vector p_vertices(vertices(p).first, vertices(p).second); + std::vector p_faces(faces(p).first, faces(p).second); + std::vector p_edges(edges(p).first, edges(p).second); + + std::size_t num_v = p_vertices.size(); + std::size_t num_e = p_edges.size(); + std::size_t num_f = p_faces.size(); + // Build a new vertices buffer has the following structure // // 0 1 ... e_begin ... f_begin ... (end_of_buffer) @@ -69,34 +82,29 @@ void PQQ_1step(Poly& p, VertexPointMap vpm, Mask mask) { // f_begin ... (end) : store the positions of the face-vertices // The index of the vertices buffer should 1-1 map to the distance // of the corresponding iterator to the begin of the iterator. - typename boost::graph_traits::vertices_size_type num_vertex = num_vertices(p); - typename boost::graph_traits::halfedges_size_type num_edge = num_halfedges(p)/2; - typename boost::graph_traits::faces_size_type num_facet = num_faces(p); // We need to reserve the memory to prevent reallocation. - reserve(p,num_vertex+num_edge+num_facet, 4*2*num_edge, 4*num_edge/2); + call_reserve(p,num_v+num_e+num_f, 4*2*num_e, 4*num_e/2); typedef typename boost::property_traits::value_type Point; - Point* vertex_point_buffer = new Point[num_vertex + num_edge + num_facet]; - Point* edge_point_buffer = vertex_point_buffer + num_vertex; - Point* face_point_buffer = edge_point_buffer + num_edge; + Point* vertex_point_buffer = new Point[num_v + num_e + num_f]; + Point* edge_point_buffer = vertex_point_buffer + num_v; + Point* face_point_buffer = edge_point_buffer + num_e; int i=0; boost::unordered_map v_index; - BOOST_FOREACH(vertex_descriptor vh, vertices(p)){ + BOOST_FOREACH(vertex_descriptor vh, p_vertices){ v_index[vh]= i++; } - std::vector v_onborder(num_vertex); - face_iterator fitr = faces(p).first; - for (size_t i = 0; i < num_facet; i++, ++fitr) + std::vector v_onborder(num_v); + typename std::vector::iterator fitr=p_faces.begin(); + for (std::size_t i = 0; i < num_f; i++, ++fitr) mask.face_node(*fitr, face_point_buffer[i]); - - { std::size_t i = 0; - BOOST_FOREACH(edge_descriptor ed, edges(p)){ + BOOST_FOREACH(edge_descriptor ed, p_edges){ if(is_border(ed,p)){ halfedge_descriptor h=halfedge(ed,p); if (is_border(h,p)) h=opposite(h,p); @@ -111,8 +119,8 @@ void PQQ_1step(Poly& p, VertexPointMap vpm, Mask mask) { } } - vertex_iterator vitr = vertices(p).first; - for (size_t i = 0; i < num_vertex; i++, ++vitr) + typename std::vector::iterator vitr = p_vertices.begin(); + for (std::size_t i = 0; i < num_v; i++, ++vitr) if (!v_onborder[v_index[*vitr]]) mask.vertex_node(*vitr, vertex_point_buffer[i]); // Build the connectivity using insert_vertex() and insert_edge() @@ -123,16 +131,16 @@ void PQQ_1step(Poly& p, VertexPointMap vpm, Mask mask) { // 4. insert_edge() between all other new inserted vertices of step 1 and // the new inserted vertex of step 3 // Step 1. - edge_iterator eitr = edges(p).first; - for (size_t i = 0; i < num_edge; i++, ++eitr) { + typename std::vector::iterator eitr = p_edges.begin(); + for (std::size_t i = 0; i < num_e; i++, ++eitr) { vertex_descriptor vh = insert_vertex(p, halfedge(*eitr,p)); put(vpm, vh, edge_point_buffer[i]); } - fitr = faces(p).first; // TODO: the topoloy modification can be done by a template function // and that gives the user a chance to create new topological masks. - for (size_t i = 0; i < num_facet; i++, ++fitr) { + fitr = p_faces.begin(); + for (std::size_t i = 0; i < num_f; i++, ++fitr) { // Step 2. Halfedge_around_facet_circulator hcir_begin(halfedge(*fitr,p),p); Halfedge_around_facet_circulator hcir = hcir_begin; @@ -161,8 +169,8 @@ void PQQ_1step(Poly& p, VertexPointMap vpm, Mask mask) { // Update the geometry data of the newly inserted vertices by the // vertices buffer - vitr = vertices(p).first; - for (size_t i = 0; i < num_vertex; i++, ++vitr) + vitr = p_vertices.begin(); + for (std::size_t i = 0; i < num_v; i++, ++vitr) put(vpm, *vitr, vertex_point_buffer[i]); // CGAL_postcondition(p.is_valid()); @@ -175,15 +183,21 @@ void PTQ_1step(Poly& p, VertexPointMap vpm, Mask mask) { typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; typedef typename boost::graph_traits::edge_descriptor edge_descriptor; - - typedef typename boost::graph_traits::vertex_iterator vertex_iterator; - typedef typename boost::graph_traits::edge_iterator edge_iterator; - typedef typename boost::graph_traits::face_iterator face_iterator; + typedef typename boost::graph_traits::face_descriptor face_descriptor; typedef Halfedge_around_face_circulator Halfedge_around_face_circulator; typedef typename boost::property_traits::value_type Point; + //First back up initial vertices/faces/edges + std::vector p_vertices(vertices(p).first, vertices(p).second); + std::vector p_faces(faces(p).first, faces(p).second); + std::vector p_edges(edges(p).first, edges(p).second); + + std::size_t num_v = p_vertices.size(); + std::size_t num_e = p_edges.size(); + std::size_t num_f = p_faces.size(); + // Build a new vertices buffer has the following structure // // 0 1 ... e_begin ... f_begin ... (end_of_buffer) @@ -191,26 +205,23 @@ void PTQ_1step(Poly& p, VertexPointMap vpm, Mask mask) { // e_begin ... (end) : store the positions of the edge-vertices // The index of the vertices buffer should 1-1 map to the distance // of the corresponding iterator to the begin of the iterator. - typename boost::graph_traits::vertices_size_type num_vertex = num_vertices(p); - typename boost::graph_traits::halfedges_size_type num_edge = num_halfedges(p)/2; - typename boost::graph_traits::faces_size_type num_facet = num_faces(p); // We need to reserve the memory to prevent reallocation. - reserve(p,num_vertex+num_edge, 2*2*num_edge, 4*num_edge/2); + call_reserve(p,num_v + num_e, 2*2*num_e, 4*num_e/2); - Point* vertex_point_buffer = new Point[num_vertex + num_edge]; - Point* edge_point_buffer = vertex_point_buffer + num_vertex; + Point* vertex_point_buffer = new Point[num_v + num_e]; + Point* edge_point_buffer = vertex_point_buffer + num_v; int i=0; boost::unordered_map v_index; - BOOST_FOREACH(vertex_descriptor vh, vertices(p)){ + BOOST_FOREACH(vertex_descriptor vh, p_vertices){ v_index[vh]= i++; } - std::vector v_onborder(num_vertex); + std::vector v_onborder(num_v); { std::size_t i = 0; - BOOST_FOREACH(edge_descriptor ed, edges(p)){ + BOOST_FOREACH(edge_descriptor ed, p_edges){ if(! is_border(ed,p)){ mask.edge_node(halfedge(ed,p), edge_point_buffer[i]); } else{ @@ -223,8 +234,8 @@ void PTQ_1step(Poly& p, VertexPointMap vpm, Mask mask) { ++i; } } - vertex_iterator vitr = vertices(p).first; - for (size_t i = 0; i < num_vertex; i++, ++vitr) + typename std::vector::iterator vitr = p_vertices.begin(); + for (std::size_t i = 0; i < num_v; i++, ++vitr) if (!v_onborder[i]) mask.vertex_node(*vitr, vertex_point_buffer[i]); // Build the connectivity using insert_vertex() and insert_edge() @@ -235,13 +246,14 @@ void PTQ_1step(Poly& p, VertexPointMap vpm, Mask mask) { // 4. insert_edge() between all other new inserted vertices of step 1 and // the new inserted vertex of step 3 // Step 1. - edge_iterator eitr = edges(p).first; - for (size_t i = 0; i < num_edge; i++, ++eitr) { - vertex_descriptor vh = insert_vertex(p, halfedge(*eitr,p)); + typename std::vector::iterator eit = p_edges.begin(); + for (std::size_t i = 0; i < num_e; i++, ++eit) { + vertex_descriptor vh = insert_vertex(p, halfedge(*eit,p)); put(vpm,vh, edge_point_buffer[i]); } - face_iterator fitr = faces(p).first; - for (size_t i = 0; i < num_facet; i++, ++fitr) { + + typename std::vector::iterator fitr = p_faces.begin(); + for (std::size_t i = 0; i < num_f; i++, ++fitr) { // Step 2. Halfedge_around_face_circulator hcir_begin(halfedge(*fitr,p),p); Halfedge_around_face_circulator hcir = hcir_begin; @@ -261,8 +273,8 @@ void PTQ_1step(Poly& p, VertexPointMap vpm, Mask mask) { // Update the geometry data of the newly inserted vertices by the // vertices buffer - vitr = vertices(p).first; - for (size_t i = 0; i < num_vertex; i++, ++vitr) + vitr = p_vertices.begin(); + for (std::size_t i = 0; i < num_v; i++, ++vitr) put(vpm, *vitr, vertex_point_buffer[i]); // CGAL_postcondition(p.is_valid()); @@ -272,54 +284,60 @@ void PTQ_1step(Poly& p, VertexPointMap vpm, Mask mask) { // ====================================================================== template -void DQQ_1step_impl(Poly& p, VertexPointMap vpm, Mask mask, CGAL::Tag_false) { +void DQQ_1step(Poly& p, VertexPointMap vpm, Mask mask) { typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; typedef typename boost::graph_traits::edge_descriptor edge_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; - typedef typename boost::graph_traits::vertex_iterator vertex_iterator; - typedef typename boost::graph_traits::edge_iterator edge_iterator; typedef typename boost::property_traits::value_type Point; - typename boost::graph_traits::vertices_size_type num_v = num_vertices(p); - typename boost::graph_traits::halfedges_size_type num_e = num_halfedges(p)/2; - typename boost::graph_traits::faces_size_type num_f = num_faces(p); + //First back up initial vertices/faces/edges + std::vector p_vertices(vertices(p).first, vertices(p).second); + std::vector p_faces(faces(p).first, faces(p).second); + std::vector p_edges(edges(p).first, edges(p).second); + + std::size_t num_v = p_vertices.size(); + std::size_t num_e = p_edges.size(); + std::size_t num_f = p_faces.size(); std::vector border_halfedges; - size_t num_be = 0 ; - BOOST_FOREACH(edge_descriptor ed, edges(p)){ + BOOST_FOREACH(edge_descriptor ed, p_edges){ if(is_border(ed,p)){ - ++num_be; border_halfedges.push_back(halfedge(ed,p)); } } Point* point_buffer = new Point[num_e*2]; // Build the point_buffer - vertex_iterator vitr, vitr_end; - boost::tie(vitr,vitr_end) = vertices(p); int pi = 0; - BOOST_FOREACH(vertex_descriptor vd, vertices(p)){ + int vi=0; + std::vector is_border_vertex(num_v, false); + BOOST_FOREACH(vertex_descriptor vd, p_vertices){ BOOST_FOREACH(halfedge_descriptor hd, halfedges_around_target(vd,p)){ if (! is_border(hd,p)){ mask.corner_node(hd, point_buffer[pi++]); } + else + is_border_vertex[vi]=true; } + ++vi; } // Reserve to avoid rellocations during insertions - reserve(p,num_v+num_e+num_f, 2*num_e, (2+4+2)*num_e); + call_reserve(p,num_v+num_e+num_f, 2*num_e, (2+4+2)*num_e); // Build the connectivity using insert_vertex() and insert_edge() pi = 0; - for (size_t i = 0; i < num_v; ++i) { + typename std::vector::iterator vitr=p_vertices.begin(); + for (std::size_t i = 0; i < num_v; ++i) { vertex_descriptor vh = *vitr; ++vitr; Halfedge_around_target_circulator vcir(vh,p); - size_t vn = degree(vh,p); - for (size_t j = 0; j < vn; ++j) { + typename boost::graph_traits::degree_size_type vn = degree(vh,p); + for (typename boost::graph_traits::degree_size_type j = 0; j < vn; ++j) { halfedge_descriptor e = *vcir; ++vcir; if (! is_border(e,p)) { @@ -329,7 +347,7 @@ void DQQ_1step_impl(Poly& p, VertexPointMap vpm, Mask mask, CGAL::Tag_false) { } vcir = Halfedge_around_target_circulator(vh,p); - for (size_t j = 0; j < vn; ++j) { + for (typename boost::graph_traits::vertices_size_type j = 0; j < vn; ++j) { if (! is_border(*vcir,p)) { halfedge_descriptor e1 = prev(*vcir, p); ++vcir; @@ -343,8 +361,8 @@ void DQQ_1step_impl(Poly& p, VertexPointMap vpm, Mask mask, CGAL::Tag_false) { } } - edge_iterator eitr = edges(p).first; - for (size_t i = 0; i < num_e; ++i) { + typename std::vector::iterator eitr = p_edges.begin(); + for (typename boost::graph_traits::edges_size_type i = 0; i < num_e; ++i) { halfedge_descriptor eh = halfedge(*eitr,p); ++eitr; if (! is_border(edge(eh,p),p)) { @@ -361,7 +379,7 @@ void DQQ_1step_impl(Poly& p, VertexPointMap vpm, Mask mask, CGAL::Tag_false) { } } - // After this point, the original border edges are in front! + BOOST_FOREACH(halfedge_descriptor eeh, border_halfedges){ halfedge_descriptor eh = eeh; if (is_border(eh,p)){ @@ -378,169 +396,17 @@ void DQQ_1step_impl(Poly& p, VertexPointMap vpm, Mask mask, CGAL::Tag_false) { Euler::remove_face(ehe,p); } - vitr = vertices(p).first; - for (size_t i = 0; i < num_v-num_be; ++i) { + vitr = p_vertices.begin(); + for (typename boost::graph_traits::vertices_size_type i = 0; i < num_v; ++i) { vertex_descriptor vh = *vitr; ++vitr; - Euler::remove_center_vertex(halfedge(vh,p),p); + if (!is_border_vertex[i]) + Euler::remove_center_vertex(halfedge(vh,p),p); } delete []point_buffer; } -template -void DQQ_1step_impl(Poly& p, VertexPointMap vpm, Mask mask, CGAL::Tag_true) { - typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; - typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; - typedef typename boost::graph_traits::face_descriptor face_descriptor; - - typedef typename boost::property_traits::value_type Point; - - // Note that for types like 'Surface_mesh', num_vertices() and other similar functions - // return the TOTAL number of elements, which may include removed vertices. - typename boost::graph_traits::vertices_size_type num_v = num_vertices(p); - typename boost::graph_traits::edges_size_type num_e = num_edges(p); - typename boost::graph_traits::faces_size_type num_f = num_faces(p); - - // Move `p` into `moved_p`, and build the subdivided mesh from scratch in `p`. - // This is done to make the algorithm work with CGAL::Surface_mesh, - // even though CGAL::Surface_mesh does not insert elements at the end (due to removed elements). - // The DooSabin subdivision of a mesh is a completely different mesh so there - // is no additional cost to rebuild from scratch (but there is a bit from - // using `copy_face_graph`). - Poly moved_p; - reserve(moved_p,num_v, num_e, num_f); - - // We must use copy_face_graph rather than an assignement operator because - // we need the correspondence between vertex_descriptors - boost::unordered_map v2v(num_v); - CGAL::copy_face_graph(p, moved_p, std::inserter(v2v, v2v.end())); - - VertexPointMap moved_vpm = get(vertex_point, moved_p); - - // Move the position information to the internal property map of moved_p - typename boost::unordered_map::iterator it = v2v.begin(), - end = v2v.end(); - for(; it!=end; ++it) { - put(moved_vpm, it->second, get(vpm, it->first)); - } - - // Temporarily change the members of the mask to `moved_p` - mask.pmesh = &moved_p; - mask.vpm = moved_vpm; - - clear(p); - reserve(p,num_v+num_e+num_f, 2*num_e, (2+4+2)*num_e); - - // Correspondence between halfedges of the original mesh and some of the - // halfedges in the subdivided mesh. Since we have halfedge_index_t, - // we can simply use a vector! - // Note: need to make sure the halfedge index map is initialized - // @fixme this overwrites previous initilizations... - helpers::init_halfedge_indices(const_cast(moved_p), - get(boost::halfedge_index, moved_p)); - std::vector old_to_new(2 * num_e); - Property_map_binder::type, - typename Pointer_property_map::type> - hmap = bind_property_maps(get(boost::halfedge_index, moved_p), - make_property_map(old_to_new)); - - // Build new n-faces - BOOST_FOREACH(face_descriptor fd, faces(moved_p)) { - halfedge_descriptor hd = halfedge(fd, moved_p); - std::list vertices_of_new_face; - - // Keep the first outside; it will be used to build the correspondence - // between old and new halfedges - Point first_pt; - mask.corner_node(hd, first_pt); - vertex_descriptor first_v = add_vertex(p); - put(vpm, first_v, first_pt); - vertices_of_new_face.push_back(first_v); - - // Loop normally and add the rest of the vertices - halfedge_descriptor done = hd; - hd = next(hd, moved_p); - while(hd != done) { - Point pt; - mask.corner_node(hd, pt); - vertex_descriptor v = add_vertex(p); - put(vpm, v, pt); - vertices_of_new_face.push_back(v); - hd = next(hd, moved_p); - } - - face_descriptor new_face = Euler::add_face(vertices_of_new_face, p); - - // Find the starting halfedge in the new face that corresponds to halfedge(fd, p) - halfedge_descriptor nf_hd = halfedge(new_face, p); - while(target(nf_hd, p) != first_v) { - nf_hd = next(nf_hd, p); - } - - // Build the correspondence between old and new halfedges - hd = halfedge(fd, moved_p); - done = nf_hd; - do { - put(hmap, hd, nf_hd); - hd = next(hd, moved_p); - nf_hd = next(nf_hd, p); - } while (nf_hd != done); - } - - // Build new edge-faces - BOOST_FOREACH(halfedge_descriptor hd, halfedges(moved_p)) { - if(is_border(hd, moved_p)) - continue; - - halfedge_descriptor hd_opp = opposite(hd, moved_p); - if(is_border(hd_opp, moved_p)) - continue; - - if(hd > hd_opp) - continue; - - halfedge_descriptor new_hd = opposite(get(hmap, hd), p); - halfedge_descriptor new_hd_opp = opposite(get(hmap, hd_opp), p); - - boost::array v = {{source(new_hd, p), - target(new_hd, p), - source(new_hd_opp, p), - target(new_hd_opp, p)}}; - - Euler::add_face(v, p); - } - - // Build new vertex-faces - BOOST_FOREACH(vertex_descriptor vd, vertices(moved_p)) { - if(is_border(vd, moved_p)) - continue; - - halfedge_descriptor hd = halfedge(vd, moved_p); - halfedge_descriptor new_hd = opposite(get(hmap, hd), p); - halfedge_descriptor new_face_hd = opposite(prev(new_hd, p), p), done = new_face_hd; - std::list vertices_of_new_faces; - do { - vertices_of_new_faces.push_back(source(new_face_hd, p)); - new_face_hd = next(new_face_hd, p); - } while(new_face_hd != done); - - Euler::add_face(vertices_of_new_faces, p); - } - - // Reset the members of the mask - mask.pmesh = &p; - mask.vpm = vpm; -} - -template -void DQQ_1step(Poly& p, VertexPointMap vpm, Mask mask) { - // Check if halfedges are index-based, which allows to use vectors instead of maps - DQQ_1step_impl(p, vpm, mask, - boost::graph_has_property()); -// CGAL_postcondition(p.is_valid()); -} - // ====================================================================== template void Sqrt3_1step(Poly& p, VertexPointMap vpm, Mask mask, @@ -556,19 +422,21 @@ void Sqrt3_1step(Poly& p, VertexPointMap vpm, Mask mask, typedef typename boost::graph_traits::edge_descriptor edge_descriptor; typedef typename boost::graph_traits::face_descriptor face_descriptor; - typedef typename boost::graph_traits::edge_iterator edge_iterator; - typedef typename boost::graph_traits::face_iterator face_iterator; - typedef typename boost::property_traits::value_type Point; - typename boost::graph_traits::vertices_size_type num_v = num_vertices(p); - typename boost::graph_traits::halfedges_size_type num_e = num_halfedges(p)/2; - typename boost::graph_traits::faces_size_type num_f = num_faces(p); + //First back up initial vertices/faces/edges + std::vector p_vertices(vertices(p).first, vertices(p).second); + std::vector p_faces(faces(p).first, faces(p).second); + std::vector p_edges(edges(p).first, edges(p).second); + + std::size_t num_v = p_vertices.size(); + std::size_t num_e = p_edges.size(); + std::size_t num_f = p_faces.size(); // reserve enough size for the new points - typename boost::graph_traits::faces_size_type new_pts_size = num_f; + std::size_t new_pts_size = num_f; if(refine_border) { - BOOST_FOREACH(edge_descriptor ed, edges(p)){ + BOOST_FOREACH(edge_descriptor ed, p_edges){ if(is_border(ed, p)) ++new_pts_size; } @@ -576,7 +444,7 @@ void Sqrt3_1step(Poly& p, VertexPointMap vpm, Mask mask, Point* cpt = new Point[new_pts_size]; // size of the subdivided mesh - reserve(p,num_v + new_pts_size, (num_e + 2*num_f + new_pts_size)*2, 3*num_f); + call_reserve(p,num_v + new_pts_size, (num_e + 2*num_f + new_pts_size)*2, 3*num_f); // keep in memory whether a face is incident to the border and, if so, which // halfedge corresponds to THE (there can only be one) border edge. @@ -587,7 +455,7 @@ void Sqrt3_1step(Poly& p, VertexPointMap vpm, Mask mask, // compute the positions of new points std::size_t i = 0; std::size_t face_id = 0; - BOOST_FOREACH (face_descriptor fd, faces(p)) { + BOOST_FOREACH (face_descriptor fd, p_faces) { //ASSERTION_MSG(circulator_size(fitr->facet_begin())==3, "(ERROR) Non-triangle facet!"); if(refine_border) { BOOST_FOREACH(halfedge_descriptor hd, halfedges_around_face(halfedge(fd, p), p)) { @@ -614,7 +482,7 @@ void Sqrt3_1step(Poly& p, VertexPointMap vpm, Mask mask, } // smooth the position of existing vertices - BOOST_FOREACH(vertex_descriptor vd, vertices(p)){ + BOOST_FOREACH(vertex_descriptor vd, p_vertices){ Point pt; if(!is_border(vd, p)) { mask.vertex_node(vd, pt); @@ -623,10 +491,9 @@ void Sqrt3_1step(Poly& p, VertexPointMap vpm, Mask mask, } // insert the new subdividing points - face_iterator b,e; - boost::tie(b,e) = faces(p); - for(std::size_t i=0, cpt_id=0; i < num_f; ++i, ++b){ - face_descriptor fd = *b; + typename std::vector::iterator fit=p_faces.begin(); + for(std::size_t i=0, cpt_id=0; i < num_f; ++i, ++fit){ + face_descriptor fd = *fit; halfedge_descriptor hd = face_halfedge_border[i]; if(refine_border && hd != boost::graph_traits::null_halfedge()) { halfedge_descriptor hd_next = next(hd, p); @@ -650,8 +517,8 @@ void Sqrt3_1step(Poly& p, VertexPointMap vpm, Mask mask, } // flip the old edges (except the border edges) - edge_iterator eitr = edges(p).first; - for (size_t i = 0; i < num_e; ++i) { + typename std::vector::iterator eitr = p_edges.begin(); + for (std::size_t i = 0; i < num_e; ++i) { halfedge_descriptor e = halfedge(*eitr,p); ++eitr; // move to next edge before flip since flip destroys current edge if (! is_border(edge(e,p),p)) { diff --git a/Subdivision_method_3/include/CGAL/Subdivision_method_3/subdivision_masks_3.h b/Subdivision_method_3/include/CGAL/Subdivision_method_3/subdivision_masks_3.h index 4fe543ad2c2..da9ad657d87 100644 --- a/Subdivision_method_3/include/CGAL/Subdivision_method_3/subdivision_masks_3.h +++ b/Subdivision_method_3/include/CGAL/Subdivision_method_3/subdivision_masks_3.h @@ -93,6 +93,8 @@ public: typedef typename Base::FT FT; typedef typename Base::Point Point; typedef typename Base::Vector Vector; + + typedef typename boost::property_traits::reference Point_ref; #endif public: @@ -117,8 +119,8 @@ public: } void edge_node(halfedge_descriptor edge, Point& pt) { - const Point& p1 = get(this->vpmap, target(edge, *(this->pmesh))); - const Point& p2 = get(this->vpmap, source(edge, *(this->pmesh))); + Point_ref p1 = get(this->vpmap, target(edge, *(this->pmesh))); + Point_ref p2 = get(this->vpmap, source(edge, *(this->pmesh))); pt = Point((p1[0]+p2[0])/2, (p1[1]+p2[1])/2, (p1[2]+p2[2])/2); } @@ -173,6 +175,8 @@ public: typedef typename Base::FT FT; typedef typename Base::Point Point; typedef typename Base::Vector Vector; + + typedef typename boost::property_traits::reference Point_ref; #endif public: @@ -197,8 +201,8 @@ public: /// computes the Catmull-Clark edge-point `pt` of the edge `edge`. void edge_node(halfedge_descriptor edge, Point& pt) { - const Point& p1 = get(this->vpmap,target(edge, *(this->pmesh))); - const Point& p2 = get(this->vpmap,source(edge, *(this->pmesh))); + Point_ref p1 = get(this->vpmap,target(edge, *(this->pmesh))); + Point_ref p2 = get(this->vpmap,source(edge, *(this->pmesh))); Point f1, f2; this->face_node(face(edge, *(this->pmesh)), f1); this->face_node(face(opposite(edge, *(this->pmesh)), *(this->pmesh)), f2); @@ -213,11 +217,11 @@ public: typename boost::graph_traits::degree_size_type n = degree(vertex, *(this->pmesh)); FT Q[] = {0.0, 0.0, 0.0}, R[] = {0.0, 0.0, 0.0}; - Point& S = get(this->vpmap,vertex); + Point_ref S = get(this->vpmap,vertex); Point q; - for (unsigned int i = 0; i < n; i++, ++vcir) { - const Point& p2 = get(this->vpmap, target(opposite(*vcir, *(this->pmesh)), *(this->pmesh))); + for (typename boost::graph_traits::degree_size_type i = 0; i < n; i++, ++vcir) { + Point_ref p2 = get(this->vpmap, target(opposite(*vcir, *(this->pmesh)), *(this->pmesh))); R[0] += (S[0] + p2[0]) / 2; R[1] += (S[1] + p2[1]) / 2; R[2] += (S[2] + p2[2]) / 2; @@ -237,15 +241,15 @@ public: /// computes the Catmull-Clark edge-point `ept` and the Catmull-Clark /// vertex-point `vpt` of the border edge `edge`. void border_node(halfedge_descriptor edge, Point& ept, Point& vpt) { - const Point& ep1 = get(this->vpmap,target(edge, *(this->pmesh))); - const Point& ep2 = get(this->vpmap,target(opposite(edge, *(this->pmesh)), *(this->pmesh))); + Point_ref ep1 = get(this->vpmap,target(edge, *(this->pmesh))); + Point_ref ep2 = get(this->vpmap,target(opposite(edge, *(this->pmesh)), *(this->pmesh))); ept = Point((ep1[0]+ep2[0])/2, (ep1[1]+ep2[1])/2, (ep1[2]+ep2[2])/2); Halfedge_around_target_circulator vcir(edge, *(this->pmesh)); - const Point& vp1 = get(this->vpmap,target(opposite(*vcir, *(this->pmesh)), *(this->pmesh))); - const Point& vp0 = get(this->vpmap, target(*vcir, *(this->pmesh))); + Point_ref vp1 = get(this->vpmap,target(opposite(*vcir, *(this->pmesh)), *(this->pmesh))); + Point_ref vp0 = get(this->vpmap, target(*vcir, *(this->pmesh))); --vcir; - const Point& vp_1 = get(this->vpmap, target(opposite(*vcir, *(this->pmesh)), *(this->pmesh))); + Point_ref vp_1 = get(this->vpmap, target(opposite(*vcir, *(this->pmesh)), *(this->pmesh))); vpt = Point((vp_1[0] + 6*vp0[0] + vp1[0])/8, (vp_1[1] + 6*vp0[1] + vp1[1])/8, (vp_1[2] + 6*vp0[2] + vp1[2])/8 ); @@ -296,6 +300,7 @@ public: typedef typename Base::FT FT; typedef typename Base::Point Point; typedef typename Base::Vector Vector; + typedef typename boost::property_traits::reference Point_ref; #endif typedef Halfedge_around_face_circulator Halfedge_around_facet_circulator; @@ -323,10 +328,10 @@ public: /// computes the Loop edge-point `pt` of the edge `edge`. void edge_node(halfedge_descriptor edge, Point& pt) { - const Point& p1 = get(this->vpmap,target(edge, *(this->pmesh))); - const Point& p2 = get(this->vpmap,target(opposite(edge, *(this->pmesh)), *(this->pmesh))); - const Point& f1 = get(this->vpmap,target(next(edge, *(this->pmesh)), *(this->pmesh))); - const Point& f2 = get(this->vpmap,target(next(opposite(edge, *(this->pmesh)), *(this->pmesh)), *(this->pmesh))); + Point_ref p1 = get(this->vpmap,target(edge, *(this->pmesh))); + Point_ref p2 = get(this->vpmap,target(opposite(edge, *(this->pmesh)), *(this->pmesh))); + Point_ref f1 = get(this->vpmap,target(next(edge, *(this->pmesh)), *(this->pmesh))); + Point_ref f2 = get(this->vpmap,target(next(opposite(edge, *(this->pmesh)), *(this->pmesh)), *(this->pmesh))); pt = Point((3*(p1[0]+p2[0])+f1[0]+f2[0])/8, (3*(p1[1]+p2[1])+f1[1]+f2[1])/8, @@ -339,10 +344,10 @@ public: size_t n = circulator_size(vcir); FT R[] = {0.0, 0.0, 0.0}; - const Point& S = get(this->vpmap,vertex); + Point_ref S = get(this->vpmap,vertex); for (size_t i = 0; i < n; i++, ++vcir) { - const Point& p = get(this->vpmap,target(opposite(*vcir, *(this->pmesh)), *(this->pmesh))); + Point_ref p = get(this->vpmap,target(opposite(*vcir, *(this->pmesh)), *(this->pmesh))); R[0] += p[0]; R[1] += p[1]; R[2] += p[2]; } if (n == 6) { @@ -361,15 +366,15 @@ public: /// computes the Loop edge-point `ept` and the Loop vertex-point `vpt` of the border edge `edge`. void border_node(halfedge_descriptor edge, Point& ept, Point& vpt) { - const Point& ep1 = get(this->vpmap,target(edge, *(this->pmesh))); - const Point& ep2 = get(this->vpmap,target(opposite(edge, *(this->pmesh)), *(this->pmesh))); + Point_ref ep1 = get(this->vpmap,target(edge, *(this->pmesh))); + Point_ref ep2 = get(this->vpmap,target(opposite(edge, *(this->pmesh)), *(this->pmesh))); ept = Point((ep1[0]+ep2[0])/2, (ep1[1]+ep2[1])/2, (ep1[2]+ep2[2])/2); Halfedge_around_vertex_circulator vcir(edge, *(this->pmesh)); - const Point& vp1 = get(this->vpmap,target(opposite(*vcir, *(this->pmesh) ), *(this->pmesh))); - const Point& vp0 = get(this->vpmap,target(*vcir, *(this->pmesh))); + Point_ref vp1 = get(this->vpmap,target(opposite(*vcir, *(this->pmesh) ), *(this->pmesh))); + Point_ref vp0 = get(this->vpmap,target(*vcir, *(this->pmesh))); --vcir; - const Point& vp_1 = get(this->vpmap,target(opposite(*vcir, *(this->pmesh)), *(this->pmesh))); + Point_ref vp_1 = get(this->vpmap,target(opposite(*vcir, *(this->pmesh)), *(this->pmesh))); vpt = Point((vp_1[0] + 6*vp0[0] + vp1[0])/8, (vp_1[1] + 6*vp0[1] + vp1[1])/8, (vp_1[2] + 6*vp0[2] + vp1[2])/8 ); @@ -457,6 +462,7 @@ public: typedef typename Base::FT FT; typedef typename Base::Point Point; typedef typename Base::Vector Vector; + typedef typename boost::property_traits::reference Point_ref; #endif public: @@ -574,12 +580,12 @@ public: /// computes the \f$ \sqrt{3}\f$ vertex-point `pt` of the vertex `vd`. void vertex_node(vertex_descriptor vertex, Point& pt) { Halfedge_around_target_circulator vcir(vertex, *(this->pmesh)); - const size_t n = degree(vertex, *(this->pmesh)); + const typename boost::graph_traits::degree_size_type n = degree(vertex, *(this->pmesh)); const FT a = (FT) ((4.0-2.0*std::cos(2.0*CGAL_PI/(double)n))/9.0); Vector cv = ((FT)(1.0-a)) * (get(this->vpmap, vertex) - CGAL::ORIGIN); - for (size_t i = 1; i <= n; ++i, --vcir) { + for (typename boost::graph_traits::degree_size_type i = 1; i <= n; ++i, --vcir) { cv = cv + (a/FT(n))*(get(this->vpmap, target(opposite(*vcir, *(this->pmesh)), *(this->pmesh)))-CGAL::ORIGIN); } diff --git a/Subdivision_method_3/test/Subdivision_method_3/test_Subdivision_method_3.cpp b/Subdivision_method_3/test/Subdivision_method_3/test_Subdivision_method_3.cpp index b0b18514ded..9250848883a 100644 --- a/Subdivision_method_3/test/Subdivision_method_3/test_Subdivision_method_3.cpp +++ b/Subdivision_method_3/test/Subdivision_method_3/test_Subdivision_method_3.cpp @@ -338,7 +338,7 @@ void test_Subdivision_surface_3_SM_NP() { // some arbitrary new coordinates (different from the internal vpm) BOOST_FOREACH(vertex_descriptor vd, vertices(P)) { - const Point& pt = get(vpm, vd); + boost::property_traits::reference pt = get(vpm, vd); Vector v = pt - Point(0., 0., -3.); put(apm, vd, pt + 0.5*v); } diff --git a/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/Surface_mesh_shortest_path.h b/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/Surface_mesh_shortest_path.h index 7ccb8bb27a3..f4b1c32fe8e 100644 --- a/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/Surface_mesh_shortest_path.h +++ b/Surface_mesh_shortest_path/include/CGAL/Surface_mesh_shortest_path/Surface_mesh_shortest_path.h @@ -1269,7 +1269,7 @@ private: propagateRight = rightSide; } - if (node->level() <= num_faces(m_graph)) + if (node->level() <= static_cast(num_faces(m_graph))) { if (propagateLeft) { diff --git a/Triangulation_2/include/CGAL/Constrained_Delaunay_triangulation_2.h b/Triangulation_2/include/CGAL/Constrained_Delaunay_triangulation_2.h index 417d0e5808f..f56addb6708 100644 --- a/Triangulation_2/include/CGAL/Constrained_Delaunay_triangulation_2.h +++ b/Triangulation_2/include/CGAL/Constrained_Delaunay_triangulation_2.h @@ -750,9 +750,9 @@ propagating_flip(Face_handle f,int i, int depth) Face_handle ni = f->neighbor(i); flip(f, i); // flip for constrained triangulations - propagating_flip(f,i); + propagating_flip(f,i, depth+1); i = ni->index(f->vertex(i)); - propagating_flip(ni,i); + propagating_flip(ni,i, depth+1); #endif } #else diff --git a/Triangulation_2/include/CGAL/Constrained_triangulation_2.h b/Triangulation_2/include/CGAL/Constrained_triangulation_2.h index 9f45f31305d..1a11c13179b 100644 --- a/Triangulation_2/include/CGAL/Constrained_triangulation_2.h +++ b/Triangulation_2/include/CGAL/Constrained_triangulation_2.h @@ -691,45 +691,53 @@ insert_constraint(Vertex_handle vaa, Vertex_handle vbb) // if a vertex vc of t lies on segment ab // or if ab intersect some constrained edges { - CGAL_triangulation_precondition( vaa != vbb); - Vertex_handle vi; + std::stack > stack; + stack.push(std::make_pair(vaa,vbb)); - Face_handle fr; - int i; - if(includes_edge(vaa,vbb,vi,fr,i)) { - mark_constraint(fr,i); - if (vi != vbb) { - insert_constraint(vi,vbb); + while(! stack.empty()){ + boost::tie(vaa,vbb) = stack.top(); + stack.pop(); + CGAL_triangulation_precondition( vaa != vbb); + Vertex_handle vi; + + Face_handle fr; + int i; + if(includes_edge(vaa,vbb,vi,fr,i)) { + mark_constraint(fr,i); + if (vi != vbb) { + stack.push(std::make_pair(vi,vbb)); + } + continue; } - return; - } - List_faces intersected_faces; - List_edges conflict_boundary_ab, conflict_boundary_ba; + List_faces intersected_faces; + List_edges conflict_boundary_ab, conflict_boundary_ba; - bool intersection = find_intersected_faces( vaa, vbb, - intersected_faces, - conflict_boundary_ab, - conflict_boundary_ba, - vi); - if ( intersection) { - if (vi != vaa && vi != vbb) { - insert_constraint(vaa,vi); - insert_constraint(vi,vbb); - } - else insert_constraint(vaa,vbb); - return; - } + bool intersection = find_intersected_faces( vaa, vbb, + intersected_faces, + conflict_boundary_ab, + conflict_boundary_ba, + vi); + if ( intersection) { + if (vi != vaa && vi != vbb) { + stack.push(std::make_pair(vaa,vi)); + stack.push(std::make_pair(vi,vbb)); + } + else{ + stack.push(std::make_pair(vaa,vbb)); + } + continue; + } - //no intersection - triangulate_hole(intersected_faces, - conflict_boundary_ab, - conflict_boundary_ba); + //no intersection + triangulate_hole(intersected_faces, + conflict_boundary_ab, + conflict_boundary_ba); - if (vi != vbb) { - insert_constraint(vi,vbb); + if (vi != vbb) { + stack.push(std::make_pair(vi,vbb)); + } } - return; } diff --git a/Triangulation_2/include/CGAL/Constrained_triangulation_plus_2.h b/Triangulation_2/include/CGAL/Constrained_triangulation_plus_2.h index b36f30363b9..0f98f9c67d4 100644 --- a/Triangulation_2/include/CGAL/Constrained_triangulation_plus_2.h +++ b/Triangulation_2/include/CGAL/Constrained_triangulation_plus_2.h @@ -781,82 +781,88 @@ insert_subconstraint(Vertex_handle vaa, // insert the subconstraint [vaa vbb] // it will eventually be splitted into several subconstraints { - CGAL_triangulation_precondition( vaa != vbb); - Vertex_handle vi; + std::stack > stack; + stack.push(std::make_pair(vaa,vbb)); - Face_handle fr; - int i; - if(this->includes_edge(vaa,vbb,vi,fr,i)) { - this->mark_constraint(fr,i); - if (vi != vbb) { - hierarchy.split_constraint(vaa,vbb,vi); - insert_subconstraint(vi,vbb, out); + while(! stack.empty()){ + boost::tie(vaa,vbb) = stack.top(); + stack.pop(); + CGAL_triangulation_precondition( vaa != vbb); + + Vertex_handle vi; + + Face_handle fr; + int i; + if(this->includes_edge(vaa,vbb,vi,fr,i)) { + this->mark_constraint(fr,i); + if (vi != vbb) { + hierarchy.split_constraint(vaa,vbb,vi); + stack.push(std::make_pair(vi,vbb)); + } + continue; } - return; - } - List_faces intersected_faces; - List_edges conflict_boundary_ab, conflict_boundary_ba; + List_faces intersected_faces; + List_edges conflict_boundary_ab, conflict_boundary_ba; - bool intersection = this->find_intersected_faces( - vaa, vbb, - intersected_faces, - conflict_boundary_ab, - conflict_boundary_ba, - vi); + bool intersection = this->find_intersected_faces( + vaa, vbb, + intersected_faces, + conflict_boundary_ab, + conflict_boundary_ba, + vi); - if ( intersection) { - if (vi != vaa && vi != vbb) { - hierarchy.split_constraint(vaa,vbb,vi); - insert_subconstraint(vaa,vi, out); - insert_subconstraint(vi,vbb, out); - } - else insert_subconstraint(vaa,vbb,out); + if ( intersection) { + if (vi != vaa && vi != vbb) { + hierarchy.split_constraint(vaa,vbb,vi); + stack.push(std::make_pair(vaa,vi)); + stack.push(std::make_pair(vi,vbb)); + } + else stack.push(std::make_pair(vaa,vbb)); - - return; - } + continue; + } - //no intersection + //no intersection - List_edges edges(conflict_boundary_ab); - std::copy(conflict_boundary_ba.begin(), conflict_boundary_ba.end(), std::back_inserter(edges)); + List_edges edges(conflict_boundary_ab); + std::copy(conflict_boundary_ba.begin(), conflict_boundary_ba.end(), std::back_inserter(edges)); - // edges may contain mirror edges. They no longer exist after triangulate_hole - // so we have to remove them before calling get_bounded_faces - if(! edges.empty()){ + // edges may contain mirror edges. They no longer exist after triangulate_hole + // so we have to remove them before calling get_bounded_faces + if(! edges.empty()){ #if defined(BOOST_MSVC) && (BOOST_VERSION == 105500) - std::set faces(intersected_faces.begin(), intersected_faces.end()); + std::set faces(intersected_faces.begin(), intersected_faces.end()); #else - boost::container::flat_set faces(intersected_faces.begin(), intersected_faces.end()); + boost::container::flat_set faces(intersected_faces.begin(), intersected_faces.end()); #endif - typename List_edges::iterator it2; - for(typename List_edges::iterator it = edges.begin(); it!= edges.end();){ - if(faces.find(it->first) != faces.end()){ - typename List_edges::iterator it2 = it; - ++it; - edges.erase(it2); - }else { - ++it; + typename List_edges::iterator it2; + for(typename List_edges::iterator it = edges.begin(); it!= edges.end();){ + if(faces.find(it->first) != faces.end()){ + typename List_edges::iterator it2 = it; + ++it; + edges.erase(it2); + }else { + ++it; + } } } + + this->triangulate_hole(intersected_faces, + conflict_boundary_ab, + conflict_boundary_ba); + + this->get_bounded_faces(edges.begin(), + edges.end(), + out); + + if (vi != vbb) { + hierarchy.split_constraint(vaa,vbb,vi); + stack.push(std::make_pair(vi,vbb)); + } } - - this->triangulate_hole(intersected_faces, - conflict_boundary_ab, - conflict_boundary_ba); - - this->get_bounded_faces(edges.begin(), - edges.end(), - out); - - if (vi != vbb) { - hierarchy.split_constraint(vaa,vbb,vi); - insert_subconstraint(vi,vbb, out); - } - return; }