diff --git a/Triangulation_2/include/CGAL/Triangulation_2/internal/Triangulation_2_projection_traits_base_3.h b/Triangulation_2/include/CGAL/Triangulation_2/internal/Triangulation_2_projection_traits_base_3.h index 9a160041dda..9a728577a4b 100644 --- a/Triangulation_2/include/CGAL/Triangulation_2/internal/Triangulation_2_projection_traits_base_3.h +++ b/Triangulation_2/include/CGAL/Triangulation_2/internal/Triangulation_2_projection_traits_base_3.h @@ -190,6 +190,7 @@ public: { CGAL_PROFILER("Projected_intersect_3::operator()") CGAL_TIME_PROFILER("Projected_intersect_3::operator()") + typedef boost::variant variant_type; const Vector_3 u1 = cross_product(s1.to_vector(), normal); if(u1 == NULL_VECTOR) return K().intersect_3_object()(s1.supporting_line(), s2); @@ -203,21 +204,29 @@ public: auto planes_intersection = intersection(plane_1, plane_2); if(! planes_intersection) { +#ifdef CGAL_T2_PTB_3_DEBUG std::cerr << "planes_intersection is empty\n"; +#endif return boost::none; } if(const Line* line = boost::get(&*planes_intersection)) { + // check if the intersection line intersects both segments by + // checking if a point on the intersection line is between + // the segments endpoints const Point& pi = line->point(0); - if(cross_product(normal, pi - s1.source()) + if( cross_product(normal, pi - s1.source()) * cross_product(normal, pi - s1.target()) > FT(0) || - cross_product(normal, pi - s2.source()) + cross_product(normal, pi - s2.source()) * cross_product(normal, pi - s2.target()) > FT(0) ) { - // the intersection of the lines is not inside the segments + // the intersection of the supporting lines is not inside both segments +#ifdef CGAL_T2_PTB_3_DEBUG std::cerr << "intersection not inside\n"; return boost::none; +#endif + return boost::none; } else { @@ -230,16 +239,67 @@ public: if(! inter){ return boost::none; } - if(const Point* point = boost::get(&*inter)){ - typedef boost::variant variant_type; + if(const Point* point = boost::get(&*inter)){ return boost::make_optional(variant_type(*point)); } } } if(boost::get(&*planes_intersection)) { - std::cerr << "coplanar lines\n"; - CGAL_error(); +#ifdef CGAL_T2_PTB_3_DEBUG + std::cerr << "coplanar supporting lines\n"; +#endif + auto is_inside_segment = [](const Segment& s, const Point& q) + { + return Vector_3(q,s.source()) * Vector_3(q,s.target()) <=0; + }; + + bool src2_in_s1 = is_inside_segment(s1, s2.source()); + bool tgt2_in_s1 = is_inside_segment(s1, s2.target()); + bool src1_in_s2 = is_inside_segment(s2, s1.source()); + bool tgt1_in_s2 = is_inside_segment(s2, s1.target()); + + if (src1_in_s2 && tgt1_in_s2) return boost::make_optional(variant_type(s1)); + if (src2_in_s1 && tgt2_in_s1) return boost::make_optional(variant_type(s2)); + + if (src1_in_s2) + { + if (src2_in_s1) + { + if (cross_product(normal, Vector_3(s1.source(), s2.source())) != NULL_VECTOR) + return boost::make_optional(variant_type(Segment(s1.source(), s2.source()))); + else + return boost::make_optional(variant_type((s1.source()))); + } + if (tgt2_in_s1) + { + if (cross_product(normal, Vector_3(s1.source(), s2.target())) != NULL_VECTOR) + return boost::make_optional(variant_type(Segment(s1.source(), s2.target()))); + else + return boost::make_optional(variant_type(s1.source())); + } + // should never get here with a Kernel with exact constructions + return boost::make_optional(variant_type(s1.source())); + } + if (tgt1_in_s2) + { + if (src2_in_s1) + { + if (cross_product(normal, Vector_3(s1.target(), s2.source())) != NULL_VECTOR) + return boost::make_optional(variant_type(Segment(s1.target(), s2.source()))); + else + return boost::make_optional(variant_type(s1.target())); + } + if (tgt2_in_s1) + { + if (cross_product(normal, Vector_3(s1.target(), s2.target())) != NULL_VECTOR) + return boost::make_optional(variant_type(Segment(s1.target(), s2.target()))); + else + return boost::make_optional(variant_type(s1.target())); + } + // should never get here with a Kernel with exact constructions + return boost::make_optional(variant_type(s1.target())); + } return boost::none; } return boost::none; diff --git a/Triangulation_2/test/Triangulation_2/test_cdt_2_projection_traits_special_case.cpp b/Triangulation_2/test/Triangulation_2/test_cdt_2_projection_traits_special_case.cpp index 26693b95a28..bf163771380 100644 --- a/Triangulation_2/test/Triangulation_2/test_cdt_2_projection_traits_special_case.cpp +++ b/Triangulation_2/test/Triangulation_2/test_cdt_2_projection_traits_special_case.cpp @@ -60,6 +60,67 @@ bool test(std::string test_name) return cdt.is_valid(); } + +template +bool test_segment_intersections(std::string test_name) +{ + std::cerr << "Testing " << test_name << std::endl; + + CDT_2_traits traits(typename K::Vector_3(0,0,1)); + typename CDT_2_traits::Intersect_2 intersect = traits.intersect_2_object(); + + typedef typename CDT_2_traits::Segment_2 Segment_2; + typedef typename CDT_2_traits::Point_2 Point_2; + + // test point intersection + Segment_2 s1( Point_2(0,0,0), Point_2(1,1,0) ), + s2( Point_2(0,1,0), Point_2(1,0,0) ); + + CGAL::Object o = intersect(s1,s2); + assert( o.is() ); + +// test segment overlap + Point_2 pts[4] = { Point_2(0,0,0), Point_2(1,0,1), Point_2(2,0,2), Point_2(4,0,3) }; + + o = intersect( Segment_2(pts[0], pts[1]), Segment_2(pts[2], pts[3]) ); + assert( o.empty() ); + + // pure overlap + o = intersect( Segment_2(pts[0], pts[2]), Segment_2(pts[1], pts[3]) ); + assert( o.is() ); + o = intersect( Segment_2(pts[0], pts[2]), Segment_2(pts[3], pts[1]) ); + assert( o.is() ); + o = intersect( Segment_2(pts[2], pts[0]), Segment_2(pts[1], pts[3]) ); + assert( o.is() ); + o = intersect( Segment_2(pts[2], pts[0]), Segment_2(pts[3], pts[1]) ); + assert( o.is() ); + // segment fully included + o = intersect( Segment_2(pts[0], pts[3]), Segment_2(pts[1], pts[2]) ); + assert( o.is() ); + assert( CGAL::object_cast(o) == Segment_2(pts[1], pts[2]) ); + // segment fully included with shared vertex + o = intersect( Segment_2(pts[0], pts[1]), Segment_2(pts[0], pts[2]) ); + assert( o.is() ); + assert( CGAL::object_cast(o) == Segment_2(pts[0], pts[1]) ); + // segment sharing a vertex + o = intersect( Segment_2(pts[0], pts[1]), Segment_2(pts[1], pts[2]) ); + assert( o.is() ); + assert( CGAL::object_cast(o) == pts[1]); + + // degenerate segment + Segment_2 sd(Point_2(1,0,6), Point_2(1,0,7)); + o = intersect( Segment_2(pts[0], pts[2]), sd ); + assert( o.is() ); + o = intersect( Segment_2(pts[0], pts[1]), sd ); + assert( o.is() ); + o = intersect( Segment_2(pts[1], pts[0]), sd ); + assert( o.is() ); + o = intersect( Segment_2(pts[2], pts[3]), sd ); + assert( o.empty() ); + + return true; +} + int main() { std::cerr.precision(17); @@ -75,5 +136,12 @@ int main() test > ("CDT_2 in a 3D plane, with Epeck"); + ok = ok && test_segment_intersections > + ("CDT_2 traits intersection with Epick"); + ok = ok && test_segment_intersections > + ("CDT_2 traits intersection with Epeck"); + return ok ? 0 : 1; }