diff --git a/src/PartSegCore_compiled_backend/triangulation/intersection.hpp b/src/PartSegCore_compiled_backend/triangulation/intersection.hpp index dbbc245..6ce8e72 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). * @@ -119,12 +125,14 @@ inline int64_t double_as_hex(double d) { * @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 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 +150,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; } /** @@ -163,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; } diff --git a/src/PartSegCore_compiled_backend/triangulation/point.hpp b/src/PartSegCore_compiled_backend/triangulation/point.hpp index 61f8143..b04fc24 100644 --- a/src/PartSegCore_compiled_backend/triangulation/point.hpp +++ b/src/PartSegCore_compiled_backend/triangulation/point.hpp @@ -208,6 +208,19 @@ 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; + 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..afda8d3 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 @@ -602,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. * @@ -615,26 +646,47 @@ struct MonotonePolygonBuilder { * @return True if the polygon is convex, false otherwise. */ inline bool _is_convex(const std::vector &polygon) { - int orientation = 0; - int triangle_orientation; - for (std::size_t i = 0; i < polygon.size() - 2; i++) { - triangle_orientation = - intersection::_orientation(polygon[i], polygon[i + 1], polygon[i + 2]); - if (triangle_orientation == 0) continue; - if (orientation == 0) + if (polygon.size() < 3) return false; + if (polygon.size() == 3) return true; + intersection::Orientation orientation = intersection::Orientation::COLLINEAR; + intersection::Orientation triangle_orientation; + 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; - else if (orientation != triangle_orientation) + break; + } + } + 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; + } } 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); + } } /** diff --git a/src/tests/test_triangulate.py b/src/tests/test_triangulate.py index f5b7286..18b6a90 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,55 @@ 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: 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 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, '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] + return np.column_stack((radius * np.cos(angles), radius * np.sin(angles))) + + +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))], + [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('n_vertex', [5, 7, 19]) +@pytest.mark.parametrize('reverse', [False, True]) +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) + + +@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_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) + assert is_convex(rotated_poly)