diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h b/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h index 269c86c3e606..471904cd97ab 100644 --- a/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h +++ b/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h @@ -39,11 +39,15 @@ struct Point_on_triangle { // triangle points are not stored in this class but are expected // to always be passed in the same order. For a triangle pqr, - // edge 0 is pq, edge 1 qr and edge 2 rp. Point 0 is p, 1 is q and 2 is r. + // edge 1 is pq, edge 2 qr and edge 3 rp. Point -1 is p, -2 is q and -3 is r. + // + // (0,-) vertex of t2 inside t1 + // (-,-) common point of t1 and t2 + // (-,0) vertex of t1 inside t2 + // (-,+) vertex of t1 on an edge of t2 + // (+,+) edge of t1 intersecting edge of t2 + // (+,-) vertex of t2 on an edge of t1 // - // (id, -1) point on t1 - // (-1, id) point on t2 - // (id1, id2) intersection of edges std::pair t1_t2_ids; boost::container::flat_set extra_t1; // store other ids of edges containing the point typename Kernel::FT alpha; // @@ -60,16 +64,16 @@ struct Point_on_triangle { switch(id) { - case 0: + case -1: return p; - case 1: + case -2: return q; default: return r; } } - Point_on_triangle(int i1, int i2=-1, typename Kernel::FT alpha = 0.) // TODO add global zero()? + Point_on_triangle(int i1, int i2, typename Kernel::FT alpha = 0.) // TODO add global zero()? : t1_t2_ids(i1,i2) , alpha(alpha) {} @@ -85,20 +89,20 @@ struct Point_on_triangle const typename Kernel::Point_3& r2, const Kernel& k) const { - if (t1_t2_ids.first!=-1) + if (t1_t2_ids.first!=0 && t1_t2_ids.second >=0) { - if (t1_t2_ids.second==-1) + if (t1_t2_ids.first<0) return (edge_id1==t1_t2_ids.first || (edge_id1+1)%3==t1_t2_ids.first) ? ZERO:POSITIVE; // it is a point on t1 // this is an intersection point if (t1_t2_ids.first==edge_id1) return ZERO; - if (t1_t2_ids.first==(edge_id1+1)%3) + if (t1_t2_ids.first-1==(edge_id1)%3) { if (alpha==0) return ZERO; return alpha>=0 ? POSITIVE:NEGATIVE; } - CGAL_assertion((t1_t2_ids.first+1)%3==edge_id1); + CGAL_assertion((t1_t2_ids.first)%3==edge_id1-1); if (alpha==1) return ZERO; return alpha<=1?POSITIVE:NEGATIVE; } @@ -111,9 +115,13 @@ struct Point_on_triangle } } + int id1() const { return t1_t2_ids.first; } int id2() const { return t1_t2_ids.second; } + void set_id1(int i) { t1_t2_ids.first=i; } + void set_id2(int i) { t1_t2_ids.second=i; } + // construct the intersection point from the info stored typename Kernel::Point_3 point(const typename Kernel::Point_3& p1, @@ -124,12 +132,13 @@ struct Point_on_triangle const typename Kernel::Point_3& r2, const Kernel& k) const { - if (t1_t2_ids.first==-1) + if (t1_t2_ids.second<0) return point_from_id(p2,q2,r2,t1_t2_ids.second); - if (t1_t2_ids.second==-1) + if (t1_t2_ids.first<0) return point_from_id(p1,q1,r1,t1_t2_ids.first); - return k.construct_barycenter_3_object()(point_from_id(p1,q1,r1,(t1_t2_ids.first+1)%3), alpha, point_from_id(p1,q1,r1,t1_t2_ids.first)) ; + return k.construct_barycenter_3_object()(point_from_id(p1,q1,r1,-(t1_t2_ids.first)%3-1), alpha, + point_from_id(p1,q1,r1,-t1_t2_ids.first)) ; } }; @@ -155,158 +164,55 @@ intersection(const Point_on_triangle& p, #ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION std::cout << " calling intersection: "; std::cout << " (" << p.id1() << "," << p.id2() << ",[" << p.alpha << "]) -"; - std::cout << " (" << q.id1() << "," << q.id2() << ",[" << q.alpha << "]) || e" << edge_id_t1; + std::cout << " (" << q.id1() << "," << q.id2() << ",[" << q.alpha << "]) || e" << edge_id_t1 << "\n"; #endif typename Kernel::Compute_alpha_for_coplanar_triangle_intersection_3 compute_alpha = k.compute_alpha_for_coplanar_triangle_intersection_3_object(); typedef Point_on_triangle Pot; - switch(p.id1()) + + auto common_edge_id = [](int pid, int qid) { - case -1: + if (pid==0 || qid==0) return 0; + if (pid<0) { - switch(q.id1()) + if (qid<0) { - case -1: // A: (-1, ip2) - (-1, iq2) - { - CGAL_assertion((p.id2()+1)%3 == q.id2() || (q.id2()+1)%3 == p.id2()); -// CGAL_assertion(p.extra_t1.empty() && q.extra_t1.empty()); // TMP to see if it's worth implementing special case -#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " -- case 1\n"; -#endif - typename Kernel::FT alpha = compute_alpha(p1, q1, - Pot::point_from_id(p2, q2, r2, p.id2()), - Pot::point_from_id(p2, q2, r2, q.id2())); - int id2 = (p.id2()+1)%3 == q.id2() ? p.id2() : q.id2(); - return Point_on_triangle(edge_id_t1, id2, alpha); // intersection with an original edge of t2 - } - default: - if (q.id2()!=-1) // B: (-1, ip2) - (iq1, iq2) - { - if (p.id2() == q.id2() || p.id2() == (q.id2()+1)%3) - { -#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " -- case 2\n"; -#endif - // points are on the same edge of t2 --> we shorten an already cut edge - typename Kernel::FT alpha = compute_alpha(p1, q1, - Pot::point_from_id(p2, q2, r2, q.id2()), - Pot::point_from_id(p2, q2, r2, (q.id2()+1)%3)); - - return Point_on_triangle(edge_id_t1, q.id2(), alpha); - } -#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " -- case 3\n"; -#endif - // point of t1: look for an edge of t1 containing both points - CGAL_assertion( p.extra_t1.count(q.id1())!=0 || p.extra_t1.count(3-q.id1()-edge_id_t1)!=0 ); - int eid1 = p.extra_t1.count(q.id1())!=0 ? q.id1() : 3-q.id1()-edge_id_t1; - return Point_on_triangle((eid1+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3, -1); // vertex of t1 - } - // C: (-1, ip2) - (iq1, -1) - //vertex of t1, special case t1 edge passed thru a vertex of t2 - CGAL_assertion(edge_id_t1 == 2); -#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " -- case 4\n"; -#endif - CGAL_assertion(q.id1()==1); - CGAL_assertion(!p.extra_t1.empty()); - return Point_on_triangle(p.extra_t1.count(0)==1?0:2,-1); + return (-pid)%3==-qid-1 ? -pid:-qid; + } + else + { + return (-pid==qid || (qid)%3==-pid-1) ? qid : 0; } } - default: + else { - switch(p.id2()) + if (qid<0) { - case -1: - { - switch(q.id1()) - { - case -1: // G: (ip1, -1) - (-1, iq2) - //vertex of t1, special case t1 edge passed thru a vertex of t2 -#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " -- case 5\n"; -#endif - CGAL_assertion(edge_id_t1 == 2); - CGAL_assertion(p.id1()==1); - CGAL_assertion(!q.extra_t1.empty()); - return Point_on_triangle(q.extra_t1.count(0)==1?0:2,-1); - default: - { -#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " -- case 6\n"; -#endif - CGAL_assertion(q.id2()!=-1); // I: (ip1, -1) - (iq2, -1) - //H: (ip1,-1), (iq1, iq2) - CGAL_assertion(edge_id_t1==2); - // p and q are on the same edge of t1 - CGAL_assertion(p.id1()==q.id1() || p.id1()==(q.id1()+1)%3); - return Point_on_triangle((q.id1()+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3 , -1); - } - } - } - default: - { - switch(q.id1()) - { - case -1: // D: (ip1, ip2) - (-1, iq2) - { - if (q.id2() == p.id2() || q.id2() == (p.id2()+1)%3) - { -#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " -- case 7\n"; -#endif - // points are on the same edge of t2 --> we shorten an already cut edge - typename Kernel::FT alpha = compute_alpha(p1, q1, - Pot::point_from_id(p2, q2, r2, p.id2()), - Pot::point_from_id(p2, q2, r2, (p.id2()+1)%3)); - - return Point_on_triangle(edge_id_t1, p.id2(), alpha); - } -#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " -- case 8\n"; -#endif - // point of t1 - //std::cout << "q.extra_t1: "; for(int qet1 : q.extra_t1) std::cout << " " << qet1; std::cout << "\n"; - CGAL_assertion( q.extra_t1.count(p.id1())!=0 || q.extra_t1.count(3-p.id1()-edge_id_t1)!=0 ); - int eid1 = q.extra_t1.count(p.id1())!=0 ? p.id1() : 3-p.id1()-edge_id_t1; - return Point_on_triangle((eid1+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3, -1); // vertex of t1 - } - default: - { - switch(q.id2()) - { - case -1: // F: (ip1, ip2) - (iq1, -1) - { - // p and q are on the same edge of t1 - CGAL_assertion(q.id1()==p.id1() || q.id1()==(p.id1()+1)%3); -#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " -- case 9\n"; -#endif - return Point_on_triangle((p.id1()+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3 , -1); - } - default: // E: (ip1, ip2) - (iq1, iq2) - { - if (p.id2()==q.id2()) - { -#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " -- case 10\n"; -#endif - typename Kernel::FT alpha = compute_alpha(p1, q1, - Pot::point_from_id(p2, q2, r2, q.id2()), - Pot::point_from_id(p2, q2, r2, (q.id2()+1)%3)); - return Point_on_triangle(edge_id_t1, q.id2(), alpha); - } - // we are intersecting an edge of t1 - CGAL_assertion(p.id1()==q.id1() || edge_id_t1==2); - int eid1 = p.id1()==q.id1() ? p.id1() : 1; - return Point_on_triangle((eid1+1)%3==edge_id_t1?edge_id_t1:(edge_id_t1+1)%3, -1); // vertex of t1 - } - } - } - } - } + return (-qid==pid || (pid)%3==-qid-1) ? pid : 0; + } + else + { + return pid==qid ? pid : 0; } } + }; + + int common_eid1 = common_edge_id(p.id1(), q.id1()); + int common_eid2 = common_edge_id(p.id2(), q.id2()); + + if (common_eid1!=0) + { + CGAL_assertion( (common_eid1)%3==edge_id_t1-1 || (edge_id_t1)%3==common_eid1-1 ); + int pid1 = (common_eid1)%3==edge_id_t1-1? -edge_id_t1:-common_eid1; + return Point_on_triangle(pid1, common_eid2); + } + else + { + CGAL_assertion(common_eid2!=0); + typename Kernel::FT alpha = compute_alpha(p1, q1, + Pot::point_from_id(p2, q2, r2, -common_eid2), + Pot::point_from_id(p2, q2, r2, -(common_eid2)%3-1)); + return Point_on_triangle(edge_id_t1, common_eid2, alpha); } } @@ -336,8 +242,17 @@ void intersection_coplanar_triangles_cutoff(const typename Kernel::Point_3& p1, for (Point_on_triangle& pot : inter_pts) { orientations[ &pot ]=pot.orientation(p1,q1,r1,edge_id,p2,q2,r2,k); - if (pot.id1()==-1 && orientations[ &pot ]==COLLINEAR) - pot.extra_t1.insert(edge_id); + if (orientations[ &pot ]==COLLINEAR) + { + if (pot.id1()==0) + pot.set_id1(edge_id); + else + { + CGAL_assertion(pot.id1()>0); + CGAL_assertion((edge_id)%3==pot.id1()-1 || (pot.id1())%3==edge_id-1); + pot.set_id1( (edge_id)%3==pot.id1()-1 ? -pot.id1(): -edge_id ); + } + } } #ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION @@ -411,32 +326,29 @@ intersection_coplanar_triangles(const typename K::Triangle_3& t1, r2 = t2.vertex(2); std::list> inter_pts; - inter_pts.push_back(Point_on_triangle(-1,0)); - inter_pts.push_back(Point_on_triangle(-1,1)); - inter_pts.push_back(Point_on_triangle(-1,2)); + inter_pts.push_back(Point_on_triangle(0,-1)); + inter_pts.push_back(Point_on_triangle(0,-2)); + inter_pts.push_back(Point_on_triangle(0,-3)); #ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION auto print_points = [&]() { - for(auto p : inter_pts) std::cout << " (" << p.id1() << "," << p.id2() << ",[" << p.alpha << "]) "; std::cout <<"\n"; + std::cout << " ipts size: " << inter_pts.size() << "\n"; + for(auto p : inter_pts) {std::cout << " (" << p.id1() << "," << p.id2() << ",[" << p.alpha << "]) ";} std::cout <<"\n"; }; - std::cout << " ipts size: " << inter_pts.size() << "\n"; print_points(); #endif //intersect t2 with the three half planes which intersection defines t1 - intersection_coplanar_triangles_cutoff(p1,q1,r1,0,p2,q2,r2,k,inter_pts); //line pq + intersection_coplanar_triangles_cutoff(p1,q1,r1,1,p2,q2,r2,k,inter_pts); //line pq #ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " ipts size: " << inter_pts.size() << "\n"; print_points(); #endif - intersection_coplanar_triangles_cutoff(q1,r1,p1,1,p2,q2,r2,k,inter_pts); //line qr + intersection_coplanar_triangles_cutoff(q1,r1,p1,2,p2,q2,r2,k,inter_pts); //line qr #ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " ipts size: " << inter_pts.size() << "\n"; print_points(); #endif - intersection_coplanar_triangles_cutoff(r1,p1,q1,2,p2,q2,r2,k,inter_pts); //line rp + intersection_coplanar_triangles_cutoff(r1,p1,q1,3,p2,q2,r2,k,inter_pts); //line rp #ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION - std::cout << " ipts size: " << inter_pts.size() << "\n"; print_points(); #endif diff --git a/Kernel_23/include/CGAL/Kernel_23/internal/Projection_traits_3.h b/Kernel_23/include/CGAL/Kernel_23/internal/Projection_traits_3.h index dd4e24005a48..1cc09421944d 100644 --- a/Kernel_23/include/CGAL/Kernel_23/internal/Projection_traits_3.h +++ b/Kernel_23/include/CGAL/Kernel_23/internal/Projection_traits_3.h @@ -20,6 +20,7 @@ #include #include +#include namespace CGAL { @@ -1210,6 +1211,14 @@ class Projection_traits_3 }; -} } //namespace CGAL::internal + +} // namespace internal + +template +struct Triangulation_structural_filtering_traits<::CGAL::internal::Projection_traits_3 > { + typedef typename Triangulation_structural_filtering_traits::Use_structural_filtering_tag Use_structural_filtering_tag; +}; + + } //namespace CGAL #endif // CGAL_INTERNAL_PROJECTION_TRAITS_3_H diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h index a3606f816a7b..8de62f6b8021 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/autorefinement.h @@ -55,6 +55,10 @@ #endif #endif +#ifdef CGAL_AUTOREF_USE_FIXED_PROJECTION_TRAITS +#include +#endif + #include //#define CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS @@ -114,7 +118,6 @@ Segment_inter_type do_coplanar_segments_intersect(std::size_t pi, std::size_t qi, std::size_t ri, std::size_t si, const std::vector& points, - const typename K::Vector_3& /* plane_normal */, const K& k = K()) { typename K::Collinear_are_ordered_along_line_3 cln_order = k.collinear_are_ordered_along_line_3_object(); @@ -249,36 +252,55 @@ do_coplanar_segments_intersect(std::size_t pi, std::size_t qi, return NO_INTERSECTION; } + // imported from Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Triangle_3_intersection.h template void coplanar_intersections(const std::array& t1, const std::array& t2, - std::vector& inter_pts) + std::vector>& inter_pts) { - const typename K::Point_3& p1 = t1[0], q1 = t1[1], r1 = t1[2]; - const typename K::Point_3& p2 = t2[0], q2 = t2[1], r2 = t2[2]; + const typename K::Point_3& p1 = t1[0], &q1 = t1[1], &r1 = t1[2]; + const typename K::Point_3& p2 = t2[0], &q2 = t2[1], &r2 = t2[2]; std::list> l_inter_pts; - l_inter_pts.push_back(Intersections::internal::Point_on_triangle(-1,0)); - l_inter_pts.push_back(Intersections::internal::Point_on_triangle(-1,1)); - l_inter_pts.push_back(Intersections::internal::Point_on_triangle(-1,2)); + l_inter_pts.push_back(Intersections::internal::Point_on_triangle(0,-1)); + l_inter_pts.push_back(Intersections::internal::Point_on_triangle(0,-2)); + l_inter_pts.push_back(Intersections::internal::Point_on_triangle(0,-3)); +#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION + auto print_points = [&]() + { + std::cout << " ipts size: " << l_inter_pts.size() << "\n"; + for(auto p : l_inter_pts) {std::cout << " (" << p.id1() << "," << p.id2() << ",[" << p.alpha << "]) ";} std::cout <<"\n"; + }; +#endif //intersect t2 with the three half planes which intersection defines t1 K k; - Intersections::internal::intersection_coplanar_triangles_cutoff(p1,q1,r1,0,p2,q2,r2,k,l_inter_pts); //line p1q1 - Intersections::internal::intersection_coplanar_triangles_cutoff(q1,r1,p1,1,p2,q2,r2,k,l_inter_pts); //line q1r1 - Intersections::internal::intersection_coplanar_triangles_cutoff(r1,p1,q1,2,p2,q2,r2,k,l_inter_pts); //line r1p1 +#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION + print_points(); +#endif + Intersections::internal::intersection_coplanar_triangles_cutoff(p1,q1,r1,1,p2,q2,r2,k,l_inter_pts); //line p1q1 +#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION + print_points(); +#endif + Intersections::internal::intersection_coplanar_triangles_cutoff(q1,r1,p1,2,p2,q2,r2,k,l_inter_pts); //line q1r1 +#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION + print_points(); +#endif + Intersections::internal::intersection_coplanar_triangles_cutoff(r1,p1,q1,3,p2,q2,r2,k,l_inter_pts); //line r1p1 +#ifdef CGAL_DEBUG_COPLANAR_T3_T3_INTERSECTION + print_points(); +#endif for (const Intersections::internal::Point_on_triangle& pot : l_inter_pts) - inter_pts.push_back( pot.point(p1,q1,r1,p2,q2,r2,k) ); + inter_pts.emplace_back(pot.point(p1,q1,r1,p2,q2,r2,k), pot.id1(), pot.id2()); } -// imported from Polygon_mesh_processing/internal/Corefinement/intersect_triangle_segment_3.h +// imported from Polygon_mesh_processing/internal/Corefinement/intersect_triangle_and_segment_3.h template -void +std::optional> find_intersection(const typename K::Point_3& p, const typename K::Point_3& q, //segment const typename K::Point_3& a, const typename K::Point_3& b, const typename K::Point_3& c, //triangle - std::vector& inter_pts, bool is_p_coplanar=false, bool is_q_coplanar=false) // note that in coref this was wrt a halfedge not p/q { Orientation ab=orientation(p,q,a,b); @@ -286,54 +308,56 @@ find_intersection(const typename K::Point_3& p, const typename K::Point_3& q, / Orientation ca=orientation(p,q,c,a); if ( ab==POSITIVE || bc==POSITIVE || ca==POSITIVE ) - return; + return std::nullopt; int nb_coplanar=(ab==COPLANAR?1:0) + (bc==COPLANAR?1:0) + (ca==COPLANAR?1:0); - if (is_p_coplanar) - { - inter_pts.push_back(p); - return; - } - if (is_q_coplanar) - { - inter_pts.push_back(q); - return; - } - - if (nb_coplanar!=2) - { - inter_pts.push_back( - typename K::Construct_plane_line_intersection_point_3()(a, b, c, p, q) - ); - } - else + if (nb_coplanar==2) { + // even if a common point is not new it is still needed to be reported so + // that the intersection segment is known. if (ab!=COPLANAR) - { // intersection is c - inter_pts.push_back(c); - return; - } - - if (bc!=COPLANAR) - { + return std::make_pair(c, -3); + else if (bc!=COPLANAR) // intersection is a - inter_pts.push_back(a); - return; + return std::make_pair(a, -1); + else + { + CGAL_assertion(ca!=COPLANAR); + // intersection is b + return std::make_pair(b, -2); } - CGAL_assertion(ca!=COPLANAR); - // intersection is b - inter_pts.push_back(b); } + + typename K::Point_3 ipt = is_p_coplanar ? p : + is_q_coplanar ? q : + typename K::Construct_plane_line_intersection_point_3() + (a, b, c, p, q); + + if (nb_coplanar == 0) + return std::make_pair(ipt, 0); + + + CGAL_assertion(nb_coplanar==1); + + if (ab==COPLANAR) + // intersection is ab + return std::make_pair(ipt, 1); + if (bc==COPLANAR) + // intersection is bc + return std::make_pair(ipt, 2); + CGAL_assertion(ca==COPLANAR); + // intersection is ca + return std::make_pair(ipt, 3); } template -void test_edge(const typename K::Point_3& p, const typename K::Point_3& q, - const typename K::Point_3& a, const typename K::Point_3& b, const typename K::Point_3& c, - const Orientation abcp, - const Orientation abcq, - std::vector& inter_pts) +std::optional> +test_edge(const typename K::Point_3& p, const typename K::Point_3& q, + const typename K::Point_3& a, const typename K::Point_3& b, const typename K::Point_3& c, + const Orientation abcp, + const Orientation abcq) { switch ( abcp ) { case POSITIVE: @@ -341,46 +365,40 @@ void test_edge(const typename K::Point_3& p, const typename K::Point_3& q, case POSITIVE: // the segment lies in the positive open halfspaces defined by the // triangle's supporting plane - break; + return std::nullopt; case NEGATIVE: // p sees the triangle in counterclockwise order - find_intersection(p,q,a,b,c,inter_pts); - break; + return find_intersection(p,q,a,b,c); //case COPLANAR: default: // q belongs to the triangle's supporting plane // p sees the triangle in counterclockwise order - find_intersection(p,q,a,b,c,inter_pts,false,true); + return find_intersection(p,q,a,b,c,false,true); } - break; case NEGATIVE: switch ( abcq ) { case POSITIVE: // q sees the triangle in counterclockwise order - find_intersection(q,p,a,b,c,inter_pts); - break; + return find_intersection(q,p,a,b,c); case NEGATIVE: // the segment lies in the negative open halfspaces defined by the // triangle's supporting plane - break; + return std::nullopt; // case COPLANAR: default: // q belongs to the triangle's supporting plane // p sees the triangle in clockwise order - find_intersection(q,p,a,b,c,inter_pts,true,false); + return find_intersection(q,p,a,b,c,true,false); } - break; default: //case COPLANAR: // p belongs to the triangle's supporting plane switch ( abcq ) { case POSITIVE: // q sees the triangle in counterclockwise order - find_intersection(q,p,a,b,c,inter_pts,false, true); - break; + return find_intersection(q,p,a,b,c,false, true); case NEGATIVE: // q sees the triangle in clockwise order - find_intersection(p,q,a,b,c,inter_pts,true); - break; + return find_intersection(p,q,a,b,c,true); //case COPLANAR: default: // the segment is coplanar with the triangle's supporting plane @@ -392,7 +410,7 @@ void test_edge(const typename K::Point_3& p, const typename K::Point_3& q, // nothing done as coplanar case handle in collect_intersections // and other intersection points will be collected with non-coplanar edges //} - break; + return std::nullopt; } } } @@ -400,7 +418,7 @@ void test_edge(const typename K::Point_3& p, const typename K::Point_3& q, template bool collect_intersections(const std::array& t1, const std::array& t2, - std::vector& inter_pts) + std::vector>& inter_pts) { // test edges of t1 vs t2 std::array ori; @@ -409,6 +427,7 @@ bool collect_intersections(const std::array& t1, if (ori[0]== COPLANAR && ori[1]==COPLANAR && ori[2]==COPLANAR) { + coplanar_intersections(t1, t2, inter_pts); #ifdef CGAL_AUTOREF_DEBUG_DEPTH for (auto p : inter_pts) @@ -421,7 +440,17 @@ bool collect_intersections(const std::array& t1, for (int i=0; i<3; ++i) { int j=(i+1)%3; - test_edge(t1[i], t1[j], t2[0], t2[1], t2[2], ori[i], ori[j], inter_pts); + std::optional > opt = + test_edge(t1[i], t1[j], t2[0], t2[1], t2[2], ori[i], ori[j]); + if (opt) + { + if (ori[i]==COPLANAR) + inter_pts.emplace_back(opt->first, -(i+1), opt->second); + else if (ori[j]==COPLANAR) + inter_pts.emplace_back(opt->first, -(j+1), opt->second); + else + inter_pts.emplace_back(opt->first, i+1 , opt->second); + } } // test edges of t2 vs t1 @@ -430,12 +459,23 @@ bool collect_intersections(const std::array& t1, for (int i=0; i<3; ++i) { int j=(i+1)%3; - test_edge(t2[i], t2[j], t1[0], t1[1], t1[2], ori[i], ori[j], inter_pts); + std::optional > opt = + test_edge(t2[i], t2[j], t1[0], t1[1], t1[2], ori[i], ori[j]); + if (opt) + { + if (ori[i]==COPLANAR) + inter_pts.emplace_back(opt->first, opt->second, -(i+1)); + else if (ori[j]==COPLANAR) + inter_pts.emplace_back(opt->first, opt->second, -(j+1)); + else + inter_pts.emplace_back(opt->first, opt->second, i+1 ); + } } +// #warning TODO get rid of sort and unique calls // because we don't handle intersection type and can have edge-edge edge-vertex duplicates - std::sort(inter_pts.begin(), inter_pts.end()); - auto last = std::unique(inter_pts.begin(), inter_pts.end()); + std::sort(inter_pts.begin(), inter_pts.end(), [](auto p, auto q){return std::get<0>(p)(q);}); + auto last = std::unique(inter_pts.begin(), inter_pts.end(), [](auto p, auto q){return std::get<0>(p)==std::get<0>(q);}); inter_pts.erase(last, inter_pts.end()); #ifdef CGAL_AUTOREF_DEBUG_DEPTH @@ -446,45 +486,75 @@ bool collect_intersections(const std::array& t1, return false; } +struct Triangle_data +{ + using Point_3 = Exact_predicates_exact_constructions_kernel::Point_3; + std::vector points; + std::vector point_locations; + std::vector> segments; + std::vector segment_input_triangle_ids; + template + std::size_t add_point(const std::tuple& tpl) + { + if (std::get(tpl) < 0) + return -std::get(tpl)-1; + points.push_back(std::get<0>(tpl)); + point_locations.push_back(std::get(tpl)); + return points.size()-1; + } + +#ifndef NDEBUG + bool on_same_edge(std::size_t i, std::size_t j) +#else + bool on_same_edge_debug(std::size_t i, std::size_t j) +#endif + { + if (point_locations[i]==0 || point_locations[j]==0) return false; + if (point_locations[i]>0) + { + if (point_locations[j]>0) + return point_locations[i]==point_locations[j]; + return -point_locations[j]==point_locations[i] || point_locations[i]%3+1==-point_locations[j]; + } + if (point_locations[j]<0) + return true; + return -point_locations[i]==point_locations[j] || point_locations[j]%3+1==-point_locations[i]; + } +#ifdef NDEBUG + bool on_same_edge(std::size_t i, std::size_t j) + { + bool on = on_same_edge_debug(i,j); + if (!on) return false; + int eid = point_locations[i]>0 ? point_locations[i] : point_locations[j]; + if (eid<0) return true; + if (!(collinear(points[eid-1], points[(eid)%3], points[i]))) + { + std::cout << point_locations[i] << " " << point_locations[j] << " " << eid << "\n"; + } + if (!(collinear(points[eid-1], points[(eid)%3], points[j]))) + { + std::cout << point_locations[i] << " " << point_locations[j] << " " << eid << "\n"; + } + CGAL_assertion(collinear(points[eid-1], points[(eid)%3], points[i])); + CGAL_assertion(collinear(points[eid-1], points[(eid)%3], points[j])); + return true; + } +#endif +}; + template void generate_subtriangles(std::size_t ti, - std::vector>& segments, - std::vector& points, - const std::vector& in_triangle_ids, + Triangle_data& triangle_data, const std::set >& intersecting_triangles, const std::set >& coplanar_triangles, const std::vector>& triangles, PointVector& new_triangles ) { - typedef CGAL::Projection_traits_3 P_traits; - typedef CGAL::No_constraint_intersection_tag Itag; - - typedef CGAL::Constrained_Delaunay_triangulation_2 CDT_2; - //typedef CGAL::Constrained_triangulation_plus_2 CDT; - typedef CDT_2 CDT; - - const std::array& t = triangles[ti]; - - // positive triangle normal - typename EK::Vector_3 n = normal(t[0], t[1], t[2]); - typename EK::Point_3 o(CGAL::ORIGIN); - - bool orientation_flipped = false; - if ( typename EK::Less_xyz_3()(o+n,o) ) - { - n=-n; - orientation_flipped = true; - } - - P_traits cdt_traits(n); - CDT cdt(cdt_traits); - cdt.insert_outside_affine_hull(t[0]); - cdt.insert_outside_affine_hull(t[1]); - typename CDT::Vertex_handle v = cdt.tds().insert_dim_up(cdt.infinite_vertex(), orientation_flipped); - v->set_point(t[2]); - #ifdef CGAL_AUTOREFINE_DEBUG_COUNTERS struct Counter { @@ -518,7 +588,12 @@ void generate_subtriangles(std::size_t ti, #define CGAL_AUTOREF_COUNTER_INSTRUCTION(X) #endif - // pre-compute segment intersections + std::vector& points=triangle_data.points; + std::vector& point_locations=triangle_data.point_locations; + std::vector>& segments=triangle_data.segments; + std::vector& in_triangle_ids=triangle_data.segment_input_triangle_ids; + +// pre-compute segment intersections if (!segments.empty()) { std::size_t nbs = segments.size(); @@ -528,7 +603,7 @@ void generate_subtriangles(std::size_t ti, CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer1.start();) std::map point_id_map; - +//TODO: we already have sorted the points while deduplicating segments! for (std::size_t pid=0; pidsecond; }; @@ -569,7 +647,7 @@ void generate_subtriangles(std::size_t ti, Segment_inter_type seg_inter_type = do_coplanar_segments_intersect(segments[i].first, segments[i].second, segments[j].first, segments[j].second, - points, n); + points); CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer5.stop();) switch(seg_inter_type) @@ -739,13 +817,147 @@ void generate_subtriangles(std::size_t ti, auto last = std::unique(segments.begin(), segments.end()); segments.erase(last, segments.end()); +// init CDT + insert points and constraints CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer3.start();) - if (segments.empty()) - cdt.insert(points.begin(), points.end()); - else - cdt.insert_constraints(points.begin(), points.end(), segments.begin(), segments.end()); +#ifdef CGAL_AUTOREF_USE_FIXED_PROJECTION_TRAITS + typedef ::CGAL::internal::Projection_traits_3 P_traits; +#else + typedef CGAL::Projection_traits_3 P_traits; +#endif + + typedef CGAL::No_constraint_intersection_tag Itag; + + typedef CGAL::Constrained_Delaunay_triangulation_2 CDT_2; + //typedef CGAL::Constrained_triangulation_plus_2 CDT; + typedef CDT_2 CDT; + + const std::array& t = triangles[ti]; + std::vector vhandles(triangle_data.points.size()); + +#ifdef CGAL_AUTOREF_USE_FIXED_PROJECTION_TRAITS + P_traits cdt_traits; + bool orientation_flipped = false; + CDT cdt(cdt_traits); + // TODO: still need to figure out why I can't make the orientation_flipped correctly + vhandles[0]=cdt.insert(t[0]); + vhandles[1]=cdt.insert(t[1]); + vhandles[2]=cdt.insert(t[2]); +#else + // positive triangle normal + typename EK::Vector_3 n = normal(t[0], t[1], t[2]); + typename EK::Point_3 o(CGAL::ORIGIN); + + bool orientation_flipped = false; + if ( typename EK::Less_xyz_3()(o+n,o) ) + { + n=-n; + orientation_flipped = true; + } + P_traits cdt_traits(n); + + CDT cdt(cdt_traits); + vhandles[0]=cdt.insert_outside_affine_hull(t[0]); + vhandles[1]=cdt.insert_outside_affine_hull(t[1]); + vhandles[2] = cdt.tds().insert_dim_up(cdt.infinite_vertex(), orientation_flipped); + vhandles[2]->set_point(t[2]); +#endif + + // insert points and fill vhandles +#if 0 + std::vector indices(triangle_data.points.size()-3); + std::iota(indices.begin(), indices.end(), 3); +#else + //start by points on edges + std::array, 3> indices_on_edges; + std::vector indices; + for (std::size_t i=3; i0 && point_locations[i]<4); + indices_on_edges[point_locations[i]-1].push_back(i); + } + } + + //sort points on edges and insert them + typename CDT::Vertex_handle vinf=cdt.infinite_vertex(); + for (int i=0; i<3; ++i) + { + if (indices_on_edges[i].empty()) continue; + std::size_t src_id=i, tgt_id=(i+1)%3; + //look for a sort axis + int coord = 0; + if (points[src_id].x()==points[tgt_id].x()) + { + coord=1; + if (points[src_id].y()==points[tgt_id].y()) + coord=2; + } + if (points[src_id][coord]>points[tgt_id][coord]) + std::swap(src_id, tgt_id); + + std::sort(indices_on_edges[i].begin(), indices_on_edges[i].end(), + [&](std::size_t id1, std::size_t id2) + { + return points[id1][coord]index(vinf)); + cdt.restore_Delaunay(vhandles[id]); // TODO maybe not each time but one global? + CGAL_assertion(cdt.is_valid()); + prev_id=id; + } + } +#endif + // then points in the interior + typedef typename Pointer_property_map::type Pmap; + typedef Spatial_sort_traits_adapter_2 Search_traits; + spatial_sort(indices.begin(), indices.end(), + Search_traits(make_property_map(points), cdt_traits)); + + typename CDT::Face_handle hint; + for (std::size_t i : indices) + { + vhandles[i] = cdt.insert(points[i], hint); + hint=vhandles[i]->face(); + } + + for (const std::pair& ids : triangle_data.segments) + { + CGAL_assertion(ids.first < vhandles.size()); + CGAL_assertion(ids.second < vhandles.size()); + CGAL_assertion( vhandles[ids.first]!= typename CDT::Vertex_handle() ); + CGAL_assertion( vhandles[ids.second]!= typename CDT::Vertex_handle() ); + cdt.insert_constraint(vhandles[ids.first], vhandles[ids.second]); + } CGAL_AUTOREF_COUNTER_INSTRUCTION(counter.timer3.stop();) +// extract new triangles for (typename CDT::Face_handle fh : cdt.finite_face_handles()) { if (orientation_flipped) @@ -895,21 +1107,30 @@ void autorefine_triangle_soup(PointRange& soup_points, // init the vector of triangles used for the autorefinement of triangles typedef CGAL::Exact_predicates_exact_constructions_kernel EK; + // even if the info is duplicated with Triangle_data::point, we keep this container + // so that we can use it in parallel calls of generate_subtriangles std::vector< std::array > triangles(tiid+1); + // vector of data for refining triangles + std::vector all_triangle_data(triangles.size()); Cartesian_converter to_exact; for(Input_TID f : intersected_faces) { - triangles[tri_inter_ids[f]]= CGAL::make_array( + std::size_t tid=tri_inter_ids[f]; + triangles[tid]= CGAL::make_array( to_exact( get(pm, soup_points[soup_triangles[f][0]]) ), to_exact( get(pm, soup_points[soup_triangles[f][1]]) ), to_exact( get(pm, soup_points[soup_triangles[f][2]]) ) ); + all_triangle_data[tid].points.resize(3); + all_triangle_data[tid].points[0]=triangles[tri_inter_ids[f]][0]; + all_triangle_data[tid].points[1]=triangles[tri_inter_ids[f]][1]; + all_triangle_data[tid].points[2]=triangles[tri_inter_ids[f]][2]; + all_triangle_data[tid].point_locations.resize(3); + all_triangle_data[tid].point_locations[0]=-1; + all_triangle_data[tid].point_locations[1]=-2; + all_triangle_data[tid].point_locations[2]=-3; } - std::vector< std::vector > > all_segments(triangles.size()); - std::vector< std::vector > all_points(triangles.size()); - std::vector< std::vector > all_in_triangle_ids(triangles.size()); - CGAL_PMP_AUTOREFINE_VERBOSE("compute intersections"); #ifdef CGAL_AUTOREF_USE_DEBUG_PARALLEL_TIMERS Real_timer t; @@ -928,7 +1149,7 @@ void autorefine_triangle_soup(PointRange& soup_points, const std::array& t1 = triangles[i1]; const std::array& t2 = triangles[i2]; - std::vector inter_pts; + std::vector> inter_pts; bool triangles_are_coplanar = autorefine_impl::collect_intersections(t1, t2, inter_pts); CGAL_assertion( @@ -941,23 +1162,54 @@ void autorefine_triangle_soup(PointRange& soup_points, switch(nbi) { case 1: - all_points[i1].push_back(inter_pts[0]); - all_points[i2].push_back(inter_pts[0]); + all_triangle_data[i1].add_point<1>(inter_pts[0]); + all_triangle_data[i2].add_point<2>(inter_pts[0]); break; case 2: - all_segments[i1].push_back(CGAL::make_array(inter_pts[0], inter_pts[1])); - all_segments[i2].push_back(CGAL::make_array(inter_pts[0], inter_pts[1])); - all_in_triangle_ids[i1].push_back(i2); - all_in_triangle_ids[i2].push_back(i1); + { + std::size_t src_id=all_triangle_data[i1].add_point<1>(inter_pts[0]), + tgt_id=all_triangle_data[i1].add_point<1>(inter_pts[1]); + if (!all_triangle_data[i1].on_same_edge(src_id, tgt_id)) + { + all_triangle_data[i1].segments.emplace_back(src_id, tgt_id); + all_triangle_data[i1].segment_input_triangle_ids.push_back(i2); + } + // + src_id=all_triangle_data[i2].add_point<2>(inter_pts[0]); + tgt_id=all_triangle_data[i2].add_point<2>(inter_pts[1]); + if (!all_triangle_data[i2].on_same_edge(src_id, tgt_id)) + { + all_triangle_data[i2].segments.emplace_back(src_id, tgt_id); + all_triangle_data[i2].segment_input_triangle_ids.push_back(i1); + } + } break; default: + { + std::vector ipt_ids1(nbi+1), ipt_ids2(nbi+1); + for (std::size_t i=0;i(inter_pts[i]); + ipt_ids2[i]=all_triangle_data[i2].add_point<2>(inter_pts[i]); } + ipt_ids1.back()=ipt_ids1.front(); + ipt_ids2.back()=ipt_ids2.front(); + + for (std::size_t i=0;i>> all_segments_ids(all_segments.size()); - auto deduplicate_inserted_segments = [&](std::size_t ti) { - if (!all_segments[ti].empty()) + if (!all_triangle_data[ti].segments.empty()) { - std::map point_id_map; + std::vector& points=all_triangle_data[ti].points; + std::vector>& segments=all_triangle_data[ti].segments; + std::vector indices(points.size()-3); + std::iota(indices.begin(), indices.end(),3); - auto get_point_id = [&](const EK::Point_3& pt) - { - auto insert_res = point_id_map.insert(std::make_pair(pt, all_points[ti].size())); - if (insert_res.second) - all_points[ti].push_back(pt); - return insert_res.first->second; - }; + std::sort(indices.begin(), indices.end(), [&points](std::size_t i, std::size_t j) + { return points[i] id_map(points.size()); + id_map[0]=0; + id_map[1]=1; + id_map[2]=2; + std::vector unique_ids; + unique_ids.reserve(indices.size()); - if (!all_points[ti].empty()) + //make points unique + create mapping between indices + for (std::size_t i=0; i tmp; - tmp.swap(all_points[ti]); - for (const EPoint_3& pt : tmp) - get_point_id(pt); + std::size_t new_id=unique_ids.size()+3; + unique_ids.push_back(indices[i]); + id_map[indices[i]]=new_id; + while(i+1!=indices.size() && points[indices[i]]==points[indices[i+1]]) + { + id_map[indices[++i]]=new_id; + } + } + + CGAL_assertion(points.size() == all_triangle_data[ti].point_locations.size()); + if (unique_ids.size()==indices.size()) + // TODO: do we want to keep points sorted? if yes always swap the 3 containers + return; // no duplicates + + // now make points unique + using EPoint_3 = EK::Point_3; // workaround for MSVC 2022 bug + std::vector tmp; + std::vector tmp_locations; + tmp.reserve(unique_ids.size()+3); + tmp_locations.reserve(unique_ids.size()+3); + tmp.push_back(points[0]); + tmp.push_back(points[1]); + tmp.push_back(points[2]); + tmp_locations.push_back(-1); + tmp_locations.push_back(-2); + tmp_locations.push_back(-3); + for(std::size_t i : unique_ids) + { + tmp.push_back(points[i]); + tmp_locations.push_back(all_triangle_data[ti].point_locations[i]); } + tmp.swap(points); + tmp_locations.swap(all_triangle_data[ti].point_locations); + CGAL_assertion(points.size() == all_triangle_data[ti].point_locations.size()); - std::size_t nbs = all_segments[ti].size(); - std::vector> filtered_segments; + // now make segments unique + std::size_t nbs = segments.size(); + std::vector> filtered_segments; std::vector filtered_in_triangle_ids; filtered_segments.reserve(nbs); + filtered_in_triangle_ids.reserve(nbs); std::set> segset; for (std::size_t si=0; sisegments[si].first); + CGAL_assertion(id_map.size()>segments[si].second); + std::pair seg_ids = + CGAL::make_sorted_pair(id_map[segments[si].first], id_map[segments[si].second]); + if (segset.insert(seg_ids).second) { - all_segments_ids[ti].emplace_back(src_id, tgt_id); - filtered_in_triangle_ids.push_back(all_in_triangle_ids[ti][si]); + filtered_segments.push_back(seg_ids); + filtered_in_triangle_ids.push_back(all_triangle_data[ti].segment_input_triangle_ids[si]); } } - if (all_segments_ids[ti].size()!=nbs) - filtered_in_triangle_ids.swap(all_in_triangle_ids[ti]); + filtered_in_triangle_ids.swap(all_triangle_data[ti].segment_input_triangle_ids); + filtered_segments.swap(segments); + + CGAL_assertion(points.size()-3==std::set(std::next(points.begin(),3), points.end()).size()); + CGAL_assertion(points.size()==std::set(points.begin(), points.end()).size()); } }; @@ -1062,11 +1350,28 @@ void autorefine_triangle_soup(PointRange& soup_points, auto refine_triangles = [&](std::size_t ti) { - if (all_segments[ti].empty() && all_points[ti].empty()) + if (all_triangle_data[ti].points.size()==3) new_triangles.push_back({triangles[ti], ti}); else { - autorefine_impl::generate_subtriangles(ti, all_segments_ids[ti], all_points[ti], all_in_triangle_ids[ti], intersecting_triangles, coplanar_triangles, triangles, new_triangles); +#ifdef CGAL_AUTOREF_USE_FIXED_PROJECTION_TRAITS + const std::array& t = triangles[ti]; + typename EK::Vector_3 orth = CGAL::normal(t[0], t[1], t[2]); // TODO::avoid construction? + int c = CGAL::abs(orth[0]) > CGAL::abs(orth[1]) ? 0 : 1; + c = CGAL::abs(orth[2]) > CGAL::abs(orth[c]) ? 2 : c; + + if (c == 0) { + autorefine_impl::generate_subtriangles(ti, all_triangle_data[ti], intersecting_triangles, coplanar_triangles, triangles, new_triangles); + } + else if (c == 1) { + autorefine_impl::generate_subtriangles(ti, all_triangle_data[ti], intersecting_triangles, coplanar_triangles, triangles, new_triangles); + } + else if (c == 2) { + autorefine_impl::generate_subtriangles(ti, all_triangle_data[ti], intersecting_triangles, coplanar_triangles, triangles, new_triangles); + } +#else + autorefine_impl::generate_subtriangles(ti, all_triangle_data[ti], intersecting_triangles, coplanar_triangles, triangles, new_triangles); +#endif } #ifdef CGAL_AUTOREF_USE_PROGRESS_DISPLAY diff --git a/Triangulation_2/include/CGAL/Constrained_Delaunay_triangulation_2.h b/Triangulation_2/include/CGAL/Constrained_Delaunay_triangulation_2.h index 52c14ea54b78..714e9295b61f 100644 --- a/Triangulation_2/include/CGAL/Constrained_Delaunay_triangulation_2.h +++ b/Triangulation_2/include/CGAL/Constrained_Delaunay_triangulation_2.h @@ -342,6 +342,25 @@ class Constrained_Delaunay_triangulation_2 return number_of_vertices() - n; } + // function usable only if no constraint has been inserted + void restore_Delaunay(Vertex_handle v) + { + if(this->dimension() <= 1) return; + + Face_handle f=v->face(); + Face_handle next; + int i; + Face_handle start(f); + do { + i = f->index(v); + next = f->neighbor(ccw(i)); // turn ccw around v + propagating_flip(f,i); + f = next; + } while(next != start); + return; + } + + #ifndef CGAL_TRIANGULATION_2_DONT_INSERT_RANGE_OF_POINTS_WITH_INFO private: //top stands for tuple-or-pair diff --git a/Triangulation_2/include/CGAL/Triangulation_2.h b/Triangulation_2/include/CGAL/Triangulation_2.h index 14e8604640f1..c0f2c7ba50e0 100644 --- a/Triangulation_2/include/CGAL/Triangulation_2.h +++ b/Triangulation_2/include/CGAL/Triangulation_2.h @@ -35,6 +35,7 @@ #include #include #include +#include #include @@ -434,6 +435,81 @@ class Triangulation_2 bool has_inexact_negative_orientation(const Point &p, const Point &q, const Point &r) const; + + +template +inline +bool +projection_traits_has_inexact_negative_orientation(const Point &p, const Point &q, const Point &r, + const T& ) const +{ + // So that this code works well with Lazy_kernel + internal::Static_filters_predicates::Get_approx get_approx; + + const double px = to_double(get_approx(p).x()); + const double py = to_double(get_approx(p).y()); + const double qx = to_double(get_approx(q).x()); + const double qy = to_double(get_approx(q).y()); + const double rx = to_double(get_approx(r).x()); + const double ry = to_double(get_approx(r).y()); + + const double pqx = qx - px; + const double pqy = qy - py; + const double prx = rx - px; + const double pry = ry - py; + + return ( determinant(pqx, pqy, prx, pry) < 0); +} + +template +inline +bool +projection_traits_has_inexact_negative_orientation(const Point &p, const Point &q, const Point &r, + const ::CGAL::internal::Projection_traits_3& ) const +{ // So that this code works well with Lazy_kernel + internal::Static_filters_predicates::Get_approx get_approx; + + const double px = to_double(get_approx(p).y()); + const double py = to_double(get_approx(p).z()); + const double qx = to_double(get_approx(q).y()); + const double qy = to_double(get_approx(q).z()); + const double rx = to_double(get_approx(r).y()); + const double ry = to_double(get_approx(r).z()); + + const double pqx = qx - px; + const double pqy = qy - py; + const double prx = rx - px; + const double pry = ry - py; + + return ( determinant(pqx, pqy, prx, pry) < 0); +} + + +template +inline +bool +projection_traits_has_inexact_negative_orientation(const Point &p, const Point &q, const Point &r, + const ::CGAL::internal::Projection_traits_3& ) const +{ // So that this code works well with Lazy_kernel + internal::Static_filters_predicates::Get_approx get_approx; + + const double px = to_double(get_approx(p).x()); + const double py = to_double(get_approx(p).z()); + const double qx = to_double(get_approx(q).x()); + const double qy = to_double(get_approx(q).z()); + const double rx = to_double(get_approx(r).x()); + const double ry = to_double(get_approx(r).z()); + + const double pqx = qx - px; + const double pqy = qy - py; + const double prx = rx - px; + const double pry = ry - py; + + return ( determinant(pqx, pqy, prx, pry) < 0); +} + + + public: Face_handle inexact_locate(const Point& p, @@ -3152,6 +3228,7 @@ inexact_locate(const Point & t, Face_handle start, int n_of_turns) const return c; } + template inline bool @@ -3159,22 +3236,8 @@ Triangulation_2:: has_inexact_negative_orientation(const Point &p, const Point &q, const Point &r) const { - // So that this code works well with Lazy_kernel - internal::Static_filters_predicates::Get_approx get_approx; - - const double px = to_double(get_approx(p).x()); - const double py = to_double(get_approx(p).y()); - const double qx = to_double(get_approx(q).x()); - const double qy = to_double(get_approx(q).y()); - const double rx = to_double(get_approx(r).x()); - const double ry = to_double(get_approx(r).y()); - - const double pqx = qx - px; - const double pqy = qy - py; - const double prx = rx - px; - const double pry = ry - py; - - return ( determinant(pqx, pqy, prx, pry) < 0); + Gt gt; + return projection_traits_has_inexact_negative_orientation(p,q,r,gt); } #endif