diff --git a/Operations_on_polyhedra/dont_submit b/Operations_on_polyhedra/dont_submit index a49429c0eba..f2c7eb7a35b 100644 --- a/Operations_on_polyhedra/dont_submit +++ b/Operations_on_polyhedra/dont_submit @@ -1,2 +1,3 @@ examples include/CGAL/Point_inside_polyhedron_3_old.h +test diff --git a/Operations_on_polyhedra/include/CGAL/intersection_of_Polyhedra_3_refinement_visitor.h b/Operations_on_polyhedra/include/CGAL/intersection_of_Polyhedra_3_refinement_visitor.h index e670794ff53..e9383a3cc87 100644 --- a/Operations_on_polyhedra/include/CGAL/intersection_of_Polyhedra_3_refinement_visitor.h +++ b/Operations_on_polyhedra/include/CGAL/intersection_of_Polyhedra_3_refinement_visitor.h @@ -30,7 +30,7 @@ #include -#include +#include #include #include #include @@ -774,7 +774,7 @@ class Node_visitor_refine_polyhedra{ } } } - #ifndef NDEBUG + #ifdef CGAL_COREFINEMENT_DEBUG else { std::cout << "X0: Found an isolated point" << std::endl; @@ -2137,17 +2137,19 @@ public: //update the info of each volume knowing about only one polyhedron: //this happens when one polyhedron has a connected component //that do not intersect the other polyhedron - typedef CGAL::Polyhedral_mesh_domain_3 Mesh_domain; + + typedef Point_inside_polyhedron_3 Inside_poly_test; + CGAL_precondition(polyhedron_to_map_node_to_polyhedron_vertex.size()==2); Polyhedron* Poly_A = polyhedron_to_map_node_to_polyhedron_vertex.begin()->first; Polyhedron* Poly_B = boost::next(polyhedron_to_map_node_to_polyhedron_vertex.begin())->first; - Mesh_domain* domain_A_ptr=NULL; - Mesh_domain* domain_B_ptr=NULL; - + Inside_poly_test* inside_A_test_ptr=NULL; + Inside_poly_test* inside_B_test_ptr=NULL; + #ifdef CGAL_COREFINEMENT_DEBUG final_map().display_characteristics(std::cout); std::cout << "\n"; #endif - + typename Combinatorial_map_3::template One_dart_per_cell_range<3> cell_range=final_map().template one_dart_per_cell<3>(); for (typename Combinatorial_map_3::template One_dart_per_cell_range<3>::iterator it = cell_range.begin(), it_end=cell_range.end(); @@ -2157,27 +2159,72 @@ public: internal_IOP::Volume_info& info=it->template attribute<3>()->info(); std::size_t inside_size=info.inside.size(); std::size_t outside_size=info.outside.size(); + + // if a volume is not classified wrt the two polyhedra, it means the component we look at does not + // is a disjoint (but maybe at a vertex TAG SL001) if ( inside_size + outside_size == 1) { bool is_inside = (inside_size==1); Polyhedron* current_poly= is_inside? (*info.inside.begin()):(*info.outside.begin()); - Polyhedron* test_poly; - Mesh_domain* domain_ptr; + Polyhedron* test_poly; + Inside_poly_test* inside_test_ptr; if ( current_poly==Poly_A) { test_poly=Poly_B; - if (domain_B_ptr == NULL) domain_B_ptr=new Mesh_domain(*Poly_B); - domain_ptr=domain_B_ptr; + if (inside_B_test_ptr == NULL) inside_B_test_ptr=new Inside_poly_test(*Poly_B); + inside_test_ptr=inside_B_test_ptr; } else { test_poly=Poly_A; - if (domain_A_ptr == NULL) domain_A_ptr=new Mesh_domain(*Poly_A); - domain_ptr=domain_A_ptr; + if (inside_A_test_ptr == NULL) inside_A_test_ptr=new Inside_poly_test(*Poly_A); + inside_test_ptr=inside_A_test_ptr; } - typename Mesh_domain::Is_in_domain is_in_domain(*domain_ptr); + + // We need to find a point of the volume that is not on the boundary of the other volume. + // Then the position of this point give the position of the volume. If all the points are on + // the bounday, we take the mid-point of an edge (which must not be on the boundary otherwise + // it contradicts the fact that volumes are disjoint + // We first use the dart we have since one_dart_per_incident_cell has a non-negligeable cost. typename Kernel::Point_3 query=it->template attribute<0>()->point(); - if ( is_in_domain(query) ) + CGAL::Bounded_side res = (*inside_test_ptr)(query); + if (res==ON_BOUNDARY) + { + typedef typename Combinatorial_map_3:: + template One_dart_per_incident_cell_range<0,3> Vertex_range; + + Vertex_range vertex_range = + final_map().template one_dart_per_incident_cell<0,3>(it); + typename Vertex_range::iterator vit = vertex_range.begin(); + + CGAL_assertion( typename Combinatorial_map_3::Dart_handle(vit) == + typename Combinatorial_map_3::Dart_handle(it) ); + ++vit; + for ( ; vit!=vertex_range.end(); ++vit) + { + query=vit->template attribute<0>()->point(); + res = (*inside_test_ptr)(query); + if ( res != ON_BOUNDARY ) break; + } + + //take edge midpoint + if (res == ON_BOUNDARY) + { + /// \todo see do a better test here. At least the construction cannot fail + /// but the mid-point can fall outside of the volume... + #ifdef CGAL_COREFINEMENT_DEBUG + #warning this is not exact!!! + #endif + typename Kernel::Point_3 p1=it->template attribute<0>()->point(); + typename Kernel::Point_3 p2=it->template beta<1>()->template attribute<0>()->point(); + query=midpoint(p1,p2); + res = (*inside_test_ptr)(query); + } + + CGAL_assertion( res!= ON_BOUNDARY ); + } + + if (res == ON_BOUNDED_SIDE ) info.inside.insert(test_poly); else info.outside.insert(test_poly); @@ -2193,9 +2240,8 @@ public: std::cout << std::endl; #endif } - if (domain_A_ptr!=NULL) delete domain_A_ptr; - if (domain_B_ptr!=NULL) delete domain_B_ptr; - + if (inside_A_test_ptr!=NULL) delete inside_A_test_ptr; + if (inside_B_test_ptr!=NULL) delete inside_B_test_ptr; } }; diff --git a/Operations_on_polyhedra/test/Operations_on_polyhedra/corefinement_point_tangency.cpp b/Operations_on_polyhedra/test/Operations_on_polyhedra/corefinement_point_tangency.cpp new file mode 100644 index 00000000000..8cf726705dd --- /dev/null +++ b/Operations_on_polyhedra/test/Operations_on_polyhedra/corefinement_point_tangency.cpp @@ -0,0 +1,112 @@ +#include +#include +#include +#include + +#include +#include + +const char* cube_a = +"OFF 8 12 0\n\ +0 0 0\n\ +0 1 0\n\ +1 1 0\n\ +1 0 0\n\ +1 0 1\n\ +1 1 1\n\ +0 1 1\n\ +0 0 1\n\ +3 0 1 2\n\ +3 3 0 2\n\ +3 4 2 5\n\ +3 4 3 2\n\ +3 1 6 5\n\ +3 2 1 5\n\ +3 0 6 1\n\ +3 0 7 6\n\ +3 4 5 6\n\ +3 7 4 6\n\ +3 3 4 7\n\ +3 0 3 7\n"; + +const char* cube_b= +"OFF\n\ +8 12 0\n\ +-1 -1 -1\n\ +-1 0 -1\n\ +0 0 -1\n\ +0 -1 -1\n\ +0 -1 0\n\ +0 0 0\n\ +-1 0 0\n\ +-1 -1 0\n\ +3 0 1 2\n\ +3 3 0 2\n\ +3 4 2 5\n\ +3 4 3 2\n\ +3 1 6 5\n\ +3 2 1 5\n\ +3 0 6 1\n\ +3 0 7 6\n\ +3 4 5 6\n\ +3 7 4 6\n\ +3 3 4 7\n\ +3 0 3 7\n"; + +const char* inside_a= +"OFF 4 4 0\n\ +1 0.5 0.5\n\ +0.5 1 0.5\n\ +0.5 0 0.5\n\ +0.5 0.5 1\n\ +3 0 2 1\n\ +3 0 1 3\n\ +3 0 3 2\n\ +3 1 2 3\n"; + +typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel ; + +int main() +{ + typedef CGAL::Polyhedron_3 Polyhedron; + typedef CGAL::Polyhedron_corefinement Corefinement; + + std::stringstream fpa(cube_a); + std::stringstream fpb(cube_b); + Polyhedron pa; + Polyhedron pb; + + fpa >> pa; + fpb >> pb; + + assert( pa.size_of_vertices() == 8); + assert( pb.size_of_vertices() == 8); + + { + Corefinement coref; + std::list > polylines; + std::vector > result; + coref( pa, pb, std::back_inserter(polylines), std::back_inserter(result), Corefinement::Intersection_tag ); + + assert( polylines.size() == 1 ); + assert( polylines.begin()->size() == 1 ); + assert( result.size() == 0 ); + } + + pb.clear(); + std::stringstream fpc(inside_a); + fpc >> pb; + assert( pb.size_of_vertices() == 4); + + { + Corefinement coref; + std::list > polylines; + std::vector > result; + coref( pa, pb, std::back_inserter(polylines), std::back_inserter(result), Corefinement::Intersection_tag ); + + assert( polylines.size() == 4 ); + assert( polylines.begin()->size() == 1 ); + assert( result.size() == 1 ); + assert( result[0].first->size_of_vertices() == 4); + } +}