diff --git a/Mesh_3/include/CGAL/Mesh_3/polylines_to_protect.h b/Mesh_3/include/CGAL/Mesh_3/polylines_to_protect.h index 19ba3f8d87d..cdbcbb956e5 100644 --- a/Mesh_3/include/CGAL/Mesh_3/polylines_to_protect.h +++ b/Mesh_3/include/CGAL/Mesh_3/polylines_to_protect.h @@ -19,12 +19,15 @@ #include #include #include // std::swap +#include // std::min + #include #include #include #include #include #include // for CGAL::Null_subdomain_index +#include #include // for boost::prior #include @@ -301,6 +304,16 @@ struct Polyline_visitor template struct Angle_tester { + const double m_angle_sq_cosine;// squared cosine of `std:min(90, angle_deg)` + + Angle_tester() + : m_angle_sq_cosine(0) + {} + + Angle_tester(const double angle_deg)//angle given in degrees for readability + : m_angle_sq_cosine(CGAL::square(std::cos((std::min)(90.,angle_deg) * CGAL_PI / 180.))) + {} + template bool operator()(vertex_descriptor& v, const Graph& g) const { @@ -320,15 +333,22 @@ struct Angle_tester const typename Kernel::Point_3& p1 = g[v1].point; const typename Kernel::Point_3& p2 = g[v2].point; - if(CGAL::angle(p1, p, p2) == CGAL::ACUTE) { - // const typename Kernel::Vector_3 e1 = p1 - p; - // const typename Kernel::Vector_3 e2 = p2 - p; - // std::cerr << "At point " << p << ": the angle is " - // << ( std::acos(e1 * e2 - // / CGAL::sqrt(e1*e1) - // / CGAL::sqrt(e2*e2)) - // * 180 / CGAL_PI ) << std::endl; + //if angle at v is acute, v must be considered as a terminal vertex + // to ensure termination + if (CGAL::angle(p1, p, p2) == CGAL::ACUTE) return true; + else if (m_angle_sq_cosine > 0.)//check angle only if angle is > 90. + { + const typename Kernel::Vector_3 e1 = p1 - p; + const typename Kernel::Vector_3 e2 = p2 - p; + + const auto scalar_product = e1 * e2; + if (CGAL::is_positive(scalar_product)) + return true; + + const auto sq_scalar_product = CGAL::square(scalar_product); + if (sq_scalar_product <= m_angle_sq_cosine * (e1 * e1) * (e2 * e2)) + return true; } } return false; @@ -997,7 +1017,8 @@ template >& polylines, PolylineInputIterator existing_polylines_begin, - PolylineInputIterator existing_polylines_end) + PolylineInputIterator existing_polylines_end, + const double& angle = 90.)//when not provided, check only for acute angles { typedef P Point_3; typedef typename Kernel_traits

::Kernel K; @@ -1035,8 +1056,9 @@ polylines_to_protect(std::vector >& polylines, Less_for_Graph_vertex_descriptors less(graph); const Graph& const_graph = graph; typedef typename Kernel_traits

::Kernel K; + Mesh_3::Angle_tester angle_tester(angle); split_graph_into_polylines(const_graph, visitor, - Mesh_3::Angle_tester(), less); + angle_tester, less); } template