diff --git a/Polygon_mesh_processing/include/CGAL/internal/Hole_filling/Fair_Polyhedron_3.h b/Polygon_mesh_processing/include/CGAL/internal/Hole_filling/Fair_Polyhedron_3.h index 3c05cf8d59d..56d6be15aaf 100644 --- a/Polygon_mesh_processing/include/CGAL/internal/Hole_filling/Fair_Polyhedron_3.h +++ b/Polygon_mesh_processing/include/CGAL/internal/Hole_filling/Fair_Polyhedron_3.h @@ -60,13 +60,17 @@ class Fair_Polyhedron_3 { typedef SparseLinearSolver Sparse_linear_solver; // members + Polyhedron& polyhedron; WeightCalculator weight_calculator; + Point_property_map ppmap; + public: - Fair_Polyhedron_3(WeightCalculator weight_calculator = WeightCalculator()) - : weight_calculator(weight_calculator) { } + Fair_Polyhedron_3(Polyhedron& polyhedron, WeightCalculator weight_calculator = WeightCalculator()) + : polyhedron(polyhedron), weight_calculator(weight_calculator), ppmap(get(CGAL::vertex_point, polyhedron)) + { } private: - double sum_weight(Vertex_handle v, Polyhedron& polyhedron) { + double sum_weight(Vertex_handle v) { double weight = 0; Halfedge_around_vertex_circulator circ(halfedge(v,polyhedron),polyhedron), done(circ); do { @@ -79,7 +83,6 @@ private: // Equation 6 in [On Linear Variational Surface Deformation Methods] void compute_row( Vertex_handle v, - Polyhedron& polyhedron, std::size_t row_id, // which row to insert in [ frees stay left-hand side ] typename Sparse_linear_solver::Matrix& matrix, double& x, double& y, double& z, // constants transfered to right-hand side @@ -93,30 +96,31 @@ private: matrix.add_coef(row_id, vertex_id_it->second, multiplier); } else { - x += multiplier * -v->point().x(); - y += multiplier * -v->point().y(); - z += multiplier * -v->point().z(); + Point_3& p = ppmap[v]; + x += multiplier * - p.x(); + y += multiplier * - p.y(); + z += multiplier * - p.z(); } } else { - double w_i = weight_calculator.w_i(v, polyhedron); + double w_i = weight_calculator.w_i(v,polyhedron); Halfedge_around_vertex_circulator circ(halfedge(v,polyhedron),polyhedron), done(circ); do { double w_i_w_ij = w_i * weight_calculator.w_ij(*circ, polyhedron) ; Vertex_handle nv = target(opposite(*circ,polyhedron),polyhedron); - compute_row(nv, polyhedron, row_id, matrix, x, y, z, -w_i_w_ij*multiplier, vertex_id_map, depth-1); + compute_row(nv, row_id, matrix, x, y, z, -w_i_w_ij*multiplier, vertex_id_map, depth-1); } while(++circ != done); - double w_i_w_ij_sum = w_i * sum_weight(v, polyhedron); - compute_row(v, polyhedron, row_id, matrix, x, y, z, w_i_w_ij_sum*multiplier, vertex_id_map, depth-1); + double w_i_w_ij_sum = w_i * sum_weight(v); + compute_row(v, row_id, matrix, x, y, z, w_i_w_ij_sum*multiplier, vertex_id_map, depth-1); } } public: template - bool fair(Polyhedron& polyhedron, InputIterator vb, InputIterator ve, Fairing_continuity fc) + bool fair(InputIterator vb, InputIterator ve, Fairing_continuity fc) { int depth = static_cast(fc) + 1; if(depth < 1 || depth > 3) { @@ -147,7 +151,7 @@ public: for(typename std::set::iterator vb = interior_vertices.begin(); vb != interior_vertices.end(); ++vb) { std::size_t v_id = vertex_id_map[*vb]; - compute_row(*vb, polyhedron, v_id, A, Bx[v_id], By[v_id], Bz[v_id], 1, vertex_id_map, depth); + compute_row(*vb, v_id, A, Bx[v_id], By[v_id], Bz[v_id], 1, vertex_id_map, depth); } CGAL_TRACE_STREAM << "**Timer** System construction: " << timer.time() << std::endl; timer.reset(); @@ -183,7 +187,7 @@ public: // update id = 0; for(typename std::set::iterator it = interior_vertices.begin(); it != interior_vertices.end(); ++it, ++id) { - (*it)->point() = Point_3(X[id], Y[id], Z[id]); + put(ppmap, *it, Point_3(X[id], Y[id], Z[id])); } return true; } @@ -221,24 +225,24 @@ The larger @a continuity gets, the more fixed vertices are required. */ template bool fair(Polyhedron& polyhedron, - InputIterator vertex_begin, - InputIterator vertex_end, - WeightCalculator weight_calculator, - Fairing_continuity continuity = FAIRING_C_1 - ) + InputIterator vertex_begin, + InputIterator vertex_end, + WeightCalculator weight_calculator, + Fairing_continuity continuity = FAIRING_C_1 + ) { - internal::Fair_Polyhedron_3 fair_functor(weight_calculator); - return fair_functor.fair(polyhedron, vertex_begin, vertex_end, continuity); + internal::Fair_Polyhedron_3 fair_functor(polyhedron, weight_calculator); + return fair_functor.fair(vertex_begin, vertex_end, continuity); } //use default SparseLinearSolver template bool fair(Polyhedron& poly, - InputIterator vb, - InputIterator ve, - WeightCalculator weight_calculator, - Fairing_continuity continuity = FAIRING_C_1 - ) + InputIterator vb, + InputIterator ve, + WeightCalculator weight_calculator, + Fairing_continuity continuity = FAIRING_C_1 + ) { typedef internal::Fair_default_sparse_linear_solver::Solver Sparse_linear_solver; return fair @@ -248,9 +252,9 @@ bool fair(Polyhedron& poly, //use default WeightCalculator template bool fair(Polyhedron& poly, - InputIterator vb, - InputIterator ve, - Fairing_continuity continuity = FAIRING_C_1 + InputIterator vb, + InputIterator ve, + Fairing_continuity continuity = FAIRING_C_1 ) { typedef internal::Cotangent_weight_with_voronoi_area_fairing Weight_calculator; @@ -261,10 +265,10 @@ bool fair(Polyhedron& poly, //use default SparseLinearSolver and WeightCalculator template bool fair(Polyhedron& poly, - InputIterator vb, - InputIterator ve, - Fairing_continuity continuity = FAIRING_C_1 - ) + InputIterator vb, + InputIterator ve, + Fairing_continuity continuity = FAIRING_C_1 + ) { typedef internal::Fair_default_sparse_linear_solver::Solver Sparse_linear_solver; return fair(poly, vb, ve, continuity); diff --git a/Polygon_mesh_processing/include/CGAL/internal/Hole_filling/Refine_Polyhedron_3.h b/Polygon_mesh_processing/include/CGAL/internal/Hole_filling/Refine_Polyhedron_3.h index 06a4927b939..1126302b3ec 100644 --- a/Polygon_mesh_processing/include/CGAL/internal/Hole_filling/Refine_Polyhedron_3.h +++ b/Polygon_mesh_processing/include/CGAL/internal/Hole_filling/Refine_Polyhedron_3.h @@ -26,8 +26,10 @@ class Refine_Polyhedron_3 { typedef Halfedge_around_target_circulator Halfedge_around_vertex_circulator; private: - - bool flippable(Polyhedron& poly, Halfedge_handle h) { + Polyhedron& poly; + Point_property_map ppmap; + + bool flippable(Halfedge_handle h) { // this check is added so that edge flip does not break manifoldness // it might happen when there is an edge where flip_edge(h) will be placed (i.e. two edges collide after flip) Vertex_handle v_tip_0 = target(next(h,poly),poly); @@ -38,24 +40,24 @@ private: } while(++v_cir != v_end); // also eliminate collinear triangle generation - if( CGAL::collinear(v_tip_0->point(), v_tip_1->point(), target(h, poly)->point()) || - CGAL::collinear(v_tip_0->point(), v_tip_1->point(), target(opposite(h, poly),poly)->point()) ) { + if( CGAL::collinear(ppmap[v_tip_0], ppmap[v_tip_1], ppmap[target(h, poly)]) || + CGAL::collinear(ppmap[v_tip_0], ppmap[v_tip_1], ppmap[target(opposite(h, poly),poly)]) ) { return false; } return true; } - bool relax(Polyhedron& poly, Halfedge_handle h) + bool relax(Halfedge_handle h) { - const Point_3& p = target(h,poly)->point(); - const Point_3& q = target(opposite(h,poly),poly)->point(); - const Point_3& r = target(next(h,poly),poly)->point(); - const Point_3& s = target(next(opposite(h,poly),poly),poly)->point(); + const Point_3& p = ppmap[target(h,poly)]; + const Point_3& q = ppmap[target(opposite(h,poly),poly)]; + const Point_3& r = ppmap[target(next(h,poly),poly)]; + const Point_3& s = ppmap[target(next(opposite(h,poly),poly),poly)]; if( (CGAL::ON_UNBOUNDED_SIDE != CGAL::side_of_bounded_sphere(p,q,r,s)) || (CGAL::ON_UNBOUNDED_SIDE != CGAL::side_of_bounded_sphere(p,q,s,r)) ) { - if(flippable(poly,h)) { + if(flippable(h)) { poly.flip_edge(h); return true; } @@ -64,8 +66,7 @@ private: } template - bool subdivide(Polyhedron& poly, - std::vector& facets, + bool subdivide(std::vector& facets, const std::set& border_edges, std::map& scale_attribute, VertexOutputIterator& vertex_out, @@ -80,11 +81,11 @@ private: Vertex_handle vi = target(halfedge(facets[i],poly),poly); Vertex_handle vj = target(next(halfedge(facets[i],poly),poly),poly); Vertex_handle vk = target(prev(halfedge(facets[i],poly),poly),poly); - Point_3 c = CGAL::centroid(vi->point(), vj->point(), vk->point()); + Point_3 c = CGAL::centroid(ppmap[vi], ppmap[vj], ppmap[vk]); double sac = (scale_attribute[vi] + scale_attribute[vj] + scale_attribute[vk])/3.0; - double dist_c_vi = std::sqrt(CGAL::squared_distance(c,vi->point())); - double dist_c_vj = std::sqrt(CGAL::squared_distance(c,vj->point())); - double dist_c_vk = std::sqrt(CGAL::squared_distance(c,vk->point())); + double dist_c_vi = std::sqrt(CGAL::squared_distance(c,ppmap[vi])); + double dist_c_vj = std::sqrt(CGAL::squared_distance(c,ppmap[vj])); + double dist_c_vk = std::sqrt(CGAL::squared_distance(c,ppmap[vk])); if((alpha * dist_c_vi > sac) && (alpha * dist_c_vj > sac) && (alpha * dist_c_vk > sac) && @@ -92,7 +93,7 @@ private: (alpha * dist_c_vj > scale_attribute[vj]) && (alpha * dist_c_vk > scale_attribute[vk])){ Halfedge_handle h = poly.create_center_vertex(halfedge(facets[i],poly)); - target(h,poly)->point() = c; + put(ppmap, target(h,poly)->point(), c); scale_attribute[target(h,poly)] = sac; *vertex_out++ = target(h,poly); @@ -107,21 +108,20 @@ private: Halfedge_handle e_jk = prev(opposite(next(h,poly),poly),poly); if(border_edges.find(e_ij) == border_edges.end()){ - relax(poly, e_ij); + relax(e_ij); } if(border_edges.find(e_ik) == border_edges.end()){ - relax(poly, e_ik); + relax(e_ik); } if(border_edges.find(e_jk) == border_edges.end()){ - relax(poly, e_jk); + relax(e_jk); } } } return facets.size() != facet_size; } - bool relax(Polyhedron& poly, - const std::vector& facets, + bool relax(const std::vector& facets, const std::set& border_edges) { int flips = 0; @@ -147,7 +147,7 @@ private: CGAL_TRACE_STREAM << "Test " << interior_edges.size() << " edges " << std::endl; //do not just use std::set (included_map) for iteration, the order effects the output (we like to make it deterministic) for(typename std::list::iterator it = interior_edges.begin(); it != interior_edges.end();++it) { - if(relax(poly,*it)) { + if(relax(*it)) { ++flips; } } @@ -156,12 +156,11 @@ private: return flips > 0; } - double average_length(Polyhedron& poly, - Vertex_handle vh, + double average_length(Vertex_handle vh, const std::set& interior_map, bool accept_internal_facets) { - const Point_3& vp = vh->point(); + const Point_3& vp = ppmap[vh]; Halfedge_around_target_circulator circ(halfedge(vh,poly),poly), done(circ); int deg = 0; double sum = 0; @@ -173,7 +172,7 @@ private: { continue; } // which means current edge is an interior edge and should not be included in scale attribute calculation } - const Point_3& vq = target(opposite(*circ,poly),poly)->point(); + const Point_3& vq = ppmap[target(opposite(*circ,poly),poly)]; sum += std::sqrt(CGAL::squared_distance(vp, vq)); ++deg; } while(++circ != done); @@ -182,8 +181,7 @@ private: return sum/deg; } - void calculate_scale_attribute(Polyhedron& poly, - const std::vector& facets, + void calculate_scale_attribute(const std::vector& facets, const std::set& interior_map, std::map& scale_attribute, bool accept_internal_facets) @@ -195,13 +193,12 @@ private: std::pair::iterator, bool> v_insert = scale_attribute.insert(std::make_pair(v, 0)); if(!v_insert.second) { continue; } // already calculated - v_insert.first->second = average_length(poly, v, interior_map, accept_internal_facets); + v_insert.first->second = average_length(v, interior_map, accept_internal_facets); } while(++circ != done); } } - bool contain_internal_facets(Polyhedron& poly, - const std::vector& facets, + bool contain_internal_facets(const std::vector& facets, const std::set& interior_map) const { for(typename std::vector::const_iterator f_it = facets.begin(); f_it != facets.end(); ++f_it) { @@ -226,9 +223,12 @@ private: } public: + Refine_Polyhedron_3(Polyhedron& poly) + : poly(poly), ppmap(get(CGAL::vertex_point, poly)) + {} + template - void refine(Polyhedron& poly, - InputIterator facet_begin, + void refine(InputIterator facet_begin, InputIterator facet_end, FacetOutputIterator& facet_out, VertexOutputIterator& vertex_out, @@ -249,18 +249,18 @@ public: } // check whether there is any need to accept internal facets - bool accept_internal_facets = contain_internal_facets(poly, facets, interior_map); + bool accept_internal_facets = contain_internal_facets(facets, interior_map); std::map scale_attribute; - calculate_scale_attribute(poly, facets, interior_map, scale_attribute, accept_internal_facets); + calculate_scale_attribute(facets, interior_map, scale_attribute, accept_internal_facets); CGAL::Timer total_timer; total_timer.start(); for(int i = 0; i < 10; ++i) { CGAL::Timer timer; timer.start(); - bool is_subdivided = subdivide(poly, facets, border_edges, scale_attribute, vertex_out, facet_out, alpha); + bool is_subdivided = subdivide(facets, border_edges, scale_attribute, vertex_out, facet_out, alpha); CGAL_TRACE_STREAM << "**Timer** subdivide() :" << timer.time() << std::endl; timer.reset(); if(!is_subdivided) { break; } - bool is_relaxed = relax(poly, facets, border_edges); + bool is_relaxed = relax(facets, border_edges); CGAL_TRACE_STREAM << "**Timer** relax() :" << timer.time() << std::endl; if(!is_relaxed) { break; } } @@ -300,16 +300,16 @@ template< class VertexOutputIterator > std::pair -refine(Polyhedron& poly, +refine(Polyhedron& poly, InputIterator facet_begin, InputIterator facet_end, FacetOutputIterator facet_out, VertexOutputIterator vertex_out, double density_control_factor = std::sqrt(2.0)) { - internal::Refine_Polyhedron_3 refine_functor; + internal::Refine_Polyhedron_3 refine_functor(poly); refine_functor.refine - (poly, facet_begin, facet_end, facet_out, vertex_out, density_control_factor); + (facet_begin, facet_end, facet_out, vertex_out, density_control_factor); return std::make_pair(facet_out, vertex_out); }