diff --git a/CMakeLists.txt b/CMakeLists.txt index dfde630..7feecf3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,7 +1,10 @@ cmake_minimum_required(VERSION 3.21) project(${SKBUILD_PROJECT_NAME} LANGUAGES CXX) -set(CMAKE_CXX_STANDARD 11) +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_FLAGS_DEBUG_INIT "-DDEBUG -g") +set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -g -DDEBUG") + # Define compiler directive add_definitions(-DNPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION) @@ -103,6 +106,18 @@ python_add_library(_napari_mapping MODULE "${napari_mapping_cython_cxx}" WITH_SO target_include_directories(_napari_mapping PRIVATE "${NUMPY_INCLUDE_DIR}") target_link_libraries(_napari_mapping PRIVATE OpenMP::OpenMP_CXX) +cython_transpile(src/PartSegCore_compiled_backend/triangulate.pyx LANGUAGE CXX OUTPUT_VARIABLE triangulate_cython_cxx) +python_add_library( + triangulate MODULE + "${triangulate_cython_cxx}" + src/PartSegCore_compiled_backend/triangulation/triangulate.hpp + src/PartSegCore_compiled_backend/triangulation/intersection.hpp + src/PartSegCore_compiled_backend/triangulation/point.hpp + src/PartSegCore_compiled_backend/triangulation/debug_util.hpp + WITH_SOABI) +target_include_directories(triangulate PRIVATE "${NUMPY_INCLUDE_DIR}") +target_include_directories(triangulate PRIVATE src/PartSegCore_compiled_backend/) + install(TARGETS euclidean_cython DESTINATION PartSegCore_compiled_backend/sprawl_utils) install(TARGETS path_sprawl_cython DESTINATION PartSegCore_compiled_backend/sprawl_utils) @@ -113,3 +128,4 @@ install(TARGETS calc_bounds DESTINATION PartSegCore_compiled_backend) install(TARGETS mso_bind DESTINATION PartSegCore_compiled_backend/multiscale_opening) install(TARGETS _fast_unique DESTINATION PartSegCore_compiled_backend) install(TARGETS _napari_mapping DESTINATION PartSegCore_compiled_backend) +install(TARGETS triangulate DESTINATION PartSegCore_compiled_backend) diff --git a/pyproject.toml b/pyproject.toml index 29c2436..be35b7c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -9,6 +9,7 @@ sdist.include = ["src/PartSegCore_compiled_backend/version.py"] cmake.version = ">=3.21" sdist.exclude = [".github", "tox.ini", "build_utils", "notebooks", ".readthedocs.yaml"] wheel.exclude = ["**.pyx"] +#cmake.build-type = "Debug" [project] @@ -139,7 +140,8 @@ quote-style = "single" indent-style = "space" docstring-code-format = true line-ending = "lf" -skip-magic-trailing-comma = true +skip-magic-trailing-comma = false + [tool.ruff.lint] select = [ diff --git a/src/PartSegCore_compiled_backend/triangulate.pyx b/src/PartSegCore_compiled_backend/triangulate.pyx new file mode 100644 index 0000000..d804520 --- /dev/null +++ b/src/PartSegCore_compiled_backend/triangulate.pyx @@ -0,0 +1,437 @@ +# distutils: define_macros=CYTHON_TRACE_NOGIL=1 +# distutils: language = c++ +# cython: language_level=3, boundscheck=False, wraparound=False, nonecheck=False, cdivision=True, embedsignature=True + +from collections.abc import Sequence + +import numpy as np +import cython + +cimport numpy as cnp + +from libcpp cimport bool +from libcpp.unordered_set cimport unordered_set +from libcpp.utility cimport pair +from libcpp.vector cimport vector +from libcpp.unordered_map cimport unordered_map +from libcpp.algorithm cimport sort + + +cdef extern from "triangulation/point.hpp" namespace "partsegcore::point": + cdef cppclass Point: + float x + float y + Point() + Point(float x, float y) + bool operator==(const Point& other) const + bool operator!=(const Point& other) const + + cdef cppclass Vector: + float x; + float y; + + cdef cppclass Segment: + Point bottom + Point top + Segment() + Segment(Point bottom, Point top) + + +cdef extern from "triangulation/intersection.hpp" namespace "partsegcore::intersection": + cdef cppclass Event: + float x + float y + int index + bool is_top + Event() + Event(float x, float y, int index, bool is_top) + + cdef cppclass OrderedPair: + int first + int second + OrderedPair() + OrderedPair(int first, int second) + + bool _on_segment(const Point& p, const Point& q, const Point& r) + int _orientation(const Point& p, const Point& q, const Point& r) + bool _do_intersect(const Segment& s1, const Segment& s2) + unordered_set[OrderedPair] _find_intersections(const vector[Segment]& segments) + vector[Point] _find_intersection(const Segment& s1, const Segment& s2) + + +cdef extern from "triangulation/triangulate.hpp" namespace "partsegcore::triangulation": + + cdef cppclass Triangle: + size_t x + size_t y + size_t z + Triangle() + Triangle(int x, int y, int z) + + cdef cppclass MonotonePolygon: + Point top + Point bottom + vector[Point] left + vector[Point] right + MonotonePolygon() + MonotonePolygon(Point top, Point bottom, vector[Point] left, vector[Point] right) + + + cdef cppclass PointTriangle: + Point p1 + Point p2 + Point p3 + + cdef cppclass PathTriangulation: + vector[Triangle] triangles + vector[Point] centers + vector[Vector] offsets + + + bool _is_convex(const vector[Point]& polygon) + vector[Triangle] _triangle_convex_polygon(const vector[Point]& polygon) + bool left_to_right(const Segment& s1, const Segment& s2) + vector[Point] find_intersection_points(const vector[Point]& segments) + vector[PointTriangle] triangulate_monotone_polygon(const MonotonePolygon& polygon) + pair[vector[Triangle], vector[Point]] triangulate_polygon_face(const vector[Point]& polygon) except + nogil + pair[vector[Triangle], vector[Point]] triangulate_polygon_face(const vector[vector[Point]]& polygon_list) except + nogil + PathTriangulation triangulate_path_edge(const vector[Point]& path, bool closed, float limit, bool bevel) except + nogil + + +ctypedef fused float_types: + cnp.float32_t + cnp.float64_t + + +def on_segment(p: Sequence[float], q: Sequence[float], r: Sequence[float]) -> bool: + """ Check if point q is on segment pr + + Parameters + ---------- + p: sequence of 2 floats + beginning of segment + q: sequence of 2 floats + point to check + r: sequence of 2 floats + end of segment + + Returns + ------- + bool: + True if q is on segment pr + """ + return _on_segment( + Point(p[0], p[1]), + Point(q[0], q[1]), + Point(r[0], r[1]) + ) + +def orientation(p: Sequence[float], q: Sequence[float], r: Sequence[float]) -> int: + """ Check orientation of 3 points + + Parameters + ---------- + p: sequence of 2 floats + first point + q: sequence of 2 floats + second point + r: sequence of 2 floats + third point + + Returns + ------- + int: + 0 if points are collinear, 1 if clockwise, 2 if counterclockwise + """ + return _orientation( + Point(p[0], p[1]), + Point(q[0], q[1]), + Point(r[0], r[1]) + ) + + +def do_intersect(s1: Sequence[Sequence[float]], s2: Sequence[Sequence[float]]) -> bool: + """ Check if two segments intersect + + Parameters + ---------- + s1: sequence of 2 sequences of 2 floats + first segment + s2: sequence of 2 sequences of 2 floats + second segment + + Returns + ------- + bool: + True if segments intersect + """ + return _do_intersect( + Segment(Point(s1[0][0], s1[0][1]), Point(s1[1][0], s1[1][1])), + Segment(Point(s2[0][0], s2[0][1]), Point(s2[1][0], s2[1][1])) + ) + + +def find_intersections(segments: Sequence[Sequence[Sequence[float]]]) -> list[tuple[int, int]]: + """ Find intersections between segments""" + cdef vector[Segment] segments_vector + cdef unordered_set[OrderedPair] intersections + cdef OrderedPair p + + segments_vector.reserve(len(segments)) + for segment in segments: + segments_vector.push_back(Segment(Point(segment[0][0], segment[0][1]), Point(segment[1][0], segment[1][1]))) + intersections = _find_intersections(segments_vector) + + return [(p.first, p.second) for p in intersections] + + +def find_intersection_point(s1: Sequence[Sequence[float]], s2: Sequence[Sequence[float]]) -> list[tuple[float, float]]: + """ Find intersection between two segments + + Parameters + ---------- + s1: sequence of 2 sequences of 2 floats + first segment + s2: sequence of 2 sequences of 2 floats + second segment + + Returns + ------- + sequence of 2 floats: + intersection point + """ + cdef vector[Point] p_li = _find_intersection( + Segment(Point(s1[0][0], s1[0][1]), Point(s1[1][0], s1[1][1])), + Segment(Point(s2[0][0], s2[0][1]), Point(s2[1][0], s2[1][1])) + ) + return [(p.x, p.y) for p in p_li] + + +def is_convex(polygon: Sequence[Sequence[float]]) -> bool: + """ Check if polygon is convex""" + cdef vector[Point] polygon_vector + cdef pair[bool, vector[int]] result + + polygon_vector.reserve(len(polygon)) + for point in polygon: + polygon_vector.push_back(Point(point[0], point[1])) + + return _is_convex(polygon_vector) + +def triangle_convex_polygon(polygon: Sequence[Sequence[float]]) -> list[tuple[int, int, int]]: + cdef vector[Point] polygon_vector + cdef vector[Triangle] result + cdef Point p1, p2 + + polygon_vector.reserve(len(polygon)) + polygon_vector.push_back(Point(polygon[0][0], polygon[0][1])) + for point in polygon[1:]: + p1 = polygon_vector[polygon_vector.size() - 1] + p2 = Point(point[0], point[1]) + if p1 != p2: + # prevent from adding polygon edge of width 0 + polygon_vector.push_back(p2) + + result = _triangle_convex_polygon(polygon_vector) + return [(triangle.x, triangle.y, triangle.z) for triangle in result] + + + +def triangulate_polygon_py(polygon: Sequence[Sequence[float]]) -> tuple[list[tuple[int, int, int]], list[tuple[float, float]]]: + """ Triangulate polygon""" + cdef vector[Point] polygon_vector + cdef Point p1, p2 + cdef pair[vector[Triangle], vector[Point]] result + + polygon_vector.reserve(len(polygon)) + polygon_vector.push_back(Point(polygon[0][0], polygon[0][1])) + for point in polygon[1:]: + p1 = polygon_vector[polygon_vector.size() - 1] + p2 = Point(point[0], point[1]) + if p1 != p2: + # prevent from adding polygon edge of width 0 + polygon_vector.push_back(p2) + + result = triangulate_polygon_face(polygon_vector) + return [(triangle.x, triangle.y, triangle.z) for triangle in result.first], [(point.x, point.y) for point in result.second] + + +def triangulate_polygon_numpy(cnp.ndarray[float_types, ndim=2] polygon: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + """ Triangulate polygon""" + cdef vector[Point] polygon_vector + cdef Point p1, p2 + cdef pair[vector[Triangle], vector[Point]] result + + polygon_vector.reserve(polygon.shape[0]) + polygon_vector.push_back(Point(polygon[0, 0], polygon[0, 1])) + for point in polygon[1:]: + p1 = polygon_vector[polygon_vector.size() - 1] + p2 = Point(point[0], point[1]) + if p1 != p2: + # prevent from adding polygon edge of width 0 + polygon_vector.push_back(p2) + + result = triangulate_polygon_face(polygon_vector) + return ( + np.array([(triangle.x, triangle.y, triangle.z) for triangle in result.first], dtype=np.uintp), + np.array([(point.x, point.y) for point in result.second], dtype=np.float32) + ) + + +def triangulate_polygon_numpy_li(list polygon_li: list[np.ndarray]) -> tuple[np.ndarray, np.ndarray]: + """ Triangulate polygon""" + cdef vector[Point] polygon_vector + cdef vector[vector[Point]] polygon_vector_list + cdef Point p1, p2 + cdef pair[vector[Triangle], vector[Point]] result + cdef size_t i; + cdef cnp.ndarray[cnp.float32_t, ndim=2] polygon32 + cdef cnp.ndarray[cnp.float64_t, ndim=2] polygon64 + + polygon_vector_list.reserve(len(polygon_li)) + for polygon in polygon_li: + polygon_vector.clear() + polygon_vector.reserve(polygon.shape[0]) + polygon_vector.push_back(Point(polygon[0, 0], polygon[0, 1])) + + if polygon.dtype == np.float32: + polygon32 = polygon + for i in range(1, polygon.shape[0]): + p1 = polygon_vector.back() + p2 = Point(polygon32[i, 0], polygon32[i, 1]) + if p1 != p2: + # prevent from adding polygon edge of width 0 + polygon_vector.push_back(p2) + else: + polygon64 = polygon + for i in range(1, polygon.shape[0]): + p1 = polygon_vector.back() + p2 = Point(polygon64[i, 0], polygon64[i, 1]) + if p1 != p2: + # prevent from adding polygon edge of width 0 + polygon_vector.push_back(p2) + + if polygon_vector.size() > 1 and polygon_vector.front() == polygon_vector.back(): + polygon_vector.pop_back() + polygon_vector_list.push_back(polygon_vector) + + result = triangulate_polygon_face(polygon_vector_list) + return ( + np.array([(triangle.x, triangle.y, triangle.z) for triangle in result.first], dtype=np.uintp), + np.array([(point.x, point.y) for point in result.second], dtype=np.float32) + ) + +def segment_left_to_right_comparator(s1: Sequence[Sequence[float]], s2: Sequence[Sequence[float]]) -> bool: + """ Compare segments by bottom point""" + return left_to_right( + Segment(Point(s1[0][0], s1[0][1]), Point(s1[1][0], s1[1][1])), + Segment(Point(s2[0][0], s2[0][1]), Point(s2[1][0], s2[1][1])) + ) + + +def find_intersection_points_py(polygon: Sequence[Sequence[float]]) -> Sequence[Sequence[float]]: + """ Find intersection points in polygon""" + cdef vector[Point] polygon_vector + cdef vector[Point] result + + polygon_vector.reserve(len(polygon)) + for point in polygon: + polygon_vector.push_back(Point(point[0], point[1])) + + result = find_intersection_points(polygon_vector) + return [(point.x, point.y) for point in result] + + +def triangulate_monotone_polygon_py(top: Sequence[float], bottom: Sequence[float], left: Sequence[Sequence[float]], right: Sequence[Sequence[float]]) -> Sequence[Sequence[Sequence[float]]]: + """ Triangulate monotone polygon""" + cdef vector[Point] left_vector, right_vector + cdef vector[PointTriangle] result + cdef MonotonePolygon mono_polygon + + left_vector.reserve(len(left)) + for point in left: + left_vector.push_back(Point(point[0], point[1])) + + right_vector.reserve(len(right)) + for point in right: + right_vector.push_back(Point(point[0], point[1])) + + mono_polygon = MonotonePolygon(Point(top[0], top[1]), Point(bottom[0], bottom[1]), left_vector, right_vector) + result = triangulate_monotone_polygon(mono_polygon) + return [ + [(triangle.p1.x, triangle.p1.y), (triangle.p2.x, triangle.p2.y), (triangle.p3.x, triangle.p3.y)] + for triangle in result + ] + + +def triangulate_path_edge_py(path: Sequence[Sequence[float]], closed: bool=False, limit: float=3.0, bevel: bool=False) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """ Triangulate path""" + cdef vector[Point] path_vector + cdef PathTriangulation result + cdef Point p1, p2 + cdef cnp.ndarray[cnp.uint32_t, ndim=2] triangles + + path_vector.reserve(len(path)) + for point in path: + path_vector.push_back(Point(point[0], point[1])) + with cython.nogil: + result = triangulate_path_edge(path_vector, closed, limit, bevel) + + if result.triangles.size() == 0: + triangles = np.zeros((0, 3), dtype=np.uint32) + else: + triangles = np.array([(triangle.x, triangle.y, triangle.z) for triangle in result.triangles], dtype=np.uint32) + + return ( + np.array([(point.x, point.y) for point in result.centers], dtype=np.float32), + np.array([(offset.x, offset.y) for offset in result.offsets], dtype=np.float32), + triangles, + ) + + +def triangulate_polygon_with_edge_numpy_li(polygon_li: list[np.ndarray]) -> tuple[tuple[np.ndarray, np.ndarray], tuple[np.ndarray, np.ndarray, np.ndarray]]: + """ Triangulate polygon""" + cdef vector[Point] polygon_vector + cdef vector[vector[Point]] polygon_vector_list + cdef Point p1, p2 + cdef pair[vector[Triangle], vector[Point]] triangulation_result + cdef vector[PathTriangulation] edge_result + cdef cnp.ndarray[cnp.uint32_t, ndim=2] triangles + + polygon_vector_list.reserve(len(polygon_li)) + for polygon in polygon_li: + polygon_vector.clear() + + polygon_vector.reserve(polygon.shape[0]) + polygon_vector.push_back(Point(polygon[0, 0], polygon[0, 1])) + + for point in polygon[1:]: + p1 = polygon_vector.back() + p2 = Point(point[0], point[1]) + if p1 != p2: + # prevent from adding polygon edge of width 0 + polygon_vector.push_back(p2) + if polygon_vector.size() > 1 and polygon_vector.front() == polygon_vector.back(): + polygon_vector.pop_back() + polygon_vector_list.push_back(polygon_vector) + with cython.nogil: + edge_result.push_back(triangulate_path_edge(polygon_vector, True, 3.0, False)) + + with cython.nogil: + triangulation_result = triangulate_polygon_face(polygon_vector_list) + + if triangulation_result.first.size() == 0: + triangles = np.zeros((0, 3), dtype=np.uint32) + else: + triangles = np.array([(triangle.x, triangle.y, triangle.z) for triangle in triangulation_result.first], dtype=np.uint32) + + return (( + triangles, + np.array([(point.x, point.y) for point in triangulation_result.second], dtype=np.float32) + ), + ( + np.array([(point.x, point.y) for res in edge_result for point in res.centers], dtype=np.float32), + np.array([(offset.x, offset.y) for res in edge_result for offset in res.offsets], dtype=np.float32), + np.array([(triangle.x, triangle.y, triangle.z) for res in edge_result for triangle in res.triangles], dtype=np.uint32), + ) + ) diff --git a/src/PartSegCore_compiled_backend/triangulation/debug_util.hpp b/src/PartSegCore_compiled_backend/triangulation/debug_util.hpp new file mode 100644 index 0000000..2bda53a --- /dev/null +++ b/src/PartSegCore_compiled_backend/triangulation/debug_util.hpp @@ -0,0 +1,60 @@ +// +// Created by Grzegorz Bokota on 24.10.24. +// + +#ifndef PARTSEGCORE_TRIANGULATION_DEBUG_UTIL_ +#define PARTSEGCORE_TRIANGULATION_DEBUG_UTIL_ + +#include + +namespace partsegcore { + +template +void print_set(std::ostream &o, const T &s, const std::string &end = "\n") { + bool first = true; + o << "{"; + for (const auto &el : s) { + if (!first) { + o << ", "; + } + o << el; + first = false; + } + o << "}" << end; +} +template +void print_vector(std::ostream &o, const T &s, const std::string &end = "\n") { + bool first = true; + o << "["; + for (const auto &el : s) { + if (!first) { + o << ", "; + } + o << el; + first = false; + } + o << "]" << end; +} +template +void print_map(std::ostream &o, const T &s, const std::string &end = "\n") { + bool first = true; + o << "{"; + for (const auto &el : s) { + if (!first) { + o << ", "; + } + o << el.first << ": "; + first = false; + print_set(o, el.second, ""); + } + o << "}" << end; +} + +template +std::ostream &operator<<(std::ostream &os, const std::pair &p) { + os << "(" << p.first << ", " << p.second << ")"; + return os; +} +} // namespace partsegcore + +#endif // PARTSEGCORE_TRIANGULATION_DEBUG_UTIL_ diff --git a/src/PartSegCore_compiled_backend/triangulation/intersection.hpp b/src/PartSegCore_compiled_backend/triangulation/intersection.hpp new file mode 100644 index 0000000..dbbc245 --- /dev/null +++ b/src/PartSegCore_compiled_backend/triangulation/intersection.hpp @@ -0,0 +1,364 @@ +// +// Created by Grzegorz Bokota on 11.10.24. +// Based on +// https://www.geeksforgeeks.org/given-a-set-of-line-segments-find-if-any-two-segments-intersect/ +// and +// Mark Berg, Otfried Cheong, Marc Kreveld, Mark Overmars. +// Computational Geometry: Algorithms and Applications. +// 3rd ed., Springer Berlin, Heidelberg, 2008. +// ISBN: 978-3-540-77974-2. +// https://link.springer.com/book/10.1007/978-3-540-77974-2 +// https://doi.org/10.1007/978-3-540-77974-2 +// + +#ifndef PARTSEGCORE_INTERSECTION_H +#define PARTSEGCORE_INTERSECTION_H + +#include +#include +#include +#include +#include + +#include "point.hpp" + +namespace partsegcore::intersection { + +struct Event { + point::Point p; + int index; + bool is_top; + + Event(point::Point::coordinate_t x, point::Point::coordinate_t y, int index, + bool is_left) + : p(x, y), index(index), is_top(is_left) {} + Event(const point::Point &p, int index, bool is_left) + : p(p), index(index), is_top(is_left) {} + Event() = default; + + bool operator<(const Event &e) const { + if (p == e.p) { + if (is_top == e.is_top) { + return index < e.index; + } + return is_top > e.is_top; + } + return p < e.p; + } +}; + +struct EventData { + std::vector tops; + std::vector bottoms; +}; + +struct OrderedPair { + std::size_t first; + std::size_t second; + OrderedPair() = default; + OrderedPair(std::size_t first, std::size_t second) { + if (first < second) { + this->first = first; + this->second = second; + } else { + this->first = second; + this->second = first; + } + } + bool operator==(const OrderedPair &pair) const { + return first == pair.first && second == pair.second; + } +}; +} // namespace partsegcore::intersection + +template <> +struct std::hash { + std::size_t operator()( + const partsegcore::intersection::OrderedPair &pair) const noexcept { + return std::hash()(pair.first) ^ + std::hash()(pair.second); + } +}; // namespace std + +namespace partsegcore::intersection { + +typedef std::map IntersectionEvents; + +/** + * Checks whether point q lies on the line segment defined by points p and r. + * + * This function determines if a given point q lies on the line segment + * that connects points p and r. It checks if the x and y coordinates + * of q are within the bounding box formed by p and r. + * It works, because we use it when orientation is 0, so we know that + * if q is in bounding box, it lies on the segment. + * + * @param p The first endpoint of the line segment. + * @param q The point to check. + * @param r The second endpoint of the line segment. + * @return True if point q lies on the segment pr, false otherwise. + */ +inline bool _on_segment(const point::Point &p, const point::Point &q, + const point::Point &r) { + if (q.x <= std::max(p.x, r.x) && q.x >= std::min(p.x, r.x) && + q.y <= std::max(p.y, r.y) && q.y >= std::min(p.y, r.y)) + return true; + return false; +} + +inline int64_t double_as_hex(double d) { + int64_t result; + std::memcpy(&result, &d, sizeof(result)); + return result; +} + +/** + * Determines the orientation of the triplet (p, q, r). + * + * @param p The first point. + * @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. + */ +inline int _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 + // and linux double val = ((q.y - p.y) * (r.x - q.x)) - ((r.y - q.y) * (q.x - + // p.x)); if (val!= 0 && val1 == val2) { + // std::ostringstream oss; + // // print coordinates and its hex representation + // oss << "val " << val << " val1 " << val1 << " val2 " << val2 << "\n"; + // oss << "p " << p << " p_hex (x=" << std::hex << double_as_hex(p.x) << ", + // y=" << double_as_hex(p.y) << ")\n"; oss << "q " << q << " q_hex (x=" << + // std::hex << double_as_hex(q.x) << ", y=" << double_as_hex(q.y) << ")\n"; + // oss << "r " << r << " r_hex (x=" << std::hex << double_as_hex(r.x) << ", + // y=" << double_as_hex(r.y) << ")\n"; throw std::runtime_error("Orientation + // is not defined for this case. " + oss.str()); + // } + // 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; +} + +/** + * Determines if two line segments intersect. + * + * This function checks whether two given line segments, s1 and s2, intersect. + * It makes use of orientation and on-segment tests to evaluate intersection + * based on the positions of the endpoints of the segments. + * + * @param s1 The first line segment, represented by two endpoints. + * @param s2 The second line segment, represented by two endpoints. + * @return True if the segments intersect, false otherwise. + */ +inline bool _do_intersect(const point::Segment &s1, const point::Segment &s2) { + const point::Point &p1 = s1.bottom; + const point::Point &q1 = s1.top; + 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); + + 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; + + return false; +} + +/** + * @brief Checks if two segments share an endpoint. + * + * This function determines whether two segments, each defined by + * two endpoints, share any endpoint. Specifically, it checks if + * the bottom or top endpoint of the first segment is equal to the + * bottom or top endpoint of the second segment. + * + * @param s1 The first segment. + * @param s2 The second segment. + * @return true if the segments share at least one endpoint, false otherwise. + */ +inline bool _share_endpoint(const point::Segment &s1, + const point::Segment &s2) { + return s1.bottom == s2.bottom || s1.bottom == s2.top || s1.top == s2.bottom || + s1.top == s2.top; +} + +template +typename T::iterator pred(T &s, typename T::iterator it) { + if (it == s.begin()) { + return s.end(); + } + return --it; +} + +template +typename T::iterator succ(T &s, typename T::iterator it) { + return ++it; +} + +inline std::unordered_set _find_intersections( + const std::vector &segments) { + std::unordered_set intersections; + IntersectionEvents intersection_events; + std::vector events; + std::map> active; + events.reserve(2 * segments.size()); + for (std::size_t i = 0; i < segments.size(); i++) { + intersection_events[segments[i].top].tops.push_back(i); + intersection_events[segments[i].bottom].bottoms.push_back(i); + } + // std::cout << "Segments "; + // print_vector(std::cout, segments, "\n"); + while (!intersection_events.empty()) { + auto event_it = --intersection_events.end(); + // std::cout << "Event " << i << ": " << event_it->first << " tops: "; + // print_vector(std::cout, event_it->second.tops, ", bottoms: "); + // print_vector(std::cout, event_it->second.bottoms, "\n"); + // i++; + auto &event_data = event_it->second; + // std::cout << "Active: "; + // print_map(std::cout, active, "\n"); + if (!event_data.tops.empty()) { + // Current implementation is not optimal, but it is was + // faster to use this to have initial working version. + // TODO based on commented code fix edge cases. + for (const auto &active_el : active) { + for (auto event_index : event_data.tops) { + for (auto index : active_el.second) { + if (_do_intersect(segments[event_index], segments[index]) && + !_share_endpoint(segments[event_index], segments[index])) { + intersections.emplace(event_index, index); + } + } + } + } + // auto next = active.lower_bound(event_it->first); + // auto prev = pred(active, next); + // if (next != active.end()) { + // std::cout << "Next: " << next->first << " "; + // print_set(std::cout, next->second, "\n"); + // for (auto event_index : event_data.tops) { + // for (auto index : next->second) { + // if (_do_intersect(segments[event_index], segments[index]) && + // !_share_endpoint(segments[event_index], + // segments[index])) { + // intersections.emplace(event_index, index); + // } + // } + // } + // } + // if (prev != active.end()) { + // std::cout << "Prev: " << prev->first << " "; + // print_set(std::cout, prev->second, "\n"); + // for (auto event_index : event_data.tops) { + // for (auto index : prev->second) { + // if (_do_intersect(segments[event_index], segments[index]) && + // !_share_endpoint(segments[event_index], + // segments[index])) { intersections.emplace(event_index, + // index); + // + // } + // } + // } + // } + active[event_it->first].insert(event_data.tops.begin(), + event_data.tops.end()); + } + if (!event_data.bottoms.empty()) { + for (auto event_index : event_data.bottoms) { + // std::cout << "Event index: " << event_index << " " << + // active.size() + // << std::endl; + auto it = active.find(segments[event_index].top); + + // auto next = succ(active, it); + // auto prev = pred(active, it); + // std::cout << "It " << it->first << " segment " << + // segments[event_index] << std::endl; std::cout << "Is next " << + // (next != active.end()) << std::endl; std::cout << "Is prev " + // << (prev != active.end()) << std::endl; if (next != + // active.end() && prev != active.end()) { + // std::cout << "Next: " << next->first << std::endl; + // for (auto n_index : next->second) { + // for (auto p_index : prev->second) { + // if (_do_intersect(segments[n_index], segments[p_index]) + // && + // !_share_endpoint(segments[n_index], + // segments[p_index])) { intersections.emplace(n_index, + // p_index); + // } + // } + // } + // } + it->second.erase(event_index); + if (it->second.empty()) { + active.erase(it); + } + } + } + intersection_events.erase(event_it); + } + + return intersections; +} + +/** + * @brief Finds the intersection point of two line segments, if it exists. + * + * This function calculates the intersection point of two given line segments. + * Each segment is defined by two endpoints. If the segments do not intersect, + * the function returns the point (0, 0). + * + * @param s1 The first line segment. + * @param s2 The second line segment. + * @return The intersection point of the two segments, or (0, 0) if they do not + * intersect. + */ +inline std::vector _find_intersection(const point::Segment &s1, + const point::Segment &s2) { + // ReSharper disable CppJoinDeclarationAndAssignment + point::Point::coordinate_t a1, b1, a2, b2, det, x, y, t; + // point::Point::coordinate_t c1, c2, u; + a1 = s1.top.y - s1.bottom.y; + b1 = s1.bottom.x - s1.top.x; + // c1 = a1 * s1.bottom.x + b1 * s1.bottom.y; + a2 = s2.top.y - s2.bottom.y; + b2 = s2.bottom.x - s2.top.x; + // c2 = a2 * s2.bottom.x + b2 * s2.bottom.y; + det = a1 * b2 - a2 * b1; + if (det == 0) { + // collinear case + std::vector res; + if (s1.point_on_line(s2.bottom)) res.push_back(s2.bottom); + if (s1.point_on_line(s2.top)) res.push_back(s2.top); + if (s2.point_on_line(s1.bottom)) res.push_back(s1.bottom); + if (s2.point_on_line(s1.top)) res.push_back(s1.top); + return res; + } + t = ((s2.top.x - s1.top.x) * (s2.bottom.y - s2.top.y) - + (s2.top.y - s1.top.y) * (s2.bottom.x - s2.top.x)) / + det; + // u = ((s2.top.x - s1.top.x) * (s1.bottom.y - s1.top.y) - + // (s2.top.y - s1.top.y) * (s2.bottom.x - s2.top.x)) / det; + if (t < 0) return {s1.top}; + if (t > 1) return {s1.bottom}; + x = s1.top.x + t * b1; + y = s1.top.y + t * (-a1); + return {{x, y}}; +} +} // namespace partsegcore::intersection + +#endif // PARTSEGCORE_INTERSECTION_H diff --git a/src/PartSegCore_compiled_backend/triangulation/point.hpp b/src/PartSegCore_compiled_backend/triangulation/point.hpp new file mode 100644 index 0000000..61f8143 --- /dev/null +++ b/src/PartSegCore_compiled_backend/triangulation/point.hpp @@ -0,0 +1,232 @@ +// +// Created by Grzegorz Bokota on 11.10.24. +// + +#ifndef PARTSEGCORE_POINT_H +#define PARTSEGCORE_POINT_H + +#include +#include +#include + +namespace partsegcore::point { + +struct Vector; + +/* Point class with x and y coordinates */ +struct Point { + typedef float coordinate_t; + + coordinate_t x; + coordinate_t y; + bool operator==(const Point &p) const { return x == p.x && y == p.y; } + bool operator!=(const Point &p) const { return !(*this == p); } + Point(coordinate_t x, coordinate_t y) : x(x), y(y) {} + Point() = default; + + bool operator<(const Point &p) const { + if (this->y == p.y) { + return this->x < p.x; + } + return this->y < p.y; + } + + // add operator + Vector operator+(const Point &p) const; + + // subtract operator + Vector operator-(const Point &p) const; + + Point operator+(const Vector &v) const; + + // Point operator/(coordinate_t f) const { + // return {this->x / f, this->y / f}; + // } + // + // Point operator*(coordinate_t f) const { + // return {this->x * f, this->y * f}; + // } + + // Overload the << operator for Point + friend std::ostream &operator<<(std::ostream &os, const Point &point) { + os << "(x=" << point.x << ", y=" << point.y << ")"; + return os; + } + + struct PointHash { + std::size_t operator()(const Point &p) const { + return std::hash()(p.x) ^ + (std::hash()(p.y) << 1); + } + }; +}; + +struct Vector { + typedef Point::coordinate_t coordinate_t; + coordinate_t x; + coordinate_t y; + Vector(coordinate_t x, coordinate_t y) : x(x), y(y) {} + Vector() = default; + explicit Vector(const Point &p) : x(p.x), y(p.y) {} + + Vector operator+(const Vector &v) const { + return {this->x + v.x, this->y + v.y}; + } + Vector operator-(const Vector &v) const { + return {this->x - v.x, this->y - v.y}; + } + Vector operator/(coordinate_t f) const { return {this->x / f, this->y / f}; } + Vector operator*(coordinate_t f) const { return {this->x * f, this->y * f}; } + Vector operator-() const { return {-this->x, -this->y}; } +}; + +inline Vector Point::operator+(const Point &p) const { + return {this->x + p.x, this->y + p.y}; +} + +inline Vector Point::operator-(const Point &p) const { + return {this->x - p.x, this->y - p.y}; +} + +inline Point Point::operator+(const Vector &v) const { + return {this->x + v.x, this->y + v.y}; +} + +/*Struct to represent edge of polygon with points ordered*/ +struct Segment { + Point top{}; + Point bottom{}; + Segment(Point p1, Point p2) { + if (p1 == p2) { + std::ostringstream oss; + oss << "Segment cannot have two identical points: " << p1; + throw std::invalid_argument(oss.str()); + } + if (p1 < p2) { + bottom = p1; + top = p2; + } else { + bottom = p2; + top = p1; + } + } + + [[nodiscard]] bool is_horizontal() const { return bottom.y == top.y; } + [[nodiscard]] bool is_vertical() const { return bottom.x == top.x; } + + Segment() = default; + + bool operator<(const Segment &s) const { + if (this->bottom == s.bottom) { + return this->top < s.top; + } + return this->bottom < s.bottom; + } + + bool operator==(const Segment &s) const { + return this->bottom == s.bottom && this->top == s.top; + } + + bool operator!=(const Segment &s) const { return !(*this == s); } + + /** + * Computes the x-coordinate of a point on the line segment at a given + * y-coordinate. This function assumes that the line segment is defined + * between the `bottom` and `top` points. + * + * If the line segment is horizontal (`bottom.y == top.y`), it returns the + * x-coordinate of `bottom`. Otherwise, it computes the x-coordinate using + * linear interpolation. + * + * @param y The y-coordinate at which to find the corresponding x-coordinate + * on the line segment. + * @return The x-coordinate of the point on the line segment at the specified + * y-coordinate. + */ + [[nodiscard]] Point::coordinate_t point_on_line_x( + Point::coordinate_t y) const { + if (bottom.y == top.y) { + return bottom.x; + } + return bottom.x + + (y - bottom.y) * ((top.x - bottom.x) / (top.y - bottom.y)); + } + + /** + * t = \frac{(x_1 - x_2)(x_3 - x_2) + (y_1 - y_2)(y_3 - y_2)}{(x_3 - x_2)^2 + + *(y_3 - y_2)^2} + **/ + [[nodiscard]] double point_projection_factor(Point p) const { + return ((p.x - this->top.x) * (this->bottom.x - this->top.x) + + (p.y - this->top.y) * (this->bottom.y - this->top.y)) / + ((this->top.x - this->bottom.x) * (this->top.x - this->bottom.x) + + (this->top.y - this->bottom.y) * (this->top.y - this->bottom.y)); + } + + [[nodiscard]] bool point_on_line(Point p) const { + if (this->is_horizontal()) { + return (this->bottom.x <= p.x && p.x <= this->top.x); + } + if (this->is_vertical()) { + return (this->bottom.y <= p.y && p.y <= this->top.y); + } + auto x_cord = this->point_on_line_x(p.y); + return (this->bottom.x <= x_cord && x_cord <= this->top.x); + } + + /** + * Overloads the << operator for the Segment structure, enabling + * segments to be directly inserted into output streams. + * + * This operator outputs a Segment in the format: + * [bottom=, top=] + * + * @param os The output stream to which the Segment is inserted. + * @param segment The Segment instance to be inserted into the stream. + * @return The output stream after the Segment has been inserted. + */ + friend std::ostream &operator<<(std::ostream &os, const Segment &segment) { + os << "[bottom=" << segment.bottom << ", top=" << segment.top << "]"; + return os; + } + + /** + * A hash functor for the Segment structure. + * + * This functor computes a hash value for a given Segment instance. + * The hash is determined by combining the hash values of the bottom + * and top points of the segment. + * + * It is required for the Segment structure to be used as a key in + * unordered containers, such as unordered_map and unordered_set. + */ + struct SegmentHash { + std::size_t operator()(const Segment &segment) const { + std::size_t h1 = Point::PointHash()(segment.bottom); + std::size_t h2 = Point::PointHash()(segment.top); + return h1 ^ (h2 << 1); + } + }; +}; +} // namespace partsegcore::point + +// overload of hash function for +// an unordered map and set +namespace std { +template <> +struct hash { + size_t operator()(const partsegcore::point::Point &point) const noexcept { + return partsegcore::point::Point::PointHash()(point); + } +}; + +template <> +struct hash { + size_t operator()(const partsegcore::point::Segment &segment) const noexcept { + return partsegcore::point::Segment::SegmentHash()(segment); + } +}; + +} // namespace std + +#endif // PARTSEGCORE_POINT_H diff --git a/src/PartSegCore_compiled_backend/triangulation/triangulate.hpp b/src/PartSegCore_compiled_backend/triangulation/triangulate.hpp new file mode 100644 index 0000000..09bdb81 --- /dev/null +++ b/src/PartSegCore_compiled_backend/triangulation/triangulate.hpp @@ -0,0 +1,1293 @@ +#ifndef PARTSEGCORE_TRIANGULATE_H +#define PARTSEGCORE_TRIANGULATE_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "intersection.hpp" +#include "point.hpp" + +namespace partsegcore::triangulation { + +enum PointType { NORMAL, SPLIT, MERGE, INTERSECTION, EMPTY }; + +struct MonotonePolygon { + point::Point top{}; + point::Point bottom{}; + std::vector left; + std::vector right; + + MonotonePolygon() = default; + explicit MonotonePolygon(point::Point top) : top(top) {} + MonotonePolygon(point::Point top, point::Point bottom, + std::vector left, + std::vector right) + : top(top), + bottom(bottom), + left(std::move(left)), + right(std::move(right)) {} +}; + +struct Interval { + point::Point last_seen{}; + point::Segment left_segment; + point::Segment right_segment; + std::vector> polygons_list{}; + Interval() = default; + explicit Interval(const point::Point &p, const point::Segment &left, + const point::Segment &right) + : last_seen(p), left_segment(left), right_segment(right) {}; + explicit Interval(const point::Point &p, const point::Segment &left, + const point::Segment &right, + std::unique_ptr &polygon) + : last_seen(p), + left_segment(left), + right_segment(right), + polygons_list() { + polygons_list.emplace_back(std::move(polygon)); + }; + + void replace_segment(const point::Segment &old_segment, + const point::Segment &new_segment) { + if (left_segment == old_segment) { + left_segment = new_segment; + return; + } else if (right_segment == old_segment) { + right_segment = new_segment; + return; + } + throw std::runtime_error("Segment not found in interval"); + }; + [[nodiscard]] point::Segment opposite_segment( + const point::Segment &segment) const { + if (segment == left_segment) { + return right_segment; + } + if (segment == right_segment) { + return left_segment; + } + throw std::runtime_error("Segment not found in interval"); + }; + + // ostream operator for Interval + friend std::ostream &operator<<(std::ostream &os, const Interval &interval) { + os << "Last Seen: " << interval.last_seen + << ", Left Segment: " << interval.left_segment + << ", Right Segment: " << interval.right_segment + << ", Polygons count: " << interval.polygons_list.size(); + return os; + } +}; + +/** + * Comparator function to determine the relative positioning of two given + * segments. + * + * Compares two segments to identify if the first segment (`s1`) is to the left + * and below the second segment (`s2`) when moving left to right. The function + * assumes that the segments do not intersect. + * + * TODO Investigate if it could be implemented in a more efficient way. + * + * @param s1 The first segment to compare. + * @param s2 The second segment to compare. + * + * @return `true` if `s1` is to the left and below `s2`, `false` otherwise. + */ +inline bool left_to_right(const point::Segment &s1, const point::Segment &s2) { + if (s1.is_horizontal() && s2.is_horizontal()) { + return s1.bottom.x < s2.bottom.x; + } + if (s1.is_horizontal()) { + int i1 = intersection::_orientation(s2.bottom, s2.top, s1.top); + int i2 = intersection::_orientation(s2.bottom, s2.top, s1.bottom); + return (i1 == 2 || i2 == 2); + } + if (s2.is_horizontal()) { + int i1 = intersection::_orientation(s1.bottom, s1.top, s2.top); + int i2 = intersection::_orientation(s1.bottom, s1.top, s2.bottom); + return (i1 == 1 || i2 == 1); + } + + if ((s1.top.y < s2.top.y)) { + if (s1.top == s2.top) { + return intersection::_orientation(s2.top, s2.bottom, s1.bottom) == 1; + } + return intersection::_orientation(s2.top, s2.bottom, s1.top) == 1; + } + if (s1.top == s2.top) { + return intersection::_orientation(s1.top, s1.bottom, s2.bottom) == 2; + } + return intersection::_orientation(s1.top, s1.bottom, s2.top) == 2; +} + +inline bool left_right_share_top(const point::Segment &s1, + const point::Segment &s2) { + return intersection::_orientation(s1.bottom, s1.top, s2.bottom) == 1; +} + +inline bool left_right_share_bottom(const point::Segment &s1, + const point::Segment &s2) { + return intersection::_orientation(s1.top, s1.bottom, s2.top) == 2; +} + +struct SegmentLeftRightComparator { + bool operator()(const point::Segment &s1, const point::Segment &s2) const { + return left_to_right(s1, s2); + } +}; + +struct Triangle { + std::size_t x; + std::size_t y; + std::size_t z; + Triangle(std::size_t x, std::size_t y, std::size_t z) : x(x), y(y), z(z) {}; + Triangle() = default; +}; + +struct PointTriangle { + point::Point p1{}; + point::Point p2{}; + point::Point p3{}; + + PointTriangle(point::Point p1, point::Point p2, point::Point p3) { + /* Sort the points in clockwise order */ + if ((p2.x - p1.x) * (p3.y - p1.y) - (p2.y - p1.y) * (p3.x - p1.x) < 0) { + this->p1 = p3; + this->p2 = p2; + this->p3 = p1; + } else { + this->p1 = p1; + this->p2 = p2; + this->p3 = p3; + } + }; + PointTriangle() = default; +}; + +typedef std::size_t EdgeIndex; + +struct PointEdgeInfo { + EdgeIndex edge_index; + point::Point opposite_point; + PointEdgeInfo(EdgeIndex edge_index, point::Point opposite_point) + : edge_index(edge_index), opposite_point(opposite_point) {} + bool operator<(const PointEdgeInfo &e) const { + return opposite_point < e.opposite_point; + } +}; + +typedef std::unordered_map> + PointToEdges; + +/** + * Get a mapping from points to the list of edges that contain those points. + * Each list of edges is sorted by the point order within the edges. + * + * @param edges A vector of Segment objects representing the edges of a polygon. + * @return A PointToEdges map where each key is a point and the corresponding + * value is a vector of edges containing that point. + */ +inline PointToEdges get_points_edges(const std::vector &edges) { + PointToEdges point_to_edges; + for (std::size_t i = 0; i < edges.size(); i++) { + point_to_edges[edges[i].bottom].emplace_back(i, edges[i].top); + point_to_edges[edges[i].top].emplace_back(i, edges[i].bottom); + } + for (auto &point_to_edge : point_to_edges) { + std::sort(point_to_edge.second.begin(), point_to_edge.second.end(), + [](const PointEdgeInfo &a, const PointEdgeInfo &b) { + return b.opposite_point < a.opposite_point; + }); + } + return point_to_edges; +} + +/** + * Determines the type of given point based on its adjacent edges. + * + * If the point has more than two edges connected to it, it is considered an + * intersection point. If it has exactly two edges, further checks categorize + * it as one of split, merge, or normal points. + * + * - If both adjacent edges have their opposite ends before the given point `p`, + * the point is classified as a merge point. + * - If both adjacent edges have their opposite ends after the given point `p`, + * the point is classified as a split point. + * - Otherwise, it is categorized as a normal point. + * + * @param p The point to classify. + * @param point_to_edges A mapping from points to their adjacent edges. + * @return The type of the point as determined by its adjacent edges. + */ +inline PointType get_point_type(point::Point p, + const PointToEdges &point_to_edges) { + if (point_to_edges.count(p) == 0) return PointType::EMPTY; + if (point_to_edges.at(p).empty()) return PointType::EMPTY; + if (point_to_edges.at(p).size() != 2) return PointType::INTERSECTION; + + const auto &edges = point_to_edges.at(p); + if (edges[0].opposite_point < p && edges[1].opposite_point < p) + return PointType::SPLIT; + if (p < edges[0].opposite_point && p < edges[1].opposite_point) + return PointType::MERGE; + return PointType::NORMAL; +} + +struct MonotonePolygonBuilder { + std::map> segment_to_line{}; + std::vector edges{}; + PointToEdges point_to_edges{}; + std::vector monotone_polygons{}; + + MonotonePolygonBuilder() = default; + explicit MonotonePolygonBuilder(const std::vector &edges) + : edges(edges) { + this->point_to_edges = get_points_edges(edges); + } + + /** + * Retrieves the left and right segments (edges) of a monotone polygon + * for a given point when edges are sharing top point. + * + * @param p The point for which the left and right edges are to be retrieved. + * @return A pair containing the left and right segments respectively. + * + * The function looks up the segments associated with the given point `p` + * in the `point_to_edges` map. Based on their mutual orientation at the + * bottom point of one of the edges, it determines which segment is on the + * left and which is on the right. This orientation check is conducted + * using a helper function to determine the relative ordering of the edges. + * + * If the orientation is counter-clockwise (value 2), the first segment is + * on the right while the second segment is on the left. Otherwise, the + * segments maintain their original order. + */ + std::pair + get_left_right_edges_top(const point::Point &p) { + auto point_info = point_to_edges.at(p); + auto fst_idx = point_info[0].edge_index; + auto snd_idx = point_info[1].edge_index; + if (intersection::_orientation(edges[fst_idx].bottom, edges[fst_idx].top, + edges[snd_idx].bottom) == 2) { + return {edges[snd_idx], edges[fst_idx]}; + } + return {edges[fst_idx], edges[snd_idx]}; + } + + /** + * Retrieves the left and right edges of a monotone polygon + * for a given point at its bottom. + * + * @param p The point for which the left and right edges are to be retrieved. + * @return A pair containing the left and right segments (edges) in order. + * + * This function looks up the edges associated with the point `p` in the + * `point_to_edges` map. Based on the orientation of the segments + * at their top and bottom points, it determines which edge is on the left + * and which is on the right. The orientation is calculated using a helper + * function to determine the relative position of one segment relative to the + * other. + * + * If the orientation indicates a counter-clockwise arrangement, the function + * swaps the order of the edges to ensure the left and right ordering. + */ + std::pair + get_left_right_edges_bottom(const point::Point &p) { + auto point_info = point_to_edges.at(p); + auto fst_idx = point_info[0].edge_index; + auto snd_idx = point_info[1].edge_index; + if (intersection::_orientation(edges[fst_idx].top, edges[fst_idx].bottom, + edges[snd_idx].top) == 1) { + return {edges[snd_idx], edges[fst_idx]}; + } + return {edges[fst_idx], edges[snd_idx]}; + } + + /** + * Processes the end point of a segment within a monotone polygon. + * + * @param p The end point to be processed. + * + * The function retrieves two segments (edges) associated with the + * given point and checks the monotone polygon information of the + * top points of these segments. Depending on whether the end point + * belongs to one or two monotone polygons, it updates the monotone + * polygon information and stores completed polygons. + * + * It also handles the cleanup of relevant mappings and data structures + * when processing is complete. + */ + + void process_end_point(const point::Point &p, const point::Segment &edge_left, + const point::Segment &edge_right, + const std::shared_ptr &interval) { + for (auto &polygon : interval->polygons_list) { + polygon->bottom = p; + monotone_polygons.push_back(*polygon); + } + + segment_to_line.erase(edge_left); + segment_to_line.erase(edge_right); + } + + void process_merge_point(const point::Point &p) { + auto [edge_left, edge_right] = get_left_right_edges_bottom(p); + + auto left_interval = segment_to_line.at(edge_left); + auto right_interval = segment_to_line.at(edge_right); + + if (left_interval != right_interval) { +#ifdef DEBUG + if (right_interval->right_segment == edge_right) { + std::ostringstream oss; + oss << "The right edge of merge point should be the left edge of the " + "right interval.\nGot interval: " + << *right_interval << " and edge: " << edge_right; + throw std::runtime_error(oss.str()); + } +#endif + segment_to_line.erase(edge_right); + segment_to_line.erase(edge_left); + left_interval->right_segment = right_interval->right_segment; +#ifdef DEBUG + if (segment_to_line.count(right_interval->right_segment) == 0) { + throw std::runtime_error("Segment not found in the map2"); + } +#endif + segment_to_line[right_interval->right_segment] = left_interval; + left_interval->last_seen = p; + left_interval->polygons_list.back()->right.push_back(p); + right_interval->polygons_list.front()->left.push_back(p); + for (auto &polygon : right_interval->polygons_list) { + left_interval->polygons_list.push_back(std::move(polygon)); + } + } else { + // This is the end point + this->process_end_point(p, edge_left, edge_right, left_interval); + } + }; + + /** + * Processes a normal point within a monotone polygon. + * + * @param p The point to be processed. + * @param edge_top The top edge segment associated with the point. + * @param edge_bottom The bottom edge segment associated with the point. + * + * This function updates the monotone polygon interval associated with + * `edge_top` to include the new point `p`. If the interval contains + * multiple polygons, it closes all but one of them by setting their + * bottom to `p` and moving them to the list of completed monotone polygons. + * + * Depending on whether `edge_top` represents the right segment of the + * interval, the function then adds the point `p` to the appropriate (left + * or right) list of points in the remaining polygon and replaces `edge_top` + * with `edge_bottom` in the segment-to-line map. + * + * The function throws a `std::runtime_error` if `edge_top` is not found in + * the `segment_to_line` map. + */ + void _process_normal_point(const point::Point &p, + const point::Segment &edge_top, + const point::Segment &edge_bottom) { + if (segment_to_line.count(edge_top) == 0) { + throw std::runtime_error("Segment not found in the map"); + } + auto interval = segment_to_line.at(edge_top); + + if (interval->polygons_list.size() > 1) { + // end all the polygons, except 1 + if (edge_top == interval->right_segment) { + for (auto i = 1; i < interval->polygons_list.size(); i++) { + interval->polygons_list[i]->bottom = p; + monotone_polygons.push_back(*interval->polygons_list[i]); + } + + } else { + for (auto i = 0; i < interval->polygons_list.size() - 1; i++) { + interval->polygons_list[i]->bottom = p; + monotone_polygons.push_back(*interval->polygons_list[i]); + } + interval->polygons_list[0] = std::move(interval->polygons_list.back()); + } + interval->polygons_list.erase(interval->polygons_list.begin() + 1, + interval->polygons_list.end()); + } + + if (edge_top == interval->right_segment) { + interval->polygons_list[0]->right.push_back(p); + } else { + interval->polygons_list[0]->left.push_back(p); + } + + segment_to_line[edge_bottom] = interval; + interval->last_seen = p; + interval->replace_segment(edge_top, edge_bottom); + segment_to_line.erase(edge_top); +#ifdef DEBUG + if (segment_to_line.count(interval->left_segment) == 0) { + throw std::runtime_error("Left segment not found in the map"); + } + if (segment_to_line.count(interval->right_segment) == 0) { + throw std::runtime_error("Right segment not found in the map"); + } +#endif + }; + + /** + * Processes a normal point within a monotone polygon. + * + * @param p The point to be processed. + * + * This function identifies the two segments (top and bottom edges) + * associated with the point `p` using its index in the `point_to_edges` map. + * It then calls the helper method `_process_normal_point` + * with the identified segments to handle the actual processing logic. + */ + void process_normal_point(const point::Point &p) { + const point::Segment &edge_top = + edges[point_to_edges.at(p).at(0).edge_index]; + const point::Segment &edge_bottom = + edges[point_to_edges.at(p).at(1).edge_index]; + + _process_normal_point(p, edge_top, edge_bottom); + } + + void process_start_point(const point::Point &p, + const point::Segment &edge_left, + const point::Segment &edge_right) { + auto new_polygon = std::make_unique(p); + auto interval = + std::make_shared(p, edge_left, edge_right, new_polygon); + segment_to_line[edge_left] = interval; + segment_to_line[edge_right] = interval; + } + + void process_split_point(const point::Point &p) { + auto [edge_left, edge_right] = get_left_right_edges_top(p); + + // We need to find to which interval the point belongs. + // If the point does not belong to any interval, we treat + // it as a start point and create a new interval + + for (auto &segment_interval : segment_to_line) { + if (segment_interval.first == segment_interval.second->right_segment) { + // scan interval only once + continue; + } + // check if the current point is inside the quadrangle defined by edges of + // the interval + auto interval = segment_interval.second; + if (interval->left_segment.point_on_line_x(p.y) < p.x && + interval->right_segment.point_on_line_x(p.y) > p.x) { + // the point is inside the interval + + // update the sweep line + auto right_segment = interval->right_segment; + interval->right_segment = edge_left; + interval->last_seen = p; + segment_to_line[edge_left] = interval; + + auto new_interval = + std::make_shared(p, edge_right, right_segment); + segment_to_line[edge_right] = new_interval; + segment_to_line[right_segment] = new_interval; + + if (interval->polygons_list.size() == 1) { + std::unique_ptr new_polygon; + if (interval->polygons_list[0]->right.empty()) { + new_polygon = std::make_unique( + interval->polygons_list[0]->top); + } else { + new_polygon = std::make_unique( + interval->polygons_list[0]->right.back()); + } + new_polygon->left.push_back(p); + new_interval->polygons_list.emplace_back(std::move(new_polygon)); + interval->polygons_list[0]->right.push_back(p); + } + + if (interval->polygons_list.size() >= 2) { + interval->polygons_list[0]->right.push_back(p); + interval->polygons_list[interval->polygons_list.size() - 1] + ->left.push_back(p); + for (auto i = 1; i < interval->polygons_list.size() - 1; i++) { + interval->polygons_list[i]->bottom = p; + monotone_polygons.push_back(*interval->polygons_list[i]); + } + new_interval->polygons_list.push_back( + std::move(interval->polygons_list.back())); + interval->polygons_list.erase(interval->polygons_list.begin() + 1, + interval->polygons_list.end()); + } + return; + } + } + // the point is not inside any interval + this->process_start_point(p, edge_left, edge_right); + }; + + void process_intersection_point(const point::Point &p) { + auto point_info = point_to_edges.at(p); + std::set processed_segments; + std::vector segments_to_normal_process; + std::vector + top_segments; // segments above the intersection point + std::vector + bottom_segments; // segments below the intersection point + + for (auto &edge_info : point_info) { + auto &edge = edges[edge_info.edge_index]; + if (processed_segments.count(edge) > 0) { + continue; + } + if (segment_to_line.count(edge) > 0) { + auto interval = segment_to_line.at(edge); + auto opposite_edge = interval->opposite_segment(edge); + if (edge.bottom == p && opposite_edge.bottom == p) { + process_end_point(p, edge, opposite_edge, interval); + processed_segments.insert(edge); + processed_segments.insert(opposite_edge); + continue; + } + } + if (edge.top == p) { + bottom_segments.push_back(edge); + } else { + top_segments.push_back(edge); + } + segments_to_normal_process.push_back(edge); + } + + // after this loop we should have at most 4 segments + + if (bottom_segments.empty() && top_segments.empty()) { + // all segments were processed + return; + } + + std::sort(top_segments.begin(), top_segments.end(), + left_right_share_bottom); + std::sort(bottom_segments.begin(), bottom_segments.end(), + left_right_share_top); + + auto bottom_begin = bottom_segments.begin(); + auto bottom_end = bottom_segments.end(); + auto top_begin = top_segments.begin(); + if (!top_segments.empty()) { + if (top_segments.front() == + segment_to_line.at(top_segments.front())->right_segment) { + ++bottom_begin; + ++top_begin; + _process_normal_point(p, top_segments.front(), bottom_segments.front()); + } + if (top_begin != top_segments.end() && + top_segments.back() == + segment_to_line.at(top_segments.back())->left_segment) { + --bottom_end; + _process_normal_point(p, top_segments.back(), bottom_segments.back()); + } + } + + for (auto it = bottom_begin; it < bottom_end; it += 2) { + process_start_point(p, *it, *(it + 1)); + } + }; +}; + +/** + * Checks if a given polygon is convex. + * + * This function takes a polygon represented as a vector of points and + * determines if it is convex. A convex polygon is one where all the interior + * angles are less than 180 degrees, meaning that the edges never turn back on + * themselves. + * + * @param polygon A vector of points representing the vertices of the polygon in + * order. + * @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) + 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) + return false; + triangle_orientation = intersection::_orientation(polygon[polygon.size() - 1], + polygon[0], polygon[1]); + if (triangle_orientation != 0 && triangle_orientation != orientation) + return false; + return true; +} + +/** + * Divides a convex polygon into triangles using the fan triangulation method. + * Starting from the first point in the polygon, each triangle is formed by + * connecting the first point with two consecutive points from the polygon list. + * + * @param polygon A vector of Point objects representing the vertices of the + * polygon. + * @return A vector of Triangle objects where each triangle is defined by + * indices corresponding to the vertices in the input polygon. + */ +inline std::vector _triangle_convex_polygon( + const std::vector &polygon) { + std::vector result; + for (std::size_t i = 1; i < polygon.size() - 1; i++) { + if (intersection::_orientation(polygon[0], polygon[i], polygon[i + 1]) != + 0) { + result.emplace_back(0, i, i + 1); + } + } + return result; +} + +/** + * Construct triangles using the current point and edges opposite to it. + * + * This function takes a current point and a stack of points, then builds + * triangles by connecting the current point to adjoining points in the stack. + * Each triangle is added to the result vector as a `PointTriangle` object. + * + * Steps performed: + * 1. Iterate over the stack to form triangles with the current point. + * 2. Replace the first element of the stack with the last element. + * 3. Set the second element of the stack to the current point. + * 4. Remove all elements from the stack except the first two. + * + * @param stack A vector of `point::Point` objects representing the stack of + * points. + * @param result A vector to store the resultant `PointTriangle` objects. + * @param current_point The current `point::Point` used to form triangles with + * points in the stack. + */ +inline void _build_triangles_opposite_edge(std::vector &stack, + std::vector &result, + point::Point current_point) { + for (std::size_t i = 0; i < stack.size() - 1; i++) { + result.emplace_back(current_point, stack[i], stack[i + 1]); + } + stack[0] = stack.back(); + stack[1] = current_point; + // remove all elements except two firsts + stack.erase(stack.begin() + 2, stack.end()); +} + +/** + * Constructs triangles from the current edge of the y-monotone and updates the + * stack and result vectors. + * + * This function processes a stack of points and an incoming point to form + * triangles that are part of the processed area. The triangles are determined + * based on the expected orientation and are added to the result vector. The + * stack is then updated to reflect the current status of the processed area. + * + * @param stack A reference to a std::vector of point::Point representing + * the current state of the stack. + * @param result A reference to a std::vector of PointTriangle where the + * generated triangles will be stored. + * @param current_point The incoming point::Point to be considered for new + * triangles and updating the stack. + * @param expected_orientation The expected orientation (clockwise or + * counterclockwise) that will guide the triangle formation. + */ +inline void _build_triangles_current_edge(std::vector &stack, + std::vector &result, + point::Point current_point, + int expected_orientation) { + auto it1 = stack.rbegin(); + auto it2 = stack.rbegin() + 1; + while (it2 != stack.rend() && + intersection::_orientation(*it2, *it1, current_point) == + expected_orientation) { + result.emplace_back(current_point, *it1, *it2); + ++it1; + ++it2; + } + stack.erase(it1.base(), stack.end()); + stack.push_back(current_point); +} + +/** + * @enum Side + * @brief An enumeration to represent different sides. + * + * This enumeration defines constants for various sides, + * in particular: + * - TOP_OR_BOTTOM: Represents the top side. + * - LEFT: Represents the left side. + * - RIGHT: Represents the right side. + */ +enum Side { + TOP_OR_BOTTOM = 0, + LEFT = 2, + RIGHT = 1, +}; + +/** + * Triangulates a given monotone polygon. + * + * The function takes a monotone polygon as input and returns a vector + * of triangles that represent the triangulation of the polygon. The + * polygon must be y-monotone, meaning that every line segment + * parallel to the y-axis intersects the polygon at most twice. + * + * The algorithm works by sorting points of the polygon, merging the + * left and right chains, and then using a stack-based approach to + * recursively build triangles. + * + * @param polygon The monotone polygon to be triangulated. + * @return A vector of PointTriangle representing the triangles of the + * triangulated polygon. + */ +inline std::vector triangulate_monotone_polygon( + const MonotonePolygon &polygon) { + std::vector result; + std::size_t left_index = 0; + std::size_t right_index = 0; + std::vector stack; + std::vector> points; + + result.reserve(polygon.left.size() + polygon.right.size()); + + points.reserve(polygon.left.size() + polygon.right.size() + 2); + points.emplace_back(polygon.top, Side::TOP_OR_BOTTOM); + while (left_index < polygon.left.size() && + right_index < polygon.right.size()) { + if (polygon.left[left_index] < polygon.right[right_index]) { + points.emplace_back(polygon.right[right_index], Side::RIGHT); + right_index++; + } else { + points.emplace_back(polygon.left[left_index], Side::LEFT); + left_index++; + } + } + while (left_index < polygon.left.size()) { + points.emplace_back(polygon.left[left_index], Side::LEFT); + left_index++; + } + while (right_index < polygon.right.size()) { + points.emplace_back(polygon.right[right_index], Side::RIGHT); + right_index++; + } + points.emplace_back(polygon.bottom, Side::TOP_OR_BOTTOM); + + stack.push_back(points[0].first); + stack.push_back(points[1].first); + Side side = points[1].second; + + for (std::size_t i = 2; i < points.size(); i++) { + if (side == points[i].second) { + _build_triangles_current_edge(stack, result, points[i].first, side); + } else { + _build_triangles_opposite_edge(stack, result, points[i].first); + } + side = points[i].second; + } + return result; +} + +/** + * Calculates the edges of a polygon represented as a sequence of points. + * + * This function takes a vector of points representing the vertices of a polygon + * and returns a vector of segments representing the edges of the polygon. + * The edges are created by connecting each point to the next, with the final + * point connecting back to the first point to close the polygon. + * + * @param polygon A vector of points representing the vertices of the polygon. + * @return A vector of segments representing the edges of the polygon. + */ +inline std::vector calc_edges( + const std::vector &polygon) { + std::vector edges; + edges.reserve(polygon.size()); + for (std::size_t i = 0; i < polygon.size() - 1; i++) { + edges.emplace_back(polygon[i], polygon[i + 1]); + } + edges.emplace_back(polygon[polygon.size() - 1], polygon[0]); + return edges; +} + +/** + * Calculates the edges of polygons from a list of polygons, provided as + * a list of points for each polygon. + * + * Each polygon is represented by a vector of points, where each point is + * defined as an instance of `point::Point`. The function iterates through + * each point in each polygon and creates edges (segments) between + * consecutive points as well as between the last point and the first point + * of each polygon. + * + * @param polygon_list A vector of polygons, with each polygon represented + * as a vector of `point::Point` instances. + * + * @return A vector of `point::Segment` instances representing the edges of + * all polygons in the input list. + */ +inline std::vector calc_edges( + const std::vector> &polygon_list) { + std::vector edges; + std::size_t points_count = 0; + for (const auto &polygon : polygon_list) { + points_count += polygon.size(); + } + edges.reserve(points_count); + for (const auto &polygon : polygon_list) { + for (std::size_t i = 0; i < polygon.size() - 1; i++) { + edges.emplace_back(polygon[i], polygon[i + 1]); + } + if (polygon.back() != polygon.front()) + edges.emplace_back(polygon.back(), polygon.front()); + } + return edges; +} + +/** + * Calculates and returns a list of unique (deduplicated) edges from a list of + * polygons. + * + * This function receives a list of polygons, where each polygon is defined by a + * series of points. It identifies the edges of the polygons and returns those + * edges that are unique (i.e., edges that appear an odd number of times across + * all polygons). Any edge that appears an even number of times is considered to + * be a duplicate and is removed. + * + * @param polygon_list A reference to a vector containing vectors of points, + * where each inner vector represents a polygon. + * @return A vector of unique edges (of type point::Segment) that are not + * duplicated across the polygons. + */ +inline std::vector calc_dedup_edges( + const std::vector> &polygon_list) { + std::set edges_set; + point::Segment edge; + + for (const auto &polygon : polygon_list) { + for (std::size_t i = 0; i < polygon.size() - 1; i++) { + edge = point::Segment(polygon[i], polygon[i + 1]); + if (edges_set.count(edge) == 0) { + edges_set.insert(edge); + } else { + edges_set.erase(edge); + } + } + if (polygon.back() == polygon.front()) continue; + edge = point::Segment(polygon.back(), polygon.front()); + if (edges_set.count(edge) == 0) { + edges_set.insert(edge); + } else { + edges_set.erase(edge); + } + } + return {edges_set.begin(), edges_set.end()}; +} + +/** + * @brief Finds intersection points in a polygon and adds mid-points for all + * intersections. + * + * The function takes a vector of points defining a polygon and finds all edge + * intersections. It then adds mid-points for all such intersections. + * + * @param polygon_list The list of polygons defined by a vector of ov vector of + * Point objects. + * @return A new vector of Point objects representing the polygon with added + * intersection points. + */ +inline std::vector> find_intersection_points( + const std::vector> &polygon_list) { + /* find all edge intersections and add mid-points for all such intersection + * places*/ + auto edges = calc_edges(polygon_list); + + auto intersections = intersection::_find_intersections(edges); + if (intersections.empty()) return polygon_list; + std::unordered_map>> + intersections_points; + for (const auto &intersection : intersections) { + auto inter_points = intersection::_find_intersection( + edges[intersection.first], edges[intersection.second]); + for (auto inter_point : inter_points) { + intersections_points[intersection.first].emplace_back( + edges[intersection.first].point_projection_factor(inter_point), + inter_point); + intersections_points[intersection.second].emplace_back( + edges[intersection.second].point_projection_factor(inter_point), + inter_point); + } + } + for (auto &intersections_point : intersections_points) { + // points_count += intersections_point.second.size() - 1; + auto edge = edges[intersections_point.first]; + intersections_point.second.emplace_back(-1, edge.top); + intersections_point.second.emplace_back(2, edge.bottom); + std::sort(intersections_point.second.begin(), + intersections_point.second.end()); + } + + std::vector> new_polygons_list; + + for (const auto &polygon : polygon_list) { + std::vector new_polygon; + new_polygon.reserve(polygon.size() * 2); + new_polygon.push_back(polygon[0]); + for (std::size_t i = 0; i < polygon.size(); i++) { + auto point = polygon[i]; + if (new_polygon.back() != point) new_polygon.push_back(point); + if (intersections_points.count(i)) { + auto &new_points = intersections_points[i]; + if (new_points[0].second == point) { + for (auto it = new_points.begin() + 1; it != new_points.end(); ++it) { + if (new_polygon.back() != it->second) + new_polygon.push_back(it->second); + } + } else { + for (auto it = new_points.rbegin() + 1; it != new_points.rend(); + ++it) { + if (new_polygon.back() != it->second) + new_polygon.push_back(it->second); + } + } + } + } + if (new_polygon.size() > 1 && new_polygon.front() == new_polygon.back()) + new_polygon.pop_back(); + new_polygons_list.push_back(new_polygon); + } + return new_polygons_list; +} + +inline std::vector find_intersection_points( + const std::vector &polygon) { + auto new_polygon = find_intersection_points( + std::vector>({polygon})); + return new_polygon[0]; +} + +inline std::vector _sorted_polygons_points( + const std::vector> &polygon_list) { + std::vector result; + std::unordered_set visited; + for (const auto &polygon : polygon_list) { + for (const auto &point : polygon) { + if (visited.count(point) == 0) { + result.push_back(point); + visited.insert(point); + } + } + } + std::sort( + result.begin(), result.end(), + [](const point::Point &p1, const point::Point &p2) { return p2 < p1; }); + return result; +} + +/* + This is an implementation of the polygon with sweeping lines. + It assumes that there are no edge intersections but may be a point with + more than 2 edges. + Described on this lecture: + https://www.youtube.com/playlist?list=PLtTatrCwXHzEqzJMaTUFgqoCNllgwk4DH + and in the book: + de Berg,M. et al. (2008) Computational Geometry Springer Berlin Heidelberg. + */ +inline std::pair, std::vector> +sweeping_line_triangulation( + const std::vector> &polygon_list) { + std::vector result; + auto edges = calc_dedup_edges(polygon_list); + MonotonePolygonBuilder builder(edges); + + PointToEdges point_to_edges = get_points_edges(edges); + + std::vector sorted_points = + _sorted_polygons_points(polygon_list); + + result.reserve(sorted_points.size() - 2); + int idx = 0; + for (auto &sorted_point : sorted_points) { + auto point_type = get_point_type(sorted_point, point_to_edges); + switch (point_type) { + case PointType::NORMAL: + // Change edge adjusted to the current sweeping line. + // Adds edge to monotone polygon. + builder.process_normal_point(sorted_point); + break; + case PointType::SPLIT: + // Split a sweeping line on two lines, + // adds edge sor cutting polygon on two parts. + builder.process_split_point(sorted_point); + break; + case PointType::MERGE: + // merge two sweeping lines into one + // merge two intervals into one + // It may be the end of the polygon, then finish the + // monotone polygon + builder.process_merge_point(sorted_point); + break; + case PointType::INTERSECTION: + // this is a merge and split point at the same time + // this is not described in the original algorithm, + // but we need it to handle self-intersecting polygons + builder.process_intersection_point(sorted_point); + break; + case PointType::EMPTY: + // this is a point without edges (removed by deduplication) + break; + } + ++idx; + } + std::unordered_map point_to_index; + point_to_index.reserve(sorted_points.size()); + for (std::size_t i = 0; i < sorted_points.size(); i++) { + point_to_index[sorted_points[i]] = i; + } + for (auto &monotone_polygon : builder.monotone_polygons) { + auto triangles = triangulate_monotone_polygon(monotone_polygon); + for (auto &triangle : triangles) { + result.emplace_back(point_to_index[triangle.p1], + point_to_index[triangle.p2], + point_to_index[triangle.p3]); + } + } + return std::make_pair(result, sorted_points); +} + +// calculate the triangulation out of a symmetric difference out of a list of +// polygons + +inline std::pair, std::vector> +triangulate_polygon_face( + const std::vector> &polygon_list) { + if (polygon_list.empty()) + // empty list + return std::make_pair(std::vector(), std::vector()); + if (polygon_list.size() == 1) { + // only one polygon in the list + std::vector polygon = polygon_list[0]; + if (polygon.size() < 3) + return std::make_pair(std::vector(), polygon); + if (polygon.size() == 3) + return std::make_pair(std::vector({Triangle(0, 1, 2)}), + polygon); + if (polygon.size() == 4) { + if (partsegcore::intersection::_orientation(polygon[0], polygon[1], + polygon[2]) != + partsegcore::intersection::_orientation(polygon[0], polygon[3], + polygon[2])) + return std::make_pair( + std::vector({Triangle(0, 1, 2), Triangle(0, 3, 2)}), + polygon); + } + + if (_is_convex(polygon)) + return std::make_pair(_triangle_convex_polygon(polygon), polygon); + } + + // Implement the sweeping line algorithm for triangulation + // described in this lecture: + // https://www.youtube.com/playlist?list=PLtTatrCwXHzEqzJMaTUFgqoCNllgwk4DH + // + return sweeping_line_triangulation(find_intersection_points(polygon_list)); +} + +inline std::pair, std::vector> +triangulate_polygon_face(const std::vector &polygon) { + // #if DDEBUG + // try{ + // #endif + return triangulate_polygon_face( + std::vector>({polygon})); + // #if DDEBUG + // } catch (const std::exception &e) { + // std::cerr << "Polygon: ["; + // for (const auto &point : polygon) { + // std::cerr << "(" << point.x << ", " << point.y << "), "; + // } + // std::cerr << "]" << std::endl; + // std::cerr << "Error: " << e.what() << std::endl; + // throw e; + // } + // #endif +} + +inline point::Point::coordinate_t vector_length(point::Point p1, + point::Point p2) { + return std::sqrt((p1.x - p2.x) * (p1.x - p2.x) + + (p1.y - p2.y) * (p1.y - p2.y)); +} + +struct PathTriangulation { + std::vector triangles; + std::vector centers; + std::vector offsets; + + void reserve(std::size_t size) { + triangles.reserve(size); + centers.reserve(size); + offsets.reserve(size); + } + + void fix_triangle_orientation() { + point::Point p1{}, p2{}, p3{}; + for (auto &triangle : triangles) { + p1 = centers[triangle.x] + offsets[triangle.x]; + p2 = centers[triangle.y] + offsets[triangle.y]; + p3 = centers[triangle.z] + offsets[triangle.z]; + if ((p2.x - p1.x) * (p3.y - p1.y) - (p2.y - p1.y) * (p3.x - p1.x) < 0) { + std::swap(triangle.x, triangle.z); + } + // if (intersection::_orientation(p1, p2, p3) == 1) { + // std::swap(triangle.x, triangle.z); + // } + } + } +}; + +inline std::pair +sign_abs(point::Point::coordinate_t x) { + if (x < 0) { + return {-1, -x}; + } + if (x > 0) { + return {1, x}; + } + return {0, 0}; +} + +inline point::Point::coordinate_t add_triangles_for_join( + PathTriangulation &triangles, point::Point p1, point::Point p2, + point::Point p3, point::Point::coordinate_t prev_length, double cos_limit, + bool bevel) { + std::size_t idx = triangles.offsets.size(); + point::Vector mitter{}; + point::Point::coordinate_t length = vector_length(p2, p3); + point::Vector p1_p2_diff_norm = (p2 - p1) / prev_length; + point::Vector p2_p3_diff_norm = (p3 - p2) / length; + + point::Point::coordinate_t cos_angle = p1_p2_diff_norm.x * p2_p3_diff_norm.x + + p1_p2_diff_norm.y * p2_p3_diff_norm.y; + point::Point::coordinate_t sin_angle = p1_p2_diff_norm.x * p2_p3_diff_norm.y - + p1_p2_diff_norm.y * p2_p3_diff_norm.x; + + triangles.centers.push_back(p2); + triangles.centers.push_back(p2); + + if (sin_angle == 0) { + mitter = {p1_p2_diff_norm.y / 2, -p1_p2_diff_norm.x / 2}; + } else { + point::Point::coordinate_t scale_factor = 1 / sin_angle; + if (bevel || cos_angle < cos_limit) { + /* Bevel join + * There is a need to check if inner vector is not to long + * See https://github.com/napari/napari/pull/7268#user-content-bevel-cut + */ + auto [sign, mag] = sign_abs(scale_factor); + scale_factor = + sign * (float)0.5 * std::min(mag, std::min(prev_length, length)); + } + mitter = (p1_p2_diff_norm - p2_p3_diff_norm) * scale_factor * 0.5; + } + if (bevel || cos_angle < cos_limit) { + triangles.centers.push_back(p2); + triangles.triangles.emplace_back(idx, idx + 1, idx + 2); + if (sin_angle < 0) { + triangles.offsets.push_back(mitter); + triangles.offsets.emplace_back(-p1_p2_diff_norm.y * 0.5, + p1_p2_diff_norm.x * 0.5); + triangles.offsets.emplace_back(-p2_p3_diff_norm.y * 0.5, + p2_p3_diff_norm.x * 0.5); + triangles.triangles.emplace_back(idx, idx + 2, idx + 3); + triangles.triangles.emplace_back(idx + 2, idx + 3, idx + 4); + } else { + triangles.offsets.emplace_back(p1_p2_diff_norm.y * 0.5, + -p1_p2_diff_norm.x * 0.5); + triangles.offsets.push_back(-mitter); + triangles.offsets.emplace_back(p2_p3_diff_norm.y * 0.5, + -p2_p3_diff_norm.x * 0.5); + triangles.triangles.emplace_back(idx + 1, idx + 2, idx + 3); + triangles.triangles.emplace_back(idx + 1, idx + 3, idx + 4); + } + } else { + triangles.offsets.push_back(mitter); + triangles.offsets.push_back(-mitter); + triangles.triangles.emplace_back(idx, idx + 1, idx + 2); + triangles.triangles.emplace_back(idx + 1, idx + 2, idx + 3); + } + + return length; +} + +inline PathTriangulation triangulate_path_edge( + const std::vector &path, bool closed, float limit, + bool bevel) { + if (path.size() < 2) + return {{{0, 1, 3}, {1, 3, 2}}, + {path[0], path[0], path[0], path[0]}, + {{0, 0}, {0, 0}, {0, 0}, {0, 0}}}; + PathTriangulation result; + point::Vector norm_diff{}; + result.reserve(path.size() * 3); + double cos_limit = 1.0 / (limit * limit / 2) - 1.0; + point::Point::coordinate_t prev_length{}; + + if (closed) { + prev_length = vector_length(path[0], path[path.size() - 1]); + prev_length = + add_triangles_for_join(result, path[path.size() - 1], path[0], path[1], + prev_length, cos_limit, bevel); + } else { + prev_length = vector_length(path[0], path[1]); + norm_diff = (path[1] - path[0]) / prev_length; + result.centers.push_back(path[0]); + result.centers.push_back(path[0]); + result.offsets.emplace_back(norm_diff.y * 0.5, -norm_diff.x * 0.5); + result.offsets.push_back(-result.offsets.back()); + result.triangles.emplace_back(0, 1, 2); + result.triangles.emplace_back(1, 2, 3); + } + for (std::size_t i = 1; i < path.size() - 1; i++) { + prev_length = + add_triangles_for_join(result, path[i - 1], path[i], path[i + 1], + prev_length, cos_limit, bevel); + } + if (closed) { + add_triangles_for_join(result, path[path.size() - 2], path[path.size() - 1], + path[0], prev_length, cos_limit, bevel); + result.centers.push_back(result.centers[0]); + result.centers.push_back(result.centers[0]); + result.offsets.push_back(result.offsets[0]); + result.offsets.push_back(result.offsets[1]); + } else { + norm_diff = (path[path.size() - 1] - path[path.size() - 2]) / prev_length; + result.centers.push_back(path[path.size() - 1]); + result.centers.push_back(path[path.size() - 1]); + result.offsets.emplace_back(norm_diff.y * 0.5, -norm_diff.x * 0.5); + result.offsets.push_back(-result.offsets.back()); + } + result.fix_triangle_orientation(); + return result; +} + +} // namespace partsegcore::triangulation + +#endif // PARTSEGCORE_TRIANGULATE_H diff --git a/src/tests/test_triangulate.py b/src/tests/test_triangulate.py new file mode 100644 index 0000000..bc092bc --- /dev/null +++ b/src/tests/test_triangulate.py @@ -0,0 +1,658 @@ +import numpy as np +import pytest + +from PartSegCore_compiled_backend.triangulate import ( + do_intersect, + find_intersection_point, + find_intersection_points_py, + find_intersections, + is_convex, + on_segment, + orientation, + segment_left_to_right_comparator, + triangle_convex_polygon, + triangulate_monotone_polygon_py, + triangulate_path_edge_py, + triangulate_polygon_numpy, + triangulate_polygon_numpy_li, + triangulate_polygon_py, + triangulate_polygon_with_edge_numpy_li, +) + + +def test_on_segment(): + assert on_segment((0, 0), (0, 1), (0, 2)) + assert not on_segment((0, 0), (1, 1), (0, 3)) + + +def test_orientation_basic(): + assert orientation((0, 0), (0, 1), (0, 2)) == 0 + assert orientation((0, 0), (0, 2), (0, 1)) == 0 + assert orientation((0, 2), (0, 0), (0, 1)) == 0 + assert orientation((0, 0), (0, 1), (1, 2)) == 1 + assert orientation((0, 0), (0, 1), (-1, 2)) == 2 + + +@pytest.mark.parametrize( + ('p1', 'p2', 'p3'), + [ + ((0, 0), (1, 1), (2, 0)), + ((0, 0), (1, 100), (2, 0)), + ((0, 0), (100, 1), (1, 0)), + ((0, 1), (1, 1), (2, 0)), + ((0, 1), (1, 1), (2000, 0)), + ], +) +def test_orientation_split(p1, p2, p3): + assert orientation(p1, p2, p3) == 1 + assert orientation(p3, p2, p1) == 2 + + +@pytest.mark.parametrize( + ('p1', 'p2', 'p3'), + [ + ((0, 1), (1, 0), (2, 1)), + ((0, 1), (1, 0), (2, 0)), + ((0, 1), (1, 0), (200, 0)), + ], +) +def test_orientation_merge(p1, p2, p3): + assert orientation(p1, p2, p3) == 2 + assert orientation(p3, p2, p1) == 1 + + +def test_do_intersect(): + assert do_intersect(((0, -1), (0, 1)), ((1, 0), (-1, 0))) + assert not do_intersect(((0, -1), (0, 1)), ((1, 0), (2, 0))) + + +def test_find_intersections(): + r""" + First test case: + (1, 0) --- (1, 1) + | | + (0, 0) --- (0, 1) + + Second test case: + (1, 0) --- (1, 1) + \ / + \ / + \ / + X + / \ + / \ + / \ + (0, 0) --- (0, 1) + """ + assert find_intersections([[(0, 0), (0, 1)], [(0, 1), (1, 1)], [(1, 1), (1, 0)], [(1, 0), (0, 0)]]) == [] + assert find_intersections([[(0, 0), (0, 1)], [(0, 1), (1, 0)], [(1, 0), (1, 1)], [(1, 1), (0, 0), (0, 1)]]) == [ + (1, 3) + ] + + +@pytest.mark.parametrize( + ('segments', 'expected'), + [ + # No intersections, simple square + ([[(0, 0), (0, 1)], [(0, 1), (1, 1)], [(1, 1), (1, 0)], [(1, 0), (0, 0)]], []), + # One intersection, crossing diagonals + ([[(0, 0), (2, 2)], [(2, 0), (0, 2)]], [(0, 1)]), + # Multiple intersections, complex shape + ( + [[(0, 0), (2, 2)], [(2, 0), (0, 2)], [(1, 0), (1, 2)], [(0, 1), (2, 1)]], + {(0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3)}, + ), + # No intersections, non-intersecting lines + ([[(0, 0), (1, 1)], [(2, 2), (3, 3)]], []), + # One intersection, T-shaped intersection + ([[(0, 0), (2, 0)], [(1, -1), (1, 1)]], [(0, 1)]), + # Multiple intersections, grid shape + ( + [ + [(0, 0), (2, 0)], + [(0, 1), (2, 1)], + [(0, 2), (2, 2)], + [(0, 0), (0, 2)], + [(1, 0), (1, 2)], + [(2, 0), (2, 2)], + ], + {(0, 4), (1, 3), (1, 4), (1, 5), (2, 4)}, + ), + ], + ids=[ + 'No intersections, simple square', + 'One intersection, crossing diagonals', + 'Multiple intersections, complex shape', + 'No intersections, non-intersecting lines', + 'One intersection, T-shaped intersection', + 'Multiple intersections, grid shape', + ], +) +def test_find_intersections_param(segments, expected): + assert set(find_intersections(segments)) == set(expected) + + +@pytest.mark.parametrize( + ('segment1', 'segment2', 'expected'), + [ + (((0, 0), (2, 2)), ((0, 2), (2, 0)), (1, 1)), # Intersecting diagonals + (((0, 0), (1, 1)), ((1, 0), (0, 1)), (0.5, 0.5)), # Intersecting diagonals + (((0, 0), (2, 2)), ((2, 2), (0, 2)), (2, 2)), # Touching at one point + (((0, 0), (2, 2)), ((1, 0), (1, 2)), (1, 1)), # Vertical line intersection + (((0, 1), (2, 1)), ((2, 0), (2, 2)), (2, 1)), # Touching in the middle + ], + ids=[ + 'Intersecting diagonals (1, 1)', + 'Intersecting diagonals (0.5, 0.5)', + 'Touching at one point (2, 2)', + 'Vertical line intersection (1, 1)', + 'Touching in the middle (2, 1)', + ], +) +def test_find_intersection_point(segment1, segment2, expected): + assert len(find_intersection_point(segment1, segment2)) == 1 + assert find_intersection_point(segment1, segment2)[0] == expected + + +@pytest.mark.parametrize( + ('segment1', 'segment2', 'expected'), + [ + (((0, 0), (2, 0)), ((1, 0), (3, 0)), [(1, 0), (2, 0)]), # Horizontal segments + (((1, 0), (3, 0)), ((0, 0), (2, 0)), [(1, 0), (2, 0)]), # Horizontal segments + (((0, 0), (0, 2)), ((0, 1), (0, 3)), [(0, 1), (0, 2)]), # Vertical segments + (((0, 1), (0, 3)), ((0, 0), (0, 2)), [(0, 1), (0, 2)]), # Vertical segments + (((0, 0), (2, 2)), ((1, 1), (3, 3)), [(1, 1), (2, 2)]), # Diagonal segments + (((1, 1), (3, 3)), ((0, 0), (2, 2)), [(1, 1), (2, 2)]), # Diagonal segments + ], +) +def test_find_intersection_points_colinear_overlap(segment1, segment2, expected): + assert sorted(find_intersection_point(segment1, segment2)) == expected + + +@pytest.mark.parametrize( + ('segment1', 'segment2', 'expected'), + [ + (((0, 0), (3, 0)), ((1, 0), (2, 0)), [(1, 0), (2, 0)]), # Horizontal segments + (((1, 0), (2, 0)), ((0, 0), (3, 0)), [(1, 0), (2, 0)]), # Horizontal segments + (((0, 0), (0, 3)), ((0, 1), (0, 2)), [(0, 1), (0, 2)]), # Vertical segments + (((0, 1), (0, 2)), ((0, 0), (0, 3)), [(0, 1), (0, 2)]), # Vertical segments + (((0, 0), (3, 3)), ((1, 1), (2, 2)), [(1, 1), (2, 2)]), # Diagonal segments + (((1, 1), (2, 2)), ((0, 0), (3, 3)), [(1, 1), (2, 2)]), # Diagonal segments + ], +) +def test_find_intersection_point_colinear_inside(segment1, segment2, expected): + assert sorted(find_intersection_point(segment1, segment2)) == expected + + +def test_is_convex(): + assert is_convex([(0, 0), (0, 1), (1, 1), (1, 0)]) + assert not is_convex([(0, 0), (0, 1), (1, 0), (1, 1)]) + + assert is_convex([(0, 0), (0, 0.5), (0, 1), (1, 1), (1, 0)]) + assert not is_convex([(0, 0), (0.1, 0.5), (0, 1), (1, 1), (1, 0)]) + + +def test_triangle_convex_polygon(): + assert triangle_convex_polygon([(0, 0), (0, 1), (1, 1), (1, 0)]) == [(0, 1, 2), (0, 2, 3)] + assert triangle_convex_polygon([(0, 0), (0, 0.5), (0, 1), (1, 1), (1, 0)]) == [(0, 2, 3), (0, 3, 4)] + + +@pytest.mark.parametrize( + ('polygon', 'expected'), + [ + ([(0, 0), (0, 1), (1, 1), (1, 0)], [(0, 1, 2), (0, 3, 2)]), + ([(0, 0), (0, 10), (1, 10), (1, 0)], [(0, 1, 2), (0, 3, 2)]), + ([(0, 0), (0, 0.5), (0, 1), (1, 1), (1, 0)], [(0, 2, 3), (0, 3, 4)]), + ], +) +def test_triangulate_polygon_py_convex(polygon, expected): + assert triangulate_polygon_py(polygon)[0] == expected + + +def _renumerate_triangles(polygon, points, triangles): + point_num = {point: i for i, point in enumerate(polygon)} + return [tuple(point_num[tuple(points[point])] for point in triangle) for triangle in triangles] + + +TEST_POLYGONS = [ + ([(0, 0), (1, 1), (0, 2), (2, 1)], [(3, 2, 1), (0, 3, 1)]), + ([(0, 0), (0, 1), (1, 2), (2, 1), (2, 0), (1, 0.5)], [(4, 3, 5), (3, 2, 1), (5, 3, 1), (5, 1, 0)]), + ([(0, 1), (0, 2), (1, 1.5), (2, 2), (2, 1), (1, 0.5)], [(4, 3, 2), (2, 1, 0), (4, 2, 0), (5, 4, 0)]), + ([(0, 1), (0, 2), (1, 0.5), (2, 2), (2, 1), (1, -0.5)], [(2, 1, 0), (2, 0, 5), (4, 3, 2), (5, 4, 2)]), + ([(0, 0), (1, 2), (2, 0), (1, 1)], [(2, 1, 3), (3, 1, 0)]), + ([(0, 0), (0, 1), (0.5, 0.5), (1, 0), (1, 1)], [(3, 4, 2), (2, 1, 0)]), + ([(0, 0), (1, 0), (0.5, 0.5), (0, 1), (1, 1)], [(2, 4, 3), (1, 2, 0)]), +] + + +@pytest.mark.parametrize(('polygon', 'expected'), TEST_POLYGONS) +def test_triangulate_polygon_py_non_convex(polygon, expected): + triangles, points = triangulate_polygon_py(polygon) + triangles_ = _renumerate_triangles(polygon, points, triangles) + assert triangles_ == expected + + +@pytest.mark.parametrize(('polygon', 'expected'), TEST_POLYGONS) +def test_triangulate_polygon_numpy_non_convex(polygon, expected): + triangles, points = triangulate_polygon_numpy(np.array(polygon).astype(np.float64)) + triangles_ = _renumerate_triangles(polygon, points, triangles) + assert triangles_ == expected + + +@pytest.mark.parametrize(('polygon', 'expected'), TEST_POLYGONS) +def test_triangulate_polygon_numpy_li_non_convex(polygon, expected): + triangles, points = triangulate_polygon_numpy_li([np.array(polygon).astype(np.float64)]) + triangles_ = _renumerate_triangles(polygon, points, triangles) + assert triangles_ == expected + + +def test_triangulate_polygon_in_polygon_numpy(): + polygons = [ + np.array([(0, 0), (10, 0), (10, 10), (0, 10)], dtype=np.float64), + np.array([(4, 4), (6, 4), (6, 6), (4, 6)], dtype=np.float32), + ] + triangles, points = triangulate_polygon_numpy_li(polygons) + assert len(triangles) == 8 + assert len(points) == 8 + + +def test_triangulate_polygon_segfault1(): + """Test on polygon that lead to segfault during test""" + polygon = [ + (205.0625, 1489.83752), + (204.212509, 1490.4751), + (204, 1491.11255), + (202.087509, 1493.45007), + (201.875, 1494.7251), + (202.300003, 1496), + (202.300003, 1498.33752), + (203.575012, 1499.82507), + (204.425003, 1500.25), + (205.0625, 1500.25), + (205.700012, 1500.67505), + (206.550003, 1500.67505), + (207.1875, 1500.25), + (208.037506, 1500.88757), + (209.3125, 1499.82507), + (209.525009, 1499.1875), + (211.012512, 1497.70007), + (210.375, 1496.42505), + (209.525009, 1495.57507), + (208.462509, 1495.15002), + (208.675003, 1494.9375), + (208.462509, 1492.8125), + (208.037506, 1491.5376), + (205.912506, 1489.83752), + ] + triangulate_polygon_py(polygon) + + +def test_triangulate_polygon_segfault2(): + polygon = [ + [1388.6875, 2744.4375], + [1388.4751, 2744.6501], + [1386.5625, 2744.6501], + [1385.925, 2744.8625], + [1385.5, 2745.2876], + [1385.2876, 2747.625], + [1385.7125, 2748.2627], + [1385.7125, 2749.1125], + [1386.1376, 2749.75], + [1389.9625, 2753.7876], + [1390.3876, 2754.6377], + [1391.025, 2754.6377], + [1392.0875, 2753.1501], + [1392.3, 2753.3625], + [1392.3, 2754.6377], + [1392.5126, 2754.2126], + [1392.3, 2754.0], + [1392.3, 2751.4502], + [1392.7251, 2750.3877], + [1391.6626, 2748.9001], + [1391.6626, 2747.4126], + [1390.8125, 2745.5], + [1390.175, 2745.2876], + [1389.3251, 2744.4375], + ] + triangulate_polygon_py(polygon) + + +def test_triangulate_polygon_segfault3(): + polygon = [ + (1066.32507, 1794.3501), + (1065.6875, 1794.77502), + (1063.77502, 1794.77502), + (1063.5625, 1794.98755), + (1063.13757, 1794.98755), + (1062.28748, 1795.83752), + (1062.28748, 1797.32507), + (1062.07507, 1797.5376), + (1062.28748, 1797.5376), + (1062.5, 1797.75), + (1063.13757, 1797.75), + (1063.5625, 1797.32507), + (1064.625, 1797.32507), + (1065.26257, 1797.96252), + (1064.83752, 1797.5376), + (1064.83752, 1796.90002), + (1065.47498, 1796.26257), + (1065.47498, 1796.05005), + (1065.6875, 1795.83752), + (1066.53748, 1795.83752), + (1066.75, 1795.625), + (1066.75, 1794.98755), + (1066.96252, 1794.77502), + (1066.75, 1794.3501), + ] + triangulate_polygon_py(polygon) + + +def test_triangulate_polygon_segfault4(): + polygon = [ + [657.6875, 2280.975], + [657.6875, 2281.6125], + [657.05005, 2282.25], + [656.2, 2284.1626], + [657.6875, 2285.8625], + [658.11255, 2286.7126], + [659.8125, 2288.4126], + [659.8125, 2288.625], + [661.9375, 2290.5376], + [662.78754, 2290.9626], + [664.7, 2292.2375], + [665.3375, 2292.2375], + [665.97504, 2291.175], + [666.61255, 2290.5376], + [666.61255, 2289.6875], + [666.1875, 2288.625], + [664.48755, 2286.925], + [664.0625, 2286.925], + [663.2125, 2286.5], + [661.9375, 2284.8], + [660.66254, 2284.8], + [660.45, 2284.5876], + [660.45, 2284.1626], + [657.6875, 2281.1875], + ] + triangulate_polygon_py(polygon) + + +def test_triangulate_polygon_segfault5(): + polygon = [ + [895.05005, 2422.5], + [894.8375, 2422.7126], + [894.41254, 2422.7126], + [893.98755, 2423.1375], + [893.35004, 2423.1375], + [892.92505, 2423.5625], + [892.7125, 2423.5625], + [891.4375, 2424.8376], + [891.22504, 2424.8376], + [892.075, 2424.8376], + [892.5, 2425.05], + [893.35004, 2425.9001], + [893.775, 2426.75], + [893.775, 2427.3875], + [894.625, 2426.75], + [895.6875, 2426.75], + [896.11255, 2426.1125], + [896.11255, 2425.05], + [895.9, 2424.8376], + [896.53754, 2423.7751], + [896.53754, 2423.35], + [896.75, 2422.925], + [896.11255, 2422.7126], + [895.9, 2422.5], + ] + triangulate_polygon_py(polygon) + + +@pytest.mark.parametrize( + ('segment1', 'segment2'), + [ + # Touching segments by top point + (((0, 0), (1, 1)), ((1, 1), (2, 0))), + # Touching segments by bottom point + (((0, 1), (1, 0)), ((1, 0), (2, 1))), + # Touching segments by bottom point with different length, longer left + (((-1, 2), (1, 0)), ((1, 0), (2, 1))), + # Touching segments by bottom point with different length, longer right + (((0, 1), (1, 0)), ((1, 0), (3, 2))), + # Parallel vertical segments + (((0, 0), (0, 1)), ((1, 0), (1, 1))), + # Parallel horizontal segments different length + (((0, 0), (0, 2)), ((1, -1), (1, 1))), + # horizontal segments + (((0, 0), (1, 0)), ((2, 0), (3, 0))), + # Two segments with top on same line + (((0, 0), (1, 1)), ((3, 1), (2, -1))), + # One horizontal segment on right and top + (((0, 0), (1, 1)), ((2, 1), (3, 1))), + # One horizontal segment on right and middle + (((0, 0), (1, 1)), ((2, 0.5), (3, 0.5))), + # One horizontal segment on right and bottom + (((0, 0), (1, 1)), ((2, 0), (3, 0))), + # One horizontal segment on left and top + (((0, 1), (1, 1)), ((2, 0), (3, 2))), + # One horizontal segment on left and middle + (((0, 0.5), (1, 0.5)), ((2, 0), (3, 2))), + # One horizontal segment on left and bottom + (((0, 0), (1, 0)), ((2, 0), (3, 2))), + # left horizontal, right oblique, bottom merge + (((0, 0), (1, 0)), ((1, 0), (2, 2))), + # left oblique, right horizontal, bottom merge + (((0, 1), (1, 0)), ((1, 0), (2, 0))), + # left horizontal, right oblique, top merge + (((0, 1), (1, 1)), ((1, 1), (2, 0))), + # left oblique, right horizontal, top merge + (((0, 0), (1, 1)), ((1, 1), (2, 1))), + ], + ids=[ + 'Touching segments by top point', + 'Touching segments by bottom point', + 'Touching segments by bottom point with different length, longer left', + 'Touching segments by bottom point with different length, longer right', + 'Parallel vertical segments', + 'Parallel horizontal segments different length', + 'horizontal segments', + 'Two segments with top on same line', + 'One horizontal segment on right and top', + 'One horizontal segment on right and middle', + 'One horizontal segment on right and bottom', + 'One horizontal segment on left and top', + 'One horizontal segment on left and middle', + 'One horizontal segment on left and bottom', + 'left horizontal, right oblique, bottom merge', + 'left oblique, right horizontal, bottom merge', + 'left horizontal, right oblique, top merge', + 'left oblique, right horizontal, top merge', + ], +) +def test_segment_left_to_right_comparator(segment1, segment2): + assert segment_left_to_right_comparator(segment1, segment2) + assert not segment_left_to_right_comparator(segment2, segment1) + + +def test_find_intersection_points_py_cross(): + r""" + (1, 0) --- (1, 1) + \ / + \ / + \ / + X + / \ + / \ + / \ + (0, 0) --- (0, 1) + """ + assert find_intersection_points_py([(0, 0), (1, 1), (1, 0), (0, 1)]) == [ + (0, 0), + (0.5, 0.5), + (1, 1), + (1, 0), + (0.5, 0.5), + (0, 1), + ] + + +def test_find_intersection_points_py_cross_intersect_in_point(): + r""" + (1, 0) --- (1, 1) + \ / + \ / + \ / + X + / \ + / \ + / \ + (0, 0) --- (0, 1) + """ + assert find_intersection_points_py([(0, 0), (0.5, 0.5), (1, 1), (1, 0), (0, 1)]) == [ + (0, 0), + (0.5, 0.5), + (1, 1), + (1, 0), + (0.5, 0.5), + (0, 1), + ] + + +@pytest.mark.parametrize( + ('polygon', 'expected'), + [ + ( + ((1, 2), (1, 0), [(0, 1)], [(2, 1)]), + [[(2.0, 1.0), (1.0, 2.0), (0.0, 1.0)], [(1.0, 0.0), (2.0, 1.0), (0.0, 1.0)]], + ), + ( + ((5, 2), (5, 0), [(0, 1)], [(2, 1)]), + [[(2.0, 1.0), (5.0, 2.0), (0.0, 1.0)], [(5.0, 0.0), (2.0, 1.0), (0.0, 1.0)]], + ), + ( + ((1, 3), (1, 0), [(0, 2), (0, 1)], [(2, 2), (2, 1)]), + [ + [(2.0, 2.0), (1.0, 3.0), (0.0, 2.0)], + [(2.0, 1.0), (2.0, 2.0), (0.0, 2.0)], + [(2.0, 1.0), (0.0, 2.0), (0.0, 1.0)], + [(1.0, 0.0), (2.0, 1.0), (0.0, 1.0)], + ], + ), + ( + ((0, 4), (0, 0), [], [(2, 3), (3, 2), (2, 1)]), + [ + [(3.0, 2.0), (2.0, 3.0), (0.0, 4.0)], + [(2.0, 1.0), (3.0, 2.0), (0.0, 4.0)], + [(2.0, 1.0), (0.0, 4.0), (0.0, 0.0)], + ], + ), + ( + ((0, 4), (0, 0), [], [(2, 3), (1, 2), (2, 1)]), + [ + [(1.0, 2.0), (2.0, 3.0), (0.0, 4.0)], + [(1.0, 2.0), (0.0, 4.0), (0.0, 0.0)], + [(2.0, 1.0), (1.0, 2.0), (0.0, 0.0)], + ], + ), + ], +) +def test_triangulate_monotone_polygon_py(polygon, expected): + assert triangulate_monotone_polygon_py(*polygon) == expected + + +@pytest.mark.parametrize( + ('path', 'closed', 'bevel', 'expected', 'exp_triangles'), + [ + ( + [[0, 0], [0, 10], [10, 10], [10, 0]], + True, + False, + 10, + [[2, 1, 0], [1, 2, 3], [4, 3, 2], [3, 4, 5], [6, 5, 4], [5, 6, 7], [8, 7, 6], [7, 8, 9]], + ), + ( + [[0, 0], [0, 10], [10, 10], [10, 0]], + False, + False, + 8, + [[2, 1, 0], [1, 2, 3], [4, 3, 2], [3, 4, 5], [6, 5, 4], [5, 6, 7]], + ), + ( + [[0, 0], [0, 10], [10, 10], [10, 0]], + True, + True, + 14, + [ + [2, 1, 0], + [3, 2, 0], + [2, 3, 4], + [5, 4, 3], + [6, 5, 3], + [5, 6, 7], + [8, 7, 6], + [9, 8, 6], + [8, 9, 10], + [11, 10, 9], + [12, 11, 9], + [11, 12, 13], + ], + ), + ( + [[0, 0], [0, 10], [10, 10], [10, 0]], + False, + True, + 10, + [[2, 1, 0], [1, 2, 3], [4, 3, 2], [5, 4, 2], [4, 5, 6], [7, 6, 5], [8, 7, 5], [7, 8, 9]], + ), + ( + [[2, 10], [0, -5], [-2, 10], [-2, -10], [2, -10]], + True, + False, + 15, + [ + [2, 1, 0], + [1, 2, 3], + [1, 3, 4], + [5, 4, 3], + [6, 5, 3], + [5, 6, 7], + [8, 7, 6], + [7, 8, 9], + [7, 9, 10], + [11, 10, 9], + [10, 11, 12], + [13, 12, 11], + [12, 13, 14], + ], + ), + ([[0, 0], [0, 10]], False, False, 4, [[2, 1, 0], [1, 2, 3]]), + ([[0, 0], [0, 10], [0, 20]], False, False, 6, [[2, 1, 0], [1, 2, 3], [4, 3, 2], [3, 4, 5]]), + ( + [[0, 0], [0, 2], [10, 1]], + True, + False, + 9, + [[2, 1, 0], [1, 2, 3], [4, 3, 2], [3, 4, 5], [6, 5, 4], [7, 6, 4], [6, 7, 8]], + ), + ([[0, 0], [10, 1], [9, 1.1]], False, False, 7, [[2, 1, 0], [1, 2, 3], [4, 3, 2], [3, 4, 5], [3, 5, 6]]), + ([[9, 0.9], [10, 1], [0, 2]], False, False, 7, [[2, 1, 0], [1, 2, 3], [4, 3, 2], [3, 4, 5], [3, 5, 6]]), + ([[0, 0], [-10, 1], [-9, 1.1]], False, False, 7, [[2, 1, 0], [1, 2, 3], [4, 3, 2], [5, 4, 2], [4, 5, 6]]), + ([[-9, 0.9], [-10, 1], [0, 2]], False, False, 7, [[2, 1, 0], [1, 2, 3], [4, 3, 2], [5, 4, 2], [4, 5, 6]]), + ], +) +def test_triangulate_path_edge_py(path, closed, bevel, expected, exp_triangles): + centers, offsets, triangles = triangulate_path_edge_py(np.array(path, dtype='float32'), closed=closed, bevel=bevel) + assert centers.shape == offsets.shape + assert centers.shape[0] == expected + assert triangles.shape[0] == expected - 2 + triangles_li = [[int(y) for y in x] for x in triangles] + assert triangles_li == exp_triangles + + +@pytest.mark.parametrize(('polygon', 'expected'), TEST_POLYGONS) +def test_triangulate_polygon_with_edge_numpy_li(polygon, expected): + (triangles, points), (centers, offsets, edge_triangles) = triangulate_polygon_with_edge_numpy_li( + [np.array(polygon)] + ) + triangles_ = _renumerate_triangles(polygon, points, triangles) + assert triangles_ == expected + assert centers.shape == offsets.shape diff --git a/tox.ini b/tox.ini index eeb4073..04b6f41 100644 --- a/tox.ini +++ b/tox.ini @@ -28,4 +28,4 @@ deps = cython commands = - pytest + pytest -v