From e2b33c4948f6488333bf552e75c0376d1b5b43e1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Sun, 29 Oct 2017 21:31:59 +0100 Subject: [PATCH] Progress on periodic optimizers - Replaced a lot of v->point() with tr.point(v) & similar - Fixed taking references to temporary: we can't do "F(construct_point_3_object()(tr.point(c, i)))" for periodic triangulations because tr.point(c,i) is _not_ a reference - Added some sanity checks to debug Lloyd optimization - Fixed some indentation --- .../IO/facets_in_complex_3_to_triangle_mesh.h | 6 +- Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h | 138 +++++---- .../Mesh_3/Cell_criteria_visitor_with_balls.h | 2 +- Mesh_3/include/CGAL/Mesh_3/Dump_c3t3.h | 19 +- .../Facet_criteria_visitor_with_balls.h | 2 +- Mesh_3/include/CGAL/Mesh_3/Lloyd_move.h | 74 +++-- .../Mesh_complex_3_in_triangulation_3_base.h | 14 +- .../CGAL/Mesh_3/Mesh_global_optimizer.h | 66 ++--- .../include/CGAL/Mesh_3/Mesh_sizing_field.h | 64 +++-- Mesh_3/include/CGAL/Mesh_3/Odt_move.h | 54 ++-- Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h | 52 ++-- .../CGAL/Mesh_3/Refine_facets_manifold_base.h | 179 ++++++------ Mesh_3/include/CGAL/Mesh_3/Sliver_perturber.h | 17 +- Mesh_3/include/CGAL/Mesh_3/Slivers_exuder.h | 62 +++-- .../CGAL/Mesh_3/Triangulation_helpers.h | 56 ++-- .../CGAL/Mesh_3/Triangulation_sizing_field.h | 60 ++-- ...tialize_triangulation_from_labeled_image.h | 18 +- .../CGAL/Mesh_3/mesh_standard_cell_criteria.h | 56 ++-- .../Mesh_3/mesh_standard_facet_criteria.h | 36 ++- .../include/CGAL/Mesh_3/vertex_perturbation.h | 262 ++++++++++-------- .../CGAL/Mesh_complex_3_in_triangulation_3.h | 23 +- .../CGAL/Polyhedral_complex_mesh_domain_3.h | 15 +- Mesh_3/include/CGAL/refine_mesh_3.h | 4 +- .../Periodic_3_mesh_3/mesh_implicit_shape.cpp | 8 +- .../CGAL/Periodic_3_mesh_triangulation_3.h | 16 ++ 25 files changed, 707 insertions(+), 596 deletions(-) diff --git a/Mesh_3/include/CGAL/IO/facets_in_complex_3_to_triangle_mesh.h b/Mesh_3/include/CGAL/IO/facets_in_complex_3_to_triangle_mesh.h index aec58ce46ff..cef36a8d688 100644 --- a/Mesh_3/include/CGAL/IO/facets_in_complex_3_to_triangle_mesh.h +++ b/Mesh_3/include/CGAL/IO/facets_in_complex_3_to_triangle_mesh.h @@ -44,14 +44,16 @@ namespace CGAL { //! @param c3t3 an instance of `C3T3`. //! @param graph an instance of `TriangleMesh`. template -void facets_in_complex_3_to_triangle_mesh(const C3T3& c3t3, TriangleMesh& graph) //complexity nlogn(number of facets on surface) +void facets_in_complex_3_to_triangle_mesh(const C3T3& c3t3, TriangleMesh& graph) { + //complexity nlogn(number of facets on surface) + + // needs a rewrite (see medit periodic) typedef typename boost::property_map::type VertexPointMap; typedef typename boost::property_traits::value_type Point_3; typedef typename C3T3::Triangulation Tr; typedef typename Tr::Vertex_handle Vertex_handle; - typedef typename C3T3::Triangulation Triangulation; typedef typename Triangulation::Vertex_handle Vertex_handle; typedef typename Triangulation::Weighted_point Weighted_point; diff --git a/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h b/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h index 3e83e7117cd..fd4ea127807 100644 --- a/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h +++ b/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h @@ -641,9 +641,6 @@ class C3T3_helpers typedef typename Gt::Plane_3 Plane_3; typedef typename Gt::Tetrahedron_3 Tetrahedron; - typedef typename Gt::Construct_point_3 Construct_point_3; - typedef typename Gt::Construct_weighted_point_3 Construct_weighted_point_3; - typedef typename Tr::Vertex_handle Vertex_handle; typedef typename Tr::Facet Facet; typedef typename Tr::Cell Cell; @@ -717,8 +714,6 @@ public: , c3t3_(c3t3) , tr_(c3t3.triangulation()) , domain_(domain) - , wp2p_(tr_.geom_traits().construct_point_3_object()) - , p2wp_(tr_.geom_traits().construct_weighted_point_3_object()) { } /** @@ -1248,13 +1243,13 @@ private: // Update facet surface center Bare_point surface_center = CGAL::cpp0x::get<0>(intersection); - facet.first->set_facet_surface_center(facet.second,surface_center); - facet_m.first->set_facet_surface_center(facet_m.second,surface_center); + facet.first->set_facet_surface_center(facet.second, surface_center); + facet_m.first->set_facet_surface_center(facet_m.second, surface_center); } else { - facet.first->set_facet_surface_center(facet.second,Bare_point()); - facet_m.first->set_facet_surface_center(facet_m.second,Bare_point()); + facet.first->set_facet_surface_center(facet.second, Bare_point()); + facet_m.first->set_facet_surface_center(facet_m.second, Bare_point()); } } @@ -2408,9 +2403,6 @@ private: C3T3& c3t3_; Tr& tr_; const MeshDomain& domain_; - // @todo purge the two members below? - Construct_point_3 wp2p_; - Construct_weighted_point_3 p2wp_; }; // class C3T3_helpers @@ -2425,9 +2417,10 @@ update_mesh(const Vertex_handle& old_vertex, OutputIterator modified_vertices, bool *could_lock_zone) { - // std::cerr << "\nupdate_mesh[v1](" << new_position << ",\n" - // << " " << (void*)(&*old_vertex) << "=" << old_vertex->point() - // << ")\n"; + // std::cerr << "\nupdate_mesh[v1](" << (void*)(&*old_vertex) + // << "=" << tr_.point(old_vertex) << ",\n" + // << " " << move << ",\n" + // << " " << new_position << ")\n"; if (could_lock_zone) *could_lock_zone = true; @@ -2459,17 +2452,19 @@ update_mesh_no_topo_change(const Vertex_handle& old_vertex, const Cell_vector& conflict_cells ) { // std::cerr << "update_mesh_no_topo_change(" - // << (void*)(&*old_vertex) << " = " << old_vertex->point() << ",\n" + // << (void*)(&*old_vertex) << " = " << tr_.point(old_vertex) << ",\n" // << " " << move << ",\n" // << " " << new_position << ")" << std::endl; + typename Gt::Construct_opposite_vector_3 cov = tr_.geom_traits().construct_opposite_vector_3_object(); + //backup metadata std::set cells_backup; fill_cells_backup(conflict_cells, cells_backup); // Get old values criterion.before_move(c3t3_cells(conflict_cells)); - Weighted_point old_position = old_vertex->point(); + Weighted_point old_position = tr_.point(old_vertex); // intentional copy // Move point reset_circumcenter_cache(conflict_cells); @@ -2496,10 +2491,6 @@ update_mesh_no_topo_change(const Vertex_handle& old_vertex, //sliver caches have been updated by valid_move reset_sliver_cache(conflict_cells); - - typename Gt::Construct_opposite_vector_3 cov = - tr_.geom_traits().construct_opposite_vector_3_object(); - move_point_no_topo_change(old_vertex, cov(move), old_position); //restore meta-data (cells should have same connectivity as before move) @@ -2507,7 +2498,7 @@ update_mesh_no_topo_change(const Vertex_handle& old_vertex, CGAL_assertion(conflict_cells.size() >= cells_backup.size()); restore_from_cells_backup(conflict_cells, cells_backup); - return std::make_pair(false,old_vertex); + return std::make_pair(false, old_vertex); } } @@ -2522,8 +2513,8 @@ update_mesh_topo_change(const Vertex_handle& old_vertex, OutputIterator modified_vertices, bool *could_lock_zone) { - // std::cerr << "update_mesh_topo_change(" - // << (void*)(&*old_vertex) << "=" << old_vertex->point() + // std::cerr << "update_mesh_topo_change(" << (void*)(&*old_vertex) + // << "=" << tr_.point(old_vertex) // << " " << move << ",\n" // << " " << new_position << ",\n" // << ")" << std::endl; @@ -2558,7 +2549,7 @@ update_mesh_topo_change(const Vertex_handle& old_vertex, CGAL_assertion(conflict_cells.size() == cells_backup.size()); criterion.before_move(c3t3_cells(conflict_cells)); - Weighted_point old_position = old_vertex->point(); + Weighted_point old_position = tr_.point(old_vertex); // intentional copy // Keep old boundary Vertex_set old_incident_surface_vertices; @@ -2581,7 +2572,7 @@ update_mesh_topo_change(const Vertex_handle& old_vertex, CGAL::Emptyset_iterator()); // If nothing changed, return - if ( old_position == new_vertex->point() ) + if ( old_position == tr_.point(new_vertex) ) { // std::cerr << "update_mesh_topo_change: no move!\n"; // check_c3t3(c3t3_); @@ -2639,9 +2630,9 @@ update_mesh(const Vertex_handle& old_vertex, OutputIterator modified_vertices, bool fill_vertices) { - // std::cerr << "\nupdate_mesh[v2](" << new_position << ",\n" - // << " " << (void*)(&*old_vertex) << "=" << old_vertex->point() - // << ")\n"; + // std::cerr << "\nupdate_mesh[v2](" << (void*)(&*old_vertex) + // << "=" << tr_.point(old_vertex) << ",\n" + // << " " << new_position << ")\n"; Cell_vector outdated_cells; Vertex_handle new_vertex = move_point(old_vertex, move, @@ -2670,7 +2661,9 @@ C3T3_helpers:: rebuild_restricted_delaunay(OutdatedCells& outdated_cells, Moving_vertices_set& moving_vertices) { + typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); typename Gt::Equal_3 equal = tr_.geom_traits().equal_3_object(); + typename OutdatedCells::iterator first_cell = outdated_cells.begin(); typename OutdatedCells::iterator last_cell = outdated_cells.end(); Update_c3t3 updater(domain_,c3t3_); @@ -2768,8 +2761,8 @@ rebuild_restricted_delaunay(OutdatedCells& outdated_cells, it != vertex_to_proj.end() ; ++it ) { - Bare_point initial_position = wp2p_((*it)->point()); - Bare_point new_pos = project_on_surface(*it, initial_position); + const Weighted_point& initial_position = tr_.point(*it); + Bare_point new_pos = project_on_surface(*it, cp(initial_position)); std::cout << "project: " << initial_position << " got: " << new_pos << std::endl; @@ -2779,7 +2772,7 @@ rebuild_restricted_delaunay(OutdatedCells& outdated_cells, // Update moving vertices (it becomes new_vertex) moving_vertices.erase(*it); - Vector_3 move(initial_position, new_pos); + Vector_3 move(cp(initial_position), new_pos); std::cout << "projection move: " << move << std::endl; CGAL_assertion(CGAL::abs(move.x()) < 0.8*(tr_.domain().xmax() - tr_.domain().xmin())); @@ -2807,8 +2800,8 @@ rebuild_restricted_delaunay(ForwardIterator first_cell, ForwardIterator last_cell, Moving_vertices_set& moving_vertices) { - typename Gt::Equal_3 equal = tr_.geom_traits().equal_3_object(); typename Gt::Construct_vector_3 vector = tr_.geom_traits().construct_vector_3_object(); + typename Gt::Equal_3 equal = tr_.geom_traits().equal_3_object(); Update_c3t3 updater(domain_,c3t3_); @@ -2886,16 +2879,17 @@ rebuild_restricted_delaunay(ForwardIterator first_cell, end = vertex_to_proj.end(); for ( ; it != end; ++it ) { - const Bare_point& initial_pos = wp2p_((it->first)->point()); - Bare_point new_pos = project_on_surface(it->first, initial_pos, it->second); + Vertex_handle vh = it->first; + const Weighted_point& initial_position = tr_.point(vh); + Bare_point new_pos = project_on_surface(vh, cp(initial_position), it->second); if ( ! equal(new_pos, Bare_point()) ) { //freezing needs 'erase' to be done before the vertex is actually destroyed // Update moving vertices (it becomes new_vertex) - moving_vertices.erase(it->first); + moving_vertices.erase(vh); - Vertex_handle new_vertex = update_mesh(it->first, vector(initial_pos, new_pos)); + Vertex_handle new_vertex = update_mesh(vh, vector(cp(initial_position), new_pos)); c3t3_.set_dimension(new_vertex, 2); moving_vertices.insert(new_vertex); @@ -2925,22 +2919,20 @@ move_point(const Vertex_handle& old_vertex, OutdatedCellsOutputIterator outdated_cells, DeletedCellsOutputIterator deleted_cells) const { - // std::cerr << "C3T3_helpers::move_point[v2](" - // << (void*)(&*old_vertex) << " = " << old_vertex->point() - // << " , move: " << move << ")\n"; + // std::cerr << "C3T3_helpers::move_point[v2](" << (void*)(&*old_vertex) + // << " = " << tr_.point(old_vertex) << ",\n" + // << " " << move << ")\n"; - typename Gt::Construct_translated_point_3 translate = - tr_.geom_traits().construct_translated_point_3_object(); - typename Gt::Construct_weighted_point_3 p2wp = - tr_.geom_traits().construct_weighted_point_3_object(); - typename Gt::Construct_point_3 wp2p = - tr_.geom_traits().construct_point_3_object(); + typename Gt::Construct_translated_point_3 translate = tr_.geom_traits().construct_translated_point_3_object(); + typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); + typename Gt::Construct_weighted_point_3 cwp = tr_.geom_traits().construct_weighted_point_3_object(); Cell_vector incident_cells_; incident_cells_.reserve(64); tr_.incident_cells(old_vertex, std::back_inserter(incident_cells_)); - const Weighted_point& new_position = p2wp(translate(wp2p(tr_.point(old_vertex)), move)); + const Weighted_point& position = tr_.point(old_vertex); + const Weighted_point& new_position = cwp(translate(cp(position), move)); if ( Th().no_topological_change(tr_, old_vertex, move, new_position, incident_cells_) ) { @@ -2964,18 +2956,16 @@ move_point(const Vertex_handle& old_vertex, Outdated_cell_set& outdated_cells_set, Moving_vertices_set& moving_vertices) const { - typename Gt::Construct_translated_point_3 translate = - tr_.geom_traits().construct_translated_point_3_object(); - typename Gt::Construct_weighted_point_3 p2wp = - tr_.geom_traits().construct_weighted_point_3_object(); - typename Gt::Construct_point_3 wp2p = - tr_.geom_traits().construct_point_3_object(); + typename Gt::Construct_translated_point_3 translate = tr_.geom_traits().construct_translated_point_3_object(); + typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); + typename Gt::Construct_weighted_point_3 cwp = tr_.geom_traits().construct_weighted_point_3_object(); Cell_vector incident_cells_; incident_cells_.reserve(64); tr_.incident_cells(old_vertex, std::back_inserter(incident_cells_)); - const Weighted_point& new_position = p2wp(translate(wp2p(tr_.point(old_vertex)), move)); + const Weighted_point& position = tr_.point(old_vertex); + const Weighted_point& new_position = cwp(translate(cp(position), move)); if ( Th().no_topological_change(tr_, old_vertex, move, new_position, incident_cells_) ) { @@ -3029,14 +3019,12 @@ move_point(const Vertex_handle& old_vertex, } //======= /Get incident cells ========== - typename Gt::Construct_translated_point_3 translate = - tr_.geom_traits().construct_translated_point_3_object(); - typename Gt::Construct_weighted_point_3 p2wp = - tr_.geom_traits().construct_weighted_point_3_object(); - typename Gt::Construct_point_3 wp2p = - tr_.geom_traits().construct_point_3_object(); + typename Gt::Construct_translated_point_3 translate = tr_.geom_traits().construct_translated_point_3_object(); + typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); + typename Gt::Construct_weighted_point_3 cwp = tr_.geom_traits().construct_weighted_point_3_object(); - const Weighted_point& new_position = p2wp(translate(wp2p(tr_.point(old_vertex)), move)); + const Weighted_point& position = tr_.point(old_vertex); + const Weighted_point& new_position = cwp(translate(cp(position), move)); if (!try_lock_point(new_position)) // LOCK { @@ -3198,7 +3186,7 @@ move_point_topo_change_conflict_zone_known( //o.w. deleted_cells will point to null pointer or so and crash const { - Weighted_point old_position = old_vertex->point(); + Weighted_point old_position = tr_.point(old_vertex); // intentional copy // make one set with conflict zone Cell_set conflict_zone; std::set_union(insertion_conflict_cells_begin, insertion_conflict_cells_end, @@ -3325,16 +3313,11 @@ project_on_surface_aux(const Bare_point& p, typedef typename Gt::Segment_3 Segment_3; // Build a segment directed as projection_direction, - typename Gt::Compute_squared_distance_3 sq_distance = - tr_.geom_traits().compute_squared_distance_3_object(); - typename Gt::Compute_squared_length_3 sq_length = - tr_.geom_traits().compute_squared_length_3_object(); - typename Gt::Construct_scaled_vector_3 scale = - tr_.geom_traits().construct_scaled_vector_3_object(); - typename Gt::Is_degenerate_3 is_degenerate = - tr_.geom_traits().is_degenerate_3_object(); - typename Gt::Construct_translated_point_3 translate = - tr_.geom_traits().construct_translated_point_3_object(); + typename Gt::Compute_squared_distance_3 sq_distance = tr_.geom_traits().compute_squared_distance_3_object(); + typename Gt::Compute_squared_length_3 sq_length = tr_.geom_traits().compute_squared_length_3_object(); + typename Gt::Construct_scaled_vector_3 scale = tr_.geom_traits().construct_scaled_vector_3_object(); + typename Gt::Is_degenerate_3 is_degenerate = tr_.geom_traits().is_degenerate_3_object(); + typename Gt::Construct_translated_point_3 translate = tr_.geom_traits().construct_translated_point_3_object(); typename MD::Construct_intersection construct_intersection = domain_.construct_intersection_object(); @@ -3388,6 +3371,8 @@ get_least_square_surface_plane(const Vertex_handle& v, Bare_point& reference_point, Surface_patch_index patch_index) const { + typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); + // Get incident facets Facet_vector facets; # ifdef CGAL_LINKED_WITH_TBB @@ -3416,7 +3401,8 @@ get_least_square_surface_plane(const Vertex_handle& v, const int& i = fit->second; // @fixme really ugly - const Bare_point& bp = tr_.get_closest_point(wp2p_(tr_.point(v)), + const Weighted_point& position = tr_.point(v); + const Bare_point& bp = tr_.get_closest_point(cp(position), cell->get_facet_surface_center(i)); surface_point_vector.push_back(bp); } @@ -3453,6 +3439,7 @@ project_on_surface(const Vertex_handle& v, { // return domain_.project_on_surface(p); + typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); typename Gt::Equal_3 equal = tr_.geom_traits().equal_3_object(); // Get plane @@ -3463,8 +3450,9 @@ project_on_surface(const Vertex_handle& v, return p; // Project - if ( ! equal(p, wp2p_(v->point())) ) - return project_on_surface_aux(p, wp2p_(v->point()), plane.orthogonal_vector()); + const Weighted_point& position = tr_.point(v); + if ( ! equal(p, cp(position)) ) + return project_on_surface_aux(p, cp(position), plane.orthogonal_vector()); else return project_on_surface_aux(p, reference_point, plane.orthogonal_vector()); } diff --git a/Mesh_3/include/CGAL/Mesh_3/Cell_criteria_visitor_with_balls.h b/Mesh_3/include/CGAL/Mesh_3/Cell_criteria_visitor_with_balls.h index e35ade55910..0dc052019fa 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Cell_criteria_visitor_with_balls.h +++ b/Mesh_3/include/CGAL/Mesh_3/Cell_criteria_visitor_with_balls.h @@ -66,7 +66,7 @@ public: : Base(ch) { - const Weighted_point& p = ch->vertex(0)->point(); + const Weighted_point& p = ch->vertex(0)->point(); // ? const Weighted_point& q = ch->vertex(1)->point(); const Weighted_point& r = ch->vertex(2)->point(); const Weighted_point& s = ch->vertex(3)->point(); diff --git a/Mesh_3/include/CGAL/Mesh_3/Dump_c3t3.h b/Mesh_3/include/CGAL/Mesh_3/Dump_c3t3.h index 9af9d966ff5..fb9a0d6dae9 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Dump_c3t3.h +++ b/Mesh_3/include/CGAL/Mesh_3/Dump_c3t3.h @@ -44,7 +44,8 @@ template ::is_specialized) > struct Dump_c3t3 { - void dump_c3t3(const C3t3& c3t3, std::string prefix) const { + void dump_c3t3(const C3t3& c3t3, std::string prefix) const + { std::clog<<"======dump c3t3===== to: " << prefix << std::endl; std::ofstream medit_file((prefix+".mesh").c_str()); medit_file.precision(17); @@ -64,7 +65,8 @@ struct Dump_c3t3 { }; // end struct template Dump_c3t3 template -struct Dump_c3t3 { +struct Dump_c3t3 +{ void dump_c3t3(const C3t3&, std::string) { std::cerr << "Warning " << __FILE__ << ":" << __LINE__ << "\n" << " the c3t3 object of following type:\n" @@ -104,8 +106,9 @@ struct Dump_c3t3 { }; // end struct template specialization Dump_c3t3 template -void dump_c3t3_edges(const C3t3& c3t3, std::string prefix) { - typename C3t3::Triangulation::Geom_traits::Construct_point_3 wp2p = +void dump_c3t3_edges(const C3t3& c3t3, std::string prefix) +{ + typename C3t3::Triangulation::Geom_traits::Construct_point_3 cp = c3t3.triangulation().geom_traits().construct_point_3_object(); std::ofstream file((prefix+".polylines.txt").c_str()); @@ -118,12 +121,14 @@ void dump_c3t3_edges(const C3t3& c3t3, std::string prefix) { const typename C3t3::Triangulation::Cell_handle c = edge_it->first; const int i = edge_it->second; const int j = edge_it->third; - file << "2 " << wp2p(c->vertex(i)->point()) - << " " << wp2p(c->vertex(j)->point()) << "\n"; + const typename C3t3::Weighted_point& ei = c3t3.triangulation().point(c, i); + const typename C3t3::Weighted_point& ej = c3t3.triangulation().point(c, j); + file << "2 " << cp(ei) << " " << cp(ej) << "\n"; } } template -void dump_c3t3(const C3t3& c3t3, std::string prefix) { +void dump_c3t3(const C3t3& c3t3, std::string prefix) +{ if(!prefix.empty()) { Dump_c3t3 dump; dump.dump_c3t3(c3t3, prefix); diff --git a/Mesh_3/include/CGAL/Mesh_3/Facet_criteria_visitor_with_balls.h b/Mesh_3/include/CGAL/Mesh_3/Facet_criteria_visitor_with_balls.h index 9eee13517b2..321c46485f6 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Facet_criteria_visitor_with_balls.h +++ b/Mesh_3/include/CGAL/Mesh_3/Facet_criteria_visitor_with_balls.h @@ -78,7 +78,7 @@ public: Squared_radius_orthogonal_sphere sq_radius_ortho_sphere = Gt().compute_squared_radius_smallest_orthogonal_sphere_3_object(); - Weighted_point p1 = fh.first->vertex ((fh.second+1)&3)->point(); + Weighted_point p1 = fh.first->vertex ((fh.second+1)&3)->point(); // ? Weighted_point p2 = fh.first->vertex ((fh.second+2)&3)->point(); Weighted_point p3 = fh.first->vertex ((fh.second+3)&3)->point(); diff --git a/Mesh_3/include/CGAL/Mesh_3/Lloyd_move.h b/Mesh_3/include/CGAL/Mesh_3/Lloyd_move.h index 8e0fb98dfb2..80802a4d649 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Lloyd_move.h +++ b/Mesh_3/include/CGAL/Mesh_3/Lloyd_move.h @@ -55,6 +55,7 @@ class Lloyd_move typedef typename Tr::Cell_handle Cell_handle; typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Weighted_point Weighted_point; typedef typename std::vector Facet_vector; typedef typename std::vector Cell_vector; @@ -79,7 +80,7 @@ public: const Sizing_field& sizing_field = Sizing_field() ) const { // std::cout << "computing move of: " << &*v -// << " pos: " << v->point().point() +// << " pos: " << c3t3.triangulation().point(v) // << " dim: " << c3t3.in_dimension(v) << std::endl; switch ( c3t3.in_dimension(v) ) @@ -260,7 +261,7 @@ private: { const Tr& tr = c3t3.triangulation(); - typename Gt::Construct_point_3 wp2p = tr.geom_traits().construct_point_3_object(); + typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object(); Facet_vector incident_facets; incident_facets.reserve(64); @@ -281,9 +282,10 @@ private: const Cell_handle& cell = fit->first; const int& i = fit->second; - // @fixme really ugly - points.push_back(tr.get_closest_point(wp2p(tr.point(v)), - cell->get_facet_surface_center(i))); + // @fixme 'get_closest_point' is really ugly + const Weighted_point& position = tr.point(v); + points.push_back(tr.get_closest_point(cp(position), + cell->get_facet_surface_center(i))); } return points; @@ -298,20 +300,19 @@ private: const C3T3& c3t3, const Sizing_field& sizing_field) const { - typename Gt::Construct_point_3 cp = - c3t3.triangulation().geom_traits().construct_point_3_object(); - typename Gt::Construct_vector_3 vector = - c3t3.triangulation().geom_traits().construct_vector_3_object(); + typename Gt::Construct_point_3 cp = c3t3.triangulation().geom_traits().construct_point_3_object(); + typename Gt::Construct_vector_3 vector = c3t3.triangulation().geom_traits().construct_vector_3_object(); - const Bare_point& p = cp(v->point()); + const Weighted_point position = c3t3.triangulation().point(v); + const Bare_point& p = cp(position); - FT da = density_1d(a,v,sizing_field); - FT db = density_1d(b,v,sizing_field); + FT da = density_1d(a, v, sizing_field); + FT db = density_1d(b, v, sizing_field); CGAL_assertion( !is_zero(da + db) ); return ( (vector(p,a)*da + vector(p,b)*db) / (da + db) ); } - + /** * Return the move from \c v to the centroid of triangle [a,b,c]. */ @@ -322,12 +323,11 @@ private: const C3T3& c3t3, const Sizing_field& sizing_field) const { - typename Gt::Construct_vector_3 vector = - c3t3.triangulation().geom_traits().construct_vector_3_object(); - typename Gt::Construct_point_3 cp = - c3t3.triangulation().geom_traits().construct_point_3_object(); + typename Gt::Construct_point_3 cp = c3t3.triangulation().geom_traits().construct_point_3_object(); + typename Gt::Construct_vector_3 vector = c3t3.triangulation().geom_traits().construct_vector_3_object(); - const Bare_point& p = cp(v->point()); + const Weighted_point& position = c3t3.triangulation().point(v); + const Bare_point& p = cp(position); FT da = density_2d(a,v,sizing_field); FT db = density_2d(b,v,sizing_field); @@ -360,7 +360,7 @@ private: // Project all points to the plane std::transform(first, last, first, Project_on_plane(plane, c3t3)); - CGAL_assertion(std::distance(first,last) > 3); + CGAL_assertion(std::distance(first, last) > 3); // Get 2D-3D transformations Aff_transformation_3 to_3d = compute_to_3d_transform(plane, *first, c3t3); @@ -403,17 +403,14 @@ private: { CGAL_precondition(std::distance(first,last) >= 3); - typename Gt::Construct_vector_3 vector = - c3t3.triangulation().geom_traits().construct_vector_3_object(); - typename Gt::Construct_centroid_3 centroid = - c3t3.triangulation().geom_traits().construct_centroid_3_object(); - typename Gt::Compute_area_3 area = - c3t3.triangulation().geom_traits().compute_area_3_object(); - typename Gt::Construct_point_3 cp = - c3t3.triangulation().geom_traits().construct_point_3_object(); + typename Gt::Compute_area_3 area = c3t3.triangulation().geom_traits().compute_area_3_object(); + typename Gt::Construct_centroid_3 centroid = c3t3.triangulation().geom_traits().construct_centroid_3_object(); + typename Gt::Construct_point_3 cp = c3t3.triangulation().geom_traits().construct_point_3_object(); + typename Gt::Construct_vector_3 vector = c3t3.triangulation().geom_traits().construct_vector_3_object(); // Vertex current position - const Bare_point& vertex_position = cp(v->point()); + const Weighted_point& vertex_weighted_position = c3t3.triangulation().point(v); + const Bare_point& vertex_position = cp(vertex_weighted_position); // Use as reference point to triangulate const Bare_point& a = *first++; @@ -457,10 +454,8 @@ private: const Bare_point& reference_point, const C3T3& c3t3) const { - typename Gt::Construct_base_vector_3 base = - c3t3.triangulation().geom_traits().construct_base_vector_3_object(); - typename Gt::Construct_orthogonal_vector_3 orthogonal_vector = - c3t3.triangulation().geom_traits().construct_orthogonal_vector_3_object(); + typename Gt::Construct_base_vector_3 base = c3t3.triangulation().geom_traits().construct_base_vector_3_object(); + typename Gt::Construct_orthogonal_vector_3 orthogonal_vector = c3t3.triangulation().geom_traits().construct_orthogonal_vector_3_object(); Vector_3 u = base(plane, 1); u = u / CGAL::sqrt(u*u); @@ -553,19 +548,20 @@ private: CGAL_assertion(c3t3.is_in_complex(current_cell)); // a & b are fixed points - const Bare_point& a = cp(v->point()); + const Weighted_point& wa = tr.point(v); + const Bare_point& a = cp(wa); Bare_point b = tr.dual(current_cell); - const Bare_point& a_b = cp(tr.point(current_cell, current_cell->index(v))); - Vector_3 ba = Vector_3(a_b, b); + const Weighted_point& a_b = tr.point(current_cell, current_cell->index(v)); + Vector_3 ba = Vector_3(cp(a_b), b); ++current_cell; CGAL_assertion(c3t3.is_in_complex(current_cell)); CGAL_assertion(current_cell != done); // c & d are moving points Bare_point c = tr.dual(current_cell); - const Bare_point& a_c = cp(tr.point(current_cell, current_cell->index(v))); - Vector_3 ca = Vector_3(a_c, c); + const Weighted_point& a_c = tr.point(current_cell, current_cell->index(v)); + Vector_3 ca = Vector_3(cp(a_c), c); ++current_cell; CGAL_assertion(current_cell != done); @@ -573,8 +569,8 @@ private: { CGAL_assertion(c3t3.is_in_complex(current_cell)); Bare_point d = tr.dual(current_cell); - const Bare_point& a_d = cp(tr.point(current_cell, current_cell->index(v))); - Vector_3 da = Vector_3(a_d, d); + const Weighted_point& a_d = tr.point(current_cell, current_cell->index(v)); + Vector_3 da = Vector_3(cp(a_d), d); Tetrahedron_3 tet = tetrahedron(a, translate(a, ba), translate(a, ca), translate(a, da)); Bare_point tet_centroid = centroid(tet); diff --git a/Mesh_3/include/CGAL/Mesh_3/Mesh_complex_3_in_triangulation_3_base.h b/Mesh_3/include/CGAL/Mesh_3/Mesh_complex_3_in_triangulation_3_base.h index 960c15bbc63..b45df2efc97 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Mesh_complex_3_in_triangulation_3_base.h +++ b/Mesh_3/include/CGAL/Mesh_3/Mesh_complex_3_in_triangulation_3_base.h @@ -291,7 +291,7 @@ public: default : #ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS std::cerr << "singular edge...\n"; - std::cerr << v->point() << std::endl; + std::cerr << tr_.point(v) << std::endl; #endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS return SINGULAR; } @@ -302,13 +302,13 @@ public: if(nb_components > 1) { #ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS std::cerr << "singular vertex: nb_components=" << nb_components << std::endl; - std::cerr << v->point() << std::endl; + std::cerr << tr_.point(v) << std::endl; #endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS return SINGULAR; } else { // REGULAR OR BOUNDARY #ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS - std::cerr << "regular or boundary: " << v->point() << std::endl; + std::cerr << "regular or boundary: " << tr_.point(v) << std::endl; #endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS if (number_of_boundary_incident_edges != 0) return BOUNDARY; @@ -933,7 +933,7 @@ Mesh_complex_3_in_triangulation_3_base::add_to_complex( if (j != i) { #ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS if(cell->vertex(j)->is_c2t3_cache_valid()) - std::cerr << "(" << cell->vertex(j)->point() << ")->invalidate_c2t3_cache()\n"; + std::cerr << "(" << tr_.point(cell, j) << ")->invalidate_c2t3_cache()\n"; #endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS cell->vertex(j)->invalidate_c2t3_cache(); } @@ -984,7 +984,7 @@ Mesh_complex_3_in_triangulation_3_base::remove_from_complex(const Facet& if (j != facet.second) { #ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS if(cell->vertex(j)->is_c2t3_cache_valid()) - std::cerr << "(" << cell->vertex(j)->point() << ")->invalidate_c2t3_cache()\n"; + std::cerr << "(" << tr_.point(cell, j) << ")->invalidate_c2t3_cache()\n"; #endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS cell->vertex(j)->invalidate_c2t3_cache(); } @@ -1008,12 +1008,12 @@ bbox() const } typename Tr::Finite_vertices_iterator vit = tr_.finite_vertices_begin(); - Bbox_3 result = (vit++)->point().bbox(); + Bbox_3 result = tr_.point(vit++).bbox(); for(typename Tr::Finite_vertices_iterator end = tr_.finite_vertices_end(); vit != end ; ++vit) { - result = result + vit->point().bbox(); + result = result + tr_.point(vit).bbox(); } return result; diff --git a/Mesh_3/include/CGAL/Mesh_3/Mesh_global_optimizer.h b/Mesh_3/include/CGAL/Mesh_3/Mesh_global_optimizer.h index 0d5003804f9..29d0b386176 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Mesh_global_optimizer.h +++ b/Mesh_3/include/CGAL/Mesh_3/Mesh_global_optimizer.h @@ -1,4 +1,4 @@ -// Copyright (c) 2009 INRIA Sophia-Antipolis (France). +// Copyright (c) 2009 INRIA Sophia-Antipolis (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). @@ -398,12 +398,9 @@ private: // operator() void operator()(const Vertex_handle& oldv) const { - typename Gt::Construct_point_3 wp2p = - m_gt.construct_point_3_object(); - typename Gt::Construct_weighted_point_3 p2wp = - m_gt.construct_weighted_point_3_object(); - typename Gt::Construct_translated_point_3 translate = - m_gt.construct_translated_point_3_object(); + typename Gt::Construct_point_3 cp = m_gt.construct_point_3_object(); + typename Gt::Construct_weighted_point_3 cwp = m_gt.construct_weighted_point_3_object(); + typename Gt::Construct_translated_point_3 translate = m_gt.construct_translated_point_3_object(); Vector_3 move = m_mgo.compute_move(oldv); if ( CGAL::NULL_VECTOR != move ) @@ -412,7 +409,8 @@ private: if(MGO::Sizing_field::is_vertex_update_needed) { - Bare_point new_position = translate(wp2p(oldv->point()), move); + const Weighted_point& position = m_gt.tr_.point(oldv); + Bare_point new_position = translate(cp(position), move); size = m_sizing_field(new_position, oldv); } @@ -464,13 +462,13 @@ private: // operator() void operator()(Vertex& v) const { - typename Gt::Construct_point_3 wp2p = - m_gt.construct_point_3_object(); + typename Gt::Construct_point_3 cp = m_gt.construct_point_3_object(); Vertex_handle vh = Tr::Triangulation_data_structure::Vertex_range::s_iterator_to(v); + const Weighted_point& position = m_mgo.tr_.point(vh); m_local_lists.local().push_back( - std::make_pair(wp2p(v.point()), m_mgo.average_circumradius_length(vh))); + std::make_pair(cp(position), m_mgo.average_circumradius_length(vh))); } }; @@ -508,7 +506,7 @@ private: void operator()( const tbb::blocked_range& r ) const { typename Gt::Construct_translated_point_3 translate = - m_gt.construct_translated_point_3_object(); + m_mgo.tr_.geom_traits().construct_translated_point_3_object(); for( size_t i = r.begin() ; i != r.end() ; ++i) { @@ -816,10 +814,8 @@ compute_moves(Moving_vertices_set& moving_vertices) else #endif // CGAL_LINKED_WITH_TBB { - typename Gt::Construct_point_3 wp2p = - tr_.geom_traits().construct_point_3_object(); - typename Gt::Construct_translated_point_3 translate = - tr_.geom_traits().construct_translated_point_3_object(); + typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); + typename Gt::Construct_translated_point_3 translate = tr_.geom_traits().construct_translated_point_3_object(); // Get move for each moving vertex typename Moving_vertices_set::iterator vit = moving_vertices.begin(); @@ -834,7 +830,8 @@ compute_moves(Moving_vertices_set& moving_vertices) FT size = 0.; if(Sizing_field::is_vertex_update_needed) { - Bare_point new_position = translate(wp2p(tr_.point(oldv)), move); + const Weighted_point& position = tr_.point(oldv); + Bare_point new_position = translate(cp(position), move); // std::cout << "new position: " << new_position << std::endl; size = sizing_field_(new_position, oldv); } @@ -870,14 +867,10 @@ typename Mesh_global_optimizer::Vector_3 Mesh_global_optimizer:: compute_move(const Vertex_handle& v) { - typename Gt::Compute_squared_length_3 sq_length = - tr_.geom_traits().compute_squared_length_3_object(); - typename Gt::Construct_vector_3 vector = - tr_.geom_traits().construct_vector_3_object(); - typename Gt::Construct_translated_point_3 translate = - tr_.geom_traits().construct_translated_point_3_object(); - typename Gt::Construct_point_3 wp2p = - tr_.geom_traits().construct_point_3_object(); + typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); + typename Gt::Compute_squared_length_3 sq_length = tr_.geom_traits().compute_squared_length_3_object(); + typename Gt::Construct_translated_point_3 translate = tr_.geom_traits().construct_translated_point_3_object(); + typename Gt::Construct_vector_3 vector = tr_.geom_traits().construct_vector_3_object(); Cell_vector incident_cells; incident_cells.reserve(64); @@ -904,16 +897,17 @@ compute_move(const Vertex_handle& v) // Project surface vertex if ( c3t3_.in_dimension(v) == 2 ) { - Bare_point new_position = translate(wp2p(tr_.point(v)), move); + const Weighted_point& position = tr_.point(v); + Bare_point new_position = translate(cp(position), move); // std::cout << "new pos before projection: " << new_position << std::endl; Bare_point projected_new_position = helper_.project_on_surface(v, new_position); // std::cout << "projected: " << projected_new_position << std::endl; - move = vector(wp2p(tr_.point(v)), projected_new_position); + move = vector(cp(position), projected_new_position); } FT local_move_sq_ratio = sq_length(move) / local_sq_size; -// std::cout << "optimizing point: " << tr_.point(v) << std::endl; +// std::cout << "optimizing point: " << position << std::endl; // std::cout << "moving with: " << move << std::endl; CGAL_assertion(CGAL::abs(move.x()) < 0.8*(tr_.domain().xmax() - tr_.domain().xmin())); @@ -1074,16 +1068,15 @@ fill_sizing_field() else #endif //CGAL_LINKED_WITH_TBB { - typename Gt::Construct_point_3 wp2p = - tr_.geom_traits().construct_point_3_object(); + typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); // Fill map with local size for(typename Tr::Finite_vertices_iterator vit = tr_.finite_vertices_begin(); vit != tr_.finite_vertices_end(); ++vit) { - value_map.insert(std::make_pair(wp2p(tr_.point(vit)), - average_circumradius_length(vit))); + const Weighted_point& position = tr_.point(vit); + value_map.insert(std::make_pair(cp(position), average_circumradius_length(vit))); } } @@ -1209,14 +1202,13 @@ typename Mesh_global_optimizer::FT Mesh_global_optimizer:: sq_circumradius_length(const Cell_handle& cell, const Vertex_handle& v) const { - typename Gt::Compute_squared_distance_3 sq_distance = - tr_.geom_traits().compute_squared_distance_3_object(); - typename Gt::Construct_point_3 wp2p = - tr_.geom_traits().construct_point_3_object(); + typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); + typename Gt::Compute_squared_distance_3 sq_distance = tr_.geom_traits().compute_squared_distance_3_object(); const Bare_point circumcenter = tr_.dual(cell); + const Weighted_point& position = tr_.point(cell, cell->index(v)); - return ( sq_distance(wp2p(tr_.point(cell, cell->index(v))), circumcenter) ); + return ( sq_distance(cp(position), circumcenter) ); } } // end namespace Mesh_3 diff --git a/Mesh_3/include/CGAL/Mesh_3/Mesh_sizing_field.h b/Mesh_3/include/CGAL/Mesh_3/Mesh_sizing_field.h index a6b1f24a0e8..a5743137290 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Mesh_sizing_field.h +++ b/Mesh_3/include/CGAL/Mesh_3/Mesh_sizing_field.h @@ -95,12 +95,13 @@ class Mesh_sizing_field typename Tr::Concurrency_tag> { // Types - typedef typename Tr::Geom_traits Gt; - typedef typename Tr::Bare_point Bare_point; - typedef typename Gt::FT FT; + typedef typename Tr::Geom_traits Gt; + typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Weighted_point Weighted_point; + typedef typename Gt::FT FT; - typedef typename Tr::Vertex_handle Vertex_handle; - typedef typename Tr::Cell_handle Cell_handle; + typedef typename Tr::Vertex_handle Vertex_handle; + typedef typename Tr::Cell_handle Cell_handle; public: // update vertices of mesh triangulation ? @@ -175,12 +176,14 @@ fill(const std::map& value_map) { typedef typename Tr::Finite_vertices_iterator Fvi; - for ( Fvi vit = tr_.finite_vertices_begin() ; - vit != tr_.finite_vertices_end() ; - ++ vit ) + typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); + + for ( Fvi vit = tr_.finite_vertices_begin(); vit != tr_.finite_vertices_end(); ++ vit ) { + const Weighted_point& position = tr_.point(vit); + typename std::map::const_iterator find_result = - value_map.find(tr_.geom_traits().construct_point_3_object()(vit->point())); + value_map.find(cp(position)); if ( find_result != value_map.end() ) { @@ -200,17 +203,16 @@ typename Mesh_sizing_field::FT Mesh_sizing_field:: operator()(const Bare_point& p, const Cell_handle& c) const { - typename Gt::Construct_weighted_point_3 p2wp = - tr_.geom_traits().construct_weighted_point_3_object(); + typename Gt::Construct_weighted_point_3 cwp = tr_.geom_traits().construct_weighted_point_3_object(); #ifdef CGAL_MESH_3_SIZING_FIELD_INEXACT_LOCATE //use the inexact locate (much faster than locate) to get a hint //and then use locate to check whether p is really inside hint // if not, an exact locate will be performed - Cell_handle hint = tr_.inexact_locate(p2wp(p),c); - const Cell_handle cell = tr_.locate(p2wp(p), hint); + Cell_handle hint = tr_.inexact_locate(cwp(p),c); + const Cell_handle cell = tr_.locate(cwp(p), hint); #else - const Cell_handle cell = tr_.locate(p2wp(p),c); + const Cell_handle cell = tr_.locate(cwp(p),c); #endif this->set_last_cell(cell); @@ -244,8 +246,8 @@ typename Mesh_sizing_field::FT Mesh_sizing_field:: interpolate_on_cell_vertices(const Bare_point& p, const Cell_handle& cell) const { + typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); typename Gt::Compute_volume_3 volume = tr_.geom_traits().compute_volume_3_object(); - typename Gt::Construct_point_3 wp2p = tr_.geom_traits().construct_point_3_object(); // Interpolate value using tet vertices values const FT& va = cell->vertex(0)->meshing_info(); @@ -253,10 +255,14 @@ interpolate_on_cell_vertices(const Bare_point& p, const Cell_handle& cell) const const FT& vc = cell->vertex(2)->meshing_info(); const FT& vd = cell->vertex(3)->meshing_info(); - const Bare_point& a = wp2p(tr_.point(cell, 0)); - const Bare_point& b = wp2p(tr_.point(cell, 1)); - const Bare_point& c = wp2p(tr_.point(cell, 2)); - const Bare_point& d = wp2p(tr_.point(cell, 3)); + const Weighted_point& wa = tr_.point(cell, 0); + const Weighted_point& wb = tr_.point(cell, 1); + const Weighted_point& wc = tr_.point(cell, 2); + const Weighted_point& wd = tr_.point(cell, 3); + const Bare_point& a = cp(wa); + const Bare_point& b = cp(wb); + const Bare_point& c = cp(wc); + const Bare_point& d = cp(wd); const FT abcp = CGAL::abs(volume(a,b,c,p)); const FT abdp = CGAL::abs(volume(a,d,b,p)); @@ -265,13 +271,12 @@ interpolate_on_cell_vertices(const Bare_point& p, const Cell_handle& cell) const // If volume is 0, then compute the average value if ( is_zero(abcp+abdp+acdp+bcdp) ) - return (va+vb+vc+vd)/4.; + return (va+vb+vc+vd) / 4.; return ( (abcp*vd + abdp*vc + acdp*vb + bcdp*va) / (abcp+abdp+acdp+bcdp) ); } - template typename Mesh_sizing_field::FT Mesh_sizing_field:: @@ -279,7 +284,7 @@ interpolate_on_facet_vertices(const Bare_point& p, const Cell_handle& cell) cons { typename Gt::Compute_area_3 area = tr_.geom_traits().compute_area_3_object(); - typename Gt::Construct_point_3 wp2p = tr_.geom_traits().construct_point_3_object(); + typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); // Find infinite vertex and put it in k0 int k0 = 0; int k1 = 1; @@ -298,13 +303,16 @@ interpolate_on_facet_vertices(const Bare_point& p, const Cell_handle& cell) cons const FT& vb = cell->vertex(k2)->meshing_info(); const FT& vc = cell->vertex(k3)->meshing_info(); - const Bare_point& a = wp2p(tr_.point(cell, k1)); - const Bare_point& b = wp2p(tr_.point(cell, k2)); - const Bare_point& c = wp2p(tr_.point(cell, k3)); + const Weighted_point& wa = tr_.point(cell, k1); + const Weighted_point& wb = tr_.point(cell, k2); + const Weighted_point& wc = tr_.point(cell, k3); + const Bare_point& a = cp(wa); + const Bare_point& b = cp(wb); + const Bare_point& c = cp(wc); - const FT abp = area(a,b,p); - const FT acp = area(a,c,p); - const FT bcp = area(b,c,p); + const FT abp = area(a, b, p); + const FT acp = area(a, c, p); + const FT bcp = area(b, c, p); CGAL_assertion(abp >= 0); CGAL_assertion(acp >= 0); diff --git a/Mesh_3/include/CGAL/Mesh_3/Odt_move.h b/Mesh_3/include/CGAL/Mesh_3/Odt_move.h index 9456a4041c4..31cf1cd57c6 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Odt_move.h +++ b/Mesh_3/include/CGAL/Mesh_3/Odt_move.h @@ -49,6 +49,7 @@ class Odt_move typedef typename Tr::Cell_handle Cell_handle; typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Weighted_point Weighted_point; typedef typename std::vector Facet_vector; typedef typename std::vector Cell_vector; @@ -73,10 +74,8 @@ public: // Compute move const Tr& tr = c3t3.triangulation(); - typename Gt::Construct_vector_3 vector = - tr.geom_traits().construct_vector_3_object(); - typename Gt::Construct_point_3 cp = - tr.geom_traits().construct_point_3_object(); + typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object(); + typename Gt::Construct_vector_3 vector = tr.geom_traits().construct_vector_3_object(); Vector_3 move = CGAL::NULL_VECTOR; FT sum_volume(0); @@ -95,8 +94,8 @@ public: Bare_point circumcenter = tr.dual(cell); // Compute move - const Bare_point& p = cp(tr.point(cell, cell->index(v))); - Vector_3 p_circum = vector(p, circumcenter); + const Weighted_point& p = tr.point(cell, cell->index(v)); + Vector_3 p_circum = vector(cp(p), circumcenter); FT volume = volume_quadrature(cell, tr, sizing_field); move = move + p_circum * volume; @@ -126,13 +125,11 @@ private: const Tr& tr, const Sizing_field& sizing_field) const { - typename Gt::Compute_volume_3 volume = - tr.geom_traits().compute_volume_3_object(); - typename Gt::Construct_centroid_3 centroid = - tr.geom_traits().construct_centroid_3_object(); + typename Gt::Construct_centroid_3 centroid = tr.geom_traits().construct_centroid_3_object(); + typename Gt::Compute_volume_3 volume = tr.geom_traits().compute_volume_3_object(); Bare_point c = centroid(tr.tetrahedron(cell)); - FT s = sizing_field(c,std::make_pair(cell,true)); + FT s = sizing_field(c, std::make_pair(cell, true)); CGAL_assertion(!is_zero(s)); // Points of cell are positively oriented @@ -206,10 +203,8 @@ private: // const Tr& tr, // const Sizing_field& sizing_field) const // { -// typename Gt::Compute_area_3 area = -// tr.geom_traits().compute_area_3_object(); -// typename Gt::Construct_centroid_3 centroid = -// tr.geom_traits().construct_centroid_3_object(); +// typename Gt::Compute_area_3 area = tr.geom_traits().compute_area_3_object(); +// typename Gt::Construct_centroid_3 centroid = tr.geom_traits().construct_centroid_3_object(); // // Bare_point c = centroid(tr.triangle(facet)); // FT s = sizing_field(c, facet.first->vertex(0)); @@ -224,15 +219,14 @@ private: // const Tr& tr, // const Sizing_field& sizing_field) const // { -// typename Gt::Construct_point_3 cp = -// tr.geom_traits().construct_point_3_object(); -// typename Gt::Compute_squared_distance_3 sq_distance = -// tr.geom_traits().compute_squared_distance_3_object(); -// typename Gt::Construct_midpoint_3 midpoint = -// tr.geom_traits().construct_midpoint_3_object(); +// typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object(); +// typename Gt::Construct_midpoint_3 midpoint = tr.geom_traits().construct_midpoint_3_object(); +// typename Gt::Compute_squared_distance_3 sq_distance = tr.geom_traits().compute_squared_distance_3_object(); // -// const Bare_point& p1 = cp(tr.point(cell, vertex_index_1)); -// const Bare_point& p2 = cp(tr.point(cell, vertex_index_2)); +// const Weighted_point& wp1 = tr.point(cell, vertex_index_1); +// const Weighted_point& wp2 = tr.point(cell, vertex_index_2); +// const Bare_point& p1 = cp(wp1); +// const Bare_point& p2 = cp(wp2); // // Bare_point c = midpoint(p1, p2); // FT s = sizing_field(c, cell->vertex(vertex_index_1)); @@ -265,10 +259,8 @@ private: // // Vector_3 normal_outside(const Facet& f, const C3T3& c3t3) const // { -// typename Gt::Construct_point_3 cp = -// c3t3.triangulation().geom_traits().construct_point_3_object(); -// typename Gt::Construct_normal_3 normal = -// c3t3.triangulation().geom_traits().construct_normal_3_object(); +// typename Gt::Construct_point_3 cp = c3t3.triangulation().geom_traits().construct_point_3_object(); +// typename Gt::Construct_normal_3 normal = c3t3.triangulation().geom_traits().construct_normal_3_object(); // // const Cell_handle& cell = f.first; // const int& i = f.second; @@ -285,11 +277,11 @@ private: // if ( ! c3t3.is_in_complex(cell) ) // std::swap(k1,k2); // -// const Bare_point& p1 = cp(c3t3.triangulation().point(cell, k1)); -// const Bare_point& p2 = cp(c3t3.triangulation().point(cell, k2)); -// const Bare_point& p3 = cp(c3t3.triangulation().point(cell, k3)); +// const Weighted_point& wp1 = c3t3.triangulation().point(cell, k1); +// const Weighted_point& wp2 = c3t3.triangulation().point(cell, k2); +// const Weighted_point& wp3 = c3t3.triangulation().point(cell, k3); // -// return normal(p1,p2,p3); +// return normal(cp(p1), cp(p2), cp(p3)); // } }; diff --git a/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h b/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h index 3be000b115f..23699189d72 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h +++ b/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h @@ -109,11 +109,11 @@ struct Get_Is_facet_bad { #ifdef _DEBUG int f1_current_erase_counter = CGAL::cpp11::get<0>(f).first->erase_counter(); int f1_saved_erase_counter = CGAL::cpp11::get<1>(f); - int f2_current_erase_counte7r = CGAL::cpp11::get<2>(f).first->erase_counter(); + int f2_current_erase_counter = CGAL::cpp11::get<2>(f).first->erase_counter(); int f2_saved_erase_counter = CGAL::cpp11::get<3>(f); //f1_current_erase_counter - f1_saved_erase_counter + f2_current_erase_counter - f2_saved_erase_counter == 1 - /*if (f1_current_erase_counter - f1_saved_erase_counter + f2_current_erase_counter - f2_saved_erase_counter == 1) + if (f1_current_erase_counter - f1_saved_erase_counter + f2_current_erase_counter - f2_saved_erase_counter == 1) { #ifdef CGAL_LINKED_WITH_TBB tbb::queuing_mutex mut; @@ -138,8 +138,8 @@ struct Get_Is_facet_bad { << "}" << std::endl; std::string s = sstr.str(); - //std::cerr << s << std::endl; - }*/ + std::cerr << s << std::endl; + } #endif return (CGAL::cpp11::get<0>(f).first->erase_counter() == CGAL::cpp11::get<1>(f) && CGAL::cpp11::get<2>(f).first->erase_counter() == CGAL::cpp11::get<3>(f) ); @@ -297,12 +297,16 @@ public: #ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS const Cell_handle c = facet.first; const int i = facet.second; - std::cerr << "Facet (" + std::cerr << "\n Facet (" << r_tr_.point(c, (i+1)&3) << " , " << r_tr_.point(c, (i+2)&3) << " , " << r_tr_.point(c, (i+3)&3) << ") : refinement point is " << get_facet_surface_center(facet) << std::endl; #endif + if(CGAL::squared_distance(get_facet_surface_center(facet), + r_tr_.point(c, (i+3)&3)) < 1e-6) + exit(0); + CGAL_assertion (this->is_facet_on_surface(facet)); this->set_last_vertex_index(get_facet_surface_center_index(facet)); return get_facet_surface_center(facet); @@ -349,9 +353,9 @@ public: { std::stringstream sstr; sstr << "Facet { " << std::endl - << " " << *facet.first->vertex((facet.second+1)%4) << std::endl - << " " << *facet.first->vertex((facet.second+2)%4) << std::endl - << " " << *facet.first->vertex((facet.second+3)%4) << std::endl + << " " << *facet.first->vertex((facet.second+1)%4) << std::endl + << " " << *facet.first->vertex((facet.second+2)%4) << std::endl + << " " << *facet.first->vertex((facet.second+3)%4) << std::endl << " 4th vertex in cell: " << *facet.first->vertex(facet.second) << std::endl << "}" << std::endl; @@ -973,9 +977,10 @@ scan_triangulation_impl() { // Cannot be const, see treat_new_facet signature Facet facet = *facet_it; - /*std::cerr << "*" << *facet.first->vertex((facet.second+1)%4) << std::endl - << " " << *facet.first->vertex((facet.second+2)%4) << std::endl - << " " << *facet.first->vertex((facet.second+3)%4) << std::endl;*/ + std::cerr << "TREAT FACET : " << std::endl + << "*" << *facet.first->vertex((facet.second+1)%4) << std::endl + << " " << *facet.first->vertex((facet.second+2)%4) << std::endl + << " " << *facet.first->vertex((facet.second+3)%4) << std::endl; this->treat_new_facet(facet); } } @@ -1388,10 +1393,10 @@ before_insertion_impl(const Facet& facet, bool source_facet_is_in_conflict = false; - /*std::cerr << "before_insertion_impl:" << std::endl - << "* " << *facet.first->vertex((facet.second+1)%4) << std::endl - << " " << *facet.first->vertex((facet.second+2)%4) << std::endl - << " " << *facet.first->vertex((facet.second+3)%4) << std::endl;*/ + std::cerr << "before_insertion_impl:" << std::endl + << "* " << *facet.first->vertex((facet.second+1)%4) << std::endl + << " " << *facet.first->vertex((facet.second+2)%4) << std::endl + << " " << *facet.first->vertex((facet.second+3)%4) << std::endl; // Iterate on conflict zone facets for (Facets_iterator facet_it = zone.internal_facets.begin(); @@ -1544,11 +1549,12 @@ treat_new_facet(Facet& facet) { insert_bad_facet(facet, *is_facet_bad); - /*std::cerr << "INSERT BAD FACET : " << std::endl - << "* " << *facet.first->vertex((facet.second+1)%4) << std::endl - << " " << *facet.first->vertex((facet.second+2)%4) << std::endl - << " " << *facet.first->vertex((facet.second+3)%4) << std::endl - << " Quality=" << is_facet_bad->second << std::endl;*/ + std::cerr << "INSERT BAD FACET : " << std::endl + << "* " << *facet.first->vertex((facet.second+1)%4) << std::endl + << " " << *facet.first->vertex((facet.second+2)%4) << std::endl + << " " << *facet.first->vertex((facet.second+3)%4) << std::endl + << " Surface center: " << get_facet_surface_center(facet) << std::endl + << " Quality = " << is_facet_bad->second << std::endl; } } else @@ -1741,9 +1747,9 @@ is_encroached_facet_refinable(Facet& facet) const ++wp_nb; } - const Weighted_point& p1 = c->vertex(k1)->point(); - const Weighted_point& p2 = c->vertex(k2)->point(); - const Weighted_point& p3 = c->vertex(k3)->point(); + const Weighted_point& p1 = r_tr_.point(c, k1); + const Weighted_point& p2 = r_tr_.point(c, k2); + const Weighted_point& p3 = r_tr_.point(c, k3); const FT min_ratio (0.16); // (0.2*2)^2 diff --git a/Mesh_3/include/CGAL/Mesh_3/Refine_facets_manifold_base.h b/Mesh_3/include/CGAL/Mesh_3/Refine_facets_manifold_base.h index ff95b2d9984..08bf217d111 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Refine_facets_manifold_base.h +++ b/Mesh_3/include/CGAL/Mesh_3/Refine_facets_manifold_base.h @@ -61,29 +61,31 @@ class Refine_facets_manifold_base Container_> Base; public: - typedef Complex3InTriangulation3 C3t3; - typedef MeshDomain Mesh_domain; - typedef typename Tr::Weighted_point Weighted_point; - typedef typename Tr::Bare_point Bare_point; - typedef typename Tr::Facet Facet; - typedef typename Tr::Vertex_handle Vertex_handle; + typedef Complex3InTriangulation3 C3t3; + typedef MeshDomain Mesh_domain; + + typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Weighted_point Weighted_point; + + typedef typename Tr::Vertex_handle Vertex_handle; + typedef typename Tr::Facet Facet; typedef typename Triangulation_mesher_level_traits_3::Zone Zone; protected: - typedef typename Tr::Geom_traits GT; - typedef typename GT::FT FT; - typedef typename GT::Construct_point_3 Construct_point_3; - typedef typename Tr::Edge Edge; - typedef typename Tr::Cell_handle Cell_handle; + typedef typename Tr::Geom_traits GT; + typedef typename GT::FT FT; - typedef typename Tr::Facet_circulator Tr_facet_circulator; - typedef std::pair EdgeVV; + typedef typename Tr::Edge Edge; + typedef typename Tr::Cell_handle Cell_handle; + + typedef typename Tr::Facet_circulator Tr_facet_circulator; + typedef std::pair EdgeVV; protected: typedef ::boost::bimap< EdgeVV, - ::boost::bimaps::multiset_of > Bad_edges; - typedef typename Bad_edges::value_type Bad_edge; + ::boost::bimaps::multiset_of > Bad_edges; + typedef typename Bad_edges::value_type Bad_edge; mutable Bad_edges m_bad_edges; mutable std::set m_bad_vertices; @@ -127,22 +129,21 @@ private: } FT compute_sq_distance_to_facet_center(const Facet& f, - const Vertex_handle v) const { - Construct_point_3 wp2p = this->r_tr_.geom_traits().construct_point_3_object(); - const Bare_point& fcenter = f.first->get_facet_surface_center(f.second); - const Bare_point& vpoint = wp2p(v->point()); + const Vertex_handle v) const + { + typename GT::Construct_point_3 cp = this->r_tr_.geom_traits().construct_point_3_object(); - return - this->r_tr_.geom_traits().compute_squared_distance_3_object()(fcenter, - vpoint) - - v->point().weight(); + const Bare_point& fcenter = f.first->get_facet_surface_center(f.second); + const Weighted_point& wp = this->r_tr_.point(v); + + return this->r_tr_.min_squared_distance(fcenter, cp(wp)) - v->point().weight(); } Facet biggest_incident_facet_in_complex(const Vertex_handle v) const { #ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS - std::cerr << "Bad vertex: " << v->point() << std::endl; + std::cerr << "Bad vertex: " << this->r_tr_.point(v) << std::endl; #endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS std::vector facets; facets.reserve(64); @@ -201,8 +202,8 @@ private: // use the list of incident facets in the complex Vertex_handle fev = edge_to_edgevv(arete).first; #ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS - std::cerr << "Bad edge: (" << fev->point() - << ", " << arete.first->vertex(arete.third)->point() + std::cerr << "Bad edge: (" << this->r_tr_.point(fev) + << ", " << this->r_tr_.point(arete.first->vertex(arete.third)) << ")\n incident facets squared sizes:\n"; #endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS Tr_facet_circulator fcirc = this->r_tr_.incident_facets(arete); @@ -290,7 +291,7 @@ private: if(m_bad_vertices.erase(c->vertex(j)) > 0) { #ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS std::cerr << "m_bad_vertices.erase(" - << c->vertex(j)->point() << ")\n"; + << this->r_tr_.point(c, j) << ")\n"; #endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS } } @@ -345,18 +346,18 @@ public: (this->r_c3t3_.face_status(*eit) == C3t3::BOUNDARY) ) ) { #ifdef CGAL_LINKED_WITH_TBB - // Parallel - if (boost::is_convertible::value) - { - this->insert_bad_facet(biggest_incident_facet_in_complex(*eit), - typename Base::Quality()); - } else + // Parallel + if (boost::is_convertible::value) + { + this->insert_bad_facet(biggest_incident_facet_in_complex(*eit), + typename Base::Quality()); + } else #endif // CGAL_LINKED_WITH_TBB - { // Sequential - m_bad_edges.insert(Bad_edge(edge_to_edgevv(*eit), - (this->r_c3t3_.face_status(*eit) == - C3t3::SINGULAR ? 0 : 1))); - } + { // Sequential + m_bad_edges.insert(Bad_edge(edge_to_edgevv(*eit), + (this->r_c3t3_.face_status(*eit) == + C3t3::SINGULAR ? 0 : 1))); + } #ifdef CGAL_MESH_3_VERBOSE ++n; #endif @@ -385,19 +386,19 @@ public: if( this->r_c3t3_.face_status(vit) == C3t3::SINGULAR ) { #ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS std::cerr << "m_bad_vertices.insert(" - << vit->point() << ")\n"; + << this->r_tr_.point(vit) << ")\n"; #endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS #ifdef CGAL_LINKED_WITH_TBB - // Parallel - if (boost::is_convertible::value) - { - this->insert_bad_facet(biggest_incident_facet_in_complex(vit), - typename Base::Quality()); - } else + // Parallel + if (boost::is_convertible::value) + { + this->insert_bad_facet(biggest_incident_facet_in_complex(vit), + typename Base::Quality()); + } else #endif // CGAL_LINKED_WITH_TBB - { // Sequential - m_bad_vertices.insert( vit ); - } + { // Sequential + m_bad_vertices.insert( vit ); + } #ifdef CGAL_MESH_3_VERBOSE ++n; #endif @@ -410,7 +411,7 @@ public: std::cerr << "Bad vertices queue:\n"; BOOST_FOREACH(Vertex_handle v2, m_bad_vertices) { - std::cerr << v2->point() << std::endl; + std::cerr << this->r_tr_.point(v2) << std::endl; } CGAL::dump_c3t3(this->r_c3t3_, "dump-at-scan-vertices"); # endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS @@ -456,9 +457,9 @@ public: #ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS const EdgeVV& edgevv = m_bad_edges.right.begin()->second; std::cerr << "Bad edge " - << edgevv.first->point() + << this->r_tr_.point(edgevv.first) << " - " - << edgevv.second->point() + << this->r_tr_.point(edgevv.second) << "\n"; #endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS return biggest_incident_facet_in_complex(first_bad_edge); @@ -469,9 +470,9 @@ public: std::cerr << "Bad vertices queue:\n"; BOOST_FOREACH(Vertex_handle v2, m_bad_vertices) { - std::cerr << v2->point() << std::endl; + std::cerr << this->r_tr_.point(v2) << std::endl; } - std::cerr << "Bad vertex " << v->point() << "\n"; + std::cerr << "Bad vertex " << this->r_tr_.point(v) << "\n"; #endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS CGAL_assertion(this->r_c3t3_.has_incident_facets_in_complex(v)); if(this->r_c3t3_.face_status(v) != C3t3::SINGULAR) { @@ -541,27 +542,27 @@ public: ) { #ifdef CGAL_LINKED_WITH_TBB - // Parallel - if (boost::is_convertible::value) - { - this->insert_bad_facet(biggest_incident_facet_in_complex(edge), - typename Base::Quality()); - } else + // Parallel + if (boost::is_convertible::value) + { + this->insert_bad_facet(biggest_incident_facet_in_complex(edge), + typename Base::Quality()); + } else #endif // CGAL_LINKED_WITH_TBB - { // Sequential - m_bad_edges.insert(Bad_edge(edge_to_edgevv(edge), - (this->r_c3t3_.face_status(edge) == - C3t3::SINGULAR ? 0 : 1))); - } + { // Sequential + m_bad_edges.insert(Bad_edge(edge_to_edgevv(edge), + (this->r_c3t3_.face_status(edge) == + C3t3::SINGULAR ? 0 : 1))); + } } else { #ifdef CGAL_LINKED_WITH_TBB - // Sequential only - if (!boost::is_convertible::value) + // Sequential only + if (!boost::is_convertible::value) #endif // CGAL_LINKED_WITH_TBB - { - m_bad_edges.left.erase( edge_to_edgevv(edge) ); // @TODO: pourquoi?! - } + { + m_bad_edges.left.erase( edge_to_edgevv(edge) ); // @TODO: pourquoi?! + } } } } @@ -590,19 +591,19 @@ public: { #ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS std::cerr << "m_bad_vertices.insert(" - << (*vit)->point() << ")\n"; + << this->r_tr_.point(*vit) << ")\n"; #endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS #ifdef CGAL_LINKED_WITH_TBB - // Parallel - if (boost::is_convertible::value) - { - this->insert_bad_facet(biggest_incident_facet_in_complex(*vit), - typename Base::Quality()); - } else + // Parallel + if (boost::is_convertible::value) + { + this->insert_bad_facet(biggest_incident_facet_in_complex(*vit), + typename Base::Quality()); + } else #endif // CGAL_LINKED_WITH_TBB - { // Sequential - m_bad_vertices.insert(*vit); - } + { // Sequential + m_bad_vertices.insert(*vit); + } } } @@ -613,19 +614,19 @@ public: { #ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS std::cerr << "m_bad_vertices.insert(" - << v->point() << ")\n"; + << this->r_tr_.point(v) << ")\n"; #endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS #ifdef CGAL_LINKED_WITH_TBB - // Parallel - if (boost::is_convertible::value) - { - this->insert_bad_facet(biggest_incident_facet_in_complex(v), - typename Base::Quality()); - } else + // Parallel + if (boost::is_convertible::value) + { + this->insert_bad_facet(biggest_incident_facet_in_complex(v), + typename Base::Quality()); + } else #endif // CGAL_LINKED_WITH_TBB - { // Sequential - m_bad_vertices.insert(v); - } + { // Sequential + m_bad_vertices.insert(v); + } } } diff --git a/Mesh_3/include/CGAL/Mesh_3/Sliver_perturber.h b/Mesh_3/include/CGAL/Mesh_3/Sliver_perturber.h index 308c3c365de..d3175ac28d8 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Sliver_perturber.h +++ b/Mesh_3/include/CGAL/Mesh_3/Sliver_perturber.h @@ -472,13 +472,15 @@ class Sliver_perturber typedef Sliver_perturber_base< typename C3T3::Triangulation, Concurrency_tag> Base; - typedef typename C3T3::Triangulation Tr; - typedef typename Tr::Geom_traits Gt; + typedef typename C3T3::Triangulation Tr; + typedef typename Tr::Geom_traits Gt; typedef typename Tr::Cell_handle Cell_handle; typedef typename Base::Vertex_handle Vertex_handle; typedef typename Tr::Vertex Vertex; - typedef typename MeshDomain::Point_3 Point_3; + + typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Weighted_point Weighted_point; typedef typename std::vector Cell_vector; typedef typename std::vector Vertex_vector; @@ -507,7 +509,7 @@ private: typedef PVertex_ PVertex; @@ -1298,7 +1300,7 @@ perturb_vertex( PVertex pv , bool *could_lock_zone ) const { - typename Gt::Construct_point_3 wp2p = tr_.geom_traits().construct_point_3_object(); + typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); #ifdef CGAL_CONCURRENT_MESH_3_PROFILING static Profile_branch_counter_3 bcounter( @@ -1313,9 +1315,10 @@ perturb_vertex( PVertex pv return; } - Point_3 p = wp2p(pv.vertex()->point()); + // and what's happening below... ? + Point_3 p = cp(pv.vertex()->point()); if (!helper_.try_lock_point_no_spin(pv.vertex()->point()) || - ! tr_.geom_traits().equal_3_object()(p, wp2p(pv.vertex()->point()))) + ! tr_.geom_traits().equal_3_object()(p, cp(pv.vertex()->point()))) { #ifdef CGAL_CONCURRENT_MESH_3_PROFILING bcounter.increment_branch_2(); // THIS is an early withdrawal! diff --git a/Mesh_3/include/CGAL/Mesh_3/Slivers_exuder.h b/Mesh_3/include/CGAL/Mesh_3/Slivers_exuder.h index d422503159d..374760eafb1 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Slivers_exuder.h +++ b/Mesh_3/include/CGAL/Mesh_3/Slivers_exuder.h @@ -99,34 +99,34 @@ namespace Mesh_3 { // Then, its operator() takes an other vertex handle v2 as input, and // returns the distance d(v1, v2). // It is used in Slivers_exuder, to constructor a transform iterator. - template + template class Min_distance_from_v : - public CGAL::unary_function + public CGAL::unary_function { const Vertex_handle * v; - const Gt& gt; + const Tr& tr; double & dist; public: Min_distance_from_v(const Vertex_handle& vh, double& dist, - const Gt& geom_traits = Gt()) - : v(&vh), gt(geom_traits), dist(dist) + const Tr& tr) + : v(&vh), tr(tr), dist(dist) { } void operator()(const Vertex_handle& vh) const { - typedef typename Gt::Compute_squared_distance_3 Compute_squared_distance_3; - Compute_squared_distance_3 distance = gt.compute_squared_distance_3_object(); + typedef typename Tr::Geom_traits::Construct_point_3 Construct_point_3; + typedef typename Tr::Weighted_point Weighted_point; - typedef typename Gt::Construct_point_3 Construct_point_3; - Construct_point_3 wp2p = gt.construct_point_3_object(); + Construct_point_3 cp = tr.geom_traits().construct_point_3_object(); - // - const double d = CGAL::to_double(distance(wp2p((*v)->point()), - wp2p(vh->point()))); + const Weighted_point& wpv = tr.point(*v); + const Weighted_point& wpvh = tr.point(vh); + + const double d = CGAL::to_double(tr.min_squared_distance(cp(wpv), cp(wpvh))); if(d < dist){ dist = d; } @@ -645,7 +645,7 @@ private: { double dist = (std::numeric_limits::max)(); - details::Min_distance_from_v min_distance_from_v(vh, dist, tr_.geom_traits()); + details::Min_distance_from_v min_distance_from_v(vh, dist, tr_); tr_.adjacent_vertices(vh, boost::make_function_output_iterator(min_distance_from_v)); @@ -1126,8 +1126,10 @@ pump_vertex(const Vertex_handle& pumped_vertex, // If best_weight < pumped_vertex weight, nothing to do if ( best_weight > pumped_vertex->point().weight() ) { - typename Gt::Construct_point_3 wp2p = tr_.geom_traits().construct_point_3_object(); - Weighted_point wp(wp2p(pumped_vertex->point()), best_weight); + typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); + + const Weighted_point& pwp = tr_.point(pumped_vertex); + Weighted_point wp(cp(pwp), best_weight); // Insert weighted point into mesh // note it can fail if the mesh is non-manifold at pumped_vertex @@ -1205,6 +1207,8 @@ expand_prestar(const Cell_handle& cell_to_add, Pre_star& pre_star, Sliver_values& criterion_values) const { + typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); + // Delete first facet of pre_star Facet start_facet = pre_star.front()->second; CGAL_assertion(tr_.mirror_facet(start_facet).first == cell_to_add); @@ -1278,14 +1282,15 @@ expand_prestar(const Cell_handle& cell_to_add, // Update ratio (ratio is needed for cells of complex only) if ( c3t3_.is_in_complex(cell_to_add) ) { - typename Gt::Construct_point_3 wp2p = tr_.geom_traits().construct_point_3_object(); - Tetrahedron_3 tet(wp2p(pumped_vertex->point()), - wp2p(cell_to_add->vertex((i+1)&3)->point()), - wp2p(cell_to_add->vertex((i+2)&3)->point()), - wp2p(cell_to_add->vertex((i+3)&3)->point())); + const Weighted_point& pwp = tr_.point(pumped_vertex); + const Weighted_point& fwp1 = tr_.point(cell_to_add, (i+1)&3); + const Weighted_point& fwp2 = tr_.point(cell_to_add, (i+2)&3); + const Weighted_point& fwp3 = tr_.point(cell_to_add, (i+3)&3); + + const Tetrahedron_3 tet(cp(pwp), cp(fwp1), cp(fwp2), cp(fwp3)); double new_value = sliver_criteria_(tet); - criterion_values.insert(std::make_pair(current_facet,new_value)); + criterion_values.insert(std::make_pair(current_facet, new_value)); } } } @@ -1369,7 +1374,7 @@ get_best_weight(const Vertex_handle& v, bool *could_lock_zone) const #ifdef CGAL_MESH_3_DEBUG_SLIVERS_EXUDER if ( best_weight > v->point().weight() ) { - Weighted_point wp(v->point(), best_weight); + Weighted_point wp(tr_.point(v), best_weight); check_pre_star(pre_star_copy, wp, v); check_ratios(ratios_copy, wp, v); } @@ -1794,7 +1799,7 @@ check_ratios(const Sliver_values& criterion_values, Facet_vector internal_facets; Facet_vector boundary_facets; - typename Gt::Construct_point_3 wp2p = tr_.geom_traits().construct_point_3_object(); + typename Gt::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object(); tr_.find_conflicts(wp, vh->cell(), @@ -1821,10 +1826,13 @@ check_ratios(const Sliver_values& criterion_values, continue; int k = it->second; - Tetrahedron_3 tet(wp2p(vh->point()), - wp2p(it->first->vertex((k+1)&3)->point()), - wp2p(it->first->vertex((k+2)&3)->point()), - wp2p(it->first->vertex((k+3)&3)->point())); + + const Weighted_point& vhwp = tr_.point(vh); + const Weighted_point& fwp1 = tr_.point(it->first, (k+1)&3)); + const Weighted_point& fwp2 = tr_.point(it->first, (k+2)&3)); + const Weighted_point& fwp3 = tr_.point(it->first, (k+3)&3)); + + Tetrahedron_3 tet(cp(vhwp), cp(fwp1), cp(fwp2), cp(fwp3)); double ratio = sliver_criteria_(tet); expected_ratios.push_back(ratio); diff --git a/Mesh_3/include/CGAL/Mesh_3/Triangulation_helpers.h b/Mesh_3/include/CGAL/Mesh_3/Triangulation_helpers.h index 0e00f245027..dbc4a6b930b 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Triangulation_helpers.h +++ b/Mesh_3/include/CGAL/Mesh_3/Triangulation_helpers.h @@ -172,7 +172,7 @@ no_topological_change(Tr& tr, tr.geom_traits().construct_opposite_vector_3_object(); bool np = true; - const Weighted_point fp = v0->point(); + const Weighted_point fp = tr.point(v0); //#define CGAL_PERIODIC_BACKUP_CHECK //#define CGAL_PERIODIC_SIDE_OF_DEBUG @@ -528,12 +528,16 @@ inside_protecting_balls(const Tr& tr, { Vertex_handle nv = tr.nearest_power_vertex(p, v->cell()); if(nv->point().weight() > 0) - return tr.geom_traits().compare_squared_distance_3_object()( - p, nv->point(), nv->point().weight()) != CGAL::LARGER; + { + typename Tr::Geom_traits::Construct_point_3 cp = tr.geom_traits().construct_point_3_object(); + const Weighted_point& nvwp = tr.point(nv); + return (tr.min_squared_distance(p, cp(nvwp)) <= nv->point().weight()); + } + return false; } - + /// This function well_oriented is called by no_topological_change after the /// position of the vertex has been (tentatively) modified. template @@ -544,7 +548,7 @@ well_oriented(const Tr& tr, { typedef typename Tr::Geom_traits Gt; typename Gt::Orientation_3 orientation = tr.geom_traits().orientation_3_object(); - typename Gt::Construct_point_3 wp2p = tr.geom_traits().construct_point_3_object(); + typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object(); typename Cell_vector::const_iterator it = cells_tos.begin(); for( ; it != cells_tos.end() ; ++it) @@ -555,17 +559,25 @@ well_oriented(const Tr& tr, int iv = c->index(tr.infinite_vertex()); Cell_handle cj = c->neighbor(iv); int mj = tr.mirror_index(c, iv); - if(orientation(wp2p(tr.point(cj, mj)), - wp2p(tr.point(c, (iv+1)&3)), - wp2p(tr.point(c, (iv+2)&3)), - wp2p(tr.point(c, (iv+3)&3))) != CGAL::NEGATIVE) + + const Weighted_point& mjwp = tr.point(cj, mj); + const Weighted_point& fwp1 = tr.point(c, (iv+1)&3); + const Weighted_point& fwp2 = tr.point(c, (iv+2)&3); + const Weighted_point& fwp3 = tr.point(c, (iv+3)&3); + + if(orientation(cp(mjwp), cp(fwp1), cp(fwp2), cp(fwp3)) != CGAL::NEGATIVE) return false; } - else if(orientation(wp2p(tr.point(c, 0)), - wp2p(tr.point(c, 1)), - wp2p(tr.point(c, 2)), - wp2p(tr.point(c, 3))) != CGAL::POSITIVE) + else + { + const Weighted_point& cwp0 = tr.point(c, 0); + const Weighted_point& cwp1 = tr.point(c, 1); + const Weighted_point& cwp2 = tr.point(c, 2); + const Weighted_point& cwp3 = tr.point(c, 3); + + if(orientation(cp(cwp0), cp(cwp1), cp(cwp2), cp(cwp3)) != CGAL::POSITIVE) return false; + } } return true; } @@ -582,7 +594,7 @@ well_oriented(const Tr& tr, { typedef typename Tr::Geom_traits Gt; typename Gt::Orientation_3 orientation = tr.geom_traits().orientation_3_object(); - typename Gt::Construct_point_3 wp2p = tr.geom_traits().construct_point_3_object(); + typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object(); typename Cell_vector::const_iterator it = cells_tos.begin(); for( ; it != cells_tos.end() ; ++it) @@ -593,16 +605,16 @@ well_oriented(const Tr& tr, int iv = c->index(tr.infinite_vertex()); Cell_handle cj = c->neighbor(iv); int mj = tr.mirror_index(c, iv); - if(orientation(wp2p(pg(cj->vertex(mj))), - wp2p(pg(c->vertex((iv+1)&3))), - wp2p(pg(c->vertex((iv+2)&3))), - wp2p(pg(c->vertex((iv+3)&3)))) != CGAL::NEGATIVE) + if(orientation(cp(pg(cj->vertex(mj))), + cp(pg(c->vertex((iv+1)&3))), + cp(pg(c->vertex((iv+2)&3))), + cp(pg(c->vertex((iv+3)&3)))) != CGAL::NEGATIVE) return false; } - else if(orientation(wp2p(pg(c->vertex(0))), - wp2p(pg(c->vertex(1))), - wp2p(pg(c->vertex(2))), - wp2p(pg(c->vertex(3)))) != CGAL::POSITIVE) + else if(orientation(cp(pg(c->vertex(0))), + cp(pg(c->vertex(1))), + cp(pg(c->vertex(2))), + cp(pg(c->vertex(3)))) != CGAL::POSITIVE) return false; } return true; diff --git a/Mesh_3/include/CGAL/Mesh_3/Triangulation_sizing_field.h b/Mesh_3/include/CGAL/Mesh_3/Triangulation_sizing_field.h index 4f9afa2b744..c5dd68c525e 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Triangulation_sizing_field.h +++ b/Mesh_3/include/CGAL/Mesh_3/Triangulation_sizing_field.h @@ -40,7 +40,6 @@ namespace CGAL { namespace Mesh_3 { - /** * @class Triangulation_sizing_field @@ -49,20 +48,21 @@ template class Triangulation_sizing_field { // Types - typedef typename Tr::Geom_traits Gt; - typedef typename Tr::Weighted_point Weighted_point; - typedef typename Gt::FT FT; - + typedef typename Tr::Geom_traits Gt; + typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Weighted_point Weighted_point; + typedef typename Gt::FT FT; + typedef Triangulation_vertex_base_with_info_3 Vb; typedef Triangulation_cell_base_3 Cb; typedef Triangulation_data_structure_3 Tds; typedef Regular_triangulation_3 Compact_triangulation; typedef Compact_triangulation Ctr; - typedef typename Tr::Vertex_handle Vertex_handle; - typedef typename Tr::Cell_handle Cell_handle; - typedef typename Ctr::Cell_handle CCell_handle; - typedef typename Ctr::Vertex_handle CVertex_handle; + typedef typename Tr::Vertex_handle Vertex_handle; + typedef typename Tr::Cell_handle Cell_handle; + typedef typename Ctr::Cell_handle CCell_handle; + typedef typename Ctr::Vertex_handle CVertex_handle; public: // Vertices of mesh triangulation do not need to be updated @@ -116,8 +116,9 @@ private: * Used by boost transform iterator */ struct Extract_point : - public CGAL::unary_function + public CGAL::unary_function { + // Weighted_point operator()(const typename Tr::Vertex& v) const { return v.point(); } }; @@ -152,7 +153,7 @@ fill(const std::map& value_map) ++ vit ) { typename std::map::const_iterator find_result = - value_map.find(vit->point()); + value_map.find(ctr_.point(vit)); if ( find_result != value_map.end() ) vit->info() = find_result->second; @@ -183,20 +184,24 @@ typename Triangulation_sizing_field::FT Triangulation_sizing_field:: interpolate_on_cell_vertices(const Weighted_point& p, const CCell_handle& cell) const { - typename Gt::Compute_volume_3 volume = - ctr_.geom_traits().compute_volume_3_object(); - + typename Gt::Construct_point_3 cp = ctr_.geom_traits().construct_point_3_object(); + typename Gt::Compute_volume_3 volume = ctr_.geom_traits().compute_volume_3_object(); + // Interpolate value using tet vertices values const FT& va = cell->vertex(0)->info(); const FT& vb = cell->vertex(1)->info(); const FT& vc = cell->vertex(2)->info(); const FT& vd = cell->vertex(3)->info(); - - const Weighted_point& a = cell->vertex(0)->point(); - const Weighted_point& b = cell->vertex(1)->point(); - const Weighted_point& c = cell->vertex(2)->point(); - const Weighted_point& d = cell->vertex(3)->point(); - + + const Weighted_point& wa = ctr_.point(cell, 0); + const Weighted_point& wb = ctr_.point(cell, 1); + const Weighted_point& wc = ctr_.point(cell, 2); + const Weighted_point& wd = ctr_.point(cell, 3); + const Bare_point& a = cp(wa); + const Bare_point& b = cp(wb); + const Bare_point& c = cp(wc); + const Bare_point& d = cp(wd); + const FT abcp = CGAL::abs(volume(a,b,c,p)); const FT abdp = CGAL::abs(volume(a,b,d,p)); const FT acdp = CGAL::abs(volume(a,c,d,p)); @@ -215,8 +220,8 @@ typename Triangulation_sizing_field::FT Triangulation_sizing_field:: interpolate_on_facet_vertices(const Weighted_point& p, const CCell_handle& cell) const { - typename Gt::Compute_area_3 area = - ctr_.geom_traits().compute_area_3_object(); + typename Gt::Construct_point_3 cp = ctr_.geom_traits().construct_point_3_object(); + typename Gt::Compute_area_3 area = ctr_.geom_traits().compute_area_3_object(); // Find infinite vertex and put it in k0 int k0 = 0; @@ -236,10 +241,13 @@ interpolate_on_facet_vertices(const Weighted_point& p, const CCell_handle& cell) const FT& vb = cell->vertex(k2)->info(); const FT& vc = cell->vertex(k3)->info(); - const Weighted_point& a = cell->vertex(k1)->point(); - const Weighted_point& b = cell->vertex(k2)->point(); - const Weighted_point& c = cell->vertex(k3)->point(); - + const Weighted_point& wa = ctr_.point(cell, k1); + const Weighted_point& wb = ctr_.point(cell, k2); + const Weighted_point& wc = ctr_.point(cell, k3); + const Bare_point& a = cp(wa); + const Bare_point& b = cp(wb); + const Bare_point& c = cp(wc); + const FT abp = area(a,b,p); const FT acp = area(a,c,p); const FT bcp = area(b,c,p); diff --git a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h index 4d64d768ad3..8d3620b1bd7 100644 --- a/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h +++ b/Mesh_3/include/CGAL/Mesh_3/initialize_triangulation_from_labeled_image.h @@ -84,21 +84,21 @@ void initialize_triangulation_from_labeled_image(C3T3& c3t3, ) { typedef typename C3T3::Triangulation Tr; + typedef typename Tr::Geom_traits Gt; typedef typename Tr::Weighted_point Weighted_point; typedef typename Tr::Bare_point Bare_point; typedef typename Tr::Segment Segment_3; - typedef typename Tr::Geom_traits::Vector_3 Vector_3; typedef typename Tr::Vertex_handle Vertex_handle; typedef typename Tr::Cell_handle Cell_handle; + typedef typename Gt::Vector_3 Vector_3; + typedef MeshDomain Mesh_domain; Tr& tr = c3t3.triangulation(); - typename Tr::Geom_traits::Construct_point_3 wp2p = - tr.geom_traits().construct_point_3_object(); - typename Tr::Geom_traits::Construct_weighted_point_3 p2wp = - tr.geom_traits().construct_weighted_point_3_object(); + typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object(); + typename Gt::Construct_weighted_point_3 cwp = tr.geom_traits().construct_weighted_point_3_object(); if(protect_features) { init_tr_from_labeled_image_call_init_features @@ -155,7 +155,8 @@ void initialize_triangulation_from_labeled_image(C3T3& c3t3, construct_intersection(Segment_3(it->first, test)); if (CGAL::cpp11::get<2>(intersect) != 0) { - Weighted_point pi = p2wp(CGAL::cpp11::get<0>(intersect)); + const Bare_point& bpi = CGAL::cpp11::get<0>(intersect); + Weighted_point pi = cwp(bpi); // This would cause trouble to optimizers // check pi will not be hidden @@ -201,8 +202,9 @@ void initialize_triangulation_from_labeled_image(C3T3& c3t3, { if (cv->point().weight() == 0.) continue; - if (CGAL::compare_squared_distance(pi.point(), wp2p(cv->point()), cv->point().weight()) - != CGAL::LARGER) + + const Weighted_point& cvwp = tr.point(cv); + if (tr.min_squared_distance(bpi, cp(cvwp) <= cv->point().weight())) { pi_inside_protecting_sphere = true; break; diff --git a/Mesh_3/include/CGAL/Mesh_3/mesh_standard_cell_criteria.h b/Mesh_3/include/CGAL/Mesh_3/mesh_standard_cell_criteria.h index a15a30904a3..78671cd9ebf 100644 --- a/Mesh_3/include/CGAL/Mesh_3/mesh_standard_cell_criteria.h +++ b/Mesh_3/include/CGAL/Mesh_3/mesh_standard_cell_criteria.h @@ -79,20 +79,26 @@ protected: virtual Badness do_is_bad(const Tr& tr, const Cell_handle& ch) const { - typedef typename Tr::Bare_point Bare_point; typedef typename Tr::Geom_traits Geom_traits; + typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Weighted_point Weighted_point; typedef typename Geom_traits::Compute_squared_distance_3 Distance; typedef typename Geom_traits::Compute_squared_radius_3 Radius; typedef typename Geom_traits::Construct_point_3 Construct_point_3; + Distance distance = tr.geom_traits().compute_squared_distance_3_object(); Radius sq_radius = tr.geom_traits().compute_squared_radius_3_object(); - Construct_point_3 wp2p = tr.geom_traits().construct_point_3_object(); + Construct_point_3 cp = tr.geom_traits().construct_point_3_object(); - const Bare_point& p = wp2p(tr.point(ch, 0)); - const Bare_point& q = wp2p(tr.point(ch, 1)); - const Bare_point& r = wp2p(tr.point(ch, 2)); - const Bare_point& s = wp2p(tr.point(ch, 3)); + const Weighted_point& wp = tr.point(ch, 0); + const Weighted_point& wq = tr.point(ch, 1); + const Weighted_point& wr = tr.point(ch, 2); + const Weighted_point& ws = tr.point(ch, 3); + const Bare_point& p = cp(wp); + const Bare_point& q = cp(wq); + const Bare_point& r = cp(wr); + const Bare_point& s = cp(ws); const FT size = sq_radius(p, q, r, s); @@ -164,18 +170,23 @@ protected: virtual Badness do_is_bad(const Tr& tr, const Cell_handle& ch) const { - typedef typename Tr::Bare_point Bare_point; - typedef typename Tr::Geom_traits Geom_traits; + typedef typename Tr::Geom_traits Geom_traits; + typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Weighted_point Weighted_point; typedef typename Geom_traits::Compute_squared_radius_3 Radius; typedef typename Geom_traits::Construct_point_3 Construct_point_3; Radius sq_radius = tr.geom_traits().compute_squared_radius_3_object(); - Construct_point_3 wp2p = tr.geom_traits().construct_point_3_object(); + Construct_point_3 cp = tr.geom_traits().construct_point_3_object(); - const Bare_point& p = wp2p(tr.point(ch, 0)); - const Bare_point& q = wp2p(tr.point(ch, 1)); - const Bare_point& r = wp2p(tr.point(ch, 2)); - const Bare_point& s = wp2p(tr.point(ch, 3)); + const Weighted_point& wp = tr.point(ch, 0); + const Weighted_point& wq = tr.point(ch, 1); + const Weighted_point& wr = tr.point(ch, 2); + const Weighted_point& ws = tr.point(ch, 3); + const Bare_point& p = cp(wp); + const Bare_point& q = cp(wq); + const Bare_point& r = cp(wr); + const Bare_point& s = cp(ws); const FT size = sq_radius(p, q, r, s); @@ -234,18 +245,23 @@ protected: virtual Badness do_is_bad(const Tr& tr, const Cell_handle& ch) const { - typedef typename Tr::Bare_point Bare_point; - typedef typename Tr::Geom_traits Geom_traits; + typedef typename Tr::Geom_traits Geom_traits; + typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Weighted_point Weighted_point; typedef typename Geom_traits::Compute_squared_radius_3 Radius; typedef typename Geom_traits::Construct_point_3 Construct_point_3; Radius sq_radius = tr.geom_traits().compute_squared_radius_3_object(); - Construct_point_3 wp2p = tr.geom_traits().construct_point_3_object(); + Construct_point_3 cp = tr.geom_traits().construct_point_3_object(); - const Bare_point& p = wp2p(tr.point(ch, 0)); - const Bare_point& q = wp2p(tr.point(ch, 1)); - const Bare_point& r = wp2p(tr.point(ch, 2)); - const Bare_point& s = wp2p(tr.point(ch, 3)); + const Weighted_point& wp = tr.point(ch, 0); + const Weighted_point& wq = tr.point(ch, 1); + const Weighted_point& wr = tr.point(ch, 2); + const Weighted_point& ws = tr.point(ch, 3); + const Bare_point& p = cp(wp); + const Bare_point& q = cp(wq); + const Bare_point& r = cp(wr); + const Bare_point& s = cp(ws); const FT size = sq_radius(p, q, r, s); const FT sq_bound = CGAL::square( size_(tr.dual(ch), 3, diff --git a/Mesh_3/include/CGAL/Mesh_3/mesh_standard_facet_criteria.h b/Mesh_3/include/CGAL/Mesh_3/mesh_standard_facet_criteria.h index 1630db6e57a..7540b4dff04 100644 --- a/Mesh_3/include/CGAL/Mesh_3/mesh_standard_facet_criteria.h +++ b/Mesh_3/include/CGAL/Mesh_3/mesh_standard_facet_criteria.h @@ -110,8 +110,9 @@ protected: CGAL_assertion (f.first->is_facet_on_surface(f.second)); CGAL_assertion (B_ != 0); - typedef typename Tr::Geom_traits Gt; - typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Geom_traits Gt; + typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Weighted_point Weighted_point; typedef typename Gt::Compute_squared_area_3 Area; typedef typename Gt::Compute_squared_distance_3 Distance; @@ -120,12 +121,15 @@ protected: Area area = tr.geom_traits().compute_squared_area_3_object(); Distance distance = tr.geom_traits().compute_squared_distance_3_object(); - Construct_point_3 wp2p = tr.geom_traits().construct_point_3_object(); + Construct_point_3 cp = tr.geom_traits().construct_point_3_object(); Construct_triangle_3 triangle = tr.geom_traits().construct_triangle_3_object(); - const Bare_point& p1 = wp2p(tr.point(f.first, (f.second+1)&3)); - const Bare_point& p2 = wp2p(tr.point(f.first, (f.second+2)&3)); - const Bare_point& p3 = wp2p(tr.point(f.first, (f.second+3)&3)); + const Weighted_point& wp1 = tr.point(f.first, (f.second+1)&3); + const Weighted_point& wp2 = tr.point(f.first, (f.second+2)&3); + const Weighted_point& wp3 = tr.point(f.first, (f.second+3)&3); + const Bare_point& p1 = cp(wp1); + const Bare_point& p2 = cp(wp2); + const Bare_point& p3 = cp(wp3); const FT triangle_area = area(triangle(p1,p2,p3)); const FT d12 = distance(p1,p2); @@ -351,12 +355,13 @@ protected: CGAL_assertion (f.first->is_facet_on_surface(f.second)); typedef typename Tr::Geom_traits Gt; - typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Weighted_point Weighted_point; - typename Gt::Construct_point_3 wp2p = - tr.geom_traits().construct_point_3_object(); + typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object(); - const Bare_point& p1 = wp2p(tr.point(f.first, (f.second+1)&3)); + const Weighted_point& wp1 = tr.point(f.first, (f.second+1)&3); + const Bare_point& p1 = cp(wp1); const Bare_point& ball_center = f.first->get_facet_surface_center(f.second); const Index& index = f.first->get_facet_surface_center_index(f.second); @@ -420,13 +425,14 @@ protected: CGAL_assertion (f.first->is_facet_on_surface(f.second)); CGAL_assertion (B_ != 0); - typedef typename Tr::Geom_traits Gt; - typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Geom_traits Gt; + typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Weighted_point Weighted_point; - typename Gt::Construct_point_3 wp2p = - tr.geom_traits().construct_point_3_object(); + typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object(); - const Bare_point& p1 = wp2p(tr.point(f.first, (f.second+1)&3)); + const Weighted_point& wp1 = tr.point(f.first, (f.second+1)&3); + const Bare_point p1 = cp(wp1); const Bare_point& ball_center = f.first->get_facet_surface_center(f.second); const FT sq_radius = tr.min_squared_distance(p1, ball_center); diff --git a/Mesh_3/include/CGAL/Mesh_3/vertex_perturbation.h b/Mesh_3/include/CGAL/Mesh_3/vertex_perturbation.h index a831e3f8874..93248cdb64f 100644 --- a/Mesh_3/include/CGAL/Mesh_3/vertex_perturbation.h +++ b/Mesh_3/include/CGAL/Mesh_3/vertex_perturbation.h @@ -123,16 +123,19 @@ namespace Mesh_3 { edge_sq_length(const typename Tr::Edge& e, const Tr& tr) { - typedef typename Tr::Geom_traits Gt; - typedef typename Tr::Bare_point Bare_point; - + typedef typename Tr::Geom_traits Gt; + typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Weighted_point Weighted_point; + typename Gt::Compute_squared_distance_3 sq_distance = tr.geom_traits().compute_squared_distance_3_object(); - typename Gt::Construct_point_3 wp2p = + typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object(); - const Bare_point& p = wp2p(tr.point(e.first, e.second)); - const Bare_point& q = wp2p(tr.point(e.first, e.third)); + const Weighted_point& wp = tr.point(e.first, e.second); + const Weighted_point& wq = tr.point(e.first, e.third); + const Bare_point& p = cp(wp); + const Bare_point& q = cp(wq); return sq_distance(p,q); } @@ -182,10 +185,11 @@ template class Abstract_perturbation { protected: - typedef typename C3T3::Vertex_handle Vertex_handle; - typedef typename C3T3::Cell_handle Cell_handle; - typedef typename C3T3::Triangulation::Geom_traits::FT FT; - + typedef typename C3T3::Vertex_handle Vertex_handle; + typedef typename C3T3::Cell_handle Cell_handle; + + typedef typename C3T3::Triangulation::Geom_traits::FT FT; + public: /** * @brief constructor @@ -393,12 +397,19 @@ class Gradient_based_perturbation : public Abstract_perturbation { protected: - typedef Abstract_perturbation Base; - typedef typename Base::Vertex_handle Vertex_handle; - typedef typename Base::Cell_handle Cell_handle; - typedef typename Base::FT FT; - typedef typename C3T3::Triangulation::Geom_traits::Vector_3 Vector_3; - + typedef Abstract_perturbation Base; + + typedef typename Base::Vertex_handle Vertex_handle; + typedef typename Base::Cell_handle Cell_handle; + + typedef typename C3T3::Triangulation Tr; + typedef typename Tr::Geom_traits Gt; + typedef typename Gt::FT FT; + typedef typename Gt::Vector_3 Vector_3; + + typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Weighted_point Weighted_point; + public: /** * @brief Constructor @@ -445,11 +456,6 @@ protected: std::vector& modified_vertices, bool *could_lock_zone = NULL) const { - typedef typename C3T3::Triangulation Tr; - typedef typename Tr::Geom_traits Gt; - typedef typename Gt::FT FT; - typedef typename Tr::Bare_point Bare_point; - typedef Triangulation_helpers Th; const Tr& tr = c3t3.triangulation(); @@ -458,9 +464,9 @@ protected: tr.geom_traits().compute_squared_length_3_object(); typename Gt::Construct_translated_point_3 translate = tr.geom_traits().construct_translated_point_3_object(); - typename Gt::Construct_point_3 wp2p = + typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object(); - typename Gt::Construct_weighted_point_3 p2wp = + typename Gt::Construct_weighted_point_3 cwp = tr.geom_traits().construct_weighted_point_3_object(); typename Gt::Construct_vector_3 vector = tr.geom_traits().construct_vector_3_object(); @@ -471,12 +477,12 @@ protected: modified_vertices.clear(); - // shenanigans with v->point() ? // norm depends on the local size of the mesh FT sq_norm = this->compute_perturbation_sq_amplitude(v, c3t3, sq_step_size_); FT step_length = CGAL::sqrt(sq_norm / sq_length(gradient_vector)); Vector_3 step_vector = step_length * gradient_vector; - Bare_point initial_loc = wp2p(v->point()); + const Weighted_point& weighted_initial_loc = c3t3.triangulation().point(v); + Bare_point initial_loc = cp(weighted_initial_loc); Bare_point new_loc = translate(initial_loc, step_vector); Bare_point final_loc = new_loc; @@ -489,7 +495,7 @@ protected: { // as long as no topological change takes place while(Th().no_topological_change__without_set_point(c3t3.triangulation(), - v, p2wp(final_loc)) + v, cwp(final_loc)) && (++i <= max_step_nb_) ) { new_loc = translate(new_loc, step_vector); @@ -504,7 +510,7 @@ protected: { Vector_3 move_vector = vector(initial_loc, final_loc); while( Th().no_topological_change(c3t3.triangulation(), v, - move_vector, p2wp(final_loc)) && + move_vector, cwp(final_loc)) && ++i <= max_step_nb_ ) { new_loc = translate(new_loc, step_vector); @@ -525,7 +531,7 @@ protected: // we know that there will be a combinatorial change return helper.update_mesh_topo_change(v, - p2wp(final_loc), + cwp(final_loc), criterion, std::back_inserter(modified_vertices), could_lock_zone); @@ -546,16 +552,21 @@ class Sq_radius_perturbation : public Gradient_based_perturbation { protected: - typedef Gradient_based_perturbation Base; - typedef typename Base::Vertex_handle Vertex_handle; - typedef typename Base::Cell_handle Cell_handle; - typedef typename Base::FT FT; - - typedef typename C3T3::Triangulation Tr; - typedef typename Tr::Geom_traits Gt; - typedef typename Gt::Vector_3 Vector_3; - typedef typename Tr::Bare_point Bare_point; - typedef typename Gt::Construct_point_3 Construct_point_3; + typedef Gradient_based_perturbation Base; + + typedef typename C3T3::Triangulation Tr; + typedef typename C3T3::Triangulation Triangulation; + + typedef typename Base::Vertex_handle Vertex_handle; + typedef typename Base::Cell_handle Cell_handle; + + typedef typename Tr::Geom_traits Gt; + typedef typename Gt::FT FT; + typedef typename Gt::Vector_3 Vector_3; + + typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Weighted_point Weighted_point; + public: /** * @brief Constructor @@ -649,22 +660,25 @@ private: const Cell_handle& cell, const Vertex_handle& v) const { - typedef typename C3T3::Triangulation Triangulation; - typedef typename Triangulation::Geom_traits Gt; - const Triangulation& tr = c3t3.triangulation(); typename Gt::Construct_translated_point_3 translate = c3t3.triangulation().geom_traits().construct_translated_point_3_object(); - typename Gt::Construct_point_3 wp2p = + typename Gt::Construct_point_3 cp = c3t3.triangulation().geom_traits().construct_point_3_object(); - // translate the tet so that cell->vertex((i+3)&3) is 0_{R^3} unsigned int index = cell->index(v); - Vector_3 translate_to_origin(CGAL::ORIGIN, wp2p(tr.point(cell, (index+3)&3))); //p4 - const Bare_point& p1 = translate(wp2p(tr.point(cell, index)), - translate_to_origin); - const Bare_point& p2 = translate(wp2p(tr.point(cell, (index+1)&3)), - translate_to_origin); - const Bare_point& p3 = translate(wp2p(tr.point(cell, (index+2)&3)), - translate_to_origin); + + const Weighted_point& wvp = tr.point(cell, index); + const Weighted_point& wp2 = tr.point(cell, (index+1)&3); + const Weighted_point& wp3 = tr.point(cell, (index+2)&3); + const Weighted_point& wp4 = tr.point(cell, (index+3)&3); + + // translate the tet so that 'wp4' is the origin + Vector_3 translate_to_origin(CGAL::ORIGIN, cp(wp4)); + const Bare_point& p1 = translate(cp(wvp), - translate_to_origin); + const Bare_point& p2 = translate(cp(wp2), - translate_to_origin); + const Bare_point& p3 = translate(cp(wp3), - translate_to_origin); // pre-compute everything FT sq_p1 = p1.x()*p1.x() + p1.y()*p1.y() + p1.z()*p1.z(); @@ -717,15 +731,19 @@ class Volume_perturbation : public Gradient_based_perturbation { protected: - typedef Gradient_based_perturbation Base; - typedef typename Base::Vertex_handle Vertex_handle; - typedef typename Base::Cell_handle Cell_handle; - typedef typename Base::FT FT; - - typedef typename C3T3::Triangulation::Geom_traits Gt; - typedef typename Gt::Vector_3 Vector_3; - typedef typename C3T3::Triangulation::Bare_point Bare_point; - + typedef Gradient_based_perturbation Base; + + typedef typename Base::Vertex_handle Vertex_handle; + typedef typename Base::Cell_handle Cell_handle; + + typedef typename C3T3::Triangulation Tr; + typedef typename Tr::Geom_traits Gt; + typedef typename Gt::FT FT; + typedef typename Gt::Vector_3 Vector_3; + + typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Weighted_point Weighted_point; + public: /** * @brief Constructor @@ -823,7 +841,7 @@ private: CGAL_assertion(cell->has_vertex(v)); const typename C3T3::Triangulation& tr = c3t3.triangulation(); - typename Gt::Construct_point_3 wp2p = + typename Gt::Construct_point_3 cp = c3t3.triangulation().geom_traits().construct_point_3_object(); const int i = cell->index(v); @@ -834,11 +852,14 @@ private: if ( (i&1) == 0 ) std::swap(k1,k3); - - const Bare_point& p1 = wp2p(tr.point(cell, k1)); - const Bare_point& p2 = wp2p(tr.point(cell, k2)); - const Bare_point& p3 = wp2p(tr.point(cell, k3)); - + + const Weighted_point& wp1 = tr.point(cell, k1); + const Weighted_point& wp2 = tr.point(cell, k2); + const Weighted_point& wp3 = tr.point(cell, k3); + const Bare_point& p1 = cp(wp1); + const Bare_point& p2 = cp(wp2); + const Bare_point& p3 = cp(wp3); + FT gx = p2.y()*p3.z() + p1.y()*(p2.z()-p3.z()) - p3.y()*p2.z() - p1.z()*(p2.y()-p3.y()); @@ -863,15 +884,18 @@ class Dihedral_angle_perturbation : public Gradient_based_perturbation { protected: - typedef Gradient_based_perturbation Base; - typedef typename Base::Vertex_handle Vertex_handle; - typedef typename Base::Cell_handle Cell_handle; - typedef typename Base::FT FT; - - typedef typename C3T3::Triangulation::Geom_traits Gt; - typedef typename Gt::Vector_3 Vector_3; - typedef typename C3T3::Triangulation::Bare_point Bare_point; - typedef typename Gt::Construct_point_3 Construct_point_3; + typedef Gradient_based_perturbation Base; + + typedef typename Base::Vertex_handle Vertex_handle; + typedef typename Base::Cell_handle Cell_handle; + + typedef typename C3T3::Triangulation Tr; + typedef typename Tr::Geom_traits Gt; + typedef typename Gt::FT FT; + typedef typename Gt::Vector_3 Vector_3; + + typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Weighted_point Weighted_point; public: /** @@ -970,13 +994,14 @@ private: CGAL_assertion(cell->has_vertex(v)); const typename C3T3::Triangulation& tr = c3t3.triangulation(); - typename Gt::Construct_point_3 wp2p = + typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object(); typename Gt::Compute_squared_distance_3 sq_distance = tr.geom_traits().compute_squared_distance_3_object(); const int i = cell->index(v); - const Bare_point& p0 = wp2p(tr.point(cell, i)); + const Weighted_point& wp0 = tr.point(cell, i); + const Bare_point& p0 = cp(wp0); // Other indices int k1 = (i+1)&3; @@ -995,9 +1020,12 @@ private: std::swap(k2,k3); // Here edge k1k2 minimizes dihedral angle - const Bare_point& p1 = wp2p(tr.point(cell, k1)); - const Bare_point& p2 = wp2p(tr.point(cell, k2)); - const Bare_point& p3 = wp2p(tr.point(cell, k3)); + const Weighted_point& wp1 = tr.point(cell, k1); + const Weighted_point& wp2 = tr.point(cell, k2); + const Weighted_point& wp3 = tr.point(cell, k3); + const Bare_point& p1 = cp(wp1); + const Bare_point& p2 = cp(wp2); + const Bare_point& p3 = cp(wp3); // grad of min dihedral angle (in cell) wrt p0 const Vector_3 p1p0 (p1,p0); @@ -1030,15 +1058,16 @@ private: typename Gt::Compute_approximate_dihedral_angle_3 approx_dihedral_angle = tr.geom_traits().compute_approximate_dihedral_angle_3_object(); - typename Gt::Construct_point_3 wp2p = + typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object(); - return CGAL::abs(approx_dihedral_angle(wp2p(p), - wp2p(tr.point(cell, k1)), - wp2p(tr.point(cell, k2)), - wp2p(tr.point(cell, k3)))); + const Weighted_point& wp1 = tr.point(cell, k1); + const Weighted_point& wp2 = tr.point(cell, k2); + const Weighted_point& wp3 = tr.point(cell, k3); + + return CGAL::abs(approx_dihedral_angle(cp(p), cp(wp1), cp(wp2), cp(wp3))); } - + /** * @brief returns the cotangent of \c value */ @@ -1055,7 +1084,7 @@ private: { const typename C3T3::Triangulation& tr = c3t3.triangulation(); - typename Gt::Construct_point_3 wp2p = tr.geom_traits().construct_point_3_object(); + typename Gt::Construct_point_3 cp = tr.geom_traits().construct_point_3_object(); int k1 = (i+1)&3; int k2 = (i+2)&3; @@ -1065,9 +1094,12 @@ private: if ( (i&1) == 1 ) std::swap(k1,k2); - const Bare_point& p1 = wp2p(tr.point(ch, k1)); - const Bare_point& p2 = wp2p(tr.point(ch, k2)); - const Bare_point& p3 = wp2p(tr.point(ch, k3)); + const Weighted_point& wp1 = tr.point(ch, k1); + const Weighted_point& wp2 = tr.point(ch, k2); + const Weighted_point& wp3 = tr.point(ch, k3); + const Bare_point& p1 = cp(wp1); + const Bare_point& p2 = cp(wp2); + const Bare_point& p3 = cp(wp3); // compute normal and return it return normal(p1, p2, p3); @@ -1087,12 +1119,18 @@ class Random_based_perturbation : public Abstract_perturbation { protected: - typedef Abstract_perturbation Base; - typedef typename Base::Vertex_handle Vertex_handle; - typedef typename Base::Cell_handle Cell_handle; - typedef typename Base::FT FT; - typedef typename C3T3::Triangulation::Geom_traits::Vector_3 Vector_3; - + typedef Abstract_perturbation Base; + + typedef typename Base::Vertex_handle Vertex_handle; + typedef typename Base::Cell_handle Cell_handle; + + typedef typename C3T3::Triangulation Tr; + typedef typename Tr::Geom_traits Gt; + typedef typename Gt::FT FT; + typedef typename Gt::Vector_3 Vector_3; + + typedef typename Tr::Bare_point Bare_point; + public: /** * @brief constructor @@ -1150,8 +1188,6 @@ protected: Vector_3 random_vector_fixed_size(const C3T3& c3t3, const FT& vector_sq_size) const { - typedef typename C3T3::Triangulation::Geom_traits Gt; - typename Gt::Compute_squared_length_3 sq_length = c3t3.triangulation().geom_traits().compute_squared_length_3_object(); @@ -1212,15 +1248,19 @@ class Li_random_perturbation : public Random_based_perturbation { protected: - typedef Random_based_perturbation Base; - typedef typename Base::Vertex_handle Vertex_handle; - typedef typename Base::Cell_handle Cell_handle; - typedef typename Base::FT FT; - - typedef typename C3T3::Triangulation::Geom_traits Gt; - typedef typename Gt::Vector_3 Vector_3; - typedef typename C3T3::Triangulation::Bare_point Bare_point; - + typedef Random_based_perturbation Base; + + typedef typename Base::Vertex_handle Vertex_handle; + typedef typename Base::Cell_handle Cell_handle; + + typedef typename C3T3::Triangulation Tr; + typedef typename Tr::Geom_traits Gt; + typedef typename Gt::FT FT; + typedef typename Gt::Vector_3 Vector_3; + + typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Weighted_point Weighted_point; + public: /** * @brief constructor @@ -1288,16 +1328,13 @@ private: { typedef Triangulation_helpers Th; - typedef typename C3T3::Triangulation Tr; - typedef typename Tr::Geom_traits Gt; - typename Gt::Construct_translated_point_3 translate = c3t3.triangulation().geom_traits().construct_translated_point_3_object(); - typename Gt::Construct_weighted_point_3 p2wp = + typename Gt::Construct_weighted_point_3 cwp = c3t3.triangulation().geom_traits().construct_weighted_point_3_object(); typename Gt::Construct_vector_3 vector = c3t3.triangulation().geom_traits().construct_vector_3_object(); - typename Gt::Construct_point_3 wp2p = + typename Gt::Construct_point_3 cp = c3t3.triangulation().geom_traits().construct_point_3_object(); modified_vertices.clear(); @@ -1308,8 +1345,9 @@ private: // norm depends on the local size of the mesh FT sq_norm = this->compute_perturbation_sq_amplitude(v, c3t3, this->sphere_sq_radius()); - const Bare_point initial_location = wp2p(v->point()); - + const Weighted_point& weighted_initial_location = c3t3.triangulation().point(v); + const Bare_point initial_location = cp(weighted_initial_location); + // Initialize loop variables bool criterion_improved = false; Vertex_handle moving_vertex = v; @@ -1338,7 +1376,7 @@ private: std::pair update = helper.update_mesh(moving_vertex, vector(initial_location, new_location), - p2wp(new_location), + cwp(new_location), criterion, std::back_inserter(tmp_mod_vertices), could_lock_zone); @@ -1352,7 +1390,7 @@ private: if ( update.first ) { criterion_improved = true; - best_location = wp2p(moving_vertex->point()); + best_location = cp(moving_vertex->point()); // this seems useless mod_vertices.insert(tmp_mod_vertices.begin(), tmp_mod_vertices.end()); diff --git a/Mesh_3/include/CGAL/Mesh_complex_3_in_triangulation_3.h b/Mesh_3/include/CGAL/Mesh_complex_3_in_triangulation_3.h index 515e49ac072..816eba41053 100644 --- a/Mesh_3/include/CGAL/Mesh_complex_3_in_triangulation_3.h +++ b/Mesh_3/include/CGAL/Mesh_complex_3_in_triangulation_3.h @@ -618,10 +618,10 @@ Mesh_complex_3_in_triangulation_3(const Self& rhs) const Vertex_handle& vb = it->left; Vertex_handle new_va; - this->triangulation().is_vertex(va->point(), new_va); + this->triangulation().is_vertex(rhs.triangulation().point(va), new_va); Vertex_handle new_vb; - this->triangulation().is_vertex(vb->point(), new_vb); + this->triangulation().is_vertex(rhs.triangulation().point(vb), new_vb); this->add_to_complex(make_internal_edge(new_va,new_vb), it->info); } @@ -631,7 +631,7 @@ Mesh_complex_3_in_triangulation_3(const Self& rhs) end = rhs.corners_.end() ; it != end ; ++it ) { Vertex_handle new_v; - this->triangulation().is_vertex(it->first->point(), new_v); + this->triangulation().is_vertex(rhs.triangulation().point(it->first), new_v); this->add_to_complex(new_v, it->second); } @@ -685,9 +685,10 @@ bool Mesh_complex_3_in_triangulation_3:: is_valid(bool verbose) const { - typedef typename Tr::Point::Point Bare_point; - typedef typename Tr::Point::Weight Weight; - typedef Weight FT; + typedef typename Tr::Bare_point Bare_point; + typedef typename Tr::Weighted_point Weighted_point; + typedef typename Tr::Weighted_point::Weight Weight; + typedef Weight FT; std::map vertex_map; @@ -712,7 +713,7 @@ is_valid(bool verbose) const { if(verbose) std::cerr << "Validity error: vertex " << (void*)(&*vit->first) - << " (" << vit->first->point() << ") " + << " (" << this->triangulation().point(vit->first) << ") " << "is not a corner (dimension " << vit->first->in_dimension() << ") but has " << vit->second << " neighbor(s)!\n"; return false; @@ -727,11 +728,13 @@ is_valid(bool verbose) const this->triangulation().geom_traits().construct_sphere_3_object(); typename Tr::Geom_traits::Do_intersect_3 do_intersect = this->triangulation().geom_traits().do_intersect_3_object(); - typename Tr::Geom_traits::Construct_point_3 wp2p = + typename Tr::Geom_traits::Construct_point_3 cp = this->triangulation().geom_traits().construct_point_3_object(); - const Bare_point& p = wp2p(it->right->point()); - const Bare_point& q = wp2p(it->left->point()); + const Weighted_point& itrwp = this->triangulation().point(it->right); + const Weighted_point& itlwp = this->triangulation().point(it->left); + const Bare_point& p = cp(itrwp); + const Bare_point& q = cp(itlwp); const FT& sq_rp = it->right->point().weight(); const FT& sq_rq = it->left->point().weight(); diff --git a/Mesh_3/include/CGAL/Polyhedral_complex_mesh_domain_3.h b/Mesh_3/include/CGAL/Polyhedral_complex_mesh_domain_3.h index b78ff0a5830..1925e99f0d3 100644 --- a/Mesh_3/include/CGAL/Polyhedral_complex_mesh_domain_3.h +++ b/Mesh_3/include/CGAL/Polyhedral_complex_mesh_domain_3.h @@ -512,9 +512,13 @@ public: std::cout.flush(); #endif CGAL::Random random(0); - typedef typename C3t3::Triangulation Tr; + + typedef typename C3t3::Triangulation Tr; + typedef typename IGT::Sphere_3 Sphere_3; + typedef typename Polyhedron::Vertex_const_handle Vertex_const_handle; + Tr& tr = c3t3.triangulation(); - typedef typename Polyhedron::Vertex_const_handle Vertex_const_handle; + typename Tr::Geom_traits::Construct_weighted_point_3 cwp = tr.geom_traits().construct_weighted_point_3_object(); @@ -538,10 +542,11 @@ public: const Patch_id patch_id = vit->halfedge()->face()->patch_id(); CGAL_assertion(std::size_t(patch_id) <= nb_of_patch_plus_one); typename Tr::Vertex_handle tr_v = tr.nearest_power_vertex(vit->point()); - if (tr_v != typename Tr::Vertex_handle()) { - typedef typename IGT::Sphere_3 Sphere_3; + if (tr_v != typename Tr::Vertex_handle()) + { const Sphere_3 sphere(tr_v->point().point(), tr_v->point().weight()); - if (!sphere.has_on_unbounded_side(vit->point())) continue; + if (!sphere.has_on_unbounded_side(vit->point())) + continue; } ++nb_of_free_vertices_on_patch[patch_id]; } diff --git a/Mesh_3/include/CGAL/refine_mesh_3.h b/Mesh_3/include/CGAL/refine_mesh_3.h index 38a4f2da50c..22b3a715926 100644 --- a/Mesh_3/include/CGAL/refine_mesh_3.h +++ b/Mesh_3/include/CGAL/refine_mesh_3.h @@ -64,7 +64,7 @@ namespace CGAL { void operator()(const Vertex& vertex) const { - typename Geom_traits::Construct_point_3 wp2p = + typename Geom_traits::Construct_point_3 cp = r_c3t3_.triangulation().geom_traits().construct_point_3_object(); typename Geom_traits::Compute_weight_3 cw = r_c3t3_.triangulation().geom_traits().compute_weight_3_object(); @@ -72,7 +72,7 @@ namespace CGAL { // Get vh properties int dimension = vertex.in_dimension(); Weight w = (dimension < 2) ? cw(vertex.point()) : 0; - Weighted_point point(wp2p(vertex.point()), w); + Weighted_point point(cp(vertex.point()), w); Index index = vertex.index(); // Insert point and restore handle properties diff --git a/Periodic_3_mesh_3/examples/Periodic_3_mesh_3/mesh_implicit_shape.cpp b/Periodic_3_mesh_3/examples/Periodic_3_mesh_3/mesh_implicit_shape.cpp index 5b5528b6113..3fe1a45cb20 100644 --- a/Periodic_3_mesh_3/examples/Periodic_3_mesh_3/mesh_implicit_shape.cpp +++ b/Periodic_3_mesh_3/examples/Periodic_3_mesh_3/mesh_implicit_shape.cpp @@ -1,4 +1,8 @@ #define CGAL_MESH_3_VERBOSE +#define _DEBUG +#define CGAL_MESHES_DEBUG_REFINEMENT_POINTS +#define CGAL_MESH_3_DEBUG_FACET_CRITERIA +#define CGAL_MESH_3_DEBUG_CELL_CRITERIA #include #include @@ -72,8 +76,8 @@ int main(int argc, char** argv) C3t3 c3t3 = CGAL::make_periodic_3_mesh_3(domain, criteria, // odt(), no_odt(), -// lloyd(), - no_lloyd(), + lloyd(), +// no_lloyd(), // exude(), no_exude(), // perturb(sliver_bound=10, time_limit=10) diff --git a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h index 3162ad102a9..52c07f56777 100644 --- a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h +++ b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h @@ -245,6 +245,11 @@ public: Bare_point canonicalize_point(const Bare_point& p) const { + if(p.x() >= domain().xmin() && p.x() < domain().xmax() && + p.y() >= domain().ymin() && p.y() < domain().ymax() && + p.z() >= domain().zmin() && p.z() < domain().zmax()) + return p; + return robust_canonicalize_point(p); } @@ -267,6 +272,11 @@ public: Weighted_point canonicalize_point(const Weighted_point& p) const { + if(p.x() >= domain().xmin() && p.x() < domain().xmax() && + p.y() >= domain().ymin() && p.y() < domain().ymax() && + p.z() >= domain().zmin() && p.z() < domain().zmax()) + return p; + return robust_canonicalize_point(p); } @@ -381,8 +391,13 @@ public: // between p and q, for all combinations of offsets FT min_squared_distance(const Bare_point& p, const Bare_point& q) const { +// std::cout << "minsqd: " << p << " // " << q << std::endl; + const Bare_point cp = canonicalize_point(p); const Bare_point cq = canonicalize_point(q); + +// std::cout << "canon: " << cp << " // " << cq << std::endl; + FT min_sq_dist = std::numeric_limits::infinity(); for(int i = 0; i < 3; ++i) { @@ -397,6 +412,7 @@ public: } } +// std::cout << "minsqdt: " << min_sq_dist << std::endl; return min_sq_dist; }