From ab9c6511bd857e9064d9d8aa4e9604d98cc2f2a3 Mon Sep 17 00:00:00 2001 From: Grzegorz Bokota Date: Wed, 12 Mar 2025 18:07:52 +0100 Subject: [PATCH 1/9] add tests, start fixing --- .../triangulation/intersection.hpp | 14 ++++-- .../triangulation/point.hpp | 10 ++++ .../triangulation/triangulate.hpp | 18 ++++++- src/tests/test_triangulate.py | 48 +++++++++++++++++++ 4 files changed, 85 insertions(+), 5 deletions(-) diff --git a/src/PartSegCore_compiled_backend/triangulation/intersection.hpp b/src/PartSegCore_compiled_backend/triangulation/intersection.hpp index dbbc245..5878f5c 100644 --- a/src/PartSegCore_compiled_backend/triangulation/intersection.hpp +++ b/src/PartSegCore_compiled_backend/triangulation/intersection.hpp @@ -112,6 +112,12 @@ inline int64_t double_as_hex(double d) { return result; } +enum Orientation { + COLLINEAR = 0, + CLOCKWISE = 1, + COUNTERCLOCKWISE = 2, +}; + /** * Determines the orientation of the triplet (p, q, r). * @@ -123,8 +129,8 @@ inline int64_t double_as_hex(double d) { * 1 if the triplet (p, q, r) is in a clockwise orientation. * 2 if the triplet (p, q, r) is in a counterclockwise orientation. */ -inline int _orientation(const point::Point &p, const point::Point &q, - const point::Point &r) { +inline Orientation _orientation(const point::Point &p, const point::Point &q, + const point::Point &r) { double val1 = ((q.y - p.y) * (r.x - q.x)); double val2 = ((r.y - q.y) * (q.x - p.x)); // This commented code if for debugging purposes of differences between macOS @@ -142,8 +148,8 @@ inline int _orientation(const point::Point &p, const point::Point &q, // } // Instead of using classical equation, we need to use two variables // to handle problem with strange behavior on macOS. - if (val1 == val2) return 0; - return (val1 > val2) ? 1 : 2; + if (val1 == val2) return Orientation::COLLINEAR; + return (val1 > val2) ? Orientation::CLOCKWISE : Orientation::COUNTERCLOCKWISE; } /** diff --git a/src/PartSegCore_compiled_backend/triangulation/point.hpp b/src/PartSegCore_compiled_backend/triangulation/point.hpp index 61f8143..a605d26 100644 --- a/src/PartSegCore_compiled_backend/triangulation/point.hpp +++ b/src/PartSegCore_compiled_backend/triangulation/point.hpp @@ -208,6 +208,16 @@ struct Segment { } }; }; +Point centroid(const std::vector &point_list) { + Point res(0, 0); + for (auto &point : point_list) { + res.x += point.x; + res.y += point.y; + } + res.x /= float(point_list.size()); + res.y /= float(point_list.size()); + return res; +} } // namespace partsegcore::point // overload of hash function for diff --git a/src/PartSegCore_compiled_backend/triangulation/triangulate.hpp b/src/PartSegCore_compiled_backend/triangulation/triangulate.hpp index e362d78..57f4a54 100644 --- a/src/PartSegCore_compiled_backend/triangulation/triangulate.hpp +++ b/src/PartSegCore_compiled_backend/triangulation/triangulate.hpp @@ -2,6 +2,7 @@ #define PARTSEGCORE_TRIANGULATE_H #include +#include #include #include // memory header is required on linux, and not on macos #include @@ -616,10 +617,25 @@ struct MonotonePolygonBuilder { */ inline bool _is_convex(const std::vector &polygon) { int orientation = 0; - int triangle_orientation; + intersection::Orientation triangle_orientation; + point::Point centroid = point::centroid(polygon); + double start_angle = + std::atan2(polygon[0].y - centroid.y, polygon[0].x - centroid.x); + double prev_angle = 0; for (std::size_t i = 0; i < polygon.size() - 2; i++) { triangle_orientation = intersection::_orientation(polygon[i], polygon[i + 1], polygon[i + 2]); + double angle = std::atan2(polygon[i + 1].y - centroid.y, + polygon[i + 1].x - centroid.x) - + start_angle; + if (angle < 0) { + angle += 2 * M_PI; + } + if (angle < prev_angle) { + return false; + } else { + prev_angle = angle; + } if (triangle_orientation == 0) continue; if (orientation == 0) orientation = triangle_orientation; diff --git a/src/tests/test_triangulate.py b/src/tests/test_triangulate.py index f5b7286..dde536a 100644 --- a/src/tests/test_triangulate.py +++ b/src/tests/test_triangulate.py @@ -208,6 +208,7 @@ def test_triangle_convex_polygon(): ], ) def test_triangulate_polygon_py_convex(polygon, expected): + assert is_convex(polygon) assert triangulate_polygon_py(polygon)[0] == expected @@ -691,3 +692,50 @@ def test_splitting_edges(split_edges, triangles): ) triangles_ = triangulate_polygon_with_edge_numpy_li([polygon], split_edges=split_edges)[1][2] assert len(triangles_) == triangles + + +def generate_regular_polygon(n, reverse, radius=1): + angles = np.linspace(0, 2 * np.pi, n, endpoint=False) + if reverse: + angles = angles[::-1] + return np.column_stack((radius * np.cos(angles), radius * np.sin(angles))) + + +def pentagram(reverse): + radius = 10 + n = 5 + angles = np.linspace(0, 4 * np.pi, n, endpoint=False) + if reverse: + angles = angles[::-1] + return np.column_stack((radius * np.cos(angles), radius * np.sin(angles))) + + +def rotation_matrix(angle): + return np.array( + [ + [np.cos(np.radians(angle)), -np.sin(np.radians(angle))], + [np.sin(np.radians(angle)), np.cos(np.radians(angle))], + ] + ) + + +ANGLES = [0, 5, 75, 95, 355] + + +@pytest.mark.parametrize('angle', ANGLES, ids=str) +@pytest.mark.parametrize('reverse', [False, True]) +def test_is_convex_self_intersection(angle, reverse): + p = pentagram(reverse) + rot = rotation_matrix(angle) + data = np.dot(p, rot) + assert not is_convex(data) + + +@pytest.mark.parametrize('angle', ANGLES, ids=str) +@pytest.mark.parametrize('n_vertex', [3, 4, 7, 12, 15, 20]) +@pytest.mark.parametrize('reverse', [False, True]) +def test_is_convex2(angle, n_vertex, reverse): + poly = generate_regular_polygon(n_vertex, reverse=reverse) + rot = rotation_matrix(angle) + rotated_poly = np.dot(poly, rot) + assert is_convex(rotated_poly) From fb47059c4785b583b18cee1f3df3064694d5c869 Mon Sep 17 00:00:00 2001 From: Grzegorz Bokota Date: Thu, 13 Mar 2025 00:38:23 +0100 Subject: [PATCH 2/9] fix check for beeing convex --- .../triangulation/triangulate.hpp | 84 +++++++++++++------ 1 file changed, 60 insertions(+), 24 deletions(-) diff --git a/src/PartSegCore_compiled_backend/triangulation/triangulate.hpp b/src/PartSegCore_compiled_backend/triangulation/triangulate.hpp index 57f4a54..afda8d3 100644 --- a/src/PartSegCore_compiled_backend/triangulation/triangulate.hpp +++ b/src/PartSegCore_compiled_backend/triangulation/triangulate.hpp @@ -603,6 +603,36 @@ struct MonotonePolygonBuilder { }; }; +/* + * Check if the polygon, that all angles have the same orientation + * do not have self-intersections + * + * @param begin Iterator to the first point of the polygon + * @param end Iterator to the end of the polygon + * @param centroid Centroid of the polygon + * + * @return True if the polygon is simple, false otherwise + */ +template +bool is_simple_polygon(Iterator begin, Iterator end, point::Point centroid) { + double start_angle = std::atan2(begin->y - centroid.y, begin->x - centroid.x); + double prev_angle = 0; + begin++; + for (; begin != end; begin++) { + double angle = + std::atan2(begin->y - centroid.y, begin->x - centroid.x) - start_angle; + if (angle < 0) { + angle += 2 * M_PI; + } + if (angle < prev_angle) { + return false; + } else { + prev_angle = angle; + } + } + return true; +} + /** * Checks if a given polygon is convex. * @@ -616,41 +646,47 @@ struct MonotonePolygonBuilder { * @return True if the polygon is convex, false otherwise. */ inline bool _is_convex(const std::vector &polygon) { - int orientation = 0; + if (polygon.size() < 3) return false; + if (polygon.size() == 3) return true; + intersection::Orientation orientation = intersection::Orientation::COLLINEAR; intersection::Orientation triangle_orientation; - point::Point centroid = point::centroid(polygon); - double start_angle = - std::atan2(polygon[0].y - centroid.y, polygon[0].x - centroid.x); - double prev_angle = 0; - for (std::size_t i = 0; i < polygon.size() - 2; i++) { - triangle_orientation = - intersection::_orientation(polygon[i], polygon[i + 1], polygon[i + 2]); - double angle = std::atan2(polygon[i + 1].y - centroid.y, - polygon[i + 1].x - centroid.x) - - start_angle; - if (angle < 0) { - angle += 2 * M_PI; + size_t idx = 0; + for (; idx < polygon.size() - 2; idx++) { + triangle_orientation = intersection::_orientation( + polygon[idx], polygon[(idx + 1)], polygon[(idx + 2)]); + if (triangle_orientation != intersection::Orientation::COLLINEAR) { + orientation = triangle_orientation; + break; } - if (angle < prev_angle) { + } + if (orientation == intersection::Orientation::COLLINEAR) { + return false; + } + for (; idx < polygon.size() - 2; idx++) { + triangle_orientation = intersection::_orientation( + polygon[idx], polygon[(idx + 1)], polygon[(idx + 2)]); + if (triangle_orientation != 0 && triangle_orientation != orientation) { return false; - } else { - prev_angle = angle; } - if (triangle_orientation == 0) continue; - if (orientation == 0) - orientation = triangle_orientation; - else if (orientation != triangle_orientation) - return false; } triangle_orientation = intersection::_orientation( polygon[polygon.size() - 2], polygon[polygon.size() - 1], polygon[0]); - if (triangle_orientation != 0 && triangle_orientation != orientation) + if (triangle_orientation != 0 && triangle_orientation != orientation) { return false; + } triangle_orientation = intersection::_orientation(polygon[polygon.size() - 1], polygon[0], polygon[1]); - if (triangle_orientation != 0 && triangle_orientation != orientation) + if (triangle_orientation != 0 && triangle_orientation != orientation) { return false; - return true; + } + + point::Point centroid = point::centroid(polygon); + + if (orientation == intersection::Orientation::COUNTERCLOCKWISE) { + return is_simple_polygon(polygon.begin(), polygon.end(), centroid); + } else { + return is_simple_polygon(polygon.rbegin(), polygon.rend(), centroid); + } } /** From 9957feb5e276211aa310a43c4ba21a1b33f8b6d5 Mon Sep 17 00:00:00 2001 From: Grzegorz Bokota Date: Thu, 13 Mar 2025 08:11:36 +0100 Subject: [PATCH 3/9] improve centroid --- src/PartSegCore_compiled_backend/triangulation/point.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/PartSegCore_compiled_backend/triangulation/point.hpp b/src/PartSegCore_compiled_backend/triangulation/point.hpp index a605d26..b04fc24 100644 --- a/src/PartSegCore_compiled_backend/triangulation/point.hpp +++ b/src/PartSegCore_compiled_backend/triangulation/point.hpp @@ -209,6 +209,9 @@ struct Segment { }; }; Point centroid(const std::vector &point_list) { + if (point_list.empty()) { + return {0, 0}; + } Point res(0, 0); for (auto &point : point_list) { res.x += point.x; From 82eed0ac386fdd8fb86eba6c1ee2aa91b3444249 Mon Sep 17 00:00:00 2001 From: Grzegorz Bokota Date: Thu, 13 Mar 2025 08:13:25 +0100 Subject: [PATCH 4/9] improve docstring --- .../triangulation/intersection.hpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/PartSegCore_compiled_backend/triangulation/intersection.hpp b/src/PartSegCore_compiled_backend/triangulation/intersection.hpp index 5878f5c..0c0beda 100644 --- a/src/PartSegCore_compiled_backend/triangulation/intersection.hpp +++ b/src/PartSegCore_compiled_backend/triangulation/intersection.hpp @@ -125,9 +125,11 @@ enum Orientation { * @param q The second point. * @param r The third point. * - * @return 0 if p, q and r are collinear. - * 1 if the triplet (p, q, r) is in a clockwise orientation. - * 2 if the triplet (p, q, r) is in a counterclockwise orientation. + * @return Orientation::COLLINEAR if p, q and r are collinear. + * Orientation::CLOCKWISE if the triplet (p, q, r) is in a clockwise + * orientation. + * Orientation::COUNTERCLOCKWISE if the triplet (p, q, r) is in a + * counterclockwise orientation. */ inline Orientation _orientation(const point::Point &p, const point::Point &q, const point::Point &r) { From f11a871985b39b51a937d36f50525229ecfc0652 Mon Sep 17 00:00:00 2001 From: Grzegorz Bokota Date: Thu, 13 Mar 2025 08:25:23 +0100 Subject: [PATCH 5/9] fix orientation usage --- .../triangulation/intersection.hpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/PartSegCore_compiled_backend/triangulation/intersection.hpp b/src/PartSegCore_compiled_backend/triangulation/intersection.hpp index 0c0beda..6ce8e72 100644 --- a/src/PartSegCore_compiled_backend/triangulation/intersection.hpp +++ b/src/PartSegCore_compiled_backend/triangulation/intersection.hpp @@ -171,17 +171,17 @@ inline bool _do_intersect(const point::Segment &s1, const point::Segment &s2) { const point::Point &p2 = s2.bottom; const point::Point &q2 = s2.top; - int o1 = _orientation(p1, q1, p2); - int o2 = _orientation(p1, q1, q2); - int o3 = _orientation(p2, q2, p1); - int o4 = _orientation(p2, q2, q1); + Orientation o1 = _orientation(p1, q1, p2); + Orientation o2 = _orientation(p1, q1, q2); + Orientation o3 = _orientation(p2, q2, p1); + Orientation o4 = _orientation(p2, q2, q1); if (o1 != o2 && o3 != o4) return true; - if (o1 == 0 && _on_segment(p1, p2, q1)) return true; - if (o2 == 0 && _on_segment(p1, q2, q1)) return true; - if (o3 == 0 && _on_segment(p2, p1, q2)) return true; - if (o4 == 0 && _on_segment(p2, q1, q2)) return true; + if (o1 == Orientation::COLLINEAR && _on_segment(p1, p2, q1)) return true; + if (o2 == Orientation::COLLINEAR && _on_segment(p1, q2, q1)) return true; + if (o3 == Orientation::COLLINEAR && _on_segment(p2, p1, q2)) return true; + if (o4 == Orientation::COLLINEAR && _on_segment(p2, q1, q2)) return true; return false; } From e0695915ed4ffaebd9c61eb397c541ebb156fae7 Mon Sep 17 00:00:00 2001 From: Grzegorz Bokota Date: Thu, 13 Mar 2025 17:04:26 +0100 Subject: [PATCH 6/9] improve tests --- src/tests/test_triangulate.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/tests/test_triangulate.py b/src/tests/test_triangulate.py index dde536a..f6e754c 100644 --- a/src/tests/test_triangulate.py +++ b/src/tests/test_triangulate.py @@ -694,16 +694,19 @@ def test_splitting_edges(split_edges, triangles): assert len(triangles_) == triangles -def generate_regular_polygon(n, reverse, radius=1): +def generate_regular_polygon(n: int, reverse: bool, radius: int = 1) -> np.ndarray: angles = np.linspace(0, 2 * np.pi, n, endpoint=False) if reverse: angles = angles[::-1] return np.column_stack((radius * np.cos(angles), radius * np.sin(angles))) -def pentagram(reverse): - radius = 10 - n = 5 +def generate_self_intersecting_polygon(n: int, reverse: bool, radius: int = 1) -> np.ndarray: + """Generate self-intersecting polygon with n vertices + + The polygon is generated by doubling the angle range of a regular polygon + """ + assert n % 2 == 1 angles = np.linspace(0, 4 * np.pi, n, endpoint=False) if reverse: angles = angles[::-1] @@ -723,9 +726,10 @@ def rotation_matrix(angle): @pytest.mark.parametrize('angle', ANGLES, ids=str) +@pytest.mark.parametrize('n_vertex', [5, 7, 19]) @pytest.mark.parametrize('reverse', [False, True]) -def test_is_convex_self_intersection(angle, reverse): - p = pentagram(reverse) +def test_is_convex_self_intersection(angle, n_vertex, reverse): + p = generate_self_intersecting_polygon(n_vertex, reverse) rot = rotation_matrix(angle) data = np.dot(p, rot) assert not is_convex(data) From ff7b775a7e7f5d89e6058a3713936069904d231d Mon Sep 17 00:00:00 2001 From: Grzegorz Bokota Date: Thu, 13 Mar 2025 17:15:14 +0100 Subject: [PATCH 7/9] improve naming and type annotation --- src/tests/test_triangulate.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/tests/test_triangulate.py b/src/tests/test_triangulate.py index f6e754c..71a78f9 100644 --- a/src/tests/test_triangulate.py +++ b/src/tests/test_triangulate.py @@ -706,14 +706,14 @@ def generate_self_intersecting_polygon(n: int, reverse: bool, radius: int = 1) - The polygon is generated by doubling the angle range of a regular polygon """ - assert n % 2 == 1 + assert n % 2 == 1, 'need odd number to generate self-intersecting polygon' angles = np.linspace(0, 4 * np.pi, n, endpoint=False) if reverse: angles = angles[::-1] return np.column_stack((radius * np.cos(angles), radius * np.sin(angles))) -def rotation_matrix(angle): +def rotation_matrix(angle: float) -> np.ndarray: return np.array( [ [np.cos(np.radians(angle)), -np.sin(np.radians(angle))], From 05bb23789137e4642899a010d051251959ff2233 Mon Sep 17 00:00:00 2001 From: Grzegorz Bokota Date: Thu, 13 Mar 2025 17:20:54 +0100 Subject: [PATCH 8/9] more improvements --- src/tests/test_triangulate.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/tests/test_triangulate.py b/src/tests/test_triangulate.py index 71a78f9..4bcd468 100644 --- a/src/tests/test_triangulate.py +++ b/src/tests/test_triangulate.py @@ -714,6 +714,7 @@ def generate_self_intersecting_polygon(n: int, reverse: bool, radius: int = 1) - def rotation_matrix(angle: float) -> np.ndarray: + """Create a 2D rotation matrix for the given angle in degrees.""" return np.array( [ [np.cos(np.radians(angle)), -np.sin(np.radians(angle))], @@ -738,7 +739,7 @@ def test_is_convex_self_intersection(angle, n_vertex, reverse): @pytest.mark.parametrize('angle', ANGLES, ids=str) @pytest.mark.parametrize('n_vertex', [3, 4, 7, 12, 15, 20]) @pytest.mark.parametrize('reverse', [False, True]) -def test_is_convex2(angle, n_vertex, reverse): +def test_is_convex_regular_polygon(angle, n_vertex, reverse): poly = generate_regular_polygon(n_vertex, reverse=reverse) rot = rotation_matrix(angle) rotated_poly = np.dot(poly, rot) From 627fb154ada21c007e229174eb86242aad6672a5 Mon Sep 17 00:00:00 2001 From: Grzegorz Bokota Date: Thu, 13 Mar 2025 17:48:55 +0100 Subject: [PATCH 9/9] improve assert information --- src/tests/test_triangulate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tests/test_triangulate.py b/src/tests/test_triangulate.py index 4bcd468..18b6a90 100644 --- a/src/tests/test_triangulate.py +++ b/src/tests/test_triangulate.py @@ -706,7 +706,7 @@ def generate_self_intersecting_polygon(n: int, reverse: bool, radius: int = 1) - The polygon is generated by doubling the angle range of a regular polygon """ - assert n % 2 == 1, 'need odd number to generate self-intersecting polygon' + assert n % 2 == 1, 'an odd number is required to generate a self-intersecting polygon' angles = np.linspace(0, 4 * np.pi, n, endpoint=False) if reverse: angles = angles[::-1]