diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h index 615ab82d7bd..af10e18962e 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/Visitor.h @@ -43,6 +43,7 @@ template struct Default_node_visitor{ typedef boost::graph_traits GT; typedef typename GT::halfedge_descriptor halfedge_descriptor; + typedef typename GT::face_descriptor face_descriptor; typedef typename GT::vertex_descriptor vertex_descriptor; void new_node_added( std::size_t /* node_id */, @@ -53,6 +54,13 @@ struct Default_node_visitor{ bool /* is_source_coplanar */ ) {} + void new_node_added_triple_face(std::size_t /* node_id */, + face_descriptor /* f1 */, + face_descriptor /* f2 */, + face_descriptor /* f3 */, + const TriangleMesh& /* tm */) + {} + void new_vertex_added(std::size_t /* node_id */, vertex_descriptor /* vh */, TriangleMesh& /*tm*/){} @@ -333,6 +341,22 @@ public: mesh_to_vertex_to_node_id[&tm].insert(std::make_pair(target(h,tm),node_id)); } + void new_node_added_triple_face(std::size_t node_id, + face_descriptor f1, + face_descriptor f2, + face_descriptor f3, + const TriangleMesh& tm) // TODO check if we need a special case if the endpoint of the intersect edge is on the third face + { + CGAL_assertion(f1!=f2 && f1!=f3 && f2!=f3); + TriangleMesh* tm_ptr = const_cast(&tm); + new_node_visitor.new_node_added_triple_face(node_id, f1, f2, f3, tm); +#ifdef CGAL_DEBUG_AUTOREFINEMENT + std::cout << "adding node " << node_id << " " << f1 << " " << f2 << " " << f3 << "\n"; +#endif + on_face[tm_ptr][f1].push_back(node_id); + on_face[tm_ptr][f2].push_back(node_id); + on_face[tm_ptr][f3].push_back(node_id); + } void new_node_added(std::size_t node_id, Intersection_type type, @@ -565,7 +589,8 @@ public: // this condition ensures to consider only graph edges that are in // the same triangle if ( !points_on_triangle || it_vh!=id_to_CDT_vh.end() ){ - CGAL_assertion(it_vh!=id_to_CDT_vh.end()); + CGAL_assertion(doing_autorefinement || it_vh!=id_to_CDT_vh.end()); + if (it_vh==id_to_CDT_vh.end()) continue; // needed for autorefinement (interior nodes) cdt.insert_constraint(vh,it_vh->second); constrained_edges.push_back(std::make_pair(id,id_n)); constrained_edges.push_back(std::make_pair(id_n,id)); diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_impl.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_impl.h index 5ecb21742aa..430b96b87b2 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_impl.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_impl.h @@ -77,10 +77,11 @@ struct Self_intersection_exception{}; // // Coplanar triangles are filtered out and handled separately. // -template +template struct Default_surface_intersection_visitor{ typedef boost::graph_traits Graph_traits; typedef typename Graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename Graph_traits::face_descriptor face_descriptor; void new_node_added( std::size_t,Intersection_type,halfedge_descriptor,halfedge_descriptor, @@ -97,7 +98,15 @@ struct Default_surface_intersection_visitor{ const TriangleMesh&, const TriangleMesh&, const VertexPointMap&, const VertexPointMap&) {} - static const bool Predicates_on_constructions_needed = false; + void new_node_added_triple_face(std::size_t /* node_id */, + face_descriptor /* f1 */, + face_descriptor /* f2 */, + face_descriptor /* f3 */, + const TriangleMesh& /* tm */) + {} + // this is required in autorefinement for the do-intersect of 3 segments. + // If we implement a predicate only test, we can get rid of it. + static const bool Predicates_on_constructions_needed = doing_autorefinement; static const bool do_need_vertex_graph = false; }; @@ -172,19 +181,19 @@ class Intersection_of_triangle_meshes // we use Face_pair_and_int and not Face_pair to handle coplanar case. // Indeed the boundary of the intersection of two coplanar triangles // may contain several segments. - typedef std::map< Face_pair_and_int, Node_id_set > Faces_to_nodes_map; - + typedef std::map< Face_pair_and_int, Node_id_set > Faces_to_nodes_map; + typedef Intersection_nodes Node_vector; // data members Edge_to_faces stm_edge_to_ltm_faces; // map edges from the triangle mesh with the smaller address to faces of the triangle mesh with the larger address Edge_to_faces ltm_edge_to_stm_faces; // map edges from the triangle mesh with the larger address to faces of the triangle mesh with the smaller address // here face descriptor are from tmi and tmj such that &tmi<&tmj Coplanar_face_set coplanar_faces; - Intersection_nodes nodes; + Node_vector nodes; Node_visitor visitor; - Faces_to_nodes_map f_to_node; //Associate a pair of triangle to their intersection points + Faces_to_nodes_map f_to_node; //Associate a pair of triangles to their intersection points CGAL_assertion_code(bool doing_autorefinement;) // member functions void filter_intersections(const TriangleMesh& tm_f, @@ -665,7 +674,7 @@ class Intersection_of_triangle_meshes Inter_type res=intersection_type(h_1,f_2,tm1,tm2,vpm1,vpm2); Intersection_type type=cpp11::get<0>(res); - //handle degenerate case: one extremity of edge belomg to f_2 + //handle degenerate case: one extremity of edge belong to f_2 std::vector all_edges; if ( cpp11::get<3>(res) ) // is edge target in triangle plane std::copy(halfedges_around_target(h_1,tm1).first, @@ -791,6 +800,140 @@ class Intersection_of_triangle_meshes } }; + /// TODO replace this by a lexical sort + struct Less_for_nodes_along_an_edge{ + Node_vector& nodes; + Node_id ref; + Less_for_nodes_along_an_edge(Node_vector& nodes, Node_id ref) + : nodes(nodes), ref(ref) + {} + bool operator()(Node_id i, Node_id j) const + { + return + compare_distance_to_point(nodes.exact_node(ref), + nodes.exact_node(i), + nodes.exact_node(j) ) == SMALLER; + } + }; + + void detect_intersections_in_the_graph(const TriangleMesh& tm, + const VertexPointMap& vpm, + Node_id& current_node) + { + boost::unordered_map > face_intersections; + for (typename Faces_to_nodes_map::iterator it=f_to_node.begin(); + it!=f_to_node.end(); + ++it){ + + face_descriptor f1 = it->first.first.first, f2 = it->first.first.second; + face_intersections[f1].push_back(f2); + face_intersections[f2].push_back(f1); + } + + std::map, std::vector > map_to_process; // TODO find a better way to avoid too many queries. + + typedef std::pair > Pair_type; + BOOST_FOREACH(Pair_type& p, face_intersections) + { + std::size_t nb_faces = p.second.size(); + // TODO handle 4 and more faces intersecting (only 3 right now) + for(std::size_t i=0; i > triple_intersections; + for(std::size_t j=i+1; jsecond[0]; + Node_id s12=f_to_node.find(std::make_pair(std::make_pair(p.first, p.second[i]), 0))->second[1]; + Node_id s21=f_to_node.find(std::make_pair(std::make_pair(p.second[i], p.second[j]), 0))->second[0]; + Node_id s22=f_to_node.find(std::make_pair(std::make_pair(p.second[i], p.second[j]), 0))->second[1]; + + if( do_intersect( + typename Node_vector::Exact_kernel::Segment_3(nodes.exact_node(s11), nodes.exact_node(s12)), + typename Node_vector::Exact_kernel::Segment_3(nodes.exact_node(s21), nodes.exact_node(s22)) + ) + ) + triple_intersections[p.second[i]].push_back(p.second[j]); + } + } + if (!triple_intersections.empty()) + { + ///TODO throw an exception when more than 2 triangles intersect along an edge + typedef std::pair > PT; + BOOST_FOREACH(PT pt, triple_intersections) + { + BOOST_FOREACH(face_descriptor fd, pt.second) + { + nodes.add_new_node(halfedge(p.first, tm), + halfedge(p.second[i], tm), + halfedge(fd, tm), + tm, vpm); + Node_id node_id=++current_node; +#ifdef CGAL_DEBUG_AUTOREFINEMENT + std::cerr << "New triple node " << node_id << "\n"; +#endif + visitor.new_node_added_triple_face(node_id,p.first, p.second[i], fd, tm); + + map_to_process[std::make_pair(p.first, p.second[i])].push_back(node_id); + map_to_process[std::make_pair(p.first, fd)].push_back(node_id); + map_to_process[std::make_pair(p.second[i], fd)].push_back(node_id); + } + } + } + } + } + +#ifdef CGAL_DEBUG_AUTOREFINEMENT + std::cout << "\nAt the end of new node creation, current_node is " << current_node << "\n\n"; +#endif + typedef std::pair, + std::vector > Faces_and_nodes; + BOOST_FOREACH(Faces_and_nodes& f_n_nids, map_to_process) + { + //get the original entry and remove it + typename Faces_to_nodes_map::iterator find_it = + f_to_node.find( std::make_pair(f_n_nids.first, 0) ); + CGAL_assertion(find_it!=f_to_node.end()); + CGAL_assertion(find_it->second.size()==2); // TODO handle case size = 1 + Node_id n1 = find_it->second[0]; + Node_id n2 = find_it->second[1]; + + // sort node ids along the edge + std::sort(f_n_nids.second.begin(), + f_n_nids.second.end(), + Less_for_nodes_along_an_edge(nodes, n1)); + + // insert new segments + Node_id prev = n1; + f_n_nids.second.push_back(n2); + int i=0; + typename Faces_to_nodes_map::iterator insert_it = find_it; +#ifdef CGAL_DEBUG_AUTOREFINEMENT + std::cout << n1 << " -> " << n2 << "\n"; +#endif + BOOST_FOREACH(Node_id id, f_n_nids.second) + { + insert_it = f_to_node.insert(insert_it, std::make_pair( + std::make_pair(f_n_nids.first,--i), Node_id_set()) ); // I have picked negative int for refined edges + insert_it->second.insert(prev); + insert_it->second.insert(id); + CGAL_assertion(insert_it->second.size()==2); +#ifdef CGAL_DEBUG_AUTOREFINEMENT + std::cerr <<" adding " << prev << " " << id << " into " + << f_n_nids.first.first << " and " << f_n_nids.first.second << "\n"; +#endif + prev=id; + } + f_to_node.erase(find_it); + } + } template void construct_polylines(Output_iterator out){ @@ -1083,6 +1226,8 @@ public: // we can remove one intersecting edge out of the two remove_duplicated_intersecting_edges(); + detect_intersections_in_the_graph(tm, vpm, current_node); + // If a pair of faces defines an isolated node, check if they share a common // vertex and create a new node in that case. add_common_vertices_for_pairs_of_faces_with_isolated_node(current_node); diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_nodes.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_nodes.h index 47b3ac1ef71..aee55abf50b 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_nodes.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Corefinement/intersection_nodes.h @@ -234,6 +234,35 @@ public: to_exact( get(vpm_a, target(h_a,tm_a)) ) ) ); } + // use to resolve intersection of 3 faces in autorefinement only + void add_new_node(halfedge_descriptor h1, + halfedge_descriptor h2, + halfedge_descriptor h3, + const TriangleMesh& tm, + const VertexPointMap& vpm) + { + // TODO Far from optimal! + typedef Exact_kernel::Plane_3 Plane_3; + Plane_3 p1(to_exact( get(vpm, source(h1,tm)) ), + to_exact( get(vpm, target(h1,tm)) ), + to_exact( get(vpm, target(next(h1,tm),tm)))), + p2(to_exact( get(vpm, source(h2,tm)) ), + to_exact( get(vpm, target(h2,tm)) ), + to_exact( get(vpm, target(next(h2,tm),tm)))), + p3(to_exact( get(vpm, source(h3,tm)) ), + to_exact( get(vpm, target(h3,tm)) ), + to_exact( get(vpm, target(next(h3,tm),tm)))); + typename cpp11::result_of< + Exact_kernel::Intersect_3(Plane_3, Plane_3, Plane_3) + >::type inter_res = intersection(p1, p2, p3); + + CGAL_assertion(inter_res != boost::none); + const Exact_kernel::Point_3* pt = + boost::get(&(*inter_res)); + CGAL_assertion(pt!=NULL); + add_new_node(*pt); + } + void add_new_node(halfedge_descriptor edge_1, face_descriptor face_2) { add_new_node(edge_1, face_2, tm1, tm2, vpm1, vpm2); diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/intersection.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/intersection.h index 9e43ceffc82..b4d97a11e95 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/intersection.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/intersection.h @@ -1730,7 +1730,9 @@ surface_self_intersection(const TriangleMesh& tm, vertex_point); // surface intersection algorithm call - Corefinement::Intersection_of_triangle_meshes + typedef Corefinement::Default_surface_intersection_visitor Visitor; + Corefinement::Intersection_of_triangle_meshes functor(tm, vpm); polyline_output=functor(polyline_output, true);