From 8e0b41bf7e81a9cb22ddd2eb00308bc643d15086 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Fri, 11 Mar 2022 08:27:03 +0100 Subject: [PATCH 1/3] add tunable angle to Angle_tester for polylines_to_protect() --- .../CGAL/Mesh_3/polylines_to_protect.h | 38 +++++++++++++------ 1 file changed, 27 insertions(+), 11 deletions(-) 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..cd103170224 100644 --- a/Mesh_3/include/CGAL/Mesh_3/polylines_to_protect.h +++ b/Mesh_3/include/CGAL/Mesh_3/polylines_to_protect.h @@ -301,6 +301,14 @@ struct Polyline_visitor template struct Angle_tester { + const double m_angle_rad;//in radian to avoid extra computations + const bool m_test_angle; + + Angle_tester(const double angle_deg)//angle given in degrees for readability + : m_angle_rad(angle_deg * CGAL_PI /180.)//stored in radian + , m_test_angle(angle_deg != 90.)//under 90 we check for acute, with exact computation + {} + template bool operator()(vertex_descriptor& v, const Graph& g) const { @@ -320,15 +328,21 @@ 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; - return true; + if (!m_test_angle) + { + if (CGAL::angle(p1, p, p2) == CGAL::ACUTE) + return true; + } + else + { + const typename Kernel::Vector_3 e1 = p1 - p; + const typename Kernel::Vector_3 e2 = p2 - p; + const double a_rad = std::acos(e1 * e2 + / CGAL::sqrt(e1*e1) + / CGAL::sqrt(e2*e2)); + // std::cerr << "At point " << p << ": the angle is " << a << std::endl; + if(a_rad < m_angle_rad) + return true; } } return false; @@ -997,7 +1011,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 +1050,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 From cfd0fd220175950b95b6d5cff500c2c25bbfb52f Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Fri, 11 Mar 2022 10:40:14 +0100 Subject: [PATCH 2/3] add missing default constructor --- Mesh_3/include/CGAL/Mesh_3/polylines_to_protect.h | 5 +++++ 1 file changed, 5 insertions(+) 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 cd103170224..80222bf0306 100644 --- a/Mesh_3/include/CGAL/Mesh_3/polylines_to_protect.h +++ b/Mesh_3/include/CGAL/Mesh_3/polylines_to_protect.h @@ -304,6 +304,11 @@ struct Angle_tester const double m_angle_rad;//in radian to avoid extra computations const bool m_test_angle; + Angle_tester() + : m_angle_rad(0.5 * CGAL_PI)//stored in radian + , m_test_angle(false)//under 90 we check for acute, with exact computation + {} + Angle_tester(const double angle_deg)//angle given in degrees for readability : m_angle_rad(angle_deg * CGAL_PI /180.)//stored in radian , m_test_angle(angle_deg != 90.)//under 90 we check for acute, with exact computation From f89a6e4ce86738351cd74fd98d7c337bfe857761 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Tue, 22 Mar 2022 15:25:43 +0100 Subject: [PATCH 3/3] test if angle is acute, and then if smaller than the input angle if angle at p is acute, then v should be considered as a terminal vertex to ensure termination --- .../CGAL/Mesh_3/polylines_to_protect.h | 35 ++++++++++--------- 1 file changed, 18 insertions(+), 17 deletions(-) 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 80222bf0306..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,17 +304,14 @@ struct Polyline_visitor template struct Angle_tester { - const double m_angle_rad;//in radian to avoid extra computations - const bool m_test_angle; + const double m_angle_sq_cosine;// squared cosine of `std:min(90, angle_deg)` Angle_tester() - : m_angle_rad(0.5 * CGAL_PI)//stored in radian - , m_test_angle(false)//under 90 we check for acute, with exact computation + : m_angle_sq_cosine(0) {} Angle_tester(const double angle_deg)//angle given in degrees for readability - : m_angle_rad(angle_deg * CGAL_PI /180.)//stored in radian - , m_test_angle(angle_deg != 90.)//under 90 we check for acute, with exact computation + : m_angle_sq_cosine(CGAL::square(std::cos((std::min)(90.,angle_deg) * CGAL_PI / 180.))) {} template @@ -333,20 +333,21 @@ struct Angle_tester const typename Kernel::Point_3& p1 = g[v1].point; const typename Kernel::Point_3& p2 = g[v2].point; - if (!m_test_angle) - { - if (CGAL::angle(p1, p, p2) == CGAL::ACUTE) - return true; - } - else + //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 double a_rad = std::acos(e1 * e2 - / CGAL::sqrt(e1*e1) - / CGAL::sqrt(e2*e2)); - // std::cerr << "At point " << p << ": the angle is " << a << std::endl; - if(a_rad < m_angle_rad) + + 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; } }